Sparse Matrix Computation is the parallel computing discipline focused on efficient storage and computation with matrices where 90-99.9% of elements are zero — using compressed storage formats (CSR, CSC, COO, ELL) and specialized algorithms that perform operations proportional to the number of nonzeros (nnz) rather than the full matrix dimensions, critical for scientific computing, graph analytics, recommendation systems, and any domain where the underlying data is naturally sparse.

Why Sparse Matrices Are Everywhere

A finite element mesh with 10 million nodes produces a 10M×10M matrix (10¹⁴ elements = 800 TB at FP64). But each node connects to only ~20 neighbors, so only 200M entries are nonzero (1.6 GB). Storing and computing with the full dense matrix is impossible; sparse formats and algorithms are mandatory.

Storage Formats

• CSR (Compressed Sparse Row): Three arrays — values[] (nonzero values), col_idx[] (column index of each nonzero), row_ptr[] (index into values[] where each row starts). Row_ptr has N+1 entries; values and col_idx have nnz entries. The default format for sparse linear algebra. Row-oriented: efficient for row-based operations (SpMV).

• CSC (Compressed Sparse Column): Transpose of CSR — column-oriented. Efficient for column-based access (sparse triangular solve, some factorization algorithms).

• COO (Coordinate): Three arrays — row[], col[], value[] — one triple per nonzero. Simplest format, easy to construct. No implicit ordering. Used as an intermediate format during matrix assembly.

• ELL (ELLPACK): Each row is padded to the same length (max nonzeros per row). Stored as two dense 2D arrays (value[N][K], col[N][K]) where K = max nnz per row. GPU-friendly due to regular access patterns but wasteful for power-law degree distributions.

• Hybrid (HYB): ELL for the regular portion + COO for overflow rows with many nonzeros. Balances GPU efficiency with storage efficiency for irregular matrices.

Sparse Matrix-Vector Multiply (SpMV)

The dominant sparse operation: y = A×x. Each row i computes a dot product of its nonzero entries with corresponding x elements. In parallel, each thread (or warp) handles one or more rows:

• CSR SpMV: Thread i iterates from row_ptr[i] to row_ptr[i+1], accumulating value[j] * x[col_idx[j]]. Performance is memory-bound: arithmetic intensity is 2 FLOP / (12-16 bytes loaded) = 0.125-0.167 FLOP/byte — deep in the memory-bound region of the roofline.
• GPU Challenge: Short rows (few nonzeros) underutilize warps. Long rows (many nonzeros) overload individual threads. Solutions: CSR-Vector (one warp per row with warp-level reduction), merge-based SpMV (load-balanced distribution of nonzeros across threads).

Sparse Linear Solvers

• Iterative Solvers: Conjugate Gradient (CG), GMRES, BiCGSTAB — dominated by SpMV and vector operations. Parallelism is straightforward (SpMV is embarrassingly parallel by rows) but convergence depends on preconditioners.
• Direct Solvers: Sparse LU/Cholesky factorization. Fill-in (new nonzeros created during factorization) must be managed. Graph-based reordering (METIS, AMD) minimizes fill-in and maximizes parallelism.

Sparse Matrix Computation is the computational backbone of scientific and data-driven applications — where the structure of the real world (physical connections, social links, molecular bonds) naturally produces sparse data that requires specialized storage and algorithms to process at scale.