specify dot product

[rveltz]

Hi,

Is it possible to pass a dot product in eigsolve ?

This is useful for me because I know a dot product for which the operator is hermitian.

Thank you for your help,

Best

[Jutho]

Not directly, but there is a simple wrapper type available in KrylovKit, though it's not really documented: https://github.com/Jutho/KrylovKit.jl/blob/master/src/innerproductvec.jl

So just call InnerProductVec(your_vector, your_dot_function) and this should work automatically. your_vector can still be any type of vector that KrylovKit accepts.

[rveltz]

Thank you for your prompt answer.

But then I pass this struct to what?

[Jutho]

Just use this new vector of the InnerProductVec type in eigsolve, together with wathever linear operator you had (matrix, sparse matrix, arbitrary function). At the end of the calculation, extract the vector(s) by saying (untested, but something like this should work)

eigvals, eigvecs, ... = eigsolve(...)
eigvecs = getfield.(eigvecs, :vec)

[rveltz]

Thank you!

[Jutho]

Let me know if you run into any problems, cause this functionality is not well tested.

[rveltz]

OK. Here is what I used

using KrylovKit
@time KrylovKit.eigsolve(x -> J(x), AF(rand(TY,size(Iext))), 5,:LR;issymmetric = false, verbosity = 2, krylovdim = 25, tol = 1e-5, maxiter = 10) #~80seconds

@time KrylovKit.eigsolve(x -> J(x), AF(rand(TY,size(Iext))), 5,:LR;issymmetric = true, verbosity = 2, krylovdim = 25, tol = 1e-5, maxiter = 10) #~65seconds

newVecDot = KrylovKit.InnerProductVec(AF(rand(TY,size(Iext))), dot)
@time KrylovKit.eigsolve(x -> KrylovKit.InnerProductVec(J(x.vec), dot) , newVecDot, 5,:LR;issymmetric = true, verbosity = 2, krylovdim = 25, tol = 1e-5, maxiter = 10) #19 ~78 seconds

I seem to lose the benefit of using issymmetric = true but the eigenvalues seem the same

[Jutho]

If the eigenvalues are reported as real numbers (i.e. Float64 instead of ComplexF64) it does use the symmetric (i.e. Lanczos) algorithm. But maybe some other problem is at play; this would require benchmarking. Norm is computed as the square root of dot for InnerProductVec. What is AF in your code? For plain Arrays as vectors, certain operations are computed more efficiently.

[rveltz]

AF = Array{Float64}

The differences in timing are due to using InnerProductVec instead of Array{Float64}. It seems a big difference. Note that size(Iext) = (256, 256, 64) so wrapping time in InnerProductVec should be negligeable.

[Jutho]

Do you have multithreading enabled? As said, Array follows a different code path and tries to apply BLAS level 2 like algorithms for doing the orthogonalisation with respect to the existing Krylov basis, and I guess this is lost when wrapping them in InnerProductVec.

[rveltz]

Yes! Thank you for the hint. So I guess it might be better (for speed reasons) to code another version of the linear operator and keep the canonical scalar product.

[Jutho]

Or make a custom wrapper, that is a subtype of AbstractArray and has linear indexing (these are the properties that determine the fast code path), i.e. it should have (again, untested)

struct MyInnerProductVec{T,N,F,A<:AbstractArray{T,N}} <: AbstractArray{T,N}
        vec::A
        dotf::F
end
# everything from innerproductvec.jl
Base.IndexStyle(::Type{<:MyInnerProductVec}) = IndexLinear()
Base.size(a::MyInnerProductVec, I...) = size(a.vec, I...)
Base.@propagate_inbounds Base.getindex(a::MyInnerProductVec, I...) = getindex(a.vec, I...)
Base.@propagate_inbounds Base.setindex!(a::MyInnerProductVec, v, I...) = setindex!(a.vec, v, I...)

[Jutho]

Did this work or did you not have time to try?

[rveltz]

Sorry, I have not found the time to try yet