Open oscardssmith opened 2 weeks ago
... This is the one functionality I missed with LinearSolve. I had a slightly different view on this: see https://github.com/SciML/LinearSolve.jl/issues/308 . Unfortunately I had no time to work on this, so thank you for this PR!
As for the functionality: Sometimes, e.g. in Newton solvers, it makes sense to update the matrix and keep the preconditioner. The "old" preconditioner may still be good enough for the "new" matrix in the sense that few additional iterations are cheaper than the re-creation of the preconditioner (thinking e.g. of Algebraic multigrid etc.)
May be a kwarg keep_precs
(with default value false) in reinit!
can help to achieve this.
@j-fu yeah, the only reason I thought about this API is that OrdinaryDiffEq has a (slightly broken very hacky) version of this API already (dolinsolve
) which someone pointed out to me a few weeks ago that it's slightly broken.
A couple of thoughts for this API extension:
What about keeping the Pl/Pr idiom (as functions) instead of precs
handling tuples ? IMHO this would cover the vast majority of use cases.
LinearSolve (and the whole SciML) has this pattern that we write e.g.UMFPACKFactorization()
to specify a factorization. In a similar way, users could expect to specify
precs=ILUZero_ILU0()
instead of
precs= (A,p=nothing)-> (ilu0(A),I)
May one could achieve both ways via a dispatch between Function
and something like AbstractPreconditioningAlgorithm
. Also this would give the possibility to have a couple of "curated" preconditioners predefined (e.g. in package extensions).
UPDATE: Well, this already works with this PR :)
struct ILUZero_ILU0 end
(::ILUZero_ILU0)(A, p=nothing)=(ilu0(A),I)
p=LinearProblem(A,b)
solve(p,KrylovJL_GMRES(precs=ILUZero_ILU0(),ldiv=true) )
Yeah, we likely would want to add aliases to create iterative algorithms with builtin preconditioners. My reasoning behind precs
rather than Pr
/Pl
is that if you want both preconditioners, it's likely that they will share computational work. We could make it so you are allowed to pass up to one of Pr
, Pl
, or precs
(where passing Pr
or Pl
would set the other one to I
)
Let's have a chat on this at JuliaCon!
This still needs tests, but I'm putting it up as a PR first to get feedback. The intended goal here is to make the
dolinsolve
functionality in OrdinaryDiffEq of being able to pass in aprecs
function that takes inA
and possibly some parameters and remakes your preconditioner for the newA
value.Here's an example of the new functionality in action: