**Purpose**

To compute an orthogonal transformation matrix Z which reduces the N-th order system (A,B,C) to the form ( Ano * ) ( Bno ) Z'*A*Z = ( ) , Z'*B = ( ) , ( 0 Ao ) ( Bo ) C*Z = ( 0 Co ) , where the NOBSV-th order system (Ao,Bo,Co) is observable. The matrix Ano of order N-NOBSV contains the unobservable eigenvalues of A. The pencil ( Ao-lambda*I ) has full column rank NOBSV for all ( Co ) lambda, and is in a staircase form, with _ _ _ _ ( Ak,k Ak,k-1 ... Ak,2 Ak,1 ) ( _ _ _ _ ) ( Ao ) = ( Ak-1,k Ak-1,k-1 ... Ak-1,2 Ak-1,1 ) , (1) ( Co ) ( : : ... _ : _ : ) ( 0 0 ... A1,2 A1,1 ) ( _ ) ( 0 0 ... 0 A0,1 ) _ where Ai-1,i is a CTAU(i-1)-by-CTAU(i) full column rank matrix (with CTAU(0) = P). The orthogonal transformation Z, performed to reduce the system matrices, can be optionally accumulated. The reduced order system (Ao,Bo,Co) has the same transfer-function matrix as the original system (A,B,C).

SUBROUTINE TB01UX( COMPZ, N, M, P, A, LDA, B, LDB, C, LDC, Z, LDZ, $ NOBSV, NLBLCK, CTAU, TOL, IWORK, DWORK, INFO ) C .. Scalar Arguments .. CHARACTER COMPZ INTEGER INFO, LDA, LDB, LDC, LDZ, M, N, NLBLCK, NOBSV, $ P DOUBLE PRECISION TOL C .. Array Arguments .. INTEGER CTAU( * ), IWORK( * ) DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ), $ DWORK( * ), Z( LDZ, * )

**Mode Parameters**

COMPZ CHARACTER*1 = 'N': do not compute Z; = 'I': Z is initialized to the unit matrix, and the orthogonal matrix Z is returned.

N (input) INTEGER The dimension of the system state vector; also the order of the square matrix A, the number of rows of the matrix B and the number of columns of the matrix C. N >= 0. M (input) INTEGER The dimension of system input vector; also the number of columns of the matrix B. M >= 0. P (input) INTEGER The dimension of system output vector; also the number of rows of the 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 N-by-N state matrix A. On exit, the leading N-by-N part of this array contains the transformed state matrix Z'*A*Z, ( Ano * ) Z'*A*Z = ( ) , ( 0 Ao ) where Ao is NOBSV-by-NOBSV and Ano is (N-NOBSV)-by-(N-NOBSV). The matrix ( Ao ) is in the observability staircase ( Co ) form (1). LDA INTEGER The leading dimension of the 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 N-by-M input matrix B. On exit, the leading N-by-M part of this array contains the transformed input matrix Z'*B. LDB INTEGER The leading dimension of the array B. LDB >= MAX(1,N) if M > 0 or P > 0; LDB >= 1 if M = 0 and P = 0. 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*Z = ( 0 Co ) , where Co is P-by-NOBSV. The matrix ( Ao ) is in the observability staircase ( Co ) form (1). LDC INTEGER The leading dimension of the array C. LDC >= MAX(1,M,P) if N > 0; LDC >= 1 if N = 0. Z (input/output) DOUBLE PRECISION array, dimension (LDZ,*) 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 Z, i.e., the product of the transformations applied to A and C on the right. LDZ INTEGER The leading dimension of the array Z. LDZ >= 1, if COMPZ = 'N'; LDZ >= MAX(1,N), if COMPZ = 'I'. NOBSV (output) INTEGER The order of the reduced matrix Ao, and the number of columns of the reduced matrix Co; also, the order of the observable part of the pair (C, A-lambda*I). NLBLCK (output) INTEGER _ The number k, of full column rank blocks Ai-1,i in the staircase form of the pencil (Ao-lambda*I) (see (1)). ( Co ) CTAU (output) INTEGER array, dimension (N) CTAU(i), for i = 1, ..., NLBLCK, is the column dimension _ of the full column rank block Ai-1,i in the staircase form (1).

TOL DOUBLE PRECISION The tolerance to be used in rank determinations when transforming the pair (A,C). 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 = N*N*EPS, is used instead, where EPS is the machine precision (see LAPACK Library routine DLAMCH). TOL < 1.

IWORK INTEGER array, dimension (P) DWORK DOUBLE PRECISION array, dimension (N+MAX(1, N, 3*P, M)) On exit, if INFO = 0, DWORK(1) returns the optimal value of LDWORK.

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

The subroutine is based on the dual of the reduction algorithms of [1].

[1] Varga, A. 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 system matrices A and C are badly scaled, it is generally recommendable to scale them with the SLICOT routine TB01ID, before calling TG01UX.

**Program Text**

None

None

None