**Purpose**

To compute orthogonal transformation matrices Q and Z which reduce the regular pole pencil A-lambda*E of the descriptor system (A-lambda*E,B,C) to the generalized real Schur form with ordered generalized eigenvalues. The pair (A,E) is reduced to the form ( A1 * * ) ( E1 * * ) Q'*A*Z = ( 0 A2 * ) , Q'*E*Z = ( 0 E2 * ) , (1) ( 0 0 A3 ) ( 0 0 E2 ) where the subpencils Ak-lambda*Ek, for k = 1, 2, 3, contain the generalized eigenvalues which belong to certain domains of interest.

SUBROUTINE TG01QD( DICO, STDOM, JOBFI, N, M, P, ALPHA, A, LDA, $ E, LDE, B, LDB, C, LDC, N1, N2, N3, ND, NIBLCK, $ IBLCK, Q, LDQ, Z, LDZ, ALPHAR, ALPHAI, BETA, $ TOL, IWORK, DWORK, LDWORK, INFO ) C .. Scalar Arguments .. CHARACTER DICO, JOBFI, STDOM INTEGER INFO, LDA, LDB, LDC, LDE, LDQ, LDWORK, LDZ, M, N, $ N1, N2, N3, ND, NIBLCK, P DOUBLE PRECISION ALPHA, TOL C .. Array Arguments .. INTEGER IBLCK( * ), IWORK(*) DOUBLE PRECISION A(LDA,*), ALPHAI(*), ALPHAR(*), B(LDB,*), $ BETA(*), C(LDC,*), DWORK(*), E(LDE,*), Q(LDQ,*), $ Z(LDZ,*)

**Mode Parameters**

DICO CHARACTER*1 Specifies the type of the descriptor system as follows: = 'C': continuous-time system; = 'D': discrete-time system. STDOM CHARACTER*1 Specifies the type of the domain of interest for the generalized eigenvalues, as follows: = 'S': stability type domain (i.e., left part of complex plane or inside of a circle); = 'U': instability type domain (i.e., right part of complex plane or outside of a circle); = 'N': whole complex domain, excepting infinity. JOBFI CHARACTER*1 Specifies the type of generalized eigenvalues in the leading diagonal block(s) as follows: = 'F': finite generalized eigenvalues are in the leading diagonal blocks (Af,Ef), and the resulting transformed pair has the form ( Af * ) ( Ef * ) Q'*A*Z = ( ) , Q'*E*Z = ( ) ; ( 0 Ai ) ( 0 Ei ) = 'I': infinite generalized eigenvalues are in the leading diagonal blocks (Ai,Ei), and the resulting transformed pair has the form ( Ai * ) ( Ei * ) Q'*A*Z = ( ) , Q'*E*Z = ( ) . ( 0 Af ) ( 0 Ef )

