### Left/right polynomial matrix representation of a given state-space representation

[Specification] [Arguments] [Method] [References] [Comments] [Example]

Purpose

```  To find a relatively prime left polynomial matrix representation
inv(P(s))*Q(s) or right polynomial matrix representation
Q(s)*inv(P(s)) with the same transfer matrix T(s) as that of a
given state-space representation, i.e.

inv(P(s))*Q(s) = Q(s)*inv(P(s)) = T(s) = C*inv(s*I-A)*B + D.

```
Specification
```      SUBROUTINE TB03AD( LERI, EQUIL, N, M, P, A, LDA, B, LDB, C, LDC,
\$                   D, LDD, NR, INDEX, PCOEFF, LDPCO1, LDPCO2,
\$                   QCOEFF, LDQCO1, LDQCO2, VCOEFF, LDVCO1, LDVCO2,
\$                   TOL, IWORK, DWORK, LDWORK, INFO )
C     .. Scalar Arguments ..
CHARACTER         EQUIL, LERI
INTEGER           INFO, LDA, LDB, LDC, LDD, LDPCO1, LDPCO2,
\$                  LDQCO1, LDQCO2, LDVCO1, LDVCO2, LDWORK, M, N,
\$                  NR, P
DOUBLE PRECISION  TOL
C     .. Array Arguments ..
INTEGER           INDEX(*), IWORK(*)
DOUBLE PRECISION  A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*),
\$                  DWORK(*), PCOEFF(LDPCO1,LDPCO2,*),
\$                  QCOEFF(LDQCO1,LDQCO2,*), VCOEFF(LDVCO1,LDVCO2,*)

```
Arguments

Mode Parameters

