**Purpose**

To find a reduced (controllable, observable, or irreducible) descriptor representation (Ar-lambda*Er,Br,Cr) for an original descriptor representation (A-lambda*E,B,C). The pencil Ar-lambda*Er is in an upper block Hessenberg form, with either Ar or Er upper triangular.

SUBROUTINE TG01JY( JOB, SYSTYP, EQUIL, CKSING, RESTOR, N, M, P, A, $ LDA, E, LDE, B, LDB, C, LDC, NR, INFRED, TOL, $ IWORK, DWORK, LDWORK, INFO ) C .. Scalar Arguments .. CHARACTER CKSING, EQUIL, JOB, RESTOR, SYSTYP INTEGER INFO, LDA, LDB, LDC, LDE, LDWORK, M, N, NR, P C .. Array Arguments .. INTEGER INFRED(*), IWORK(*) DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), DWORK(*), $ E(LDE,*), TOL(*)

**Mode Parameters**

JOB CHARACTER*1 Indicates whether the user wishes to remove the uncontrollable and/or unobservable parts as follows: = 'I': Remove both the uncontrollable and unobservable parts to get an irreducible descriptor representation; = 'C': Remove the uncontrollable part only to get a controllable descriptor representation; = 'O': Remove the unobservable part only to get an observable descriptor representation. SYSTYP CHARACTER*1 Indicates the type of descriptor system algorithm to be applied according to the assumed transfer-function matrix as follows: = 'R': Rational transfer-function matrix; = 'S': Proper (standard) transfer-function matrix; = 'P': Polynomial transfer-function matrix. EQUIL CHARACTER*1 Specifies whether the user wishes to preliminarily scale the system (A-lambda*E,B,C) as follows: = 'S': Perform scaling; = 'N': Do not perform scaling. CKSING CHARACTER*1 Specifies whether the user wishes to check if the pencil (A-lambda*E) is singular as follows: = 'C': Check singularity; = 'N': Do not check singularity. If the pencil is singular, the reduced system computed for CKSING = 'N' can be wrong. RESTOR CHARACTER*1 Specifies whether the user wishes to save the system matrices before each phase and restore them if no order reduction took place as follows: = 'R': Save and restore; = 'N': Do not save the matrices.

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 original state matrix A. On exit, the leading NR-by-NR part of this array contains the reduced order state matrix Ar of an irreducible, controllable, or observable realization for the original system, depending on the value of JOB, JOB = 'I', JOB = 'C', or JOB = 'O', respectively. The matrix Ar is upper triangular if SYSTYP = 'P'. If SYSTYP = 'S' and JOB = 'C', the matrix [Br Ar] is in a controllable staircase form (see SLICOT Library routine TG01HD). If SYSTYP = 'S' and JOB = 'I' or 'O', the matrix ( Ar ) ( Cr ) is in an observable staircase form (see TG01HD). The resulting Ar has INFRED(5) nonzero sub-diagonals. The block structure of staircase forms is contained in the leading INFRED(7) elements of IWORK. 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 original descriptor matrix E. On exit, the leading NR-by-NR part of this array contains the reduced order descriptor matrix Er of an irreducible, controllable, or observable realization for the original system, depending on the value of JOB, JOB = 'I', JOB = 'C', or JOB = 'O', respectively. The resulting Er has INFRED(6) nonzero sub-diagonals. If at least for one k = 1,...,4, INFRED(k) >= 0, then the resulting Er is structured being either upper triangular or block Hessenberg, in accordance to the last performed order reduction phase (see METHOD). The block structure of staircase forms is contained in the leading INFRED(7) elements of IWORK. LDE INTEGER The leading dimension of array E. LDE >= MAX(1,N). B (input/output) DOUBLE PRECISION array, dimension (LDB,M), if JOB = 'C', or (LDB,MAX(M,P)), otherwise. On entry, the leading N-by-M part of this array must contain the original input matrix B; if JOB = 'I', or JOB = 'O', the remainder of the leading N-by-MAX(M,P) part is used as internal workspace. On exit, the leading NR-by-M part of this array contains the reduced input matrix Br of an irreducible, controllable, or observable realization for the original system, depending on the value of JOB, JOB = 'I', JOB = 'C', or JOB = 'O', respectively. If JOB = 'C', only the first IWORK(1) rows of B are nonzero. 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 original output matrix C; if JOB = 'I', or JOB = 'O', the remainder of the leading MAX(M,P)-by-N part is used as internal workspace. On exit, the leading P-by-NR part of this array contains the transformed state/output matrix Cr of an irreducible, controllable, or observable realization for the original system, depending on the value of JOB, JOB = 'I', JOB = 'C', or JOB = 'O', respectively. If JOB = 'I', or JOB = 'O', only the last IWORK(1) columns (in the first NR columns) of C are nonzero. LDC INTEGER The leading dimension of array C. LDC >= MAX(1,M,P) if N > 0. LDC >= 1 if N = 0. NR (output) INTEGER The order of the reduced descriptor representation (Ar-lambda*Er,Br,Cr) of an irreducible, controllable, or observable realization for the original system, depending on JOB = 'I', JOB = 'C', or JOB = 'O', respectively. INFRED (output) INTEGER array, dimension 7 This array contains information on performed reduction and on structure of resulting system matrices as follows: INFRED(k) >= 0 (k = 1, 2, 3, or 4) if Phase k of reduction (see METHOD) has been performed. In this case, INFRED(k) is the achieved order reduction in Phase k. INFRED(k) < 0 (k = 1, 2, 3, or 4) if Phase k was not performed. INFRED(5) - the number of nonzero sub-diagonals of A. INFRED(6) - the number of nonzero sub-diagonals of E. INFRED(7) - the number of blocks in the resulting staircase form at last performed reduction phase. The block dimensions are contained in the first INFRED(7) elements of IWORK.

