A Systematic Survey of General Sparse Matrix-Matrix Multiplication

• SpGEMM is an expensive kernel in scientific computing and graph algorithms, and its optimization has potential to affect many applications.Source: §1 · p. 1
• The authors present the first survey, to their knowledge, covering SpGEMM developments and available algorithms, data structures, and libraries.It follows systematic literature review guidelines and covers applications, formats, formulations, challenges, architecture-oriented optimizations, and performance evaluation.Source: §1 · pp. 1–2
• The survey addresses application and formulation, research status, implementation performance, and challenges and future directions through four research questions.Source: §1 · p. 2
• The literature search iteratively queried multiple digital libraries and included 92 SpGEMM-related papers.The paper reports a rapid upward trend in articles over the prior three years.Source: §1 · p. 2
• A fully systematic taxonomy is difficult because topics such as load balance can be understood differently across computing contexts.The survey therefore selects frequently covered topics and presents paper-topic correlations in section tables or figures.Source: §1 · p. 2

2.1.1 Notation.

• The paper distinguishes matrices, vectors, and scalars typographically and provides a table of symbolic representations.Source: §2.1.1 · p. 3, Table 1 · p. 3

2.1.2 Sparse Matrix.

• Sparse matrices have no strict formal definition; one quantitative definition requires NNZ to be O(n).Source: §2.1.2 · p. 3

2.1.3 Compression Format.

• Sparse storage formats reduce redundant storage and useless calculations; the survey covers basic formats and later specialized variants.COO, CSR, CSC, ELL, and DIA are basic formats; later examples include DCSC for hypersparse matrices and BCSR for distributed multiplication.Source: §2.1.3 · pp. 3–4

2.2 Typical Applications

• The section introduces application contexts and explains how SpGEMM is formulated and used in them.Source: §2.2 · p. 4

2.2.1 Multi-source BFS.

• Multi-source BFS represents concurrent searches from multiple vertices as repeated products of the graph adjacency matrix with a source or frontier matrix.Concurrent BFSs can share computation in one kernel without synchronization costs compared with sequential runs.Source: §2.2.1 · pp. 4–5, Figure 2 · p. 5

2.2.2 Markov Clustering (MCL).

• Markov clustering alternates matrix-square expansion, threshold pruning, and inflation until graph points gradually cluster into groups.Source: §2.2.2 · p. 5

2.2.3 Algebraic Multigrid Solvers.

• In algebraic multigrid setup, two SpGEMMs construct the coarse-grid system through the Galerkin product; these products take more than 80% of construction time.The setup constructs interpolation and restriction operators; the restriction operator is the transpose of interpolation.Source: §2.2.3 · p. 5

2.2.4 Others.

• The survey also identifies SpGEMM as an important component in genome assembly, databases, graph algorithms, and molecular dynamics.Source: §2.2.4 · p. 6

3.1 Overview

• RbR is the most-used formulation, followed by OP and IP, while CbC is least explored.The four approaches are row-by-row, outer-product, inner-product, and column-by-column, respectively.Source: §3.1 · p. 6

3.2 Row-by-row

• RbR computes each output row by summing products of nonzeros in an A row with corresponding rows of B.It uses row-wise partitioning of both inputs.Source: §3.2 · p. 6

3.3 Inner-product

• IP forms each output entry as the inner product of an A row and a B column, matching only shared nonzero indices.It partitions A by rows and B by columns.Source: §3.3 · pp. 6–7

3.4 Outer-product

• OP sums outer products of corresponding A columns and B rows to form C.It partitions A by columns and B by rows.Source: §3.4 · p. 7

3.5 Column-by-column

• CbC computes each output column by summing products of nonzeros in a B column with corresponding columns of A.It partitions both matrices by columns.Source: §3.5 · p. 7

3.6 Discussion

