Linear System Solvers

DefaultLinearSolver(;safetyfallback=true)

The default linear solver. This is the algorithm chosen when solve(prob) is called. It's a polyalgorithm that detects the optimal method for a given A, b and hardware (Intel, AMD, GPU, etc.).

Keyword Arguments

• safetyfallback: determines whether to fallback to a column-pivoted QR factorization when an LU factorization fails. This can be required if A is rank-deficient. Defaults to true.

RFLUFactorization{P, T}(; pivot = Val(true), thread = Val(true))

A fast pure Julia LU-factorization implementation using RecursiveFactorization.jl. This is by far the fastest LU-factorization implementation, usually outperforming OpenBLAS and MKL for smaller matrices (<500x500), but currently optimized only for Base Array with Float32 or Float64. Additional optimization for complex matrices is in the works.

Type Parameters

• P: Pivoting strategy as Val{Bool}. Val{true} enables partial pivoting for stability.
• T: Threading strategy as Val{Bool}. Val{true} enables multi-threading for performance.

Constructor Arguments

• pivot = Val(true): Enable partial pivoting. Set to Val{false} to disable for speed at the cost of numerical stability.
• thread = Val(true): Enable multi-threading. Set to Val{false} for single-threaded execution.
• throwerror = true: Whether to throw an error if RecursiveFactorization.jl is not loaded.

Performance Notes

• Fastest for dense matrices with dimensions roughly < 500×500
• Optimized specifically for Float32 and Float64 element types
• Recursive blocking strategy provides excellent cache performance
• Multi-threading can provide significant speedups on multi-core systems

Requirements

Using this solver requires that RecursiveFactorization.jl is loaded: using RecursiveFactorization

Example

using RecursiveFactorization
# Fast, stable (with pivoting)
alg1 = RFLUFactorization()
# Fastest (no pivoting), less stable
alg2 = RFLUFactorization(pivot=Val(false))  

LUFactorization(pivot=LinearAlgebra.RowMaximum())

Julia's built in lu. Equivalent to calling lu!(A)

• On dense matrices, this uses the current BLAS implementation of the user's computer, which by default is OpenBLAS but will use MKL if the user does using MKL in their system.
• On sparse matrices, this will use UMFPACK from SparseArrays. Note that this will not cache the symbolic factorization.
• On CuMatrix, it will use a CUDA-accelerated LU from CuSolver.
• On BandedMatrix and BlockBandedMatrix, it will use a banded LU.

Positional Arguments

• pivot: The choice of pivoting. Defaults to LinearAlgebra.RowMaximum(). The other choice is LinearAlgebra.NoPivot().

GenericLUFactorization(pivot=LinearAlgebra.RowMaximum())

Julia's built in generic LU factorization. Equivalent to calling LinearAlgebra.generic_lufact!. Supports arbitrary number types but does not achieve as good scaling as BLAS-based LU implementations. Has low overhead and is good for small matrices.

Positional Arguments

• pivot: The choice of pivoting. Defaults to LinearAlgebra.RowMaximum(). The other choice is LinearAlgebra.NoPivot().

QRFactorization(pivot=LinearAlgebra.NoPivot(),blocksize=16)

Julia's built in qr. Equivalent to calling qr!(A).

• On dense matrices, this uses the current BLAS implementation of the user's computer which by default is OpenBLAS but will use MKL if the user does using MKL in their system.
• On sparse matrices, this will use SPQR from SparseArrays
• On CuMatrix, it will use a CUDA-accelerated QR from CuSolver.
• On BandedMatrix and BlockBandedMatrix, it will use a banded QR.

SVDFactorization(full=false,alg=LinearAlgebra.DivideAndConquer())

Julia's built in svd. Equivalent to svd!(A).

• On dense matrices, this uses the current BLAS implementation of the user's computer which by default is OpenBLAS but will use MKL if the user does using MKL in their system.

CholeskyFactorization(; pivot = nothing, tol = 0.0, shift = 0.0, perm = nothing)

Julia's built in cholesky. Equivalent to calling cholesky!(A).

Keyword Arguments

• pivot: defaluts to NoPivot, can also be RowMaximum.
• tol: the tol argument in CHOLMOD. Only used for sparse matrices.
• shift: the shift argument in CHOLMOD. Only used for sparse matrices.
• perm: the perm argument in CHOLMOD. Only used for sparse matrices.

