Exception in tgsyl! subfunction

[chjaesch]

The subfunction tgsyl! throw an exception while applying the function

fiblkdiag(A, E, missing, missing; fast=false, finite_infinite = true)

to the regular pencil

A = [-8.828446764485715e7 0.0 0.0 0.0 1.0 0.0; 0.0 -599202.6439090811 0.0 0.0 1.0 0.0; 0.0 0.0 -44599.09517911089 3.1081161629028176e6 2.0 0.0; 0.0 0.0 -3.1081161629028176e6 -44599.09517911089 0.0 0.0; 0.0 0.0 0.0 0.0 -1.0 0.0; 1.4376687293449677e10 33249.50029419384 29200.436106193243 22920.479599602455 -161.90134290419934 -1.0]

E = [1.0 0.0 0.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 -1.8422498299893581e-6 0.0]

The same happens in the sylvsys function (package MatrixEquations) when computing the Weierstrass canonical form as in

http://www.control.isy.liu.se/research/reports/2004/2602.pdf

while using the fischur function for diagonalisation.

Do you have any suggestions for solving this problem?

[andreasvarga]

Yes, I checked and tgsyl! fails on this example. This function is a wrapper to the LAPACK subroutine dtgsyl.f, so it is out my control. Some scaling of data could be helpful, but unfortunately I did not implemented any of them.

However, just for my curiosity: your data are very badly scaled. Such data usually originate from a particular (wrong) setting of tolerances in previous computations, which are probably atol1 (the absolute tolerance for the elements of A) and atol2 (the same for E). I would suggest to try with other settings if possible. Just looking to the eigenvalues, (4 finite, 2 infinite), I can imagine that with another setting of tolerances this could easily become (3 finite, 3 infinite) which is quite different.

Note: I tried the same computatios in MATLAB using the DescriptorSystemTools. The result was the same:

[At,Et,Bt,Ct,DIMS,TAU,Q,Z] = sl_gsep(4,A,E,zeros(6,0),zeros(0,6),1,1.e-10)
Error using sl_gsep
The diagonal blocks have close eigenvalues

[chjaesch]

Okay, thank you for the quick feedback.
I try to transform the stats of the descriptor system to get better scaled matrices.

If you like to extend your package, it would be nice if you also provide a function which compute the Weierstrass canonical form.

[andreasvarga]

Weierstrass canonical form, as well as Jordan and Kronecker canonical forms are not necessary for numerical computations. Their computations are extremely sensitive numerical problems and should be generally avoided.