• Under the stated uniform-distribution assumptions, OP has the highest data reuse and on-chip memory requirement, while IP has the lowest of both.RbR and CbC are equal on both measures in the comparison.Source: §3.6 · p. 7
• Intermediate-result volume is lowest for IP and highest for OP; RbR and CbC are equal between them.IP generates no intermediate results, while OP generates O(N × nnz_C); RbR and CbC have O(nnz × nnz_C / N^2) as stated for RbR.Source: §3.6 · p. 8
• RbR favors row-major storage and CbC column-major storage; IP alone requires index matching, while semirings allow custom multiplication and addition.IP and OP prefer opposite input storage layouts, and neither dictates the output layout.Source: §3.6 · p. 8

4.1 Overview

• The workflow predicts output size, allocates memory, partitions and balances work, performs numeric multiplication, and accumulates results.Multiplication generates intermediate results that accumulation reduces to the final output.Source: §4.1 · p. 8
• Because output sparsity is unknown, precise output-size prediction can be expensive, and NNZ largely determines CSR/CSC memory footprint.SpGEMM also complicates parallelization by involving three sparse matrices; irregular nonzero distributions make accumulator design challenging.Source: §4.1 · p. 8

4.2.1 Precise Prediction.

• Precise prediction computes each output row or column's exact NNZ in a symbolic phase before calculating values in the numeric phase.Kokkos Kernels, cuSPARSE, MKL, and RMerge are listed as representative implementations.Source: §4.2.1 · p. 9

4.2.2 Probabilistic Method.

• Probabilistic methods estimate output NNZ from a three-layer graph representation or input-row sampling.The cited graph-based method gives a 1 ± ε approximation in O(n/ε^2) time; a later variant has expected O(n) time for some ε.Source: §4.2.2 · pp. 9–10

4.2.3 Upper-Bound Prediction.

• Upper-bound prediction allocates for a maximum output NNZ, often by summing the NNZ of B rows selected by each nonzero in an A row.ESC is a representative; another approach sizes reusable per-thread hash tables from a maximum row count.Source: §4.2.3 · p. 10

4.2.4 Progressive Method.

• Progressive prediction starts with an allocation and grows it when needed; Matlab is a representative implementation.Matlab typically makes a new block 1.5 times the current space when more is required; hybrid approaches use upper bounds for short rows and progressive allocation for long rows.Source: §4.2.4 · p. 10

4.2.5 Discussion.

• Precise prediction saves memory and permits reuse of C's structure for repeated products, but traverses inputs twice; upper bounds can overallocate, while failed estimates or initial allocations require further allocation.The preferred method depends on the problems being addressed.Source: §4.2.5 · p. 10

4.3.1 Block Partition.

• The workcube represents input matrices and output as three views, with each voxel corresponding to a scalar multiplication contributing to an output cell.Block partitioning divides this workcube among computing units.Source: §4.3.1 · pp. 10–11
• In 1D partitioning, each unit receives a workcube layer; T, F, and S variants divide along planes parallel to different views.T assigns rows of A and all of B to produce corresponding rows of C; cited T-partition schemes also distribute rows by estimated or precise NNZ to balance work.Source: §4.3.1 · pp. 11–12, Figure 9 · p. 11
• When large row or column NNZ makes 1D communication substantial, 2D partitions divide two workcube dimensions; 3D partitions divide all three.2D methods can bottleneck as process counts rise, motivating 3D methods to reduce communication cost.Source: §4.3.1 · pp. 12–13, Figure 10 · p. 13
• SpSUMMA and SpCannon are 2D TF-partition examples: processors compute C blocks using row and column blocks of A and B.SpSUMMA broadcasts within processor rows and columns; SpCannon uses point-to-point exchanges.Source: §4.3.1 · p. 12
• Traditional block partitions may ignore per-block NNZ workload, and SpSUMMA assumes a uniform random distribution that may not hold.The survey also describes 3D approaches, including batching B blocks to fit available memory.Source: §4.3.2 · p. 13, §4.3.1 · p. 13

4.3.2 Graph Partition.