N (input) INTEGER The number of rows of the matrix B, the number of columns of the matrix C and the order of the square matrices A and E. N >= 0. M (input) INTEGER The number of columns of the matrix B. M >= 0. P (input) INTEGER The number of rows of the matrix C. P >= 0. ALPHA (input) DOUBLE PRECISION The boundary of the domain of interest for the finite generalized eigenvalues of the pair (A,E). For a continuous-time system (DICO = 'C'), ALPHA is the boundary value for the real parts of the generalized eigenvalues, while for a discrete-time system (DICO = 'D'), ALPHA >= 0 represents the boundary value for the moduli of the generalized eigenvalues. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the leading N-by-N part of this array must contain the state dynamics matrix A. On exit, the leading N-by-N part of this array contains the matrix Q'*A*Z in real Schur form, with the elements below the first subdiagonal set to zero. If JOBFI = 'I', the N1-by-N1 pair (A1,E1) contains the infinite spectrum, the N2-by-N2 pair (A2,E2) contains the finite spectrum in the domain of interest, and the N3-by-N3 pair (A3,E3) contains the finite spectrum ouside of the domain of interest. If JOBFI = 'F', the N1-by-N1 pair (A1,E1) contains the finite spectrum in the domain of interest, the N2-by-N2 pair (A2,E2) contains the finite spectrum ouside of the domain of interest, and the N3-by-N3 pair (A3,E3) contains the infinite spectrum. Ai has a block structure as in (2), where A0,0 is ND-by-ND and Ai,i is IBLCK(i)-by-IBLCK(i), for i = 1, ..., NIBLCK. The domain of interest for the pair (Af,Ef), containing the finite generalized eigenvalues, is defined by the parameters ALPHA, DICO and STDOM as follows: For DICO = 'C': Real(eig(Af,Ef)) < ALPHA if STDOM = 'S'; Real(eig(Af,Ef)) > ALPHA if STDOM = 'U'. For DICO = 'D': Abs(eig(Af,Ef)) < ALPHA if STDOM = 'S'; Abs(eig(Af,Ef)) > ALPHA if STDOM = 'U'. LDA INTEGER The leading dimension of the array A. LDA >= MAX(1,N). E (input/output) DOUBLE PRECISION array, dimension (LDE,N) On entry, the leading N-by-N part of this array must contain the descriptor matrix E. On exit, the leading N-by-N part of this array contains the matrix Q'*E*Z in upper triangular form, with the elements below the diagonal set to zero. Its structure corresponds to the block structure of the matrix Q'*A*Z (see description of A). LDE INTEGER The leading dimension of the array E. LDE >= MAX(1,N). B (input/output) DOUBLE PRECISION array, dimension (LDB,M) On entry, the leading N-by-M part of this array must contain the input matrix B. On exit, the leading N-by-M part of this array contains the transformed input matrix Q'*B. LDB INTEGER The leading dimension of the array B. LDB >= MAX(1,N). C (input/output) DOUBLE PRECISION array, dimension (LDC,N) On entry, the leading P-by-N part of this array must contain the output matrix C. On exit, the leading P-by-N part of this array contains the transformed output matrix C*Z. LDC INTEGER The leading dimension of the array C. LDC >= MAX(1,P). N1 (output) INTEGER N2 (output) INTEGER N3 (output) INTEGER The number of the generalized eigenvalues of the pairs (A1,E1), (A2,E2) and (A3,E3), respectively. ND (output) INTEGER. The number of non-dynamic infinite eigenvalues of the pair (A,E). Note: N-ND is the rank of the matrix E. NIBLCK (output) INTEGER If ND > 0, the number of infinite blocks minus one. If ND = 0, then NIBLCK = 0. IBLCK (output) INTEGER array, dimension (N) IBLCK(i) contains the dimension of the i-th block in the staircase form (2), where i = 1,2,...,NIBLCK. Q (output) DOUBLE PRECISION array, dimension (LDQ,N) The leading N-by-N part of this array contains the orthogonal matrix Q, where Q' is the product of orthogonal transformations applied to A, E, and B on the left. LDQ INTEGER The leading dimension of the array Q. LDQ >= MAX(1,N). Z (output) DOUBLE PRECISION array, dimension (LDZ,N) The leading N-by-N part of this array contains the orthogonal matrix Z, which is the product of orthogonal transformations applied to A, E, and C on the right. LDZ INTEGER The leading dimension of the array Z. LDZ >= MAX(1,N). ALPHAR (output) DOUBLE PRECISION array, dimension (N) ALPHAI (output) DOUBLE PRECISION array, dimension (N) BETA (output) DOUBLE PRECISION array, dimension (N) On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j = 1, ..., N, are the generalized eigenvalues. ALPHAR(j) + ALPHAI(j)*i, and BETA(j), j = 1, ..., N, are the diagonals of the complex Schur form (S,T) that would result if the 2-by-2 diagonal blocks of the real Schur form of (A,E) were further reduced to triangular form using 2-by-2 complex unitary transformations. If ALPHAI(j) is zero, then the j-th eigenvalue is real; if positive, then the j-th and (j+1)-st eigenvalues are a complex conjugate pair, with ALPHAI(j+1) negative.

