Parallel Sparse Matrix Computation is the high-performance computing discipline focused on efficient parallel algorithms for sparse matrices — matrices where the vast majority of elements are zero (>95% for typical scientific problems) — where specialized storage formats (CSR, CSC, COO, ELL), sparse matrix-vector multiplication (SpMV), and sparse direct/iterative solvers are the computational workhorses of scientific simulation, graph analytics, and machine learning, and where the irregular memory access patterns of sparse data make efficient parallelization fundamentally harder than dense linear algebra.

Why Sparse Matrices Are Hard to Parallelize

Dense matrix operations (GEMM) have regular, predictable memory access patterns that achieve >90% of peak FLOPS. Sparse matrices have indexed, indirect access patterns — for CSR format, computing row i requires loading column indices from col_idx[row_ptr[i]:row_ptr[i+1]] and then gathering values from the input vector at those indices. The indirect access causes random memory reads with near-zero cache hit rate on large problems.

Storage Formats

Parallel SpMV (Sparse Matrix-Vector Multiply)

SpMV (y = A·x) is the dominant kernel in iterative solvers (CG, GMRES, BiCGSTAB). Parallelization approaches:

• Row-per-thread (CSR): Each thread computes one row's dot product. Load imbalance when row lengths vary (power-law graphs).
• Warp-per-row (GPU): A full warp cooperatively computes one row using shuffle-based reduction. Better load balance for medium-length rows.
• Merge-based (CSR-Adaptive): Balances work by evenly distributing NON-ZEROS across threads (not rows). Each thread processes an equal share of the values array. Requires binary search to determine row boundaries.
• Segmented Reduction: Treat SpMV as a segmented reduction over the values array, where segments correspond to rows. GPU-friendly with balanced work distribution.

Sparse Solvers

• Iterative: Krylov methods (CG, GMRES) repeatedly apply SpMV + preconditioning. Parallel SpMV + parallel preconditioner (ILU, AMG) = parallel solver. Communication: one all-reduce per iteration for global dot products.
• Direct: Sparse LU/Cholesky factorization (SuperLU, CHOLMOD, MUMPS). Fill-in creates new nonzeros during factorization. Supernodal methods group dense subblocks for BLAS-3 efficiency. Parallel scalability limited by the elimination tree structure.

Parallel Sparse Matrix Computation is where the elegance of parallel algorithms meets the harsh reality of irregular memory access — requiring creative data structures and load-balancing techniques to extract parallelism from the inherently unstructured access patterns of sparse data.