Estimation-Based SpGEMM Workflow

• Estimation-based SpGEMM workflow is a family of techniques that predict sparse matrix product structures to optimize memory allocation and execution scheduling.
• It leverages lightweight predictors—such as sampled compression ratios, work estimates, and HLL sketches—to accurately forecast nonzero distributions and guide algorithm dispatch.
• By integrating estimation into memory pre-allocation and kernel selection, these workflows enhance performance across diverse architectures including CPUs, GPUs, vector, and accelerator-based systems.

Searching arXiv for relevant papers on estimation-based SpGEMM workflows and related algorithm-selection work. Estimation-based SpGEMM workflow denotes a family of sparse general matrix multiplication procedures in which structural properties of the product matrix are predicted or estimated before, or in lieu of, exact symbolic computation. In sparse matrix multiplication $C = A \cdot B$, the principal target of estimation is typically the output structure of $C$, especially the number of nonzeros per row or per column, because memory allocation, load balancing, and accumulator selection depend on it. Across CPU, GPU, vector, and accelerator settings, the workflow replaces or supplements exact symbolic phases with lightweight predictors derived from intermediate-product counts, sampled compression ratios, HyperLogLog sketches, or work metrics such as per-column operation counts. The approach is motivated by the cost of exact symbolic traversal and by the fact that SpGEMM performance is strongly conditioned on irregular sparsity, accumulator pressure, and communication patterns (Du et al., 2022, Li et al., 21 Apr 2026, Liu et al., 2015, Fèvre et al., 2023, Ballard et al., 2016, Yadav et al., 2024).

1. Definition and motivating problem

Sparse general matrix multiplication (SpGEMM) computes $C = A \cdot B$ for sparse matrices and is a fundamental building block in scientific computing, graph analytics, and machine learning (Du et al., 2022, Fèvre et al., 2023, Li et al., 21 Apr 2026, Liu et al., 2015). A central difficulty is that the output sparsity pattern is unknown a priori. In CSR-based row-wise formulations, the critical quantity is the number of nonzeros per output row of $C$; in CSC-based column-wise vector formulations, the analogous quantity is the work or structural complexity per output column (Du et al., 2022, Fèvre et al., 2023).

The motivation for estimation is operational rather than merely analytical. Accurate prediction enables memory allocation for CSR row pointers and per-row buffers, load balancing across CPU threads or GPU warps and CTAs, and kernel selection among hash-based, merge-based, sorting-based, dense, or hybrid accumulators (Du et al., 2022, Li et al., 21 Apr 2026). In vector architectures, work estimation supports dynamic switching among SPA, SPARS, and HASH, while in communication-centric settings it guides the selection of 1D, 2D, or 3D algorithm families (Fèvre et al., 2023, Ballard et al., 2016). In accelerator-oriented dataflow selection, estimation extends beyond output size to operational and traffic surrogates that distinguish inner-product, outer-product, and row-wise dataflows (Yadav et al., 2024).

The immediate alternative to estimation is exact symbolic computation. Exact symbolic methods determine the output structure precisely but often perform a traversal similar in complexity to the numeric phase, thereby adding substantial overhead when used strictly for sizing or scheduling (Du et al., 2022). On GPUs, the symbolic pass alone accounts for roughly $28\%$ of total runtime on average in state-of-the-art baselines such as spECK, even though it produces only one number per row (Li et al., 21 Apr 2026). Upper bounds based on intermediate products are cheap, but can be extremely loose, particularly when the compression ratio is high (Du et al., 2022). Estimation-based workflows occupy the design space between these two extremes.

2. Structural quantities and estimator targets

The quantities estimated in SpGEMM workflows vary by execution model, but they all describe latent structure or work in the product.

For row-wise CSR SpGEMM, let $A \in \mathbb{R}^{m \times k}$, $B \in \mathbb{R}^{k \times n}$, and $C = A \cdot B$. The row-wise dataflow forms

$$C_{i,:} = \sum_{k \in \mathrm{NZ}(A_{i,:})} A_{i,k} \cdot B_{k,:}.$$

A fundamental quantity is the row-level intermediate-product count

$$\mathrm{FLOP}_i = \sum_{k \in \mathrm{NZ}(A_{i,:})} \mathrm{nnz}(B_{k,:}),$$

with total

$C$0

The corresponding output compression ratio is

