**Purpose**

To compute orthogonal transformation matrices Q and Z which reduce the N-th order descriptor system (A-lambda*E,B,C) to the form ( Ac * ) ( Ec * ) ( Bc ) Q'*A*Z = ( ) , Q'*E*Z = ( ) , Q'*B = ( ) , ( 0 Anc ) ( 0 Enc ) ( 0 ) C*Z = ( Cc Cnc ) , where the NCONT-th order descriptor system (Ac-lambda*Ec,Bc,Cc) is a finite and/or infinite controllable. The pencil Anc - lambda*Enc is regular of order N-NCONT and contains the uncontrollable finite and/or infinite eigenvalues of the pencil A-lambda*E. For JOBCON = 'C' or 'I', the pencil ( Bc Ec-lambda*Ac ) has full row rank NCONT for all finite lambda and is in a staircase form with _ _ _ _ ( E1,0 E1,1 ... E1,k-1 E1,k ) ( _ _ _ ) ( Bc Ec ) = ( 0 E2,1 ... E2,k-1 E2,k ) , (1) ( ... _ _ ) ( 0 0 ... Ek,k-1 Ek,k ) _ _ _ ( A1,1 ... A1,k-1 A1,k ) ( _ _ ) Ac = ( 0 ... A2,k-1 A2,k ) , (2) ( ... _ ) ( 0 ... 0 Ak,k ) _ where Ei,i-1 is an rtau(i)-by-rtau(i-1) full row rank matrix _ (with rtau(0) = M) and Ai,i is an rtau(i)-by-rtau(i) upper triangular matrix. For JOBCON = 'F', the pencil ( Bc Ac-lambda*Ec ) has full row rank NCONT 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 ) , (3) ( ... _ _ ) ( 0 0 ... Ak,k-1 Ak,k ) _ _ _ ( E1,1 ... E1,k-1 E1,k ) ( _ _ ) Ec = ( 0 ... E2,k-1 E2,k ) , (4) ( ... _ ) ( 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. For JOBCON = 'C', the (N-NCONT)-by-(N-NCONT) regular pencil Anc - lambda*Enc has the form ( Ainc - lambda*Einc * ) Anc - lambda*Enc = ( ) , ( 0 Afnc - lambda*Efnc ) where: 1) the NIUCON-by-NIUCON regular pencil Ainc - lambda*Einc, with Ainc upper triangular and nonsingular, contains the uncontrollable infinite eigenvalues of A - lambda*E; 2) the (N-NCONT-NIUCON)-by-(N-NCONT-NIUCON) regular pencil Afnc - lambda*Efnc, with Efnc upper triangular and nonsingular, contains the uncontrollable finite eigenvalues of A - lambda*E. Note: The significance of the two diagonal blocks can be interchanged by calling the routine with the arguments A and E interchanged. In this case, Ainc - lambda*Einc contains the uncontrollable zero eigenvalues of A - lambda*E, while Afnc - lambda*Efnc contains the uncontrollable nonzero finite and infinite eigenvalues of A - lambda*E. For JOBCON = 'F', the pencil Anc - lambda*Enc has the form Anc - lambda*Enc = Afnc - lambda*Efnc , where the regular pencil Afnc - lambda*Efnc, with Efnc upper triangular and nonsingular, contains the uncontrollable finite eigenvalues of A - lambda*E. For JOBCON = 'I', the pencil Anc - lambda*Enc has the form Anc - lambda*Enc = Ainc - lambda*Einc , where the regular pencil Ainc - lambda*Einc, with Ainc upper triangular and nonsingular, contains the uncontrollable nonzero finite and infinite eigenvalues of A - lambda*E. The left and/or right orthogonal transformations Q and Z performed to reduce the system matrices can be optionally accumulated. The reduced order descriptor system (Ac-lambda*Ec,Bc,Cc) has the same transfer-function matrix as the original system (A-lambda*E,B,C).

SUBROUTINE TG01HD( JOBCON, COMPQ, COMPZ, N, M, P, A, LDA, E, LDE, $ B, LDB, C, LDC, Q, LDQ, Z, LDZ, NCONT, NIUCON, $ NRBLCK, RTAU, TOL, IWORK, DWORK, INFO ) C .. Scalar Arguments .. CHARACTER COMPQ, COMPZ, JOBCON INTEGER INFO, LDA, LDB, LDC, LDE, LDQ, LDZ, $ M, N, NCONT, NIUCON, 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**

JOBCON CHARACTER*1 = 'C': separate both finite and infinite uncontrollable eigenvalues; = 'F': separate only finite uncontrollable eigenvalues: = 'I': separate only nonzero finite and infinite uncontrollable eigenvalues. 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.

