Spada: Accelerating Sparse Matrix Multiplication with Adaptive Dataflow

Which loop order wins at sparse-times-sparse matrix multiply depends on the matrix — and even on the region inside the matrix — so Spada makes the dataflow a dial, sliding a reshapeable window over the compressed rows and letting a few profiled passes per region choose the window's shape.

Why this problem is worth solving

Accelerators answered, but each froze one dataflow — one ordering of the three loops over M, N, and the reduction dimension K — into silicon:

• Inner product (reduction innermost): each output element is a dot product. Perfect output reuse, but the two operands must be intersected — most fetched pairs share no index in sparse data, so bandwidth is wasted on misses.
• Outer product (reduction outermost): each column of A times each row of B yields a partial matrix. Each B row is fetched exactly once — best input reuse — but the flood of
  partial sums — intermediate products that must later be merged into final output values
  can outgrow on-chip memory and spill.
• Row-based / Gustavson (reduction in the middle): build one output row at a time from the B rows matching that A row’s nonzeros. Psums stay small and on-chip, everything stays in CSR — but B rows are gathered irregularly, with reuse only if nearby A rows happen to share column indices.

The paper’s measurements make the consequence concrete: the inner-product design trails by 38× on average across the suite; the outer-product design wins compressed fully-connected layers while the row-based design wins convolutions and attention; and within single matrices, different regions prefer different shapes. A fixed dataflow is a standing bet against half your workload.

What people did before, and where it broke

SIGMA, OuterSPACE, SpArch, MatRaptor, GAMMA — each a strong point design, each rigid. (This site’s Flexagon explainer covers the closest prior spirit: supporting all three dataflows on one substrate, switched per DNN layer by an offline mapper.) Three gaps remained. First, the three classic dataflows have drastically different execution behaviors — inner product walks A by rows and B by columns, outer product the opposite — so naively unioning them in hardware means either replicated datapaths or heavyweight reconfiguration. Second, choosing a dataflow requires knowing the sparse pattern, which for arbitrary matrices (unlike DNN layers with stable statistics) is unknown until runtime — and even then, patterns shift within a matrix. Third, whatever mechanism adapts must not cost much: SpGEMM is memory-bound, and area belongs to the cache.

The idea

Don’t union the dataflows — parametrize one. A naive multi-dataflow machine would implement inner-, outer-, and row-based execution separately and switch — expensive, slow, and mostly redundant. Spada instead starts from the row-based scheme (already the best all-rounder, with consistent CSR formats and small psums) and applies the oldest trick in dense linear algebra: loop blocking. Blocking the row and reduction loops turns the per-element traversal into an α×β window over A’s compressed nonzeros. The endpoints of the shape space are the classic dataflows: α=1 recovers row-based behavior (many psums merged immediately — output reuse), β=1 recovers outer-product behavior (a column of A elements sharing fetched B rows — input reuse). The inner-product scheme is deliberately left out: its ineffectual index intersections are a pure liability on sparse data.

Make reconfiguration nearly free by making modes similar. Because every shape runs the same template — fetch window, fetch matching B rows, multiply, merge groups of lanes — switching shapes only regroups lanes: P-queue thresholds, the bitonic sorting network’s output level, and the flexible reduction tree change; fetchers and multipliers don’t. Reconfiguration happens between passes (a few stall cycles, hidden behind B prefetching), not per element.

Balance load where sparsity makes it uneven. Lanes process B rows of wildly different lengths. Each lane pair shares a sort array that each cycle picks the two globally smallest column indices from either lane’s stream — so a lane that finished early transparently helps its neighbor, keeping merge order intact. Empty A rows are skipped outright via the CSR offsets.

Let the matrix teach you its shape — cheaply. Analytical prediction of reuse would require scanning A’s index structure (as expensive as doing the work) and modeling cache interference. Spada instead exploits two empirical regularities: within a band of contiguous rows, index overlap between neighbors is stable, and band boundaries correlate with jumps in row length — detectable from the CSR offsets array alone. So the window adapter partitions rows into bands by row-length discontinuities, profiles the first rows of each large band under all four shapes, and applies the winner to the band’s remainder. Small bands use cost-descent: stop profiling the moment a shape performs worse than its predecessor.

band 1: rows 0–799 (short, uniform)
band 2: rows 800+ (long)
thresholds: T_abs=5, T_rel=2, T_band=128 rows

Why it works

• The classic dataflows differ, at core, in loop order → blocking the row-based nest exposes that difference as two block sizes → one hardware template spans the reuse spectrum, with the classic schemes as endpoints of a dial rather than separate machines.
• All shapes share the fetch-multiply-merge pipeline → mode switches are lane regroupings, not datapath changes → reconfiguration costs cycles, not area, so adapting per region is affordable.
• Sparse patterns are locally stable (bands) and band boundaries announce themselves in row lengths → profiling a band’s first rows predicts its remainder → the system learns the right shape from the matrix itself, with no offline analysis and no model of cache behavior.
• Load imbalance is the tax on parallel sparsity → pairwise sort arrays and task-aware P queues absorb row-length variance at O(n²) cost only within pairs → multipliers, the scarce resource, stay busy.
• Memory dominates SpGEMM → the shared cache plus RidxLRU keeps both B rows and psums useful across window shapes → adaptivity comes at GAMMA-equivalent area (6.32 vs 6.13 mm²) instead of SpArch’s 2×.

What it costs

• Profiling passes run suboptimal shapes. Each large band spends its first rows trying all four shapes; small bands may settle on a suboptimal shape via early-stopped cost descent.
• Reconfiguration stalls. Draining P queues between passes of different shapes costs cycles, partly hidden by prefetching.
• The window is small. Shapes are limited to the 8-lane budget (1×8…8×1); patterns that would reward much larger windows are out of reach, and sort-array balancing only spans lane pairs.
• Two FC layers still lose. On compressed fully-connected layers the gap to SpArch’s massive merger remains — adaptation inherits the row-based template’s limits rather than transcending them.
• Simulation + partial RTL. Performance comes from a cycle-accurate simulator with synthesized key components, not a full chip.

The evaluation

Twenty-seven workloads stress the diversity claim: 18 SuiteSparse matrices spanning densities from 2.5×10⁻⁷ to 0.49 and 4 to 5,162 nonzeros per row (squared, or multiplied by their transpose, following GAMMA), plus compressed ResNet-50, AlexNet, and pruned BERT layers. Spada is compared against re-implementations of the three fixed-dataflow flagships — SIGMA (inner product), SpArch (outer product), GAMMA (row-based) — with equal multiplier counts.