Skip to content

Commit 6a174ed

Browse files
ChrisRackauckasChrisRackauckas
committed
Start formalizing the preconditioner interface
1 parent d52408b commit 6a174ed

File tree

4 files changed

+147
-62
lines changed

4 files changed

+147
-62
lines changed

‎docs/src/basics/Preconditioners.md

Lines changed: 86 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,86 @@
1+
# Preconditioners
2+
3+
Many linear solvers can be accelerated by using what is known as a **preconditioner**,
4+
an approximation to the matrix inverse action which is cheap to evaluate. These
5+
can improve the numerical conditioning of the solver process and in turn improve
6+
the performance. LinearSolve.jl provides an interface for the definition of
7+
preconditioners which works with the wrapped packages.
8+
9+
## Using Preconditioners
10+
11+
### Mathematical Definition
12+
13+
Preconditioners are specified in the keyword arguments of `init` or `solve`. The
14+
right preconditioner, `Pr` transforms the linear system ``Au = b`` into the
15+
form:
16+
17+
```math
18+
AP_r^{-1}(Pu) = AP_r^{-1}y = b
19+
```
20+
21+
to add the solving step ``P_r u = y``. The left preconditioner, `Pl`, transforms
22+
the linear system into the form:
23+
24+
```math
25+
P_l^{-1}(Au - b) = 0
26+
```
27+
28+
A two-sided preconditioned system is of the form:
29+
30+
```math
31+
P_l A P_r^{-1} (P_r u) = P_l b
32+
```
33+
34+
By default, if no preconditioner is given the preconditioner is assumed to be
35+
the identity ``I``.
36+
37+
### Using Preconditioners
38+
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.
43+
44+
```julia
45+
s = rand(n)
46+
47+
# Pr applies 1 ./ s .* vec
48+
Pr = LinearSolve.DiagonalPreconditioner(s)
49+
# Pl applies s .* vec
50+
Pl = LinearSolve.DiagonalPreconditioner(s)
51+
52+
A = rand(n,n)
53+
b = rand(n)
54+
55+
prob = LinearProblem(A,b)
56+
sol = solve(prob,IterativeSolvers_GMRES(),Pl=Pl,Pr=Pr)
57+
```
58+
59+
## Pre-Defined Preconditioners
60+
61+
To simplify the usage of preconditioners, LinearSolve.jl comes with many standard
62+
preconditioners written to match the required interface.
63+
64+
- `DiagonalPreconditioner(s::Union{Number,AbstractVector})`: the diagonal
65+
preconditioner, defined as a diagonal matrix `Diagonal(s)`.
66+
67+
## Preconditioner Interface
68+
69+
To define a new preconditioner you define a Julia type which satisfies the
70+
following interface:
71+
72+
### General
73+
74+
- `Base.eltype(::Preconditioner)`
75+
- `Base.adjoint(::Preconditioner)`
76+
- `Base.inv(::Preconditioner)` (Optional?)
77+
78+
### Required for Right Preconditioners
79+
80+
- `Base.\(::Preconditioner,::AbstractVector)`
81+
- `LinearAlgebra.ldiv!(::AbstractVector,::Preconditioner,::AbstractVector)`
82+
83+
### Required for Left Preconditioners
84+
85+
- `Base.*(::Preconditioner,::AbstractVector)`
86+
- `LinearAlgebra.mul!(::AbstractVector,::Preconditioner,::AbstractVector)`

‎src/LinearSolve.jl

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,8 @@ abstract type AbstractKrylovSubspaceMethod <: SciMLLinearSolveAlgorithm end
2929

3030
include("common.jl")
3131
include("factorization.jl")
32-
include("wrappers.jl")
32+
include("iterative_wrappers.jl")
33+
include("preconditioners.jl")
3334
include("default.jl")
3435

3536
const IS_OPENBLAS = Ref(true)

‎src/wrappers.jl renamed to ‎src/iterative_wrappers.jl

