Unified handling of preconditoners and factorizations ?

[j-fu]

Hi,

I see that 2.0 is underway, so this may be an occasion to have some discussion about the API. The current API handles preconditioners and direct solvers very differently. IMHO they both are (possibly incomplete) factorizations and could be handeled in a similar way.

OTOH, the current API puts UMFPACKFactorization() and the like and iterative solvers into one basket which seems to be not very logical to me either.

So my view is the following:

• A complete factorization provides a solver: Pardiso, Umfpack, KLU, Cholesky etc.

• whereas an approximate factorization provides a preconditioner: ILU0, Jacobi, AlgebraicMultigrid etc.

Iterative algorithms are CG, GMRES etc. Direct solution essentially corresponds to a "OneStep" algorithm which works together with a complete factorization aka solver.

In that sense, a solve API could look as follows:

solve(P::LinearProblem, F::Factorization, I::Iteration)

and one could have something like

solve(P::LinearProblem, F::Factorization)= solve(P,F,OneStep())

(with a possible warning that an incomplete factorization doesn't solve)

As factorization one would provide a preconditioner constructor instead of a the preconditioner itself. So we could have

solve(p, KLUFactorization())
solve(p, ILU0Factorization(), CG())

Repeated solve after updating the matrix could update the preconditioner in a similar way as a factorization.

May be this proposal comes too late or doesn't fit for some other reasons, I can live either way (and already coded around this...).

[ChrisRackauckas]

This sounds interesting, though v2 is already late and hopefully going out today, so unless you plan on making a PR soon it sounds like a 3.0. I'm not sure it's breaking either since it looks like an extension: we do not handle preconditioners in the interface right now other than passing them along to some solvers, so elevating that to some new preconditioner interface would be possible without even breaking any of the "onestep" method interfaces.

Couldn't we just have solve(p, ILU0Factorization()) as something that solves with ILU0 and then you can pass ILU0Factorization() as a preconditioner choice and it'll know to do that?

[j-fu]

Unfortunately I have no uninterrupted time to make a reasonable PR immediately (which also would involve some thinking).

But I see your interest and would work on a PR in the next couple of weeks.

I coded something like this for VoronoiFVM and could try to lift these ideas over here.

Couldn't we just have solve(p, ILU0Factorization()) as something that solves with ILU0 and then you can pass ILU0Factorization() as a preconditioner choice and it'll know to do that?

Basically this would be the idea. But it would have to be solve(p, ILU0Factorization(), KrylovJL_CG()) unless there is an automatic way to choose the iterative solver algorithm (among cg, bicgstab, gmres at least...).

[oscardssmith]

Closing in favor of https://github.com/SciML/LinearSolve.jl/pull/514