Iterative Linear Solvers

Why Iterative Linear Solvers Matters

A temperature model on a large grid can produce millions of linear equations, while each equation involves only a few neighboring grid cells. A dense factorization would ignore that structure and consume far more memory than storing the nonzero coefficients. Iterative solvers improve a candidate solution through repeated updates or matrix-vector products, often reaching useful accuracy without a full factorization.

This lesson solves real square systems $A\mathbf x=\mathbf b$ using Jacobi, Gauss–Seidel, successive overrelaxation (SOR), conjugate gradient (CG), and generalized minimal residual (GMRES). You will distinguish stationary iterations from Krylov methods, trace a small system exactly, implement all five, and check the true residual before declaring convergence.

The teaching programs use dense nested arrays so every update is visible. Sparse storage and matrix-free products are the natural next step for large problems. An iterative method is not automatically cheaper or more reliable than a direct solve: convergence depends on the matrix, the starting guess, tolerances, and often a preconditioner.

Core Concepts and Terminology

Given $A\in\mathbb R^{n\times n}$, $\mathbf b\in\mathbb R^n$, and an initial vector $\mathbf x_0$, the output is an approximate solution, its residual norm, the number of completed updates, and a termination reason. The programs use

$$ \lVert\mathbf b-A\mathbf x\rVert_2 \le \text{atol}+\text{rtol}\,\lVert\mathbf b\rVert_2, \qquad \text{atol}\gt0,\quad\text{rtol}\ge0. $$

The positive absolute tolerance handles a zero right-hand side without dividing by its norm. This additive threshold is the contract used here; other libraries may use a maximum of absolute and relative thresholds instead.

Use $A=\begin{bmatrix}4&1\\1&3\end{bmatrix}$ and $\mathbf b=[1,2]^T$. The exact solution is $[1/11,7/11]^T$. Starting from zero gives residual $[1,2]^T$. Each method changes the candidate while leaving the original equations fixed.

A small residual is not automatically a small solution error. For nonsingular $A$ and nonzero $\mathbf b$, $\lVert\mathbf e\rVert_2/\lVert\mathbf x_\star\rVert_2\le\kappa_2(A)\lVert\mathbf r\rVert_2/\lVert\mathbf b\rVert_2$. Ill-conditioning can magnify the effect of residual error and data perturbations.

How Iterative Linear Solvers Works, Step by Step

1. Validate the matrix and the method's assumptions. Require a finite, nonempty square matrix and matching finite vectors. Jacobi, Gauss–Seidel, and SOR need usable diagonal entries; the code rejects $\lvert a_{ii}\rvert\le10^{-14}\max_j\lvert a_{ij}\rvert$. Supported SOR parameters satisfy $0\lt\omega\lt2$. GMRES needs a positive restart length. CG checks exact stored symmetry and relies on the caller to supply positive definiteness; checking symmetry and encountered positive curvatures does not certify SPD on every direction.

2. Compute the initial residual and test it. If the starting guess already meets the residual threshold, return converged with zero updates. Otherwise proceed within the iteration budget. All method-specific input requirements are checked before this shortcut. An iteration budget of zero permits the initial check but makes no updates.

3. For stationary methods, choose which values a row update uses. Write $A=D+L+U$, with $L$ and $U$ containing the actual strictly lower and upper entries of $A$.

Jacobi uses only the previous iterate:

$$ x_i^{(k+1)}=\frac{b_i-\sum_{j\ne i}a_{ij}x_j^{(k)}}{a_{ii}}. $$

Gauss–Seidel uses values already computed in the current sweep:

$$ x_i^{(k+1)}=\frac{b_i-\sum_{j\lt i}a_{ij}x_j^{(k+1)}-\sum_{j\gt i}a_{ij}x_j^{(k)}}{a_{ii}}. $$

SOR first forms that same row's Gauss–Seidel candidate $\widehat x_i$, then relaxes it using $x_i^{(k+1)}=(1-\omega)x_i^{(k)}+\omega\widehat x_i$. With $\omega=1$ this is exactly Gauss–Seidel. A larger $\omega$ is not automatically faster, and $0\lt\omega\lt2$ does not guarantee convergence on an arbitrary nonsymmetric matrix.

Each completed sweep is followed by a fresh residual calculation. Reusing the updated array for Jacobi would silently turn it into a different iteration. The implementations keep the old and new vectors separate and select the appropriate source for each off-diagonal term.

4. For CG, update conjugate search directions. Initialize $\mathbf r=\mathbf b-A\mathbf x$, $\mathbf p=\mathbf r$, and $\gamma=\mathbf r^T\mathbf r$. On each step,

$$ \alpha=\frac{\gamma}{\mathbf p^TA\mathbf p},\qquad \mathbf x_{\mathrm{new}}=\mathbf x+\alpha\mathbf p,\qquad \mathbf r_{\mathrm{new}}=\mathbf r-\alpha A\mathbf p, $$

then $\beta=(\mathbf r_{\mathrm{new}}^T\mathbf r_{\mathrm{new}})/\gamma$ and $\mathbf p_{\mathrm{new}}=\mathbf r_{\mathrm{new}}+\beta\mathbf p$. Under SPD assumptions, directions are mutually $A$-conjugate in exact arithmetic. The code tests the true residual after changing $\mathbf x$, while keeping the recursive residual for the next CG direction. A nonpositive curvature $\mathbf p^TA\mathbf p$, or an unusable zero residual dot product while the true residual still fails, produces breakdown.

