## TB01PD

### Minimal, controllable or observable block Hessenberg realization for a given state-space representation

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

Purpose

```  To find a reduced (controllable, observable, or minimal) state-
space representation (Ar,Br,Cr) for any original state-space
representation (A,B,C). The matrix Ar is in upper block
Hessenberg form.

```
Specification
```      SUBROUTINE TB01PD( JOB, EQUIL, N, M, P, A, LDA, B, LDB, C, LDC,
\$                   NR, TOL, IWORK, DWORK, LDWORK, INFO )
C     .. Scalar Arguments ..
CHARACTER         EQUIL, JOB
INTEGER           INFO, LDA, LDB, LDC, LDWORK, M, N, NR, P
DOUBLE PRECISION  TOL
C     .. Array Arguments ..
INTEGER           IWORK(*)
DOUBLE PRECISION  A(LDA,*), B(LDB,*), C(LDC,*), DWORK(*)

```
Arguments

Mode Parameters

```  JOB     CHARACTER*1
Indicates whether the user wishes to remove the
uncontrollable and/or unobservable parts as follows:
= 'M':  Remove both the uncontrollable and unobservable
parts to get a minimal state-space representation;
= 'C':  Remove the uncontrollable part only to get a
controllable state-space representation;
= 'O':  Remove the unobservable part only to get an
observable state-space representation.

EQUIL   CHARACTER*1
Specifies whether the user wishes to preliminarily balance
the triplet (A,B,C) as follows:
= 'S':  Perform balancing (scaling);
= 'N':  Do not perform balancing.

```
Input/Output Parameters
```  N       (input) INTEGER
The order of the original state-space representation, i.e.
the order of the 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 Ar of a
minimal, controllable, or observable realization for the
original system, depending on the value of JOB, JOB = 'M',
JOB = 'C', or JOB = 'O', respectively.

LDA     INTEGER
The leading dimension of array A.  LDA >= 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/state matrix B; if JOB = 'M',
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 transformed input/state matrix Br of a minimal,
controllable, or observable realization for the original
system, depending on the value of JOB, JOB = 'M',
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 state/output matrix C; if JOB = 'M',
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 a minimal,
controllable, or observable realization for the original
system, depending on the value of JOB, JOB = 'M',
JOB = 'C', or JOB = 'O', respectively.
If JOB = 'M', 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 state-space representation
(Ar,Br,Cr) of a minimal, controllable, or observable
realization for the original system, depending on
JOB = 'M', JOB = 'C', or JOB = 'O'.

```
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)).
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.

```
Method
```  If JOB = 'M', the matrices A and B are operated on by orthogonal
similarity transformations (made up of products of Householder
transformations) so as to produce an upper block Hessenberg matrix
A1 and a matrix B1 with all but its first rank(B) rows zero; this
separates out the controllable part of the original system.
Applying the same algorithm to the dual of this subsystem,
therefore separates out the controllable and observable (i.e.
minimal) part of the original system representation, with the
final Ar upper block Hessenberg (after using pertransposition).
If JOB = 'C', or JOB = 'O', only the corresponding part of the
above procedure is applied.

```
References
```   Van Dooren, P.
The Generalized Eigenstructure Problem in Linear System
Theory. (Algorithm 1)
IEEE Trans. Auto. Contr., AC-26, pp. 111-129, 1981.

```
Numerical Aspects
```                            3
The algorithm requires 0(N ) operations and is backward stable.

```
```  None
```
Example

Program Text

```*     TB01PD 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
PARAMETER        ( LDA = NMAX, LDB = NMAX, LDC = MAXMP )
INTEGER          LIWORK
PARAMETER        ( LIWORK = NMAX+MAXMP )
INTEGER          LDWORK
PARAMETER        ( LDWORK = NMAX+MAX( NMAX, 3*MAXMP ) )
*     .. Local Scalars ..
DOUBLE PRECISION TOL
INTEGER          I, INFO, J, M, N, NR, P
CHARACTER        JOB, EQUIL
*     .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), B(LDB,MAXMP), C(LDC,NMAX),
\$                 DWORK(LDWORK)
INTEGER          IWORK(LIWORK)
*     .. External Subroutines ..
EXTERNAL         TB01PD
*     .. 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, JOB, EQUIL
IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99990 ) 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 = 99989 ) 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 = 99988 ) P
ELSE
READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,N ), I = 1,P )
*              Find a minimal ssr for (A,B,C).
CALL TB01PD( JOB, EQUIL, N, M, P, A, LDA, B, LDB, C, LDC,
\$                      NR, TOL, IWORK, DWORK, LDWORK, INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
WRITE ( NOUT, FMT = 99997 ) NR
WRITE ( NOUT, FMT = 99996 )
DO 20 I = 1, NR
WRITE ( NOUT, FMT = 99995 ) ( A(I,J), J = 1,NR )
20             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 60 I = 1, P
WRITE ( NOUT, FMT = 99995 ) ( C(I,J), J = 1,NR )
60             CONTINUE
END IF
END IF
END IF
END IF
STOP
*
99999 FORMAT (' TB01PD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from TB01PD = ',I2)
99997 FORMAT (' The order of the minimal realization = ',I2)
99996 FORMAT (/' The transformed state dynamics matrix of a minimal re',
\$       'alization is ')
99995 FORMAT (20(1X,F8.4))
99993 FORMAT (/' The transformed input/state matrix of a minimal reali',
\$       'zation is ')
99992 FORMAT (/' The transformed state/output matrix of a minimal real',
\$       'ization is ')
99990 FORMAT (/' N is out of range.',/' N = ',I5)
99989 FORMAT (/' M is out of range.',/' M = ',I5)
99988 FORMAT (/' P is out of range.',/' P = ',I5)
END
```
Program Data
``` TB01PD EXAMPLE PROGRAM DATA
3     1     2     0.0     M     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
```
Program Results
``` TB01PD EXAMPLE PROGRAM RESULTS

The order of the minimal realization =  3

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

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

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