TOL DOUBLE PRECISION A tolerance used in rank decisions to determine the effective rank, which is defined as the order of the largest leading (or trailing) triangular submatrix in the QR factorization with column pivoting whose estimated condition number is less than 1/TOL. If the user sets TOL <= 0, then an implicitly computed, default tolerance TOLDEF = N**2*EPS, is used instead, where EPS is the machine precision (see LAPACK Library routine DLAMCH). TOL < 1.

IWORK INTEGER array, dimension (N) DWORK DOUBLE PRECISION array, dimension (LDWORK) On exit, if INFO = 0, DWORK(1) returns the optimal value of LDWORK. LDWORK INTEGER The length of the array DWORK. LDWORK >= 1, and if N = 0, LDWORK >= 4*N, if STDOM = 'N'; LDWORK >= 4*N+16, if STDOM = 'S' or 'U'. For optimum performance LDWORK should be larger. If LDWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the DWORK array, returns this value as the first entry of the DWORK array, and no error message related to LDWORK is issued by XERBLA.

INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value; = 1: the pencil A-lambda*E is not regular; = 2: the QZ algorithm failed to compute all generalized eigenvalues of the pair (A,E); = 3: a failure occured during the ordering of the generalized real Schur form of the pair (A,E).

The separation of the finite and infinite parts is based on the reduction algorithm of [1]. If JOBFI = 'F', the matrices of the pair (Ai,Ei), containing the infinite generalized eigenvalues, have the form ( A0,0 A0,k ... A0,1 ) ( 0 E0,k ... E0,1 ) Ai = ( 0 Ak,k ... Ak,1 ) , Ei = ( 0 0 ... Ek,1 ) ; (2) ( : : . : ) ( : : . : ) ( 0 0 ... A1,1 ) ( 0 0 ... 0 ) if JOBFI = 'I', the matrices Ai and Ei have the form ( A1,1 ... A1,k A1,0 ) ( 0 ... E1,k E1,0 ) Ai = ( : . : : ) , Ei = ( : . : : ) , (3) ( : ... Ak,k Ak,0 ) ( : ... 0 Ek,0 ) ( 0 ... 0 A0,0 ) ( 0 ... 0 0 ) where Ai,i , for i = 0, 1, ..., k, are nonsingular upper triangular matrices, and A0,0 corresponds to the non-dynamic infinite modes of the system.

[1] Misra, P., Van Dooren, P., and Varga, A. Computation of structural invariants of generalized state-space systems. Automatica, 30, pp. 1921-1936, 1994.

3 The algorithm requires about 25N floating point operations.

The number of infinite poles is computed as NIBLCK NINFP = Sum IBLCK(i) = N - ND - NF, i=1 where NF is the number of finite generalized eigenvalues. The multiplicities of infinite poles can be computed as follows: there are IBLCK(k)-IBLCK(k+1) infinite poles of multiplicity k, for k = 1, ..., NIBLCK, where IBLCK(NIBLCK+1) = 0. Note that each infinite pole of multiplicity k corresponds to an infinite eigenvalue of multiplicity k+1.

**Program Text**