• Hypergraph partitioning models outer products and output nonzeros as vertices connected by nets; bipartite models represent the same dependencies with edges.Partitioning seeks roughly equal-sized parts while minimizing the cost of nets spanning parts.Source: §4.3.2 · p. 13, Figure 11 · p. 14
• Graph-based models have been used to reduce message volume and communication latency, and to partition work for MapReduce as well as 2D and 3D SpGEMM.The survey cites both hypergraph and bipartite graph approaches.Source: §4.3.2 · p. 14

4.3.3 Discussion.

• For irregular computation, graph partitioning can yield more balanced loads and lower communication cost than block partitioning, but at higher partitioning overhead.In 1D T-partition, the longest A row can dominate performance when B is regular; irregular B makes imbalance more serious.Source: §4.3.3 · p. 14

4.4 Result Accumulating

• The survey classifies accumulators by how they store intermediate results: dense, sparse, and hybrid.Figure 12 presents the classification.Source: §4.4 · p. 14, Figure 12 · p. 15

4.4.1 Dense Accumulator.

• Dense accumulators use three vectors to store values, mark occupied output positions, and record column indices.They offer random access but require large, often private per-thread storage, limiting scalability in highly concurrent settings.Source: §4.4.1 · p. 15

4.4.2 Sparse Accumulator.

• Sparse accumulators store row intermediates compactly, using list-based sorting or hash-based merging.List methods include radix, merge, and heap sorting; hash-based methods use column indices as keys, then shrink and sort the table to form the compressed output.Source: §4.4.2 · pp. 15–16

4.4.3 Hybrid Accumulator.

• Hybrid accumulators select dense or sparse strategies according to row, block, or tile characteristics.Examples choose among accumulators using per-row upper-bound non-zeros, expected block performance, or a preset tile-density threshold.Source: §4.4.3 · p. 16

4.4.4 Discussion.

• Dense accumulators suit short output rows, CPU caches, or rows with many non-zeros, while sparse accumulators suit wide outputs and GPUs.Dense storage can constrain concurrency; hash methods save memory but handle collisions, and radix sort is preferred for large arrays.Source: §4.4.4 · p. 17

5.1.1 Optimization for Memory Access.

• CPU techniques improve memory behavior through matrix partitioning, cache-aware allocation, propagation blocks, and coalesced DMA with local reuse of B.Patwary et al. target cache misses; Chen et al. [29] reserve a minimum partition of B in local data memory.Source: §5.1.1 · p. 17

5.1.2 Optimization for Load Balance.

• CPU load balancing uses dynamic partitions; Sunway's three-level partitioning divides A and B by computational load rather than rows.Patwary et al. report better results when partition count is 6–10 times the thread count.Source: §5.1.2 · pp. 17–18

5.1.3 Optimization for Data Structure.

• Patwary et al. convert B to blocked CSR only when over 30% of rows exceed the L2-cache-sized non-zero limit.The threshold avoids format-conversion overhead when the estimated proportion is lower.Source: §5.1.3 · p. 18

5.2.1 Optimization for Memory Access.

• GPU memory optimizations reduce global-memory traffic by staging or merging data in shared memory and choosing data placement based on input matrices.Reported methods also use shared memory for A or intermediate C results and hybrid pre-allocation for output storage.Source: §5.2.1 · p. 18, §5.2 · p. 18

5.2.2 Optimization for Load Balance.

• GPU load-balancing methods assign work across threads, warps, and blocks or group rows and tiles to handle irregular non-zero distributions.Examples include reorganizing column-row pairs, processing multiple rows per block, and selecting task assignment from input non-zero distributions [99].Source: §5.2.2 · pp. 18–19

5.2.3 Optimization for Data Structure.

• GPU data-structure optimizations parallelize merging and sorting, use shared-memory sorting when feasible, and exploit registers for sort, merge, and hash accumulators.When shared memory is insufficient, Dalton et al. launch global-memory sorting.Source: §5.2.3 · p. 19

5.3 Field-Programmable Gate Array (FPGA)

• FPGA designs exploit customizable data paths and parallelism, including tunable power-delay or energy-delay trade-offs and index comparisons of at least 8×8 per clock cycle.One reported framework is scalable and reconfigurable, but does not provide SpGEMM implementation details.Source: §5.3 · p. 19
• Small on-chip memories, lower operating frequency, and costly data transmission constrain FPGA SpGEMM use.The survey identifies these as reasons FPGA is not widely employed for SpGEMM.Source: §5.3 · p. 19

