**Purpose**

Given the descriptor system (A-lambda*E,B,C) with the system matrices A, E and B of the form ( A1 X1 ) ( E1 Y1 ) ( B1 ) A = ( ) , E = ( ) , B = ( ) , ( 0 X2 ) ( 0 Y2 ) ( 0 ) where - B is an L-by-M matrix, with B1 an N1-by-M submatrix, - A is an L-by-N matrix, with A1 an N1-by-N1 submatrix, - E is an L-by-N matrix, with E1 an N1-by-N1 submatrix with LBE nonzero sub-diagonals, this routine reduces the pair (A1-lambda*E1,B1) to the form Qc'*[ B1 A1-lambda*E1 ]*diag(I,Zc) = ( Bc Ac-lambda*Ec * ) ( ) , ( 0 0 Anc-lambda*Enc ) where: 1) the pencil ( Bc Ac-lambda*Ec ) has full row rank NR for all finite lambda and is in a staircase form with _ _ _ _ ( A1,0 A1,1 ... A1,k-1 A1,k ) ( _ _ _ ) ( Bc Ac ) = ( 0 A2,1 ... A2,k-1 A2,k ) , (1) ( ... _ _ ) ( 0 0 ... Ak,k-1 Ak,k ) _ _ _ ( E1,1 ... E1,k-1 E1,k ) ( _ _ ) Ec = ( 0 ... E2,k-1 E2,k ) , (2) ( ... _ ) ( 0 ... 0 Ek,k ) _ where Ai,i-1 is an rtau(i)-by-rtau(i-1) full row rank _ matrix (with rtau(0) = M) and Ei,i is an rtau(i)-by-rtau(i) upper triangular matrix. 2) the pencil Anc-lambda*Enc is regular of order N1-NR with Enc upper triangular; this pencil contains the uncontrollable finite eigenvalues of the pencil (A1-lambda*E1). The transformations are applied to the whole matrices A, E, B and C. The left and/or right orthogonal transformations Qc and Zc, performed to reduce the pencil, can be optionally accumulated in the matrices Q and Z, respectively. The reduced order descriptor system (Ac-lambda*Ec,Bc,Cc) has no uncontrollable finite eigenvalues and has the same transfer- function matrix as the original system (A-lambda*E,B,C).

SUBROUTINE TG01HY( COMPQ, COMPZ, L, N, M, P, N1, LBE, A, LDA, $ E, LDE, B, LDB, C, LDC, Q, LDQ, Z, LDZ, NR, $ NRBLCK, RTAU, TOL, IWORK, DWORK, LDWORK, INFO ) C .. Scalar Arguments .. CHARACTER COMPQ, COMPZ INTEGER INFO, L, LBE, LDA, LDB, LDC, LDE, LDQ, LDWORK, $ LDZ, M, N, N1, NR, NRBLCK, P DOUBLE PRECISION TOL C .. Array Arguments .. INTEGER IWORK( * ), RTAU( * ) DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ), $ DWORK( * ), E( LDE, * ), Q( LDQ, * ), $ Z( LDZ, * )

**Mode Parameters**

COMPQ CHARACTER*1 = 'N': do not compute Q; = 'I': Q is initialized to the unit matrix, and the orthogonal matrix Q is returned; = 'U': Q must contain an orthogonal matrix Q1 on entry, and the product Q1*Q is returned. COMPZ CHARACTER*1 = 'N': do not compute Z; = 'I': Z is initialized to the unit matrix, and the orthogonal matrix Z is returned; = 'U': Z must contain an orthogonal matrix Z1 on entry, and the product Z1*Z is returned.