BunchKaufmanFactorization(; rook = false)

Julia's built in bunchkaufman. Equivalent to calling bunchkaufman(A). Only for Symmetric matrices.

Keyword Arguments

• rook: whether to perform rook pivoting. Defaults to false.

CHOLMODFactorization(; shift = 0.0, perm = nothing)

A wrapper of CHOLMOD's polyalgorithm, mixing Cholesky factorization and ldlt. Tries cholesky for performance and retries ldlt if conditioning causes Cholesky to fail.

Only supports sparse matrices.

Keyword Arguments

• shift: the shift argument in CHOLMOD.
• perm: the perm argument in CHOLMOD

NormalCholeskyFactorization(pivot = RowMaximum())

A fast factorization which uses a Cholesky factorization on A * A'. Can be much faster than LU factorization, but is not as numerically stable and thus should only be applied to well-conditioned matrices.

Positional Arguments

• pivot: Defaults to RowMaximum(), but can be NoPivot()

NormalBunchKaufmanFactorization(rook = false)

A fast factorization which uses a BunchKaufman factorization on A * A'. Can be much faster than LU factorization, but is not as numerically stable and thus should only be applied to well-conditioned matrices.

Positional Arguments

• rook: whether to perform rook pivoting. Defaults to false.

DiagonalFactorization()

A special implementation only for solving Diagonal matrices fast.

SimpleGMRES(; restart::Bool = true, blocksize::Int = 0, warm_start::Bool = false,
    memory::Int = 20)

A simple GMRES implementation for square non-Hermitian linear systems.

This implementation handles Block Diagonal Matrices with Uniformly Sized Square Blocks with specialized dispatches.

Arguments

• restart::Bool: If true, then the solver will restart after memory iterations.

• memory::Int = 20: The number of iterations before restarting. If restart is false, this value is used to allocate memory and later expanded if more memory is required.

• blocksize::Int = 0: If blocksize is > 0, the solver assumes that the matrix has a uniformly sized block diagonal structure with square blocks of size blocksize. Misusing this option will lead to incorrect results.

  • If this is set ≤ 0 and during runtime we get a Block Diagonal Matrix, then we will check if the specialized dispatch can be used.

FastLUFactorization()

A high-performance LU factorization using the FastLapackInterface.jl package. This provides an optimized interface to LAPACK routines with reduced overhead compared to the standard LinearAlgebra LAPACK wrappers.

Features

• Reduced function call overhead compared to standard LAPACK wrappers
• Optimized for performance-critical applications
• Uses partial pivoting (no choice of pivoting method available)
• Suitable for dense matrices where maximum performance is required

Limitations

• Does not allow customization of pivoting strategy (always uses partial pivoting)
• Requires FastLapackInterface.jl to be loaded
• Limited to dense matrix types supported by LAPACK

Requirements

Using this solver requires that FastLapackInterface.jl is loaded: using FastLapackInterface

Performance Notes

This factorization is optimized for cases where the overhead of standard LAPACK function calls becomes significant, typically for moderate-sized dense matrices or when performing many factorizations.

Example

using FastLapackInterface
alg = FastLUFactorization()
sol = solve(prob, alg)

FastQRFactorization{P}(; pivot = ColumnNorm(), blocksize = 36)

A high-performance QR factorization using the FastLapackInterface.jl package. This provides an optimized interface to LAPACK QR routines with reduced overhead compared to the standard LinearAlgebra LAPACK wrappers.

Type Parameters

• P: The type of pivoting strategy used

Fields

• pivot::P: Pivoting strategy (e.g., ColumnNorm() for column pivoting, nothing for no pivoting)
• blocksize::Int: Block size for the blocked QR algorithm (default: 36)

Features

• Reduced function call overhead compared to standard LAPACK wrappers
• Supports various pivoting strategies for numerical stability
• Configurable block size for optimal performance
• Suitable for dense matrices, especially overdetermined systems

Performance Notes

The block size can be tuned for optimal performance depending on matrix size and architecture. The default value of 36 is generally good for most cases, but experimentation may be beneficial for specific applications.

Requirements

Using this solver requires that FastLapackInterface.jl is loaded: using FastLapackInterface

Example

using FastLapackInterface
# QR with column pivoting
alg1 = FastQRFactorization()  
# QR without pivoting for speed
alg2 = FastQRFactorization(pivot=nothing)
# Custom block size
alg3 = FastQRFactorization(blocksize=64)

