JuliaMath / DoubleFloats.jl

math with more good bits
MIT License
153 stars 33 forks source link

Matrix exponential throws error #127

Open judober opened 3 years ago

judober commented 3 years ago

I tried calculating a matrix exponential but

julia> exp([1 2; 3 4])
2×2 Matrix{Float64}:
  51.969   74.7366
 112.105  164.074

julia> exp(Double64.([1 2; 3 4]))
ERROR: MethodError: no method matching svdvals!(::Matrix{Complex{Double64}})

I thought it worked with an older version, so I tried (randomly) v1.1.15 and it worked (v1.1.21 did not work):

julia> exp(Double64.([1 2; 3 4]))
2×2 Matrix{Complex{Double64}}:
    51.968956198705 + 0.0im   74.73656456700321 + 0.0im
 112.10484685050481 + 0.0im  164.07380304920983 + 0.0im

However, for some reason it tourned out complex.

Update: I did some more tests. v1.1.18 ist the last version that works (with complex output). In v1.1.19 and v1.1.20 eigen! seems missing. v1.1.18 has Generic SVD in its dependencies that, given the name, might be connected.

JeffreySarnoff commented 3 years ago

fixed-ish


julia> using DoubleFloats

julia> using LinearAlgebra

julia> S  = convert.(Double64, [1 2 0; 2 1 3; 0 3 1]) #Symmetric Matrix
3×3 Matrix{Double64}:
 1.0  2.0  0.0
 2.0  1.0  3.0
 0.0  3.0  1.0

julia> eigvals(S) #Works
3-element Vector{Complex{Double64}}:
 -2.605551275463989 + 0.0im
                1.0 + 0.0im
   4.60555127546399 + 0.0im

julia> eigvecs(S) 
3×3 Matrix{Complex{Double64}}:
 -0.3922322702763681 + 0.0im     0.8320502943378437 + 0.0im  0.3922322702763681 + 0.0im
  0.7071067811865476 + 0.0im  1.288877052411983e-32 + 0.0im  0.7071067811865476 + 0.0im
 -0.5883484054145521 + 0.0im    -0.5547001962252291 + 0.0im  0.5883484054145521 + 0.0im
JeffreySarnoff commented 3 years ago

thanks for the information .. I will give this more attention later

judober commented 3 years ago

I tried again with the new releases. With add DoubleFloats #v1.1.22 I get the same error for my original example. With add DoubleFloats #v1.1.23 I get

ERROR: MethodError: no method matching eigen!(::Matrix{Double64}; permute=true, scale=true, sortby=LinearAlgebra.eigsortby)

probably the same issue as #126. Adding GenericSchur to v1.1.23 fixes the problem.

Update: Actually, import GenericSchur.eigen! is sufficient and might be a solution for the other problems with GenericSchur that are referenced in #126?