* TG01QD EXAMPLE PROGRAM TEXT * Copyright (c) 2002-2017 NICONET e.V. * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D0 ) INTEGER NIN, NOUT PARAMETER ( NIN = 5, NOUT = 6 ) INTEGER NMAX, MMAX, PMAX PARAMETER ( NMAX = 20, MMAX = 20, PMAX = 20 ) INTEGER LDA, LDB, LDC, LDE, LDQ, LDZ PARAMETER ( LDA = NMAX, LDB = NMAX, LDC = PMAX, $ LDE = NMAX, LDQ = NMAX, LDZ = NMAX ) INTEGER LDWORK PARAMETER ( LDWORK = 4*NMAX+16 ) * .. Local Scalars .. CHARACTER*1 DICO, JOBFI, STDOM INTEGER I, INFO, J, M, N, N1, N2, N3, ND, NIBLCK, P DOUBLE PRECISION ALPHA, TOL * .. Local Arrays .. INTEGER IBLCK(NMAX), IWORK(NMAX) DOUBLE PRECISION A(LDA,NMAX), ALPHAI(NMAX), ALPHAR(NMAX), $ B(LDB,MMAX), BETA(NMAX), C(LDC,NMAX), $ DWORK(LDWORK), E(LDE,NMAX), Q(LDQ,NMAX), $ Z(LDZ,NMAX) * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL TG01QD * .. Intrinsic Functions .. INTRINSIC DCMPLX * .. Executable Statements .. * WRITE ( NOUT, FMT = 99999 ) * Skip the heading in the data file and read the data. READ ( NIN, FMT = '()' ) READ ( NIN, FMT = * ) N, M, P, DICO, STDOM, JOBFI, ALPHA, TOL IF ( N.LT.0 .OR. N.GT.NMAX ) THEN WRITE ( NOUT, FMT = 99988 ) N ELSE READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N ) READ ( NIN, FMT = * ) ( ( E(I,J), J = 1,N ), I = 1,N ) IF ( M.LT.0 .OR. M.GT.MMAX ) THEN WRITE ( NOUT, FMT = 99987 ) M ELSE READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,M ), I = 1,N ) IF ( P.LT.0 .OR. P.GT.PMAX ) THEN WRITE ( NOUT, FMT = 99986 ) P ELSE READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,N ), I = 1,P ) * Find the reduced descriptor system * (A-lambda E,B,C). CALL TG01QD( DICO, STDOM, JOBFI, N, M, P, ALPHA, A, LDA, $ E, LDE, B, LDB, C, LDC, N1, N2, N3, ND, $ NIBLCK, IBLCK, Q, LDQ, Z, LDZ, ALPHAR, $ ALPHAI, BETA, TOL, IWORK, DWORK, LDWORK, $ INFO ) * IF ( INFO.NE.0 ) THEN WRITE ( NOUT, FMT = 99998 ) INFO ELSE WRITE ( NOUT, FMT = 99994 ) N1, N2, N3 WRITE ( NOUT, FMT = 99983 ) ND WRITE ( NOUT, FMT = 99989 ) NIBLCK + 1 IF ( NIBLCK.GT.0 ) THEN WRITE ( NOUT, FMT = 99985 ) $ ( IBLCK(I), I = 1, NIBLCK ) END IF WRITE ( NOUT, FMT = 99997 ) DO 10 I = 1, N WRITE ( NOUT, FMT = 99995 ) ( A(I,J), J = 1,N ) 10 CONTINUE WRITE ( NOUT, FMT = 99996 ) DO 20 I = 1, N WRITE ( NOUT, FMT = 99995 ) ( E(I,J), J = 1,N ) 20 CONTINUE WRITE ( NOUT, FMT = 99993 ) DO 30 I = 1, N WRITE ( NOUT, FMT = 99995 ) ( B(I,J), J = 1,M ) 30 CONTINUE WRITE ( NOUT, FMT = 99992 ) DO 40 I = 1, P WRITE ( NOUT, FMT = 99995 ) ( C(I,J), J = 1,N ) 40 CONTINUE WRITE ( NOUT, FMT = 99991 ) DO 50 I = 1, N WRITE ( NOUT, FMT = 99995 ) ( Q(I,J), J = 1,N ) 50 CONTINUE WRITE ( NOUT, FMT = 99990 ) DO 60 I = 1, N WRITE ( NOUT, FMT = 99995 ) ( Z(I,J), J = 1,N ) 60 CONTINUE WRITE ( NOUT, FMT = 99985 ) DO 70 I = 1, N IF ( BETA(I).EQ.ZERO .OR. ALPHAI(I).EQ.ZERO ) THEN WRITE ( NOUT, FMT = 99984 ) $ ALPHAR(I)/BETA(I) ELSE WRITE ( NOUT, FMT = 99984 ) $ DCMPLX( ALPHAR(I), ALPHAI(I) )/BETA(I) END IF 70 CONTINUE END IF END IF END IF END IF STOP * 99999 FORMAT (' TG01QD EXAMPLE PROGRAM RESULTS',/1X) 99998 FORMAT (' INFO on exit from TG01QD = ',I2) 99997 FORMAT (/' The transformed state dynamics matrix Q''*A*Z is ') 99996 FORMAT (/' The transformed descriptor matrix Q''*E*Z is ') 99995 FORMAT (20(1X,F8.4)) 99994 FORMAT (' Number of eigenvalues of the three blocks =', 3I5) 99993 FORMAT (/' The transformed input/state matrix Q''*B is ') 99992 FORMAT (/' The transformed state/output matrix C*Z is ') 99991 FORMAT (/' The left transformation matrix Q is ') 99990 FORMAT (/' The right transformation matrix Z is ') 99989 FORMAT ( ' Number of infinite blocks = ',I5) 99988 FORMAT (/' N is out of range.',/' N = ',I5) 99987 FORMAT (/' M is out of range.',/' M = ',I5) 99986 FORMAT (/' P is out of range.',/' P = ',I5) 99985 FORMAT (/' The finite generalized eigenvalues are '/ $ ' real part imag part ') 99984 FORMAT (1X,F9.4,SP,F9.4,S,'i ') 99983 FORMAT ( ' Number of non-dynamic infinite eigenvalues = ',I5) END

