# 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