TOL DOUBLE PRECISION array, dimension 3 TOL(1) is the tolerance to be used in rank determinations when transforming (A-lambda*E,B,C). If the user sets TOL(1) > 0, then the given value of TOL(1) 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(1) is considered to be of full rank. If the user sets TOL(1) <= 0, then an implicitly computed, default tolerance, defined by TOLDEF1 = N*N*EPS, is used instead, where EPS is the machine precision (see LAPACK Library routine DLAMCH). TOL(1) < 1. TOL(2) is the tolerance to be used for checking pencil singularity when CKSING = 'C', or singularity of the matrices A and E when CKSING = 'N'. If the user sets TOL(2) > 0, then the given value of TOL(2) is used. If the user sets TOL(2) <= 0, then an implicitly computed, default tolerance, defined by TOLDEF2 = 10*EPS, is used instead. TOL(2) < 1. TOL(3) is the threshold value for magnitude of the matrix elements, if EQUIL = 'S': elements with magnitude less than or equal to TOL(3) are ignored for scaling. If the user sets TOL(3) >= 0, then the given value of TOL(3) is used. If the user sets TOL(3) < 0, then an implicitly computed, default threshold, defined by THRESH = c*EPS, where c = MAX(norm_1(A,E,B,C)) is used instead. TOL(3) = 0 is not always a good choice. TOL(3) < 1. TOL(3) is not used if EQUIL = 'N'.