KLUFactorization(;reuse_symbolic=true, check_pattern=true)

A fast sparse LU-factorization which specializes on sparsity patterns with “less structure”.

UMFPACKFactorization(;reuse_symbolic=true, check_pattern=true)

A fast sparse multithreaded LU-factorization which specializes on sparsity patterns with “more structure”.

SparspakFactorization(reuse_symbolic = true)

This is the translation of the well-known sparse matrix software Sparspak (Waterloo Sparse Matrix Package), solving large sparse systems of linear algebraic equations. Sparspak is composed of the subroutines from the book "Computer Solution of Large Sparse Positive Definite Systems" by Alan George and Joseph Liu. Originally written in Fortran 77, later rewritten in Fortran 90. Here is the software translated into Julia.

The Julia rewrite is released under the MIT license with an express permission from the authors of the Fortran package. The package uses multiple dispatch to route around standard BLAS routines in the case e.g. of arbitrary-precision floating point numbers or ForwardDiff.Dual. This e.g. allows for Automatic Differentiation (AD) of a sparse-matrix solve.

CliqueTreesFactorization(
    alg = nothing,
    snd = nothing,
    reuse_symbolic = true,
)

The sparse Cholesky factorization algorithm implemented in CliqueTrees.jl. The implementation is pure-Julia and accepts arbitrary numeric types. It is somewhat slower than CHOLMOD.

KrylovJL_CG(args...; kwargs...)

A generic CG implementation for Hermitian and positive definite linear systems

KrylovJL_MINRES(args...; kwargs...)

A generic MINRES implementation for Hermitian linear systems

KrylovJL_GMRES(args...; gmres_restart = 0, window = 0, kwargs...)

A generic GMRES implementation for square non-Hermitian linear systems

KrylovJL_BICGSTAB(args...; kwargs...)

A generic BICGSTAB implementation for square non-Hermitian linear systems

KrylovJL_LSMR(args...; kwargs...)

A generic LSMR implementation for least-squares problems

KrylovJL_CRAIGMR(args...; kwargs...)

A generic CRAIGMR implementation for least-norm problems

KrylovJL(args...; KrylovAlg = Krylov.gmres!,
    Pl = nothing, Pr = nothing,
    gmres_restart = 0, window = 0,
    kwargs...)

A generic wrapper over the Krylov.jl krylov-subspace iterative solvers.

MKLLUFactorization()

A wrapper over Intel's Math Kernel Library (MKL). Direct calls to MKL in a way that pre-allocates workspace to avoid allocations and does not require libblastrampoline.

MKL32MixedLUFactorization()

A mixed precision LU factorization using Intel MKL that performs factorization in Float32 precision while maintaining Float64 interface. This can provide significant speedups for large matrices when reduced precision is acceptable.

Performance Notes

• Converts Float64 matrices to Float32 for factorization
• Uses optimized MKL routines for the factorization
• Can be 2x faster than full precision for memory-bandwidth limited problems
• May have reduced accuracy compared to full Float64 precision

Requirements

This solver requires MKL to be available through MKL_jll.

Example

alg = MKL32MixedLUFactorization()
sol = solve(prob, alg)

OpenBLASLUFactorization()

A direct wrapper over OpenBLAS's LU factorization (getrf! and getrs!). This solver makes direct calls to OpenBLAS_jll without going through Julia's libblastrampoline, which can provide performance benefits in certain configurations.

Performance Characteristics

• Pre-allocates workspace to avoid allocations during solving
• Makes direct ccalls to OpenBLAS routines
• Can be faster than LUFactorization when OpenBLAS is well-optimized for the hardware
• Supports Float32, Float64, ComplexF32, and ComplexF64 element types

When to Use

• When you want to ensure OpenBLAS is used regardless of the system BLAS configuration
• When benchmarking shows better performance than LUFactorization on your specific hardware
• When you need consistent behavior across different systems (always uses OpenBLAS)

Example

using LinearSolve, LinearAlgebra

A = rand(100, 100)
b = rand(100)
prob = LinearProblem(A, b)
sol = solve(prob, OpenBLASLUFactorization())

AppleAccelerateLUFactorization()

A wrapper over Apple's Accelerate Library. Direct calls to Acceelrate in a way that pre-allocates workspace to avoid allocations and does not require libblastrampoline.