TG01QD EXAMPLE PROGRAM DATA 4 2 2 C S F -1.E-7 0.0 -1 0 0 3 0 0 1 2 1 1 0 4 0 0 0 0 1 2 0 0 0 1 0 1 3 9 6 3 0 0 2 0 1 0 0 0 0 1 1 1 -1 0 1 0 0 1 -1 1

TG01QD EXAMPLE PROGRAM RESULTS Number of eigenvalues of the three blocks = 1 2 1 Number of non-dynamic infinite eigenvalues = 1 Number of infinite blocks = 1 The transformed state dynamics matrix Q'*A*Z is 1.6311 2.1641 -0.5848 -3.6517 0.0000 -0.4550 1.0935 1.6851 0.0000 0.0000 0.0000 1.5770 0.0000 0.0000 0.0000 2.2913 The transformed descriptor matrix Q'*E*Z is -0.4484 9.6340 5.1601 -2.6183 0.0000 -3.3099 -1.6050 -0.2756 0.0000 0.0000 2.3524 -0.6008 0.0000 0.0000 0.0000 0.0000 The transformed input/state matrix Q'*B is 0.0232 -0.9413 0.7251 0.2478 -0.4336 -0.9538 1.1339 0.3780 The transformed state/output matrix C*Z is 0.8621 0.3754 -0.8847 0.5774 0.1511 -1.1192 1.1795 0.5774 The left transformation matrix Q is 0.0232 0.7251 0.3902 0.5669 -0.3369 -0.6425 0.3902 0.5669 -0.9413 0.2478 -0.1301 -0.1890 0.0000 0.0000 -0.8238 0.5669 The right transformation matrix Z is -0.8621 -0.3754 -0.0843 -0.3299 0.4258 -0.9008 -0.0211 -0.0825 0.0000 0.0000 -0.9689 0.2474 -0.2748 -0.2184 0.2317 0.9073 The finite generalized eigenvalues are real part imag part -3.6375 0.1375 0.0000 Infinity