use of AbstractArray

[rveltz]

Hi,

In the matrix free context, I would like to apply an operator to a Matrix representing the solution of a PDE. In order to use Krylov.jl, I have to pass AbstractVectors hence reshape / views. Is there a better workflow?

Thank you a lot

[amontoison]

Hi Romain!
Is this solution not good enough for you:
https://jso.dev/Krylov.jl/dev/matrix_free/#Example-with-discretized-PDE ?

I have another alternative if you want but I would like to know the type of the matrix before (just Matrix or something on GPU / custom type)?

[rveltz]

You are using vec basically hence a reshape. It means I have to dispatch my operator on AbstractVectors

like

function myop(x::AbstractVector)
    xc = reshape(x, N, N)
    return vec(myop(xc))
end

I was hoping for a solution ala KrylovKit where you dont need those but I guess this is not a big deal. IterativeSolvers has the same requirements

[amontoison]

You allocate 32 bits at each vec / reshape.
So it's quite cheap but you could remove these allocations with a custom workspace: https://jso.dev/Krylov.jl/dev/custom_workspaces/

You just need to create your own AbstractVector that wraps your matrix.
The drawback is that you will lose by default the dispatch to BLAS / CUBLAS / rocBLAS for the linear algebra operations.

I only recommend it for an advanced use like multi-GPU with MPI.

[rveltz]

Ok it just halfed my computation time and I did not optimise it.

With bifurcationKit, I regret now not having tried earlier, the fact that you can reuse the krylov space is a massive boost in perfs, I thought KrylovKit was good enough

[amontoison]

If you use it on a GPU, the speed-up will be even greater because allocations and deallocations are more expensive, and the garbage collection system works differently.

KrylovKit.jl is probably more flexible but was not developed for performances and also doesn't support preconditioners which is quite relevant for Krylov methods.

[rveltz]

The only con is the lack of general eigensolver like in KrylovKit

[amontoison]

Yes, you're totally right. I should implement some eigensolvers one day.

Off-topic: Romain, I will probably come to see JB-Caillau at INRIA the second week of January.

[rveltz]

Great! Please ping me the week before so I can free some time.

[rveltz]

By the way, thanks for your detailed answers, I am typing on my phone so I am not really verbose