Internal API Documentation

SciMLLinearSolveAlgorithm <: SciMLBase.AbstractLinearAlgorithm

The root abstract type for all linear solver algorithms in LinearSolve.jl. All concrete linear solver implementations should inherit from one of the specialized subtypes rather than directly from this type.

This type integrates with the SciMLBase ecosystem, providing a consistent interface for linear algebra operations across the Julia scientific computing ecosystem.

Interface

A concrete MyAlg <: SciMLLinearSolveAlgorithm must implement SciMLBase.solve!(cache::LinearCache, alg::MyAlg; kwargs...) and needs_concrete_A(alg::MyAlg)::Bool. It may implement init_cacheval, default_alias_A, default_alias_b, needs_square_A, and update_tolerances_internal!; each of these has a documented default.

Subtyping one of the categorized abstract types below supplies needs_concrete_A and the aliasing defaults. Direct subtypes must define needs_concrete_A themselves. The complete contract, including cache lifecycle rules and extension boundaries, is documented on the Linear Solver Algorithm Interface page.

Extension rules

Define the four traits needs_concrete_A, needs_square_A, default_alias_A, and default_alias_b in the package that defines MyAlg, next to the algorithm type. Downstream solvers query them before an optional backend is loaded. Methods that actually call the backend, such as solve! and init_cacheval, may be defined in the package extension instead.

Examples

struct MyAlg <: LinearSolve.AbstractKrylovSubspaceMethod end
LinearSolve.needs_square_A(::MyAlg) = true

Use algorithm_interface_issues to check a custom algorithm before passing it to init or solve.

AbstractFactorization <: SciMLLinearSolveAlgorithm