5. For GMRES, build an orthonormal Krylov basis and solve a small least-squares problem. At a cycle's base vector $\mathbf x_b$, let $\mathbf r_b=\mathbf b-A\mathbf x_b$, $\beta=\lVert\mathbf r_b\rVert_2$, and $\mathbf v_1=\mathbf r_b/\beta$. Arnoldi orthogonalizes $A\mathbf v_j$ against the existing basis to obtain

$$ AV_j=V_{j+1}\overline H_j, $$

where $\overline H_j$ is upper Hessenberg. GMRES chooses

$$ \mathbf z_j=\operatorname{arg\,min}_{\mathbf z} \lVert\beta\mathbf e_1-\overline H_j\mathbf z\rVert_2, \qquad \mathbf x_j=\mathbf x_b+V_j\mathbf z_j. $$

The implementations use two passes of modified Gram–Schmidt for orthogonalization. Previously computed Givens rotations are applied to each new Hessenberg column, followed by a new rotation that eliminates its last subdiagonal entry. The rotated small system is upper triangular, so back-substitution yields the coefficient vector without forming normal equations.

The triangular solve explicitly checks its pivots against $10^{-14}$ times the largest magnitude in its current upper triangle. This is a guarded QR-based small solve, not Gaussian elimination on the original matrix. No row pivoting is applied to $A$: arbitrary row changes can alter symmetry, stationarity, and the intended Krylov problem. A diagonal problem in a stationary method should prompt a deliberate reordering or another method, rather than an undocumented permutation.

6. Verify the candidate and handle restart or breakdown. GMRES forms a candidate and recomputes $\mathbf b-A\mathbf x$ at every completed inner step, regardless of the rotated residual estimate. If Arnoldi's next vector has norm at most $10^{-14}$ times the original $\lVert A\mathbf v_j\rVert$, the code treats the basis as closed. If the resulting candidate passes the true residual check, this is a successful termination, often called happy breakdown. If it fails, return breakdown. A zero rotation or rejected triangular pivot also returns breakdown with the latest completed candidate.

After $m$ inner steps without convergence, restart at the latest candidate with a new residual and basis. The global iteration budget counts completed candidate updates across all cycles. A breakdown attempt can perform work before producing another completed update, so iteration count is not an exact count of matrix-vector products.

flowchart TD
    A["Validate system, controls, and method requirements"] --> B["Compute true initial residual"]
    B --> C{"Residual within threshold?"}
    C -->|Yes| R["Return converged"]
    C -->|No| D{"Update budget remains?"}
    D -->|No| L["Return limit"]
    D -->|Yes| E["Build sweep, CG direction, or GMRES candidate"]
    E --> F{"Required diagonal, curvature, or small solve usable?"}
    F -->|No| X["Return breakdown or input error"]
    F -->|Yes| G["Complete update; recompute true residual"]
    G --> H{"Residual within threshold?"}
    H -->|Yes| R
    H -->|No| I{"GMRES basis closed?"}
    I -->|Yes| X
    I -->|No| J["Continue recurrence or restart GMRES when due"]
    J --> D

Invalid inputs and checked nonfinite arithmetic raise an exception. Ordinary budget exhaustion returns limit; it is not labeled convergence. The code stores one true residual norm per completed update, including the initial norm, making history length equal to iterations + 1.

Worked Example with a Full Trace

Solve $A\mathbf x=\mathbf b$ for $A=\begin{bmatrix}4&1\\1&3\end{bmatrix}$, $\mathbf b=[1,2]^T$, and $\mathbf x_0=[0,0]^T$. The exact solution is $[1/11,7/11]^T$. For the stationary-method trace, impose a budget of three sweeps with a very tight residual threshold; the complete budgeted runs are:

Both capped runs return limit, because their residuals remain above the tight target. On the first Gauss–Seidel sweep, the second row uses the newly computed $1/4$, rather than Jacobi's old zero. SOR with $\omega=1.1=11/10$ gives first-sweep values $[11/40,253/400]$; its first residual is $[-293/400,-69/400]$. This particular relaxation initially overshoots and is not guaranteed to beat Gauss–Seidel.

CG reaches the exact solution in two steps in exact arithmetic. Its complete recurrence trace is:

Between the two steps, $\beta=(5/16)/5=1/16$. The CG residual norms are $\sqrt5$, $\sqrt5/4$, and $0$.

For GMRES with restart length $m=2$, the complete Arnoldi construction is especially small:

1. Start $\beta=\sqrt5$ and $\mathbf v_1=[1,2]^T/\sqrt5$.
2. $A\mathbf v_1=[6,7]^T/\sqrt5$. Its projection onto $\mathbf v_1$ is $h_{11}=4$. The remaining vector is $[2,-1]^T/\sqrt5$, with norm $h_{21}=1$, giving $\mathbf v_2$.
3. $A\mathbf v_2=[7,-1]^T/\sqrt5$. Its basis coefficients are $h_{12}=1$ and $h_{22}=3$; the remainder is zero, so $h_{32}=0$.

Thus $\overline H_2=\begin{bmatrix}4&1\\1&3\\0&0\end{bmatrix}$. The first GMRES coefficient is $4\sqrt5/17$, giving $\mathbf x_1=[4/17,8/17]^T$ and residual $[-7/17,6/17]^T$. At the second step the least-squares solution is $\mathbf z=[3\sqrt5/11,-\sqrt5/11]^T$, giving $\mathbf x_2=[1/11,7/11]^T$ and zero residual. The closed basis is a successful termination after the true check.

