**Purpose**

To compute for the descriptor system (A-lambda E,B,C) the orthogonal transformation matrices Q and Z such that the transformed system (Q'*A*Z-lambda Q'*E*Z, Q'*B, C*Z) is in an SVD (singular value decomposition) coordinate form with the system matrices Q'*A*Z and Q'*E*Z in the form ( A11 A12 ) ( Er 0 ) Q'*A*Z = ( ) , Q'*E*Z = ( ) , ( A21 A22 ) ( 0 0 ) where Er is an invertible diagonal matrix having on the diagonal the decreasingly ordered nonzero singular values of E. Optionally, the A22 matrix can be further reduced to the SVD form ( Ar 0 ) A22 = ( ) , ( 0 0 ) where Ar is an invertible diagonal matrix having on the diagonal the decreasingly ordered nonzero singular values of A22. The left and/or right orthogonal transformations performed to reduce E and A22 are accumulated.

SUBROUTINE TG01ED( JOBA, L, N, M, P, A, LDA, E, LDE, B, LDB, $ C, LDC, Q, LDQ, Z, LDZ, RANKE, RNKA22, TOL, $ DWORK, LDWORK, INFO ) C .. Scalar Arguments .. CHARACTER JOBA INTEGER INFO, L, LDA, LDB, LDC, LDE, LDQ, LDWORK, $ LDZ, M, N, P, RNKA22, RANKE DOUBLE PRECISION TOL C .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ), $ DWORK( * ), E( LDE, * ), Q( LDQ, * ), $ Z( LDZ, * )

**Mode Parameters**

JOBA CHARACTER*1 = 'N': do not reduce A22; = 'R': reduce A22 to an SVD form.

L (input) INTEGER The number of rows of matrices A, B, and E. L >= 0. N (input) INTEGER The number of columns of matrices A, E, and C. N >= 0. M (input) INTEGER The number of columns of matrix B. M >= 0. P (input) INTEGER The number of rows of matrix C. 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, the leading L-by-N part of this array contains the transformed matrix Q'*A*Z. If JOBA = 'R', this matrix is in the form ( A11 * * ) Q'*A*Z = ( * Ar 0 ) , ( * 0 0 ) where A11 is a RANKE-by-RANKE matrix and Ar is a RNKA22-by-RNKA22 invertible diagonal matrix, with decresingly ordered positive diagonal elements. LDA INTEGER The leading dimension of 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, the leading L-by-N part of this array contains the transformed matrix Q'*E*Z. ( Er 0 ) Q'*E*Z = ( ) , ( 0 0 ) where Er is a RANKE-by-RANKE invertible diagonal matrix having on the diagonal the decreasingly ordered positive singular values of E. LDE INTEGER The leading dimension of 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/state matrix B. On exit, the leading L-by-M part of this array contains the transformed matrix Q'*B. LDB INTEGER The leading dimension of array B. LDB >= MAX(1,L) if M > 0 or LDB >= 1 if M = 0. 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 (output) DOUBLE PRECISION array, dimension (LDQ,L) The leading L-by-L part of this array contains the orthogonal matrix Q, which is the accumulated product of transformations applied to A, E, and B on the left. LDQ INTEGER The leading dimension of array Q. LDQ >= MAX(1,L). 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 accumulated product of transformations applied to A, E, and C on the right. LDZ INTEGER The leading dimension of array Z. LDZ >= MAX(1,N). RANKE (output) INTEGER The effective rank of matrix E, and thus also the order of the invertible diagonal submatrix Er. RANKE is computed as the number of singular values of E greater than TOL*SVEMAX, where SVEMAX is the maximum singular value of E. RNKA22 (output) INTEGER If JOBA = 'R', then RNKA22 is the effective rank of matrix A22, and thus also the order of the invertible diagonal submatrix Ar. RNKA22 is computed as the number of singular values of A22 greater than TOL*SVAMAX, where SVAMAX is an estimate of the maximum singular value of A. If JOBA = 'N', then RNKA22 is not referenced.

TOL DOUBLE PRECISION The tolerance to be used in determining the rank of E and of A22. If TOL > 0, then singular values less than TOL*SVMAX are treated as zero, where SVMAX is the maximum singular value of E or an estimate of it for A and E. If TOL <= 0, the default tolerance TOLDEF = EPS*L*N is used instead, where EPS is the machine precision (see LAPACK Library routine DLAMCH). TOL < 1.

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,MIN(L,N) + MAX(3*MIN(L,N)+MAX(L,N), 5*MIN(L,N), M, P)).

INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value; > 0: the QR algorithm has failed to converge when computing singular value decomposition. In this case INFO specifies how many superdiagonals did not converge. This failure is not likely to occur.

The routine computes the singular value decomposition (SVD) of E, in the form ( Er 0 ) E = Q * ( ) * Z' ( 0 0 ) and finds the largest RANKE-by-RANKE leading diagonal submatrix Er whose condition number is less than 1/TOL. RANKE defines thus the effective rank of matrix E. If JOBA = 'R' the same reduction is performed on A22 in the partitioned matrix ( A11 A12 ) Q'*A*Z = ( ) , ( A21 A22 ) to obtain it in the form ( Ar 0 ) A22 = ( ) , ( 0 0 ) with Ar an invertible diagonal matrix. The accumulated transformations are also applied to the rest of matrices B <- Q' * B, C <- C * Z.

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

None

**Program Text**

* TG01ED 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, MIN( LMAX, NMAX ) + $ MAX( MMAX, PMAX, 3*MIN( LMAX, NMAX ) + $ MAX( LMAX, NMAX ), $ 5*MIN( LMAX, NMAX ) ) ) ) * .. Local Scalars .. CHARACTER*1 JOBA INTEGER I, INFO, J, L, M, N, P, RANKE, RNKA22 DOUBLE PRECISION TOL * .. Local Arrays .. 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 TG01ED * .. 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, TOL JOBA = 'R' 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 ) * Find the transformed descriptor system * (A-lambda E,B,C). CALL TG01ED( JOBA, L, N, M, P, A, LDA, E, LDE, B, LDB, $ C, LDC, Q, LDQ, Z, LDZ, RANKE, RNKA22, $ TOL, DWORK, LDWORK, INFO ) * IF ( INFO.NE.0 ) THEN WRITE ( NOUT, FMT = 99998 ) INFO ELSE WRITE ( NOUT, FMT = 99994 ) RANKE, RNKA22 WRITE ( NOUT, FMT = 99997 ) DO 10 I = 1, L WRITE ( NOUT, FMT = 99995 ) ( A(I,J), J = 1,N ) 10 CONTINUE WRITE ( NOUT, FMT = 99996 ) DO 20 I = 1, L WRITE ( NOUT, FMT = 99995 ) ( E(I,J), J = 1,N ) 20 CONTINUE WRITE ( NOUT, FMT = 99993 ) DO 30 I = 1, L 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, L WRITE ( NOUT, FMT = 99995 ) ( Q(I,J), J = 1,L ) 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 END IF STOP * 99999 FORMAT (' TG01ED EXAMPLE PROGRAM RESULTS',/1X) 99998 FORMAT (' INFO on exit from TG01ED = ',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 (' Rank of matrix E =', I5/ $ ' Rank of matrix A22 =', 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) END

TG01ED EXAMPLE PROGRAM DATA 4 4 2 2 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

TG01ED EXAMPLE PROGRAM RESULTS Rank of matrix E = 3 Rank of matrix A22 = 1 The transformed state dynamics matrix Q'*A*Z is 2.1882 -0.8664 -3.5097 -2.1353 -0.4569 -0.2146 1.9802 0.3531 -0.5717 -0.5245 -0.4591 0.4696 -0.4766 -0.5846 2.1414 0.3086 The transformed descriptor matrix Q'*E*Z is 11.8494 0.0000 0.0000 0.0000 0.0000 2.1302 0.0000 0.0000 0.0000 0.0000 1.0270 0.0000 0.0000 0.0000 0.0000 0.0000 The transformed input/state matrix Q'*B is -0.2396 -1.0668 -0.2656 -0.8393 -0.7657 -0.1213 1.1339 0.3780 The transformed state/output matrix C*Z is -0.2499 -1.0573 0.3912 -0.8165 -0.5225 1.3958 0.8825 0.0000 The left transformation matrix Q is -0.1534 0.5377 -0.6049 0.5669 -0.0872 0.2536 0.7789 0.5669 -0.9805 -0.0360 0.0395 -0.1890 -0.0863 -0.8033 -0.1608 0.5669 The right transformation matrix Z is -0.2612 0.2017 -0.4737 0.8165 -0.7780 0.4718 -0.0738 -0.4082 -0.5111 -0.8556 -0.0826 0.0000 -0.2556 0.0684 0.8737 0.4082