Krylov.jl Documentation

Krylov.jl provides implementations for several classes of linear systems, making it suitable for scenarios where matrix factorization is impossible or inefficient (e.g., when the operator $A$ is not explicitly available, is too large to materialize, or is a 'fast' operator).

Supported problem types include:

1. Square or rectangular full-rank systems ($Ax = b$): Used when $b$ lies in the range space of $A$ (e.g., $A$ is square/nonsingular or tall/full column rank).
2. Linear least-squares problems ($\min |b - Ax|$): Used for inconsistent systems where $b$ is not in the range of $A$ (e.g., $A$ is square/singular or tall/thin). If $A$ is column rank-deficient, it identifies the minimum norm solution.
3. Linear least-norm problems ($\min |x|$ subject to $Ax = b$): Used when $A$ is column rank-deficient but $b$ is in the range of $A$ (e.g., $A$ is square/singular or short/wide).
4. Adjoint systems ($Ax = b$ and $A^H y = c$): $A$ can have any shape.
5. Saddle-point and Hermitian quasi-definite systems: Solves systems of the form $\begin{bmatrix} M & A \ A^H & -N \ \end{bmatrix} \begin{bmatrix} x \ y \end{bmatrix} = \dots$
6. Generalized saddle-point and non-Hermitian partitioned systems: Solves systems of the form $\begin{bmatrix} M & A \ B & N \end{bmatrix} \begin{bmatrix} x \ y \end{bmatrix} = \begin{bmatrix} b \ c \end{bmatrix}$, where $B$ has the shape of $A^H$.

Key Features:

• All solvers include in-place versions.
• GPU compatible.
• Supports any floating-point data type.

A HaloVector is a specialized AbstractVector designed for finite difference methods that use halo regions (ghost cells). By embedding an OffsetArray (from OffsetArrays.jl) within the HaloVector, you can enable "if-less" stencil applications, where boundary conditions are handled by the data layout rather than explicit conditional checks in the core loop.

To implement a custom HaloVector, you must provide:

• Base.similar(v::HaloVector): Returns a new HaloVector with the same element type and similar data structure.
• Base.length(v::HaloVector): Returns the number of interior (non-halo) elements.
• Base.size(v::HaloVector): Returns the dimensions of the interior vector.
• Base.getindex(v::HaloVector, idx): Maps the linear index to the underlying 2D/3D grid coordinates.

using OffsetArrays

struct HaloVector{FC, D} <: AbstractVector{FC}
    data::D

    function HaloVector(data::D) where {D}
        FC = eltype(data)
        return new{FC, D}(data)
    end
end

function Base.similar(v::HaloVector)
    data = similar(v.data)
    fill!(data, zero(eltype(data)))
    return HaloVector(data)
end

function Base.length(v::HaloVector)
    m, n = size(v.data)
    l = (m - 2) * (n - 2)
    return l
end

function Base.size(v::HaloVector)
    l = length(v)
    return (l,)
end

function Base.getindex(v::HaloVector, idx)
    m, n = size(v.data)
    row = div(idx - 1, n - 2) + 1
    col = mod(idx - 1, n - 2) + 1
    return v.data[row, col]
end

Selecting a solver from KrylovSolverType requires understanding three dimensions of compatibility:

1. Adjoint Requirements (matvec_At)

Some solvers require the adjoint operator. If your solver is in the 'Provide matvec_At' list, you must pass a valid callback to the third argument of krylov_solve.

• No Adjoint Needed: CG, CR, CAR, MINRES, MINRES-QLP, MINARES, SYMMLQ, GMRES, FGMRES, FOM, DIOM, DQGMRES, BiCGSTAB, CGS.
• Adjoint Required: BiLQ, QMR, BiLQR, TriLQR, USYMLQ, USYMLQR, TriCG, TriMR, LSLQ, LSQR, LSMR, CGLS, CRLS, CGNE, CRMR, CRAIG, CRAIGMR, LNLQ, GPMR (uses $B$ instead of $A^H$).

