An Introduction to Preconditioners Victor Eijkhout

594, March 2005

594 Dongarra/Eijkhout 2005/03/16

Introduction Algebraic preconditioners Incomplete factorization (ILU) preconditioners Block methods Approximations of the inverse Domain partitioning methods Optimal solvers Multigrid

594 Dongarra/Eijkhout 2005/03/16

‘Simple’ preconditioners

Preconditioners constructed from matrix elements: Jacobi: M = DA Gauss-Seidel: M = DA + LA SOR: M = DA + ωLA SSOR: M = (ω −1 DA + LA ) (2ω −1 − 1)DA


−1

(ω −1 DA + UA )

594 Dongarra/Eijkhout 2005/03/16

Classical theory of simple preconditioners

Convergence condition: ρ(I − M −1 A) < 1 Convergence guaranteed only for simple problems: M-matrices Jacobi and G-S: #it ∼ h−2 SOR: #it ∼ h−1 for optimal omega

594 Dongarra/Eijkhout 2005/03/16

Current state of simple preconditioners

Stationary iterative methods are not used: convergence theory too limited Problem with G-S and SOR: nonsymmetry Only Jacobi still used with non-stationary methods, sometimes SSOR but only with ω = 1

594 Dongarra/Eijkhout 2005/03/16

Preconditioners in non-stationary methods

Convergence theorypis incomplete: only bounds known for instance #it ∼ κ(A) Possible criterium κ(M −1 A) < κ(A) either order of magnitude or constant

594 Dongarra/Eijkhout 2005/03/16

Incomplete factorization (ILU) preconditioners

594 Dongarra/Eijkhout 2005/03/16

Inspiration from direct methods

Direct methods: Gaussian elimination A = LU, Ax = b ⇒ x = A−1 b = U −1 (L−1 b) Problem with LU is fill-in: discarding fill-in gives approximate solution Aim: let LU take storage similar to A

594 Dongarra/Eijkhout 2005/03/16

Discarding fill-in

Exact factorization: −1 ∀i,j>k : aij ← aij − aik akk akj .

Approximate factorization: ∀i,j>k : if (i, j) ∈ S

−1 akj . aij ← aij − aik akk

Limit storage by defining S −1 Alternatively: only fill in zero locations if aik akk akj large enough

594 Dongarra/Eijkhout 2005/03/16

Error analysis of ILU Laplacian: 

 4 −1 −1  −1 4 −1 −1     . . . . . . . .   . . . .   A=  −1 4 −1 −1    −1 −1 4 −1 −1   .. .. .. .. .. . . . . . One elimination step:   4 −1 −1 −.25 −1 A =  3.75 −1 −1 4 −1 One more elimination step:  4 −1 A =  3.75 −1 3.73 −1

 −1  −.25 −1 −.066 −.2666 −1 594 Dongarra/Eijkhout 2005/03/16

Limit analysis

Because of the sign pattern, pivots decrease ⇒ inverses increase ⇒ limit pivot wanted

594 Dongarra/Eijkhout 2005/03/16

Limit analysis’

One-dimensional case: 

 2 −1  −1 2 −1   A=  −1 2 −1   .. .. .. . . .

Pivot equation di+1 = 2 − 1/di Limit equation d = 2 − 1/d, solution d = 1 Two-dimensional case, pivot is altered by previous point and previous line: dij = 4 − 1/di−1,j − 1/di,j−1 √ Limit equation d = 4 − 2/d, solution 2 + 2

594 Dongarra/Eijkhout 2005/03/16

Fourier analysis

Pivots converge quickly ⇒ pretend they are constant Matrix and factorization have constant diagonal ⇒ difference operators, eigenfunctions are sines

594 Dongarra/Eijkhout 2005/03/16

Modified ILU

Instead of discarding fill-in, add to the diagonal ( −1 aij ← aij − aik akk akj accept fill −1 aii ← aii − aik akk akj discard fill Physical meaning: conservation of mass Theory: possible reduction κ(M −1 A) = O(h−1 )

594 Dongarra/Eijkhout 2005/03/16

ILU-D

If the only fill-in is on the diagonal, reuse storage of LA and UA : only one extra vector for the preconditioner. ILU-D: only allow fill-in on the diagonal Lemma: ILU-D ≡ ILU(0) if there are no triangles in the matrix graph Proof: homework

594 Dongarra/Eijkhout 2005/03/16

Computational stuff Assume ILU-D, so L = LA , U = UA . Different ways of writing the factorization: M = (D + L)D −1 (D + U) = (I + LD −1 )(D + U) = (D + L)(I + D −1 U) = (I + LD −1 )D(I + D −1 U)

Computational cost? (2) and (3) cheaper if scaled factors are stored; otherwise all the same if we store D or D −1 as needed Storage cost? (1) and (4): two vectors; (2) and (3) only one ease of implementation? 594 Dongarra/Eijkhout 2005/03/16

Computational stuff, cont’d Solving (D + L)x = y or (D + U)x = y simple: X xi ← di−1 (yi − `ij xj ) i = 1...n ji

However, (I + D −1 U)x = y xi ← yi − di−1

X

uij xj

j>i

while (I + LD −1 )x = y xi ← yi −

X

dj−1 `ij xj

j