Edit or run this notebook
Scientific Computing TU Berlin Winter 2021/22 © Jürgen Fuhrmann
Notebook 18
3.0 ms
4.9 s

Practical iterative methods

4.9 ms

Incomplete LU (ILU) preconditioning

4.6 ms

Idea (Varga, Buleev, 1960 : derive a preconditioner not from an additive decomposition but from the LU factorization.

  • LU factorization has large fill-in. For a preconditioner, just limit the fill-in to a fixed pattern.

  • Apply the standard LU factorization method, but calculate only a part of the entries, e.g. only those which are larger than a certain threshold value, or only those which correspond to certain predefined pattern.

  • Result: incomplete LU factors L, U, remainder R: A=LUR

  • What about zero pivots which prevent such an algoritm from being computable ?

11.0 ms

Theorem (Saad, Th. 10.2): If A is an M-Matrix, then the algorithm to compute the incomplete LU factorization with a given pattern is stable, i.e. does not detriorate due to zero pivots (main diagonal elements) Moreover, A=LUR=MN where M=LU and N=R is a regular splitting.

86.7 ms

Discussion

  • Generally better convergence properties than Jacobi, Gauss-Seidel, though we mostly cannot apply the comparison theorem for regular splittings

  • Block variants are possible

  • ILU Variants:

    • ILUM: ("modified"): add ignored off-diagonal entries to main diagonal

    • ILUT: ("threshold"): zero pattern calculated dynamically based on drop tolerance

    • ILU0: Drop all fill-in

    • Incomplete Cholesky: symmetric variant of ILU

  • Dependence on ordering

  • Can be parallelized using graph coloring

  • Not much theory: experiment for particular systems and see if it works well

  • I recommend it as the default initial guess for a sensible preconditioner

5.0 ms

Further approaches to preconditioning

These are based on ideas which are best explained and developed with multidimensional PDEs in mind.

  • Multigrid: gives indeed O(N) optimal solver complexity in some situations. This is the holy grail method... I will try to discuss this later in the course.

  • Domain decomposition - based on the idea the subdivision of the computational domain into a number of subdomains and subsequent repeated solution of the smaller subdomain problems

323 μs

Iterative methods in Julia

Julia has some well maintained packages for iterative methods and preconditioning.

439 μs

Random sparse M-Matrices

141 μs

We will test the methods with random sparse M matrices, so we define a function which gives us a random, strictly diagonally dominant M-Matrix which is not necessarily irreducible. For skew=0 it is also symmetric:

179 μs
sprandm (generic function with 1 method)
1.7 ms

Test the method a bit...

150 μs
N
5
22.3 μs
634 ms
x1x2x3x4x5
Float64Float64Float64Float64Float64
1
2.6006
-0.922413
-0.165199
-0.840878
-0.462586
2
-0.966885
1.1229
0.0
-0.937833
-0.0462712
3
-0.566097
0.0
0.165822
0.0
0.0
4
-0.703205
0.0
0.0
2.32171
0.0
5
-0.363993
-0.20023
0.0
-0.542756
0.509066
244 ms

Up to rounding errors, the inverse is nonnegative, as predicted by the theory. There are zero entries because it is not necessarily irreducible. Invertibility is guaranteed by strict diagonal dominance.

143 μs
Ainv
5×5 Matrix{Float64}:
  308.942    308.846   307.78    308.839    308.807
  363.27     364.062   361.903   363.534    363.193
 1054.69    1054.36   1056.75   1054.34    1054.23
   93.5729    93.544    93.221    93.9725    93.532
  463.55     463.762   461.806   464.006    465.343
62.9 μs
93.22100392714542
9.6 μs
ρ_jacobi (generic function with 1 method)
543 μs
0.9994562900651217
79.0 μs
Loading...
ii