Linear Solver Algorithm Interface

Every algorithm that LinearSolve.jl can be handed as solve(prob, alg) is a subtype of LinearSolve.SciMLLinearSolveAlgorithm. This page is the contract that such a type has to satisfy. It applies equally to algorithms defined inside LinearSolve.jl, in one of its package extensions, and in downstream packages.

The interface is checked in LinearSolve.jl's test suite for every algorithm the package knows about, using LinearSolve.algorithm_interface_issues; downstream packages can run the same check against their own algorithms.

Choosing a supertype

SciMLLinearSolveAlgorithm
├── AbstractFactorization
│   ├── AbstractDenseFactorization
│   └── AbstractSparseFactorization
├── AbstractKrylovSubspaceMethod
└── AbstractSolveFunction

Subtype the most specific abstract type that fits. Doing so supplies several of the methods below; subtyping SciMLLinearSolveAlgorithm directly supplies none of them, so such an algorithm must define needs_concrete_A itself. Every trait default that comes from the categorized types is listed in the tables below.

Required methods

solve! is the algorithm proper. It reads cache.A, cache.b, cache.u, cache.Pl/cache.Pr, cache.abstol/cache.reltol/cache.maxiters, and writes the solution into cache.u. cache.isfresh is true when A has changed since the previous solve; algorithms that build a factorization or a solver object should do so only when it is set, store the result in cache.cacheval, and then set cache.isfresh = false.

needs_concrete_A is a pure function of the algorithm type. Downstream solvers — OrdinaryDiffEq.jl and NonlinearSolve.jl in particular — call it on a user-supplied linsolve to decide whether to assemble a concrete Jacobian.

Optional methods

Each of these has a default, so an algorithm only defines the ones whose default does not fit.

update_tolerances_internal! needs no definition when the algorithm has no tolerances at all. Define it as nothing when the algorithm reads cache.abstol/cache.reltol at solve time (Krylov methods do), and define it to write into cache.cacheval when the algorithm snapshots the tolerances into its own solver object at init time.

The four Bool-valued traits are called while init builds the cache, so they must be inferable as Bool — write them as literal returns (needs_concrete_A(::MyAlg) = true) rather than as runtime computations.

A complete example

using LinearSolve, LinearAlgebra, SciMLBase

struct MyLUFactorization <: LinearSolve.AbstractDenseFactorization end

# Optional: give the cache its final type up front so repeated solves stay
# type-stable. `ArrayInterface.lu_instance(A)` is the cheap way to do this.
function LinearSolve.init_cacheval(
        alg::MyLUFactorization, A, b, u, Pl, Pr, maxiters::Int, abstol, reltol,
        verbose, assump::LinearSolve.OperatorAssumptions
    )
    return lu(convert(AbstractMatrix, A))
end

function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::MyLUFactorization; kwargs...)
    if cache.isfresh
        cache.cacheval = lu!(convert(AbstractMatrix, cache.A))
        cache.isfresh = false
    end
    y = ldiv!(cache.u, cache.cacheval, cache.b)
    return SciMLBase.build_linear_solution(alg, y, nothing, cache)
end

Subtyping AbstractDenseFactorization supplies needs_concrete_A(alg) = true and the false aliasing defaults, which is what a destructive dense factorization wants. An algorithm that subtypes SciMLLinearSolveAlgorithm directly would additionally need

LinearSolve.needs_concrete_A(::MyLUFactorization) = true

Matrix-free implementation

An iterative algorithm can promise that it only needs the generic mul! operation. The operator does not need to expose a factorization or a dense matrix representation; the algorithm still receives the standard LinearCache and returns the standard LinearSolution.

struct DocIdentityOperator{T} <: AbstractMatrix{T}
    n::Int
end

Base.size(A::DocIdentityOperator) = (A.n, A.n)
Base.getindex(A::DocIdentityOperator{T}, i::Int, j::Int) where {T} =
    i == j ? one(T) : zero(T)

function LinearAlgebra.mul!(y::AbstractVector, ::DocIdentityOperator, x::AbstractVector)
    copyto!(y, x)
    return y
end

struct DocMatVecAlg <: LinearSolve.AbstractKrylovSubspaceMethod end

function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::DocMatVecAlg; kwargs...)
    mul!(cache.u, cache.A, cache.b)
    return SciMLBase.build_linear_solution(
        alg, cache.u, nothing, cache; retcode = SciMLBase.ReturnCode.Success
    )
end

A = DocIdentityOperator{Float64}(3)
b = [1.0, 2.0, 3.0]
cache = init(LinearProblem(A, b), DocMatVecAlg())
@assert solve!(cache).u == b

The example is intentionally limited to the generic contract: a real iterative algorithm must also define its convergence, stopping, and preconditioning rules and should exercise those rules with dedicated tests.

Checking compliance

using Test
@test isempty(LinearSolve.algorithm_interface_issues(MyLUFactorization()))

algorithm_interface_issues accepts either an instance or the type, and returns one message per violation. Pass check_solve = false to check only the traits, which is what LinearSolve.jl's own test suite does for algorithms whose solve! lives in an extension whose backend is not loaded in that test group.

To sweep every algorithm a session knows about:

for T in LinearSolve.concrete_algorithm_types()
    @test isempty(LinearSolve.algorithm_interface_issues(T))
end

LinearSolve.algorithm_interface_issues(alg; check_solve = true) -> Vector{String}

struct MyLUFactorization <: LinearSolve.SciMLLinearSolveAlgorithm end
LinearSolve.needs_concrete_A(::MyLUFactorization) = true
function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::MyLUFactorization; kwargs...)
    # ...
end

@test isempty(LinearSolve.algorithm_interface_issues(MyLUFactorization()))

LinearSolve.concrete_algorithm_types(T = SciMLLinearSolveAlgorithm)

algorithms = LinearSolve.concrete_algorithm_types()
@assert all(T -> T <: LinearSolve.SciMLLinearSolveAlgorithm, algorithms)

needs_square_A(alg)

LinearSolve.update_tolerances!(cache; abstol = nothing, reltol = nothing)

LinearSolve.update_tolerances_internal!(cache, alg, abstol, reltol)