Accuracy loss for small magnitude eigenvalues

[andreasvarga]

This is basically a similar example to that in #4, but for a symplectic product of K = 300 8x8 matrices (instead K = 200) (the underlying Riccati equation is the same, but the number of discretization steps K is different, i.e. 300 vs. 200).

The first computation, to serve as reference, has been performed with the SLICOT wrapper pschur available in the PeriodicSystems and produces the following eigenvalues and their reciprocals:

using PeriodicSystems
using JLD
hpd = load("test1.jld","hpd");
S, Z, ev, sind, = PeriodicSystems.pschur(hpd); 
[sort(ev,by=abs) sort(1 ./ev,by=abs)] 

julia> [sort(ev,by=abs) sort(1 ./ev,by=abs)]
8×2 Matrix{ComplexF64}:
 -1.92309e-56+4.89783e-56im  -1.92309e-56-4.89783e-56im
 -1.92309e-56-4.89783e-56im  -1.92309e-56+4.89783e-56im
  2.43966e-47+0.0im           2.43966e-47-0.0im
   2.21356e-7+0.0im            2.21356e-7-0.0im
    4.51761e6+0.0im             4.51761e6-0.0im
   4.09892e46+0.0im            4.09892e46-0.0im
  -6.94581e54+1.769e55im      -6.94581e54-1.769e55im
  -6.94581e54-1.769e55im      -6.94581e54+1.769e55im

which clearly illustrates the symplectic nature of the eigenproblem.

The same computation performed with the PeriodicSchurDecompositions produces slightly inaccurate eigenvalues, which do not fullfil the symplectic requirement, as can be seen below:

using PeriodicSchurDecompositions
using JLD
hpd = load("test1.jld","hpd");
PSF = PeriodicSchurDecompositions.pschur(hpd,:L);
ev1 = PSF.values;
[sort(ev1,by=abs) sort(1 ./ev1,by=abs)]

julia> [sort(ev1,by=abs) sort(1 ./ev1,by=abs)]
8×2 Matrix{ComplexF64}:
  1.07037e-56+0.0im       -1.92309e-56-4.89783e-56im
 -1.65442e-56+0.0im       -1.92309e-56+4.89783e-56im
  2.36581e-48-0.0im        2.43966e-47-0.0im
   2.21356e-7-0.0im         2.21356e-7-0.0im
    4.51761e6+0.0im          4.51761e6+0.0im
   4.09892e46+0.0im         4.22688e47+0.0im
  -6.94581e54+1.769e55im   -6.04442e55-0.0im
  -6.94581e54-1.769e55im    9.34258e55-0.0im

This loss of accuracy in the computed small size eigenvalues prevents the success of some tests for solving periodic Riccati differential equations. The computations for K = 100 are OK on Windows and Linux with Julia 1.7, but fail for K = 200 and K = 300 on both Windows and Linux.

I am attaching the data in the zip-file
test1.zip

[RalphAS]

The standard real decomposition code here is mainly a translation of MB03WD, which seems to have a similar problem. I'll try to investigate further.

I have recently done a translation of the generalized code (MB03BD), which I gather you use in PeriodicSystems. The new generalized code also appears to be more robust than my standard one, when acting on all positive signature cases, but testing is still in progress. Look for that in the next release of the package.

[andreasvarga]

Indeed, there is an issue with the accuracy of MB03WD. However, I did some tests with your code and the problem has been apparently corrected by you. If your fix can be easily moved to the Fortran code, Sima would be (probably) delighted to include in MB03WD.

[RalphAS]

I have made deflation criteria more strict, so the real PSD should work for this case as of v0.1.2. I haven't seen failures with the new version, but please reopen the issue if you encounter them.