```  LERI    CHARACTER*1
Indicates whether the left polynomial matrix
representation or the right polynomial matrix
representation is required as follows:
= 'L':  A left matrix fraction is required;
= 'R':  A right matrix fraction is required.

EQUIL   CHARACTER*1
Specifies whether the user wishes to balance the triplet
(A,B,C), before computing a minimal state-space
representation, as follows:
= 'S':  Perform balancing (scaling);
= 'N':  Do not perform balancing.

```
Input/Output Parameters
```  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 Amin of a
minimal realization for the original system.

LDA     INTEGER
The leading dimension of array A.  LDA >= MAX(1,N).

B       (input/output) DOUBLE PRECISION array, dimension
(LDB,MAX(M,P))
On entry, the leading N-by-M part of this array must
contain the original input/state matrix B; 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 Bmin.

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; 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 Cmin.

LDC     INTEGER
The leading dimension of array C.  LDC >= MAX(1,M,P).

D       (input) DOUBLE PRECISION array, dimension (LDD,MAX(M,P))
The leading P-by-M part of this array must contain the
original direct transmission matrix D; 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,M,P).

NR      (output) INTEGER
The order of the minimal state-space representation
(Amin,Bmin,Cmin).

INDEX   (output) INTEGER array, dimension (P), if LERI = 'L', or
dimension (M), if LERI = 'R'.
If LERI = 'L', INDEX(I), I = 1,2,...,P, contains the
maximum degree of the polynomials in the I-th row of the
denominator matrix P(s) of the left polynomial matrix
representation.
These elements are ordered so that
INDEX(1) >= INDEX(2) >= ... >= INDEX(P).
If LERI = 'R', INDEX(I), I = 1,2,...,M, contains the
maximum degree of the polynomials in the I-th column of
the denominator matrix P(s) of the right polynomial
matrix representation.
These elements are ordered so that
INDEX(1) >= INDEX(2) >= ... >= INDEX(M).

PCOEFF  (output) DOUBLE PRECISION array, dimension
(LDPCO1,LDPCO2,N+1)
If LERI = 'L' then porm = P, otherwise porm = M.
The leading porm-by-porm-by-kpcoef part of this array
contains the coefficients of the denominator matrix P(s),
where kpcoef = MAX(INDEX(I)) + 1.
PCOEFF(I,J,K) is the coefficient in s**(INDEX(iorj)-K+1)
of polynomial (I,J) of P(s), where K = 1,2,...,kpcoef; if
LERI = 'L' then iorj = I, otherwise iorj = J.
Thus for LERI = 'L', P(s) =
diag(s**INDEX(I))*(PCOEFF(.,.,1)+PCOEFF(.,.,2)/s+...).

LDPCO1  INTEGER
The leading dimension of array PCOEFF.
LDPCO1 >= MAX(1,P), if LERI = 'L';
LDPCO1 >= MAX(1,M), if LERI = 'R'.

LDPCO2  INTEGER
The second dimension of array PCOEFF.
LDPCO2 >= MAX(1,P), if LERI = 'L';
LDPCO2 >= MAX(1,M), if LERI = 'R'.

QCOEFF  (output) DOUBLE PRECISION array, dimension
(LDQCO1,LDQCO2,N+1)
If LERI = 'L' then porp = M, otherwise porp = P.
If LERI = 'L', the leading porm-by-porp-by-kpcoef part
of this array contains the coefficients of the numerator
matrix Q(s).
If LERI = 'R', the leading porp-by-porm-by-kpcoef part
of this array contains the coefficients of the numerator
matrix Q(s).
QCOEFF(I,J,K) is defined as for PCOEFF(I,J,K).

LDQCO1  INTEGER
The leading dimension of array QCOEFF.
LDQCO1 >= MAX(1,P),   if LERI = 'L';
LDQCO1 >= MAX(1,M,P), if LERI = 'R'.

LDQCO2  INTEGER
The second dimension of array QCOEFF.
LDQCO2 >= MAX(1,M),   if LERI = 'L';
LDQCO2 >= MAX(1,M,P), if LERI = 'R'.

VCOEFF  (output) DOUBLE PRECISION array, dimension
(LDVCO1,LDVCO2,N+1)
The leading porm-by-NR-by-kpcoef part of this array
contains the coefficients of the intermediate matrix V(s).
VCOEFF(I,J,K) is defined as for PCOEFF(I,J,K).

LDVCO1  INTEGER
The leading dimension of array VCOEFF.
LDVCO1 >= MAX(1,P), if LERI = 'L';
LDVCO1 >= MAX(1,M), if LERI = 'R'.

LDVCO2  INTEGER
The second dimension of array VCOEFF.  LDVCO2 >= MAX(1,N).

```
Tolerances
```  TOL     DOUBLE PRECISION
The tolerance to be used in rank determination when
transforming (A, B, C). If the user sets TOL > 0, then
the given value of TOL 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/TOL is considered to be of full rank.  If the user sets
TOL <= 0, then an implicitly computed, default tolerance
(determined by the SLICOT routine TB01UD) is used instead.

```
Workspace
```  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 + MAX(N, 3*M, 3*P), PM*(PM + 2))
where  PM = P, if LERI = 'L';
PM = M, if LERI = 'R'.
For optimum performance LDWORK should be larger.

```
Error Indicator
```  INFO    INTEGER
= 0:  successful exit;
< 0:  if INFO = -i, the i-th argument had an illegal
value;
= 1:  if a singular matrix was encountered during the
computation of V(s);
= 2:  if a singular matrix was encountered during the
computation of P(s).

```
Method
```  The method for a left matrix fraction will be described here:
right matrix fractions are dealt with by constructing a left
fraction for the dual of the original system. The first step is to
obtain, by means of orthogonal similarity transformations, a
minimal state-space representation (Amin,Bmin,Cmin,D) for the
original system (A,B,C,D), where Amin is lower block Hessenberg
with all its superdiagonal blocks upper triangular and Cmin has
all but its first rank(C) columns zero.  The number and dimensions
of the blocks of Amin now immediately yield the row degrees of
P(s) with P(s) row proper: furthermore, the P-by-NR polynomial
matrix V(s) (playing a similar role to S(s) in Wolovich's
Structure Theorem) can be calculated a column block at a time, in
reverse order, from Amin. P(s) is then found as if it were the
O-th column block of V(s) (using Cmin as well as Amin), while
Q(s) = (V(s) * Bmin) + (P(s) * D). Finally, a special similarity
transformation is used to put Amin in an upper block Hessenberg
form.

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

 Patel, R.V.
On Computing Matrix Fraction Descriptions and Canonical
Forms of Linear Time-Invariant Systems.
UMIST Control Systems Centre Report 489, 1980.
(Algorithms 1 and 2, extensively modified).

```
Numerical Aspects
```                            3
The algorithm requires 0(N ) operations.

```
```  None
```
Example

