Sparse matrices (CSR)

Deep dive

How it works

In: COO triplets (row, col, value) and later a vector x · Out: the matrix as row_ptr/col/val, then y = A·x · State: the three arrays (read-only after the build) · Cost: build O(nnz + n); SpMV O(nnz), about 12 bytes read per nonzero for 2 flops

COO → CSR in O(nnz + n) with a counting sort:

1. Count entries per row into row_ptr[r+1].
2. Prefix-sum: row_ptr[i+1] += row_ptr[i].
3. Scatter each triplet to next[r]++ (a copy of row_ptr).
4. Sort each row's columns (short rows: insertion sort) and merge duplicates by summing, as FEM assembly needs.

SpMV: for each row i, y[i] = Σ val[k]·x[col[k]] for k in row_ptr[i]..row_ptr[i+1]. val and col stream sequentially; x is gathered, so cache hits on x depend on how well-banded the matrix is (reorder with reverse Cuthill-McKee to improve). Aᵀ·x with the same CSR scatters instead: y[col[k]] += val[k]·x[i], which is a race when parallel; store CSC (the transpose) for parallel transpose products.

A = [10  0  0  2]     row_ptr = [0, 2, 3, 5, 6]
    [ 0  3  0  0]     col     = [0, 3, 1, 0, 2, 3]
    [ 1  0  4  0]     val     = [10, 2, 3, 1, 4, 7]
    [ 0  0  0  7]

Code

pub struct Csr { pub nrows: usize, pub ncols: usize, pub row_ptr: Vec<usize>, pub col: Vec<u32>, pub val: Vec<f64> }

impl Csr {
    /// Build from COO triplets; duplicates are summed, columns sorted within each row.
    pub fn from_coo(nrows: usize, ncols: usize, t: &[(u32, u32, f64)]) -> Csr {
        let mut row_ptr = vec![0usize; nrows + 1];
        for &(r, _, _) in t { row_ptr[r as usize + 1] += 1; }
        for i in 0..nrows { row_ptr[i + 1] += row_ptr[i]; }
        let mut next = row_ptr.clone();
        let (mut col, mut val) = (vec![0u32; t.len()], vec![0.0; t.len()]);
        for &(r, c, v) in t {
            let k = next[r as usize];
            col[k] = c; val[k] = v; next[r as usize] += 1;
        }
        // sort each row by column and merge duplicates, compacting in place
        let mut out = 0;
        let mut new_ptr = vec![0usize; nrows + 1];
        for i in 0..nrows {
            let (s, e) = (row_ptr[i], row_ptr[i + 1]);
            let mut row: Vec<(u32, f64)> = (s..e).map(|k| (col[k], val[k])).collect();
            row.sort_unstable_by_key(|p| p.0);
            for (c, v) in row {
                if out > new_ptr[i] && col[out - 1] == c { val[out - 1] += v; } else { col[out] = c; val[out] = v; out += 1; }
            }
            new_ptr[i + 1] = out;
        }
        col.truncate(out); val.truncate(out);
        Csr { nrows, ncols, row_ptr: new_ptr, col, val }
    }

    pub fn nnz(&self) -> usize { self.val.len() }

    /// y = A x
    pub fn spmv(&self, x: &[f64], y: &mut [f64]) {
        assert_eq!(x.len(), self.ncols);
        for i in 0..self.nrows {
            let (s, e) = (self.row_ptr[i], self.row_ptr[i + 1]);
            let mut acc = 0.0;
            for (v, &c) in self.val[s..e].iter().zip(&self.col[s..e]) { acc += v * x[c as usize]; }
            y[i] = acc;
        }
    }

    /// y = A^T x (scatter; sequential only)
    pub fn spmv_t(&self, x: &[f64], y: &mut [f64]) {
        y.iter_mut().for_each(|v| *v = 0.0);
        for i in 0..self.nrows {
            for k in self.row_ptr[i]..self.row_ptr[i + 1] { y[self.col[k] as usize] += self.val[k] * x[i]; }
        }
    }

    /// Explicit transpose (= CSC of A), O(nnz + n); columns come out sorted.
    pub fn transpose(&self) -> Csr {
        let mut ptr = vec![0usize; self.ncols + 1];
        for &c in &self.col { ptr[c as usize + 1] += 1; }
        for j in 0..self.ncols { ptr[j + 1] += ptr[j]; }
        let mut next = ptr.clone();
        let (mut col, mut val) = (vec![0u32; self.nnz()], vec![0.0; self.nnz()]);
        for i in 0..self.nrows {
            for k in self.row_ptr[i]..self.row_ptr[i + 1] {
                let j = self.col[k] as usize;
                col[next[j]] = i as u32; val[next[j]] = self.val[k]; next[j] += 1;
            }
        }
        Csr { nrows: self.ncols, ncols: self.nrows, row_ptr: ptr, col, val }
    }

    /// Parallel y = A x: rows split so each thread gets ~nnz/threads nonzeros.
    pub fn spmv_par(&self, x: &[f64], y: &mut [f64], threads: usize) {
        let nnz = self.nnz();
        let mut bounds = vec![0usize];
        for t in 1..threads {
            let target = nnz * t / threads;
            bounds.push(self.row_ptr.partition_point(|&p| p < target).min(self.nrows));
        }
        bounds.push(self.nrows);
        std::thread::scope(|s| {
            let mut rest = &mut y[..];
            for w in bounds.windows(2) {
                let (r0, r1) = (w[0], w[1].max(w[0]));
                let (mine, tail) = std::mem::take(&mut rest).split_at_mut(r1 - r0);
                rest = tail;
                s.spawn(move || {
                    for (yi, i) in mine.iter_mut().zip(r0..r1) {
                        let (a, b) = (self.row_ptr[i], self.row_ptr[i + 1]);
                        *yi = self.val[a..b].iter().zip(&self.col[a..b]).map(|(v, &c)| v * x[c as usize]).sum();
                    }
                });
            }
        });
    }
}

