Multigrid

Deep dive

How it works

In: operator A on a grid hierarchy (or the matrix, for algebraic multigrid), right-hand side f, initial guess, smoothing counts ν₁, ν₂ · Out: u with the residual reduced by a fixed factor per cycle · State: u, f and the residual on every level (coarse operators built once) · Cost: about 50 flops per unknown per V(2,2) cycle in 2D

V-cycle on level ℓ for A_ℓ u = f:

1. Pre-smooth ν₁ times: red-black Gauss-Seidel updates all "red" points (i + j even) from their black neighbours, then all black ones.
2. Residual r = f − A_ℓ u.
3. Restrict r to level ℓ+1 by full weighting: weights 1/4 centre, 1/8 edge neighbours, 1/16 corners.
4. Recurse on A_{ℓ+1} e = r with e = 0 initial guess (or solve exactly on the coarsest level).
5. Prolong e back with bilinear interpolation and add: u += P e.
6. Post-smooth ν₂ times.

Why it works: a smoother reduces error components with wavelength ~ h by a factor μ ≈ 1/4 per sweep (red-black GS on Poisson) but barely touches wavelength ~ L. After restriction, the same smooth component has wavelength ~ 2h relative to the coarse spacing, and is again "high frequency" a few levels down. Full multigrid (FMG) starts from the coarsest grid and interpolates upward, reaching discretization accuracy in about one V-cycle per level.

fine   h   S ─R─┐                      ┌─P─ S
           2h    S ─R─┐          ┌─P─ S
           4h          S ─R─┐┌─P─ S
coarse 8h                  solve

Code

/// Grid of (n+1) x (n+1) points on the unit square, boundary values fixed at 0. n = 2^k.
#[derive(Clone)]
pub struct Grid { pub n: usize, pub v: Vec<f64> }
impl Grid {
    pub fn zeros(n: usize) -> Grid { Grid { n, v: vec![0.0; (n + 1) * (n + 1)] } }
    #[inline] fn at(&self, i: usize, j: usize) -> f64 { self.v[j * (self.n + 1) + i] }
    #[inline] fn set(&mut self, i: usize, j: usize, x: f64) { let n = self.n; self.v[j * (n + 1) + i] = x; }
}

/// Red-black Gauss-Seidel for -Laplace(u) = f with the 5-point stencil.
pub fn smooth(u: &mut Grid, f: &Grid, sweeps: usize) {
    let n = u.n;
    let h2 = 1.0 / (n * n) as f64;
    for _ in 0..sweeps {
        for color in 0..2 {
            for j in 1..n {
                let start = 1 + (j + 1 + color) % 2; // points with (i + j) % 2 == color
                for i in (start..n).step_by(2) {
                    let s = u.at(i - 1, j) + u.at(i + 1, j) + u.at(i, j - 1) + u.at(i, j + 1);
                    u.set(i, j, 0.25 * (s + h2 * f.at(i, j)));
                }
            }
        }
    }
}

pub fn residual(u: &Grid, f: &Grid) -> Grid {
    let n = u.n;
    let inv_h2 = (n * n) as f64;
    let mut r = Grid::zeros(n);
    for j in 1..n {
        for i in 1..n {
            let lap = 4.0 * u.at(i, j) - u.at(i - 1, j) - u.at(i + 1, j) - u.at(i, j - 1) - u.at(i, j + 1);
            r.set(i, j, f.at(i, j) - inv_h2 * lap);
        }
    }
    r
}

/// Full weighting: 1/4 centre, 1/8 edges, 1/16 corners.
pub fn restrict(r: &Grid) -> Grid {
    let nc = r.n / 2;
    let mut c = Grid::zeros(nc);
    for j in 1..nc {
        for i in 1..nc {
            let (fi, fj) = (2 * i, 2 * j);
            let edges = r.at(fi - 1, fj) + r.at(fi + 1, fj) + r.at(fi, fj - 1) + r.at(fi, fj + 1);
            let corners = r.at(fi - 1, fj - 1) + r.at(fi + 1, fj - 1) + r.at(fi - 1, fj + 1) + r.at(fi + 1, fj + 1);
            c.set(i, j, 0.25 * r.at(fi, fj) + 0.125 * edges + 0.0625 * corners);
        }
    }
    c
}

