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.

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.

Characteristics

• Requires concrete matrix representation (needs_concrete_A() = true)
• Typically efficient for multiple solves with the same matrix
• Generally provides high accuracy for well-conditioned problems
• Memory requirements depend on the specific factorization type

Subtypes

• AbstractDenseFactorization: For dense matrix factorizations
• AbstractSparseFactorization: For sparse matrix factorizations

Examples of concrete subtypes

• LUFactorization, QRFactorization, CholeskyFactorization
• UMFPACKFactorization, KLUFactorization

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.

Characteristics

• Optimized for matrices with few zeros or no sparsity structure
• Leverage highly optimized BLAS/LAPACK routines when available
• Generally provide excellent performance for moderately-sized dense problems
• Memory usage scales as O(n²) with matrix size

Examples of concrete subtypes

• LUFactorization: Dense LU with partial pivoting (via LAPACK)
• QRFactorization: Dense QR factorization for overdetermined systems
• CholeskyFactorization: Dense Cholesky for symmetric positive definite matrices
• BunchKaufmanFactorization: For symmetric indefinite matrices

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.

Characteristics

• Optimized for matrices with many zero entries
• Often use specialized pivoting strategies to preserve sparsity
• May reorder rows/columns to minimize fill-in during factorization
• Typically more memory-efficient than dense methods for sparse problems

Examples of concrete subtypes

• UMFPACKFactorization: General sparse LU with partial pivoting
• KLUFactorization: Sparse LU optimized for circuit simulation
• CHOLMODFactorization: Sparse Cholesky for positive definite systems
• SparspakFactorization: Envelope/profile method for sparse systems

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.

Characteristics

• Does not require concrete matrix representation (needs_concrete_A() = false)
• Only needs matrix-vector products A*v (can work with operators/functions)
• Memory usage typically O(n) or O(kn) where k << n
• Convergence depends on matrix properties (condition number, eigenvalue distribution)
• Often benefits significantly from preconditioning

Advantages

• Low memory requirements for large sparse systems
• Can handle matrix-free operators (functions that compute A*v)
• Often the only feasible approach for very large systems
• Can exploit matrix structure through specialized operators

Examples of concrete subtypes

• GMRESIteration: Generalized Minimal Residual method
• CGIteration: Conjugate Gradient (for symmetric positive definite systems)
• BiCGStabLIteration: Bi-Conjugate Gradient Stabilized
• Wrapped external iterative solvers (KrylovKit.jl, IterativeSolvers.jl)

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.

Characteristics

• Does not require concrete matrix representation (needs_concrete_A() = false)
• Provides maximum flexibility for custom solving strategies
• Can wrap external solver libraries or implement specialized algorithms
• Performance and stability depend entirely on the wrapped implementation

Examples of concrete subtypes

• LinearSolveFunction: Wraps arbitrary user-defined solve functions
• DirectLdiv!: Direct application of the \ operator

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

The core cache structure used by LinearSolve for storing and managing the state of linear solver computations. This mutable struct acts as the primary interface for iterative solving and caching of factorizations and intermediate results.

Fields

• A::TA: The matrix operator of the linear system.
• b::Tb: The right-hand side vector of the linear system.
• u::Tu: The solution vector (preallocated storage for the result).
• p::Tp: Parameters passed to the linear solver algorithm.
• alg::Talg: The linear solver algorithm instance.
• cacheval::Tc: Algorithm-specific cache storage for factorizations and intermediate computations.
• isfresh::Bool: Cache validity flag for the matrix A. false means cacheval is up-to-date with respect to A, true means cacheval needs to be updated.
• precsisfresh::Bool: Cache validity flag for preconditioners. false means Pl and Pr are up-to-date with respect to A, true means they need to be updated.
• Pl::Tl: Left preconditioner operator.
• Pr::Tr: Right preconditioner operator.
• abstol::Ttol: Absolute tolerance for iterative solvers.
• reltol::Ttol: Relative tolerance for iterative solvers.
• maxiters::Int: Maximum number of iterations for iterative solvers.
• verbose::LinearVerbosity: Whether to print verbose output during solving.
• assumptions::OperatorAssumptions{issq}: Assumptions about the operator properties.
• sensealg::S: Sensitivity analysis algorithm for automatic differentiation.

Usage

The LinearCache is typically created via init(::LinearProblem, ::SciMLLinearSolveAlgorithm) and then used with solve!(cache) for efficient repeated solves with the same matrix structure but potentially different right-hand sides or parameter values.

Cache Management