IWORK INTEGER array, dimension (2*N+MAX(M,P)) On exit, if INFO = 0, the leading INFRED(7) elements of IWORK contain the orders of the diagonal blocks of Ar-lambda*Er. 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,x,y,8*N), if EQUIL = 'S', LDWORK >= MAX(1,x,y), if EQUIL = 'N', where x = MAX(2*(z+MAX(M,P)+N-1),N*N+4*N), if RESTOR = 'R' x = MAX( 2*(MAX(M,P)+N-1),N*N+4*N), if RESTOR = 'N' y = 2*N*N+10*N+MAX(N,23), if CKSING = 'C', y = 0, if CKSING = 'N', z = 2*N*N+N*M+N*P, if JOB = 'I', z = 0, if JOB <> 'I'. For good performance, LDWORK should be generally larger. If RESTOR = 'R', or LDWORK >= MAX(1,2*N*N+N*M+N*P+2*(MAX(M,P)+N-1), more accurate results are to be expected by considering only those reductions phases (see METHOD), where effective order reduction occurs. This is achieved by saving the system matrices before each phase and restoring them if no order reduction took place. Actually, if JOB = 'I' and RESTOR = 'N', then the saved matrices are those obtained after orthogonally triangularizing the matrix A (if SYSTYP = 'R' or 'P'), or the matrix E (if SYSTYP = 'R' or 'S'). 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 workspace includes the extra space for improving the accuracy.

INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value; = 1: the given pencil A - lambda*E is numerically singular and the reduced system is not computed. This error can be returned only if CKSING = 'C'.

The subroutine is based on the reduction algorithms of [1], but with a different ordering of the phases. The order reduction is performed in 4 phases: Phase 1: Eliminate all infinite and finite nonzero uncontrollable eigenvalues. The resulting matrix ( Br Er ) is in a controllable staircase form (see TG01HD), and Ar is upper triangular. This phase is performed if JOB = 'I' or 'C' and SYSTYP = 'R' or 'P'. Phase 2: Eliminate all infinite and finite nonzero unobservable eigenvalues. The resulting matrix ( Er ) is in an ( Cr ) observable staircase form (see SLICOT Library routine TG01ID), and Ar is upper triangular. This phase is performed if JOB = 'I' or 'O' and SYSTYP = 'R' or 'P'. Phase 3: Eliminate all finite uncontrollable eigenvalues. The resulting matrix ( Br Ar ) is in a controllable staircase form (see TG01HD), and Er is upper triangular. This phase is performed if JOB = 'I' or 'C' and SYSTYP = 'R' or 'S'. Phase 4: Eliminate all finite unobservable eigenvalues. The resulting matrix ( Ar ) is in an observable ( Cr ) staircase form (see TG01ID), and Er is upper triangular. This phase is performed if JOB = 'I' or 'O' and SYSTYP = 'R' or 'S'. The routine checks the singularity of the matrices A and/or E (depending on JOB and SYSTYP) and skips the unnecessary phases. See FURTHER COMMENTS.

[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 pencil A-lambda*E has no zero eigenvalues, then an irreducible realization is computed skipping Phases 3 and 4 (equivalent to setting: JOB = 'I' and SYSTYP = 'P').

**Program Text**

* TG01JY EXAMPLE PROGRAM TEXT * Copyright (c) 2002-2017 NICONET e.V. * * .. Parameters .. INTEGER NIN, NOUT PARAMETER ( NIN = 5, NOUT = 6 ) INTEGER NMAX, MMAX, PMAX PARAMETER ( NMAX = 20, MMAX = 20, PMAX = 20 ) INTEGER LDA, LDB, LDC, LDE PARAMETER ( LDA = NMAX, LDB = NMAX, LDC = PMAX, $ LDE = NMAX ) INTEGER LDWORK, LIWORK PARAMETER ( LDWORK = 2*NMAX*NMAX + $ MAX( 2*( NMAX*( NMAX + MMAX + PMAX ) + $ MAX( MMAX, PMAX ) + NMAX - 1 ), $ 10*NMAX + MAX( NMAX, 23 ) ), $ LIWORK = 2*NMAX + MAX( MMAX, PMAX ) ) * .. Local Scalars .. CHARACTER CKSING, EQUIL, JOB, RESTOR, SYSTYP INTEGER I, INFO, J, M, N, NR, P * .. Local Arrays .. INTEGER INFRED(7), IWORK(LIWORK) DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,NMAX), $ DWORK(LDWORK), E(LDE,NMAX), TOL(3) * .. External Subroutines .. EXTERNAL TG01JY * .. 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(1), TOL(2), TOL(3), JOB, $ SYSTYP, EQUIL, CKSING, RESTOR 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 irreducible descriptor system (Ar-lambda Er,Br,Cr). CALL TG01JY( JOB, SYSTYP, EQUIL, CKSING, RESTOR, N, M, P, $ A, LDA, E, LDE, B, LDB, C, LDC, NR, INFRED, $ TOL, IWORK, DWORK, LDWORK, INFO ) * IF ( INFO.NE.0 ) THEN WRITE ( NOUT, FMT = 99998 ) INFO ELSE WRITE ( NOUT, FMT = 99994 ) NR WRITE ( NOUT, FMT = 99991 ) DO 10 I = 1, 4 IF( INFRED(I).GE.0 ) $ WRITE ( NOUT, FMT = 99990 ) I, INFRED(I) 10 CONTINUE WRITE ( NOUT, FMT = 99997 ) DO 20 I = 1, NR WRITE ( NOUT, FMT = 99995 ) ( A(I,J), J = 1,NR ) 20 CONTINUE WRITE ( NOUT, FMT = 99996 ) DO 30 I = 1, NR WRITE ( NOUT, FMT = 99995 ) ( E(I,J), J = 1,NR ) 30 CONTINUE WRITE ( NOUT, FMT = 99993 ) DO 40 I = 1, NR WRITE ( NOUT, FMT = 99995 ) ( B(I,J), J = 1,M ) 40 CONTINUE WRITE ( NOUT, FMT = 99992 ) DO 50 I = 1, P WRITE ( NOUT, FMT = 99995 ) ( C(I,J), J = 1,NR ) 50 CONTINUE END IF END IF END IF END IF STOP * 99999 FORMAT (' TG01JY EXAMPLE PROGRAM RESULTS',/1X) 99998 FORMAT (' INFO on exit from TG01JY = ',I2) 99997 FORMAT (/' The reduced state dynamics matrix Ar is ') 99996 FORMAT (/' The reduced descriptor matrix Er is ') 99995 FORMAT (20(1X,F8.4)) 99994 FORMAT (' Order of reduced system =', I5 ) 99993 FORMAT (/' The reduced input/state matrix Br is ') 99992 FORMAT (/' The reduced state/output matrix Cr is ') 99991 FORMAT (/' Achieved order reductions in different phases') 99990 FORMAT (' Phase',I2,':', I3, ' elliminated eigenvalue(s)' ) 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

TG01JY EXAMPLE PROGRAM DATA 9 2 2 0.0 0.0 0.0 I R N N N -2 -3 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 -2 -3 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 1 0 0 -1 0 0 0 0 -1 0 0 0 0 1 0 1 -3 0 1 0 2 0 0 1 1 3 0 1 0 0 1

TG01JY EXAMPLE PROGRAM RESULTS Order of reduced system = 7 Achieved order reductions in different phases Phase 1: 0 elliminated eigenvalue(s) Phase 2: 2 elliminated eigenvalue(s) The reduced state dynamics matrix Ar is 1.0000 -0.0393 -0.0980 0.1066 -0.0781 0.2330 -0.0777 0.0000 1.0312 0.2717 -0.2609 0.1533 -0.6758 0.3553 0.0000 0.0000 1.3887 -0.6699 0.4281 -1.6389 0.7615 0.0000 0.0000 0.0000 1.2147 -0.2423 0.9792 -0.4788 0.0000 0.0000 0.0000 0.0000 1.0545 -0.5035 0.2788 0.0000 0.0000 0.0000 0.0000 0.0000 1.6355 -0.4323 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 The reduced descriptor matrix Er is 0.4100 0.2590 0.5080 0.3109 -0.0705 -0.1429 0.1477 -0.7629 -0.3464 0.0992 0.3007 -0.0619 -0.2483 0.0152 0.1120 -0.2124 -0.4184 0.1288 -0.0569 0.4213 0.6182 0.0000 0.1122 -0.0039 -0.2771 0.0758 -0.0975 -0.3923 0.0000 0.0000 0.3708 0.4290 -0.1006 -0.1402 0.2699 0.0000 0.0000 0.0000 0.0000 0.9458 -0.2211 0.2378 0.0000 0.0000 0.0000 0.5711 0.2648 0.5948 -0.5000 The reduced input/state matrix Br is 0.5597 -0.2363 0.4843 0.0498 0.4727 0.1491 -0.1802 -1.1574 -0.5995 -0.1556 -0.1729 -0.3999 0.0000 0.2500 The reduced state/output matrix Cr is 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 4.0000 0.0000 0.0000 0.0000 0.0000 0.0000 3.1524 -1.7500