The first Givens rotation has cosine $4/\sqrt{17}$ and sine $1/\sqrt{17}$. The final rotated upper triangle is $\begin{bmatrix}\sqrt{17}&7/\sqrt{17}\\0&11/\sqrt{17}\end{bmatrix}$ with right-hand side $[4\sqrt5/\sqrt{17},-\sqrt5/\sqrt{17}]^T$. This exposes the small back-substitution carried out by the programs.

Implementation in Python, C++, and Java

Each panel is a complete program implementing all five methods on dense nested arrays. The common solver takes a method name, $A$, $\mathbf b$, an initial guess, relative and absolute tolerances, a maximum number of completed updates, SOR's relaxation parameter, and GMRES's restart length. Unused method-specific parameters do not affect the other methods.

Returned fields are the latest candidate x, completed iterations, its true residual norm, status, and the residual history. Status is converged, limit, or breakdown. Inputs are not modified. No solution-vector history is stored, keeping stationary and CG workspace small; the scalar residual history uses one extra number per update.

For CG, exact stored symmetry is enforced; the routines do not silently symmetrize an input. Positive definiteness remains a caller assumption. For GMRES, the small QR triangular pivots are explicitly checked, while the original matrix is not row-pivoted. Reordering a stationary problem or applying a preconditioner requires a deliberate extension, rather than inserting a direct-elimination step into these recurrences.

Use the worked matrix and zero starting guess, $\text{rtol}=10^{-10}$, $\text{atol}=10^{-12}$, budget $200$, $\omega=1.1$, and GMRES restart $2$. The driver verifies both the true residual criterion and error relative to the independently known exact solution.

import math
from dataclasses import dataclass

GUARD = 1e-14


def finite(v):
    if not math.isfinite(v):
        raise ValueError("nonfinite input or arithmetic result")
    return v


def dot(a,b):
    return finite(sum(x*y for x,y in zip(a,b)))


def norm(v):
    total = 0.0
    for x in v:
        total = math.hypot(total,x)
    return finite(total)


def mv(a,x):
    return [dot(row,x) for row in a]


def residual(a,b,x):
    return [finite(y-z) for y,z in zip(b,mv(a,x))]


@dataclass
class Result:
    x: list
    iterations: int
    residual: float
    status: str
    history: list


def solve(method,a,b,x0,rtol=1e-10,atol=1e-12,limit=200,omega=1.1,restart=20):
    n = len(b)
    if (n == 0 or len(a) != n or len(x0) != n or any(len(row) != n for row in a)
            or method not in ("jacobi","gauss-seidel","sor","cg","gmres")
            or not isinstance(limit,int) or limit < 0
            or not math.isfinite(rtol) or rtol < 0 or not math.isfinite(atol) or atol <= 0):
        raise ValueError("invalid dimensions, method, tolerances, or budget")
    for row in a:
        for v in row:
            finite(v)
    for v in b+x0:
        finite(v)
    if method == "sor" and (not math.isfinite(omega) or not 0 < omega < 2):
        raise ValueError("SOR requires omega in (0, 2)")
    if method == "gmres" and (not isinstance(restart,int) or restart < 1):
        raise ValueError("positive GMRES restart required")
    stationary = method in ("jacobi","gauss-seidel","sor")
    if stationary:
        for i in range(n):
            if abs(a[i][i]) <= GUARD*max(abs(v) for v in a[i]):
                raise ValueError("zero or relatively small diagonal")
    if method == "cg" and any(a[i][j] != a[j][i] for i in range(n) for j in range(i)):
        raise ValueError("CG requires a symmetric matrix; positive definiteness is a caller assumption")
    x = x0.copy()
    r = residual(a,b,x)
    target = finite(atol+rtol*norm(b))
    history = [norm(r)]
    def finish(status):
        return Result(x.copy(),len(history)-1,history[-1],status,history.copy())
    if history[-1] <= target:
        return finish("converged")
    if stationary:
        for _ in range(limit):
            old, new = x.copy(), x.copy()
            for i in range(n):
                total = 0.0
                for j in range(n):
                    if j != i:
                        value = new[j] if method != "jacobi" and j < i else old[j]
                        total += a[i][j]*value
                candidate = finite((b[i]-total)/a[i][i])
                new[i] = finite((1-omega)*old[i]+omega*candidate) if method == "sor" else candidate
            x = new
            history.append(norm(residual(a,b,x)))
            if history[-1] <= target:
                return finish("converged")
        return finish("limit")
    if method == "cg":
        p, rr = r.copy(), dot(r,r)
        for _ in range(limit):
            ap = mv(a,p)
            curvature = dot(p,ap)
            if rr <= 0 or curvature <= 0:
                return finish("breakdown")
            alpha = finite(rr/curvature)
            x = [finite(v+alpha*d) for v,d in zip(x,p)]
            history.append(norm(residual(a,b,x)))
            if history[-1] <= target:
                return finish("converged")
            r = [finite(v-alpha*w) for v,w in zip(r,ap)]
            next_rr = dot(r,r)
            if next_rr <= 0:
                return finish("breakdown")
            beta = finite(next_rr/rr)
            p = [finite(v+beta*d) for v,d in zip(r,p)]
            rr = next_rr
        return finish("limit")
    # Restarted GMRES: two-pass modified Gram-Schmidt and Givens rotations.
    m = min(restart,n)
    while len(history)-1 < limit:
        base = x.copy()
        r = residual(a,b,base)
        beta = norm(r)
        v = [[z/beta for z in r]]
        h = [[0.0]*m for _ in range(m+1)]
        cs, sn, g = [0.0]*m, [0.0]*m, [beta]+[0.0]*m
        for j in range(min(m,limit-(len(history)-1))):
            w = mv(a,v[j])
            original_norm = norm(w)
            for _ in range(2):
                for i in range(j+1):
                    projection = dot(v[i],w)
                    h[i][j] = finite(h[i][j]+projection)
                    w = [finite(z-projection*q) for z,q in zip(w,v[i])]
            next_norm = norm(w)
            h[j+1][j] = next_norm
            closed = next_norm <= GUARD*original_norm
            if not closed:
                v.append([z/next_norm for z in w])
            for i in range(j):
                top = cs[i]*h[i][j]+sn[i]*h[i+1][j]
                h[i+1][j] = finite(-sn[i]*h[i][j]+cs[i]*h[i+1][j])
                h[i][j] = finite(top)
            rho = finite(math.hypot(h[j][j],h[j+1][j]))
            if rho == 0.0:
                return finish("breakdown")
            cs[j], sn[j] = h[j][j]/rho, h[j+1][j]/rho
            h[j][j], h[j+1][j] = rho, 0.0
            g[j+1], g[j] = finite(-sn[j]*g[j]), finite(cs[j]*g[j])
            # The small upper-triangular problem has explicit relative pivot tests.
            scale = max(abs(h[ii][jj]) for ii in range(j+1) for jj in range(ii,j+1))
            z = [0.0]*(j+1)
            for i in range(j,-1,-1):
                if h[i][i] == 0.0 or abs(h[i][i]) <= GUARD*scale:
                    return finish("breakdown")
                total = g[i]
                for k in range(i+1,j+1):
                    total = finite(total-h[i][k]*z[k])
                z[i] = finite(total/h[i][i])
            x = base.copy()
            for i in range(j+1):
                x = [finite(value+z[i]*q) for value,q in zip(x,v[i])]
            history.append(norm(residual(a,b,x)))
            if history[-1] <= target:
                return finish("converged")
            if closed:
                return finish("breakdown")
    return finish("limit")


