Solving linear systems
Linear systems are of the form A*x=b where A should be a linear map that has the same type of output as input, i.e. the solution x should be of the same type as the right hand side b. They can be solved using the function linsolve:
KrylovKit.linsolve — Function.linsolve(A::AbstractMatrix, b::AbstractVector, T::Type = promote_type(eltype(A), eltype(b)); kwargs...)
linsolve(f, b, T::Type = eltype(b); kwargs...)
linsolve(f, b, x₀; kwargs...)
linsolve(f, b, x₀, algorithm)Compute a solution x to the linear system A*x = b or 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 fill!(similar(b, T), 0). The argument T acts as a hint in which Number type the computation should be performed, but is not restrictive. If the linear map automatically produces complex values, complex arithmetic will be used even though T<:Real was specified.
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 sidebbut possibly with a differenteltypeinfo: an object of type [ConvergenceInfo], which has the following fieldsinfo.converged::Int: takes value 0 or 1 depending on whether the solution was converged up to the requested toleranceinfo.residual: residualb - f(x)of the approximate solutionxinfo.normres::Real: norm of the residual, i.e.norm(info.residual)info.numops::Int: number of times the linear map was applied, i.e. number of timesfwas called, or a vector was multiplied withAinfo.numiter::Int: number of times the Krylov subspace was restarted (see below)
No warning is printed if not all requested eigenvalues were converged, so always check if info.converged == 1.
Keyword arguments:
Keyword arguments and their default values are given by:
krylovdim = 30: the maximum dimension of the Krylov subspace that will be constructed. Note that the dimension of the vector space is not known or checked, e.g.x₀should not necessarily support theBase.lengthfunction. If you know the actual problem dimension is smaller than the default value, it is useful to reduce the value ofkrylovdim, though in principle this should be detected.tol = 1e-12: the requested accuracy (corresponding to the 2-norm of the residual for Schur vectors, not the eigenvectors). If you work in e.g. single precision (Float32), you should definitely change the default value.maxiter = 100: the number of times the Krylov subspace can be rebuilt; see below for further details on the algorithms.issymmetric: if the linear map is symmetric, only meaningful ifT<:Realishermitian: if the linear map is hermitian
The default value for the last two depends on the method. If an AbstractMatrix is used, issymmetric and ishermitian are checked for that matrix, ortherwise the default values are issymmetric = false and ishermitian = T <: Real && issymmetric.
The keyword argument tol represents the actual tolerance, i.e. the solution x is deemed converged if norm(b - f(x)) < tol. If you want to use some relative tolerance rtol, use tol = rtol*norm(b) as keyword argument.
Algorithm
The last 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 GMRES is implemented. Note that we do not use the standard GMRES terminology, where every new application of the linear map (matrix vector multiplications) counts as a new iteration, and there is a restart parameter that indicates the largest Krylov subspace that is being built. Indeed, the keyword argument krylovdim corresponds to the traditional restart parameter, whereas maxiter counts the number of outer iterations, i.e. the number of times this Krylov subspace is built. So the maximal number of operations of the linear map is roughly maxiter * krylovdim.