Source-linked AI summary
Parallel Sparse Matrix-Matrix Multiplication and Indexing: Implementation and Experiments
Aydin Buluc, John Gilbert
TL;DR
The paper addresses the lack of scalable, general sparse-matrix indexing in distributed memory and the need for flexible SpGEMM implementations. It uses two-dimensional block distributions with serial hypersparse kernels to build SpGEMM, SpRef, and SpAsgn. The implementations scale to thousands of processors, achieve performance predicted by the analyses, and expose communication and hierarchical-parallelism limits at high concurrency.
Problem
General sparse indexing with arbitrary row and column vectors lacked a parallel algorithm, while prior 1D SpGEMM algorithms were not scalable to thousands of processors.
Method
The paper develops flexible distributed-memory SpGEMM using two-dimensional block distributions and serial hypersparse kernels, then uses it to implement SpRef and SpAsgn.
Results
The implementations scale to thousands of processors, and experiments confirm performance predicted by the analyses; a one-sided implementation can be up to 46% faster for a specified R-MAT workload.
Takeaways & Limitations
SpGEMM serves as a building block for sparse indexing, assignment, and graph operations including contraction, subgraph extraction, breadth-first search, and batch updates.
Takeaways & Limitations
Beyond 2000 processors, more than 75% of runtime is spent in inter-node communication, motivating communication-avoiding and hierarchical algorithms.
Abstract
from arXiv · showhide
Generalized sparse matrix-matrix multiplication (or SpGEMM) is a key primitive for many high performance graph algorithms as well as for some linear solvers, such as algebraic multigrid. Here we show that SpGEMM also yields efficient algorithms for general sparse-matrix indexing in distributed memory, provided that the underlying SpGEMM implementation is sufficiently flexible and scalable. We demonstrate that our parallel SpGEMM methods, which use two-dimensional block data distributions with serial hypersparse kernels, are indeed highly flexible, scalable, and memory-efficient in the general case. This algorithm is the first to yield increasing speedup on an unbounded number of processors; our experiments show scaling up to thousands of processors in a variety of test scenarios.
1. Introduction.
The paper develops scalable parallel implementations of SpGEMM and generalized sparse indexing, using SpGEMM to support SpRef and extending the approach to assignment and graph operations.
- Sparse kernels: SpGEMM computes sparse matrix products over a general semiring, while SpRef extracts A(I, J) for arbitrary row and column index vectors.The SpRef approach uses SpGEMM as its key subroutine, regularizing computation and data access patterns.
- Main contributions: SpGEMM enables a simple and efficient implementation of SpRef, alongside a distributed-memory implementation supporting broader applications and processor layouts.The paper also reports extensive performance experiments with SpGEMM and SpRef.
- Main contributions: SpAsgn assigns a sparse matrix B to a submatrix A(I, J), supporting streaming batch updates to a graph.The operation is presented as an additional sparse-matrix primitive with reported parallel performance.
- Paper organization: The paper presents parallel algorithms and theoretical performance analyses for SpGEMM and SpRef, followed by large-scale performance evaluations and implementation studies.The evaluation includes comparisons with similar Trilinos primitives and examines how implementation decisions affect performance.
- Notation: The notation defines sparse matrices over a semiring, nonzero counts, Matlab-style indexing, one-based indices, array lengths, and nonzero-operation counts for products.The notation distinguishes flops((AB)·C) from flops(A·(BC)) because multiplication order changes the required work.
2. Notation.
The analysis measures message latency and inverse bandwidth as operation-time costs and uses Θ notation for asymptotically tight bounds.
- Communication costs: Message latency α and inverse bandwidth β are expressed as time for a floating-point operation, including cache-miss and memory-indirection costs.These parameters are used in parallel running-time analyses.
- Asymptotic notation: The notation f(x) = Θ(g(x)) means that f is bounded asymptotically above and below by g.
3. Sparse matrix-matrix multiplication.
SpGEMM computes sparse matrix products over general semirings using distributed two-dimensional decompositions, hypersparse storage, and Sparse SUMMA. These choices address scalability and memory costs while supporting flexible matrix shapes, processor grids, and blocking parameters.
- SpGEMM supports graph algorithms, scientific computing applications, and general computing tasks such as parsing and intersection searching.
- O(flops + nnz + n) is the classical serial column-wise SpGEMM cost when flops ≥ max{nnz, n}, using CSC storage.
- 2D decomposition assigns submatrices to a pr × pc processor grid, with processor P(i, j) storing an (m/pr) × (n/pc) block.
- 1D SpGEMM algorithms do not scale to thousands of processors because auxiliary data structures cannot be loaded and unloaded quickly enough to amortize their costs.
- Hypersparse local submatrices make CSC storage asymptotically wasteful: total memory becomes O(n√p + nnz), versus O(n + nnz) for the whole matrix on one processor.
- 3.2. DCSC Data Structure: DCSC stores only nonempty columns and their indices, providing O(nnz) storage and enabling fast column indexing through an AUX array.
- 3.3. Sparse SUMMA algorithm: Sparse SUMMA broadcasts operand blocks along processor rows and columns, then updates each local result with HyperSparseGEMM.
- 3.3. Sparse SUMMA algorithm: Sparse SUMMA handles rectangular matrices, rectangular processor grids, and adjustable blocking parameters through a generalized 2D algorithm.
4. Sparse matrix indexing and subgraph selection.
SpRef generalizes sparse indexing by extracting ordered row-column submatrices through SpGEMM, while SpAsgn assigns sparse submatrices using related products. Their parallel implementations exploit hypersparsity and 2D communication patterns to support scalable distributed execution.
- SpRef: SpRef extracts the ordered submatrix A(I, J), and selecting identical row and column vectors yields an induced subgraph.It also supports random row and column permutations for load balancing in parallel matrix computations.
- SpRef: A general parallel algorithm for arbitrary row and column index vectors was previously absent; the proposed method uses parallel SpGEMM.Simple row, column, and element indexing were commonly handled by special-purpose routines.
- SpAsgn: SpAsgn implements A(I, J) = B through sparse products that embed B and remove the prior A(I, J) region.When additive inverses are unavailable, subtraction is replaced by generalized sparse elementwise multiplication.
- SpRef: O(nnz(A)) sequential complexity makes SpRef asymptotically optimal for arbitrary I and J.The two products satisfy flops(R·A) ≤ nnz(A) and flops((R·A)·Q) ≤ nnz(A), while writing a general permutation requires Ω(nnz(A)).
- SpAsgn: O(nnz(A) + nnz(B) + len(I) + len(J)) is the sequential SpAsgn complexity using an optimal Θ(flops) SpGEMM routine.The bound follows from the input nonzeros, index-vector lengths, and sparse intermediate products.
- Parallel SpRef: Θ(√p) is the asymptotic speedup limit for parallel SpRef because SpGEMM costs dominate its execution.The parallel design respects hypersparsity and seeks communication along one processor dimension; vector distribution across all processors preserves asymptotic costs.
5. Experimental Results.
Experiments show that the sparse kernels scale effectively across diverse matrix operations and processor counts, while communication becomes the dominant cost at high concurrency. The implementation also outperforms Trilinos substantially and remains scalable and memory-efficient on large problems.
- Experimental platform: The experiments used NERSC’s Franklin system, a 9660-node Cray XT4 with quad-core 2.3 GHz AMD Opteron processors and a 3D torus interconnect.The algorithms also performed similarly on a fat-tree topology in earlier Ranger experiments.
- Indexing and subgraph operations: Linear speedup persisted to about 64 processors for random permutation and then followed √p scaling; extracting 10 subgraphs scaled slightly worse because smaller extractions reduce available parallelism.The permutation used a scale 22 R-MAT matrix with approximately 32 million nonzeros.
- Sparse matrix multiplication: SpGEMM exceeded 2 billion useful double-precision flops per second on 8100 cores when multiplying scale 24 R-MAT matrices.Useful-flop rates depend on matrix structure and sparsity.
- Communication and blocking: At 8100 cores, communication consumed 83% of SpGEMM time, despite computation decreasing by over 90x from 36 cores; doubling blocking stages slowed computation by only 1–6%.Communication includes synchronization and load-imbalance overheads; its cost was unaffected by the tested blocking sizes.
- Restriction operations: Restriction operations scaled to 4096 processors, with order changes affecting performance by less than 5%; full order-8 restriction achieved 110× speedup on 1024 processors and 163× on 4096.Different parenthesizations of the triple product did not affect performance.
- Comparison with Trilinos: Sparse SUMMA was 66× faster than Trilinos at 256 processors and used less memory, while restriction multiplication was 65X faster at 121 processors and scaled to 4096 cores.Trilinos ran out of memory beyond scale 21 or 256 processors, whereas EpetraExt slowed beyond 16 cores.
6. Conclusions and Future Work.
The paper presents Sparse SUMMA as a flexible distributed-memory SpGEMM algorithm that supports scalable sparse indexing, assignment, and graph operations. Experiments validate the predicted performance, while future work targets communication, matrix-chain, asynchronous, and hierarchical extensions.
- Conclusions: Sparse SUMMA scales to thousands of processors in distributed memory and serves as a building block for scalable SpRef and SpAsgn routines.These routines support graph coarsening by edge contraction, subgraph extraction, multi-source breadth-first search, and batch graph updates.
- Conclusions: Using SpGEMM as a building block enabled the most general analysis of SpRef, and experiments confirmed performance predicted by the analyses.
- Future Work: Sparse matrix chain products remain a future extension that could eliminate temporary intermediates and enable structure prediction and sparsity-based multiplication ordering.Triple products arise in algebraic multigrid coarsening, sparse matrix indexing, and parallel graph contraction.
- Future Work: More than 75% of execution time is spent in inter-node communication after 2000 processors, motivating communication-avoiding sparse matrix algorithms.
- Future Work: A truly one-sided implementation can run up to 46% faster than synchronous alternatives on scale-20 R-MAT matrices using 4000 processors.The experiments suggest synchronous SpGEMM algorithms have considerably higher load imbalance than asynchronous ones.
- Future Work: Increasing cores per node intensify network-interface contention, making hierarchical SpGEMM and SpRef parallelism an important future direction.The paper identifies flat MPI parallelism as unscalable without exploiting faster on-chip networks.