## TG01HX

### Orthogonal reduction of a descriptor system to a system with the same transfer-function matrix and with no uncontrollable finite eigenvalues

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

Purpose

```  Given the descriptor system (A-lambda*E,B,C) with the system
matrices A, E and B of the form

( A1 X1 )        ( E1 Y1 )        ( B1 )
A = (       ) ,  E = (       ) ,  B = (    ) ,
( 0  X2 )        ( 0  Y2 )        ( 0  )

where
- B is an L-by-M matrix, with B1 an N1-by-M  submatrix
- A is an L-by-N matrix, with A1 an N1-by-N1 submatrix
- E is an L-by-N matrix, with E1 an N1-by-N1 submatrix
with LBE nonzero sub-diagonals,
this routine reduces the pair (A1-lambda*E1,B1) to the form

Qc'*[B1 A1-lambda*E1]*diag(I,Zc) =

( Bc Ac-lambda*Ec      *         )
(                                ) ,
( 0     0         Anc-lambda*Enc )

where:
1) the pencil ( Bc Ac-lambda*Ec ) has full row rank NR for
all finite lambda and is in a staircase form with
_      _          _        _
( A1,0   A1,1  ...  A1,k-1   A1,k  )
(        _          _        _     )
( Bc Ac ) = (  0     A2,1  ...  A2,k-1   A2,k  ) ,  (1)
(              ...  _        _     )
(  0       0   ...  Ak,k-1   Ak,k  )

_          _        _
( E1,1  ...  E1,k-1   E1,k  )
(            _        _     )
Ec      = (   0   ...  E2,k-1   E2,k  ) ,         (2)
(       ...           _     )
(   0   ...    0      Ek,k  )
_
where Ai,i-1 is an rtau(i)-by-rtau(i-1) full row rank
_
matrix (with rtau(0) = M) and Ei,i is an rtau(i)-by-rtau(i)
upper triangular matrix.

2) the pencil Anc-lambda*Enc is regular of order N1-NR with Enc
upper triangular; this pencil contains the uncontrollable
finite eigenvalues of the pencil (A1-lambda*E1).

The transformations are applied to the whole matrices A, E, B
and C. The left and/or right orthogonal transformations Qc and Zc
performed to reduce the pencil can be optionally accumulated
in the matrices Q and Z, respectively.

The reduced order descriptor system (Ac-lambda*Ec,Bc,Cc) has no
uncontrollable finite eigenvalues and has the same
transfer-function matrix as the original system (A-lambda*E,B,C).

```
Specification
```      SUBROUTINE TG01HX( COMPQ, COMPZ, L, N, M, P, N1, LBE, A, LDA,
\$                   E, LDE, B, LDB, C, LDC, Q, LDQ, Z, LDZ, NR,
\$                   NRBLCK, RTAU, TOL, IWORK, DWORK, INFO )
C     .. Scalar Arguments ..
CHARACTER          COMPQ, COMPZ
INTEGER            INFO, L, LBE, LDA, LDB, LDC, LDE, LDQ, LDZ, M,
\$                   N, N1, NR, NRBLCK, P
DOUBLE PRECISION   TOL
C     .. Array Arguments ..
INTEGER            IWORK( * ), RTAU( * )
DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( LDC, * ),
\$                   DWORK( * ), E( LDE, * ), Q( LDQ, * ),
\$                   Z( LDZ, * )

```
Arguments

Mode Parameters