$C$1

These definitions underpin sampled compression-ratio estimation on CPUs and upper-bound workflows on GPUs (Du et al., 2022, Liu et al., 2015, Li et al., 21 Apr 2026).

In the GPU formulation of Ocean, the structural target for row $C$2 is

$C$3

the number of distinct column indices reachable by two-hop expansion through the nonzeros of row $C$4 of $C$5 into rows of $C$6 (Li et al., 21 Apr 2026). This is precisely the quantity needed for per-row capacity and bin assignment in Gustavson-style numeric kernels.

In the vector-architecture workflow using CSC storage, the principal estimator is per-column work:

$C$7

where $C$8 is the number of nonzeros in column $C$9 of $C = A \cdot B$0 (Fèvre et al., 2023). Here the estimation target is not directly $C = A \cdot B$1 but the intermediate-product count governing whether SPA, SPARS, or HASH is preferable.

The 2015 GPU framework uses

$C = A \cdot B$2

which serves simultaneously as a work estimate and an upper bound on $C = A \cdot B$3 via $C = A \cdot B$4 (Liu et al., 2015). This bound is structurally loose but easy to compute in $C = A \cdot B$5 time.

A broader implication is that estimation-based workflows can be organized around three estimator families: direct size prediction, work prediction, and communication prediction. That taxonomy is not stated explicitly in the sources, but it is consistent with the sampled-$C = A \cdot B$6 method (Du et al., 2022), the $C = A \cdot B$7 vector work model (Fèvre et al., 2023), and the hypergraph communication framework (Ballard et al., 2016).

3. Ratio-based estimation of output structure

A prominent estimation-based workflow is the sampled compression-ratio method for CSR SpGEMM proposed in "Predicting the Output Structure of Sparse Matrix Multiplication with Sampled Compression Ratio" (Du et al., 2022). Its central idea is to estimate $C = A \cdot B$8 from a sampled result computed on a small uniform random subset of rows of $C = A \cdot B$9, while keeping all rows of $C$0 visible. For a sampled row set $C$1,

$C$2

and the sample-based compression-ratio estimator is

$C$3

Predicted row lengths then follow as

$C$4

with total prediction

$C$5

The workflow computes $C$6 by precomputing the row degrees $C$7 and summing them over the nonzeros in each row of $C$8 (Du et al., 2022).

The distinguishing argument is error neutralization. If the sampled FLOP and sampled nonzero counts are perturbed by relative errors $C$9 and $28\%$0, then the induced relative error in the estimated total is

$28\%$1

When $28\%$2, the estimator error becomes small. The paper reports a positive correlation between sampled FLOP and sampled nnz with correlation coefficient $28\%$3 across $28\%$4 test cases, which supports the ratio-based design (Du et al., 2022).

The sampling strategy uses a small, uniform random subset of rows of $28\%$5 with

$28\%$6

and computes $28\%$7 exactly on the sample using a per-row hash-based accumulator with linear probing (Du et al., 2022). Because both matrices are typically in CSR, row sampling is preferred over column sampling, and row-wise dataflow is described as faster and simpler than inner-product dataflow for CSR inputs (Du et al., 2022).

The empirical results are specific. On $28\%$8 test cases derived from $28\%$9 representative SuiteSparse matrices, the absolute relative error of the proposed method is $A \in \mathbb{R}^{m \times k}$0 on average and $A \in \mathbb{R}^{m \times k}$1 in the worst case, versus $A \in \mathbb{R}^{m \times k}$2 average and $A \in \mathbb{R}^{m \times k}$3 worst-case for a reference sampling design that scales sampled nnz only and ignores sampled FLOP (Du et al., 2022). The parallel implementation adds on average $A \in \mathbb{R}^{m \times k}$4 of execution time relative to BRMerge-Precise; on matrix-square benchmarks, computing $A \in \mathbb{R}^{m \times k}$5 costs $A \in \mathbb{R}^{m \times k}$6 on average, up to $A \in \mathbb{R}^{m \times k}$7, and the prediction step costs $A \in \mathbb{R}^{m \times k}$8 on average, up to $A \in \mathbb{R}^{m \times k}$9 (Du et al., 2022).

The workflow then integrates prediction into downstream execution. It allocates CSR with

$B \in \mathbb{R}^{k \times n}$0

