Struct Ilu0Preconditioner 

pub struct Ilu0Preconditioner { /* private fields */ }

ILU(0) preconditioner: incomplete LU factorisation with no fill-in.

The factors L (unit lower-triangular) and U (upper-triangular) are computed with the standard IKJ Gaussian-elimination algorithm but every update that would land outside the sparsity pattern of A is discarded — so L\U occupy exactly the same nonzero positions as A. Applying M^{-1} = (LU)^{-1} is one forward substitution (L y = r) followed by one back substitution (U z = y).

This is a genuine ILU(0), not a Jacobi fallback

The factorisation is real: for a symmetric problem it reduces to an incomplete Cholesky-like factor and typically beats Jacobi substantially; for a nonsymmetric problem it captures both off-diagonal directions. The only safeguard is a small-pivot floor (ILU_PIVOT_FLOOR): a pivot smaller than that in magnitude is bumped up (keeping its sign) so the triangular solves stay finite. In the degenerate limit where every off-diagonal contribution vanishes, ILU(0) coincides exactly with Jacobi — that is a property of the method, not a substitution of one for the other.

Storage (CSR, built internally from the LDU coefficients)

The LDU (face-addressed) matrix is expanded once, at construction, into a compressed-sparse-row (CSR) view: for each row i the diagonal diag[i], plus upper[f] at column neighbour[f] for every face owned by i, plus lower[f] at column owner[f] for every face neighboured by i. Columns within a row are stored in ascending order. The factorisation overwrites those values with the combined L\U (unit L diagonal is implicit and not stored).

Implementations

impl Ilu0Preconditioner

pub fn new(a: &LduMatrix) -> Self

Build the ILU(0) factorisation of a.

Constructs the CSR view then runs the in-place IKJ incomplete-LU elimination restricted to existing nonzeros. Runs in roughly O(nnz · avg_row) time. Pivots below ILU_PIVOT_FLOOR are floored to keep the factor invertible.

pub fn apply_on(&self, r: &[f64], z: &mut [f64], _backend: ComputeBackend)

Apply z = M^{-1} r — always sequentially, on every backend.

This kernel does not parallelise, and the code does not pretend it does

ILU(0) is sequential by construction in both of its phases:

• The factorisation (Self::factor) has a loop-carried dependence: eliminating row i reads the already-eliminated rows k < i through self.lu[self.diag_ptr[k]] and find(k, j), so row i cannot start until every row it touches has finished.
• The triangular solves (Self::apply) have the same shape: the forward solve reads z[j] for j < i while computing z[i], and the back solve reads z[j] for j > i.

There are parallel variants in the literature — level-scheduling, which exposes only as much parallelism as the sparsity graph’s critical path allows, and block-Jacobi-over-subdomains, which is a different preconditioner: it drops the couplings between subdomains, so it changes the iteration count and needs its own validation. Bead op-yvj.4.4 names that explicitly and says not to substitute one silently for the other, so this implementation does neither. It keeps the exact standard ILU(0) and runs it on one thread.

The practical consequence is Amdahl’s law: in a solve preconditioned with ILU(0), the fraction of time spent here is a hard serial floor on the achievable speed-up. That fraction is measured, not assumed: on an asymmetric 7-point-stencil system, 4 logical cores, release, --features parallel, 2026-08-13, these triangular solves were 29-46% of the serial solve time (rising with problem size), which caps an ILU(0)-preconditioned speed-up at 1.7-2.1x on 4 cores and 2.2-3.4x on infinitely many. The measured end-to-end figure is ~1.5x. Full derivation, all four repeat runs and the limitations are on the ilu0_serial_fraction_benchmark in crate::krylov::hybrid_tests.

The factorisation (Self::new) is a further sequential cost on top: at 262 144 cells it measured 124-151 ms against a 522-664 ms serial solve, i.e. roughly a quarter of one solve — 23.8% in the fourth run, the only one whose factorisation and solve figures are quoted here as a matched pair — paid once per coefficient update rather than once per iteration.

Arguments

• r — residual to precondition, length n, dimensionless.
• z — output buffer, length n. Fully overwritten.
• _backend — accepted so this matches JacobiPreconditioner::apply_on and the enum can dispatch uniformly. It is ignored; the kernel runs serially whatever is requested. That also means the result is trivially bitwise identical across backends.