clarify Krylov.jl scope

[dpo]

No description provided.

[lkdvos]

I guess the any Datatype is meant with respect to the operator, and not so much the vectors, which still have to be AbstractVector, or am I wrong in this? (Not that I think it matters too much, but that would be the main difference between these packages I think?)

[amontoison]

Yes, you are right @lkdvos.
By "any datatype", Dominique means any real or complex precision (Float32, Float64, ComplexF32, Float64, Float128, ...).

The workspace of each Krylov method is still composed of abstract vectors (Vector{Float64}, CuVector{Float32}, rocVector{Float64}, ...).

The operator A can be anything.

If the user wants a specific type in the workspace, they need to wrap it in a custom wrapper and define a few routines for this type:

struct FakeVector{T, V} <: AbstractVector{T}
    field::V
end

Krylov.knorm(n::Integer, x::FakeVector{T, V}) where {T, V} = ...

Krylov.kdot(n::Integer, x::FakeVector{T, V}, y::FakeVector{T, V}) where {T, V} = ...

Krylov.kscal!(n::Integer, s::T, x::FakeVector{T, V}) where {T, V} = ...

Krylov.kaxpy!(n::Integer, s::T, x::FakeVector{T, V}, y::FakeVector{T, V}) where {T, V} = ...

Krylov.kaxpby!(n::Integer, s::T, x::FakeVector{T, V}, t::T, y::FakeVector{T, V}) where {T, V} = ...

Krylov.kcopy!(n::Integer, y::FakeVector{T, V}, x::FakeVector{T, V}) where {T, V} = ...

Krylov.kfill!(x::FakeVector{T, V}, val::T) where {T, V} = ...

It's the best compromise I've found so far between flexibility and performance.

[Jutho]

This seems very similar (at least in terms of functionality) to the (mutating) interface imposed by VectorInterface.jl, on which also KrylovKit.jl is built and which is a very light-weight dependency. So shamelessly self-promoting, would you not be willing to consider using that interface instead? That would make it easy for users to just implement just one interface, and then gives them the flexibility to try out both Krylov.jl and KrylovKit.jl

Some questions also emerge from your comment:

1. What is the argument n in most of the methods? Is it just the total vector length
2. If FakeVector (or any other type) adheres to this interface, why should it still be a subtype of AbstractVector? Related to 1: does it need to have a Base.length method?
3. What do you use the kfill! method for in Krylov algorithms? Just to fill with zeros? So then all you need is the equivalent of VectorInterface.zerovector!?
4. If that interface is sufficient, does this mean that also in Krylov.jl, you don't use BLAS level 2 methods in the manipulation of the Krylov space vectors (e.g. orthonormalizing them or constructing linear combinations), i.e. you do not store the Krylov basis vectors as the columns of a matrix? I considered this to be a major downside of the KrylovKit approach, that might affect performance, although in practice it probably is not. I guess this is typically not the dominating part anyway, and BLAS level 2 performance can still be reasonably well approached with simple handwritten implementations.

[amontoison]

This seems very similar (at least in terms of functionality) to the (mutating) interface provided by VectorInterface.jl, on which KrylovKit.jl is also built, and which is a very lightweight dependency. So, shamelessly self-promoting, would you consider using that interface instead? This would allow users to implement just one interface and give them the flexibility to try out both Krylov.jl and KrylovKit.jl.

Very interesting package! I was not aware of it. Initially, I designed this interface with routines knorm and kdot to improve dispatch for various arrays. For GPU arrays, dispatch was inefficient for a long time because a generic fallback was used that involved indexing. I managed to target broadcast temporary or better functions. Now, each GPU package has a proper dispatch for dot, norm, etc., but that wasn't the case a few years ago. The issue I might have with VectorInterface.jl is that I also need to work with AbstractMatrix in block m-Krylov methods (like block_gmres in Krylov.jl), and I require another set of routines:
link.

Some questions also arise from your comment:

1. What is the argument n in most of the methods? Is it just the total vector length?

The argument n represents the length of the array. It's only used for Vector{T} where T is a BLASFloat to call the BLAS routines directly, which require this parameter. It avoids recomputing the length at every call.

2. If FakeVector (or any other type) adheres to this interface, why should it still be a subtype of AbstractVector? Related to 1: does it need to have a Base.length method?

Base.length is mostly used for consistency checks, but some methods for rectangular problems require two variants of vectors. Some vectors should be of size n and others of size m, where (m, n) is the size of your abstract linear operator A. As input, you only have one kind of vector, so how can you instantiate the other one from the first? Currently, I retrieve the type of vector b to have a constructor that depends on a size argument. I’m not sure if you’ve encountered this issue here because you seem to handle only square systems.

3. What do you use the kfill! method for in Krylov algorithms? Just to fill with zeros? So then all you need is the equivalent of?

Yes, kfill! dispatches to fill!, but I use it for both vectors and matrices in Krylov and block-Krylov subspaces, respectively. VectorInterface.zerovector! will work for non-block Krylov methods.

4. If that interface is sufficient, does this mean that in Krylov.jl you don't use BLAS level 2 methods for manipulating the Krylov space vectors (e.g., orthonormalizing them or constructing linear combinations), i.e., you don’t store the Krylov basis vectors as the columns of a matrix? I considered this a major downside of the KrylovKit approach, which might affect performance, although in practice it probably doesn't. I suppose this is typically not the dominating factor anyway, and BLAS level 2 performance can still be reasonably well approximated with simple handwritten implementations.

For non-block Krylov methods, you can't use BLAS level 2 operations in many solvers (except for the operator-vector products). The only methods where we could use them are GMRES, FOM, FGMRES, and GPMR. However, the gain is limited because we need to rely on views for other operations, which don’t work well on GPUs. Furthermore, storing the Krylov vectors as a vector of vectors allows us to easily extend the basis if we don’t restart GMRES and similar methods, whereas using a matrix would require resize!. resize! is also inefficient on GPUs.