**Purpose**

To find a reduced descriptor representation (Ar-lambda*Er,Br,Cr) without non-dynamic modes for a descriptor representation (A-lambda*E,B,C). Optionally, the reduced descriptor system can be put into a standard form with the leading diagonal block of Er identity.

SUBROUTINE TG01GD( JOBS, L, N, M, P, A, LDA, E, LDE, B, LDB, $ C, LDC, D, LDD, LR, NR, RANKE, INFRED, TOL, $ IWORK, DWORK, LDWORK, INFO ) C .. Scalar Arguments .. CHARACTER JOBS INTEGER INFO, INFRED, L, LDA, LDB, LDC, LDD, LDE, $ LDWORK, LR, M, N, NR, P, RANKE DOUBLE PRECISION TOL C .. Array Arguments .. INTEGER IWORK(*) DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*), $ DWORK(*), E(LDE,*)

**Mode Parameters**

JOBS CHARACTER*1 Indicates whether the user wishes to transform the leading diagonal block of Er to an identity matrix, as follows: = 'S': make Er with leading diagonal identity; = 'D': keep Er unreduced or upper triangular.

L (input) INTEGER The number of rows of the matrices A, E, and B; also the number of differential equations. L >= 0. N (input) INTEGER The number of columns of the matrices A, E, and C; also the dimension of descriptor state vector. N >= 0. M (input) INTEGER The number of columns of the matrix B; also the dimension of the input vector. M >= 0. P (input) INTEGER The number of rows of the matrix C. also the dimension of the output vector. P >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the leading L-by-N part of this array must contain the state dynamics matrix A. On exit, if NR < N, the leading LR-by-NR part of this array contains the reduced order state matrix Ar of a descriptor realization without non-dynamic modes. Array A contains the original state dynamics matrix if INFRED < 0. 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 descriptor matrix E. On exit, if INFRED >= 0, the leading LR-by-NR part of this array contains the reduced order descriptor matrix Er of a descriptor realization without non-dynamic modes. In this case, only the leading RANKE-by-RANKE submatrix of Er is nonzero and this submatrix is nonsingular and upper triangular. Array E contains the original descriptor matrix if INFRED < 0. If JOBS = 'S', then the leading RANKE-by-RANKE submatrix results in an identity matrix. 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 input matrix B. On exit, the leading LR-by-M part of this array contains the reduced order input matrix Br of a descriptor realization without non-dynamic modes. Array B contains the original input matrix if INFRED < 0. 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 output matrix C. On exit, the leading P-by-NR part of this array contains the reduced order output matrix Cr of a descriptor realization without non-dynamic modes. Array C contains the original output matrix if INFRED < 0. LDC INTEGER The leading dimension of the array C. LDC >= MAX(1,P). D (input/output) DOUBLE PRECISION array, dimension (LDD,M) On entry, the leading P-by-M part of this array must contain the original feedthrough matrix D. On exit, the leading P-by-M part of this array contains the feedthrough matrix Dr of a reduced descriptor realization without non-dynamic modes. LDD INTEGER The leading dimension of the array D. LDD >= MAX(1,P). LR (output) INTEGER The number of reduced differential equations. NR (output) INTEGER The dimension of the reduced descriptor state vector. RANKE (output) INTEGER The estimated rank of the matrix E. INFRED (output) INTEGER This parameter contains information on performed reduction and on structure of resulting system matrices, as follows: INFRED >= 0 the reduced system is in an SVD-like coordinate form with Er upper triangular; INFRED is the achieved order reduction. INFRED < 0 no reduction achieved and the original system has been restored.

