GitHub

In-place methods

All solvers in Krylov.jl have an in-place variant implemented in a method whose name ends with !. A workspace (KrylovSolver), which contains the storage needed by a Krylov method, can be used to solve multiple linear systems with the same dimensions and the same floating-point precision. The section storage requirements specifies the memory needed for each Krylov method. Each KrylovSolver has three constructors:

XyzSolver(A, b)
XyzSolver(m, n, S)
XyzSolver(kc::KrylovConstructor)

Xyz represents the name of the Krylov method, written in lowercase except for its first letter (e.g., Cg, Minres, Lsmr, Bicgstab, etc.). If the name of the Krylov method contains an underscore (e.g., minres_qlp or cgls_lanczos_shift), the workspace constructor transforms it by capitalizing each word and removing underscores, resulting in names like MinresQlpSolver or CglsLanczosShiftSolver.

Given an operator A and a right-hand side b, you can create a KrylovSolver based on the size of A and the type of b, or explicitly provide the dimensions (m, n) and the storage type S. We assume that S(undef, 0), S(undef, n), and S(undef, m) are well-defined for the storage type S. For more advanced vector types, workspaces can also be created with the help of a KrylovConstructor.

Krylov.KrylovConstructor — Type
KrylovConstructor(vm; vm_empty=vm)
KrylovConstructor(vm, vn; vm_empty=vm, vn_empty=vn)

Krylov methods require a workspace containing vectors of length m and n to solve linear problems of size m × n. The KrylovConstructor facilitates the allocation of these vectors using similar.

For square problems (m == n), use the first constructor with a single vector vm. For rectangular problems (m ≠ n), use the second constructor with vm and vn.

Input arguments

  • vm: a vector of length m;
  • vn: a vector of length n.

Keyword arguments

  • vm_empty: an empty vector that may be replaced with a vector of length m;
  • vn_empty: an empty vector that may be replaced with a vector of length n.

Note

Empty vectors vm_empty and vn_empty reduce storage requirements when features such as warm-start or preconditioners are unused. These empty vectors will be replaced within a KrylovSolver only if required, such as when preconditioners are provided.

source

See the section custom workspaces for an example where this constructor is the only applicable option.

For example, use S = Vector{Float64} if you want to solve linear systems in double precision on the CPU and S = CuVector{Float32} if you want to solve linear systems in single precision on an Nvidia GPU.

Note

DiomSolver, FomSolver, DqgmresSolver, GmresSolver, BlockGmresSolver, FgmresSolver, GpmrSolver, CgLanczosShiftSolver and CglsLanczosShiftSolver require an additional argument (memory or nshifts).

The workspace is always the first argument of the in-place methods:

minres_solver = MinresSolver(n, n, Vector{Float64})
minres!(minres_solver, A1, b1)

dqgmres_solver = DqgmresSolver(n, n, memory, Vector{BigFloat})
dqgmres!(dqgmres_solver, A2, b2)

lsqr_solver = LsqrSolver(m, n, CuVector{Float32})
lsqr!(lsqr_solver, A3, b3)

A generic function solve! is also available and dispatches to the appropriate Krylov method.

Krylov.solve! — Function
solve!(solver, args...; kwargs...)

Generic function that dispatches to the appropriate in-place Krylov method based on the type of solver.

source

In-place methods return an updated solver workspace. Solutions and statistics can be recovered via solver.x, solver.y and solver.stats. Functions solution, statistics and results can be also used.

Krylov.issolved — Function
issolved(solver)

Return a boolean that determines whether the Krylov method associated to solver succeeded.

source
Krylov.solution — Function
solution(solver)

Return the solution(s) stored in the solver. Optionally you can specify which solution you want to recover, solution(solver, 1) returns x and solution(solver, 2) returns y.

source
Krylov.results — Function
results(solver)

Return a tuple containing the solution(s) and the statistics associated with the solver. Allows retrieving the output arguments of an out-of-place method from the in-place method.

For example, instead of x, stats = cg(A, b), you can use:

    solver = CgSolver(A, b)
    cg!(solver, A, b)
    x, stats = results(solver)
source

Examples

We illustrate the use of in-place Krylov solvers with two well-known optimization methods. The details of the optimization methods are described in the section about Factorization-free operators.

Example 1: Newton's method for convex optimization without linesearch

using Krylov

function newton(∇f, ∇²f, x₀; itmax = 200, tol = 1e-8)

    n = length(x₀)
    x = copy(x₀)
    gx = ∇f(x)
    
    iter = 0
    S = typeof(x)
    solver = CgSolver(n, n, S)
    Δx = solver.x

    solved = false
    tired = false

    while !(solved || tired)
 
        Hx = ∇²f(x)           # Compute ∇²f(xₖ)
        cg!(solver, Hx, -gx)  # Solve ∇²f(xₖ)Δx = -∇f(xₖ)
        x = x + Δx            # Update xₖ₊₁ = xₖ + Δx
        gx = ∇f(x)            # ∇f(xₖ₊₁)
        
        iter += 1
        solved = norm(gx) ≤ tol
        tired = iter ≥ itmax
    end
    return x
end

Example 2: The Gauss-Newton method for nonlinear least squares without linesearch

using Krylov

function gauss_newton(F, JF, x₀; itmax = 200, tol = 1e-8)

    n = length(x₀)
    x = copy(x₀)
    Fx = F(x)
    m = length(Fx)
    
    iter = 0
    S = typeof(x)
    solver = LsmrSolver(m, n, S)
    Δx = solver.x

    solved = false
    tired = false

    while !(solved || tired)
 
        Jx = JF(x)              # Compute J(xₖ)
        lsmr!(solver, Jx, -Fx)  # Minimize ‖J(xₖ)Δx + F(xₖ)‖
        x = x + Δx              # Update xₖ₊₁ = xₖ + Δx
        Fx_old = Fx             # F(xₖ)
        Fx = F(x)               # F(xₖ₊₁)
        
        iter += 1
        solved = norm(Fx - Fx_old) / norm(Fx) ≤ tol
        tired = iter ≥ itmax
    end
    return x
end