Complexity and constants

• Memory: f64 values + u32 columns = 12 bytes per nonzero, plus 8 bytes per row pointer. A 5-point Laplacian on 1000 × 1000 has n = 10⁶ rows and about 5·10⁶ nonzeros: 68 MB.
• SpMV traffic per call: 12·nnz + 8·n (pointers) + 8·n (y) + x reads (≥ 8·n): about 84 MB for that matrix.
• Measured with the code above on that Laplacian (Apple M4 Pro, -O): COO → CSR build 46 ms; serial SpMV 2.1 ms (about 40 GB/s effective, 4.8 GFLOP/s); spmv_par 1.27 ms on 2 threads, 0.71 ms on 4, 0.68 ms on 8 (about 120 GB/s): it saturates memory bandwidth at a few cores, more threads do not help.
• Conversion COO → CSR: O(nnz) plus per-row sorting; the triplet array is 16 bytes per entry, so assembly peaks at about 2.3× the final CSR size.

Efficiency

Cost: one SpMV on the 5-point Laplacian with 10⁶ rows and about 5·10⁶ nonzeros moves about 84 MB. Measured with this page's code on an Apple M4 Pro: 2.1 ms serial (about 40 GB/s effective, 4.8 GFLOP/s), 0.68 ms on 8 threads (about 120 GB/s). Bound: memory bandwidth, about 2 flops per 12 bytes.

Levers:

• Threads — 3.1× at 4 threads, no gain beyond (measured above) — saturates bandwidth within a few cores.
• u32 instead of u64 column indices — 12 instead of 16 bytes per nonzero, a quarter less traffic — fewer than 2³² columns.
• Equal-nnz row split — removes the one-thread-does-most case — power-law rows.
• Reverse Cuthill-McKee reordering — better cache hits on x gathers — unstructured meshes.
• BSR or matrix-free — one index per block, or none — multi-component PDEs, stencils.

Naive → standard → state of the art:

• COO triplets, 16 bytes per entry, scatter per entry · slower, about 1.3× the bytes of CSR (16 vs 12)
• CSR, u32 columns, nnz-balanced threads · bandwidth-bound, 12 B/nnz
• Blocked or SELL-C-σ layouts, matrix-free stencils · fewer index bytes, SIMD/GPU friendly

Implementation notes

• Index width: u32 columns and u32 (or u64 when nnz ≥ 2³²) row pointers; the code above uses usize pointers for simplicity.
• Block CSR (BSR): for elasticity or multi-component PDEs, store dense 3×3 or 4×4 blocks; one column index per block cuts index traffic and enables SIMD.
• SELL-C-σ / ELLPACK: pad rows to equal length within chunks of C rows for SIMD and GPUs; sorting rows by length within windows of σ limits padding.
• Symmetric storage (upper triangle only) halves memory but makes SpMV a scatter; usually not worth it for parallel code.
• Reordering (RCM, METIS nested dissection) improves locality of x gathers and is mandatory for sparse direct solvers.
• Matrix-free: for stencils and high-order FEM, compute A·x on the fly; no index traffic at all.
• Compressed indices: delta-encode column indices within rows (bit packing) when bandwidth is the bottleneck and decode is cheap; for 0/1 matrices, bitmaps (roaring) or pattern-only CSR.

Pitfalls

• Splitting rows evenly across threads when a few rows are dense (power-law graphs): one thread does most of the work.
• Parallel Aᵀx on CSR by scattering without atomics: data races; build the transpose once.
• Assuming sorted, duplicate-free columns from external data: validate or canonicalize, many kernels (merging, triangular solves) depend on it.
• Storing explicit zeros from assembly (cancellations) and then counting them as structural nonzeros: harmless for SpMV, wasteful for direct solvers; prune after assembly if they are structural.
• Using u32 for row pointers with more than 4 billion nonzeros: silent overflow.

Variants and alternatives

In production

SciPy scipy.sparse.csr_matrix, PETSc MATAIJ (CSR) and MATBAIJ (BSR), Intel MKL sparse BLAS, cuSPARSE (CSR, COO, BSR, SELL), Eigen SparseMatrix (CSC by default), SuiteSparse (CSC), Rust sprs and faer, PyTorch torch.sparse_csr_tensor.

Implementations

Sources

• Yousef Saad, Iterative Methods for Sparse Linear Systems, 2nd ed., SIAM, 2003, chapter 3. book
• Timothy A. Davis, Direct Methods for Sparse Linear Systems, SIAM, 2006. book
• Moritz Kreutzer, Georg Hager, Gerhard Wellein, Holger Fehske, Alan R. Bishop, A Unified Sparse Matrix Data Format for Efficient General Sparse Matrix-Vector Multiplication on Modern Processors with Wide SIMD Units, SIAM Journal on Scientific Computing 36(5), 2014. paper
• Samuel Williams, Leonid Oliker, Richard Vuduc, John Shalf, Katherine Yelick, James Demmel, Optimization of Sparse Matrix-Vector Multiplication on Emerging Multicore Platforms, SC 2007. paper