Skip to content

Commit 9649e7e

Browse files
Merge pull request #551 from j-fu/jf/document-precs
Document `precs`
2 parents d3390ec + bd0d087 commit 9649e7e

File tree

2 files changed

+56
-15
lines changed

2 files changed

+56
-15
lines changed

docs/src/basics/Preconditioners.md

Lines changed: 55 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -4,15 +4,13 @@ Many linear solvers can be accelerated by using what is known as a **preconditio
44
an approximation to the matrix inverse action which is cheap to evaluate. These
55
can improve the numerical conditioning of the solver process and in turn improve
66
the performance. LinearSolve.jl provides an interface for the definition of
7-
preconditioners which works with the wrapped packages.
7+
preconditioners which works with the wrapped iterative solver packages.
88

99
## Using Preconditioners
1010

1111
### Mathematical Definition
1212

13-
Preconditioners are specified in the keyword arguments to `init` or `solve`: `Pl` for left
14-
and `Pr` for right preconditioner, respectively.
15-
The right preconditioner, ``P_r`` transforms the linear system ``Au = b`` into the form:
13+
A right preconditioner, ``P_r`` transforms the linear system ``Au = b`` into the form:
1614

1715
```math
1816
AP_r^{-1}(P_r u) = AP_r^{-1}y = b
@@ -31,30 +29,73 @@ A two-sided preconditioned system is of the form:
3129
P_l^{-1}A P_r^{-1} (P_r u) = P_l^{-1}b
3230
```
3331

34-
By default, if no preconditioner is given, the preconditioner is assumed to be
35-
the identity ``I``.
32+
### Specifying Preconditioners
3633

37-
### Using Preconditioners
34+
One way to specify preconditioners uses the `Pl` and `Pr` keyword arguments to `init` or `solve`: `Pl` for left
35+
and `Pr` for right preconditioner, respectively. By default, if no preconditioner is given, the preconditioner is assumed to be
36+
the identity ``I``.
3837

39-
In the following, we will use the `DiagonalPreconditioner` to define a two-sided
40-
preconditioned system which first divides by some random numbers and then
41-
multiplies by the same values. This is commonly used in the case where if, instead
42-
of random, `s` is an approximation to the eigenvalues of a system.
38+
In the following, we will use a left sided diagonal (Jacobi) preconditioner.
4339

44-
```@example precon
40+
```@example precon1
4541
using LinearSolve, LinearAlgebra
4642
n = 4
47-
s = rand(n)
48-
Pl = Diagonal(s)
4943
5044
A = rand(n, n)
5145
b = rand(n)
5246
47+
Pl = Diagonal(A)
48+
5349
prob = LinearProblem(A, b)
5450
sol = solve(prob, KrylovJL_GMRES(), Pl = Pl)
5551
sol.u
5652
```
5753

54+
Alternatively, preconditioners can be specified via the `precs` argument to the constructor of
55+
an iterative solver specification. This argument shall deliver a factory method mapping `A` and a
56+
parameter `p` to a tuple `(Pl,Pr)` consisting a left and a right preconditioner.
57+
58+
```@example precon2
59+
using LinearSolve, LinearAlgebra
60+
n = 4
61+
62+
A = rand(n, n)
63+
b = rand(n)
64+
65+
prob = LinearProblem(A, b)
66+
sol = solve(prob, KrylovJL_GMRES(precs = (A, p) -> (Diagonal(A), I)))
67+
sol.u
68+
```
69+
70+
This approach has the advantage that the specification of the preconditioner is possible without
71+
the knowledge of a concrete matrix `A`. It also allows to specify the preconditioner via a callable object
72+
and to pass parameters to the constructor of the preconditioner instances. The example below also shows how
73+
to reuse the preconditioner once constructed for the subsequent solution of a modified problem.
74+
75+
```@example precon3
76+
using LinearSolve, LinearAlgebra
77+
78+
Base.@kwdef struct WeightedDiagonalPreconBuilder
79+
w::Float64
80+
end
81+
82+
(builder::WeightedDiagonalPreconBuilder)(A, p) = (builder.w * Diagonal(A), I)
83+
84+
n = 4
85+
A = n * I - rand(n, n)
86+
b = rand(n)
87+
88+
prob = LinearProblem(A, b)
89+
sol = solve(prob, KrylovJL_GMRES(precs = WeightedDiagonalPreconBuilder(w = 0.9)))
90+
sol.u
91+
92+
B = A .+ 0.1
93+
cache = sol.cache
94+
reinit!(cache, A = B, reuse_precs = true)
95+
sol = solve!(cache, KrylovJL_GMRES(precs = WeightedDiagonalPreconBuilder(w = 0.9)))
96+
sol.u
97+
```
98+
5899
## Preconditioner Interface
59100

60101
To define a new preconditioner you define a Julia type which satisfies the

test/resolve.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@ using LinearSolve, LinearAlgebra, SparseArrays, InteractiveUtils, Test
22
using LinearSolve: AbstractDenseFactorization, AbstractSparseFactorization
33

44
for alg in vcat(InteractiveUtils.subtypes(AbstractDenseFactorization),
5-
InteractiveUtils.subtypes(AbstractSparseFactorization))
5+
InteractiveUtils.subtypes(AbstractSparseFactorization))
66
if alg in [PardisoJL]
77
## Pardiso has extra tests in test/pardiso/pardiso.jl
88
continue

0 commit comments

Comments
 (0)