search

JuliaLang / julia

The Julia Programming Language
https://julialang.org/
MIT License
45.47k stars 5.46k forks source link

eigvals fails on a pair of complex matrices of order 40 #40492

Closed andreasvarga closed 1 month ago

andreasvarga commented 3 years ago

I encountered the following error in computing the eigenvalues of a pair of complex matrices A and B, arising by a similarity transformation of a pair of Hamiltonian-skew Hamiltonian matrices. The eigenvalues should have a certain symmetry: ifλ is a generalized eigenvalue of the pair(A,B), then -conj(λ) is a generalized eigenvalue as well. 

julia> eigvals(A,B)
ERROR: LAPACKException(34)
Stacktrace:
 [1] chklapackerror(ret::Int64)
   @ LinearAlgebra.LAPACK C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\lapack.jl:38
 [2] ggev!(jobvl::Char, jobvr::Char, A::Matrix{ComplexF64}, B::Matrix{ComplexF64})
   @ LinearAlgebra.LAPACK C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\lapack.jl:2321
 [3] eigvals!(A::Matrix{ComplexF64}, B::Matrix{ComplexF64}; sortby::typeof(LinearAlgebra.eigsortby))
   @ LinearAlgebra C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\eigen.jl:551
 [4] eigvals!
   @ C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\eigen.jl:550 [inlined]
 [5] #eigvals#82
   @ C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\eigen.jl:580 [inlined]
 [6] eigvals(A::Matrix{ComplexF64}, B::Matrix{ComplexF64})
   @ LinearAlgebra C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\eigen.jl:579
 [7] top-level scope
   @ REPL[77]:1

To reproduce this computation, I am attaching the saved matrices in the file test_ggev.jld (packed in the zip-file test_ggev.zip) and the following code produces the error with Julia v1.6.0 running on a machine with Windows 10: 

using JLD
A, B = load("test_ggev.jld","A","B");
eigvals(A,B)

For checking purposes, I transfered the matrices to MATLAB and obtained the following eigenvalues, which exhibit the expected symmetry: 

eig(A,B)

ans =

  -7.8949 - 3.1677i
  -7.8950 - 3.1674i
   7.8949 - 3.1677i
   7.8950 - 3.1674i
   7.8949 - 3.1675i
   7.8949 - 3.1675i
  -7.8949 - 3.1675i
  -7.8949 - 3.1675i
  -2.6435 - 1.3271i
  -2.6434 - 1.3269i
  -2.1045 + 0.3119i
  -2.1048 + 0.3122i
   2.1045 + 0.3119i
   2.1048 + 0.3122i
  -0.5757 - 2.2930i
  -0.5760 - 2.2949i
   0.5757 - 2.2930i
   0.5760 - 2.2949i
   2.1046 + 0.3121i
   2.1046 + 0.3121i
   2.6435 - 1.3271i
   2.6434 - 1.3269i
   2.6434 - 1.3270i
   2.6434 - 1.3270i
   0.5758 - 2.2940i
   0.5758 - 2.2940i
   0.1491 + 1.3045i
   0.1491 + 1.3045i
  -0.1491 + 1.3045i
  -0.1491 + 1.3045i
   0.1491 + 1.3045i
   0.1491 + 1.3045i
  -2.6434 - 1.3270i
  -2.6434 - 1.3270i
  -2.1046 + 0.3121i
  -2.1046 + 0.3121i
  -0.5758 - 2.2940i
  -0.5758 - 2.2940i
  -0.1491 + 1.3045i
  -0.1491 + 1.3045i

The following computation in Julia produces approximately the same result:

julia> eigvals(A/B)
40-element Vector{ComplexF64}:
  -7.894987786971432 - 3.1673553694785017im
  -7.894924711201389 - 3.167535553514198im
  -7.894924711198504 - 3.1675355535182783im
  -7.894861512653495 - 3.1677157151387565im
  -2.643461272885884 - 1.3271419108713263im
 -2.6434480473706086 - 1.3270066949327255im
 -2.6434480473705486 - 1.3270066949316872im
  -2.643434600429585 - 1.326871448734545im
 -2.1047552691676503 + 0.31216749748540107im
 -2.1046238787850124 + 0.3120568420550187im
 -2.1046238787848073 + 0.3120568420537384im
                     ⋮
  2.1046238787850733 + 0.31205684205442474im
   2.104623878785536 + 0.31205684205421114im
   2.104755269167242 + 0.3121674974850989im
  2.6434346004304508 - 1.3268714487339714im
  2.6434480473693047 - 1.3270066949322732im
  2.6434480473698527 - 1.3270066949328774im
  2.6434612728870044 - 1.3271419108711606im
   7.894861512656981 - 3.167715715138961im
  7.8949247111950225 - 3.1675355535141736im
   7.894924711201617 - 3.167535553515905im
   7.894987786971231 - 3.167355369480694im
antoine-levitt commented 3 years ago

Looks similar to https://github.com/JuliaLang/julia/issues/40260

jarlebring commented 3 years ago

I cannot reproduce this 

julia> A, B = load("/tmp/test_ggev.jld","A","B");
julia> eigvals(A,B)
40-element Vector{ComplexF64}:
   -7.894987786961872 - 3.167355369477298im
   -7.894924711217828 - 3.16753555352056im
  -7.8949247112013765 - 3.1675355535141883im
   -7.894861512643777 - 3.1677157151376583im
