Crate faer_precond

Expand description

Numerical preconditioners for iterative linear solvers, built on faer.

Every preconditioner in this crate implements faer’s matrix_free::{LinOp, Precond, BiLinOp, BiPrecond} traits, so they plug directly into faer’s Krylov solvers (CG, GMRES, BiCGSTAB, LSMR, …).

§Available preconditioners

Type	Source	Apply cost	Notes
`JacobiPrecond`	diagonal of `A`	`O(n)`	Diagonally-dominant problems.
`BlockJacobiPrecond`	dense diagonal blocks of `A`	`O(sum b_k²)`	Arbitrary block partition; LU per block.
`Ssor`	`A`’s `D`/`L`/`U` split + relaxation	`O(nnz(A))`	Stationary; two triangular solves, no fill.
`Ilu0`	CSC sparsity of `A`	`O(nnz(A))`	Zero-fill incomplete LU.
`Iluk`	CSC of `A` + fill level `k`	`O(nnz_LU)`	Level-of-fill incomplete LU; between ILU(0) and ILUTP.
`Ic0`	CSC lower triangle of `A`	`O(nnz_L)`	Zero-fill incomplete Cholesky for HPD `A`.
`Ict`	CSC of `A` + threshold/fill params	`O(nnz_L)`	Threshold incomplete Cholesky; SPD analogue of ILUTP.
`Ilutp`	CSC of `A` + threshold/fill params	`O(nnz_LU)`	Threshold ILU with partial pivoting; general nonsymmetric workhorse.
`Poly`	`A` + polynomial degree/bounds	`O(degree · nnz(A))`	Matvec-only (Neumann/Chebyshev); no triangular solves.
`Fsai`	CSC of `A` + pattern	`O(nnz_G)`	Factorised approximate inverse for HPD `A`; matvec apply.
`Spai`	CSC of `A` + pattern	`O(nnz_M)`	Sparse approximate inverse (nonsymmetric); matvec apply.
`SolvePrecond`	any faer factorisation (`Llt`, `Lu`, `Qr`, …)	factorisation-dependent	Adapter, not a factorisation.

§Choosing a preconditioner

There is no single best preconditioner; the right one depends on the structure of A and how much work you can afford per iteration.

Start with JacobiPrecond. It is almost free to build and apply, and it helps whenever A’s rows differ in scale — diagonally dominant systems, variable-coefficient PDEs, badly-scaled unknowns. If A has a constant diagonal it does nothing, so move on.
Ic0 for symmetric positive-definite A. The standard choice for SPD problems from PDE discretisations (Laplacians, diffusion, elasticity) solved with conjugate gradient. It cuts iteration counts sharply; each apply costs two sparse triangular solves, so it always wins on iteration count and wins on wall-clock time once the problem is ill-conditioned enough to need many iterations.
Ilu0 for general (nonsymmetric) sparse A. The nonsymmetric counterpart to IC(0), paired with GMRES or BiCGSTAB. Same zero-fill idea: cheap to build and stores nothing beyond A’s sparsity pattern.
Ilutp when Ilu0 is too weak. Threshold ILU with partial pivoting: it adds fill where the factor needs it (tuned by a drop tolerance and a fill budget) and pivots for stability — the robust choice for hard nonsymmetric problems, badly-scaled operators, or matrices with small/zero diagonal entries. Costs more to build and apply than Ilu0, and its pattern is value-dependent (no zero-allocation refactorisation).
BlockJacobiPrecond when unknowns cluster into small dense groups. Several fields per mesh node, coupled species, or tightly-coupled sub-systems. Inverting those blocks exactly captures the strong local coupling that point-Jacobi misses.
Ssor for a cheap, fill-free stationary preconditioner. Built straight from A’s own triangles with one relaxation knob; a step up from point-Jacobi when you do not want to store a factorisation. Symmetric for SPD A, so it pairs with CG.
Iluk when Ilu0 is too weak but you want a fixed pattern. Level-of-fill ILU: more accurate than ILU(0), but with a value-independent pattern (allocation-free refactorisation), unlike Ilutp.
Ict for hard SPD problems. The threshold incomplete Cholesky: IC(0)’s adaptive cousin, for ill-conditioned SPD systems where zero-fill stalls. The SPD counterpart to Ilutp.
Poly, Fsai or Spai when triangular solves are the bottleneck. These apply through matrix-vector products only — no sequential forward/back substitution — so they parallelise far better on many cores or accelerators. Poly (Neumann/Chebyshev) and Fsai are for SPD A; Spai is the nonsymmetric approximate inverse. The trade-off is a weaker approximation per unit of single-core work; Poly’s Chebyshev form also needs a spectral-interval estimate.
SolvePrecond to reuse an exact factorisation. Factorise a cheaper approximation of A once and let the Krylov method correct the rest.

When A’s values change between iterations but its sparsity pattern does not — the inner loop of a nonlinear solver — build the symbolic factor once and call refactorize (see Ilu0 and Ic0) to avoid reallocating.

