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.
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_modConjugate 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_iterand 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.
Practice what you learned
1. When is a direct solver (e.g., LU decomposition) generally preferable to an iterative solver?
2. What property must a matrix A have for Conjugate Gradient to be guaranteed to converge?
3. What is the primary purpose of a preconditioner in an iterative solver like CG?
4. Why does the example CG implementation only need matrix-vector products (`matmul(a, p)`) rather than computing A^-1?
5. What is a key risk of a solver loop that stops after `max_iter` iterations without explicitly checking a `converged` flag?
Was this page helpful?
You May Also Like
Fortran Best Practices
Practical guidelines for writing clean, portable, and performant modern Fortran code, from type safety to array handling.
Fortran vs Modern Languages
How Fortran compares to C++, Python, Rust, and Julia for numerical and scientific computing workloads.
Fortran Quick Reference
A condensed cheat sheet of core modern Fortran syntax, intrinsics, and conventions for day-to-day coding.
Related Reading
Related Study Notes in Programming
Browse all study notesApache Spark Study Notes
Programming · 30 topics
ProgrammingApache Flink Study Notes
Programming · 30 topics
ProgrammingHadoop Study Notes
Programming · 30 topics
ProgrammingSnowflake Study Notes
Programming · 30 topics
ProgrammingApache Airflow Study Notes
Programming · 30 topics
Programmingdbt (Data Build Tool) Study Notes
Programming · 30 topics