TG01BD

Orthogonal reduction of a descriptor system to the generalized Hessenberg form

Purpose

  To reduce the matrices A and E of the system pencil

          S =  ( A  B ) - lambda ( E  0 ) ,
               ( C  0 )          ( 0  0 )

  corresponding to the descriptor triple (A-lambda E,B,C),
  to generalized upper Hessenberg form using orthogonal
  transformations,

       Q' * A * Z = H,   Q' * E * Z = T,

  where H is upper Hessenberg, T is upper triangular, Q and Z
  are orthogonal, and ' means transpose. The corresponding
  transformations, written compactly as diag(Q',I) * S * diag(Z,I),
  are also applied to B and C, getting Q' * B and C * Z.

  The orthogonal matrices Q and Z are determined as products of
  Givens rotations. They may either be formed explicitly, or they
  may be postmultiplied into input matrices Q1 and Z1, so that

       Q1 * A * Z1' = (Q1*Q) * H * (Z1*Z)'
       Q1 * E * Z1' = (Q1*Q) * T * (Z1*Z)'.

      SUBROUTINE TG01BD( JOBE, COMPQ, COMPZ, N, M, P, ILO, IHI, A, LDA,
     $                   E, LDE, B, LDB, C, LDC, Q, LDQ, Z, LDZ, DWORK,
     $                   LDWORK, INFO )
C     .. Scalar Arguments ..
      CHARACTER          COMPQ, COMPZ, JOBE
      INTEGER            IHI, ILO, INFO, LDA, LDB, LDC, LDE, LDQ,
     $                   LDWORK, LDZ, M, N, P
C     .. Array Arguments ..
      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( LDC, * ),
     $                   DWORK(  * ), E( LDE, * ), Q( LDQ, * ),
     $                   Z( LDZ, * )

Mode Parameters

  JOBE    CHARACTER*1
          Specifies whether E is a general square or an upper
          triangular matrix, as follows:
          = 'G':  E is a general square matrix;
          = 'U':  E is an upper triangular matrix.

  COMPQ   CHARACTER*1
          Indicates what should be done with matrix Q, as follows:
          = 'N':  do not compute Q;
          = 'I':  Q is initialized to the unit matrix, and the
                  orthogonal matrix Q is returned;
          = 'V':  Q must contain an orthogonal matrix Q1 on entry,
                  and the product Q1*Q is returned.

  COMPZ   CHARACTER*1
          Indicates what should be done with matrix Z, as follows:
          = 'N':  do not compute Z;
          = 'I':  Z is initialized to the unit matrix, and the
                  orthogonal matrix Z is returned;
          = 'V':  Z must contain an orthogonal matrix Z1 on entry,
                  and the product Z1*Z is returned.

  N       (input) INTEGER
          The order of the matrices A, E, and the number of rows of
          the matrix B.  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.

  ILO     (input) INTEGER
  IHI     (input) INTEGER
          It is assumed that A and E are already upper triangular in
          rows and columns 1:ILO-1 and IHI+1:N.  ILO and IHI could
          normally be set by a previous call to LAPACK Library
          routine DGGBAL; otherwise they should be set to 1 and N,
          respectively.
          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
          If JOBE = 'U', the matrix E is assumed upper triangular.

  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
          On entry, the leading N-by-N part of this array must
          contain the state dynamics matrix A.
          On exit, the leading N-by-N part of this array contains
          the upper Hessenberg matrix H = Q' * A * Z. The elements
          below the first subdiagonal are set to zero.

  LDA     INTEGER
          The leading dimension of 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 descriptor matrix E. If JOBE = 'U', this
          matrix is assumed upper triangular.
          On exit, the leading N-by-N part of this array contains
          the upper triangular matrix T = Q' * E * Z. The elements
          below the diagonal are set to zero.

  LDE     INTEGER
          The leading dimension of 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 input/state matrix B.
          On exit, if M > 0, the leading N-by-M part of this array
          contains the transformed matrix Q' * B.
          The array B is not referenced if M = 0.

  LDB     INTEGER
          The leading dimension of array B.
          LDB >= MAX(1,N) if M > 0;  LDB >= 1 if M = 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, if P > 0, the leading P-by-N part of this array
          contains the transformed matrix C * Z.
          The array C is not referenced if P = 0.

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

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

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

  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, and on exit
                           it contains the orthogonal matrix Z,
                           which is the product of the Givens
                           transformations applied to A, E, and C
                           on the right;
          If COMPZ = 'V':  on entry, Z must contain an orthogonal
                           matrix Z1, and on exit this is
                           overwritten by Z1*Z.

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

  DWORK   DOUBLE PRECISION array, dimension (LDWORK)
          On exit, if INFO = 0, DWORK(1) contains the optimal value
          of LDWORK.

  LDWORK  INTEGER
          The dimension of the array DWORK.
          LDWORK >= 1,                          if JOBE = 'U';
          LDWORK >= MAX(1,IHI+1-ILO+MAX(NI,M)), if JOBE = 'G', where
          NI = N+1-ILO, if COMPQ = 'N', and NI = N, otherwise.
          For good performance, if JOBE = 'G', LDWORK must generally
          be larger, LDWORK >= MAX(1,IHI+1-ILO+MAX(NI,M)*NB), where
          NB is the optimal block size.

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

  First, this routine computes the QR factorization of E and applies
  the transformations to A, B, and possibly Q. Then, the routine
  reduces A to upper Hessenberg form, preserving E triangular, by
  an unblocked reduction [1], using two sequences of plane rotations
  applied alternately from the left and from the right. The
  corresponding transformations may be accumulated and/or applied
  to the matrices B and C. If JOBE = 'U', the initial reduction of E
  to upper triangular form is skipped.

  This routine is a modification and extension of the LAPACK Library
  routine DGGHRD [2].

  [1] Golub, G.H. and van Loan, C.F.
      Matrix Computations. Third Edition.
      M. D. Johns Hopkins University Press, Baltimore, 1996.

  [2] Anderson, E., Bai, Z., Bischof, C., Demmel, J., Dongarra, J.,
      Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A.,
      Ostrouchov, S., and Sorensen, D.
      LAPACK Users' Guide: Second Edition.
      SIAM, Philadelphia, 1995.

  None

Program Text

  None

  None

  None

Return to index