L (input) INTEGER The number of descriptor state equations; also the number of rows of the matrices A, E and B. L >= 0. N (input) INTEGER The dimension of the descriptor state vector; also the number of columns of the matrices A, E and C. N >= 0. M (input) INTEGER The dimension of descriptor system input vector; also the number of columns of the matrix B. M >= 0. P (input) INTEGER The dimension of descriptor system output; also the number of rows of the matrix C. P >= 0. N1 (input) INTEGER The order of the subsystem (A1-lambda*E1,B1,C1) to be reduced. MIN(L,N) >= N1 >= 0. LBE (input) INTEGER The number of nonzero sub-diagonals of the submatrix E1. MAX(0,N1-1) >= LBE >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the leading L-by-N part of this array must contain the L-by-N state matrix A in the partitioned form ( A1 X1 ) A = ( ) , ( 0 X2 ) where A1 is N1-by-N1. On exit, the leading L-by-N part of this array contains the transformed state matrix, ( Ac * * ) Qc'*A*diag(Zc,I) = ( 0 Anc * ) , ( 0 0 * ) where Ac is NR-by-NR and Anc is (N1-NR)-by-(N1-NR). The matrix ( Bc Ac ) is in the controllability staircase form (1). LDA INTEGER The leading dimension of the array A. LDA >= MAX(1,L). E (input/output) DOUBLE PRECISION array, dimension (LDE,N) On entry, the leading L-by-N part of this array must contain the L-by-N descriptor matrix E in the partitioned form ( E1 Y1 ) E = ( ) , ( 0 Y2 ) where E1 is an N1-by-N1 matrix with LBE nonzero sub-diagonals. On exit, the leading L-by-N part of this array contains the transformed descriptor matrix ( Ec * * ) Qc'*E*diag(Zc,I) = ( 0 Enc * ) , ( 0 0 * ) where Ec is NR-by-NR and Enc is (N1-NR)-by-(N1-NR). Both Ec and Enc are upper triangular. LDE INTEGER The leading dimension of the array E. LDE >= MAX(1,L). B (input/output) DOUBLE PRECISION array, dimension (LDB,M) On entry, the leading L-by-M part of this array must contain the L-by-M input matrix B in the partitioned form ( B1 ) B = ( ) , ( 0 ) where B1 is N1-by-M. On exit, the leading L-by-M part of this array contains the transformed input matrix ( Bc ) Qc'*B = ( ) , ( 0 ) where Bc is NR-by-M. The matrix ( Bc Ac ) is in the controllability staircase form (1). LDB INTEGER The leading dimension of the array B. LDB >= MAX(1,L). C (input/output) DOUBLE PRECISION array, dimension (LDC,N) On entry, the leading P-by-N part of this array must contain the state/output matrix C. On exit, the leading P-by-N part of this array contains the transformed matrix C*Zc. LDC INTEGER The leading dimension of the array C. LDC >= MAX(1,P). Q (input/output) DOUBLE PRECISION array, dimension (LDQ,L) If COMPQ = 'N': Q is not referenced. If COMPQ = 'I': on entry, Q need not be set; on exit, the leading L-by-L part of this array contains the orthogonal matrix Qc, where Qc' is the product of the transformations applied to A, E, and B on the left. If COMPQ = 'U': on entry, the leading L-by-L part of this array must contain an orthogonal matrix Q; on exit, the leading L-by-L part of this array contains the orthogonal matrix Q*Qc. LDQ INTEGER The leading dimension of the array Q. LDQ >= 1, if COMPQ = 'N'; LDQ >= MAX(1,L), if COMPQ = 'U' or 'I'. Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N) If COMPZ = 'N': Z is not referenced. If COMPZ = 'I': on entry, Z need not be set; on exit, the leading N-by-N part of this array contains the orthogonal matrix Zc, i.e., the product of the transformations applied to A, E, and C on the right. If COMPZ = 'U': on entry, the leading N-by-N part of this array must contain an orthogonal matrix Z; on exit, the leading N-by-N part of this array contains the orthogonal matrix Z*Zc. LDZ INTEGER The leading dimension of the array Z. LDZ >= 1, if COMPZ = 'N'; LDZ >= MAX(1,N), if COMPZ = 'U' or 'I'. NR (output) INTEGER The order of the reduced matrices Ac and Ec, and the number of rows of the reduced matrix Bc; also the order of the controllable part of the pair (B, A-lambda*E). NRBLCK (output) INTEGER _ The number k, of full row rank blocks Ai,i in the staircase form of the pencil (Bc Ac-lambda*Ec) (see (1) and (2)). RTAU (output) INTEGER array, dimension (N1) RTAU(i), for i = 1, ..., NRBLCK, is the row dimension of _ the full row rank block Ai,i-1 in the staircase form (1).

TOL DOUBLE PRECISION The tolerance to be used in rank determinations when transforming (A-lambda*E, B). If the user sets TOL > 0, then the given value of TOL is used as a lower bound for reciprocal condition numbers in rank determinations; a (sub)matrix whose estimated condition number is less than 1/TOL is considered to be of full rank. If the user sets TOL <= 0, then an implicitly computed, default tolerance, defined by TOLDEF = L*N*EPS, is used instead, where EPS is the machine precision (see LAPACK Library routine DLAMCH). TOL < 1.

IWORK INTEGER array, dimension (M) 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 >= MAX(1,N,L,2*(M+N1-1)) For good performance, LDWORK should be generally 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.

The subroutine is based on the reduction algorithm of [1]. If suitable, it uses block algorithms for the reduction of the matrix E and for the corresponding updates of the matrices A, B, and Q. Moreover, for large systems, the row transformations are applied on panels of columns of the matrices A, B, and E.

[1] Varga, A. Computation of Irreducible Generalized State-Space Realizations. Kybernetika, vol. 26, pp. 89-106, 1990.

The algorithm is numerically backward stable and requires 0( N*N1**2 ) floating point operations.

If INFO > 0 on entry, that value is used as block size for the block algorithms. Otherwise, the block size is chosen internally.

**Program Text**

None

None

None