/// Bilinear interpolation of the coarse correction, added into u.
pub fn prolong_add(e: &Grid, u: &mut Grid) {
    let nc = e.n;
    for j in 0..nc {
        for i in 0..nc {
            let (a, b, c, d) = (e.at(i, j), e.at(i + 1, j), e.at(i, j + 1), e.at(i + 1, j + 1));
            let (fi, fj) = (2 * i, 2 * j);
            let add = |u: &mut Grid, x: usize, y: usize, v: f64| { let w = u.at(x, y) + v; u.set(x, y, w); };
            if i > 0 && j > 0 { add(u, fi, fj, a); }
            if j > 0 { add(u, fi + 1, fj, 0.5 * (a + b)); }
            if i > 0 { add(u, fi, fj + 1, 0.5 * (a + c)); }
            add(u, fi + 1, fj + 1, 0.25 * (a + b + c + d));
        }
    }
}

pub fn v_cycle(u: &mut Grid, f: &Grid, nu1: usize, nu2: usize) {
    if u.n <= 2 { smooth(u, f, 1); return; } // one interior point: GS is exact
    smooth(u, f, nu1);
    let rc = restrict(&residual(u, f));
    let mut ec = Grid::zeros(rc.n);
    v_cycle(&mut ec, &rc, nu1, nu2);
    prolong_add(&ec, u);
    smooth(u, f, nu2);
}

pub fn norm(g: &Grid) -> f64 { (g.v.iter().map(|x| x * x).sum::<f64>() / g.v.len() as f64).sqrt() }

/// Solve to a relative residual `rtol`; returns the number of V-cycles.
pub fn solve(u: &mut Grid, f: &Grid, rtol: f64, max_cycles: usize) -> usize {
    let r0 = norm(&residual(u, f)).max(f64::MIN_POSITIVE);
    for c in 1..=max_cycles {
        v_cycle(u, f, 2, 2);
        if norm(&residual(u, f)) <= rtol * r0 { return c; }
    }
    max_cycles
}

Complexity and constants

• Work per V(2,2) cycle in 2D: 4 smoothing sweeps + 1 residual + transfers on the fine level, times 4/3 for all levels: on the order of 50 flops per unknown. 3D factor 8/7.
• Measured with the code above (−Δu = f, V(2,2), rtol 10⁻¹⁰, Apple M4 Pro, one thread): 9 cycles for every n from 64 to 2048; steady residual reduction 0.058-0.060 per cycle; about 9-13 ns per unknown per cycle, so 2047² ≈ 4.2·10⁶ unknowns solve in 0.38 s. V(1,1) gives 0.124 per cycle, V(3,3) 0.044. Discretization error falls 8.0·10⁻⁴ → 2.0·10⁻⁴ → 5.0·10⁻⁵ for n = 32, 64, 128 (second order), so in practice a few cycles (or FMG) suffice.
• Memory: about 4/3 × (u, f, r) on the finest level in 2D.
• Coarsening ratio 2 per dimension; AMG typically coarsens by 2-4× unknowns per level, with operator complexity (total nnz over fine nnz) of 1.2-2 a common target.

Efficiency

