Replace IterativeSolvers.jl by Krylov.jl

[amontoison]

No description provided.

[amontoison]

The error with the test set CPU hydrostatic free surface model is probably related to a preconditioner that is not symmetric positive definite.

[simone-silvestri]

Hello, @amontoison. Nice work. Do you see a performance improvement when switching to this package?
There should be some benchmarks in the benchmark folder that we can test (probably we need to update that folder a bit, let me know if you have problems with it)

[amontoison]

Hello, @amontoison. Nice work. Do you see a performance improvement when switching to this package? There should be some benchmarks in the benchmark folder that we can test (probably we need to update that folder a bit, let me know if you have problems with it)

Hi @simone-silvestri, cg is not too hard to implement, so I don't think we will see a significant performance improvement on CPU.
However, for other methods like gmres, Krylov.jl easily outperforms IterativeSolvers.jl.
On GPU, though, we might see a difference because I try to dispatch to the BLAS/LAPACK routines of the GPU vendors as much as possible, whereas IterativeSolvers.jl relies on some broadcast.

Also, to the best of my knowledge, only cg works on (NVIDIA) GPUs for IterativeSolvers.jl, while all solvers in Krylov.jl work on the GPUs of any vendor.

I'll try to run some benchmarks before the end of the week.

[glwagner]

@simone-silvestri can we remove the matrix solver? I don't think we use it anymore and it will help us move forward.

[amontoison]

We can remove the dependency on IterativeSolvers.jl if we remove HeptadiagonalIterativeSolver.
For Krylov.jl, I developed a dedicated solver in #4041.
I don't need to replace IterativeSolvers.jl by Krylov.jl if we remove the HeptadiagonalIterativeSolver.

[simone-silvestri]

@simone-silvestri can we remove the matrix solver? I don't think we use it anymore and it will help us move forward.

I can remove it. I ll open a PR. I can merge it next week, however, it is still our default implicit solver for non-rectilinear grids, so I suggest a bit of benchmarking before merging the PR

[amontoison]

@simone-silvestri I close the PR.
Maybe you can keep the solver but switch to your own implementation of CG or the one of Krylov.jl.
I wanted to reduce the number of dependencies.

[simone-silvestri]

Well, we can also try switching the implementation to Krylov over here and seeing the results.
I am not 100% sure we should remove the solver if it is faster than our pcg implementation.

[amontoison]

I benchmarked a few times Krylov.jl versus the other Julia packages for Krylov methods and it was never slower.
One big mistake that you do in your implementation is that you compute norm(r) explicitly while you have this information for free.

[amontoison]

-> https://github.com/CliMA/Oceananigans.jl/blob/main/src/Solvers/conjugate_gradient_solver.jl#L194
If you take the square root of this quantity you have your residual norm and don't need to compute anything at each iteration for your stopping conditions.

[simone-silvestri]

that's a good idea! Thanks for the suggestion, let's try it out. I think there are still performance problems with our custom pcg solver, that is why I am a bit cautious in removing completely the other implementation. If we do so before fixing our pcg, we will not be able to compare our custom implementation anymore.

[glwagner]

@simone-silvestri I think it is a bad infrastructure strategy to keep around lots of methods. Performance at the present is not that important, what matters far more is the direction we want to go with the code.