acroy / Expokit.jl

Julia implementation of EXPOKIT routines
Other
27 stars 11 forks source link

Similar calls to expmv and expm return very different results #4

Closed marcusps closed 9 years ago

marcusps commented 9 years ago

This happens even for 2x2 matrices, and affects complex and real matrices equally

julia> expm(pi/4*[0 1; 1 0])*[1.;0.]
2-element Array{Float64,1}:
 1.32461 
 0.868671

julia> expmv(pi/4,[0 1; 1 0]|>sparse, [1.;0.])
2-element Array{Float64,1}:
 1.19636
 0.65672

julia> expm(-pi/4*1im*[0 1; 1 0])*[1.;0.]
2-element Array{Complex{Float64},1}:
 0.707107+0.0im     
      0.0-0.707107im

julia> expmv(pi/4,-1im*[0 1; 1 0]|>sparse, [1.;0.])
2-element Array{Complex{Float64},1}:
 0.815705+0.0im     
      0.0-0.578469im
acroy commented 9 years ago

There seems to be a problem with tau?

julia> expmv(1.,pi/4*[0 1; 1 0]|>sparse, [1.;0.])
2-element Array{Float64,1}:
 1.32461 
 0.868671

julia> expmv(1.,pi/4*(-1im)*[0 1; 1 0]|>sparse, complex([1.;0.]))
2-element Array{Complex{Float64},1}:
 0.707107+0.0im     
      0.0-0.707107im
acroy commented 9 years ago

I have somehow introduced this in the process of optimizing the code - for aeb3c24b8a483eb92c7a7c17426c654880dcf074 it still works.

acroy commented 9 years ago

Ok, the mutating scale! in line 85 seems to be problematic; changing it to scale makes it work again. Can't investigate further now ...

marcusps commented 9 years ago

I guess the issue is that you don't want to modify hm in that call, and using scale! in combination with slice does just that. I'll submit a pull request with the fix.