TB01UX

Observable-unobservable decomposition of a standard system

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

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

Specification
      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, * )

Arguments

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.

Input/Output Parameters
  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).

Tolerances
  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.

Workspace
  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.

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 dual of the reduction
  algorithms 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**3 )  floating point operations.

Further Comments
  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.

Example

Program Text

  None
Program Data
  None
Program Results
  None

Return to Supporting Routines index