def main():
    a, b, x0 = [[4.,1.],[1.,3.]], [1.,2.], [0.,0.]
    exact, ok = [1/11,7/11], True
    for method in ["jacobi","gauss-seidel","sor","cg","gmres"]:
        r = solve(method,a,b,x0,omega=1.1,restart=2)
        good = r.status == "converged" and norm(residual(a,b,r.x)) <= 1e-12+1e-10*norm(b)
        good = good and max(abs(x-y) for x,y in zip(r.x,exact)) <= 1e-9
        ok = ok and good
        print(f"{method}: x = [{r.x[0]:.6f}, {r.x[1]:.6f}], residual check = {'PASS' if good else 'FAIL'}")
    print(f"all methods = {'PASS' if ok else 'FAIL'}")


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <string>
#include <vector>
using Vec=std::vector<double>;using Matrix=std::vector<Vec>;
const double GUARD=1e-14;
double finite_value(double v){if(!std::isfinite(v))throw std::invalid_argument("nonfinite input or arithmetic result");return v;}
double dot(const Vec& a,const Vec& b){double v=0;for(std::size_t i=0;i<a.size();++i)v+=a[i]*b[i];return finite_value(v);}
double norm(const Vec& v){double total=0;for(double x:v)total=std::hypot(total,x);return finite_value(total);}
Vec mv(const Matrix& a,const Vec& x){Vec out;for(const auto& row:a)out.push_back(dot(row,x));return out;}
Vec residual(const Matrix& a,const Vec& b,const Vec& x){Vec v=mv(a,x);for(std::size_t i=0;i<v.size();++i)v[i]=finite_value(b[i]-v[i]);return v;}
struct Result{Vec x;int iterations;double residual;std::string status;Vec history;};
Result solve(const std::string& method,const Matrix& a,const Vec& b,const Vec& x0,
             double rtol=1e-10,double atol=1e-12,int limit=200,double omega=1.1,int restart=20){
    int n=static_cast<int>(b.size());
    bool known=method=="jacobi" || method=="gauss-seidel" || method=="sor" || method=="cg" || method=="gmres";
    if(n==0 || static_cast<int>(a.size())!=n || static_cast<int>(x0.size())!=n || !known || limit<0
       || !std::isfinite(rtol) || rtol<0 || !std::isfinite(atol) || atol<=0)
        throw std::invalid_argument("invalid dimensions, method, tolerances, or budget");
    for(const auto& row:a){if(static_cast<int>(row.size())!=n)throw std::invalid_argument("nonsquare matrix");for(double v:row)finite_value(v);}
    for(double v:b)finite_value(v);
    for(double v:x0)finite_value(v);
    if(method=="sor" && (!std::isfinite(omega) || omega<=0 || omega>=2))throw std::invalid_argument("SOR requires omega in (0, 2)");
    if(method=="gmres" && restart<1)throw std::invalid_argument("positive GMRES restart required");
    bool stationary=method=="jacobi" || method=="gauss-seidel" || method=="sor";
    if(stationary)for(int i=0;i<n;++i){double scale=0;for(double v:a[i])scale=std::max(scale,std::abs(v));if(std::abs(a[i][i])<=GUARD*scale)throw std::invalid_argument("zero or relatively small diagonal");}
    if(method=="cg")for(int i=0;i<n;++i)for(int j=0;j<i;++j)if(a[i][j]!=a[j][i])throw std::invalid_argument("CG requires symmetry; positive definiteness is a caller assumption");
    Vec x=x0,r=residual(a,b,x);double target=finite_value(atol+rtol*norm(b));Vec history={norm(r)};
    auto finish=[&](const std::string& status){return Result{x,static_cast<int>(history.size())-1,history.back(),status,history};};
    if(history.back()<=target)return finish("converged");
    if(stationary){
        for(int step=0;step<limit;++step){
            Vec old=x,next=x;
            for(int i=0;i<n;++i){
                double total=0;
                for(int j=0;j<n;++j)if(j!=i){double value=method!="jacobi" && j<i ? next[j] : old[j];total+=a[i][j]*value;}
                double candidate=finite_value((b[i]-total)/a[i][i]);
                next[i]=method=="sor" ? finite_value((1-omega)*old[i]+omega*candidate) : candidate;
            }
            x=next;history.push_back(norm(residual(a,b,x)));if(history.back()<=target)return finish("converged");
        }
        return finish("limit");
    }
    if(method=="cg"){
        Vec p=r;double rr=dot(r,r);
        for(int step=0;step<limit;++step){
            Vec ap=mv(a,p);double curvature=dot(p,ap);
            if(rr<=0 || curvature<=0)return finish("breakdown");
            double alpha=finite_value(rr/curvature);for(int i=0;i<n;++i)x[i]=finite_value(x[i]+alpha*p[i]);
            history.push_back(norm(residual(a,b,x)));if(history.back()<=target)return finish("converged");
            for(int i=0;i<n;++i)r[i]=finite_value(r[i]-alpha*ap[i]);
            double next_rr=dot(r,r);if(next_rr<=0)return finish("breakdown");
            double beta=finite_value(next_rr/rr);for(int i=0;i<n;++i)p[i]=finite_value(r[i]+beta*p[i]);rr=next_rr;
        }
        return finish("limit");
    }
    // Restarted GMRES: two-pass modified Gram-Schmidt and Givens rotations.
    int m=std::min(restart,n);
    while(static_cast<int>(history.size())-1<limit){
        Vec base=x;r=residual(a,b,base);double beta=norm(r);Vec first=r;for(double& value:first)value/=beta;
        Matrix v={first},h(m+1,Vec(m));Vec cs(m),sn(m),g(m+1);g[0]=beta;
        int steps=std::min(m,limit-(static_cast<int>(history.size())-1));
        for(int j=0;j<steps;++j){
            Vec w=mv(a,v[j]);double original_norm=norm(w);
            for(int pass=0;pass<2;++pass)for(int i=0;i<=j;++i){
                double projection=dot(v[i],w);h[i][j]=finite_value(h[i][j]+projection);
                for(int k=0;k<n;++k)w[k]=finite_value(w[k]-projection*v[i][k]);
            }
            double next_norm=norm(w);h[j+1][j]=next_norm;bool closed=next_norm<=GUARD*original_norm;
            if(!closed){Vec next=w;for(double& value:next)value/=next_norm;v.push_back(next);}
            for(int i=0;i<j;++i){
                double top=cs[i]*h[i][j]+sn[i]*h[i+1][j];
                h[i+1][j]=finite_value(-sn[i]*h[i][j]+cs[i]*h[i+1][j]);h[i][j]=finite_value(top);
            }
            double rho=finite_value(std::hypot(h[j][j],h[j+1][j]));if(rho==0)return finish("breakdown");
            cs[j]=h[j][j]/rho;sn[j]=h[j+1][j]/rho;h[j][j]=rho;h[j+1][j]=0;
            g[j+1]=finite_value(-sn[j]*g[j]);g[j]=finite_value(cs[j]*g[j]);
            // Explicit relative pivot tests in the small upper-triangular solve.
            double scale=0;for(int ii=0;ii<=j;++ii)for(int jj=ii;jj<=j;++jj)scale=std::max(scale,std::abs(h[ii][jj]));
            Vec z(j+1);
            for(int i=j;i>=0;--i){
                if(h[i][i]==0 || std::abs(h[i][i])<=GUARD*scale)return finish("breakdown");
                double total=g[i];for(int k=i+1;k<=j;++k)total=finite_value(total-h[i][k]*z[k]);z[i]=finite_value(total/h[i][i]);
            }
            x=base;for(int i=0;i<=j;++i)for(int k=0;k<n;++k)x[k]=finite_value(x[k]+z[i]*v[i][k]);
            history.push_back(norm(residual(a,b,x)));if(history.back()<=target)return finish("converged");
            if(closed)return finish("breakdown");
        }
    }
    return finish("limit");
}
int main(){
    Matrix a={ {4.,1.}, {1.,3.} };Vec b={1.,2.},x0={0.,0.},exact={1./11,7./11};bool ok=true;
    std::cout<<std::fixed<<std::setprecision(6);
    for(const auto& method:std::vector<std::string>{"jacobi","gauss-seidel","sor","cg","gmres"}){
        auto r=solve(method,a,b,x0,1e-10,1e-12,200,1.1,2);
        bool good=r.status=="converged" && norm(residual(a,b,r.x))<=1e-12+1e-10*norm(b);
        good=good && std::max(std::abs(r.x[0]-exact[0]),std::abs(r.x[1]-exact[1]))<=1e-9;ok=ok && good;
        std::cout<<method<<": x = ["<<r.x[0]<<", "<<r.x[1]<<"], residual check = "<<(good?"PASS":"FAIL")<<'\n';
    }
    std::cout<<"all methods = "<<(ok?"PASS":"FAIL")<<'\n';return 0;
}

