Feat/015 gmres restarted resolver

docs/SUMMARY.md

@@ -32,6 +32,8 @@
     - [Power Iteration](book/08-krylov/power-iteration.md)
     - [Conjugate Gradient](book/08-krylov/conjugate-gradient.md)
     - [Lanczos Iteration](book/08-krylov/lanczos.md)
+    - [Arnoldi Iteration](book/08-krylov/arnoldi.md)
+    - [GMRES(m)](book/08-krylov/gmres.md)
   - [Eigenvalues & Eigenvectors](book/09-eigenvalues.md)
   - [Numerical Stability](book/12-numerical-stability/README.md)
     - [NaN/Inf Policy](book/12-numerical-stability/nan-inf-policy.md)

docs/book/08-krylov/README.md

@@ -7,6 +7,11 @@
 - **Conjugate Gradient (CG)** for solving symmetric positive-definite linear systems
 - **Lanczos iteration**, which builds an orthonormal basis of a Krylov subspace and the
   symmetric tridiagonal matrix a symmetric operator projects onto within it
+- **Arnoldi iteration**, the non-symmetric counterpart to Lanczos, which builds an orthonormal
+  basis of a Krylov subspace and the upper Hessenberg matrix a general operator projects onto
+  within it
+- **GMRES(m)**, a restarted iterative solver for general (possibly non-symmetric) linear
+  systems, built on top of Arnoldi iteration
 
 Unlike the direct decompositions in [Decompositions](../06-decompositions/README.md), these
 refine an estimate (or a basis) over many iterations and can fail to converge — or, for

docs/book/08-krylov/arnoldi.md

@@ -0,0 +1,48 @@
+# Arnoldi Iteration
+
+`arnoldi` reduces a general, row-major `n x n` matrix `a` to upper Hessenberg form over a
+`K`-dimensional Krylov subspace: starting from a normalized `v0`, it builds an orthonormal
+basis `Q` of `span{v0, a*v0, ..., a^{K-1}*v0}` and returns the projection `H = Qᵗ * a * Q`,
+which is upper Hessenberg. Unlike [Lanczos Iteration](lanczos.md), `a` need not be symmetric —
+there is no three-term recurrence to exploit, so every step orthogonalizes the candidate
+vector against the *entire* basis built so far, by modified Gram-Schmidt, rather than just the
+two previous vectors.
+
+```rust
+{{#include ../../../examples/krylov/arnoldi.rs}}
+```
+
+## Orthogonalization
+
+Modified Gram-Schmidt (subtracting each projection immediately, rather than computing all
+projections against the original candidate and subtracting them at the end) is the standard
+trade for Arnoldi: markedly more stable than classical Gram-Schmidt at the same `O(K * n)`
+per-step cost, though still less stable than Householder Arnoldi, which trades that extra
+stability for losing the explicit basis vectors the Krylov projection needs. No
+reorthogonalization pass is added on top, unlike Lanczos's full reorthogonalization — Arnoldi's
+every-vector orthogonalization does not erode as quickly as Lanczos's three-term recurrence
+does, so a second pass doesn't earn its keep the same way.
+
+## Breakdown is success, not failure
+
+When the candidate vector's norm falls to (numerically) zero relative to `‖a * q_j‖` after
+orthogonalization, `q_0, ..., q_j` already span an invariant subspace of `a`: there's no new
+direction to extend the basis with, but the vectors and the leading block of `H` already built
+are exact and useful. This is reported as `Ok((h, reached))` with `reached < K`, not an error —
+unlike Lanczos, which reports the analogous condition as `ConvergenceError::Breakdown`. The
+difference is what the caller does next: a Lanczos caller that requested `K` vectors and got
+fewer has nothing it can use without changing its request, while GMRES, built on top of
+Arnoldi, can solve directly in the smaller subspace `Ok` reports — the exact subspace is often
+exactly where the true solution already lives. Callers that do need the full `K` vectors
+distinguish this case from a complete run by checking `reached < K`.
+
+## Gotchas
+
+- `K` is a `const` generic on the caller's `Basis` buffer, not a runtime parameter — see
+  [Krylov Basis-Size Const-Generic Convention](../../specs/krylov-basis-size-const-generics.md).
+  `K > n` is a `DimensionMismatch`.
+- Both `ConvergenceError::ZeroVector` and `ConvergenceError::NonFinite` on `v0` are checked up
+  front, even when `K == 0` means no basis vector is ever written.
+- `tol` has no auto-computed default — see
+  [Krylov Tolerance and Convergence Criteria](../../specs/krylov-tolerance-and-convergence.md).
+  A `tol` of `0` detects only exact breakdown.

docs/book/08-krylov/gmres.md

@@ -0,0 +1,57 @@
+# GMRES(m)
+
+`gmres` solves the general (possibly non-symmetric) linear system `A x = b` via restarted
+GMRES, GMRES(`M`): unlike [Conjugate Gradient](conjugate-gradient.md), `A` need not be
+symmetric positive-definite. Each restart cycle runs [Arnoldi Iteration](arnoldi.md) 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.
+
+```rust
+{{#include ../../../examples/krylov/gmres.rs}}
+```
+
+## 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](arnoldi.md).
+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](../../specs/krylov-basis-size-const-generics.md).
+  `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](../../specs/krylov-tolerance-and-convergence.md).
+- 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`.

examples/krylov/gmres.rs

@@ -0,0 +1,57 @@
+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:?}");
+}

examples/krylov/main.rs

@@ -1,14 +1,19 @@
 //! Tour `rustebra::krylov`'s subspace methods: Lanczos iteration, which builds an orthonormal
 //! basis of a Krylov subspace and the symmetric tridiagonal matrix a symmetric operator
-//! projects onto within it, and Arnoldi iteration, its non-symmetric counterpart producing an
-//! upper Hessenberg projection.
+//! projects onto within it, Arnoldi iteration, its non-symmetric counterpart producing an
+//! upper Hessenberg projection, and GMRES(m), the restarted linear solver built on top of it.
 //!
-//! Run with: `cargo run --example krylov`
+//! Run with: `cargo run --example krylov` (add `--features alloc` for the GMRES section, since
+//! it takes its operator as a `SparseLinearOp`).
 
 mod arnoldi;
+#[cfg(feature = "alloc")]
+mod gmres;
 mod lanczos;
 
 fn main() {
     lanczos::run();
     arnoldi::run();
+    #[cfg(feature = "alloc")]
+    gmres::run();
 }