§Design contract

No heap allocation during apply. Every preconditioner returns either StackReq::EMPTY or a precise pre-computed scratch size from apply_scratch / apply_in_place_scratch. All temporary memory flows through the MemStack faer’s trait interface provides.
Refactorisation reuses storage. Ilu0 and Ic0 expose a refactorize(&mut self, a) method for the steady-state case (nonlinear Krylov drivers where the sparsity pattern is fixed but the values change). It allocates nothing.
In-place semantics match faer. apply performs out = M^{-1} rhs; apply_in_place overwrites rhs with M^{-1} rhs. Transpose and adjoint variants follow the same convention against M^{-T} and M^{-H}.

§Example: point-Jacobi

use dyn_stack::MemStack;
use faer::{mat, Par};
use faer::matrix_free::Precond;
use faer_precond::JacobiPrecond;

let pc = JacobiPrecond::try_from_diagonal(&[4.0_f64, 2.0, 8.0]).unwrap();

let mut x = mat![[8.0_f64], [6.0], [16.0]];
pc.apply_in_place(x.as_mut(), Par::Seq, MemStack::new(&mut []));

assert!((*x.as_ref().get(0, 0) - 2.0).abs() < 1e-12);
assert!((*x.as_ref().get(1, 0) - 3.0).abs() < 1e-12);
assert!((*x.as_ref().get(2, 0) - 2.0).abs() < 1e-12);

§Example: incomplete LU on a sparse matrix

use dyn_stack::MemStack;
use faer::sparse::{SparseColMat, Triplet};
use faer::{mat, Par};
use faer::matrix_free::Precond;
use faer_precond::Ilu0;

// Build a 5x5 SPD tridiagonal: diag = 4, off = -1.
let mut triplets = Vec::new();
for i in 0..5 {
    triplets.push(Triplet::new(i, i, 4.0_f64));
    if i > 0 {
        triplets.push(Triplet::new(i, i - 1, -1.0));
        triplets.push(Triplet::new(i - 1, i, -1.0));
    }
}
let a = SparseColMat::<usize, f64>::try_new_from_triplets(5, 5, &triplets).unwrap();

let pc = Ilu0::try_new(a.as_ref()).expect("non-singular pattern");

// For a tridiagonal there is no fill, so ILU(0) is the exact LU.
let mut b = mat![[1.0_f64], [0.0], [0.0], [0.0], [0.0]];
pc.apply_in_place(b.as_mut(), Par::Seq, MemStack::new(&mut []));

§Cargo features

There are no opt-in feature flags in 0.1. The crate inherits faer’s default feature set, which includes sparse linear algebra (sparse-linalg) required by Ilu0 and Ic0.

§License

Dual-licensed under MIT OR Apache-2.0.

Re-exports§

pub use adapters::SolvePrecond;
pub use block_jacobi::BlockJacobiError;
pub use block_jacobi::BlockJacobiPrecond;
pub use fsai::Fsai;
pub use fsai::FsaiError;
pub use fsai::FsaiPattern;
pub use ic0::Ic0;
pub use ic0::Ic0Error;
pub use ic0::SymbolicIc0;
pub use ict::Ict;
pub use ict::IctError;
pub use ict::IctParams;
pub use ilu0::Ilu0;
pub use ilu0::Ilu0Error;
pub use ilu0::SymbolicIlu0;
pub use iluk::Iluk;
pub use iluk::IlukError;
pub use iluk::IlukParams;
pub use iluk::SymbolicIluk;
pub use ilutp::FillControl;
pub use ilutp::Ilutp;
pub use ilutp::IlutpError;
pub use ilutp::IlutpParams;
pub use ilutp::RowNorm;
pub use jacobi::JacobiError;
pub use jacobi::JacobiPrecond;
pub use poly::BoundEstimate;
pub use poly::Poly;
pub use poly::PolyError;
pub use poly::PolyKind;
pub use poly::PolyParams;
pub use spai::Spai;
pub use spai::SpaiError;
pub use spai::SpaiPattern;
pub use ssor::Ssor;
pub use ssor::SsorError;
pub use ssor::SsorParams;

Modules§

adapters
block_jacobi: Block-Jacobi preconditioner.
fsai: Factorized sparse approximate inverse (FSAI) preconditioner.
ic0: Zero-fill incomplete Cholesky preconditioner.
ict: Threshold incomplete Cholesky preconditioner (ICT).
ilu0: Zero-fill incomplete LU preconditioner.
iluk: Level-of-fill incomplete LU preconditioner, ILU(k).
ilutp: Threshold ILU with partial pivoting (ILUTP).
jacobi: Point-Jacobi (diagonal) preconditioner.
poly: Polynomial preconditioner (Neumann series and Chebyshev).
spai: Sparse approximate inverse (SPAI) preconditioner.
ssor: SSOR / symmetric Gauss-Seidel preconditioner.