Introduction
Installing
Install KrylovKit.jl via the package manager:
using Pkg
Pkg.add("KrylovKit")KrylovKit.jl is a pure Julia package; no dependencies (aside from the Julia standard library) are required.
Getting started
After installation, start by loading KrylovKit
using KrylovKitThe help entry of the KrylovKit module states
KrylovKit — ModuleKrylovKitA Julia package collecting a number of Krylov-based algorithms for linear problems, singular value and eigenvalue problems and the application of functions of linear maps or operators to vectors.
KrylovKit accepts general functions or callable objects as linear maps, and general Julia objects with vector like behavior as vectors.
The high level interface of KrylovKit is provided by the following functions:
linsolve: solve linear systemseigsolve: find a few eigenvalues and corresponding eigenvectorsgeneigsolve: find a few generalized eigenvalues and corresponding vectorssvdsolve: find a few singular values and corresponding left and right singular vectorsexponentiate: apply the exponential of a linear map to a vectorexpintegrator: exponential integrator for a linear non-homogeneous ODE, computes a linear combination of theϕⱼfunctions which generalizeϕ₀(z) = exp(z).
Common interface
The for high-level function linsolve, eigsolve, geneigsolve, svdsolve, exponentiate and expintegrator follow a common interface
results..., info = problemsolver(A, args...; kwargs...)where problemsolver is one of the functions above. Here, A is the linear map in the problem, which could be an instance of AbstractMatrix, or any function or callable object that encodes the action of the linear map on a vector. In particular, one can write the linear map using Julia's do block syntax as
results..., info = problemsolver(args...; kwargs...) do x
y = # implement linear map on x
return y
endRead the documentation for problems that require both the linear map and its adjoint to be implemented, e.g. svdsolve, or that require two different linear maps, e.g. geneigsolve.
Furthermore, args is a set of additional arguments to specify the problem. The keyword arguments kwargs contain information about the linear map (issymmetric, ishermitian, isposdef) and about the solution strategy (tol, krylovdim, maxiter). Finally, there is a keyword argument verbosity that determines how much information is printed to STDOUT. The default value verbosity = SILENT_LEVEL means that no information will be printed. With verbosity = WARN_LEVEL, a single message at the end of the algorithm will be displayed, which is a warning if the algorithm did not succeed in finding the solution, or some information if it did. For verbosity = STARTSTOP_LEVEL, information about the current state is displayed after every iteration of the algorithm. Finally, for verbosity > STARTSTOP_LEVEL, information about the individual Krylov expansion steps is displayed.
The return value contains one or more entries that define the solution, and a final entry info of type ConvergeInfo that encodes information about the solution, i.e. whether it has converged, the residual(s) and the norm thereof, the number of operations used:
KrylovKit.ConvergenceInfo — Typestruct ConvergenceInfo{S,T}
converged::Int
residual::T
normres::S
numiter::Int
numops::Int
endUsed to return information about the solution found by the iterative method.
converged: the number of solutions that have converged according to an appropriate error measure and requested tolerance for the problem. Its value can be zero or one forlinsolve,exponentiateandexpintegrator, or any integer>= 0foreigsolve,schursolveorsvdsolve.residual:the (list of) residual(s) for the problem, ornothingfor problems without the concept of a residual (i.e.exponentiate,expintegrator). This is a single vector (of the same type as the type of vectors used in the problem) forlinsolve, or aVectorof such vectors foreigsolve,schursolveorsvdsolve.normres: the norm of the residual(s) (in the previous field) or the value of any other error measure that is appropriate for the problem. This is aRealforlinsolveandexponentiate, and aVector{<:Real}foreigsolve,schursolveandsvdsolve. The number of values innormresthat are smaller than a predefined tolerance corresponds to the numberconvergedof solutions that have converged.numiter: the number of iterations (sometimes called restarts) used by the algorithm.numops: the number of times the linear map or operator was applied
There is also an expert interface where the user specifies the algorithm that should be used explicitly, i.e.
results..., info = problemsolver(A, args..., algorithm(; kwargs...))Most algorithm constructions take the same keyword arguments (tol, krylovdim, maxiter and verbosity) discussed above.
While KrylovKit.jl does currently not provide a general interface for including preconditioners, it is possible to e.g. use a modified inner product. KrylovKit.jl provides a specific type for this purpose:
KrylovKit.InnerProductVec — Typev = InnerProductVec(vec, dotf)Create a new vector v from an existing vector dotf with a modified inner product given by inner. The vector vec, which can be any type (not necessarily Vector) that supports the basic vector interface required by KrylovKit, is wrapped in a custom struct v::InnerProductVec. All vector space functionality such as addition and multiplication with scalars (both out of place and in place using mul!, rmul!, axpy! and axpby!) applied to v is simply forwarded to v.vec. The inner product between vectors v1 = InnerProductVec(vec1, dotf) and v2 = InnerProductVec(vec2, dotf) is computed as dot(v1, v2) = dotf(v1.vec, v2.vec) = dotf(vec1, vec2). The inner product between vectors with different dotf functions is not defined. Similarly, The norm of v::InnerProductVec is defined as v = sqrt(real(dot(v, v))) = sqrt(real(dotf(vec, vec))).
In a (linear) map applied to v, the original vector can be obtained as v.vec or simply as v[].