Bordered solvers [BifurcationKit jacobian-free methods optimization]

[rveltz]

Hi,

In BifurcationKit, we repeatedly solve bordered linear solvers (see https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/borderedlinearsolver/), also with Jacobian-free solvers. At the moment, BifurcationKit does not take advantage of LinearSolve.

I could combine the bordered linear solvers of BK with LinearSolve.jl for example to take advantage of the Krylov space re-use between continuation steps.

Should it be another package like BorderedLinearSolve.jl, should I PR to LinearSolve.jl or keep the infrastructure in BifurcationKit.jl?

[ChrisRackauckas]

I think this would need a separate package because it would be a substantial amount of code for the new solvers? I would be fine maintaining it in SciML if someone put it together but I haven't had a particular need for it so it just never came about. It does seem a bit specialized though.

[j-fu]

I have something like this on my task list, outside of the bifurcation context, coming from a semiconductor laser model. I would be interested in wide bordered systems.

I did some experiments with BlockArrays.jl to handle mixed sparse/dense situations which played out quite well. Used GMRES directly from Krylov.jl for that.

There could also be an implementation of
Govaerts/Pryce: Mixed block elimination for linear systems with wider borders, IMA Journal of Numerical Analysis (1993) 13, 161-180

So yes, such a package would be an interesting effort. I might participate in that.

[rveltz]

For BifurcationKit, this must run on GPU. Hence, in order to avoid using scalar indexing to solve a system bordered by a scalar, I introduced a type BorderedArray. Otherwise, all this is relatively straightforward.