Module krylov 

Layer 2 — asymmetric Krylov iterative solvers (BiCGStab, restarted GMRES) and preconditioners (Jacobi, ILU(0)) for the sparse LduMatrix. Asymmetric Krylov iterative solvers and preconditioners for sparse A x = b.

This module complements the crate’s existing SPD-only machinery (DIC-PCG and GAMG in crate::ldu_matrix::solvers) with the nonsymmetric iterative solvers a Newton–Krylov subsurface-flow solver needs, where the Jacobian is not symmetric:

• bicgstab — preconditioned BiCGStab: fixed work/storage per iteration, breakdown-guarded.
• gmres — restarted, right-preconditioned GMRES(m): residual-minimising, robust, O(m) storage.

and three preconditioners dispatched by the Preconditioner enum (never trait objects, per the workspace design rules):

• Preconditioner::identity — no preconditioning (M = I).
• Preconditioner::jacobi — diagonal scaling; always applicable.
• Preconditioner::ilu0 — genuine ILU(0) incomplete factorisation.

Matrix representation and conventions

All solvers act on crate::ldu_matrix::LduMatrix, the crate’s face-addressed sparse matrix, and use only its multiply (SpMV) and residual kernels. The system size n is LduMatrix::n_cells; all right-hand-side, guess, and solution slices have length n. Every quantity here is a dimensionless f64: a Krylov subspace mixes residuals, search directions and increments that share no single physical dimension, so no uom typing is applied — apply units at the field/equation layer that assembles the matrix.

Execution backend

Both solvers run on the hybrid ComputeBackend, driving the kernels in crate::ldu_matrix::parallel. Each has one implementation with the backend as a parameter — bicgstab_prepared and gmres_prepared, which take a HybridLdu so the cell-gather index is built once per mesh rather than once per solve — plus a convenience adapter (bicgstab, gmres) for a caller holding a bare LduMatrix, which builds the index and runs on ComputeBackend::Serial. They are not a serial/parallel pair: they differ in who owns the index, and the backend is a parameter of both.

Whole-solve wall clock on 4 logical cores, release, --features parallel, 512 000 cells, measured 2026-08-13 over five independent runs: 2.4-2.7x with Jacobi preconditioning, ~1.5x with ILU(0) — the gap being ILU(0)’s inherently sequential triangular solves, measured at 29-46% of the serial solve, which caps it at 1.7-2.1x on 4 cores. Below roughly 13 000 cells the parallel path loses (0.6-0.8x at 4 096 cells). Full tables, methodology and limitations are on the benchmarks in hybrid_tests and summarised on bicgstab_prepared.

Backend parity is bitwise, not tolerance-based: Serial and CpuMulti produce identical iterates, identical residual histories and identical iteration counts at any thread count, because every kernel underneath carries that guarantee individually.

What is not bitwise is the comparison against the solvers as they stood before they took the hybrid path, because crate::ldu_matrix::parallel::dot sums in blocks of 1 024 where crate::krylov::vecops::dot sums flat, and floating-point addition is not associative. Below 1 024 elements the two are identical; above it they differ in the last bits, which in principle can move an iteration count by one when a residual crosses the tolerance within those bits. Measured rather than assumed: over a 216-solve sweep of 3 mesh sizes x 3 diagonal dominances x 3 right-hand sides x 4 tolerances (1e-6 to 1e-12) x both solvers, the iteration count changed zero times, for GMRES — the more exposed of the two — as well as for BiCGStab. Converged residuals differ by up to 2.2e-4 relative, which at the tightest tolerance swept is an absolute difference of ~1.9e-16 on a residual of 8.6e-13. See hybrid_tests for the full grid and its limitations.

Convergence

The stopping test for both solvers is the relative residual ||b − A x||₂ / ||b||₂ ≤ tolerance. The reported final_residual is always the true residual of the returned iterate (recomputed from A and b), not an internal estimate. A right-hand side that is exactly zero returns x = 0, converged = true, 0 iterations.

Example

use outram_foam_basic_lib::ldu_matrix::LduMatrix;
use outram_foam_basic_lib::krylov::{bicgstab, Preconditioner, KrylovSettings};

// 3-cell chain: cells 0-1 and 1-2 share a face each.
let mut a = LduMatrix::new(3, vec![0, 1], vec![1, 2]);
a.diag = vec![4.0, 4.0, 4.0];
a.lower = vec![-1.0, -1.0];
a.upper = vec![-1.0, -1.0];
let b = vec![1.0, 2.0, 3.0];

let precond = Preconditioner::jacobi(&a);
let settings = KrylovSettings::default();
let (x, result) = bicgstab(&a, &b, None, &precond, &settings);

assert!(result.converged);
// Verify: A x ≈ b.
let ax = a.multiply(&x);
for i in 0..3 {
    assert!((ax[i] - b[i]).abs() < 1e-6);
}

Modules

vecops

Dense BLAS-1 vector primitives for the Krylov solvers.

Structs

Ilu0Preconditioner

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

JacobiPreconditioner

Diagonal (Jacobi) preconditioner: M = diag(A), so M^{-1} r = r / diag.

KrylovResult

Outcome of a Krylov solve.

KrylovSettings

Iteration controls shared by bicgstab and gmres.

Enums

Preconditioner

Preconditioner M^{-1} ≈ A^{-1}, dispatched by enum (no trait objects).

Functions

bicgstab

Solve A x = b with preconditioned BiCGStab, serially, from a bare LduMatrix.

bicgstab_prepared

Solve A x = b with preconditioned BiCGStab on a chosen ComputeBackend.

gmres

Solve A x = b with restarted right-preconditioned GMRES(m), serially, from a bare LduMatrix.

gmres_prepared

Solve A x = b with restarted right-preconditioned GMRES(m) on a chosen ComputeBackend.