AppleAccelerate32MixedLUFactorization()

A mixed precision LU factorization using Apple's Accelerate framework that performs factorization in Float32 precision while maintaining Float64 interface. This can provide significant speedups on Apple hardware when reduced precision is acceptable.

Performance Notes

• Converts Float64 matrices to Float32 for factorization
• Uses optimized Accelerate routines for the factorization
• Particularly effective on Apple Silicon with unified memory
• May have reduced accuracy compared to full Float64 precision

Requirements

This solver is only available on Apple platforms and requires the Accelerate framework.

Example

alg = AppleAccelerate32MixedLUFactorization()
sol = solve(prob, alg)

MetalLUFactorization()

A wrapper over Apple's Metal GPU library for LU factorization. Direct calls to Metal in a way that pre-allocates workspace to avoid allocations and automatically offloads to the GPU. This solver is optimized for Metal-capable Apple Silicon Macs.

Requirements

Using this solver requires that Metal.jl is loaded: using Metal

Performance Notes

• Most efficient for large dense matrices where GPU acceleration benefits outweigh transfer costs
• Automatically manages GPU memory and transfers
• Particularly effective on Apple Silicon Macs with unified memory

Example

using Metal
alg = MetalLUFactorization()
sol = solve(prob, alg)

MetalOffload32MixedLUFactorization()

A mixed precision Metal GPU-accelerated LU factorization that converts matrices to Float32 before offloading to Metal GPU for factorization, then converts back for the solve. This can provide speedups on Apple Silicon when reduced precision is acceptable.

Performance Notes

• Converts Float64 matrices to Float32 for GPU factorization
• Can be significantly faster for large matrices where memory bandwidth is limiting
• Particularly effective on Apple Silicon Macs with unified memory architecture
• May have reduced accuracy compared to full precision methods

Requirements

Using this solver requires that Metal.jl is loaded: using Metal

Example

using Metal
alg = MetalOffload32MixedLUFactorization()
sol = solve(prob, alg)

MKLPardisoFactorize(; nprocs::Union{Int, Nothing} = nothing,
    matrix_type = nothing,
    cache_analysis = false,
    iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing,
    dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing)

A sparse factorization method using MKL Pardiso.

Keyword Arguments

Setting cache_analysis = true disables Pardiso's scaling and matching defaults and caches the result of the initial analysis phase for all further computations with this solver.

For the definition of the other keyword arguments, see the Pardiso.jl documentation. All values default to nothing and the solver internally determines the values given the input types, and these keyword arguments are only for overriding the default handling process. This should not be required by most users.

MKLPardisoIterate(; nprocs::Union{Int, Nothing} = nothing,
    matrix_type = nothing,
    cache_analysis = false,
    iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing,
    dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing)

A mixed factorization+iterative method using MKL Pardiso.

Keyword Arguments

Setting cache_analysis = true disables Pardiso's scaling and matching defaults and caches the result of the initial analysis phase for all further computations with this solver.

For the definition of the other keyword arguments, see the Pardiso.jl documentation. All values default to nothing and the solver internally determines the values given the input types, and these keyword arguments are only for overriding the default handling process. This should not be required by most users.

PardisoJL(; nprocs::Union{Int, Nothing} = nothing,
    solver_type = nothing,
    matrix_type = nothing,
    iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing,
    dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing,
    vendor::Union{Symbol, Nothing} = nothing
)

A generic method using Pardiso. Specifying solver_type is required.

Keyword Arguments

The vendor keyword allows to choose between Panua pardiso (former pardiso-project.org; vendor=:Panua) and MKL Pardiso (vendor=:MKL). If vendor==nothing, Panua pardiso is preferred over MKL Pardiso.

For the definition of the other keyword arguments, see the Pardiso.jl documentation. All values default to nothing and the solver internally determines the values given the input types, and these keyword arguments are only for overriding the default handling process. This should not be required by most users.

CudaOffloadLUFactorization()

An offloading technique used to GPU-accelerate CPU-based computations using LU factorization. Requires a sufficiently large A to overcome the data transfer costs.

CudaOffloadQRFactorization()

An offloading technique used to GPU-accelerate CPU-based computations using QR factorization. Requires a sufficiently large A to overcome the data transfer costs.

CUDAOffload32MixedLUFactorization()

