Fast implementations of linear solving algorithms in Julia that satisfy the SciML common interface. LinearSolve.jl makes it easy to define high level algorithms which allow for swapping out the linear solver that is used while maintaining maximum efficiency. Specifically, LinearSolve.jl includes:

• Fast pure Julia LU factorizations which outperform standard BLAS
• KLU for faster sparse LU factorization on unstructured matrices
• UMFPACK for faster sparse LU factorization on matrices with some repeated structure
• MKLPardiso wrappers for handling many sparse matrices faster than SuiteSparse (KLU, UMFPACK) methods
• Sparspak.jl for sparse LU factorization in pure Julia for generic number types and for non-GPL distributions
• GPU-offloading for large dense matrices
• Wrappers to all of the Krylov implementations (Krylov.jl, IterativeSolvers.jl, KrylovKit.jl) for easy testing of all of them. LinearSolve.jl handles the API differences, especially with the preconditioner definitions
• A polyalgorithm that smartly chooses between these methods
• A caching interface which automates caching of symbolic factorizations and numerical factorizations as optimally as possible

For information on using the package, see the stable documentation. Use the in-development documentation for the version of the documentation which contains the unreleased features.

n = 4
A = rand(n, n)
b1 = rand(n);
b2 = rand(n);
prob = LinearProblem(A, b1)

linsolve = init(prob)
sol1 = solve(linsolve)

sol1.u
#=
4-element Vector{Float64}:
 -0.9247817429364165
 -0.0972021708185121
  0.6839050402960025
  1.8385599677530706
=#

linsolve = LinearSolve.set_b(linsolve, b2)
sol2 = solve(linsolve)

sol2.u
#=
4-element Vector{Float64}:
  1.0321556637762768
  0.49724400693338083
 -1.1696540870182406
 -0.4998342686003478
=#

linsolve = LinearSolve.set_b(linsolve, b2)
sol2 = solve(linsolve, IterativeSolversJL_GMRES()) # Switch to GMRES
sol2.u
#=
4-element Vector{Float64}:
  1.0321556637762768
  0.49724400693338083
 -1.1696540870182406
 -0.4998342686003478
=#

A2 = rand(n, n)
linsolve = LinearSolve.set_A(linsolve, A2)
sol3 = solve(linsolve)

sol3.u
#=
4-element Vector{Float64}:
 -6.793605395935224
  2.8673042300837466
  1.1665136934977371
 -0.4097250749016653
=#