|
| 1 | +# Developing New Linear Solvers |
| 2 | + |
| 3 | +Developing new or custom linear solvers for the SciML interface can be done in |
| 4 | +one of two ways: |
| 5 | + |
| 6 | +1. You can either create a completely new set of dispatches for `init` and `solve`. |
| 7 | +2. You can extend LinearSolve.jl's internal mechanisms. |
| 8 | + |
| 9 | +For developer ease, we highly recommend (2) as that will automatically make the |
| 10 | +caching API work. Thus this is the documentation for how to do that. |
| 11 | + |
| 12 | +## Developing New Linear Solvers with LinearSolve.jl Primitives |
| 13 | + |
| 14 | +Let's create a new wrapper for a simple LU-factorization which uses only the |
| 15 | +basic machinery. A simplified version is: |
| 16 | + |
| 17 | +```julia |
| 18 | +struct MyLUFactorization{P} <: SciMLBase.AbstractLinearAlgorithm end |
| 19 | + |
| 20 | +init_cacheval(alg::MyLUFactorization, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose) = lu!(convert(AbstractMatrix,A)) |
| 21 | + |
| 22 | +function SciMLBase.solve(cache::LinearCache, alg::MyLUFactorization; kwargs...) |
| 23 | + if cache.isfresh |
| 24 | + A = convert(AbstractMatrix,A) |
| 25 | + fact = lu!(A) |
| 26 | + cache = set_cacheval(cache, fact) |
| 27 | + end |
| 28 | + y = ldiv!(cache.u, cache.cacheval, cache.b) |
| 29 | + SciMLBase.build_linear_solution(alg,y,nothing,cache) |
| 30 | +end |
| 31 | +``` |
| 32 | + |
| 33 | +The way this works is as follows. LinearSolve.jl has a `LinearCache` that everything |
| 34 | +shares (this is what gives most of the ease of use). However, many algorithms |
| 35 | +need to cache their own things, and so there's one value `cacheval` that is |
| 36 | +for the algorithms to modify. The function: |
| 37 | + |
| 38 | +```julia |
| 39 | +init_cacheval(alg::MyLUFactorization, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose) |
| 40 | +``` |
| 41 | + |
| 42 | +is what is called at `init` time to create the first `cacheval`. Note that this |
| 43 | +should match the type of the cache later used in `solve` as many algorithms, like |
| 44 | +those in OrdinaryDiffEq.jl, expect type-groundedness in the linear solver definitions. |
| 45 | +While there are cheaper ways to obtain this type for LU factorizations (specifically, |
| 46 | +`ArrayInterface.lu_instance(A)`), for a demonstration this just performs an |
| 47 | +LU-factorization to get an `LU{T, Matrix{T}}` which it puts into the `cacheval` |
| 48 | +so its typed for future use. |
| 49 | + |
| 50 | +After the `init_cacheval`, the only thing left to do is to define |
| 51 | +`SciMLBase.solve(cache::LinearCache, alg::MyLUFactorization)`. Many algorithms |
| 52 | +may use a lazy matrix-free representation of the operator `A`. Thus if the |
| 53 | +algorithm requires a concrete matrix, like LU-factorization does, the algorithm |
| 54 | +should `convert(AbstractMatrix,cache.A)`. The flag `cache.isfresh` states whether |
| 55 | +`A` has changed since the last `solve`. Since we only need to factorize when |
| 56 | +`A` is new, the factorization part of the algorithm is done in a `if cache.isfresh`. |
| 57 | +`cache = set_cacheval(cache, fact)` puts the new factorization into the cache |
| 58 | +so it's updated for future solves. Then `y = ldiv!(cache.u, cache.cacheval, cache.b)` |
| 59 | +performs the solve and a linear solution is returned via |
| 60 | +`SciMLBase.build_linear_solution(alg,y,nothing,cache)`. |
0 commit comments