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.

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

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.

CudaOffloadFactorization()

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

DiagonalFactorization()

A special implementation only for solving Diagonal matrices fast.

FastLUFactorization()

The FastLapackInterface.jl version of the LU factorization. Notably, this version does not allow for choice of pivoting method.

FastQRFactorization()

The FastLapackInterface.jl version of the QR factorization.

GenericFactorization(;fact_alg=LinearAlgebra.factorize): Constructs a linear solver from a generic factorization algorithm fact_alg which complies with the Base.LinearAlgebra factorization API. Quoting from Base:

  * If `A` is upper or lower triangular (or diagonal), no factorization of `A` is
    required. The system is then solved with either forward or backward substitution.
    For non-triangular square matrices, an LU factorization is used.
    For rectangular `A` the result is the minimum-norm least squares solution computed by a
    pivoted QR factorization of `A` and a rank estimate of `A` based on the R factor.
    When `A` is sparse, a similar polyalgorithm is used. For indefinite matrices, the `LDLt`
    factorization does not use pivoting during the numerical factorization and therefore the
    procedure can fail even for invertible matrices.

Keyword Arguments

• fact_alg: the factorization algorithm to use. Defaults to LinearAlgebra.factorize, but can be swapped to choices like lu, qr

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

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)

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

A generic wrapper over the IterativeSolvers.jl solvers.

KLUFactorization(;reuse_symbolic=true, check_pattern=true)

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

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.

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

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

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

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.

MetalLUFactorization()

A wrapper over Apple's Metal GPU library. Direct calls to Metal in a way that pre-allocates workspace to avoid allocations and automatically offloads to the GPU.

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.

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

OperatorAssumptions(issquare = nothing; condition::OperatorCondition.T = IllConditioned)

Sets the operator A assumptions used as part of the default algorithm

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)

A generic method using MKL Pardiso. Specifying solver_type is required.

Keyword Arguments

For the definition of the 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.

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.

RFLUFactorization()

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.

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.

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.

SimpleLUFactorization(pivot::Bool = true)

A simple LU-factorization implementation without BLAS. Fast for small matrices.

Positional Arguments

• pivot::Bool: whether to perform pivoting. Defaults to true

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.

UMFPACKFactorization(;reuse_symbolic=true, check_pattern=true)

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

if alg.alg === DefaultAlgorithmChoice.LUFactorization SciMLBase.solve!(cache, LUFactorization(), args...; kwargs...)) else ... end

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

A wrapper over the IterativeSolvers.jl BICGSTAB.

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_IDRS(args...; Pl = nothing, kwargs...)

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

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

A wrapper over the IterativeSolvers.jl MINRES.

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

A generic BICGSTAB implementation for square non-Hermitian linear systems

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

A generic CG implementation for Hermitian and positive definite linear systems

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

A generic CRAIGMR implementation for least-norm problems

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

A generic GMRES implementation for square non-Hermitian linear systems

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

A generic LSMR implementation for least-squares problems

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

A generic MINRES implementation for Hermitian linear systems

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.

MKLPardisoFactorize(; nprocs::Union{Int, Nothing} = nothing,
                    matrix_type = nothing,
                    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

For the definition of the 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,
                    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

For the definition of the 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.

(c, s, ρ) = _sym_givens(a, b)

Numerically stable symmetric Givens reflection. Given a and b reals, return (c, s, ρ) such that

[ c  s ] [ a ] = [ ρ ]
[ s -c ] [ b ] = [ 0 ].

elseif DefaultAlgorithmChoice.LUFactorization === cache.alg
    (cache.cacheval.LUFactorization)' \ dy
else
    ...
end

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

needs_square_A(alg)

Returns true if the algorithm requires a square matrix.

set_A(cache, A)

set_b(cache, b)

set_p(cache, p)

set_prec(cache, Pl, Pr)

set_u(cache, u)

OperatorCondition

Specifies the assumption of matrix conditioning for the default linear solver choices. Condition number is defined as the ratio of eigenvalues. The numerical stability of many linear solver algorithms can be dependent on the condition number of the matrix. The condition number can be computed as:

using LinearAlgebra
cond(rand(100,100))

However, in practice this computation is very expensive and thus not possible for most practical cases. Therefore, OperatorCondition lets one share to LinearSolve the expected conditioning. The higher the expected condition number, the safer the algorithm needs to be and thus there is a trade-off between numerical performance and stability. By default the method assumes the operator may be ill-conditioned for the standard linear solvers to converge (such as LU-factorization), though more extreme ill-conditioning or well-conditioning could be the case and specified through this assumption.

OperatorCondition.IllConditioned

The default assumption of LinearSolve. Assumes that the operator can have minor ill-conditioning and thus needs to use safe algorithms.

OperatorCondition.SuperIllConditioned

Assumes that the operator can have fairly extreme ill-conditioning and thus the most stable algorithm is used.

OperatorCondition.VeryIllConditioned

Assumes that the operator can have fairly major ill-conditioning and thus the standard linear algebra algorithms cannot be used.

OperatorCondition.WellConditioned

Assumes that the operator can have fairly contained conditioning and thus the fastest algorithm is used.