Lines changed: 0 additions & 61 deletions
Original file line numberDiff line numberDiff line change
@@ -1,64 +1,3 @@
1-
2-
## Preconditioners
3-
4-
scaling_preconditioner(s::Number) = I * (1/s), I * s
5-
scaling_preconditioner(s::AbstractVector) = Diagonal(inv.(s)),Diagonal(s)
6-
7-
struct ComposePreconditioner{Ti,To}
8-
inner::Ti
9-
outer::To
10-
end
11-
12-
Base.eltype(A::ComposePreconditioner) = promote_type(eltype(A.inner), eltype(A.outer))
13-
Base.adjoint(A::ComposePreconditioner) = ComposePreconditioner(A.outer', A.inner')
14-
Base.inv(A::ComposePreconditioner) = InvComposePreconditioner(A)
15-
16-
function LinearAlgebra.ldiv!(A::ComposePreconditioner, x)
17-
@unpack inner, outer = A
18-
19-
ldiv!(inner, x)
20-
ldiv!(outer, x)
21-
end
22-
23-
function LinearAlgebra.ldiv!(y, A::ComposePreconditioner, x)
24-
@unpack inner, outer = A
25-
26-
ldiv!(y, inner, x)
27-
ldiv!(outer, y)
28-
end
29-
30-
struct InvComposePreconditioner{Tp <: ComposePreconditioner}
31-
P::Tp
32-
end
33-
34-
InvComposePreconditioner(inner, outer) = InvComposePreconditioner(ComposePreconditioner(inner, outer))
35-
36-
Base.eltype(A::InvComposePreconditioner) = Base.eltype(A.P)
37-
Base.adjoint(A::InvComposePreconditioner) = InvComposePreconditioner(A.P')
38-
Base.inv(A::InvComposePreconditioner) = deepcopy(A.P)
39-
40-
function LinearAlgebra.mul!(y, A::InvComposePreconditioner, x)
41-
@unpack P = A
42-
ldiv!(y, P, x)
43-
end
44-
45-
function get_preconditioner(Pi, Po)
46-
47-
ifPi = Pi !== Identity()
48-
ifPo = Po !== Identity()
49-
50-
P =
51-
if ifPi & ifPo
52-
ComposePreconditioner(Pi, Po)
53-
elseif ifPi | ifPo
54-
ifPi ? Pi : Po
55-
else
56-
Identity()
57-
end
58-
59-
return P
60-
end
61-
621
## Krylov.jl
632

643
struct KrylovJL{F,Tl,Tr,I,A,K} <: AbstractKrylovSubspaceMethod

‎src/preconditioners.jl

Lines changed: 59 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,59 @@
1+
## Preconditioners
2+
3+
scaling_preconditioner(s::Number) = I * (1/s), I * s
4+
scaling_preconditioner(s::AbstractVector) = Diagonal(inv.(s)),Diagonal(s)
5+
6+
struct ComposePreconditioner{Ti,To}
7+
inner::Ti
8+
outer::To
9+
end
10+
11+
Base.eltype(A::ComposePreconditioner) = promote_type(eltype(A.inner), eltype(A.outer))
12+
Base.adjoint(A::ComposePreconditioner) = ComposePreconditioner(A.outer', A.inner')
13+
Base.inv(A::ComposePreconditioner) = InvComposePreconditioner(A)
14+
15+
function LinearAlgebra.ldiv!(A::ComposePreconditioner, x)
16+
@unpack inner, outer = A
17+
18+
ldiv!(inner, x)
19+
ldiv!(outer, x)
20+
end
21+
22+
function LinearAlgebra.ldiv!(y, A::ComposePreconditioner, x)
23+
@unpack inner, outer = A
24+
25+
ldiv!(y, inner, x)
26+
ldiv!(outer, y)
27+
end
28+
29+
struct InvComposePreconditioner{Tp <: ComposePreconditioner}
30+
P::Tp
31+
end
32+
33+
InvComposePreconditioner(inner, outer) = InvComposePreconditioner(ComposePreconditioner(inner, outer))
34+
35+
Base.eltype(A::InvComposePreconditioner) = Base.eltype(A.P)
36+
Base.adjoint(A::InvComposePreconditioner) = InvComposePreconditioner(A.P')
37+
Base.inv(A::InvComposePreconditioner) = deepcopy(A.P)
38+
39+
function LinearAlgebra.mul!(y, A::InvComposePreconditioner, x)
40+
@unpack P = A
41+
ldiv!(y, P, x)
42+
end
43+
44+
function get_preconditioner(Pi, Po)
45+
46+
ifPi = Pi !== Identity()
47+
ifPo = Po !== Identity()
48+
49+
P =
50+
if ifPi & ifPo
51+
ComposePreconditioner(Pi, Po)
52+
elseif ifPi | ifPo
53+
ifPi ? Pi : Po
54+
else
55+
Identity()
56+
end
57+
58+
return P
59+
end

0 commit comments

Comments
 (0)