The cache automatically tracks when matrix A or parameters p change by setting the appropriate freshness flags. When solve! is called, stale cache entries are automatically recomputed as needed.

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, ...)

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()

Display what algorithm choices are actually made by the default solver for representative matrix sizes. Shows current preferences and system information.

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.

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.

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)

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.

Specializations

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

LinearSolveFunction{F} <: AbstractSolveFunction

A flexible wrapper that allows using custom functions as linear solver algorithms. This provides a way to integrate user-defined solving strategies into the LinearSolve.jl framework while maintaining compatibility with the caching and interface systems.

Fields

• solve_func::F: A callable that implements the custom linear solving logic

Function Signature

The wrapped function should have the signature:

solve_func(A, b, u, p, isfresh, Pl, Pr, cacheval; kwargs...)

Arguments to wrapped function

• A: The matrix operator of the linear system
• b: The right-hand side vector
• u: Pre-allocated solution vector (can be used as working space)
• p: Parameters passed to the solver
• isfresh: Boolean indicating if the matrix A has changed since last solve
• Pl: Left preconditioner operator
• Pr: Right preconditioner operator
• cacheval: Algorithm-specific cache storage
• kwargs...: Additional keyword arguments

Returns

The wrapped function should return a solution vector.

Example

function my_custom_solver(A, b, u, p, isfresh, Pl, Pr, cacheval; kwargs...)
    # Custom solving logic here
    return A \ b  # Simple example
end

alg = LinearSolveFunction(my_custom_solver)
sol = solve(prob, alg)

DirectLdiv! <: AbstractSolveFunction

A simple linear solver that directly applies the left-division operator (\) to solve the linear system. This algorithm calls ldiv!(u, A, b) which computes u = A \ b in-place.

Usage

alg = DirectLdiv!()
sol = solve(prob, alg)

Notes

• This is essentially a direct wrapper around Julia's built-in ldiv! function
• Suitable for cases where the matrix A has a natural inverse or factorization
• Performance depends on the specific matrix type and its ldiv! implementation
• No preconditioners or advanced numerical techniques are applied
• Best used for small to medium problems or when A has special structure

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)

SimpleLUFactorization(pivot::Bool = true)

A pure Julia LU factorization implementation without BLAS dependencies. This solver is optimized for small matrices and situations where BLAS is not available or desirable.

Constructor Arguments

• pivot::Bool = true: Whether to perform partial pivoting for numerical stability. Set to false for slightly better performance at the cost of stability.

Features

• Pure Julia implementation (no BLAS dependencies)
• Partial pivoting support for numerical stability
• In-place matrix modification for memory efficiency
• Fast for small matrices (typically < 100×100)
• Educational value for understanding LU factorization

Performance Characteristics

• Optimal for small dense matrices
• No overhead from BLAS calls
• Linear scaling with problem size (O(n³) operations)
• Memory efficient due to in-place operations

Use Cases

• Small matrices where BLAS overhead is significant
• Systems without optimized BLAS libraries
• Educational and prototyping purposes
• Embedded systems with memory constraints

Example

# Stable version with pivoting (default)
alg1 = SimpleLUFactorization()
# Faster version without pivoting
alg2 = SimpleLUFactorization(false)

prob = LinearProblem(A, b)
sol = solve(prob, alg1)

Notes

For larger matrices (> 100×100), consider using BLAS-based factorizations like LUFactorization() for better performance.

LUSolver{T}

A mutable workspace for performing LU factorization and solving linear systems. This struct maintains all necessary arrays and state information for the factorization and solve phases, allowing for efficient reuse when solving multiple systems with the same matrix structure.

Fields

• n::Int: Dimension of the square matrix
• A::Matrix{T}: Working copy of the matrix to be factorized (modified in-place)
• b::Vector{T}: Right-hand side vector storage
• x::Vector{T}: Solution vector storage
• pivots::Vector{Int}: Pivot indices from the factorization
• perms::Vector{Int}: Permutation vector tracking row exchanges
• info::Int: Status information (0 = success, >0 indicates singularity)

Constructor

LUSolver{T}(n)  # Create solver for n×n matrix with element type T

Usage

The solver is typically created from a matrix using the convenience constructors:

solver = LUSolver(A)        # From matrix A
solver = LUSolver(A, b)     # From matrix A and RHS b

Then factorized and solved:

simplelu_factorize!(solver)    # Perform LU factorization
simplelu_solve!(solver)        # Solve for the stored RHS

Notes

This is a pure Julia implementation primarily for educational purposes and small matrices. For production use, prefer optimized LAPACK-based factorizations.