Abstract type for linear solvers that work by computing a matrix factorization. These algorithms typically decompose the matrix A into a product of simpler matrices (e.g., A = LU, A = QR, A = LDL') and then solve the system using forward/backward substitution.

Interface

Factorization algorithms receive a concrete representation of A, normally factorize it when cache.isfresh is true, store the factorization in cache.cacheval, and return a SciMLBase.build_linear_solution from solve!. The concrete subtypes are AbstractDenseFactorization and AbstractSparseFactorization.

Examples

LUFactorization, QRFactorization, CholeskyFactorization, UMFPACKFactorization, and KLUFactorization are concrete subtypes.

AbstractDenseFactorization <: AbstractFactorization

Abstract type for factorization-based linear solvers optimized for dense matrices. These algorithms assume the matrix has no particular sparsity structure and use dense linear algebra routines (typically from BLAS/LAPACK) for optimal performance.

Interface

Use this supertype when the algorithm needs the entries of a dense A to build a factorization. The default needs_concrete_A is true and the default aliasing traits are false, preserving the caller's matrix while a factorization is built.

Examples

LUFactorization, QRFactorization, CholeskyFactorization, and BunchKaufmanFactorization are concrete subtypes.

AbstractSparseFactorization <: AbstractFactorization

Abstract type for factorization-based linear solvers optimized for sparse matrices. These algorithms take advantage of sparsity patterns to reduce memory usage and computational cost compared to dense factorizations.

Interface

Use this supertype when the algorithm can preserve and exploit sparse structure. The default needs_concrete_A, default_alias_A, and default_alias_b traits are appropriate for the usual non-mutating sparse factorization workflow. An algorithm that mutates its input must override the aliasing traits and document that requirement.

Examples

UMFPACKFactorization, KLUFactorization, CHOLMODFactorization, SparspakFactorization, and ParUFactorization are concrete subtypes.

AbstractKrylovSubspaceMethod <: SciMLLinearSolveAlgorithm

Abstract type for iterative linear solvers based on Krylov subspace methods. These algorithms solve linear systems by iteratively building an approximation from a sequence of Krylov subspaces, without requiring explicit matrix factorization.

Interface

Use this supertype for algorithms that can solve from matrix-vector products without materializing A. needs_concrete_A defaults to false. The implementation reads the right-hand side and tolerances from LinearCache, uses cache.Pl and cache.Pr when it supports preconditioning, writes into cache.u, and returns a SciMLBase.build_linear_solution.

Examples

Krylov wrappers such as KrylovJL_GMRES, KrylovJL_CG, and IterativeSolversJL_GMRES are concrete subtypes. A matrix-free implementation must support the mul! operations required by its algorithm.

AbstractSolveFunction <: SciMLLinearSolveAlgorithm

Abstract type for linear solvers that wrap custom solving functions or provide direct interfaces to specific solve methods. These provide flexibility for integrating custom algorithms or simple solve strategies.

Interface

Use this supertype for an algorithm that delegates solving to a callable or a specialized direct operation. needs_concrete_A defaults to false, but the wrapped implementation is responsible for accepting the operator types that the algorithm advertises. The callable must return a solution compatible with cache.u.

Examples

LinearSolveFunction wraps a user callable, while DirectLdiv! delegates to Julia's ldiv! implementation.

LinearCache{TA, Tb, Tu, Tp, Talg, Tc, Tl, Tr, Ttol, issq, S}

The mutable state passed to a linear solver algorithm by init and reused by solve!. Construct it with SciMLBase.init(::LinearProblem, alg) rather than calling the constructor directly.

Fields

• A::TA: Operator or matrix for the system.
• b::Tb: Right-hand side. It may be a vector or a matrix of right-hand sides.
• u::Tu: Preallocated solution storage written by solve!.
• p::Tp: Problem parameters forwarded to the algorithm.
• alg::Talg: Algorithm instance used by this cache.
• cacheval::Tc: Algorithm-owned factorization, workspace, or solver object.
• isfresh::Bool: Whether cacheval must be rebuilt because A changed.
• precsisfresh::Bool: Whether the preconditioners must be refreshed.
• Pl::Tl: Left preconditioner, or nothing.
• Pr::Tr: Right preconditioner, or nothing.
• abstol::Ttol: Absolute convergence tolerance.
• reltol::Ttol: Relative convergence tolerance.
• maxiters::Int: Maximum iteration count for iterative algorithms.
• verbose::Tlv: Verbosity specification.
• assumptions::OperatorAssumptions{issq}: Properties promised about A.
• sensealg::S: Sensitivity algorithm associated with the solve.
• sparse_reduction::Tred: State for persistent sparse-pattern reduction, or nothing.
• alias_A::Bool: Whether the caller permits replacing or mutating A.

Interface rules

An algorithm's solve! method may update u, cacheval, and the freshness flags, but must preserve the meaning of the other fields. When isfresh is true, rebuild any factorization or backend object that depends on A; after doing so, set it to false. Algorithms that read tolerances at solve time use abstol, reltol, and maxiters directly. Algorithms that copy tolerances into cacheval must implement update_tolerances_internal!.

Examples

prob = LinearProblem(A, b)
cache = init(prob, LUFactorization())
sol = solve!(cache)
cache.b = b2
sol2 = solve!(cache)

init_cacheval(alg::SciMLLinearSolveAlgorithm, args...)

Initialize algorithm-specific cache values for the given linear solver algorithm. This function returns nothing by default and is intended to be overloaded by specific algorithm implementations that need to store intermediate computations or factorizations.

Arguments

• alg: The linear solver algorithm instance
• args...: Additional arguments passed to the cache initialization

Returns

Algorithm-specific cache value or nothing for algorithms that don't require caching.

cache.cacheval = NamedTuple(LUFactorization = cache of LUFactorization, ...)

defaultalg(A, b, assumptions::OperatorAssumptions)

Select a default linear solver algorithm for the operator A, right-hand side b, and operator assumptions. This is the dispatch point used by solve(::LinearProblem) when no algorithm is supplied explicitly.

Arguments

• A: Matrix, factorization, or abstract operator to solve with.
• b: Right-hand side vector or matrix.
• assumptions: OperatorAssumptions describing whether A is square, its conditioning, and its structural-zero behavior.

Returns

A concrete SciMLLinearSolveAlgorithm, usually wrapped in a DefaultLinearSolver, selected for the input representation and available extensions.

Examples

A = rand(100, 100)
b = rand(100)
alg = defaultalg(A, b, OperatorAssumptions(true))
solve(LinearProblem(A, b), alg)

For an abstract matrix-free operator, the default is a Krylov algorithm when the operator does not provide a direct solve method:

A = SciMLOperators.MatrixOperator(rand(10, 10))
b = rand(10)
alg = defaultalg(A, b, OperatorAssumptions(true))

Notes

The two-argument form assumes a square system. Use an explicit algorithm when the application requires a particular factorization or iterative method.

get_tuned_algorithm(::Type{eltype_A}, ::Type{eltype_b}, matrix_size) where {eltype_A, eltype_b}

Get the tuned algorithm preference for the given element type and matrix size. Returns nothing if no preference exists. Uses preloaded constants for efficiency. Fast path when no preferences are set.

is_algorithm_available(alg::DefaultAlgorithmChoice.T)

Check if the given algorithm is currently available (extensions loaded, etc.).

show_algorithm_choices()

Print a report of the dense default algorithm selection to stdout and return nothing. Takes no arguments. Use it to check whether autotune preferences are set and took effect, or to see which LU variant solve(prob) will pick on this machine.

The report has three parts, followed by a short legend of the size categories:

• Current Preferences: the Preferences.jl entries best_algorithm_<eltype>_<size> and best_always_loaded_<eltype>_<size> stored for LinearSolve (typically written to LocalPreferences.toml by LinearSolveAutotune.jl) for the eltypes Float32, Float64, ComplexF32, ComplexF64 and the size categories tiny (n <= 20), small (21-100), medium (101-300), large (301-1000), big (> 1000), or "No autotune preferences currently set." if there are none.
• Default Algorithm Choices: a table of the algorithm defaultalg returns for a dense random square matrix (with OperatorAssumptions(true)) of each eltype at one representative size per category: 8, 50, 200, 500, and 1500. Sizes of 10 or less always resolve to GenericLUFactorization.
• System Information: whether MKL and Apple Accelerate are available and whether RecursiveFactorization.jl is enabled.

Preferences are read into constants when LinearSolve is compiled, so a preference set in the current session shows up in the first section but only changes the second section (and actual solves) after Julia is restarted.

make_preferences_dynamic!()

Internal function for testing only. Makes preferences dynamic by redefining gettunedalgorithm to check preferences at runtime instead of using compile-time constants. This allows tests to verify that the preference system works correctly.

supernodal_panel_solve!(W, B, np; operation, algorithm = :auto)

Apply a supernodal triangular-panel operation using the requested backend.

Arguments

• W: Matrix containing the factored diagonal block.
• B: Panel or right-hand side to update in place.
• np: Width of the diagonal block in W.

Keywords

• operation: which triangular operation to apply. One of

  • :factor_right_upper: B := B / U11, a right solve with the non-unit upper triangle (used during numeric factorization to form L21).
  • :factor_lower: B := L11 \ B, a left solve with the unit-lower triangle (used during numeric factorization to form U12).
  • :lower: B := L11 \ B, the unit-lower forward substitution of the solve phase.
  • :upper: B := U11 \ B, the non-unit upper back substitution of the solve phase.

  Any other symbol throws an ArgumentError. :factor_lower and :lower compute the same result; they differ only in which backend code path each algorithm routes them to (see below).

• algorithm: backend to dispatch to through supernodal_panel_solve_backend!. Defaults to :auto. Accepted values:

  • :kernel: the in-tree column-oriented kernels (generic over the element type, allocation-free) for :lower/:upper; the :factor_* operations are forwarded to :triangularsolve.
  • :blas: BLAS.trsm! for :lower/:upper when W and B are strided BlasFloat matrices, the kernels for other element types; the :factor_* operations are forwarded to :triangularsolve.
  • :triangularsolve: LinearAlgebra.ldiv!/rdiv! on the triangular wrappers, replaced by TriangularSolve.jl kernels for strided Float32 and Float64 panels when RecursiveFactorization.jl and TriangularSolve.jl are loaded.
  • :auto: :factor_right_upper and :factor_lower always go to :triangularsolve. :lower and :upper use :kernel when np <= PANEL_KERNEL_MAX_NP (256), when B has a single column, or when the element type is not a BlasFloat; otherwise :blas when np > PANEL_BLAS_MIN_NP (1792) and :triangularsolve in between.

Returns

The updated B.

supernodal_panel_solve_backend!(algorithm, W, B, np; operation)

Backend extension hook for supernodal_panel_solve!.

Interface rules

An extension may specialize algorithm::Val for supported operand types. The method must apply operation to B in place, return B, and treat the factored diagonal block W[1:np, 1:np] as read-only. B may be a view into the remainder of W. Unsupported operations must throw ArgumentError.

The built-in backends use Val(:kernel), Val(:blas), and Val(:triangularsolve). Extensions should add methods only for backend and operand combinations they implement; the generic methods provide the fallback behavior.

Built-in backends

• Val(:kernel): runs :lower/:upper through the in-tree column-oriented kernels; :factor_right_upper/:factor_lower are forwarded to Val(:triangularsolve).
• Val(:blas): runs :lower/:upper through BLAS.trsm! when W and B are strided BlasFloat matrices, and through the kernels otherwise; :factor_right_upper/:factor_lower are forwarded to Val(:triangularsolve).
• Val(:triangularsolve): the in-tree method runs :factor_right_upper through LinearAlgebra.rdiv! with UpperTriangular and :factor_lower through LinearAlgebra.ldiv! with UnitLowerTriangular, and forwards :lower/:upper to Val(:blas).

This is the extension hook: a package can add a more specific method for Val(:triangularsolve) and supported panel types. The LinearSolveRecursiveFactorizationExt extension defines

supernodal_panel_solve_backend!(::Val{:triangularsolve}, W::StridedMatrix{Tv},
    B::StridedMatrix{Tv}, np::Int; operation::Symbol) where {Tv <: Union{Float32, Float64}}

which runs all four operations through TriangularSolve.jl and is active once RecursiveFactorization.jl and TriangularSolve.jl are both loaded (using RecursiveFactorization is enough, since it depends on TriangularSolve.jl).

AMD

Direct pure-Julia port of SuiteSparse's AMD (Approximate Minimum Degree) ordering algorithm. Matches the SuiteSparse implementation closely enough to produce identical permutations on the same inputs.

The entry point is amd_order!; supporting routines mirror amd_aat, amd_1, amd_2, amd_postorder and amd_post_tree from SuiteSparse.

amd_order!(n, Ap, Ai, P; dense_alpha=10.0, aggressive=true) -> (status, lnz)

Compute the AMD ordering of the symmetric pattern of A+A'. P is the output permutation. Mirrors amd_order.c. Returns (AMD_OK, lnz) on success, where lnz is the number of off-diagonal nonzeros in L (SuiteSparse's Info[AMD_LNZ]). On failure returns (status, 0.0).

