Julia users are building a wide variety of applications in the SciML ecosystem,
often requiring problem-specific handling of their linear solves. As existing solvers in LinearSolve.jl
may not
be optimally suited for novel applications, it is essential for the linear solve
interface to be easily extendable by users. To that end, the linear solve algorithm
LinearSolveFunction()
accepts a user-defined function for handling the solve. A
user can pass in their custom linear solve function, say my_linsolve
, to
LinearSolveFunction()
. A contrived example of solving a linear system with a custom solver is below.
using LinearSolve, LinearAlgebra
function my_linsolve(A, b, u, p, newA, Pl, Pr, solverdata; verbose = true, kwargs...)
if verbose == true
println("solving Ax=b")
end
u = A \ b
return u
end
prob = LinearProblem(Diagonal(rand(4)), rand(4))
alg = LinearSolveFunction(my_linsolve)
sol = solve(prob, alg)
sol.u
The inputs to the function are as follows:
u
, the solution initialized aszero(b)
,A
, the linear operatorb
, the right-hand-sidep
, a set of parametersnewA
, aBool
which istrue
ifA
has been modified since last solvePl
, left-preconditionerPr
, right-preconditionersolverdata
, solver cache set tonothing
if solver hasn't been initializedkwargs
, standard SciML keyword arguments such asverbose
,maxiters
,abstol
,reltol
The function my_linsolve
must accept the above specified arguments and modify the, and return
the solution, u
. As memory for u
is already allocated, the user may choose
to modify u
in place as follows:
function my_linsolve!(A, b, u, p, newA, Pl, Pr, solverdata; verbose = true, kwargs...)
if verbose == true
println("solving Ax=b")
end
u .= A \ b # in place
return u
end
alg = LinearSolveFunction(my_linsolve!)
sol = solve(prob, alg)
sol.u