dhgeqz - implement a single-/double-shift version of the QZ method for finding the generalized eigenvalues w(j)=(ALPHAR(j) + i*ALPHAI(j))/BETAR(j) of the equation det( A-w(i) B ) = 0 In addition, the pair A,B may be reduced to generalized Schur form
SUBROUTINE DHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, * ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO) CHARACTER * 1 JOB, COMPQ, COMPZ INTEGER N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK, INFO DOUBLE PRECISION A(LDA,*), B(LDB,*), ALPHAR(*), ALPHAI(*), BETA(*), Q(LDQ,*), Z(LDZ,*), WORK(*)
SUBROUTINE DHGEQZ_64( JOB, COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, * LDB, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO) CHARACTER * 1 JOB, COMPQ, COMPZ INTEGER*8 N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK, INFO DOUBLE PRECISION A(LDA,*), B(LDB,*), ALPHAR(*), ALPHAI(*), BETA(*), Q(LDQ,*), Z(LDZ,*), WORK(*)
SUBROUTINE HGEQZ( JOB, COMPQ, COMPZ, [N], ILO, IHI, A, [LDA], B, * [LDB], ALPHAR, ALPHAI, BETA, Q, [LDQ], Z, [LDZ], [WORK], [LWORK], * [INFO]) CHARACTER(LEN=1) :: JOB, COMPQ, COMPZ INTEGER :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK, INFO REAL(8), DIMENSION(:) :: ALPHAR, ALPHAI, BETA, WORK REAL(8), DIMENSION(:,:) :: A, B, Q, Z
SUBROUTINE HGEQZ_64( JOB, COMPQ, COMPZ, [N], ILO, IHI, A, [LDA], B, * [LDB], ALPHAR, ALPHAI, BETA, Q, [LDQ], Z, [LDZ], [WORK], [LWORK], * [INFO]) CHARACTER(LEN=1) :: JOB, COMPQ, COMPZ INTEGER(8) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK, INFO REAL(8), DIMENSION(:) :: ALPHAR, ALPHAI, BETA, WORK REAL(8), DIMENSION(:,:) :: A, B, Q, Z
#include <sunperf.h>
void dhgeqz(char job, char compq, char compz, int n, int ilo, int ihi, double *a, int lda, double *b, int ldb, double *alphar, double *alphai, double *beta, double *q, int ldq, double *z, int ldz, int *info);
void dhgeqz_64(char job, char compq, char compz, long n, long ilo, long ihi, double *a, long lda, double *b, long ldb, double *alphar, double *alphai, double *beta, double *q, long ldq, double *z, long ldz, long *info);
dhgeqz implements a single-/double-shift version of the QZ method for finding the generalized eigenvalues B is upper triangular, and A is block upper triangular, where the diagonal blocks are either 1-by-1 or 2-by-2, the 2-by-2 blocks having complex generalized eigenvalues (see the description of the argument JOB.)
If JOB='S', then the pair (A,B) is simultaneously reduced to Schur
form by applying one orthogonal tranformation (usually called Q) on
the left and another (usually called Z) on the right. The 2-by-2
upper-triangular diagonal blocks of B corresponding to 2-by-2 blocks
of A will be reduced to positive diagonal matrices. (I.e.,
if A(j+1,j) is non-zero, then B(j+1,j)=B(j,j+1)=0 and B(j,j) and
B(j+1,j+1) will be positive.)
If JOB='E', then at each iteration, the same transformations are computed, but they are only applied to those parts of A and B which are needed to compute ALPHAR, ALPHAI, and BETAR.
If JOB='S' and COMPQ and COMPZ are 'V' or 'I', then the orthogonal transformations used to reduce (A,B) are accumulated into the arrays Q and Z s.t.:
(in) A(in) Z(in)* = Q(out) A(out) Z(out)*
(in) B(in) Z(in)* = Q(out) B(out) Z(out)*
Ref: C.B. Moler & G.W. Stewart, ``An Algorithm for Generalized Matrixigenvalue Problems'', SIAM J. Numer. Anal., 10(1973),p. 241--256.
= 'E': compute only ALPHAR, ALPHAI, and BETA. A and B will not necessarily be put into generalized Schur form. = 'S': put A and B into generalized Schur form, as well as computing ALPHAR, ALPHAI, and BETA.
= 'N': do not modify Q.
= 'V': multiply the array Q on the right by the transpose of the orthogonal tranformation that is applied to the left side of A and B to reduce them to Schur form. = 'I': like COMPQ ='V', except that Q will be initialized to the identity first.
= 'N': do not modify Z.
= 'V': multiply the array Z on the right by the orthogonal tranformation that is applied to the right side of A and B to reduce them to Schur form. = 'I': like COMPZ ='V', except that Z will be initialized to the identity first.
A(j+1,j) is non-zero,
then B(j+1,j) =B(j,j+1) =0 and B(j,j) and B(j+1,j+1) will be
positive.)
If JOB ='S', then on exit A and B will have been
simultaneously reduced to Schur form.
If JOB ='E', then on exit B will have been destroyed.
Elements corresponding to diagonal blocks of A will be
correct, but the off-diagonal portion will be meaningless.
ALPHAR(1:N) will be set to real parts of the diagonal
elements of A that would result from reducing A and B to
Schur form and then further reducing them both to triangular
form using unitary transformations s.t. the diagonal of B
was non-negative real. Thus, if A(j,j) is in a 1-by-1 block
(i.e., A(j+1,j) =A(j,j+1) =0), then ALPHAR(j) =A(j,j).
Note that the (real or complex) values
(ALPHAR(j) + i*ALPHAI(j))/BETA(j), j =1,...,N, are the
generalized eigenvalues of the matrix pencil A - wB.
ALPHAI(1:N) will be set to imaginary parts of the diagonal
elements of A that would result from reducing A and B to
Schur form and then further reducing them both to triangular
form using unitary transformations s.t. the diagonal of B
was non-negative real. Thus, if A(j,j) is in a 1-by-1 block
(i.e., A(j+1,j) =A(j,j+1) =0), then ALPHAR(j) =0.
Note that the (real or complex) values
(ALPHAR(j) + i*ALPHAI(j))/BETA(j), j =1,...,N, are the
generalized eigenvalues of the matrix pencil A - wB.
BETA(1:N) will be set to the (real) diagonal elements of B
that would result from reducing A and B to Schur form and
then further reducing them both to triangular form using
unitary transformations s.t. the diagonal of B was
non-negative real. Thus, if A(j,j) is in a 1-by-1 block
(i.e., A(j+1,j) =A(j,j+1) =0), then BETA(j) =B(j,j).
Note that the (real or complex) values
(ALPHAR(j) + i*ALPHAI(j))/BETA(j), j =1,...,N, are the
generalized eigenvalues of the matrix pencil A - wB.
(Note that BETA(1:N) will always be non-negative, and no
BETAI is necessary.)
WORK(1) returns the optimal LWORK.
If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA.
= 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal value
= 1,...,N: the QZ iteration did not converge. (A,B) is not in Schur form, but ALPHAR(i), ALPHAI(i), and BETA(i), i =INFO+1,...,N should be correct. = N+1,...,2*N: the shift calculation failed. (A,B) is not in Schur form, but ALPHAR(i), ALPHAI(i), and BETA(i), i =INFO-N+1,...,N should be correct. > 2*N: various "impossible" errors.
Iteration counters:
JITER -- counts iterations.
IITER -- counts iterations run since ILAST was last
changed. This is therefore reset only when a 1-by-1 or
2-by-2 block deflates off the bottom.