Building a Numerical Solver

Analogy categoryAnalogies

Scoping a Solver: Direct vs Iterative

Before writing a line of code, decide whether the linear system Ax = b calls for a direct method (Gaussian elimination, LU decomposition via LAPACK's dgesv) or an iterative method (Jacobi, Gauss-Seidel, or Conjugate Gradient). Direct methods give an exact answer up to floating-point error in a predictable number of operations, O(n^3), but become impractical in memory and time once n exceeds roughly 10,000-50,000 for dense systems. Iterative methods exploit sparsity — common in discretized PDEs where each row of A only has a handful of nonzero entries — converging to a chosen tolerance in far fewer operations per iteration, at the cost of needing a convergence criterion and, often, a preconditioner.

🏏

Cricket analogy: Choosing a direct solver is like sending in a specialist to bat out every ball of a chase for a guaranteed, exact result, while an iterative solver is like a chasing team using DLS-style incremental targets, converging toward the win in stages when the exact full chase isn't feasible.

Implementing Conjugate Gradient

Conjugate Gradient (CG) is the workhorse iterative method for symmetric positive-definite systems, common when discretizing elliptic PDEs like the Poisson equation with finite differences or finite elements. Each iteration computes a search direction conjugate to previous directions, guaranteeing convergence to the exact solution in at most n steps in exact arithmetic, though in practice far fewer iterations suffice once the residual drops below a tolerance like 1e-8. The implementation needs only matrix-vector products (never the full matrix inverse), making it ideal for sparse matrices stored in compressed formats where you never materialize the dense A explicitly.

🏏

Cricket analogy: CG's conjugate search directions are like a bowling attack that never repeats the same line and length twice against a set batsman, each new delivery deliberately orthogonal to the last plan so the batsman can't settle, converging on a wicket faster than random bowling.

fortran

module cg_solver_mod
  use, intrinsic :: iso_fortran_env, only: rk => real64
  implicit none
  private
  public :: conjugate_gradient

contains

  subroutine conjugate_gradient(a, b, x, tol, max_iter, iters, converged)
    real(rk), intent(in)    :: a(:,:), b(:)
    real(rk), intent(inout) :: x(:)
    real(rk), intent(in)    :: tol
    integer,  intent(in)    :: max_iter
    integer,  intent(out)   :: iters
    logical,  intent(out)   :: converged

    real(rk), allocatable :: r(:), p(:), ap(:)
    real(rk) :: rs_old, rs_new, alpha, beta
    integer  :: n, k

    n = size(b)
    allocate(r(n), p(n), ap(n))

    r = b - matmul(a, x)
    p = r
    rs_old = dot_product(r, r)
    converged = .false.

    do k = 1, max_iter
      ap = matmul(a, p)
      alpha = rs_old / dot_product(p, ap)
      x = x + alpha * p
      r = r - alpha * ap
      rs_new = dot_product(r, r)

      if (sqrt(rs_new) < tol) then
        converged = .true.
        iters = k
        return
      end if

      beta = rs_new / rs_old
      p = r + beta * p
      rs_old = rs_new
    end do

    iters = max_iter
  end subroutine conjugate_gradient

end module cg_solver_mod

Conjugate Gradient only converges reliably for symmetric positive-definite matrices. Applying it to a non-symmetric or indefinite system either fails to converge or converges to a nonsensical result; use GMRES or BiCGSTAB for non-symmetric systems, or symmetrize the normal equations (A^T A x = A^T b) only as a last resort, since it squares the condition number.

Preconditioning and Convergence

CG's convergence rate depends on the condition number of A: a poorly conditioned matrix (large ratio of largest to smallest eigenvalue) converges painfully slowly, sometimes stalling in floating-point arithmetic before reaching tolerance. A preconditioner M approximates A^-1 cheaply and transforms the system into one with a much tighter eigenvalue spread — even a simple diagonal (Jacobi) preconditioner, dividing each row by its diagonal entry, often cuts iteration counts substantially for diagonally dominant systems, while incomplete Cholesky (IC(0)) preconditioning is a common stronger choice for finite-element stiffness matrices.

🏏

Cricket analogy: A preconditioner is like a pitch report given to a bowler before the match — instead of bowling blind on an unknown surface (ill-conditioned), knowing the pitch tends to seam early lets the bowler adjust length immediately and take wickets (converge) faster.

A practical convergence check combines a relative residual test (norm(r) / norm(b) < tol) with a hard iteration cap (max_iter), and logs the residual history — a solver that silently stops at max_iter without flagging non-convergence is a common source of silently wrong downstream results in production pipelines.

Choose direct solvers for small/dense systems needing exact answers; choose iterative solvers for large sparse systems.
Conjugate Gradient requires a symmetric positive-definite matrix and needs only matrix-vector products, not the full inverse.
Convergence speed depends on the condition number of A; preconditioning tightens the eigenvalue spread to speed it up.
Diagonal (Jacobi) preconditioning is cheap and effective for diagonally dominant systems; incomplete Cholesky is stronger for FEM stiffness matrices.
Never apply CG to a non-symmetric system directly; use GMRES or BiCGSTAB instead.
Always cap iterations with max_iter and explicitly report non-convergence rather than silently returning a stale result.
Store sparse matrices in compressed formats (CSR/CSC) so matrix-vector products avoid ever materializing the dense matrix.