Repository navigation
Bordered solvers [BifurcationKit jacobian-free methods optimization] #1
Description
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?
Activity
ChrisRackauckas commentedon Dec 1, 2024
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 commentedon Dec 1, 2024
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 commentedon Dec 3, 2024
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.AJ0070 commentedon Sep 1, 2026
I put together a package for this: BorderedLinearSolve.jl. It does block elimination on top of LinearSolve, so both inner solves share one cache and it carries between continuation steps, which is the Krylov reuse this issue was about.
@rveltz it runs under
CUDA.allowscalar(false)with a CuArray or CuSparseMatrixCSR J, at 7e-16 with LU and 7e-10 with GMRES on a T4 at n=200. @j-fu wide borders take the same path, m inner solves and then a small m by m system, and BEC+k is in there: on a J with a 1e-11 singular value, k=0 gives 6.8e-5 and k=3 gives 1.9e-15.Not registered yet. The interface has only had my eyes on it so far, so feedback on that would be welcome before it gets fixed in a registry.
rveltz commentedon Sep 2, 2026
For now, the only thing I see is that line https://github.com/AJ0070/BorderedLinearSolve.jl/blob/main/src/problem.jl#L29 where you seem to imply an abstract Matrix.
I dont think you need this. Also, I would not require to have
AbstractVectorfor the rhs. When you solve 2d-pde, you have a matrix to reprensent the solution. Forcing it to be anAbstractVector, requires the user to usevec,reshape. I know thatKrylovmakes that requirements but I dont think it need to.KrylovKitis was more user friendly in this case.AJ0070 commentedon Sep 2, 2026
Both fixed, thanks. The border width now comes from an explicit tuple of states instead of being inferred from an array's shape, so a matrix-shaped state is no longer mistaken for a wide border. The
size(J, ...)checks are gone as well, sinceJis only ever applied and solved against and may be matrix-free.The state can be any array now. A 4x4 state with a matrix-free operator solves to 8e-16 and comes back 4x4 rather than flattened.
You are right about the underlying limitation though. LinearSolve reads a matrix right-hand side as multiple right-hand sides, so what I do is flatten at the inner solve boundary and restore the shape on the way out. The bordering algebra itself never sees a flattened state, but that is a workaround on my side rather than LinearSolve handling array-valued states the way KrylovKit does.
AJ0070 commentedon Sep 4, 2026
@rveltz
rveltz commentedon Sep 4, 2026
That looks pretty good to me.