Sparse direct solver

Deep dive

How it works

Parts:

• Ordering — permutation P that limits fill (AMD, nested dissection).
• Elimination tree — pattern-only parent links; gives column counts, nnz(L) and flops before any arithmetic.
• Supernodes / fronts — columns sharing a pattern, grouped into dense blocks.
• Factor L (or L, U) — the stored result, reused for every right-hand side.
• Triangular solves — forward and back substitution on L.

A, pattern ──order──▶ PAPᵀ ──symbolic──▶ etree, nnz(L) ──numeric──▶ L ──solve──▶ x
 (once per pattern)                                (per values)           (per b)

1. Order. Choose a permutation P so that PAPᵀ has little fill: minimum degree/AMD (small to medium), nested dissection (METIS, large, parallel). See fill-reducing ordering. For LU, a column ordering (COLAMD) plus numerical row pivoting.
2. Symbolic analysis (pattern only). Elimination tree: parent(j) = first k > j with L[k][j] ≠ 0, computed in near-O(nnz(A)) with path compression. Row k of L is the union of tree paths from the nonzeros of A[:,k] up to k, so column counts, nnz(L) and flops are known before any arithmetic. Postorder the tree and merge columns with identical structure into supernodes.
3. Numeric factorization. Up-looking (code below, one row of L at a time: a sparse triangular solve), left-looking (supernodal, CHOLMOD) or multifrontal (a dense frontal matrix per tree node; children's update matrices are summed in by "extend-add"; MUMPS, UMFPACK). Supernodal and multifrontal codes spend nearly all time in dense GEMM/TRSM.
4. Solve. Forward L y = Pb, back Lᵀ z = y, x = Pᵀz. SPD needs no pivoting; symmetric indefinite uses LDLᵀ with 1×1/2×2 pivots; general LU uses threshold partial pivoting, which can add fill beyond the symbolic prediction (static pivoting plus refinement avoids it).

Independent subtrees of the tree factor in parallel; a node waits for its children.

Code

const NONE: usize = usize::MAX;

/// Upper triangle of a symmetric matrix in CSC: column k holds rows i <= k.
pub struct Csc { pub n: usize, pub p: Vec<usize>, pub i: Vec<usize>, pub x: Vec<f64> }

/// Build the upper-triangular CSC of P A P^T from (row, col, value) triplets of one triangle.
/// perm[new] = old.
pub fn permuted_upper(n: usize, trip: &[(usize, usize, f64)], perm: &[usize]) -> Csc {
    let mut pinv = vec![0; n];
    for (new, &old) in perm.iter().enumerate() { pinv[old] = new; }
    let mut p = vec![0; n + 1];
    let e: Vec<_> = trip.iter().map(|&(i, j, v)| {
        let (a, b) = (pinv[i], pinv[j]);
        let (r, c) = if a <= b { (a, b) } else { (b, a) };
        p[c + 1] += 1;
        (r, c, v)
    }).collect();
    for c in 0..n { p[c + 1] += p[c]; }
    let mut next = p.clone();
    let (mut i, mut x) = (vec![0; e.len()], vec![0.0; e.len()]);
    for (r, c, v) in e { i[next[c]] = r; x[next[c]] = v; next[c] += 1; }
    Csc { n, p, i, x }
}

/// Elimination tree: parent[j] = first column k > j with L[k][j] != 0.
pub fn etree(a: &Csc) -> Vec<usize> {
    let n = a.n;
    let (mut parent, mut anc) = (vec![NONE; n], vec![NONE; n]);
    for k in 0..n {
        for &r in &a.i[a.p[k]..a.p[k + 1]] {
            let mut i = r;
            while i != NONE && i < k { // walk up with path compression
                let next = anc[i];
                anc[i] = k;
                if next == NONE { parent[i] = k; }
                i = next;
            }
        }
    }
    parent
}

/// Nonzero pattern of row k of L (excluding the diagonal), in topological order.
fn ereach(a: &Csc, k: usize, parent: &[usize], mark: &mut [usize], s: &mut Vec<usize>) {
    s.clear();
    mark[k] = k;
    for &r in &a.i[a.p[k]..a.p[k + 1]] {
        let (mut i, start) = (r, s.len());
        if i >= k { continue; }
        while mark[i] != k { s.push(i); mark[i] = k; i = parent[i]; }
        s[start..].reverse();
    }
    s.reverse(); // later paths first, each path leaf to ancestor: descendants before ancestors
}

pub struct Chol { pub n: usize, pub lp: Vec<usize>, pub li: Vec<usize>, pub lx: Vec<f64> }

/// Symbolic pass + numeric up-looking Cholesky; diagonal is first in each column. None if not SPD.
pub fn cholesky(a: &Csc) -> Option<Chol> {
    let n = a.n;
    let parent = etree(a);
    let (mut mark, mut s) = (vec![NONE; n], Vec::new());
    let mut cnt = vec![1usize; n]; // symbolic: column counts of L
    for k in 0..n { ereach(a, k, &parent, &mut mark, &mut s); for &i in &s { cnt[i] += 1; } }
    let mut lp = vec![0; n + 1];
    for j in 0..n { lp[j + 1] = lp[j] + cnt[j]; }
    let (mut li, mut lx) = (vec![0; lp[n]], vec![0.0; lp[n]]);
    let mut fill: Vec<usize> = lp[..n].iter().map(|&p| p + 1).collect(); // next free slot per column
    let (mut x, mut mark) = (vec![0.0; n], vec![NONE; n]);
    for k in 0..n {
        ereach(a, k, &parent, &mut mark, &mut s);
        let mut d = 0.0;
        for p in a.p[k]..a.p[k + 1] { // scatter column k of A into the dense work vector
            if a.i[p] < k { x[a.i[p]] = a.x[p]; } else { d = a.x[p]; }
        }
        for &i in &s { // row k of L: sparse triangular solve L[0..k,0..k] y = A[0..k,k]
            let lki = x[i] / lx[lp[i]];
            x[i] = 0.0;
            for p in lp[i] + 1..fill[i] { x[li[p]] -= lx[p] * lki; }
            d -= lki * lki;
            li[fill[i]] = k; lx[fill[i]] = lki; fill[i] += 1;
        }
        if !(d > 0.0) { return None; }
        li[lp[k]] = k; lx[lp[k]] = d.sqrt();
    }
    Some(Chol { n, lp, li, lx })
}

impl Chol {
    /// Solve (P A P^T) x = b in the permuted space.
    pub fn solve(&self, b: &mut [f64]) {
        for j in 0..self.n {
            b[j] /= self.lx[self.lp[j]];
            for p in self.lp[j] + 1..self.lp[j + 1] { b[self.li[p]] -= self.lx[p] * b[j]; }
        }
        for j in (0..self.n).rev() {
            for p in self.lp[j] + 1..self.lp[j + 1] { b[j] -= self.lx[p] * b[self.li[p]]; }
            b[j] /= self.lx[self.lp[j]];
        }
    }
    pub fn nnz(&self) -> usize { self.lp[self.n] }
    /// Multiply-add count of the numeric factorization: sum of (column count)^2.
    pub fn flops(&self) -> usize { (0..self.n).map(|j| (self.lp[j + 1] - self.lp[j]).pow(2)).sum() }
}

Tested on a k×k 5-point Laplacian in natural and nested-dissection order (max error about 10⁻¹³; an indefinite 2×2 returns None).

Complexity and constants

• Measured with this code (scalar, one Apple-silicon core, release build, symbolic included) on the 5-point Laplacian, n = k²:

At k = 127 ordering alone cuts nnz(L) 4.5× and flops 7.9× (natural order is the banded case: nnz(L) ~ n^1.5, flops ~ n²).

• Asymptotics for model meshes under nested dissection: 2D nnz(L) O(n log n), flops O(n^1.5); 3D nnz(L) O(n^(4/3)), flops O(n²).
• Memory: 8 B per value + 4-8 B per row index per nonzero of L (supernodes store one index per block row, much less). The solution phase streams all of L once for 2·nnz(L) multiply-adds: memory-bandwidth bound, order of 0.1-0.2 flop/byte.
• Cost drivers: nnz(L) (bytes), Σ colcount² (flops), supernode size (closeness to GEMM peak), number of right-hand sides.

Efficiency

Cost: factorization of the 5-point Laplacian with this page's scalar code on one Apple-silicon core: 10.8 ms at n = 16,129 with nested dissection, 73 ms in natural order (table above). Supernodal and multifrontal codes are compute-bound on dense GEMM when fronts are large, otherwise memory-bound; the solve phase streams all of L once (2·nnz(L) multiply-adds) and is always memory-bound, order of 0.1-0.2 flop/byte.

Levers:

• Fill-reducing ordering — 4.5× less nnz(L), 7.9× fewer flops at n = 16k (measured above) — always, before factoring.
• Reuse ordering and symbolic analysis — skips them on every refactorization — when only values change (Newton, time stepping).
• Batch right-hand sides into an n×m block — L streamed once instead of m times (BLAS-3) — many solves per factor.
• Mixed precision, factor in f32 and refine in f64 — halves factor memory and traffic (see LU and Cholesky) — well-conditioned systems.
• Incomplete factorization or multigrid — O(nnz) memory — when the factor will not fit.

Naive → standard → state of the art:

• Natural order, simplicial code · nnz(L) ~ n^1.5 in 2D (table above)
• Nested dissection or AMD plus symbolic reuse · nnz(L) O(n log n) in 2D, one pattern analysis per matrix family
• Supernodal/multifrontal with threaded BLAS (CHOLMOD, MUMPS) · dense-block speed on big fronts

Implementation notes

• Build the permuted upper-triangular CSC of PAPᵀ once; store L in CSC with the diagonal first in each column.
• The symbolic pass allocates L exactly (no reallocation while factoring). Davis's cs_counts gets column counts in near-O(nnz(A)); the pass here costs O(nnz(L)).
• One dense work vector (zeroed as consumed) and one pattern stack; the pattern must be in topological order (descendants first) or the triangular solve reads unfinished entries.
• Threads follow the tree: subtrees in parallel, the few big top fronts with threaded BLAS.

Pitfalls

• Skipping the ordering (natural mesh numbering): several times more fill in 2D (4.5× above), far worse in 3D.
• General LU on an SPD matrix: about twice the flops and memory of Cholesky.
• Re-running ordering and symbolic analysis on every Newton step.
• Sizing memory from nnz(A): read nnz(L) from the symbolic phase and check it against RAM first.
• Static pivoting without checking the residual and refining.
• Forgetting the permutation: x is in permuted space until Pᵀ is applied.

Variants and alternatives

In production

SuiteSparse (CHOLMOD, UMFPACK; behind MATLAB and Julia sparse \; scikit-sparse wraps CHOLMOD), MUMPS, SuperLU, Intel oneMKL PARDISO, Eigen, faer; FEM codes and interior-point optimizers call these as black boxes.

Implementations

Sources

• Timothy A. Davis, Direct Methods for Sparse Linear Systems, SIAM 2006. book · code
• Joseph W. H. Liu, The role of elimination trees in sparse factorization, SIAM J. Matrix Anal. Appl. 11(1), 1990. paper
• Iain S. Duff, John K. Reid, The multifrontal solution of indefinite sparse symmetric linear equations, ACM TOMS 9(3), 1983. paper
• Alan George, Nested dissection of a regular finite element mesh, SIAM J. Numer. Anal. 10(2), 1973. paper
• Yanqing Chen, Timothy A. Davis, William W. Hager, Sivasankaran Rajamanickam, Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization and update/downdate, ACM TOMS 35(3), 2008. paper
• Patrick R. Amestoy, Iain S. Duff, Jean-Yves L'Excellent, Jacko Koster, A fully asynchronous multifrontal solver using distributed dynamic scheduling, SIAM J. Matrix Anal. Appl. 23(1), 2001. paper