Incomplete factorization preconditioners

Deep dive

How it works

Parts:

• Pattern rule — which entries of L, U may be nonzero.
• Factorization — row-by-row elimination in CSR that drops entries outside the rule.
• Ordering — permutation applied first (fill-reducing ordering).
• Triangular solves — forward and back substitution every iteration.
• Outer Krylov loop — CG with M = LLᵀ (SPD), GMRES otherwise.

1. Reorder A symmetrically (RCM, minimum degree).
2. For each row i and allowed entry (i, j), j ≤ i: subtract the products over shared columns of rows i and j, divide by l_jj.
3. If a pivot is ≤ 0 (or tiny for ILU), refactor A + α·diag(A) with larger α (Manteuffel shift).
4. Each iteration: z = U⁻¹(L⁻¹ r), plus one SpMV with A.

Pattern rules:

• IC(0)/ILU(0) — A's pattern; memory equals A.
• ILU(k) — keep fill of level ≤ k, levels counted along fill paths, computed symbolically once.
• ILUT(p, τ) (Saad 1994) — drop entries below τ·‖row‖, keep the p largest in each of L and U per row.
• MILU — add dropped entries to the diagonal so row sums match A; faster on Poisson, more breakdown-prone.

A -> reorder -> factor with dropping -> L,U -> loop: r -> L⁻¹ -> U⁻¹ -> z ; q = A p

Code

IC(0) and PCG with M = LLᵀ against plain CG on 5-point Poisson. Compiled and run.

use std::cmp::Ordering::*;
pub struct Csr { pub n: usize, pub row_ptr: Vec<usize>, pub col: Vec<u32>, pub val: Vec<f64> }

impl Csr {
    /// Build from per-row (col, value) lists with ascending columns.
    pub fn from_rows(rows: Vec<Vec<(usize, f64)>>) -> Csr {
        let (mut row_ptr, mut col, mut val) = (vec![0], vec![], vec![]);
        for r in &rows {
            for &(c, v) in r { col.push(c as u32); val.push(v); }
            row_ptr.push(col.len());
        }
        Csr { n: rows.len(), row_ptr, col, val }
    }
    pub fn spmv(&self, x: &[f64], y: &mut [f64]) {
        for i in 0..self.n { y[i] = (self.row_ptr[i]..self.row_ptr[i + 1]).map(|k| self.val[k] * x[self.col[k] as usize]).sum(); }
    }
}
/// 5-point Poisson matrix on a k x k grid (4 on the diagonal, -1 per neighbour).
pub fn poisson2d(k: usize) -> Csr {
    Csr::from_rows((0..k * k).map(|id| {
        let (i, j) = (id % k, id / k);
        let nb = [(j > 0, id.wrapping_sub(k), -1.0), (i > 0, id.wrapping_sub(1), -1.0), (true, id, 4.0), (i + 1 < k, id + 1, -1.0), (j + 1 < k, id + k, -1.0)];
        nb.into_iter().filter(|t| t.0).map(|t| (t.1, t.2)).collect()
    }).collect())
}
/// IC(0): A ~ L L^T with L on the pattern of tril(A); None on a non-positive pivot.
pub fn ic0(a: &Csr) -> Option<Csr> {
    let mut l = Csr { n: a.n, row_ptr: vec![0], col: vec![], val: vec![] };
    for i in 0..a.n {
        for k in a.row_ptr[i]..a.row_ptr[i + 1] { if a.col[k] as usize <= i { l.col.push(a.col[k]); l.val.push(a.val[k]); } }
        l.row_ptr.push(l.col.len());
    }
    for i in 0..l.n {
        let (si, ei) = (l.row_ptr[i], l.row_ptr[i + 1]);
        for p in si..ei {
            let j = l.col[p] as usize;
            // l_ij = (a_ij - sum_{m<j} l_im l_jm) / l_jj; the sum runs over columns present in both rows
            let (mut s, mut x, mut y) = (l.val[p], si, l.row_ptr[j]);
            let ey = l.row_ptr[j + 1] - 1; // last entry of row j is its diagonal
            while x < p && y < ey {
                match l.col[x].cmp(&l.col[y]) {
                    Less => x += 1,
                    Greater => y += 1,
                    Equal => { s -= l.val[x] * l.val[y]; x += 1; y += 1; }
                }
            }
            if j == i && s <= 0.0 { return None; } // breakdown: shift the diagonal and retry
            l.val[p] = if j == i { s.sqrt() } else { s / l.val[ey] };
        }
    }
    Some(l)
}
/// z = (L L^T)^-1 r: forward solve L y = r, then back solve L^T z = y by scattering rows of L.
pub fn solve_llt(l: &Csr, r: &[f64], z: &mut [f64]) {
    for i in 0..l.n {
        let (s, e) = (l.row_ptr[i], l.row_ptr[i + 1]);
        let acc: f64 = (s..e - 1).map(|k| l.val[k] * z[l.col[k] as usize]).sum();
        z[i] = (r[i] - acc) / l.val[e - 1];
    }
    for i in (0..l.n).rev() {
        let (s, e) = (l.row_ptr[i], l.row_ptr[i + 1]);
        z[i] /= l.val[e - 1];
        for k in s..e - 1 { z[l.col[k] as usize] -= l.val[k] * z[i]; }
    }
}
fn dot(a: &[f64], b: &[f64]) -> f64 { a.iter().zip(b).map(|(x, y)| x * y).sum() }

