/
linsolve.jl
178 lines (157 loc) · 7.48 KB
/
linsolve.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
"""
linsolve(A::AbstractMatrix, b::AbstractVector, [x₀, a₀::Number = 0, a₁::Number = 1];
kwargs...)
linsolve(f, b, [x₀, a₀::Number = 0, a₁::Number = 1]; kwargs...)
# expert version:
linsolve(f, b, x₀, algorithm, [a₀::Number = 0, a₁::Number = 1])
Compute a solution `x` to the linear system `(a₀ + a₁ * A)*x = b` or
`a₀ * x + a₁ * f(x) = b`, possibly using a starting guess `x₀`. Return the approximate
solution `x` and a `ConvergenceInfo` structure.
### Arguments:
The linear map can be an `AbstractMatrix` (dense or sparse) or a general function or
callable object. If no initial guess is specified, it is chosen as `(zero(a₀)*zero(a₁))*b`
which should generate an object similar to `b` but initialized with zeros. The numbers `a₀`
and `a₁` are optional arguments; they are applied implicitly, i.e. they do not contribute
the computation time of applying the linear map or to the number of operations on vectors of
type `x` and `b`.
### Return values:
The return value is always of the form `x, info = linsolve(...)` with
- `x`: the approximate solution to the problem, similar type as the right hand side `b`
but possibly with a different `eltype`
- `info`: an object of type [`ConvergenceInfo`], which has the following fields
+ `info.converged::Int`: takes value 0 or 1 depending on whether the solution was
converged up to the requested tolerance
+ `info.residual`: residual `b - f(x)` of the approximate solution `x`
+ `info.normres::Real`: norm of the residual, i.e. `norm(info.residual)`
+ `info.numops::Int`: total number of times that the linear map was applied, i.e. the
number of times that `f` was called, or a vector was multiplied with `A`
+ `info.numiter::Int`: number of times the Krylov subspace was restarted (see below)
!!! warning "Check for convergence"
No warning is printed if no converged solution was found, so always check if
`info.converged == 1`.
### Keyword arguments:
Keyword arguments are given by:
- `verbosity::Int = 0`: verbosity level, i.e. 0 (no messages), 1 (single message
at the end), 2 (information after every iteration), 3 (information per Krylov step)
- `atol::Real`: the requested accuracy, i.e. absolute tolerance, on the norm of the
residual.
- `rtol::Real`: the requested accuracy on the norm of the residual, relative to the norm
of the right hand side `b`.
- `tol::Real`: the requested accuracy on the norm of the residual that is actually used by
the algorithm; it defaults to `max(atol, rtol*norm(b))`. So either use `atol` and `rtol`
or directly use `tol` (in which case the value of `atol` and `rtol` will be ignored).
- `krylovdim::Integer`: the maximum dimension of the Krylov subspace that will be
constructed.
- `maxiter::Integer`: the number of times the Krylov subspace can be rebuilt; see below for
further details on the algorithms.
- `orth::Orthogonalizer`: the orthogonalization method to be used, see
[`Orthogonalizer`](@ref)
- `issymmetric::Bool`: if the linear map is symmetric, only meaningful if `T<:Real`
- `ishermitian::Bool`: if the linear map is hermitian
- `isposdef::Bool`: if the linear map is positive definite
The default values are given by `atol = KrylovDefaults.tol`, `rtol = KrylovDefaults.tol`,
`tol = max(atol, rtol*norm(b))`, `krylovdim = KrylovDefaults.krylovdim`,
`maxiter = KrylovDefaults.maxiter`, `orth = KrylovDefaults.orth`;
see [`KrylovDefaults`](@ref) for details.
The default value for the last three parameters depends on the method. If an
`AbstractMatrix` is used, `issymmetric`, `ishermitian` and `isposdef` are checked for that
matrix, ortherwise the default values are `issymmetric = false`,
`ishermitian = T <: Real && issymmetric` and `isposdef = false`.
### Algorithms
The final (expert) method, without default values and keyword arguments, is the one that is
finally called, and can also be used directly. Here, one specifies the algorithm explicitly.
Currently, only [`CG`](@ref) and [`GMRES`](@ref) are implemented, where `CG` is chosen if
`isposdef == true`. Note that in standard `GMRES` terminology, our parameter `krylovdim` is
referred to as the *restart* parameter, and our `maxiter` parameter counts the number of
outer iterations, i.e. restart cycles. In `CG`, the Krylov subspace is only implicit
because short recurrence relations are being used, and therefore no restarts are required.
Therefore, we pass `krylovdim*maxiter` as the maximal number of CG iterations that can be
used by the `CG` algorithm.
"""
function linsolve end
linsolve(A::AbstractMatrix, b::AbstractVector, a₀::Number = 0, a₁::Number = 1; kwargs...) =
linsolve(A, b, (zero(a₀) * zero(a₁)) * b, a₀, a₁; kwargs...)
linsolve(f, b, a₀::Number = 0, a₁::Number = 1; kwargs...) =
linsolve(f, b, (zero(a₀) * zero(a₁)) * b, a₀, a₁; kwargs...)
function linsolve(f, b, x₀, a₀::Number = 0, a₁::Number = 1; kwargs...)
Tx = promote_type(typeof(x₀))
Tb = typeof(b)
Tfx = Core.Compiler.return_type(apply, Tuple{typeof(f),Tx})
T = promote_type(Core.Compiler.return_type(dot, Tuple{Tb,Tfx}), typeof(a₀), typeof(a₁))
alg = linselector(f, b, T; kwargs...)
return linsolve(f, b, x₀, alg, a₀, a₁)
end
function linselector(
f,
b,
T::Type;
issymmetric::Bool = false,
ishermitian::Bool = T <: Real && issymmetric,
isposdef::Bool = false,
krylovdim::Int = KrylovDefaults.krylovdim,
maxiter::Int = KrylovDefaults.maxiter,
rtol::Real = KrylovDefaults.tol,
atol::Real = KrylovDefaults.tol,
tol::Real = max(atol, rtol * norm(b)),
orth = KrylovDefaults.orth,
verbosity::Int = 0
)
if (T <: Real && issymmetric) || ishermitian
isposdef &&
return CG(maxiter = krylovdim * maxiter, tol = tol, verbosity = verbosity)
# TODO: implement MINRES for symmetric but not posdef; for now use GRMES
# return MINRES(krylovdim*maxiter, tol=tol)
return GMRES(
krylovdim = krylovdim,
maxiter = maxiter,
tol = tol,
orth = orth,
verbosity = verbosity
)
else
return GMRES(
krylovdim = krylovdim,
maxiter = maxiter,
tol = tol,
orth = orth,
verbosity = verbosity
)
end
end
function linselector(
A::AbstractMatrix,
b,
T::Type;
issymmetric::Bool = T <: Real && LinearAlgebra.issymmetric(A),
ishermitian::Bool = issymmetric || LinearAlgebra.ishermitian(A),
isposdef::Bool = ishermitian ? LinearAlgebra.isposdef(A) : false,
krylovdim::Int = KrylovDefaults.krylovdim,
maxiter::Int = KrylovDefaults.maxiter,
rtol::Real = KrylovDefaults.tol,
atol::Real = KrylovDefaults.tol,
tol::Real = max(atol, rtol * norm(b)),
orth = KrylovDefaults.orth,
verbosity::Int = 0
)
if (T <: Real && issymmetric) || ishermitian
isposdef &&
return CG(maxiter = krylovdim * maxiter, tol = tol, verbosity = verbosity)
# TODO: implement MINRES for symmetric but not posdef; for now use GRMES
# return MINRES(krylovdim*maxiter, tol=tol, verbosity = verbosity)
return GMRES(
krylovdim = krylovdim,
maxiter = maxiter,
tol = tol,
orth = orth,
verbosity = verbosity
)
else
return GMRES(
krylovdim = krylovdim,
maxiter = maxiter,
tol = tol,
orth = orth,
verbosity = verbosity
)
end
end