Erroneous eigenvalue computation

[andreasvarga]

The following eigenvalue computation is part of solving periodic Riccati equations. The undelying periodic matrix hpd, is a 200 vector of 8x8 state transition matrices, whose left product is a simplectic matrix. Thus, the eigenvalues are symmetrical with respect to the unit circle.

The first computation has been performed with the SLICOT wrapper pschur available in the PeriodicSystems and produces the following eigenvalues:

using PeriodicSystems
using JLD
hpd = load("test.jld","hpd")
S, Z, ev, = PeriodicSystems.pschur(hpd); ev

8-element Vector{ComplexF64}:
   9.473709584260244e71 + 6.144171599316869e71im
   9.473709584260244e71 - 6.144171599316869e71im
   7.833585365370729e44 + 0.0im
    4.644872743033464e6 + 0.0im
  2.1529115119486817e-7 + 0.0im
 1.2765546724226738e-45 + 0.0im
  7.430242860588899e-73 + 4.8188818491875555e-73im
  7.430242860588899e-73 - 4.8188818491875555e-73im

It can be easily checked the symmetry of eigenvalues

julia> [1. ./ev ev]
8×2 Matrix{ComplexF64}:
 7.43024e-73-4.81888e-73im   9.47371e71+6.14417e71im
 7.43024e-73+4.81888e-73im   9.47371e71-6.14417e71im
 1.27655e-45-0.0im           7.83359e44+0.0im
  2.15291e-7-0.0im            4.64487e6+0.0im
   4.64487e6-0.0im           2.15291e-7+0.0im
  7.83359e44-0.0im          1.27655e-45+0.0im
  9.47371e71-6.14417e71im   7.43024e-73+4.81888e-73im
  9.47371e71+6.14417e71im   7.43024e-73-4.81888e-73im

The same computation performed with the PeriodicSchurDecompositions produces wrong eigenvalues:

using PeriodicSchurDecompositions
using JLD
hpd = load("test.jld","hpd")
PSF = PeriodicSchurDecompositions.pschur(hpd,:L); 

julia> PSF.values
8-element Vector{ComplexF64}:
   9.473709584260287e71 + 6.144171599318682e71im
   9.473709584260287e71 - 6.144171599318682e71im
                    0.0 + 0.0im
   3.184085626273608e45 + 0.0im
  -9.313225746154785e-9 + 0.0im
   4.0418154850413655e6 + 0.0im
 2.3165926597911904e-45 + 0.0im
                    0.0 + 0.0im

The next computation would be the reordering of eigenvalues, but this is probably another story. I am attaching the data in a separate file. test.zip

[RalphAS]

Thanks for reporting. This example is challenging the convergence logic; I will need to study it.

[RalphAS]

This appears to be fixed by the revisions leading up to v0.1.1.