eigvals fails on a pair of complex matrices of order 40

[andreasvarga]

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]

Looks similar to #40260

[jarlebring]

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]

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]

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

[andreasvarga]

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 <#40492 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ALJDHECJ3HMYOJOWV5FLGCDTMK7J5ANCNFSM427DMQTQ> .

[stevengj]

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]

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]

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]

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]

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.