faer-precond

Numerical preconditioners for iterative linear solvers, built on top of the faer linear algebra crate.

Every preconditioner implements faer's LinOp, Precond, BiLinOp, and BiPrecond traits from the matrix_free module, so they plug directly into faer's Krylov solvers (CG, GMRES, BiCGSTAB, LSMR, ...).

Preconditioners

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. Almost free to build and apply. 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 is two sparse triangular solves, so it always wins on iteration count and wins on wall-clock time once the problem is ill-conditioned enough.
• Ilu0 for general (nonsymmetric) sparse A. The nonsymmetric counterpart to IC(0), paired with GMRES or BiCGSTAB. Cheap to build, 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 diagonals. Costs more than Ilu0, and its pattern is value-dependent (no zero-allocation refactorisation). See examples/ilutp.rs for the fill-vs-iterations tradeoff.
• BlockJacobiPrecond when unknowns cluster into small dense groups. Several fields per mesh node, coupled species, tightly-coupled sub-systems. Inverting those blocks exactly captures the strong local coupling point-Jacobi misses.
• SolvePrecond to reuse an exact factorisation. Factorise a cheaper approximation of A once and let the Krylov method correct the rest.

The examples/speedup.rs example measures all of this end-to-end on a badly-conditioned diffusion problem — run it with cargo run --release --example speedup.

Design contract

• No heap allocation during apply. Every apply_scratch returns either StackReq::EMPTY or a precise pre-computed size; all temporary memory flows through the MemStack provided by faer's trait interface.
• Refactorisation reuses storage. Ilu0::refactorize and Ic0::refactorize mutate the existing value buffer for the next iteration of a nonlinear Krylov driver — no allocation occurs after the first factorisation.
• In-place semantics match faer. apply performs out = M⁻¹ rhs; apply_in_place overwrites rhs.

Install

[dependencies]
faer-precond = "0.1"
faer = "0.24"
dyn-stack = "0.13"

Quick start

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 []));
// x is now [2.0, 3.0, 2.0]

Block-Jacobi

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

// 5x5 matrix; two diagonal blocks of size 2 and 3.
let a = mat![
    [4.0_f64, 1.0, 0.0, 0.0, 0.0],
    [2.0,     3.0, 0.0, 0.0, 0.0],
    [0.0,     0.0, 6.0, 1.0, 2.0],
    [0.0,     0.0, 3.0, 5.0, 1.0],
    [0.0,     0.0, 2.0, 1.0, 4.0],
];
let pc = BlockJacobiPrecond::try_new(a.as_ref(), &[0, 2, 5]).unwrap();

let mut x = mat![[1.0_f64], [2.0], [3.0], [-1.0], [0.5]];
pc.apply_in_place(x.as_mut(), Par::Seq, MemStack::new(&mut []));

ILU(0)

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

// 5x5 tridiagonal SPD: diag 4, off -1 (no fill — ILU(0) is exact).
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");

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 []));

IC(0)

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

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();

// Hermitian PD input — IC(0) silently ignores the strict upper triangle.
let pc = Ic0::try_new(a.as_ref()).expect("matrix is positive definite");

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 []));

ILUTP (threshold + partial pivoting)

use dyn_stack::{MemBuffer, MemStack};
use faer::sparse::{SparseColMat, Triplet};
use faer::{mat, Par};
use faer::matrix_free::Precond;
use faer_precond::{FillControl, Ilutp, IlutpParams};

// A non-symmetric tridiagonal: diag 4, sub -2, super -1.
let mut triplets = Vec::new();
for i in 0..5usize {
    triplets.push(Triplet::new(i, i, 4.0_f64));
    if i > 0 {
        triplets.push(Triplet::new(i, i - 1, -2.0));
        triplets.push(Triplet::new(i - 1, i, -1.0));
    }
}
let a = SparseColMat::<usize, f64>::try_new_from_triplets(5, 5, &triplets).unwrap();

// Tune the drop tolerance, fill budget, and pivot aggressiveness.
let params = IlutpParams {
    drop_tol: 1e-3,
    fill: FillControl::PerRow(10),
    pivot_tol: 0.1,
    ..Default::default()
};
let pc = Ilutp::try_new_with_params(a.as_ref(), params).unwrap();

// Apply needs a length-n scratch column (it may apply the pivot permutation).
let mut b = mat![[1.0_f64], [0.0], [0.0], [0.0], [0.0]];
let mut buf = MemBuffer::new(pc.apply_in_place_scratch(1, Par::Seq));
pc.apply_in_place(b.as_mut(), Par::Seq, MemStack::new(&mut buf));

Wrapping a faer factorisation as a preconditioner

use dyn_stack::MemStack;
use faer::{mat, Par, Side};
use faer::linalg::solvers::Llt;
use faer::matrix_free::Precond;
use faer_precond::SolvePrecond;

let a = mat![[4.0_f64, 1.0], [1.0, 3.0]];
let llt = Llt::new(a.as_ref(), Side::Lower).expect("matrix is SPD");
let pc = SolvePrecond::new(llt);

let mut x = mat![[1.0_f64], [2.0]];
pc.apply_in_place(x.as_mut(), Par::Seq, MemStack::new(&mut []));
// x is now [1/11, 7/11]

Repeated factorisation (nonlinear Krylov)

When A's sparsity pattern is fixed but its values change between Krylov iterations, build the symbolic factor once and refactorise repeatedly with zero allocation:

use faer::sparse::{SparseColMat, Triplet};
use faer_precond::{Ilu0, SymbolicIlu0};

# let triplets: Vec<Triplet<usize, usize, f64>> = (0..3).map(|i| Triplet::new(i, i, 1.0_f64)).collect();
# let a0 = SparseColMat::<usize, f64>::try_new_from_triplets(3, 3, &triplets).unwrap();
let symbolic = SymbolicIlu0::try_new(a0.as_ref().symbolic()).unwrap();
let mut pc = Ilu0::<usize, f64>::new_with_symbolic(symbolic);

// Hot loop — no allocation:
pc.refactorize(a0.as_ref()).unwrap();
// pc.refactorize(a1.as_ref()).unwrap();
// pc.refactorize(a2.as_ref()).unwrap();

Trait coverage

Every preconditioner implements the full trait hierarchy:

• LinOp<T> — apply, conj_apply
• Precond<T> — apply_in_place, conj_apply_in_place
• BiLinOp<T> — transpose_apply, adjoint_apply
• BiPrecond<T> — transpose_apply_in_place, adjoint_apply_in_place

Minimum supported Rust version

rustc 1.88 (edition 2024 + let-chains).

License

Dual-licensed under either of

• Apache License, Version 2.0 — LICENSE-APACHE
• MIT license — LICENSE-MIT

at your option.

Contribution

Unless you explicitly state otherwise, any contribution intentionally submitted for inclusion in the work by you, as defined in the Apache-2.0 license, shall be dual-licensed as above, without any additional terms or conditions.