**Purpose**

To reduce the pair (A,E) to a real generalized Schur form by using an orthogonal equivalence transformation (A,E) <-- (Q'*A*Z,Q'*E*Z) and to apply the transformation to the matrices B and C: B <-- Q'*B and C <-- C*Z.

SUBROUTINE TG01WD( N, M, P, A, LDA, E, LDE, B, LDB, C, LDC, $ Q, LDQ, Z, LDZ, ALPHAR, ALPHAI, BETA, DWORK, $ LDWORK, INFO ) C .. Scalar Arguments .. INTEGER INFO, LDA, LDB, LDC, LDE, LDQ, LDWORK, LDZ, $ M, N, P C .. Array Arguments .. DOUBLE PRECISION A(LDA,*), ALPHAI(*), ALPHAR(*), B(LDB,*), $ BETA(*), C(LDC,*), DWORK(*), E(LDE,*), $ Q(LDQ,*), Z(LDZ,*)

**Input/Output Parameters**

N (input) INTEGER The order of the original state-space representation, i.e., the order of the matrices A and E. N >= 0. M (input) INTEGER The number of system inputs, or of columns of B. M >= 0. P (input) INTEGER The number of system outputs, or of rows of 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 original state dynamics matrix A. On exit, the leading N-by-N part of this array contains the matrix Q' * A * Z in an upper quasi-triangular form. 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 original descriptor matrix E. On exit, the leading N-by-N part of this array contains the matrix Q' * E * Z in an upper triangular form. 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 matrix B. On exit, the leading N-by-M part of this array contains the transformed input matrix Q' * B. LDB INTEGER The leading dimension of 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 output matrix C. On exit, the leading P-by-N part of this array contains the transformed output matrix C * Z. LDC INTEGER The leading dimension of array C. LDC >= MAX(1,P). Q (output) DOUBLE PRECISION array, dimension (LDQ,N) The leading N-by-N part of this array contains the left orthogonal transformation matrix used to reduce (A,E) to the real generalized Schur form. The columns of Q are the left generalized Schur vectors of the pair (A,E). LDQ INTEGER The leading dimension of array Q. LDQ >= max(1,N). Z (output) DOUBLE PRECISION array, dimension (LDZ,N) The leading N-by-N part of this array contains the right orthogonal transformation matrix used to reduce (A,E) to the real generalized Schur form. The columns of Z are the right generalized Schur vectors of the pair (A,E). LDZ INTEGER The leading dimension of array Z. LDZ >= max(1,N). ALPHAR (output) DOUBLE PRECISION array, dimension (N) ALPHAI (output) DOUBLE PRECISION array, dimension (N) BETA (output) DOUBLE PRECISION array, dimension (N) On exit, if INFO = 0, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will be the generalized eigenvalues. ALPHAR(j) + ALPHAI(j)*i, and BETA(j), j=1,...,N, are the diagonals of the complex Schur form that would result if the 2-by-2 diagonal blocks of the real Schur form of (A,E) were further reduced to triangular form using 2-by-2 complex unitary transformations. If ALPHAI(j) is zero, then the j-th eigenvalue is real; if positive, then the j-th and (j+1)-st eigenvalues are a complex conjugate pair, with ALPHAI(j+1) negative.

DWORK DOUBLE PRECISION array, dimension (LDWORK) On exit, if INFO = 0, DWORK(1) returns the optimal value of LDWORK. LDWORK INTEGER The dimension of working array DWORK. LDWORK >= 8*N+16. For optimum performance LDWORK should be larger.

INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value; > 0: if INFO = i, the QZ algorithm failed to compute the generalized real Schur form; elements i+1:N of ALPHAR, ALPHAI, and BETA should be correct.

The pair (A,E) is reduced to a real generalized Schur form using an orthogonal equivalence transformation (A,E) <-- (Q'*A*Z,Q'*E*Z) and the transformation is applied to the matrices B and C: B <-- Q'*B and C <-- C*Z.

3 The algorithm requires about 25N floating point operations.

None

**Program Text**

None

None

None