• Cost: about 50 flops per unknown per V(2,2) cycle in 2D (4/3 of the finest level). Measured with the code above (Apple M4 Pro, one thread): 9 cycles for every n from 64 to 2048, 9-13 ns per unknown per cycle, 2047² ≈ 4.2·10⁶ unknowns solved in 0.38 s. Memory about 4/3 × (u, f, r) on the finest level. Bound: memory bandwidth of the fine-level sweeps.
• Levers:
  • V(2,2) smoothing — residual factor 0.058-0.060 per cycle, against 0.124 for V(1,1) and 0.044 for V(3,3) (measured above) — 2D Poisson, a good default
  • full multigrid start — discretization accuracy in about one V-cycle per level (page above) — when the discretization error is the target
  • multigrid as the CG preconditioner — nearly mesh-independent iteration counts (CG) — robust production use
  • fused smoothing and residual sweeps with red-black or damped Jacobi — fewer memory passes, parallel — fine levels, GPUs
  • Galerkin coarse operators or AMG — right convergence with jumping coefficients (page above) — variable-coefficient and unstructured problems
• Naive → standard → state of the art:
  • Gauss-Seidel alone · iterations grow as h⁻² with the grid
  • geometric V(2,2) · 9 cycles for any n, O(n) work (measured above)
  • FMG or AMG-preconditioned CG on MPI/GPU (hypre, AmgX) · near-constant iterations at scale

Implementation notes

• Red-black ordering makes Gauss-Seidel parallel (each colour is independent) and is a better smoother for the 5-point stencil than lexicographic order.
• Damped Jacobi (ω = 4/5 for 2D 5-point) is the GPU-friendly alternative; slightly weaker smoothing, embarrassingly parallel.
• Galerkin coarse operator A_c = R A P (used by AMG) handles variable coefficients correctly; rediscretizing on the coarse grid is cheaper but can fail for jumping coefficients.
• As a preconditioner: one V-cycle with symmetric smoothing (forward then backward, or red-black then black-red) keeps M symmetric for PCG.
• Neumann / periodic problems: the coarse problem is singular; fix one value or project out the mean on each level.
• Anisotropy and stretched grids break pointwise smoothers: use line relaxation or semi-coarsening.
• Cache blocking: fuse smoothing sweeps and residual computation; the fine-level sweeps are memory-bound.

Pitfalls

• Grid sizes that are not compatible with coarsening (n not divisible by 2 enough times): the coarse grid misaligns.
• Smoothing on the boundary nodes of Dirichlet problems, or restricting boundary residuals: corrupts the correction.
• Using only 1 smoothing sweep with a weak smoother: convergence factor degrades toward that of the smoother.
• Measuring convergence by the error norm without a reference solution: use the residual norm; converge to discretization error, not machine precision.
• Interpolation order too low for the problem (injection for prolongation): convergence factor near 1. Rule: interpolation order + restriction order > PDE order.

Variants and alternatives

In production

Solver libraries: hypre (BoomerAMG), PETSc (PCGAMG, PCMG), Trilinos ML/MueLu, AmgX (NVIDIA), AMGCL. Benchmarks: HPGMG. Simulation: pressure solvers in large CFD and climate codes; fluid animation in graphics uses multigrid pressure solves (Chentanez and Müller) or MG-PCG. Machine learning: multigrid-style neural architectures (Ke et al.).

Implementations

Sources

• William L. Briggs, Van Emden Henson, Steve F. McCormick, A Multigrid Tutorial, 2nd ed., SIAM, 2000. book
• Achi Brandt, Multi-Level Adaptive Solutions to Boundary-Value Problems, Mathematics of Computation 31(138), 1977. paper
• Ulrich Trottenberg, Cornelius W. Oosterlee, Anton Schüller, Multigrid, Academic Press, 2001.
• John W. Ruge, Klaus Stüben, Algebraic Multigrid, in S. F. McCormick (ed.), Multigrid Methods, SIAM, 1987. chapter
• Van Emden Henson, Ulrike Meier Yang, BoomerAMG: A Parallel Algebraic Multigrid Solver and Preconditioner, Applied Numerical Mathematics 41(1), 2002. paper · code
• Nuttapong Chentanez, Matthias Müller, A Multigrid Fluid Pressure Solver Handling Separating Solid Boundary Conditions, SCA 2011. paper
• Tsung-Wei Ke et al., Multigrid Neural Architectures, CVPR 2017. paper