**Purpose**

To find the transfer matrix T(s) of a given state-space representation (A,B,C,D). T(s) is expressed as either row or column polynomial vectors over monic least common denominator polynomials.

SUBROUTINE TB04AD( ROWCOL, N, M, P, A, LDA, B, LDB, C, LDC, D, $ LDD, NR, INDEX, DCOEFF, LDDCOE, UCOEFF, LDUCO1, $ LDUCO2, TOL1, TOL2, IWORK, DWORK, LDWORK, $ INFO ) C .. Scalar Arguments .. CHARACTER ROWCOL INTEGER INFO, LDA, LDB, LDC, LDD, LDDCOE, LDUCO1, $ LDUCO2, LDWORK, M, N, NR, P DOUBLE PRECISION TOL1, TOL2 C .. Array Arguments .. INTEGER INDEX(*), IWORK(*) DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*), $ DCOEFF(LDDCOE,*), DWORK(*), $ UCOEFF(LDUCO1,LDUCO2,*)

**Mode Parameters**

ROWCOL CHARACTER*1 Indicates whether the transfer matrix T(s) is required as rows or columns over common denominators as follows: = 'R': T(s) is required as rows over common denominators; = 'C': T(s) is required as columns over common denominators.

N (input) INTEGER The order of the state-space representation, i.e. the order of the original state dynamics matrix A. N >= 0. M (input) INTEGER The number of system inputs. M >= 0. P (input) INTEGER The number of system outputs. 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 dynamics matrix A. On exit, the leading NR-by-NR part of this array contains the upper block Hessenberg state dynamics matrix A of a transformed representation for the original system: this is completely controllable if ROWCOL = 'R', or completely observable if ROWCOL = 'C'. LDA INTEGER The leading dimension of array A. LDA >= MAX(1,N). B (input/output) DOUBLE PRECISION array, dimension (LDB,M), if ROWCOL = 'R', and (LDB,MAX(M,P)) if ROWCOL = 'C'. On entry, the leading N-by-M part of this array must contain the original input/state matrix B; if ROWCOL = 'C', 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 transformed input/state matrix B. 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 state/output matrix C; if ROWCOL = 'C', 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 C. LDC INTEGER The leading dimension of array C. LDC >= MAX(1,P) if ROWCOL = 'R'; LDC >= MAX(1,M,P) if ROWCOL = 'C'. D (input) DOUBLE PRECISION array, dimension (LDD,M), if ROWCOL = 'R', and (LDD,MAX(M,P)) if ROWCOL = 'C'. The leading P-by-M part of this array must contain the original direct transmission matrix D; if ROWCOL = 'C', this array is modified internally, but restored on exit, and the remainder of the leading MAX(M,P)-by-MAX(M,P) part is used as internal workspace. LDD INTEGER The leading dimension of array D. LDD >= MAX(1,P) if ROWCOL = 'R'; LDD >= MAX(1,M,P) if ROWCOL = 'C'. NR (output) INTEGER The order of the transformed state-space representation. INDEX (output) INTEGER array, dimension (porm), where porm = P, if ROWCOL = 'R', and porm = M, if ROWCOL = 'C'. The degrees of the denominator polynomials. DCOEFF (output) DOUBLE PRECISION array, dimension (LDDCOE,N+1) The leading porm-by-kdcoef part of this array contains the coefficients of each denominator polynomial, where kdcoef = MAX(INDEX(I)) + 1. DCOEFF(I,K) is the coefficient in s**(INDEX(I)-K+1) of the I-th denominator polynomial, where K = 1,2,...,kdcoef. LDDCOE INTEGER The leading dimension of array DCOEFF. LDDCOE >= MAX(1,P) if ROWCOL = 'R'; LDDCOE >= MAX(1,M) if ROWCOL = 'C'. UCOEFF (output) DOUBLE PRECISION array, dimension (LDUCO1,LDUCO2,N+1) If ROWCOL = 'R' then porp = M, otherwise porp = P. The leading porm-by-porp-by-kdcoef part of this array contains the coefficients of the numerator matrix U(s). UCOEFF(I,J,K) is the coefficient in s**(INDEX(iorj)-K+1) of polynomial (I,J) of U(s), where K = 1,2,...,kdcoef; if ROWCOL = 'R' then iorj = I, otherwise iorj = J. Thus for ROWCOL = 'R', U(s) = diag(s**INDEX(I))*(UCOEFF(.,.,1)+UCOEFF(.,.,2)/s+...). LDUCO1 INTEGER The leading dimension of array UCOEFF. LDUCO1 >= MAX(1,P) if ROWCOL = 'R'; LDUCO1 >= MAX(1,M) if ROWCOL = 'C'. LDUCO2 INTEGER The second dimension of array UCOEFF. LDUCO2 >= MAX(1,M) if ROWCOL = 'R'; LDUCO2 >= MAX(1,P) if ROWCOL = 'C'.