/// CG from x = 0; `l = None` is plain CG, `Some(L)` preconditions with M = L L^T.
pub fn pcg(a: &Csr, b: &[f64], l: Option<&Csr>, rtol: f64, max_it: usize) -> (Vec<f64>, usize) {
    let n = b.len();
    let (mut x, mut r, mut z, mut q) = (vec![0.0; n], b.to_vec(), vec![0.0; n], vec![0.0; n]);
    let apply = |r: &[f64], z: &mut [f64]| match l { Some(l) => solve_llt(l, r, z), None => z.copy_from_slice(r) };
    apply(&r, &mut z);
    let (mut p, mut rz, tol) = (z.clone(), dot(&r, &z), rtol * dot(b, b).sqrt());
    for it in 1..=max_it {
        a.spmv(&p, &mut q);
        let alpha = rz / dot(&p, &q);
        for i in 0..n { x[i] += alpha * p[i]; r[i] -= alpha * q[i]; }
        if dot(&r, &r).sqrt() <= tol { return (x, it); }
        apply(&r, &mut z);
        let rz_new = dot(&r, &z);
        let beta = rz_new / rz;
        rz = rz_new;
        for i in 0..n { p[i] = z[i] + beta * p[i]; }
    }
    (x, max_it)
}

fn rel_res(a: &Csr, x: &[f64], b: &[f64]) -> f64 {
    let mut r = vec![0.0; b.len()];
    a.spmv(x, &mut r);
    let d: Vec<f64> = b.iter().zip(&r).map(|(b, r)| b - r).collect();
    dot(&d, &d).sqrt() / dot(b, b).sqrt()
}