import java.util.ArrayList;
import java.util.Locale;
public class Main {
    static final double GUARD=1e-14;
    static double finite(double v){if(!Double.isFinite(v))throw new IllegalArgumentException("nonfinite input or arithmetic result");return v;}
    static double dot(double[] a,double[] b){double value=0;for(int i=0;i<a.length;++i)value+=a[i]*b[i];return finite(value);}
    static double norm(double[] v){double total=0;for(double x:v)total=Math.hypot(total,x);return finite(total);}
    static double[] mv(double[][] a,double[] x){double[] out=new double[a.length];for(int i=0;i<a.length;++i)out[i]=dot(a[i],x);return out;}
    static double[] residual(double[][] a,double[] b,double[] x){double[] out=mv(a,x);for(int i=0;i<out.length;++i)out[i]=finite(b[i]-out[i]);return out;}
    static class Result {
        double[] x;int iterations;double residual;String status;ArrayList<Double> history;
        Result(double[] x,ArrayList<Double> history,String status){this.x=x.clone();this.history=new ArrayList<>(history);this.iterations=history.size()-1;this.residual=history.get(history.size()-1);this.status=status;}
    }
    static Result solve(String method,double[][] a,double[] b,double[] x0,double rtol,double atol,int limit,double omega,int restart){
        int n=b.length;boolean known=method.equals("jacobi") || method.equals("gauss-seidel") || method.equals("sor") || method.equals("cg") || method.equals("gmres");
        if(n==0 || a.length!=n || x0.length!=n || !known || limit<0 || !Double.isFinite(rtol) || rtol<0 || !Double.isFinite(atol) || atol<=0)
            throw new IllegalArgumentException("invalid dimensions, method, tolerances, or budget");
        for(double[] row:a){if(row.length!=n)throw new IllegalArgumentException("nonsquare matrix");for(double v:row)finite(v);}
        for(double v:b)finite(v);for(double v:x0)finite(v);
        if(method.equals("sor") && (!Double.isFinite(omega) || omega<=0 || omega>=2))throw new IllegalArgumentException("SOR requires omega in (0, 2)");
        if(method.equals("gmres") && restart<1)throw new IllegalArgumentException("positive GMRES restart required");
        boolean stationary=method.equals("jacobi") || method.equals("gauss-seidel") || method.equals("sor");
        if(stationary)for(int i=0;i<n;++i){double scale=0;for(double v:a[i])scale=Math.max(scale,Math.abs(v));if(Math.abs(a[i][i])<=GUARD*scale)throw new IllegalArgumentException("zero or relatively small diagonal");}
        if(method.equals("cg"))for(int i=0;i<n;++i)for(int j=0;j<i;++j)if(a[i][j]!=a[j][i])throw new IllegalArgumentException("CG requires symmetry; positive definiteness is a caller assumption");
        double[] x=x0.clone(),r=residual(a,b,x);double target=finite(atol+rtol*norm(b));ArrayList<Double> history=new ArrayList<>();history.add(norm(r));
        if(history.get(0)<=target)return new Result(x,history,"converged");
        if(stationary){
            for(int step=0;step<limit;++step){
                double[] old=x.clone(),next=x.clone();
                for(int i=0;i<n;++i){
                    double total=0;for(int j=0;j<n;++j)if(j!=i){double value=!method.equals("jacobi") && j<i ? next[j] : old[j];total+=a[i][j]*value;}
                    double candidate=finite((b[i]-total)/a[i][i]);next[i]=method.equals("sor") ? finite((1-omega)*old[i]+omega*candidate) : candidate;
                }
                x=next;history.add(norm(residual(a,b,x)));if(history.get(history.size()-1)<=target)return new Result(x,history,"converged");
            }
            return new Result(x,history,"limit");
        }
        if(method.equals("cg")){
            double[] p=r.clone();double rr=dot(r,r);
            for(int step=0;step<limit;++step){
                double[] ap=mv(a,p);double curvature=dot(p,ap);if(rr<=0 || curvature<=0)return new Result(x,history,"breakdown");
                double alpha=finite(rr/curvature);for(int i=0;i<n;++i)x[i]=finite(x[i]+alpha*p[i]);
                history.add(norm(residual(a,b,x)));if(history.get(history.size()-1)<=target)return new Result(x,history,"converged");
                for(int i=0;i<n;++i)r[i]=finite(r[i]-alpha*ap[i]);double nextRr=dot(r,r);if(nextRr<=0)return new Result(x,history,"breakdown");
                double beta=finite(nextRr/rr);for(int i=0;i<n;++i)p[i]=finite(r[i]+beta*p[i]);rr=nextRr;
            }
            return new Result(x,history,"limit");
        }
        // Restarted GMRES: two-pass modified Gram-Schmidt and Givens rotations.
        int m=Math.min(restart,n);
        while(history.size()-1<limit){
            double[] base=x.clone();r=residual(a,b,base);double beta=norm(r);double[] first=r.clone();for(int i=0;i<n;++i)first[i]/=beta;
            ArrayList<double[]> v=new ArrayList<>();v.add(first);double[][] h=new double[m+1][m];double[] cs=new double[m],sn=new double[m],g=new double[m+1];g[0]=beta;
            int steps=Math.min(m,limit-(history.size()-1));
            for(int j=0;j<steps;++j){
                double[] w=mv(a,v.get(j));double originalNorm=norm(w);
                for(int pass=0;pass<2;++pass)for(int i=0;i<=j;++i){
                    double projection=dot(v.get(i),w);h[i][j]=finite(h[i][j]+projection);for(int k=0;k<n;++k)w[k]=finite(w[k]-projection*v.get(i)[k]);
                }
                double nextNorm=norm(w);h[j+1][j]=nextNorm;boolean closed=nextNorm<=GUARD*originalNorm;
                if(!closed){double[] next=w.clone();for(int i=0;i<n;++i)next[i]/=nextNorm;v.add(next);}
                for(int i=0;i<j;++i){
                    double top=cs[i]*h[i][j]+sn[i]*h[i+1][j];h[i+1][j]=finite(-sn[i]*h[i][j]+cs[i]*h[i+1][j]);h[i][j]=finite(top);
                }
                double rho=finite(Math.hypot(h[j][j],h[j+1][j]));if(rho==0)return new Result(x,history,"breakdown");
                cs[j]=h[j][j]/rho;sn[j]=h[j+1][j]/rho;h[j][j]=rho;h[j+1][j]=0;
                g[j+1]=finite(-sn[j]*g[j]);g[j]=finite(cs[j]*g[j]);
                // Explicit relative pivot tests in the small upper-triangular solve.
                double scale=0;for(int ii=0;ii<=j;++ii)for(int jj=ii;jj<=j;++jj)scale=Math.max(scale,Math.abs(h[ii][jj]));double[] z=new double[j+1];
                for(int i=j;i>=0;--i){
                    if(h[i][i]==0 || Math.abs(h[i][i])<=GUARD*scale)return new Result(x,history,"breakdown");
                    double total=g[i];for(int k=i+1;k<=j;++k)total=finite(total-h[i][k]*z[k]);z[i]=finite(total/h[i][i]);
                }
                x=base.clone();for(int i=0;i<=j;++i)for(int k=0;k<n;++k)x[k]=finite(x[k]+z[i]*v.get(i)[k]);
                history.add(norm(residual(a,b,x)));if(history.get(history.size()-1)<=target)return new Result(x,history,"converged");
                if(closed)return new Result(x,history,"breakdown");
            }
        }
        return new Result(x,history,"limit");
    }
    public static void main(String[] args){
        Locale.setDefault(Locale.ROOT);double[][] a={ {4.,1.}, {1.,3.} };double[] b={1.,2.},x0={0.,0.},exact={1./11,7./11};boolean ok=true;
        for(String method:new String[]{"jacobi","gauss-seidel","sor","cg","gmres"}){
            Result r=solve(method,a,b,x0,1e-10,1e-12,200,1.1,2);boolean good=r.status.equals("converged") && norm(residual(a,b,r.x))<=1e-12+1e-10*norm(b);
            good=good && Math.max(Math.abs(r.x[0]-exact[0]),Math.abs(r.x[1]-exact[1]))<=1e-9;ok=ok && good;
            System.out.printf("%s: x = [%.6f, %.6f], residual check = %s%n",method,r.x[0],r.x[1],good?"PASS":"FAIL");
        }
        System.out.println("all methods = "+(ok?"PASS":"FAIL"));
    }
}

