specify dot product about krylovkit.jl HOT 13 CLOSED

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.

from krylovkit.jl.

Thank you for your prompt answer.

But then I pass this struct to what?

from krylovkit.jl.

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)

from krylovkit.jl.

Thank you!

from krylovkit.jl.

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

from krylovkit.jl.

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

from krylovkit.jl.

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.

from krylovkit.jl.

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.

from krylovkit.jl.

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.

from krylovkit.jl.

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.

from krylovkit.jl.

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...)

from krylovkit.jl.

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

from krylovkit.jl.

Sorry, I have not found the time to try yet

from krylovkit.jl.