fn main() {
    for k in [32usize, 64] {
        let a = poisson2d(k);
        let l = ic0(&a).unwrap();
        let b = vec![1.0; a.n];
        let lower = (0..a.n).map(|i| a.col[a.row_ptr[i]..a.row_ptr[i + 1]].iter().filter(|&&c| c as usize <= i).count()).sum::<usize>();
        assert_eq!(l.val.len(), lower); // same pattern as tril(A)
        let (x0, it0) = pcg(&a, &b, None, 1e-8, 10_000);
        let (x1, it1) = pcg(&a, &b, Some(&l), 1e-8, 10_000);
        assert!(rel_res(&a, &x0, &b) < 1e-8 && rel_res(&a, &x1, &b) < 1e-8 && it1 * 3 < it0 * 2);
        println!("k={k}: CG {it0} it, IC(0)-PCG {it1} it");
    }
    // 1-D case: tridiagonal A has no fill, so IC(0) is the exact Cholesky factor
    let t = Csr::from_rows((0..200).map(|i: usize| {
        [(i.wrapping_sub(1), -1.0), (i, 2.0), (i + 1, -1.0)].into_iter().filter(|&(c, _)| c < 200).collect()
    }).collect());
    let b: Vec<f64> = (0..200).map(|i| (i as f64).sin()).collect();
    let (x, it) = pcg(&t, &b, Some(&ic0(&t).unwrap()), 1e-8, 100);
    assert_eq!(it, 1);
    assert!(rel_res(&t, &x, &b) < 1e-12);
    println!("tridiagonal n=200: IC(0)-PCG {it} iteration");
}

Output: k=32 CG 59 / PCG 29 iterations; k=64 119 / 52; tridiagonal 1. At k = 128, 256: CG 239 / 470, PCG 100 / 176.

Complexity and constants

• Factor: one pass over the lower triangle; memory nnz(tril A).
• Apply: two sequential sweeps over L, about one extra SpMV of traffic.
• Poisson: CG doubles per grid doubling (59, 119, 239, 470); IC(0)-PCG grew 29, 52, 100, 176 (1.8-1.9x): a constant-factor gain, not a new order.

Efficiency

• Cost: setup is one sparse pass; each iteration adds about nnz(L) multiply-adds in two sequential sweeps (memory- and latency-bound). One noisy run: k = 64 took about 2.4 ms for CG and 3.3 ms for IC(0)-PCG; sequential solves ate the saving.
• Levers:
  • factor once, reuse across right-hand sides and Newton steps — setup amortized — slowly changing A
  • RCM or minimum-degree ordering — less fill for the same drop rule — ILU(k), ILUT
  • level scheduling or multicolor ordering — rows of one level or color solved in parallel (Saad) — GPUs, many cores
  • Chow and Patel fixed-point sweeps — factor entries in parallel — GPU setup
• Naive → standard → state of the art:
  • Jacobi · one multiply per row
  • IC(0)/ILU(0) · about 2-2.7x fewer iterations on Poisson (measured), sequential solves
  • ILUT/ILU(k) with ordering and parallel solves · fewer iterations, memory grows with fill

Implementation notes

• The diagonal is last in each row; the back solve scatters rows of L, so no transpose. Sort columns per row.
• Shift: start with a small α and grow it geometrically until no pivot fails; IC(0) needs none for M-matrices (Meijerink and van der Vorst).
• Level scheduling: compute dependency levels once; block-Jacobi with ILU per block is the cheap parallel option.

Pitfalls

• A bad pivot appears as NaN or a stalled PCG; check pivots.
• Weak diagonal dominance (convection) can make L⁻¹U⁻¹ huge.
• Ordering changes iteration counts.

Variants and alternatives

In production

Structural finite-element codes, reservoir and circuit simulators and CFD pressure solves use these as preconditioners or smoothers.

Implementations

Sources

• Yousef Saad, Iterative Methods for Sparse Linear Systems, 2nd ed., SIAM, 2003. book · free PDF
• Yousef Saad, ILUT: A dual threshold incomplete LU factorization, Numerical Linear Algebra with Applications 1(4), 1994. paper
• J. A. Meijerink, H. A. van der Vorst, An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix, Mathematics of Computation 31(137), 1977. paper
• Thomas A. Manteuffel, An incomplete factorization technique for positive definite linear systems, Mathematics of Computation 34(150), 1980. paper
• Edmond Chow, Aftab Patel, Fine-Grained Parallel Incomplete LU Factorization, SIAM Journal on Scientific Computing 37(2), 2015. paper