2. Second Right-Hand Side (c)

Some solvers solve for two right-hand sides simultaneously. For these, you must provide a second vector c to krylov_solve.

• Requires c: TriCG, TriMR, BiLQR, TriLQR, USYMLQ, USYMLQR, and GPMR.

3. Dual Solution Retrieval (krylov_get_y)

If a solver produces a second solution $y$, you can retrieve it using krylov_get_y. For all other solvers, this function returns -2.

• Produces y: TriCG, TriMR, USYMLQR, GPMR, BiLQR, TriLQR, CRAIG, CRAIGMR, and LNLQ.

The libkrylov interface is designed around three core principles:

1. Opaque workspace handle: You manage memory using a void * workspace handle. You create it with krylov_workspace_create and must manually release it with krylov_workspace_free. The library handles the internal Julia-rooted memory.
2. Matrix-free via callbacks: Solvers do not take matrix objects. Instead, you provide C function pointers (callbacks) that perform matrix-vector products (e.g., A * x, Aᴴ * x, or M⁻¹ * x). This allows the library to work with dense, sparse, or stencil-based operators.
3. Static typing: Workspaces are specialized for a specific (solver, precision) pair at creation time, ensuring high performance without dynamic dispatch during the solve phase.

Standard Solver Lifecycle

Regardless of the solver chosen, the workflow is:

1. Create: krylov_workspace_create to allocate the workspace.
2. Solve: krylov_solve passing your matrix-vector product callbacks.
3. Retrieve: krylov_get_x (or krylov_get_y for solvers with two solutions) to get the result.
4. Free: krylov_workspace_free to clean up.

Krylov.jl solvers are matrix-free, meaning they only require operator-vector products rather than explicit matrices. To use a custom object as a linear operator (for A or B arguments), your object must implement:

• mul!(y, A, v): Multiplication with a vector.
• size(A): Returns the dimensions of the operator.
• eltype(A): Returns the element type.

Some methods (like CGLS, LSQR, GMRES, etc.) also require the adjoint operator. In these cases, you must implement adjoint(A) or provide a method for mul!(y, transpose(A), w)/mul!(y, A', u).

It is highly recommended to use LinearOperators.jl to model these operators. Using LinearOperator(type, nrows, ncols, symmetric, hermitian, prod, tprod, ctprod), you can define:

• prod(y, v): for mul!(y, A, v)
• tprod(y, w): for mul!(y, transpose(A), w)
• ctprod(u, w): for mul!(y, A', u)

# Example using LinearOperators.jl
A = LinearOperator(Float64, 4, 4, true, true, (y, v) -> H(y, v))

A preconditioner in Krylov.jl is a matrix-free callback that computes the action of the inverse: $y = M^{-1}x$. It must solve the system $My = x$, not compute $Mx$.

• Symmetric solvers (e.g., CG, CR, CAR, MINRES): Use a single centered preconditioner $M$. In krylov_solve, pass this to the third argument (matvec_M).
• Non-symmetric solvers (e.g., GMRES, BiCGSTAB): Use a left preconditioner $M$ and a right preconditioner $N$. In krylov_solve, pass $M$ to the third argument and $N$ to the fourth argument.

/* y = M⁻¹ x with M = diag(A): apply the INVERSE, i.e. divide by the diagonal. */
static void precond_M(const void *xv, void *yv, void *userdata)
{
  const double *x = (const double *)xv;
  double       *y = (double *)yv;
  const Tri    *A = (const Tri *)userdata;
  for (int i = 0; i < A->n; i++) y[i] = x[i] / A->diag[i];
}

// ... in main ...
// For symmetric solvers, matvec_M is the centered preconditioner; matvec_N is NULL.
krylov_solve(ws, matvec_A, NULL, precond_M, NULL, b, NULL, &A, &opts);