TOL1 DOUBLE PRECISION The tolerance to be used in determining the i-th row of T(s), where i = 1,2,...,porm. If the user sets TOL1 > 0, then the given value of TOL1 is used as an absolute tolerance; elements with absolute value less than TOL1 are considered neglijible. If the user sets TOL1 <= 0, then an implicitly computed, default tolerance, defined in the SLICOT Library routine TB01ZD, is used instead. TOL2 DOUBLE PRECISION The tolerance to be used to separate out a controllable subsystem of (A,B,C). If the user sets TOL2 > 0, then the given value of TOL2 is used as a lower bound for the reciprocal condition number (see the description of the argument RCOND in the SLICOT routine MB03OD); a (sub)matrix whose estimated condition number is less than 1/TOL2 is considered to be of full rank. If the user sets TOL2 <= 0, then an implicitly computed, default tolerance, defined in the SLICOT Library routine TB01UD, is used instead.

IWORK INTEGER array, dimension (N+MAX(M,P)) On exit, if INFO = 0, the first nonzero elements of IWORK(1:N) return the orders of the diagonal blocks of A. 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*(N + 1) + MAX(N*MP + 2*N + MAX(N,MP), 3*MP, PM)), where MP = M, PM = P, if ROWCOL = 'R'; MP = P, PM = M, if ROWCOL = 'C'. For optimum performance LDWORK should be larger.

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

The method for transfer matrices factorized by rows will be described here: T(s) factorized by columns is dealt with by operating on the dual of the original system. Each row of T(s) is simply a single-output relatively left prime polynomial matrix representation, so can be calculated by applying a simplified version of the Orthogonal Structure Theorem to a minimal state-space representation for the corresponding row of the given system. A minimal state-space representation is obtained using the Orthogonal Canonical Form to first separate out a completely controllable one for the overall system and then, for each row in turn, applying it again to the resulting dual SIMO (single-input multi-output) system. Note that the elements of the transformed matrix A so calculated are individually scaled in a way which guarantees a monic denominator polynomial.

[1] Williams, T.W.C. An Orthogonal Structure Theorem for Linear Systems. Control Systems Research Group, Kingston Polytechnic, Internal Report 82/2, 1982.

3 The algorithm requires 0(N ) operations.

None

**Program Text**