The identical printed output is:

• jacobi: x = [0.090909, 0.636364], residual check = PASS
• gauss-seidel: x = [0.090909, 0.636364], residual check = PASS
• sor: x = [0.090909, 0.636364], residual check = PASS
• cg: x = [0.090909, 0.636364], residual check = PASS
• gmres: x = [0.090909, 0.636364], residual check = PASS
• all methods = PASS

For this setup the completed update counts are $19$, $10$, $11$, $2$, and $2$, respectively. The printed vectors are rounded to six decimals; the convergence check uses the unrounded values. SOR taking eleven updates here versus Gauss–Seidel's ten illustrates why a fixed relaxation choice should be measured rather than assumed optimal.

Complexity, Correctness, and Invariants

Let $n$ be the dimension, $k$ the completed update count, $K$ its budget, and $m\le n$ the GMRES restart length. The programs store $A$ densely, so a matrix-vector product or sweep costs $\Theta(n^2)$. With a suitable sparse representation, let $q$ denote the number of stored nonzeros; the analogous product or sweep is typically $O(q+n)$. The table lists extra workspace beyond the matrix and input vectors, including the stored residual history.

The $m^2$ term per worst-case inner step includes the repeated small back-substitution; the $nm$ term covers basis orthogonalization and candidate construction. GMRES with no effective restart uses $m=n$ and a growing basis; its cumulative full-cycle orthogonalization cost is $O(nk^2)$ for $k\le n$, plus matrix-vector and small-solve work. Sparse products replace the dense $n^2$ product term by approximately $q+n$, but the supplied loops themselves remain dense.