TOL DOUBLE PRECISION The tolerance to be used in rank determinations when transforming (A-lambda*E). 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 (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, if MIN(L,N) = 0; otherwise, LDWORK >= MAX( N+P, MIN(L,N)+MAX(3*N-1,M,L) ). If LDWORK >= 2*L*N+L*M+N*P+ MAX( 1, N+P, MIN(L,N)+MAX(3*N-1,M,L) ) then the original matrices are restored if no order reduction is possible. This is achieved by saving system matrices before reduction and restoring them if no order reduction took place. 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. The optimal size does not necessarily include the space needed for saving the original system matrices.

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

The subroutine elliminates the non-dynamics modes in two steps: Step 1: Reduce the system to the SVD-like coordinate form (Q'*A*Z-lambda*Q'*E*Z, Q'*B, C*Z) , where ( A11 A12 A13 ) ( E11 0 0 ) ( B1 ) Q'*A*Z = ( A21 A22 0 ), Q'*E*Z = ( 0 0 0 ), Q'*B = ( B2 ), ( A31 0 0 ) ( 0 0 0 ) ( B3 ) C*Z = ( C1 C2 C3 ), where E11 and A22 are upper triangular invertible matrices. Step 2: Compute the reduced system as (Ar-lambda*Er,Br,Cr,Dr), where ( A11 - A12*inv(A22)*A21, A13 ) ( E11 0 ) Ar = ( ), Er = ( ), ( A31 0 ) ( 0 0 ) ( B1 - A12*inv(A22)*B2 ) Br = ( ), Cr = ( C1 - C2*inv(A22)*A21, C3 ), ( B3 ) Dr = D - C2*inv(A22)*B2. Step 3: If desired (JOBS = 'S'), reduce the descriptor system to the standard form Ar <- diag(inv(Er),I)*Ar; Br <- diag(inv(Er),I)*Br; Er = diag(I,0). If L = N and LR = NR = RANKE, then if Step 3 is performed, the resulting system is a standard state space system.

If L = N, the algorithm requires 0( N**3 ) floating point operations.

None

**Program Text**

* TG01GD 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, LDD, LDE PARAMETER ( LDA = LMAX, LDB = LMAX, LDC = PMAX, $ LDD = PMAX, LDE = LMAX ) INTEGER LDWORK PARAMETER ( LDWORK = MIN( LMAX, NMAX ) + $ MAX( 3*NMAX - 1, MMAX, LMAX ) + $ 2*LMAX*NMAX + LMAX*MMAX + PMAX*NMAX ) * .. Local Scalars .. CHARACTER*1 JOBS INTEGER I, INFO, INFRED, J, L, LR, M, N, NR, P, RANKE DOUBLE PRECISION TOL * .. Local Arrays .. INTEGER IWORK(NMAX) DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,NMAX), $ D(LDD,MMAX), DWORK(LDWORK), E(LDE,NMAX) * .. External Subroutines .. EXTERNAL TG01GD * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. Executable Statements .. * WRITE ( NOUT, FMT = 99999 ) * Skip the heading in the data file and read the data. READ ( NIN, FMT = '()' ) READ ( NIN, FMT = * ) L, N, M, P, JOBS, TOL IF ( L.LT.0 .OR. L.GT.LMAX ) THEN WRITE ( NOUT, FMT = 99989 ) L ELSE 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,L ) READ ( NIN, FMT = * ) ( ( E(I,J), J = 1,N ), I = 1,L ) 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,L ) 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 ) READ ( NIN, FMT = * ) ( ( D(I,J), J = 1,M ), I = 1,P ) * Find the reduced descriptor system * (A-lambda E,B,C,D). CALL TG01GD( JOBS, L, N, M, P, A, LDA, E, LDE, B, LDB, $ C, LDC, D, LDD, LR, NR, RANKE, INFRED, $ TOL, IWORK, DWORK, LDWORK, INFO ) * IF ( INFO.NE.0 ) THEN WRITE ( NOUT, FMT = 99998 ) INFO ELSE WRITE ( NOUT, FMT = 99994 ) RANKE WRITE ( NOUT, FMT = 99997 ) DO 10 I = 1, LR WRITE ( NOUT, FMT = 99995 ) ( A(I,J), J = 1,NR ) 10 CONTINUE WRITE ( NOUT, FMT = 99996 ) DO 20 I = 1, LR WRITE ( NOUT, FMT = 99995 ) ( E(I,J), J = 1,NR ) 20 CONTINUE WRITE ( NOUT, FMT = 99993 ) DO 30 I = 1, LR 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,NR ) 40 CONTINUE WRITE ( NOUT, FMT = 99991 ) DO 50 I = 1, P WRITE ( NOUT, FMT = 99995 ) ( D(I,J), J = 1,M ) 50 CONTINUE END IF END IF END IF END IF END IF STOP * 99999 FORMAT (' TG01GD EXAMPLE PROGRAM RESULTS',/1X) 99998 FORMAT (' INFO on exit from TG01GD = ',I2) 99997 FORMAT (/' The reduced state dynamics matrix is ') 99996 FORMAT (/' The reduced descriptor matrix is ') 99995 FORMAT (20(1X,F8.4)) 99994 FORMAT (' Rank of matrix E =', I5) 99993 FORMAT (/' The reduced input/state matrix is ') 99992 FORMAT (/' The reduced state/output matrix is ') 99991 FORMAT (/' The transformed feedthrough matrix 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) END

TG01GD EXAMPLE PROGRAM DATA 4 4 2 2 D 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 1 0 1 1

TG01GD EXAMPLE PROGRAM RESULTS Rank of matrix E = 3 The reduced state dynamics matrix is 2.5102 -3.8550 -11.4533 -0.0697 0.0212 0.7015 0.3798 -0.1156 -3.8250 The reduced descriptor matrix is 10.1587 5.8230 1.3021 0.0000 -2.4684 -0.1896 0.0000 0.0000 1.0338 The reduced input/state matrix is 7.7100 1.6714 0.7678 1.1070 2.5428 0.6935 The reduced state/output matrix is 0.5477 -2.5000 -6.2610 -1.0954 1.0000 -0.8944 The transformed feedthrough matrix is 4.0000 1.0000 1.0000 1.0000