TG01NX

Block-diagonal decomposition of a descriptor system in generalized real Schur form

Purpose

  To compute equivalence transformation matrices Q and Z which
  reduce the regular pole pencil A-lambda*E of the descriptor system
  (A-lambda*E,B,C), with (A,E) in a generalized real Schur form, to
  the block-diagonal form

             ( A1  0  )             ( E1  0  )
     Q*A*Z = (        ) ,   Q*E*Z = (        ) ,                 (1)
             ( 0   A2 )             ( 0   E2 )

  where the pair (Q*A*Z,Q*E*Z) is in a generalized real Schur form,
  with (A1,E1) and (A2,E2) having no common generalized eigenvalues.
  This decomposition corresponds to an additive spectral
  decomposition of the transfer-function matrix of the descriptor
  system as the sum of two terms containing the generalized
  eigenvalues of (A1,E1) and (A2,E2), respectively.

      SUBROUTINE TG01NX( JOBT, N, M, P, NDIM, A, LDA, E, LDE, B, LDB,
     $                   C, LDC, Q, LDQ, Z, LDZ, IWORK, INFO )
C     .. Scalar Arguments ..
      CHARACTER          JOBT
      INTEGER            INFO, LDA, LDB, LDC, LDE, LDQ, LDZ, M, N, NDIM,
     $                   P
C     .. Array Arguments ..
      INTEGER            IWORK(*)
      DOUBLE PRECISION   A(LDA,*), B(LDB,*), C(LDC,*), E(LDE,*),
     $                   Q(LDQ,*), Z(LDZ,*)

Mode Parameters

  JOBT    CHARACTER*1
          = 'D':  compute the direct transformation matrices;
          = 'I':  compute the inverse transformation matrices
                  inv(Q) and inv(Z).

  N       (input) INTEGER
          The number of rows of the matrix B, the number of columns
          of the matrix C and the order of the square matrices A
          and E.  N >= 0.

  M       (input) INTEGER
          The number of columns of the matrix B.  M >= 0.

  P       (input) INTEGER
          The number of rows of the matrix C.  P >= 0.

  NDIM    (input) INTEGER
          The dimension of the leading diagonal blocks of (A,E)
          having generalized eigenvalues distinct from those of the
          trailing diagonal block.  0 <= NDIM <= N.

  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 in a real Schur form.
          On exit, the leading N-by-N part of this array contains
          the transformed state matrix Q*A*Z (if JOBT = 'D') or
          inv(Q)*A*inv(Z) (if JOBT = 'I'), in the form (1), where
          A1 is a NDIM-by-NDIM matrix.

  LDA     INTEGER
          The leading dimension of the array A.  LDA >= MAX(1,N).

  E       (input/output) DOUBLE PRECISION array, dimension (LDE,N)
          On entry, the leading N-by-N part of this array must
          contain the N-by-N descriptor matrix E in upper triangular
          form.
          On exit, the leading N-by-N part of this array contains
          the transformed descriptor matrix Q*E*Z (if JOBT = 'D') or
          inv(Q)*E*inv(Z) (if JOBT = 'I'), in the form (1), where
          E1 is an NDIM-by-NDIM matrix.

  LDE     INTEGER
          The leading dimension of the array E.  LDE >= MAX(1,N).

  B       (input/output) DOUBLE PRECISION array, dimension (LDB,M)
          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 Q*B (if JOBT = 'D') or
          inv(Q)*B (if JOBT = 'I').

  LDB     INTEGER
          The leading dimension of the 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 state/output matrix C.
          On exit, the leading P-by-N part of this array contains
          the transformed matrix C*Z (if JOBT = 'D') or C*inv(Z)
          (if JOBT = 'I').

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

  Q       (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
          On entry, the leading N-by-N part of this array contains
          Q1, the orthogonal left transformation matrix Q used to
          reduce the pair (A,E) to the generalized real Schur form.
          On exit, the leading N-by-N part of this array contains
          the left transformation matrix Q = Q2*Q1, if JOBT = 'D',
          or its inverse inv(Q), if JOBT = 'I'.

  LDQ     INTEGER
          The leading dimension of the array Q.  LDQ >= MAX(1,N).

  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
          On entry, the leading N-by-N part of this array contains
          the orthogonal right transformation matrix Z1 used to
          reduce the pair (A,E) to the generalized real Schur form.
          On exit, the leading N-by-N part of this array contains
          the right transformation matrix Z = Z1*Z2, if JOBT = 'D',
          or its inverse inv(Z), if JOBT = 'I'.

  LDZ     INTEGER
          The leading dimension of the array Z.  LDZ >= MAX(1,N).

  IWORK   INTEGER array, dimension (N+6)

  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value;
          = 1:  the separation of the two diagonal blocks failed
                because of very close eigenvalues.

  For the separation, transformation matrices Q2 and Z2 of the form

          ( I -X )          ( I  Y )
     Q2 = (      ) ,   Z2 = (      )
          ( 0  I )          ( 0  I )

  are determined, such that Q2*A*Z2 and Q2*E*Z2 are block diagonal
  as in (1). X and Y are computed by solving generalized Sylvester
  equations.

  If we partition Q2*B and C*Z2 according to (1) in the form ( B1 )
                                                             ( B2 )
  and ( C1 C2 ), then (A1-lambda*E1,B1,C1) and (A2-lambda*E2,B2,C2)
  represent an additive spectral decomposition of the system
  transfer-function matrix.

  [1] Kagstrom, B. and Van Dooren, P.
      Additive decomposition of a transfer function with respect
      to a specified region.
      Proc. MTNS Symp., Brussels, 1989.

  The algorithm is numerically backward stable and requires
  0( N**3 )  floating point operations.

  None

Program Text

  None

  None

  None