```  COMPQ   CHARACTER*1
= 'N':  do not compute Q;
= 'I':  Q is initialized to the unit matrix, and the
orthogonal matrix Q is returned;
= 'U':  Q must contain an orthogonal matrix Q1 on entry,
and the product Q1*Q is returned.

COMPZ   CHARACTER*1
= 'N':  do not compute Z;
= 'I':  Z is initialized to the unit matrix, and the
orthogonal matrix Z is returned;
= 'U':  Z must contain an orthogonal matrix Z1 on entry,
and the product Z1*Z is returned.

```
Input/Output Parameters
```  L       (input) INTEGER
The number of descriptor state equations; also the number
of rows of matrices A, E and B.  L >= 0.

N       (input) INTEGER
The dimension of the descriptor state vector; also the
number of columns of matrices A, E and 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; also the
number of rows of matrix C.  P >= 0.

N1      (input) INTEGER
The order of subsystem (A1-lambda*E1,B1,C1) to be reduced.
MIN(L,N) >= N1 >= 0.

LBE     (input) INTEGER
The number of nonzero sub-diagonals of submatrix E1.
MAX(0,N1-1) >= LBE >= 0.

A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the leading L-by-N part of this array must
contain the L-by-N state matrix A in the partitioned
form
( A1 X1 )
A = (       ) ,
( 0  X2 )

where A1 is N1-by-N1.
On exit, the leading L-by-N part of this array contains
the transformed state matrix,

( Ac  *   * )
Qc'*A*Zc = ( 0  Anc  * ) ,
( 0   0   * )

where Ac is NR-by-NR and Anc is (N1-NR)-by-(N1-NR).
The matrix ( Bc Ac ) is in the controllability
staircase form (1).

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 L-by-N descriptor matrix E in the partitioned
form
( E1 Y1 )
E = (       ) ,
( 0  Y2 )

where E1 is N1-by-N1 matrix with LBE nonzero
sub-diagonals.
On exit, the leading L-by-N part of this array contains
the transformed descriptor matrix

( Ec  *   * )
Qc'*E*Zc = ( 0  Enc  * ) ,
( 0   0   * )

where Ec is NR-by-NR and Enc is (N1-NR)-by-(N1-NR).
Both Ec and Enc are upper triangular and Enc is
nonsingular.

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 L-by-M input matrix B in the partitioned
form
( B1 )
B = (    ) ,
( 0  )

where B1 is N1-by-M.
On exit, the leading L-by-M part of this array contains
the transformed input matrix

( Bc )
Qc'*B = (    ) ,
( 0  )

where Bc is NR-by-M.
The matrix ( Bc Ac ) is in the controllability
staircase form (1).

LDB     INTEGER
The leading dimension of array B.  LDB >= MAX(1,L).

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*Zc.

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

Q       (input/output) DOUBLE PRECISION array, dimension (LDQ,L)
If COMPQ = 'N': Q is not referenced.
If COMPQ = 'I': on entry, Q need not be set;
on exit, the leading L-by-L part of this
array contains the orthogonal matrix Qc,
where Qc' is the product of transformations
which are applied to A, E, and B on
the left.
If COMPQ = 'U': on entry, the leading L-by-L part of this
array must contain an orthogonal matrix Q;
on exit, the leading L-by-L part of this
array contains the orthogonal matrix
Q*Qc.

LDQ     INTEGER
The leading dimension of array Q.
LDQ >= 1,        if COMPQ = 'N';
LDQ >= MAX(1,L), if COMPQ = 'U' or 'I'.

Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
If COMPZ = 'N': Z is not referenced.
If COMPZ = 'I': on entry, Z need not be set;
on exit, the leading N-by-N part of this
array contains the orthogonal matrix Zc,
which is the product of transformations
applied to A, E, and C on the right.
If COMPZ = 'U': on entry, the leading N-by-N part of this
array must contain an orthogonal matrix Z;
on exit, the leading N-by-N part of this
array contains the orthogonal matrix
Z*Zc.

LDZ     INTEGER
The leading dimension of array Z.
LDZ >= 1,        if COMPZ = 'N';
LDZ >= MAX(1,N), if COMPZ = 'U' or 'I'.

NR      (output) INTEGER
The order of the reduced matrices Ac and Ec, and the
number of rows of the reduced matrix Bc; also the order of
the controllable part of the pair (B, A-lambda*E).

NRBLCK  (output) INTEGER                      _
The number k, of full row rank blocks Ai,i in the
staircase form of the pencil (Bc Ac-lambda*Ec) (see (1)
and (2)).

RTAU    (output) INTEGER array, dimension (N1)
RTAU(i), for i = 1, ..., NRBLCK, is the row dimension of
_
the full row rank block Ai,i-1 in the staircase form (1).

```
Tolerances
```  TOL     DOUBLE PRECISION
The tolerance to be used in rank determinations when
transforming (A-lambda*E, B). If the user sets TOL > 0,
then the given value of TOL 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 is considered to be of full rank.  If the user sets
TOL <= 0, then an implicitly computed, default tolerance,
defined by  TOLDEF = L*N*EPS,  is used instead, where
EPS is the machine precision (see LAPACK Library routine
DLAMCH).  TOL < 1.

```
Workspace
```  IWORK   INTEGER array, dimension (M)

DWORK   DOUBLE PRECISION array, dimension (MAX(N,L,2*M))

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

```
Method
```  The subroutine is based on the reduction algorithm of [1].

```
References
```  [1] Varga, A.
Computation of Irreducible Generalized State-Space
Realizations.
Kybernetika, vol. 26, pp. 89-106, 1990.

```
Numerical Aspects
```  The algorithm is numerically backward stable and requires
0( N*N1**2 )  floating point operations.

```
```  None
```
Example

Program Text

```  None
```
Program Data
```  None
```
Program Results
```  None
```