All methods must validate the dense matrix and form an initial residual, so even an accepted initial guess costs $O(n^2)$ here. For a given budget, worst-case work is obtained by substituting $k=K$, unless a failure exits sooner. There is no matrix-independent average iteration count, and stationary methods or restarted GMRES have no unconditional iteration count to a target tolerance. Omitting the scalar history removes its $O(k)$ storage term.

Stationary convergence. In exact arithmetic, error follows $\mathbf e_{k+1}=T\mathbf e_k$. Convergence to the solution for every starting guess is equivalent to $\rho(T)\lt1$. For the splitting convention used here,

$$ T_J=-D^{-1}(L+U),\qquad T_{GS}=-(D+L)^{-1}U,\qquad T_{SOR}=(D+\omega L)^{-1}[(1-\omega)D-\omega U]. $$

Diagonal dominance and SPD results are useful sufficient conditions for the appropriate methods, not interchangeable guarantees for all stationary schemes.

CG's invariant and bound. Under SPD assumptions in exact arithmetic, $\mathbf x_k$ lies in $\mathbf x_0+\mathcal K_k(A,\mathbf r_0)$ and minimizes the error in the $A$-norm there. Search directions are $A$-conjugate and residuals are mutually orthogonal. If $\kappa=\kappa_2(A)$, the classical bound is

