**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