@@ -18,19 +18,21 @@ julia> Pkg.checkout("IterativeSolvers")
18
18
## Interface
19
19
20
20
All linear-algebraic routines will take as input a linear operator ` A ` that maps
21
- vectors to vectors. ` A ` is not explicitly typed, but must either be a
22
- [ ` KrylovSubspace ` ] ( @ref ) or support multiplication ` * ` or function composition (apply)
23
- that behave as necessary to produce the correct mapping on the vector space.
21
+ vectors to vectors. Typically ` A ` is a ` Matrix ` or a ` SparseMatrixCSC ` , but since
22
+ ` A ` is not explicitly typed, any linear operator that supports matrix operations
23
+ can be used as well. This makes it possible to apply solvers * matrix-free* . In
24
+ IterativeSolvers.jl we strongly recommend [ LinearMaps.jl] ( https://github.com/Jutho/LinearMaps.jl )
25
+ for non-matrix types of ` A ` .
24
26
25
- A custom type for ` A ` may be specified. The following interface is expected to
26
- be defined on ` A ` :
27
+ For matrix-free types of ` A ` the following interface is expected to be defined:
27
28
28
- ` A*v ` is defined and computes the matrix-vector product on a ` v::Vector ` .
29
+ ` A*v ` computes the matrix-vector product on a ` v::AbstractVector ` .
29
30
30
- ` eltype(A) ` is defined and returns the element type implicit in the equivalent
31
- matrix representation of ` A ` .
31
+ ` A_mul_B!(y, A, v) ` computes the matrix-vector product on a ` v::AbstractVector ` in-place.
32
32
33
- ` size(A, d) ` is defined and returns the nominal dimensions along the dth axis
33
+ ` eltype(A) ` returns the element type implicit in the equivalent matrix representation of ` A ` .
34
+
35
+ ` size(A, d) ` is defined and returns the nominal dimensions along the ` d ` th axis
34
36
in the equivalent matrix representation of ` A ` .
35
37
36
38
### Solvers
@@ -158,109 +160,3 @@ plot(ch, :resnorm, sep = :blue)
158
160
* Plot additional keywords*
159
161
160
162
` sep::Symbol = :white ` : color of the line separator in restarted methods.
161
-
162
- ## KrylovSubspace
163
-
164
- When [ ` KrylovSubspace ` ] ( @ref ) is supported by the method, ` A ` can be replaced
165
- by an instance ` K ` of it. To check if a certain function supports it use
166
- ` methods ` or ` ? ` to find out, if there is a ` K ` in the argument list then it does.
167
-
168
- ``` julia
169
- julia> ?cg!
170
-
171
- cg! (x, A, b)
172
- cg! (x, K, b)
173
-
174
- ...
175
- ```
176
-
177
- ` KrylovSubspace ` is only allowed on functions that accept ` x ` as an argument,
178
- that's why ` cg ` doesn't allow it. This is also true when ` x ` is a keyword as
179
- is the case of ` powm ` .
180
-
181
- ``` julia
182
- julia> ?cg!
183
-
184
- powm (A)
185
- powm (K)
186
-
187
- ...
188
- ```
189
-
190
- The ` KrylovSubspace ` type collects information on the Krylov subspace generated
191
- over the course of an iterative Krylov solver.
192
-
193
- Recall that the Krylov subspace of order r given a starting vector ` b ` and a
194
- linear operator ` A ` is spanned by the vectors ` [b, A*b, A^2*b,... A^(r-1)*b] ` .
195
- Many modern iterative solvers work on Krylov spaces which expand until they span
196
- enough of the range of ` A ` for the solution vector to converge. Most practical
197
- algorithms, however, will truncate the order of the Krylov subspace to save
198
- memory and control the accumulation of roundoff errors. Additionally, they do
199
- not work directly with the raw iterates ` A^n*b ` , but will orthogonalize
200
- subsequent vectors as they are computed so as to improve numerical stability.
201
- ` KrylovSubspace ` s provide a natural framework to express operations such as a
202
- (possibly non-orthonormalized) basis for the Krylov subspace, retrieving the
203
- next vector in the subspace, and orthogonalizing an arbitrary vector against
204
- (part or all of) the subspace.
205
-
206
- The implementation of ` KrylovSubspace ` in this package differs from standard
207
- textbook presentations of iterative solvers. First, the ` KrylovSubspace ` type
208
- shows clearly the relationship between the linear operator A and the sequence of
209
- basis vectors for the Krylov subspace that is generated by each new iteration.
210
- Second, the grouping together of basis vectors also makes the orthogonalization
211
- steps in each iteration routine clearer. Third, unlike in other languages, the
212
- powerful type system of Julia allows for a simplified notation for iterative
213
- algorithms without compromising on performance, and enhances code reuse across
214
- multiple solvers.
215
-
216
- ### Constructors
217
-
218
- A [ ` KrylovSubspace ` ] ( @ref ) can be initialized by three constructors, depending on the type
219
- of the linear operator:
220
-
221
- ```
222
- KrylovSubspace{T}(A::AbstractMatrix{T}, [order::Int, v::Vector{Vector{T}}])
223
-
224
- KrylovSubspace{T}(A::KrylovSubspace{T}, [order::Int, v::Vector{Vector{T}}])
225
-
226
- KrylovSubspace(A, n::Int, [order::Int, v::Vector{Vector{T}}])
227
- ```
228
-
229
- ` A ` : The linear operator associated with the ` KrylovSubspace ` .
230
-
231
- ` order ` : the order of the ` KrylovSubspace ` , i.e. the maximal number of Krylov
232
- vectors to remember.
233
-
234
- ` n ` : the dimensionality of the underlying vector space ` T^n ` .
235
-
236
- ` v ` : the iterable collection of Krylov vectors (of maximal length order).
237
-
238
- The dimensionality of the underlying vector space is automatically inferred where
239
- possible, i.e. when the linear operator is an ` AbstractMatrix ` or ` KrylovSubspace ` .
240
-
241
- * ` Note: ` * second constructor destroys the old ` KrylovSubspace ` .
242
-
243
- ### Orthogonalization
244
-
245
- Orthogonalizing the basis vectors for a ` KrylovSubspace ` is crucial for numerical
246
- stability, and is a core low-level operation for many iterative solvers.
247
-
248
- ```
249
- orthogonalize{T}(v::Vector{T}, K::KrylovSubspace{T}, [p::Int]; [method::Symbol], [normalize::Bool])
250
- ```
251
-
252
- ` v ` : vector to orthogonalize.
253
-
254
- ` K ` : ` KrylovSubspace ` to orthogonalize against.
255
-
256
- ` p ` : number of Krylov vectors to orthogonalize against. (Default is all available)
257
-
258
- ` method ` : Orthogonalization method. Currently supported methods are:
259
-
260
- * ` :GramSchmidt `
261
- * ` :ModifiedGramSchmidt ` (default)
262
- * ` :Householder ` .
263
-
264
- ## Define function as matrix
265
-
266
- This package supports the use of [ LinearMaps.jl] ( https://github.com/Jutho/LinearMaps.jl ) for creating matrices using custom functions.
0 commit comments