Solving linear systems Ax = b is at the heart of numerical computing — from fitting models to simulating physics. Equana offers three solver backends, each suited to different problem sizes and matrix types:

1. NumWasm — dense direct solvers (LU, QR, Cholesky, least-squares) for small-to-moderate systems
2. SciWasm — sparse iterative solvers (GMRES, BiCGSTAB, PCG, LSQR) for large sparse systems
3. mfemwasm — object-oriented FEM solvers with preconditioner support for finite element problems

This tutorial walks through each backend, explains when to use which solver, and provides runnable examples. Run the cells in order — variables persist across the notebook.

For sparse matrix creation and manipulation, see the Sparse Matrices tutorial. For dense matrix basics, see the Working with Matrices tutorial.

Choosing the Right Solver

The right solver depends on your matrix properties and problem size:

Direct vs. iterative: Direct solvers (LU, QR, Cholesky) compute an exact answer in a fixed number of operations — O(n³) for dense matrices. Iterative solvers approximate the solution through successive refinement and are much faster for large sparse systems where most entries are zero.

Dense Direct Solvers (NumWasm)

Dense solvers compute the exact solution (up to floating-point precision) by factoring the matrix. They are best for systems up to a few thousand unknowns.

Backslash Operator and linsolve

The backslash operator A \ b, linsolve(A, b), and mldivide(A, b) are three equivalent ways to solve Ax = b directly:

LU Decomposition

LU decomposition factors A into a lower triangular matrix L, an upper triangular matrix U, and a permutation matrix P such that P·A = L·U. Once factored, solving for different right-hand sides is cheap (forward and back substitution):

QR Decomposition

QR decomposes A into an orthogonal matrix Q (where Q'·Q = I) and an upper triangular matrix R. QR is numerically more stable than LU and is the basis for least-squares solvers:

Cholesky Factorization

For symmetric positive definite (SPD) matrices, Cholesky computes an upper triangular R such that R'·R = A. It requires half the work of LU and guarantees numerical stability:

Least-Squares Solutions

When a system is overdetermined (more equations than unknowns, m > n), no exact solution exists. lstsq(A, b) finds the x that minimizes ||Ax - b||₂. With multiple outputs, it also returns residuals, rank, and singular values:

Sparse Iterative Solvers (SciWasm)

For large sparse systems, direct solvers become impractical — a 10,000×10,000 dense LU costs ~667 billion operations. Iterative solvers approximate the solution through successive refinement, exploiting the sparsity to perform each iteration in O(nnz) time.

Available Solvers

All solvers support optional tolerance and iteration-limit parameters. With multiple outputs, they return convergence diagnostics:

(x, flag, relres, iter) = solver(A, b, ...)

Setting Up a Test Problem

We build a classic tridiagonal Laplacian — symmetric positive definite and well-suited for benchmarking:

GMRES — General Systems

GMRES (Generalized Minimum Residual) works for any square nonsingular matrix. The optional restart parameter limits memory usage by restarting the Krylov subspace after restart iterations:

BiCGSTAB — Nonsymmetric Systems

BiCGSTAB (Bi-Conjugate Gradient Stabilized) is another choice for general nonsymmetric systems. It uses less memory than GMRES (no restart needed) and often has smoother convergence:

PCG — Symmetric Positive Definite Systems

Preconditioned Conjugate Gradient is the most efficient solver when the matrix is symmetric and positive definite (SPD). It uses half the work per iteration of BiCGSTAB and is guaranteed to converge for SPD matrices. Many physical problems — diffusion, elasticity, Poisson — produce SPD systems:

LSQR — Least-Squares Problems

lsqr solves overdetermined or underdetermined sparse systems in the least-squares sense (minimizes ||Ax - b||₂). It works on rectangular sparse matrices where direct methods and square solvers cannot be applied:

CGS — Conjugate Gradient Squared

CGS is another Krylov method for general nonsymmetric systems. It converges faster than BiCGSTAB in some cases but can have irregular convergence behavior:

Sparse Direct Solvers

For moderate-sized sparse systems where iterative convergence is problematic, sparse direct solvers factorize the matrix. On desktop, Equana uses SuiteSparse (UMFPACK/CHOLMOD) for native-speed factorization:

MFEM Solvers

For finite element problems, MFEM provides an object-oriented solver framework. You create a solver object, configure tolerances and iteration limits, optionally attach a preconditioner, then solve. The solver operates on the assembled system matrix and vectors.

Setting Up a FEM Test Problem

We build a small 2D Poisson problem (-∇²u = 1 with u = 0 on boundary) as a test system:

CG Solver (Object-Oriented API)

Create a CG solver, configure it, set the operator (system matrix), and call mult(b, x) to solve. After solving, query iteration count and residual:

GMRES Solver

For nonsymmetric FEM systems (e.g., advection-diffusion), use GMRES. The kdim parameter controls the Krylov subspace dimension (restart value):

Available MFEM Solver and Preconditioner Classes

Iterative Solvers:

Direct Solvers:

Preconditioners:

Solver Configuration Methods:

Quick Reference

Dense Solvers (NumWasm)

Sparse Iterative Solvers (SciWasm)

Sparse Direct Solvers

MFEM Solvers

Summary

This tutorial covered all three solver backends in Equana:

• Dense direct solvers (NumWasm) — linsolve, mldivide, the backslash operator A \ b, plus LU, QR, and Cholesky factorizations, and lstsq for overdetermined systems
• Sparse iterative solvers (SciWasm) — gmres, bicgstab, pcg, lsqr, cgs, bicg, minres, and tfqmr with convergence diagnostics (flag, relres, iter)
• Sparse direct solvers — spsolve, spchol, splu for moderate-sized sparse systems
• MFEM solvers — object-oriented CG, GMRES, and other Krylov solvers with configurable tolerances, preconditioners, and iteration limits for finite element systems

Next Steps

• Explore the Sparse Matrices tutorial for creating and manipulating sparse matrices
• Try the Finite Element Method tutorial for a complete FEM workflow with solvers
• See the Linear Elasticity tutorial for solving structural problems with CG