Program Text

```*     TB03AD 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, LDPCO1, LDPCO2, LDQCO1,
\$                 LDQCO2, LDVCO1, LDVCO2, NMAXP1
PARAMETER        ( LDA = NMAX, LDB = NMAX, LDC = MAXMP,
\$                   LDD = MAXMP, LDPCO1 = MAXMP, LDPCO2 = MAXMP,
\$                   LDQCO1 = MAXMP, LDQCO2 = MAXMP, LDVCO1 = MAXMP,
\$                   LDVCO2 = NMAX, NMAXP1 = NMAX+1 )
INTEGER          LIWORK
PARAMETER        ( LIWORK = NMAX + MAXMP )
INTEGER          LDWORK
PARAMETER        ( LDWORK = MAX( NMAX + MAX( NMAX, 3*MAXMP ),
\$                                 MAXMP*( MAXMP + 2 ) ) )
*     .. Local Scalars ..
DOUBLE PRECISION TOL
INTEGER          I, INDBLK, INFO, J, K, KPCOEF, M, N, NR, P, PORM,
\$                 PORP
CHARACTER*1      EQUIL, LERI
LOGICAL          LLERI
*     .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), B(LDB,MAXMP), C(LDC,NMAX),
\$                 D(LDD,MAXMP), DWORK(LDWORK),
\$                 PCOEFF(LDPCO1,LDPCO2,NMAXP1),
\$                 QCOEFF(LDQCO1,LDQCO2,NMAXP1),
\$                 VCOEFF(LDVCO1,LDVCO2,NMAXP1)
INTEGER          INDEX(MAXMP), IWORK(LIWORK)
*     .. External Functions ..
LOGICAL          LSAME
EXTERNAL         LSAME
*     .. External Subroutines ..
*     .. 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, LERI, EQUIL
LLERI = LSAME( LERI, 'L' )
IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99987 ) 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 = 99986 ) 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 = 99985 ) 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 right pmr which is equivalent to the ssr
*              C*inv(sI-A)*B+D.
CALL TB03AD( LERI, EQUIL, N, M, P, A, LDA, B, LDB, C,
\$                      LDC, D, LDD, NR, INDEX, PCOEFF, LDPCO1,
\$                      LDPCO2, QCOEFF, LDQCO1, LDQCO2, VCOEFF,
\$                      LDVCO1, LDVCO2, TOL, IWORK, DWORK, LDWORK,
\$                      INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
WRITE ( NOUT, FMT = 99997 ) NR
DO 20 I = 1, NR
WRITE ( NOUT, FMT = 99996 ) ( A(I,J), J = 1,NR )
20             CONTINUE
INDBLK = 0
DO 40 I = 1, N
IF ( IWORK(I).NE.0 ) INDBLK = INDBLK + 1
40             CONTINUE
WRITE ( NOUT, FMT = 99995 ) ( IWORK(I), I = 1,INDBLK )
WRITE ( NOUT, FMT = 99994 )
DO 60 I = 1, NR
WRITE ( NOUT, FMT = 99996 ) ( B(I,J), J = 1,M )
60             CONTINUE
WRITE ( NOUT, FMT = 99993 )
DO 80 I = 1, P
WRITE ( NOUT, FMT = 99996 ) ( C(I,J), J = 1,NR )
80             CONTINUE
IF ( LLERI ) THEN
PORM = P
PORP = M
WRITE ( NOUT, FMT = 99992 ) INDBLK
ELSE
PORM = M
PORP = P
WRITE ( NOUT, FMT = 99991 ) INDBLK
END IF
WRITE ( NOUT, FMT = 99990 ) ( INDEX(I), I = 1,PORM )
KPCOEF = 0
DO 100 I = 1, PORM
KPCOEF = MAX( KPCOEF, INDEX(I) )
100             CONTINUE
KPCOEF = KPCOEF + 1
WRITE ( NOUT, FMT = 99989 )
DO 140 I = 1, PORM
DO 120 J = 1, PORM
WRITE ( NOUT, FMT = 99996 )
\$                        ( PCOEFF(I,J,K), K = 1,KPCOEF )
120                CONTINUE
140             CONTINUE
WRITE ( NOUT, FMT = 99988 )
IF ( LLERI ) THEN
DO 180 I = 1, PORM
DO 160 J = 1, PORP
WRITE ( NOUT, FMT = 99996 )
\$                           ( QCOEFF(I,J,K), K = 1,KPCOEF )
160                   CONTINUE
180                CONTINUE
ELSE
DO 220 I = 1, PORP
DO 200 J = 1, PORM
WRITE ( NOUT, FMT = 99996 )
\$                           ( QCOEFF(I,J,K), K = 1,KPCOEF )
200                   CONTINUE
220                CONTINUE
END IF
END IF
END IF
END IF
END IF
STOP
*
99999 FORMAT (' TB03AD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from TB03AD = ',I2)
99997 FORMAT (' The order of the minimal state-space representation = ',
\$       I2,//' The transformed state dynamics matrix of a minimal',
\$       ' realization is ')
99996 FORMAT (20(1X,F8.4))
99995 FORMAT (/' and the dimensions of its diagonal blocks are ',/20(I5)
\$       )
99994 FORMAT (/' The transformed input/state matrix of a minimal reali',
\$       'zation is ')
99993 FORMAT (/' The transformed state/output matrix of a minimal real',
\$       'ization is ')
99992 FORMAT (/' The observability index of the transformed minimal sy',
\$       'stem representation = ',I2)
99991 FORMAT (/' The controllability index of the transformed minimal ',
\$       'system representation = ',I2)
99990 FORMAT (/' INDEX is ',/20(I5))
99989 FORMAT (/' The denominator matrix P(s) is ')
99988 FORMAT (/' The numerator matrix Q(s) is ')
99987 FORMAT (/' N is out of range.',/' N = ',I5)
99986 FORMAT (/' M is out of range.',/' M = ',I5)
99985 FORMAT (/' P is out of range.',/' P = ',I5)
END
```
Program Data
``` TB03AD EXAMPLE PROGRAM DATA
3     1     2     0.0     R     N
1.0   2.0   0.0
4.0  -1.0   0.0
0.0   0.0   1.0
1.0   0.0   1.0
0.0   1.0  -1.0
0.0   0.0   1.0
0.0   1.0
```
Program Results
``` TB03AD EXAMPLE PROGRAM RESULTS

The order of the minimal state-space representation =  3

The transformed state dynamics matrix of a minimal realization is
1.0000  -1.4142   0.0000
-2.8284  -1.0000   2.8284
0.0000   1.4142   1.0000

and the dimensions of its diagonal blocks are
1    1    1

The transformed input/state matrix of a minimal realization is
-1.4142
0.0000
0.0000

The transformed state/output matrix of a minimal realization is
0.7071   1.0000   0.7071
-0.7071   0.0000  -0.7071

The controllability index of the transformed minimal system representation =  3

INDEX is
3

The denominator matrix P(s) is
0.1768  -0.1768  -1.5910   1.5910

The numerator matrix Q(s) is
0.0000  -0.1768   0.7071   0.8839
0.1768   0.0000  -1.5910   0.0000
```