where $B \in \mathbb{R}^{k \times n}$1 is a small guard factor such as $B \in \mathbb{R}^{k \times n}$2–$B \in \mathbb{R}^{k \times n}$3, and uses predicted row sizes to choose hash-based, merge-based, or sorting-based accumulation and to distribute rows proportionally to $B \in \mathbb{R}^{k \times n}$4 or $B \in \mathbb{R}^{k \times n}$5 (Du et al., 2022).

4. Upper-bound and hybrid pre-allocation workflows

An earlier and influential GPU framework estimates structure using per-row work bounds rather than direct output-size prediction (Liu et al., 2015). In "A Framework for General Sparse Matrix-Matrix Multiplication on GPUs and Heterogeneous Processors" (Liu et al., 2015), all stages use CSR, and the symbolic estimation phase computes $B \in \mathbb{R}^{k \times n}$6 on the GPU with one thread per output row. Because $B \in \mathbb{R}^{k \times n}$7 depends only on row lengths read from $B \in \mathbb{R}^{k \times n}$8’s row pointers, the estimation cost is $B \in \mathbb{R}^{k \times n}$9 (Liu et al., 2015).

The framework then performs host-side binning into $C = A \cdot B$0 bins organized into five groups:

Bin group	Range of $C = A \cdot B$1	Allocation policy
1	$C = A \cdot B$2	Empty rows
2	$C = A \cdot B$3	Single-entry rows
3	$C = A \cdot B$4	Exact-size upper bound per row
4	$C = A \cdot B$5	Exact-size upper bound per row
5	$C = A \cdot B$6	Initial capacity $C = A \cdot B$7, then progressive growth

For bin groups $C = A \cdot B$8–$C = A \cdot B$9, the temporary CSR matrix $$C_{i,:} = \sum_{k \in \mathrm{NZ}(A_{i,:})} A_{i,k} \cdot B_{k,:}.$$0 allocates $$C_{i,:} = \sum_{k \in \mathrm{NZ}(A_{i,:})} A_{i,k} \cdot B_{k,:}.$$1. For bin group $$C_{i,:} = \sum_{k \in \mathrm{NZ}(A_{i,:})} A_{i,k} \cdot B_{k,:}.$$2, capacity is initialized to $$C_{i,:} = \sum_{k \in \mathrm{NZ}(A_{i,:})} A_{i,k} \cdot B_{k,:}.$$3 and grown by doubling when needed (Liu et al., 2015). This is described as a hybrid strategy interpolating between upper-bound allocation and progressive allocation. The rationale is explicit: full upper-bound allocation at the intermediate size $$C_{i,:} = \sum_{k \in \mathrm{NZ}(A_{i,:})} A_{i,k} \cdot B_{k,:}.$$4 can exceed device memory by orders of magnitude on irregular matrices, whereas precise pre-computation doubles the work and time (Liu et al., 2015).

The estimator is integral to algorithm dispatch. Rows with $$C_{i,:} = \sum_{k \in \mathrm{NZ}(A_{i,:})} A_{i,k} \cdot B_{k,:}.$$5 use a heap-based shared-memory accumulator. Rows with $$C_{i,:} = \sum_{k \in \mathrm{NZ}(A_{i,:})} A_{i,k} \cdot B_{k,:}.$$6 use bitonic ESC (Expansion–Sorting–Compression). Rows with $$C_{i,:} = \sum_{k \in \mathrm{NZ}(A_{i,:})} A_{i,k} \cdot B_{k,:}.$$7 use merge-based insertion with GPU merge path, with checkpointing, dumping to global memory, re-allocation by $$C_{i,:} = \sum_{k \in \mathrm{NZ}(A_{i,:})} A_{i,k} \cdot B_{k,:}.$$8, and reload for long rows (Liu et al., 2015). Load balancing follows the same estimate: short rows use one thread per row, moderate and long rows use one thread group per row, and per-bin counters avoid launching empty kernels (Liu et al., 2015).

This workflow does not estimate $$C_{i,:} = \sum_{k \in \mathrm{NZ}(A_{i,:})} A_{i,k} \cdot B_{k,:}.$$9 tightly. Instead, it treats arithmetic work as a safe upper bound and uses bin-specific accumulation strategies to absorb overestimation. A plausible implication is that it belongs to the same estimation-based lineage as later ratio-based and HLL-based methods, but emphasizes robust scheduling and memory control over prediction accuracy. That interpretation is consistent with the later critique that upper bounds can be “extremely loose” when compression ratios are large (Du et al., 2022).

