Comments (13)
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 Array
s 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.
Related Issues (20)
- eigsolve not working with howmany option HOT 8
- "Access to undefined reference" when using a `Vector` in exponentiate HOT 5
- Free vibration problem fails to converge HOT 1
- [Question] Eigenvalue from eigsolve by Lanczos is not correct HOT 2
- geneigsolve not correctly identifying positive definite and hermitian matrices HOT 1
- `DomainError` in `svdsolve` when running multithreaded code
- Access to undefined reference when using KrylovKit on a non-concrete vector HOT 4
- GMRES Fails Silently From Stagnation HOT 4
- Feature request: Multiple timesteps in Exponentiate HOT 1
- Convergence issue of `exponentiate` function HOT 1
- GEP of free vibration still fails; different error
- Multithreaded GPU Thread Errors
- Standard evp (non-symm, real matrix): too many iterations HOT 2
- eigsolve fails to converge at times HOT 11
- citation HOT 2
- Feature request: Support multiple input vectors for `eigsolve`
- WARNING: using ChainRulesCore.add!! in module KrylovKit conflicts with an existing identifier. HOT 3
- Feature request: block lanczos for solving degenerate ground states HOT 2
- re-export VectorInterface? HOT 10
- Weird behavior of exponentiate when trying to use inplace behavior HOT 4
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from krylovkit.jl.