needs_concrete_A(alg) -> Bool

Trait function that determines whether a linear solver algorithm requires a concrete matrix representation or can work with abstract operators.

Arguments

• alg: A linear solver algorithm instance

Returns

• true: Algorithm requires a concrete matrix (e.g., for factorization)
• false: Algorithm can work with abstract operators (e.g., matrix-free methods)

Usage

This trait is used internally by LinearSolve.jl to optimize algorithm dispatch and determine when matrix operators need to be converted to concrete arrays. It is also queried by downstream solvers such as OrdinaryDiffEq.jl and NonlinearSolve.jl to decide whether to assemble a concrete Jacobian, which is why every algorithm must implement it, and why it must be implemented next to the algorithm struct rather than in a package extension: the callers run before the backend package is necessarily loaded.

Algorithm-Specific Behavior

• AbstractFactorization: true (needs explicit matrix entries for factorization)
• AbstractKrylovSubspaceMethod: false (only needs matrix-vector products)
• AbstractSolveFunction: false (depends on the wrapped function's requirements)
• Direct subtypes of SciMLLinearSolveAlgorithm: no default; defining this is required

Example

needs_concrete_A(LUFactorization())  # true
needs_concrete_A(GMRESIteration())   # false

default_tol(T)

Compute the default tolerance for iterative linear solvers based on the element type. The tolerance is typically set as the square root of the machine epsilon for the given floating point type, ensuring numerical accuracy appropriate for that precision.

Arguments

• T: The element type of the linear system

Returns

• For floating point types: √(eps(T))
• For exact types (Rational, Integer): 0 (exact arithmetic)
• For Any type: 0 (conservative default)

default_alias_A(alg, A, b) -> Bool

Determine the default aliasing behavior for the matrix A given the algorithm type. Aliasing allows the algorithm to modify the original matrix in-place for efficiency, but this may not be desirable or safe for all algorithm types.

Arguments

• alg: The linear solver algorithm
• A: The matrix operator
• b: The right-hand side vector

Returns

• false: Safe default, algorithm will not modify the original matrix A
• true: Algorithm may modify A in-place for efficiency

Algorithm-Specific Behavior

• Dense factorizations: false (destructive, need to preserve original)
• Krylov methods: true (non-destructive, safe to alias)
• Sparse factorizations: true (typically preserve sparsity structure)

default_alias_b(alg, A, b) -> Bool

Determine the default aliasing behavior for the right-hand side vector b given the algorithm type. Similar to default_alias_A but for the RHS vector.

Returns

• false: Safe default, algorithm will not modify the original vector b
• true: Algorithm may modify b in-place for efficiency

__init_u0_from_Ab(A, b)

Initialize the solution vector u0 with appropriate size and type based on the matrix A and right-hand side b. The solution vector is allocated with the same element type as b and sized to match the number of columns in A.

Arguments

• A: The matrix operator (determines solution vector size)
• b: The right-hand side vector (determines element type)

Returns

A zero-initialized vector of size (size(A, 2),) with element type matching b. For a matrix (batched) right-hand side b of size (size(A, 1), k), returns a zero-initialized matrix of size (size(A, 2), k) so that each column of u0 corresponds to a column of b.

Specializations

• For static matrices (SMatrix): Returns a static vector (SVector)
• For regular matrices: Returns a similar vector to b with appropriate size

solve!(cache::LinearCache, args...; adjoint = false, kwargs...)

adjoint = true solves adjoint(A) x = b instead of A x = b, reusing the factorization the cache already holds rather than factorizing the adjoint afresh:

cache = init(LinearProblem(A, b))
x = solve!(cache).u
cache.b = c
lambda = solve!(cache; adjoint = true).u   # adjoint(A) lambda = c, same factorization

Algorithms that have not opted into reusing their factorization for an adjoint solve fall back to factorizing adjoint(A), so the keyword is always answerable. The solution is written into cache.u as usual.

ComposePreconditioner{Ti, To}

A preconditioner that composes two preconditioners by applying them sequentially. The inner preconditioner is applied first, followed by the outer preconditioner. This allows for building complex preconditioning strategies by combining simpler ones.

Fields

• inner::Ti: The inner (first) preconditioner to apply
• outer::To: The outer (second) preconditioner to apply

Usage

# Compose a diagonal preconditioner with an ILU preconditioner
inner_prec = DiagonalPreconditioner(diag(A))
outer_prec = ILUFactorization()  
composed = ComposePreconditioner(inner_prec, outer_prec)

The composed preconditioner applies: outer(inner(x)) for any vector x.

Mathematical Interpretation

For a linear system Ax = b, if P₁ is the inner and P₂ is the outer preconditioner, then the composed preconditioner effectively applies P₂P₁ as the combined preconditioner.

InvPreconditioner{T}

A preconditioner wrapper that treats a matrix or operator as if it represents the inverse of the actual preconditioner. Instead of solving Px = y, it computes P*y where P is stored as the "inverse" preconditioner matrix.

Fields

• P::T: The stored preconditioner matrix/operator (representing P⁻¹)

Usage

This is useful when you have a matrix that approximates the inverse of your desired preconditioner. For example, if you have computed an approximate inverse matrix Ainv ≈ A⁻¹, you can use:

prec = InvPreconditioner(Ainv)

Mathematical Interpretation

For a linear system Ax = b with preconditioner M, normally we solve M⁻¹Ax = M⁻¹b. With InvPreconditioner, the stored matrix P represents M⁻¹ directly, so applying the preconditioner becomes a matrix-vector multiplication rather than a linear solve.

Methods

• ldiv!(A::InvPreconditioner, x): Computes x ← P*x (in-place)
• ldiv!(y, A::InvPreconditioner, x): Computes y ← P*x
• mul!(y, A::InvPreconditioner, x): Computes y ← P⁻¹*x (inverse operation)