The framework also includes a heterogeneous-processor variant using re-allocatable shared virtual memory on APUs. For matrices that need re-allocation, this yields average $$\mathrm{FLOP}_i = \sum_{k \in \mathrm{NZ}(A_{i,:})} \mathrm{nnz}(B_{k,:}),$$0 and up to $$\mathrm{FLOP}_i = \sum_{k \in \mathrm{NZ}(A_{i,:})} \mathrm{nnz}(B_{k,:}),$$1 speedup on AMD A10-7850K by enabling in-place growth and shrink without copy (Liu et al., 2015).

5. Sketch-based estimation and workflow selection on GPUs

The GPU system Ocean introduces a more elaborate estimation-based symbolic workflow that replaces most exact symbolic accumulation with HyperLogLog (HLL) sketches (Li et al., 21 Apr 2026). Ocean constructs one HLL sketch per row of $$\mathrm{FLOP}_i = \sum_{k \in \mathrm{NZ}(A_{i,:})} \mathrm{nnz}(B_{k,:}),$$2, using $$\mathrm{FLOP}_i = \sum_{k \in \mathrm{NZ}(A_{i,:})} \mathrm{nnz}(B_{k,:}),$$3 registers with $$\mathrm{FLOP}_i = \sum_{k \in \mathrm{NZ}(A_{i,:})} \mathrm{nnz}(B_{k,:}),$$4 bits per register. It uses $$\mathrm{FLOP}_i = \sum_{k \in \mathrm{NZ}(A_{i,:})} \mathrm{nnz}(B_{k,:}),$$5 adaptively: $$\mathrm{FLOP}_i = \sum_{k \in \mathrm{NZ}(A_{i,:})} \mathrm{nnz}(B_{k,:}),$$6 for $$\mathrm{FLOP}_i = \sum_{k \in \mathrm{NZ}(A_{i,:})} \mathrm{nnz}(B_{k,:}),$$7 when the input expansion ratio $$\mathrm{FLOP}_i = \sum_{k \in \mathrm{NZ}(A_{i,:})} \mathrm{nnz}(B_{k,:}),$$8, otherwise $$\mathrm{FLOP}_i = \sum_{k \in \mathrm{NZ}(A_{i,:})} \mathrm{nnz}(B_{k,:}),$$9 for $C$00 (Li et al., 21 Apr 2026).

For each row $C$01 of $C$02, Ocean merges the sketches of the rows of $C$03 referenced by $C$04 using elementwise maxima and estimates the union cardinality with the standard raw HLL estimator

$C$05

where $C$06 (Li et al., 21 Apr 2026). Linear counting is used when $C$07, and large-range corrections are enabled when $C$08 (Li et al., 21 Apr 2026).

Ocean embeds HLL into a matrix-level workflow selector driven by two ratios. The input expansion ratio is

$C$09

and the output compression ratio is

$C$10

Ocean estimates $C$11 cheaply via sampling: it builds HLL sketches for all rows of $C$12 and randomly samples $C$13 of rows of $C$14, with minimum $C$15 and maximum $C$16, to merge their sketches and estimate approximate $C$17 values (Li et al., 21 Apr 2026).

The workflow then chooses among three modes:

Workflow	Selection criterion
Upper-bound estimation	Average intermediate products per row $C$18
HLL estimation	Average intermediate products per row $C$19, $C$20, and sampled $C$21
Exact symbolic	Otherwise

The rationale is architectural. Each HLL register is a byte, each CSR column index is $C$22 bytes, and $C$23 makes merging $C$24 registers per referenced row about as cheap as reading approximately $C$25 symbolic indices (Li et al., 21 Apr 2026). Likewise, $C$26 is used because estimation-based workflows delay the final CSR layout, and late compaction overhead grows when compression is low (Li et al., 21 Apr 2026).

Ocean treats estimator uncertainty through allocation slack. Per-row capacity is

$C$27

with default $C$28, and $C$29 may be raised to $C$30 when $C$31 (Li et al., 21 Apr 2026). The relative standard error of HLL is approximated by $C$32, yielding about $C$33 for $C$34 and $C$35 for $C$36 (Li et al., 21 Apr 2026). Measured average relative per-row errors are about $C$37, $C$38, and $C$39 for $C$40, with overflow rates after binning with $C$41 expansion and upward rounding of $C$42, $C$43, and less than $C$44, respectively (Li et al., 21 Apr 2026).