A mixed precision GPU-accelerated LU factorization that converts matrices to Float32 before offloading to CUDA GPU for factorization, then converts back for the solve. This can provide speedups when the reduced precision is acceptable and memory bandwidth is a bottleneck.

Performance Notes

• Converts Float64 matrices to Float32 for GPU factorization
• Can be significantly faster for large matrices where memory bandwidth is limiting
• May have reduced accuracy compared to full precision methods
• Most beneficial when the condition number of the matrix is moderate

AMDGPUOffloadLUFactorization()

An offloading technique using LU factorization to GPU-accelerate CPU-based computations on AMD GPUs. Requires a sufficiently large A to overcome the data transfer costs.

AMDGPUOffloadQRFactorization()

An offloading technique using QR factorization to GPU-accelerate CPU-based computations on AMD GPUs. Requires a sufficiently large A to overcome the data transfer costs.

CUSOLVERRFFactorization(; symbolic = :RF, reuse_symbolic = true)

A GPU-accelerated sparse LU factorization using NVIDIA's cusolverRF library. This solver is specifically designed for sparse matrices on CUDA GPUs and provides high-performance factorization and solve capabilities.

Keyword Arguments

• symbolic: The symbolic factorization method to use. Options are:
  • :RF (default): Use cusolverRF's built-in symbolic analysis
  • :KLU: Use KLU for symbolic analysis
• reuse_symbolic: Whether to reuse the symbolic factorization when the sparsity pattern doesn't change (default: true)

IterativeSolversJL_CG(args...; Pl = nothing, Pr = nothing, kwargs...)

A wrapper over the IterativeSolvers.jl CG.

IterativeSolversJL_GMRES(args...; Pl = nothing, Pr = nothing, gmres_restart = 0, kwargs...)

A wrapper over the IterativeSolvers.jl GMRES.

IterativeSolversJL_BICGSTAB(args...; Pl = nothing, Pr = nothing, kwargs...)

A wrapper over the IterativeSolvers.jl BICGSTAB.

IterativeSolversJL_MINRES(args...; Pl = nothing, Pr = nothing, kwargs...)

A wrapper over the IterativeSolvers.jl MINRES.

IterativeSolversJL_IDRS(args...; Pl = nothing, kwargs...)

A wrapper over the IterativeSolvers.jl IDR(S).

IterativeSolversJL(args...;
    generate_iterator = IterativeSolvers.gmres_iterable!,
    Pl = nothing, Pr = nothing,
    gmres_restart = 0, kwargs...)

A generic wrapper over the IterativeSolvers.jl solvers.

KrylovKitJL_CG(args...; Pl = nothing, Pr = nothing, kwargs...)

A generic CG implementation for Hermitian and positive definite linear systems

KrylovKitJL_GMRES(args...; Pl = nothing, Pr = nothing, gmres_restart = 0, kwargs...)

A generic GMRES implementation.

KrylovKitJL(args...; KrylovAlg = Krylov.gmres!, kwargs...)

A generic iterative solver implementation allowing the choice of KrylovKit.jl solvers.

HYPREAlgorithm(solver; Pl = nothing)

HYPRE.jl is an interface to hypre and provide iterative solvers and preconditioners for sparse linear systems. It is mainly developed for large multi-process distributed problems (using MPI), but can also be used for single-process problems with Julias standard sparse matrices.

If you need more fine-grained control over the solver/preconditioner options you can alternatively pass an already created solver to HYPREAlgorithm (and to the Pl keyword argument). See HYPRE.jl docs for how to set up solvers with specific options.

Positional Arguments

The single positional argument solver has the following choices:

• HYPRE.BiCGSTAB
• HYPRE.BoomerAMG
• HYPRE.FlexGMRES
• HYPRE.GMRES
• HYPRE.Hybrid
• HYPRE.ILU
• HYPRE.ParaSails (as preconditioner only)
• HYPRE.PCG

Keyword Arguments

• Pl: A choice of left preconditioner.

Example

For example, to use HYPRE.PCG as the solver, with HYPRE.BoomerAMG as the preconditioner, the algorithm should be defined as follows:

A, b = setup_system(...)
prob = LinearProblem(A, b)
alg = HYPREAlgorithm(HYPRE.PCG)
prec = HYPRE.BoomerAMG
sol = solve(prob, alg; Pl = prec)