Extended decompositions (time-varying matrix size)

Thanks for fixing this issue.

I performed some tests including also the reordering of eigenvalues. I implemented a wrapper to the SLICOT routine MB03KD, which served as basis for tests. In solving small dimension periodic Sylvester equations, this routine fully exploits the block structure of the Kronecker expansion using a custom QR-decomposition and solve algorithm. Therefore, it is very fast (probably at the limit of achievable performance). For your information, here are the timing results obtained running the SLICOT wrappers based code to compute and reorder the PSD and the codes in the PeriodicSchurDecompositions package.

using PeriodicSystems 
using PeriodicSchurDecompositions
using JLD
hpd = load("test.jld","hpd");

@time S, Z, ev, sind, = PeriodicSystems.pschur(hpd);
select1 = abs.(ev) .< 1
@time PeriodicSystems.psordschur!(S, Z, select1; sind);

@time PSF = PeriodicSchurDecompositions.pschur(hpd,:L); 
select = abs.(PSF.values) .< 1
@time PeriodicSchurDecompositions.ordschur!(PSF,select);

julia> @time S, Z, ev, sind, = PeriodicSystems.pschur(hpd);
  0.004264 seconds (5.85 k allocations: 620.297 KiB)

julia> select1 = abs.(ev) .< 1
8-element BitVector:
 0
 0
 0
 0
 1
 1
 1
 1

julia> @time PeriodicSystems.psordschur!(S, Z, select1; sind);
  0.002052 seconds (27.17 k allocations: 1.520 MiB)

julia> @time PSF = PeriodicSchurDecompositions.pschur(hpd,:L);
  0.029634 seconds (216.12 k allocations: 21.850 MiB)

julia> select = abs.(PSF.values) .< 1
8-element BitVector:
 0
 0
 0
 0
 1
 1
 1
 1

julia> @time PeriodicSchurDecompositions.ordschur!(PSF,select);
  2.313631 seconds (194.97 k allocations: 6.110 GiB, 21.71% gc time)

The requested reordering corresponds to the worst case (i.e., moving the 4 trailing eigenvalues to the front).

I included your code in the continuous-time periodic Riccati solvers as an option to the SLICOT-based code to compute the PSD. Here are the timing results:

julia>     @time X, EVALS, F = prcric(psysc.A, psysc.B, r, q; K = 200, solver = "symplectic", reltol = 1.e-10, abstol = 1.e-10, fast = false, PSD_SLICOT = true);
  0.626887 seconds (3.61 M allocations: 522.883 MiB, 11.34% gc time)

julia>     @time X, EVALS, F = prcric(psysc.A, psysc.B, r, q; K = 200, solver = "symplectic", reltol = 1.e-10, abstol = 1.e-10, fast = false, PSD_SLICOT = false);
  2.981597 seconds (3.98 M allocations: 6.639 GiB, 18.63% gc time)

I implemented several structure exploiting solvers for small periodic Sylvester-type equations. A typical gain of speed for a problem with 8x8 matrices of period 200 was about 10 times. I wonder if it would be possible to use a similar approach for speeding up the Sylvester solvers used for reordering.

RalphAS / PeriodicSchurDecompositions.jl

Extended decompositions (time-varying matrix size) #6