5.4 ASIC

• Six of seven ASIC-optimization papers were reported after 2020, while only one presented dedicated SpGEMM hardware.ASIC designs include outer-product and row-wise approaches, on-chip partial-result merging, and data-reuse structures; most reported work was evaluated with simulators.Source: §5.4 · p. 20

5.5 Heterogeneous Platform

• Heterogeneous methods divide work across CPU and GPU, stream CPU-prepared sparse data to FPGA, or place data across DRAM and persistent memory.CPU-GPU approaches use runtime or heuristics for work division; REAP performs symbolic analysis and packs an intermediate representation on CPU.Source: §5.5 · p. 20

5.6 Distributed Platform

• Distributed SpGEMM targets balanced node workloads and low communication, using 1D and 2D partitions, broadcasts, scheduling, and graph-based partitioning.SpCannon uses nearest-neighbor exchanges, whereas SpSUMMA uses row-wise and column-wise broadcasts; other approaches reduce communication or balance map and reduce loads.Source: §5.6 · pp. 21–22

5.7 Discussion

• Platform-specific priorities include cache use and load balance on CPUs, reduced global-memory access on GPUs, and balanced work with low communication on distributed systems.Heterogeneous platforms also seek work partitioning that uses both CPU and accelerator capability.Source: §5.7 · p. 22

6 PROGRAMMING MODEL

• The survey presents programming models as bridges between programmers’ logical views and physical execution, and lists popular models and frameworks.Table 4 summarizes different programming models.Source: §6 · p. 22, Table 4 · p. 22

6.1 Task-based programming model

• Task-based models recursively partition workloads into fork-join tasks, with implementations including Cilk, TBB, and Chunks and Tasks.Cilk uses work stealing; reported SpGEMM frameworks use tasks or blocks for load balance, while some customized libraries appear only in their reported work.Source: §6.1 · pp. 22–23

6.2 OpenMP and MPI

• OpenMP parallelizes threaded code, while MPI supports communication across distributed computers; distributed systems commonly combine MPI between nodes with OpenMP within nodes.OpenMP is also used for CPU-accelerator work division through heuristics.Source: §6.2 · p. 23

6.3 MapReduce

• A MapReduce SpGEMM implementation schedules map and reduce tasks statically to improve data locality and load balance.The reported implementation uses MR-MPI on distributed machines.Source: §6.3 · p. 23

6.4 CUDA and OpenCL

• CUDA and OpenCL support GPU programming with SIMT, but high performance requires architectural understanding and correct code remains challenging on newer GPUs.CUDA streams can overlap executions and data transfers; OpenCL is an open standard supporting devices from multiple vendors.Source: §6.4 · p. 23

7.1 Overview

• The study compares state-of-the-art SpGEMM implementations across CPU, GPU, and distributed-system platforms.Kokkos Kernels implementations are labeled KK-OpenMP on CPU and KK-CUDA on GPU.Source: §7.1 · p. 24

7.2 System Setup and Benchmark

• CPU and GPU tests use 1,880 square and 678 rectangular SuiteSparse matrices; distributed tests use 10 matrices with more than 5 million nonzeros.The benchmarks are A×A for square matrices and A×A^T for rectangular matrices, with single and double precision tested.Source: §7.2 · p. 24

7.3 Evaluation Results on CPU

• MKL-CSR and MKL-BSR outperform KK-OpenMP on most CPU matrices, averaging over 13x and 3x speedups, respectively; MKL-BSR wins on more than 900 matrices but fails on about 53%.MKL-CSR runs successfully on more than 90% of matrices and is fastest on about half; KK-OpenMP is comparable on large matrices. Single-to-double precision speedups average 1.05x for MKL-CSR, 1.17x for MKL-BSR, and 1.09x for KK-OpenMP.Source: Table 7 · p. 25, Figure 17 · p. 25, §7.3 · pp. 25–26

