Direct against iterative

Factorisation is robust and memory-hungry; iteration is cheap and depends entirely on the preconditioner.

A terminal showing a solver residual trace on a dark screen

The two ways to solve a large sparse system

Every finite element run ends at the same place: a linear system Ku = f, where K is the stiffness matrix, f is the load vector, and u is what you came for. The matrix is large, sparse, and structured by the mesh topology. How you solve it determines whether the run finishes in seconds or days, and whether the answer is trustworthy or subtly wrong.

Two families of solver exist. Direct methods factorise K into triangular factors — classically L and U, or L and Lᵀ for symmetric positive-definite problems — and then solve two triangular systems by forward and back substitution. The factorisation is expensive to compute and to store, but once it exists the system can be solved for any right-hand side at low additional cost. Iterative methods start with a guess and refine it, performing a sequence of matrix-vector products until the residual falls below a tolerance. They never form an explicit inverse, and their memory footprint is a fraction of a direct solver's, but their convergence depends on the eigenvalue distribution of K — which depends on the mesh, the physics, and what you do to help them.

The choice is not aesthetic. It follows from problem size, conditioning, how many right-hand sides you need, and whether the physics produces a matrix that iterative methods can handle without heroic intervention.

A convergence plot printed and pinned to a partition

What factorisation actually costs

For a direct solver operating on a two-dimensional problem, fill-in — the nonzero entries that appear in the factors L and U but not in K — grows roughly as O(n log n) in the number of unknowns n, and the operation count is O(n^(3/2)). In three dimensions those numbers become O(n^(4/3)) and O(n²) respectively. These are not constants someone invented to warn you off: they fall out of the nested-dissection reordering theory developed by Alan George in 1973, which is documented in detail in the SIAM literature. For a three-dimensional mesh with a million degrees of freedom, the fill-in alone can exhaust the RAM of a workstation.

Reordering algorithms — nested dissection, approximate minimum degree — reduce fill-in dramatically by reordering the equations before factorising, but cannot eliminate the fundamental scaling. Codes such as MUMPS, PARDISO and SuperLU implement sparse direct solvers with sophisticated reordering, and they are the default choice in many commercial FE packages for problems up to roughly a few hundred thousand to a million unknowns. Beyond that, memory becomes the hard limit.

How you solve it determines whether the run finishes in seconds or days, and whether the answer is trustworthy or subtly wrong.

The robustness of a direct solver is its defining virtue. Given a non-singular K with sufficient floating-point precision, it will produce a solution. It does not converge — it either succeeds or it runs out of memory and time. For ill-conditioned problems where iterative methods struggle to converge at all, a direct solver is often the only reliable path. It is also the right choice when you need to solve for many right-hand sides: the factorisation is computed once, and each additional solve costs only the two triangular back-substitutions.

What iteration actually costs — and demands

An iterative solver replaces factorisation with repeated multiplication of K by a vector. The workhorse for symmetric positive-definite problems is the conjugate gradient method, first published by Magnus Hestenes and Eduard Stiefel in 1952. For non-symmetric or indefinite systems — as arise in incompressible flow or certain coupled problems — GMRES and its variants, introduced by Youcef Saad and Martin Schultz in 1986, are standard. Each iteration is cheap: one sparse matrix-vector product, some dot products, vector updates. Memory scales with problem size and the Krylov subspace dimension, not with fill-in.

The catch is convergence. The number of iterations required grows with the spectral condition number of K — the ratio of the largest to the smallest eigenvalue. A mesh of poor element quality drives the condition number up; so does a wide range of material stiffnesses, or a fine mesh without any special treatment. For a raw finite element stiffness matrix, unpreconditioned CG can require thousands of iterations to reach engineering accuracy, and GMRES can stagnate.

An aisle of compute racks in a machine room

The preconditioner M is what makes iterative solvers practical. Instead of solving Ku = f, the solver works on M⁻¹Ku = M⁻¹f (or a symmetric equivalent), where M is chosen so that M⁻¹K has a tightly clustered spectrum. The ideal M is K itself — which would converge in one step but costs exactly as much as a direct solve. The practical M is something cheaper that nonetheless captures the dominant structure of K.

Diagonal (Jacobi) preconditioning — dividing each equation by its diagonal entry — is trivially cheap and almost always insufficient. Incomplete LU (ILU) factorisation computes a sparse approximate factorisation by dropping fill-in below a threshold; it is far more effective but introduces its own complexity and can fail for indefinite problems. For elliptic PDEs on structured meshes, multigrid preconditioners are optimal: they achieve iteration counts that are essentially independent of problem size, because they correct error on multiple length scales simultaneously. Algebraic multigrid (AMG) extends this to unstructured meshes without geometric information and is the basis of solvers such as BoomerAMG, part of the hypre library developed at Lawrence Livermore National Laboratory, whose technical documentation is publicly available ↗.

The difficulty is that no single preconditioner works for all problems. AMG performs well on scalar elliptic problems but struggles with highly anisotropic systems or saddle-point problems. Block preconditioners for mixed formulations — pressure-velocity coupling in incompressible Navier-Stokes, for instance — require problem-specific knowledge to construct. The iterative solver user must understand the physics to choose a preconditioner, and an inappropriate choice can mean non-convergence with no obvious warning beyond a stagnating residual.

Choosing in practice

The decision tree used by most practitioners is roughly as follows. For small to medium problems — up to a few hundred thousand unknowns in 3D — a sparse direct solver is the default. It is robust, its memory cost is acceptable on modern hardware, and it removes uncertainty about convergence. For large 3D problems where memory is the constraint, or for problems run in parallel on many cores, iterative solvers with appropriate preconditioners become necessary. Sparse direct factorisation parallelises poorly beyond a moderate number of cores because fill-in creates dense subblocks that are difficult to distribute; iterative methods, whose inner loop is a sparse matrix-vector product, scale far more naturally.

For time-dependent problems where K changes modestly between steps, a direct solver can re-use a previously computed ordering and update the factorisation at reduced cost — an advantage that partially offsets the memory burden. For nonlinear problems solved by Newton iteration, the system changes at every step, which levels the field: a direct solve per Newton step versus preconditioned iterative solve per Newton step, with the comparison depending on how the preconditioner handles the changing matrix.

Verification of the linear solver — confirming that it actually solves the algebraic system to the required tolerance, separate from whether the discretisation itself is correct — is often overlooked. The residual norm ||Ku - f|| / ||f|| should be checked explicitly, not assumed. An iterative solver that exits on an iteration limit rather than a residual tolerance has not solved the problem; it has stopped trying.