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 two constructors:

XyzSolver(A, b)
XyzSolver(m, n, S)

Xyz is the name of the Krylov method with lowercase letters except its first one (Cg, Minres, Lsmr, Bicgstab, ...). 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 give the dimensions (m, n) and the storage type S.

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