* TB04AD 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 MAXMP PARAMETER ( MAXMP = MAX( MMAX, PMAX ) ) INTEGER LDA, LDB, LDC, LDD, LDDCOE, LDUCO1, LDUCO2, $ NMAXP1 PARAMETER ( LDA = NMAX, LDB = NMAX, LDC = MAXMP, $ LDD = MAXMP, LDDCOE = MAXMP, LDUCO1 = MAXMP, $ LDUCO2 = MAXMP, NMAXP1 = NMAX+1 ) INTEGER LIWORK PARAMETER ( LIWORK = NMAX + MAXMP ) INTEGER LDWORK PARAMETER ( LDWORK = NMAX*( NMAX + 1 ) + $ MAX( NMAX*MAXMP + 2*NMAX + $ MAX( NMAX, MAXMP ), 3*MAXMP ) ) * .. Local Scalars .. DOUBLE PRECISION TOL1, TOL2 INTEGER I, II, IJ, INDBLK, INFO, J, JJ, KDCOEF, M, N, $ NR, P, PORM, PORP CHARACTER*1 ROWCOL CHARACTER*132 ULINE LOGICAL LROWCO * .. Local Arrays .. DOUBLE PRECISION A(LDA,NMAX), B(LDB,MAXMP), C(LDC,NMAX), $ D(LDD,MAXMP), DCOEFF(LDDCOE,NMAXP1), $ DWORK(LDWORK), UCOEFF(LDUCO1,LDUCO2,NMAXP1) INTEGER INDEX(MAXMP), IWORK(LIWORK) * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL TB04AD * .. 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, TOL1, TOL2, ROWCOL LROWCO = LSAME( ROWCOL, 'R' ) ULINE(1:20) = ' ' DO 20 I = 21, 132 ULINE(I:I) = '-' 20 CONTINUE IF ( N.LT.0 .OR. N.GT.NMAX ) THEN WRITE ( NOUT, FMT = 99986 ) N ELSE READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N ) IF ( M.LT.0 .OR. M.GT.MMAX ) THEN WRITE ( NOUT, FMT = 99985 ) M ELSE READ ( NIN, FMT = * ) ( ( B(I,J), I = 1,N ), J = 1,M ) IF ( P.LT.0 .OR. P.GT.PMAX ) THEN WRITE ( NOUT, FMT = 99984 ) 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 transfer matrix T(s) of (A,B,C,D). CALL TB04AD( ROWCOL, N, M, P, A, LDA, B, LDB, C, LDC, D, $ LDD, NR, INDEX, DCOEFF, LDDCOE, UCOEFF, $ LDUCO1, LDUCO2, TOL1, TOL2, IWORK, DWORK, $ LDWORK, INFO ) * IF ( INFO.NE.0 ) THEN WRITE ( NOUT, FMT = 99998 ) INFO ELSE WRITE ( NOUT, FMT = 99997 ) NR DO 40 I = 1, NR WRITE ( NOUT, FMT = 99996 ) ( A(I,J), J = 1,NR ) 40 CONTINUE WRITE ( NOUT, FMT = 99995 ) DO 60 I = 1, NR WRITE ( NOUT, FMT = 99996 ) ( B(I,J), J = 1,M ) 60 CONTINUE WRITE ( NOUT, FMT = 99994 ) DO 80 I = 1, P WRITE ( NOUT, FMT = 99996 ) ( C(I,J), J = 1,NR ) 80 CONTINUE INDBLK = 0 DO 100 I = 1, N IF ( IWORK(I).NE.0 ) INDBLK = INDBLK + 1 100 CONTINUE IF ( LROWCO ) THEN PORM = P PORP = M WRITE ( NOUT, FMT = 99993 ) INDBLK, $ ( IWORK(I), I = 1,INDBLK ) ELSE PORM = M PORP = P WRITE ( NOUT, FMT = 99992 ) INDBLK, $ ( IWORK(I), I = 1,INDBLK ) END IF WRITE ( NOUT, FMT = 99991 ) ( INDEX(I), I = 1,PORM ) WRITE ( NOUT, FMT = 99990 ) KDCOEF = 0 DO 120 I = 1, PORM KDCOEF = MAX( KDCOEF, INDEX(I) ) 120 CONTINUE KDCOEF = KDCOEF + 1 DO 160 II = 1, PORM DO 140 JJ = 1, PORP WRITE ( NOUT, FMT = 99989 ) II, JJ, $ ( UCOEFF(II,JJ,IJ), IJ = 1,KDCOEF ) WRITE ( NOUT, FMT = 99988 ) ULINE(1:7*KDCOEF+21) WRITE ( NOUT, FMT = 99987 ) $ ( DCOEFF(II,IJ), IJ = 1,KDCOEF ) 140 CONTINUE 160 CONTINUE END IF END IF END IF END IF STOP * 99999 FORMAT (' TB04AD EXAMPLE PROGRAM RESULTS',/1X) 99998 FORMAT (' INFO on exit from TB04AD = ',I2) 99997 FORMAT (' The order of the transformed state-space representatio', $ 'n = ',I2,//' The transformed state dynamics matrix A is ') 99996 FORMAT (20(1X,F8.4)) 99995 FORMAT (/' The transformed input/state matrix B is ') 99994 FORMAT (/' The transformed state/output matrix C is ') 99993 FORMAT (/' The controllability index of the transformed state-sp', $ 'ace representation = ',I2,//' The dimensions of the diag', $ 'onal blocks of the transformed A are ',/20(I5)) 99992 FORMAT (/' The observability index of the transformed state-spac', $ 'e representation = ',I2,//' The dimensions of the diagon', $ 'al blocks of the transformed A are ',/20(I5)) 99991 FORMAT (/' The degrees of the denominator polynomials are',/20(I5) $ ) 99990 FORMAT (/' The coefficients of polynomials in the transfer matri', $ 'x T(s) are ') 99989 FORMAT (/' element (',I2,',',I2,') is ',20(1X,F6.2)) 99988 FORMAT (1X,A) 99987 FORMAT (20X,20(1X,F6.2)) 99986 FORMAT (/' N is out of range.',/' N = ',I5) 99985 FORMAT (/' M is out of range.',/' M = ',I5) 99984 FORMAT (/' P is out of range.',/' P = ',I5) END

TB04AD EXAMPLE PROGRAM DATA 3 2 2 0.0 0.0 R -1.0 0.0 0.0 0.0 -2.0 0.0 0.0 0.0 -3.0 0.0 1.0 -1.0 1.0 1.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0 1.0

TB04AD EXAMPLE PROGRAM RESULTS The order of the transformed state-space representation = 3 The transformed state dynamics matrix A is -2.5000 -0.2887 -0.4082 -0.2887 -1.5000 -0.7071 -0.4082 -0.7071 -2.0000 The transformed input/state matrix B is -1.4142 -0.7071 0.0000 1.2247 0.0000 0.0000 The transformed state/output matrix C is 0.0000 0.8165 1.1547 0.0000 1.6330 0.5774 The controllability index of the transformed state-space representation = 2 The dimensions of the diagonal blocks of the transformed A are 2 1 The degrees of the denominator polynomials are 2 3 The coefficients of polynomials in the transfer matrix T(s) are element ( 1, 1) is 1.00 5.00 7.00 0.00 ----------------------------- 1.00 5.00 6.00 0.00 element ( 1, 2) is 0.00 1.00 3.00 0.00 ----------------------------- 1.00 5.00 6.00 0.00 element ( 2, 1) is 0.00 0.00 1.00 1.00 ----------------------------- 1.00 6.00 11.00 6.00 element ( 2, 2) is 1.00 8.00 20.00 15.00 ----------------------------- 1.00 6.00 11.00 6.00