imaginary t is not working
Roger-luo opened this issue · 2 comments
Roger-luo commented
MWE
using ExponentialUtilities, LinearAlgebra
A = Hermitian(rand(ComplexF64, 4, 4))
v = rand(ComplexF64, 4)
Ks = arnoldi(A, v);
w = similar(v)
expv!(w, -1e2 * im, Ks)Got error:
ERROR: InexactError: Float64(0.0 + 56.73590474494471im)
Stacktrace:
[1] setindex!(::Array{Float64,1}, ::Complex{Float64}, ::Int64) at ./complex.jl:37
[2] macro expansion at /Users/roger/Documents/Repos/julia/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/generic.jl:101 [inlined]
[3] macro expansion at ./simdloop.jl:73 [inlined]
[4] lmul! at /Users/roger/Documents/Repos/julia/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/generic.jl:100 [inlined]
[5] #expv!#20(::Nothing, ::Function, ::Array{Complex{Float64},1}, ::Complex{Float64}, ::KrylovSubspace{Float64,Complex{Float64},Float64}) at /Users/roger/.julia/dev/ExponentialUtilities/src/krylov_phiv.jl:86
[6] expv!(::Array{Complex{Float64},1}, ::Complex{Float64}, ::KrylovSubspace{Float64,Complex{Float64},Float64}) at /Users/roger/.julia/dev/ExponentialUtilities/src/krylov_phiv.jl:73
[7] top-level scope at none:0
This is because when A is hermitian, F is always real, but it uses lmul! on F.vectors here:
Is there anyway to avoid the allocation while supporting complex valued t ? Or just dispatch expv! to a special version for complex t?
Roger-luo commented
A much simpler one is:
julia> A = rand(4, 4);
julia> v = rand(4);
julia> Ks = arnoldi(A, v);
julia> w = similar(v);
julia> expv!(w, -1e2 * im, Ks)
ERROR: InexactError: Float64(-0.0 - 107.84113054914846im)
Stacktrace:
[1] setindex!(::Array{Float64,2}, ::Complex{Float64}, ::Int64) at ./complex.jl:37
[2] macro expansion at /Users/roger/Documents/Repos/julia/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/generic.jl:101 [inlined]
[3] macro expansion at ./simdloop.jl:73 [inlined]
[4] lmul! at /Users/roger/Documents/Repos/julia/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/generic.jl:100 [inlined]
[5] #expv!#20(::Nothing, ::Function, ::Array{Float64,1}, ::Complex{Float64}, ::KrylovSubspace{Float64,Float64,Float64}) at /Users/roger/.julia/dev/ExponentialUtilities/src/krylov_phiv.jl:88
[6] expv!(::Array{Float64,1}, ::Complex{Float64}, ::KrylovSubspace{Float64,Float64,Float64}) at /Users/roger/.julia/dev/ExponentialUtilities/src/krylov_phiv.jl:73
[7] top-level scope at none:0