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.

Further Comments
  None
Example

Program Text

  None
Program Data
  None
Program Results
  None

Return to Supporting Routines index