make preconditioners part of the solver rather than a random extra

[oscardssmith]

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 a precs function that takes in A and possibly some parameters and remakes your preconditioner for the new A value.

Here's an example of the new functionality in action:

using LinearSolve, Krylov, KrylovPreconditioners, LinearAlgebra, SparseArrays
A, b = sprand(10,10, .5), rand(10);
prob = LinearProblem(A,b)
solver = KrylovJL_CG(precs=(A,p=nothing) -> (BlockJacobiPreconditioner(A, 2), I), ldiv=false)
cache = init(prob, solver)
solve!(cache)
reinit!(cache, A = sprand(10,10, .2), b = rand(10))
solve!(cache)

[j-fu]

... This is the one functionality I missed with LinearSolve. I had a slightly different view on this: see #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.

[oscardssmith]

@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.

[j-fu]

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

[oscardssmith]

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)

[j-fu]

Let's have a chat on this at JuliaCon!

[j-fu]

Was nice to talk to you and to hear that there is quite some interest in this PR in the SciML org. How can I help do move this on ?

[oscardssmith]

This is mostly one that I think I just need to put a couple more hours into to get the Pr/Pl/precs handling uniform across everything. This is my next priority after SciML/OrdinaryDiffEq.jl#229.

[oscardssmith]

This is now working and tested. I'm planning on merging soon assuming CI is clear.

[oscardssmith]

This is failing a lot of tests, but so is master. is this ready to merge?

[ChrisRackauckas]

@fredrikekre can you figure out what LoadError: mpiexec() is deprecated, use the non-do-block form is?

[fredrikekre]

#501 (comment) + #501 (comment) :)