The numeric phase uses the estimator to guide accumulator choice. Ocean provides hash-based accumulators with hybrid memory, dense accumulators, and ESC accumulators, each with multiple bins (Li et al., 21 Apr 2026). Dense is preferred over hash when both fit the same bin, except in the largest bin where hybrid hash is favored; ESC is used only in the upper-bound workflow (Li et al., 21 Apr 2026). Overflow rows are handled by a fallback kernel using the largest dense accumulator configuration and a global buffer sized by the row’s intermediate product count (Li et al., 21 Apr 2026).

The performance data are end-to-end. On $C$45 square matrices and $C$46 rectangular matrices from SuiteSparse, Ocean achieves on A100 a geometric mean of $C$47 GFLOP/s for square cases and $C$48 GFLOP/s for rectangular cases, and on H100 $C$49 GFLOP/s for square cases (Li et al., 21 Apr 2026). It wins on $C$50 square matrices on A100 and $C$51 on H100, with average speedups of $C$52–$C$53 over leading GPU baselines across the study; more specifically, on A100 square matrices the average speedups over cuSPARSE, spECK, opSparse, TileSpGEMM, and HSMU are $C$54, $C$55, $C$56, $C$57, and $C$58, respectively (Li et al., 21 Apr 2026). The estimation step costs approximately $C$59 of runtime for selected matrices, compared with symbolic-step costs of about $C$60–$C$61 in spECK or the baseline Ocean variant (Li et al., 21 Apr 2026).

6. Estimation for kernel, dataflow, and communication selection

Estimation-based SpGEMM workflows are not confined to output-size prediction. In several settings, the estimated quantity is used to select among algorithmic dataflows, accumulator families, or communication decompositions.

On RISC-V vector architectures, the workflow in "Optimization of SpGEMM with Risc-V vector instructions" (Fèvre et al., 2023) computes $C$62 for every output column and sorts columns by decreasing work. This work estimate drives dynamic switching between SPA, SPARS, and HASH, yielding hybrid algorithms H-SPA$C$63 and H-HASH$C$64 governed by a threshold $C$65 (Fèvre et al., 2023). If $C$66, the system uses SPARS or HASH inside blocks formed to homogenize column work; otherwise it uses SPA (Fèvre et al., 2023). The threshold guidance comes from sensitivity experiments: $C$67 in $C$68 yields the best average, and $C$69 is used in the main evaluation (Fèvre et al., 2023). H-SPA$C$70 and H-HASH$C$71 obtain average speedups of $C$72 and $C$73 over SPA on $C$74 SuiteSparse matrices, and $C$75 and $C$76 on the $C$77 most sparse matrices (Fèvre et al., 2023).

In Misam, estimation is fused with machine learning for dataflow selection among inner-product, outer-product, and row-wise SpGEMM (Yadav et al., 2024). The system extracts light-weight structural features such as sparsityA, sparsityB, avg_row_lengthA, avg_row_lengthA_var, avg_col_lengthB, and blocks_accessed, then predicts the best dataflow per block using a decision tree or DQN (Yadav et al., 2024). The decision tree, limited to depth $C$78, achieves $C$79 selection accuracy; the DQN achieves $C$80 (Yadav et al., 2024). Relative to any single fixed dataflow, the decision tree produces average speedups of $C$81 over IP, $C$82 over OP, and $C$83 over RW, with $C$84 improvement over the heuristic baseline (Yadav et al., 2024). The abstract reports gains of up to $C$85 times (Yadav et al., 2024). Although Misam frames the problem as dataflow choice rather than symbolic estimation, it exemplifies the same principle: low-cost structural surrogates are used to avoid committing to a single static execution strategy.

At a larger scale, the hypergraph framework for SpGEMM estimates communication and memory-hierarchy traffic from sparsity patterns (Ballard et al., 2016). The fine-grained model associates computation vertices with nontrivial products and data vertices with nonzeros of $C$86, $C$87, and $C$88, then encodes dependencies as hyperedges (Ballard et al., 2016). Estimation in this setting concerns communication volume along the critical path, total connectivity, and sequential I/O lower bounds. The workflow computes sparsity features, estimates flops $C$89, constructs coarse 1D or 2D hypergraph models, partitions or approximates them, and selects the algorithm family minimizing predicted communication subject to load balance (Ballard et al., 2016). The paper’s empirical conclusion is that at least one coarse-grained model is as communication efficient as the fine-grained model in each of the three application domains studied, but different applications favor different families (Ballard et al., 2016).

Taken together, these works suggest that “estimation-based SpGEMM workflow” is best understood as a general control architecture: estimate a structural proxy, map it to a resource or algorithm decision, then execute the numeric phase under that guidance. The proxy may be $C$90, $C$91, $C$92, feature vectors for ML classifiers, or hypergraph cuts.

7. Accuracy limits, guardrails, and open directions

The principal limitation of estimation-based workflows is that their reliability depends on how well the estimator captures overlap and skew in the product structure. In the sampled compression-ratio method, when the correlation between sampled FLOP and sampled nnz weakens, error neutralization diminishes; the observed worst-case error reaches $C$93 in extreme cases (Du et al., 2022). The paper identifies rows with extreme degree skew and duplicate aggregation patterns as contributors to residual error, and suggests adaptive sampling, importance sampling proportional to $C$94, and per-cluster compression-ratio estimation as future directions (Du et al., 2022).

Ocean incorporates explicit guardrails for such uncertainty (Li et al., 21 Apr 2026). It uses slack factors, upward bin rounding, and a robust overflow kernel. Under-estimation is therefore treated as a performance issue rather than a correctness issue: overflow rows are rerun with safe upper-bound capacity, while over-estimation increases compaction cost modestly but does not violate CSR correctness (Li et al., 21 Apr 2026). This suggests a design principle common to modern workflows: approximate sizing is acceptable when backed by bounded fallback mechanisms.

Upper-bound workflows have the opposite problem. They are safe but can over-allocate severely when the compression ratio is large, consuming memory and skewing binning or scheduling (Du et al., 2022, Liu et al., 2015). Ocean formalizes part of this trade-off with $C$95 and $C$96 thresholds, choosing exact symbolic whenever late compaction would not amortize or when HLL merging would not be traffic-efficient (Li et al., 21 Apr 2026).

Work-estimation workflows for vector architectures face different failure modes. The synthetic cross-overs reported for SPA versus SPARS and HASH depend on assumptions about uniform sparsity and sorted column degrees, and the paper notes that the general case deviates depending on column nnz distributions (Fèvre et al., 2023). Hence, $C$97 is a useful but not exhaustive predictor. Sorting and blocking reduce imbalance but do not eliminate it (Fèvre et al., 2023).

ML-based selectors such as Misam face hardware dependence and distribution shift (Yadav et al., 2024). The learned efficiencies and feature-action mappings are tied to buffer sizes, interconnect, and PE scheduling, so retraining per accelerator configuration is advisable (Yadav et al., 2024). That caveat is consistent with the broader literature: estimated structure is platform-agnostic only at a coarse level, while thresholding, accumulator choice, and occupancy interactions are platform-specific.

Several future directions recur across the sources. These include adaptive or stratified sampling for nonuniform compression ratios (Du et al., 2022); alternative sketches with lower relative standard error at equal bandwidth, hierarchical sketches, and better memory-footprint modeling for multi-GPU settings (Li et al., 21 Apr 2026); predictive models for kernel efficiency beyond shared-memory capacity on vector hardware (Fèvre et al., 2023); and tighter theoretical error bounds or uncertainty-aware fallbacks that bridge statistical estimation and resource provisioning (Du et al., 2022, Li et al., 21 Apr 2026). A plausible implication is that future estimation-based SpGEMM systems will combine multiple estimators—sampling, sketches, analytical bounds, and learned predictors—inside a unified controller rather than relying on a single proxy.

In aggregate, the literature defines an estimation-based SpGEMM workflow as a multistage procedure in which symbolic exactness is selectively replaced by predictive structure analysis, with the resulting estimates steering memory layout, load balancing, accumulator design, kernel dispatch, dataflow selection, and, in distributed settings, communication decomposition. The technical diversity of these workflows reflects the fact that the unknown output structure of sparse matrix products is not merely a combinatorial nuisance; it is the principal control variable of high-performance SpGEMM.