Contents

2021-10-22 Preconditioners

Last time

  • Limitations of stationary iterative methods

  • Preconditioning as a concept

  • Krylov methods (GMRES)

Today

  • Experiments with preconditioners and diagnostics in PETSc

  • Basic preconditioners

  • Incomplete factorization and fill

using Plots
using LinearAlgebra
using SparseArrays

default(linewidth=4)

function advdiff_matrix(n; kappa=1, wind=[0, 0])
    "Advection-diffusion with Dirichlet boundary conditions eliminated"
    h = 2 / (n + 1)
    rows = Vector{Int64}()
    cols = Vector{Int64}()
    vals = Vector{Float64}()
    idx((i, j),) = (i-1)*n + j
    in_domain((i, j),) = 1 <= i <= n && 1 <= j <= n
    stencil_advect = [-wind[1], -wind[2], 0, wind[1], wind[2]] / h
    stencil_diffuse = [-1, -1, 4, -1, -1] * kappa / h^2
    stencil = stencil_advect + stencil_diffuse
    for i in 1:n
        for j in 1:n
            neighbors = [(i-1, j), (i, j-1), (i, j), (i+1, j), (i, j+1)]
            mask = in_domain.(neighbors)
            append!(rows, idx.(repeat([(i,j)], 5))[mask])
            append!(cols, idx.(neighbors)[mask])
            append!(vals, stencil[mask])
        end
    end
    sparse(rows, cols, vals)
end

function my_spy(A)
    cmax = norm(vec(A), Inf)
    s = max(1, ceil(120 / size(A, 1)))
    spy(A, marker=(:square, s), c=:diverging_rainbow_bgymr_45_85_c67_n256, clims=(-cmax, cmax))
end

grad_descent (generic function with 1 method)

Experiments in PETSc

  • PETSc = Portable Extensible Toolkit for Scientific computing

  • ./configure, make

    • Depends on BLAS and LAPACK (default system or package manager)

      • or ./configure --download-f2cblaslapack --download-blis

    • Depends on MPI for parallelism (package manager)

      • or ./configure --download-mpich (or --download-openmpi)

  • docker run -it --rm jedbrown/petsc

  • Lots of examples (mostly C and Fortran, some Python)

  • Experimental bindings in Rust, Julia

  • We’ll use src/snes/tutorials/ex15.c

    • “p-Bratu”: combines p-Laplacian with Bratu nonlinearity

      \[ -\nabla\cdot\big(\kappa(\nabla u) \nabla u\big) - \lambda e^u = f \]

    • Newton solver with Krylov on each iteration

Simple preconditioners

Jacobi -pc_type jacobi

\[ P_{\text{Jacobi}}^{-1} = D^{-1} \]

where \(D\) is the diagonal of \(A\).

Gauss-Seidel -pc_type sor

\[ P_{GS}^{-1} = (L+D)^{-1} \]

where \(L\) is the (strictly) lower triangular part of \(A\). The upper triangular part may be used instead, or a symmetric form

\[ P_{SGS}^{-1} = (L+U)^{-1} A \Big( I - (L+D)^{-1} \Big) . \]

Over-relaxation

-pc_sor_omega 1. is default (Gauss-Seidel)

Run p-Laplacian example

$ cd src/snes/tutorials
$ make ex15
$ ./ex15 -da_refine 2 -dm_view
$ ./ex15 -ksp_monitor -pc_type jacobi
$ ./ex15 -snes_view

Experiments

  • How does iteration count vary under grid refinement?

  • How sensitive is it to parameters

    • p-Laplacian -p \(> 1\) and -epsilon \(> 0\)

    • Bratu -lambda \(< 6.8\)

  • How sensitive to -ksp_gmres_restart?

  • -ksp_monitor_true_residual

Symmetric problems

Lanczos iteration: like GMRES for symmetric problems

If \(A\) is symmetric, then \(H = Q^T A Q\) is also symmetric. Since \(H\) is Hessenberg, this means \(H\) is tridiagonal. Instead of storing \(Q_n\), it is sufficient to store only the last two columns since the iteration satisfies a 3-term recurrence. The analog of GMRES for the symmetric case is called MINRES and is also useful for solving symmetric indefinite problems.

Conjugate Gradients: changing the norm

Instead of minimizing the \(A^T A\) norm of the error, the Conjugate Gradient method minimizes the \(A\) norm of the error. For \(A\) to induce a norm, it must be symmetric positive definite. Jeremy Shewchuck’s guide to CG is an excellent resource.

Incomplete factorization

  • Start factoring like in a sparse direct solver

    • ILU(0) Discard fill outside the sparsity pattern of \(A\)

    • ILU(\(k\)): only allow \(k\) levels of fill

    • ILUT(\(\epsilon\)): only keep fill values larger than \(\epsilon\)

Experiments

  • Try -pc_factor_levels 2

  • Impact on cost (check -log_view)

  • Impact on scalability (-da_refine 6)

Kershaw (1978) matrix

Incomplete Cholesky can break down for SPD matrices.

K = [3 -2 0 2; -2 3 -2 0; 0 -2 3 -2; 2 0 -2 3]
eigvals(K)

4-element Vector{Float64}:
 0.17157287525380904
 0.17157287525381085
 5.828427124746186
 5.82842712474619