N (input) INTEGER The dimension of the descriptor state vector; also the order of square matrices A and E, the number of rows of matrix B, and the number of columns of matrix C. N >= 0. M (input) INTEGER The dimension of descriptor system input vector; also the number of columns of matrix B. M >= 0. P (input) INTEGER The dimension of descriptor system output vector; also the number of rows of matrix C. P >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the leading N-by-N part of this array must contain the N-by-N state matrix A. On exit, the leading N-by-N part of this array contains the transformed state matrix Q'*A*Z, ( Ac * ) Q'*A*Z = ( ) , ( 0 Anc ) where Ac is NCONT-by-NCONT and Anc is (N-NCONT)-by-(N-NCONT). If JOBCON = 'F', the matrix ( Bc Ac ) is in the controllability staircase form (3). If JOBCON = 'C' or 'I', the submatrix Ac is upper triangular. If JOBCON = 'C', the Anc matrix has the form ( Ainc * ) Anc = ( ) , ( 0 Afnc ) where the NIUCON-by-NIUCON matrix Ainc is nonsingular and upper triangular. If JOBCON = 'I', Anc is nonsingular and upper triangular. LDA INTEGER The leading dimension of 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 N-by-N descriptor matrix E. On exit, the leading N-by-N part of this array contains the transformed descriptor matrix Q'*E*Z, ( Ec * ) Q'*E*Z = ( ) , ( 0 Enc ) where Ec is NCONT-by-NCONT and Enc is (N-NCONT)-by-(N-NCONT). If JOBCON = 'C' or 'I', the matrix ( Bc Ec ) is in the controllability staircase form (1). If JOBCON = 'F', the submatrix Ec is upper triangular. If JOBCON = 'C', the Enc matrix has the form ( Einc * ) Enc = ( ) , ( 0 Efnc ) where the NIUCON-by-NIUCON matrix Einc is nilpotent and the (N-NCONT-NIUCON)-by-(N-NCONT-NIUCON) matrix Efnc is nonsingular and upper triangular. If JOBCON = 'F', Enc is nonsingular and upper triangular. LDE INTEGER The leading dimension of 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 N-by-M input matrix B. On exit, the leading N-by-M part of this array contains the transformed input matrix ( Bc ) Q'*B = ( ) , ( 0 ) where Bc is NCONT-by-M. For JOBCON = 'C' or 'I', the matrix ( Bc Ec ) is in the controllability staircase form (1). For JOBCON = 'F', the matrix ( Bc Ac ) is in the controllability staircase form (3). LDB INTEGER The leading dimension of 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 state/output matrix C. On exit, the leading P-by-N part of this array contains the transformed matrix C*Z. LDC INTEGER The leading dimension of array C. LDC >= MAX(1,P). Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N) If COMPQ = 'N': Q is not referenced. If COMPQ = 'I': on entry, Q need not be set; on exit, the leading N-by-N part of this array contains the orthogonal matrix Q, where Q' is the product of transformations which are applied to A, E, and B on the left. If COMPQ = 'U': on entry, the leading N-by-N part of this array must contain an orthogonal matrix Qc; on exit, the leading N-by-N part of this array contains the orthogonal matrix Qc*Q. LDQ INTEGER The leading dimension of array Q. LDQ >= 1, if COMPQ = 'N'; LDQ >= MAX(1,N), 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 Z, which is the product of 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 Zc; on exit, the leading N-by-N part of this array contains the orthogonal matrix Zc*Z. LDZ INTEGER The leading dimension of array Z. LDZ >= 1, if COMPZ = 'N'; LDZ >= MAX(1,N), if COMPZ = 'U' or 'I'. NCONT (output) INTEGER The order of the reduced matrices Ac and Ec, and the number of rows of reduced matrix Bc; also the order of the controllable part of the pair (A-lambda*E,B). NIUCON (output) INTEGER For JOBCON = 'C', the order of the reduced matrices Ainc and Einc; also the number of uncontrollable infinite eigenvalues of the pencil A - lambda*E. For JOBCON = 'F' or 'I', NIUCON has no significance and is set to zero. NRBLCK (output) INTEGER For JOBCON = 'C' or 'I', the number k, of full row rank _ blocks Ei,i in the staircase form of the pencil (Bc Ec-lambda*Ac) (see (1) and (2)). For JOBCON = 'F', the number k, of full row rank blocks _ Ai,i in the staircase form of the pencil (Bc Ac-lambda*Ec) (see (3) and (4)). RTAU (output) INTEGER array, dimension (N) RTAU(i), for i = 1, ..., NRBLCK, is the row dimension of _ _ the full row rank block Ei,i-1 or Ai,i-1 in the staircase form (1) or (3) for JOBCON = 'C' or 'I', or for JOBCON = 'F', respectively.

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 = N*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 (MAX(N,2*M))

INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value.

The subroutine is based on the reduction algorithms of [1].

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

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

If the system matrices A, E and B are badly scaled, it is generally recommendable to scale them with the SLICOT routine TG01AD, before calling TG01HD.

**Program Text**