....
    7.894924711197777 - 3.167535553515384im
   7.8949247112438705 - 3.1675355535069936im
    7.894987786948488 - 3.167355369484353im
julia> versioninfo()
Julia Version 1.7.0-DEV.931
Commit 29d5158d27* (2021-04-15 20:13 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-8550U CPU @ 1.80GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, skylake)
Environment:
  JULIA_LOAD_PATH = .:/home/jarl/archive_synced/scripts:
andreasvarga commented 3 years ago

This is the version on which the failure takes place:

julia> versioninfo()
Julia Version 1.6.0
Commit f9720dc2eb (2021-03-24 12:55 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-10700 CPU @ 2.90GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, skylake)
jarlebring commented 3 years ago

Would you mind trying it also on 1.7 nightly? There were some changes in MKL / BLAS cf MKL.jl.

andreasvarga commented 3 years ago

Unfortunately I have no possibility to test on Julia 1.7. Let's hope the error will disapear in 1.7.

Elias Jarlebring @.***> schrieb am Do., 6. Mai 2021, 18:26:

Would you mind trying it also on 1.7 nightly? There were some changes in MKL / BLAS cf MKL.jl.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaLang/julia/issues/40492#issuecomment-833661243, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALJDHECJ3HMYOJOWV5FLGCDTMK7J5ANCNFSM427DMQTQ .

stevengj commented 3 years ago

I can reproduce the failure in Julia 1.6, but it works for me in Julia 1.7:

julia> eigvals(A,B)
40-element Vector{ComplexF64}:
   -7.894987786963967 - 3.1673553694860406im
   -7.894924711213602 - 3.1675355535031704im
   -7.894924711201402 - 3.1675355535142016im
  -7.8948615126459245 - 3.1677157151462993im
  -2.6434612728853852 - 1.3271419108695268im
   -2.643448047371438 - 1.327006694935356im
  -2.6434480473702506 - 1.3270066949335426im
  -2.6434346004295546 - 1.3268714487318745im
   -2.104755269168861 + 0.3121674974841139im
   -2.104623878784103 + 0.31205684205486894im
  -2.1046238787832072 + 0.31205684205592776im
  -2.1044925457954906 + 0.31194640700815734im
  -0.5759667502484053 - 2.2949494777856017im
  -0.5758304237483377 - 2.293954362150927im
   -0.575830423748109 - 2.2939543621510876im
   -0.575692391915598 - 2.2929595029807692im
 -0.14914950669611507 + 1.304475457103618im
 -0.14914876709856803 + 1.3044904884775734im
 -0.14914876709733502 + 1.3044904884781316im
 -0.14914802717631268 + 1.3045055032652189im
   0.1491480271754271 + 1.30450550326327im
  0.14914876709901037 + 1.3044904884778397im
  0.14914876709930827 + 1.3044904884809951im
  0.14914950669459048 + 1.3044754571024302im
   0.5756923919155538 - 2.292959502980632im
   0.5758304237482132 - 2.293954362151071im
   0.5758304237482588 - 2.2939543621510494im
   0.5759667502484193 - 2.2949494777856083im
   2.1044925457939017 + 0.31194640700901366im
    2.104623878784881 + 0.3120568420547012im
   2.1046238787852563 + 0.31205684205454237im
   2.1047552691676366 + 0.3121674974848107im
   2.6434346004317373 - 1.3268714487313003im
   2.6434480473683797 - 1.327006694935784im
   2.6434480473700557 - 1.3270066949329222im
   2.6434612728864564 - 1.3271419108702713im
   7.8948615126366795 - 3.1677157151337028im
    7.894924711199914 - 3.1675355535202785im
    7.894924711236199 - 3.1675355535178897im
    7.894987786952084 - 3.1673553694778707im
andreasvarga commented 2 years ago

I reran this example with Julia LTS version 1.6.5. On Windows, this example continues to fail. But I can confirm that with 1.7.1 it runs correctly.

In the meantime I encountered other examples of complex matrices of various sizes on which eigvals (and of course schur) fail with version 1.6.5. My dilemma is that, since I included the LTS version as compulsory test platform for four packages on which I am working, I am not sure if this is a wise option for me, unless such failures can still happen. I would appreciate hearing your opinions on this issue and especially knowing what perspectives exist to correct this error in 1.6.5.

andreasvarga commented 2 years ago

I think I found the issue where the convergence problems of generalized complex eigenvalues have been discussed in the LAPACK/Julia context #477 (see also #475). I tried several of the test matrices in #475, and yes, they fail on 1.6.5.

Backport issues to Julia 1.6 are discussed in #39216 in relation to Upgrade to OpenBLAS 0.3.13, but for certain reasons, the idea has been rejected. However I had the impression from that discussion that backporting would still be possible.

vchuravy commented 2 years ago

The fix turned out to cause multiple other issues and so being conservative with backporting seemed to be the right choice. Now that these patches have been used for a while an no other issues have cropped up backporting them should be okay.

giordano commented 1 month ago

It's unlikely there will be a new release in the v1.6 series at this point, and v1.10 will become the new LTS soon™️ . I'm going to close this ticket as won't fix.