7.4 Evaluation Results on GPU

• On both GPUs, spECK and NSparse achieve the best performance for most matrices; TileSpGEMM wins on 10%–20% of successful runs, while KK-CUDA’s results differ sharply between P100 and V100.KK-CUDA averages 1.25x and 1.40x speedups over CUSP in single and double precision on P100, respectively, but is slower overall than the other algorithms on V100. CUSP and TileSpGEMM run successfully for about 72% and 68% of matrices on both GPUs.Source: Figure 18 · p. 26, Figure 19 · p. 26, §7.4 · pp. 26–27
• Single precision is only slightly faster than double precision: average speedups range from 1.00 to 1.23 on Tesla P100 and 1.02 to 1.20 on Tesla V100.The symbolic stage performs no floating-point calculation, so when it occupies a large share of runtime, changing precision has less effect.Source: Table 10 · p. 27, §7.4 · pp. 27–28

7.5 Evaluation Results on Distributed System

• CTF has higher parallel efficiency than PETSc in almost all tested process settings, but PETSc has better average performance; CTF’s average performance approaches or exceeds PETSc’s at higher process counts in double precision.CTF’s maximum speedup is higher beyond 32 processes for single precision and at every tested setting for double precision.Source: Figure 20 · p. 28, Table 11 · p. 28, §7.5 · p. 28

7.6 Evaluation Results of Preprocessing Overhead

• Format conversion dominates preprocessing: MKL-BSR conversion takes about 8 SpGEMM runs on average, while TileSpGEMM’s CSR2Tile conversion can cost tens of runs for some matrices.Tile2CSR costs about 5 SpGEMM runs on average; MKL export is negligible, and import takes less than one SpGEMM on average.Source: §7.6 · p. 29

Optimization based on machine learning

• Machine-learning optimization of SpGEMM is challenging because it must account for three matrices’ sparsity, matrix partitioning, and load balancing.The survey notes that machine-learning models have been used for parameter tuning and format or algorithm selection in SpMV.Source: §8 · p. 29

Size prediction

• More efficient size prediction is needed because allocation strategies can over-allocate memory or trigger costly reallocation, and the symbolic stage can exceed numerical-stage runtime for some matrices.On-device memory limits and data transfers make over- and re-allocation especially costly on accelerators.Source: §8 · p. 29

Heterogeneous architecture-oriented optimization

• Heterogeneous SpGEMM must balance work and minimize communication while matching denser blocks to accelerators and ultra-sparse blocks to CPUs.The paper contrasts CPUs’ strength in complicated logic with GPUs’ strength in dense computations.Source: §8 · p. 29

Application specific optimization

• Application-specific tuning can improve performance at the cost of extensibility; for classical AMG, P^T(AP) tends to require less per-row intermediate accumulation than (P^TA)P.The paper gives AmgX’s custom kernel for unsmoothed aggregation-based AMG as another example of exploiting application characteristics.Source: §8 · p. 30

SpGEMM in AI and HPC

• SpGEMM may suit convolution when both the input and model are sparse, and the survey expects future HPC devices to support it as AI techniques develop.The paper contrasts this prospect with dense conversion and tensor-core execution, which it says waste work and add overhead for sparse matrices.Source: §8 · p. 30

9 CONCLUSION

• Row-by-row is the most commonly used formulation, while CSR is the most frequently used storage format.The survey attributes row-by-row's use to high parallelism and low cache requirement, and CSR's to broad adoption across fields that avoids expensive format conversion.Source: §9 · p. 30
• On CPU, MKL performs best; on GPU, spECK and NSparse lead on many matrices, while TileSpGEMM excels on regular square matrices.Source: §9 · p. 30
• On distributed systems, PETSc leads with few processes, while CTF performs better as the process count increases.Source: §9 · p. 30
• Further work remains on supporting different formats, application- and architecture-oriented optimization, and efficient size prediction.Heterogeneous architectures may offer new optimization opportunities, and architecture-specific implementations are expected to continue.Source: §9 · p. 30