* TG01HD EXAMPLE PROGRAM TEXT * Copyright (c) 2002-2017 NICONET e.V. * * .. Parameters .. INTEGER NIN, NOUT PARAMETER ( NIN = 5, NOUT = 6 ) INTEGER LMAX, NMAX, MMAX, PMAX PARAMETER ( LMAX = 20, NMAX = 20, MMAX = 20, PMAX = 20 ) INTEGER LDA, LDB, LDC, LDE, LDQ, LDZ PARAMETER ( LDA = LMAX, LDB = LMAX, LDC = PMAX, $ LDE = LMAX, LDQ = LMAX, LDZ = NMAX ) INTEGER LDWORK PARAMETER ( LDWORK = MAX( 1, NMAX, 2*MMAX ) ) * .. Local Scalars .. CHARACTER*1 COMPQ, COMPZ, JOBCO INTEGER I, INFO, J, M, N, NCONT, NIUCON, NRBLCK, P DOUBLE PRECISION TOL * .. Local Arrays .. INTEGER IWORK(MMAX), RTAU(NMAX) DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,NMAX), $ DWORK(LDWORK), E(LDE,NMAX), Q(LDQ,LMAX), $ Z(LDZ,NMAX) * .. External Subroutines .. EXTERNAL TG01HD * .. Intrinsic Functions .. INTRINSIC MAX * .. 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, TOL, JOBCO COMPQ = 'I' COMPZ = 'I' 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 transformed descriptor system (A-lambda E,B,C). CALL TG01HD( JOBCO, COMPQ, COMPZ, N, M, P, A, LDA, $ E, LDE, B, LDB, C, LDC, Q, LDQ, Z, LDZ, $ NCONT, NIUCON, NRBLCK, RTAU, TOL, IWORK, $ DWORK, INFO ) * IF ( INFO.NE.0 ) THEN WRITE ( NOUT, FMT = 99998 ) INFO ELSE WRITE ( NOUT, FMT = 99994 ) NCONT, NIUCON WRITE ( NOUT, FMT = 99985 ) WRITE ( NOUT, FMT = 99984 ) ( RTAU(I), I = 1,NRBLCK ) 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 END IF END IF END IF END IF STOP * 99999 FORMAT (' TG01HD EXAMPLE PROGRAM RESULTS',/1X) 99998 FORMAT (' INFO on exit from TG01HD = ',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 (' Dimension of controllable part =', I5/ $ ' Number of uncontrollable infinite eigenvalues =', I5) 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 (/' L is out of range.',/' L = ',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 staircase form row dimensions are ' ) 99984 FORMAT (10I5) END

TG01HD EXAMPLE PROGRAM DATA 7 3 2 0.0 C 2 0 2 0 -1 3 1 0 1 0 0 1 0 0 0 0 0 1 0 0 1 0 0 2 0 -1 3 1 0 0 0 1 0 0 1 0 1 0 0 1 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 1 0 -1 0 0 1 3 0 2 0 0 0 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 3 1 0 0 1 0 0 1 0 -1 1 0 -1 1 0

TG01HD EXAMPLE PROGRAM RESULTS Dimension of controllable part = 3 Number of uncontrollable infinite eigenvalues = 1 The staircase form row dimensions are 2 1 The transformed state dynamics matrix Q'*A*Z is 0.0000 0.0000 0.0000 0.0000 -1.2627 0.4334 0.4666 0.0000 2.0000 0.0000 -3.7417 -0.8520 0.2924 -0.4342 0.0000 0.0000 1.7862 0.3780 -0.2651 -0.7723 0.0000 0.0000 0.0000 0.0000 3.7417 0.8520 -0.2924 0.4342 0.0000 0.0000 0.0000 0.0000 -1.5540 0.5334 0.5742 0.0000 0.0000 0.0000 0.0000 -0.6533 0.2242 0.2414 0.0000 0.0000 0.0000 0.0000 -0.5892 0.2022 0.2177 The transformed descriptor matrix Q'*E*Z is -1.8325 1.0000 2.3752 0.0000 -0.8214 0.2819 1.8016 0.4887 0.0000 0.3770 -0.5345 0.1874 0.5461 0.0000 -0.1728 0.0000 -0.1333 -1.1339 0.1325 0.3861 0.0000 0.0000 0.0000 0.0000 0.0000 0.8520 -0.2924 0.4342 0.0000 0.0000 0.0000 0.0000 -1.0260 -0.1496 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.1937 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 The transformed input/state matrix Q'*B is 1.0000 2.0000 3.0000 2.0000 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 The transformed state/output matrix C*Z is 0.0000 1.0000 0.0000 0.0000 -1.2627 0.4334 0.4666 0.3665 0.0000 -0.9803 -1.6036 0.1874 0.5461 0.0000 The left transformation matrix Q is 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.7071 0.0000 0.2740 -0.6519 0.0000 0.0000 0.0000 0.0000 0.0000 0.8304 0.3491 -0.4342 0.0000 0.0000 0.0000 -1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.4003 0.1683 0.9008 0.0000 0.0000 0.7071 0.0000 -0.2740 0.6519 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 The right transformation matrix Z is 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.6108 0.0000 0.7917 0.0000 0.0000 0.0000 0.0000 0.4887 0.0000 0.3770 -0.5345 0.1874 0.5461 0.0000 0.0000 0.0000 0.0000 0.0000 -0.4107 0.1410 0.9008 0.6108 0.0000 0.4713 0.2673 -0.1874 -0.5461 0.0000 -0.1222 0.0000 -0.0943 -0.8018 -0.1874 -0.5461 0.0000 0.0000 0.0000 0.0000 0.0000 -0.8520 0.2924 -0.4342