GMRES(m)

gmres solves the general (possibly non-symmetric) linear system A x = b via restarted GMRES, GMRES(M): unlike Conjugate Gradient, A need not be symmetric positive-definite. Each restart cycle runs Arnoldi Iteration from the current residual to build an M-dimensional Krylov basis, solves the resulting small least-squares problem via Givens rotations, and updates x — restarting from the improved iterate until either the residual meets tol or max_restarts cycles are exhausted.

A is supplied as a sparse linear operator rather than a dense matrix: applying it never allocates, so restart cycles reuse the same workspace (out_x, basis, scratch) throughout.

#![allow(unused)]
fn main() {
use rustebra::krylov::gmres;
use rustebra::sparse::CsrMatrix;
use rustebra::storage::Basis;

pub(crate) fn run() {
    println!("\n== GMRES(m) ==");

    // Non-symmetric, non-SPD system, unlike Conjugate Gradient's requirements:
    // [[4, 1], [2, 3]] x = [1, 2]. Solution: x = [0.1, 0.6].
    let a = CsrMatrix::new(
        2,
        2,
        vec![0, 2, 4],
        vec![0, 1, 0, 1],
        vec![4.0_f64, 1.0, 2.0, 3.0],
    )
    .expect("valid CSR");
    let b = [1.0, 2.0];
    let x0 = [0.0, 0.0];
    let mut out_x = [0.0; 2];
    let mut buffer = [0.0; 4];
    let mut basis = Basis::<f64, 2>::new(&mut buffer, 2).unwrap();
    let mut scratch = [0.0; 2];

    gmres(&a, &b, &x0, 10, 1e-10, &mut out_x, &mut basis, &mut scratch).expect("converges");
    println!("full-basis solve: x = {out_x:?}");

    // Restart size smaller than the problem dimension (M = 1 < n = 3): the solution still
    // emerges, just carried forward across several restart cycles instead of one.
    let a3 = CsrMatrix::new(
        3,
        3,
        vec![0, 2, 5, 7],
        vec![0, 1, 0, 1, 2, 1, 2],
        vec![5.0_f64, 1.0, 1.0, 4.0, 1.0, 1.0, 3.0],
    )
    .expect("valid CSR");
    let b3 = [6.0, 6.0, 4.0];
    let x0_3 = [0.0, 0.0, 0.0];
    let mut out_x3 = [0.0; 3];
    let mut buffer3 = [0.0; 3];
    let mut basis3 = Basis::<f64, 1>::new(&mut buffer3, 3).unwrap();
    let mut scratch3 = [0.0; 3];

    gmres(
        &a3,
        &b3,
        &x0_3,
        500,
        1e-10,
        &mut out_x3,
        &mut basis3,
        &mut scratch3,
    )
    .expect("converges across restarts");
    println!("restarted GMRES(1) solve: x = {out_x3:?}");
}
}

Algorithm

Each restart cycle:

1. Computes the residual r = b - A x and its norm β = ‖r‖, returning Ok immediately if β <= tol.
2. Runs Arnoldi iteration from q_0 = r / β, building an orthonormal basis Q of up to M vectors and the upper Hessenberg projection H, stopping early (before M steps) on breakdown — an invariant subspace found before the basis filled up, the same “success, not failure” case documented on Arnoldi Iteration.
3. Solves min_y ‖β e_1 - H y‖ via incremental Givens rotations, then updates x <- x + Q y.

Convergence

GMRES’s residual norm decreases monotonically within a cycle (each additional basis vector can only improve the least-squares fit) and never increases across a restart, because restarting recomputes the same residual the next cycle continues from. It is not guaranteed to decrease strictly every cycle, though: a starting vector aligned with an invariant subspace the operator doesn’t expand (breakdown on the very first Arnoldi step) leaves x unchanged, and the iteration stagnates. Restarting also discards the larger Krylov subspace full (non-restarted) GMRES would have kept building, so GMRES(M) can converge slower, or stagnate on problems full GMRES would resolve — the restart budget max_restarts bounds the cost of that risk rather than eliminating it.

Gotchas

• The restart size M is a const generic on the caller’s Basis buffer, not a runtime parameter — see Krylov Basis-Size Const-Generic Convention. M > n (the operator’s dimension) is a DimensionMismatch.
• M == 0 never reduces a nonzero residual: no basis vector can be built, so every restart cycle is a no-op check against tol, and a nonzero residual exhausts max_restarts without ever moving x.
• tol has no auto-computed default — see Krylov Tolerance and Convergence Criteria.
• Non-finite residuals or Arnoldi iterates (NaN or infinite, including values that overflow f64 mid-computation) return ConvergenceError::NonFinite rather than silently producing a wrong answer.
• A zero pivot in the small least-squares system built from H — a coincidental exact singularity in the projected system that Arnoldi’s own breakdown test didn’t already catch — returns ConvergenceError::Breakdown.