$$ \lVert\mathbf e_k\rVert_A\le 2\left(\frac{\sqrt\kappa-1}{\sqrt\kappa+1}\right)^k\lVert\mathbf e_0\rVert_A. $$

Finite precision can weaken orthogonality and separate the recursive residual from the true residual. The extra true-residual calculation prevents declaring success merely because the recurrence reports a misleadingly small value; it does not restore every exact-arithmetic invariant.

GMRES's invariant. Arnoldi constructs an orthonormal basis of the current Krylov space and a small Hessenberg representation of its action under $A$. Minimizing the small least-squares residual is equivalent to minimizing the full Euclidean residual in that affine search space. For a nonsingular $n$-dimensional system, full GMRES reaches a solution in at most $n$ steps in exact arithmetic. Restart discards basis information and loses that finite-step guarantee. Reorthogonalization improves numerical basis quality, but the final original-system residual remains the acceptance check.

Variants, Pitfalls, and Optimizations

SPD alone does not guarantee Jacobi convergence. Consider the $3\times3$ matrix with diagonal $1$ and every off-diagonal entry $0.6$, with $\mathbf b=[2.2,2.2,2.2]^T$ and zero starting guess. Its eigenvalues are $2.2,0.4,0.4$, so it is SPD. Jacobi's iteration matrix has an eigenvalue $-1.2$, and all components follow $u_{k+1}=2.2-1.2u_k$: $0,2.2,-0.44,2.728,-1.0736,\ldots$. Errors grow in magnitude. Gauss–Seidel converges on this SPD system, while CG solves this particular eigenvector-aligned right-hand side in one step in exact arithmetic.

Restart can cause GMRES stagnation. For $A=\begin{bmatrix}0&-1\\1&0\end{bmatrix}$, $\mathbf b=[1,0]^T$, and zero start, GMRES with restart length $1$ repeatedly minimizes along $[1,0]^T$. Its best coefficient is zero, so the iterate stays zero even though the matrix is nonsingular. Full GMRES with restart length $2$ reaches $[0,-1]^T$. Increase restart length or choose a useful preconditioner based on evidence, while accounting for memory and orthogonalization cost.

Choose a compatible preconditioner. Diagonal scaling is cheap; incomplete factorizations, block methods, algebraic multigrid, and domain decomposition can be more effective for suitable systems. Standard preconditioned CG needs an SPD-compatible preconditioner. GMRES permits more general left or right preconditioning, but its internally monitored preconditioned residual may differ from the original residual. Flexible GMRES is appropriate when the preconditioner changes between iterations; ordinary fixed-preconditioner formulas should not be assumed valid for that case.

Respect the matrix class. Nonsymmetric or indefinite systems are not justified CG inputs merely because they sometimes run without breakdown. MINRES can be useful for symmetric indefinite problems; GMRES handles general nonsymmetric systems. Forming $A^TA$ to force SPD squares the two-norm condition number when $A$ has full column rank, so least-squares methods such as LSQR or LSMR can be preferable to explicitly forming normal equations.

Do not hide numerical failure. A zero diagonal can defeat a stationary update even when a system is nonsingular. A tiny update can reflect stagnation rather than a small residual. A nearly closed Arnoldi basis can result from finite-precision dependence; verify the candidate instead of assuming happy breakdown. The code's relative guard is a numerical policy, so very ill-conditioned valid problems may return breakdown. Overflow in a diverging iteration raises an exception rather than reporting a solution.

Exploit structure deliberately. Jacobi updates can be parallelized readily, whereas Gauss–Seidel has row dependencies; colored or block variants change the implementation and convergence behavior. SOR parameter tuning depends on the problem and ordering. Sparse CSR storage avoids processing zeros, and matrix-free Krylov methods need only an operator product. For many right-hand sides with the same matrix, reusing a direct factorization or a preconditioner may be more efficient than repeated unprepared solves.

Interpret residuals with scale and conditioning. Compare the achieved true residual with the requested threshold, and consider a normwise backward-error measure such as $\lVert\mathbf r\rVert/(\lVert A\rVert\lVert\mathbf x\rVert+\lVert\mathbf b\rVert)$ when meaningful. A small backward error can coexist with a larger forward error for an ill-conditioned matrix. The printed sample's known solution supplies a separate forward-error check that a general application may not have.

Practice Problems and Real-World Use

Start with exact small systems, then vary one structural property at a time. Record completed updates, matrix-vector work where measured, true residual histories, and solution error when an independent reference is known. Comparing iteration counts alone hides GMRES's growing basis cost and repeated residual verification.

For further reading, the authors' Templates for the Solution of Linear Systems develops stationary methods, Krylov methods, preconditioners, and stopping criteria. SciPy's CG documentation and GMRES documentation state their matrix assumptions and API-specific convergence controls. Compare restart, iteration-budget, and residual conventions before treating two implementations as identical.