Contents
shgeqz - 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 SHGEQZ(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
REAL A(LDA,*), B(LDB,*), ALPHAR(*), ALPHAI(*), BETA(*),
Q(LDQ,*), Z(LDZ,*), WORK(*)
SUBROUTINE SHGEQZ_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
REAL A(LDA,*), B(LDB,*), ALPHAR(*), ALPHAI(*), BETA(*),
Q(LDQ,*), Z(LDZ,*), WORK(*)
F95 INTERFACE
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, DIMENSION(:) :: ALPHAR, ALPHAI, BETA, WORK
REAL, 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, DIMENSION(:) :: ALPHAR, ALPHAI, BETA, WORK
REAL, DIMENSION(:,:) :: A, B, Q, Z
C INTERFACE
#include <sunperf.h>
void shgeqz(char job, char compq, char compz, int n, int
ilo, int ihi, float *a, int lda, float *b, int
ldb, float *alphar, float *alphai, float *beta,
float *q, int ldq, float *z, int ldz, int *info);
void shgeqz_64(char job, char compq, char compz, long n,
long ilo, long ihi, float *a, long lda, float *b,
long ldb, float *alphar, float *alphai, float
*beta, float *q, long ldq, float *z, long ldz,
long *info);
shgeqz 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 diag-
onal 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 posi-
tive 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 accumu-
lated into the arrays Q and Z s.t.:
(in) A(in) Z(in)* = Q(out) A(out) Z(out)*
Ref: C.B. Moler & G.W. Stewart, "An Algorithm for General-
ized Matrixigenvalue Problems", SIAM J. Numer. Anal.,
10(1973),p. 241--256.
JOB (input)
= '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.
COMPQ (input)
= '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.
COMPZ (input)
= '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.
N (input) The order of the matrices A, B, Q, and Z. N >= 0.
ILO (input)
It is assumed that A is already upper triangular
in rows and columns 1:ILO-1 and IHI+1:N. 1 <= ILO
<= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
IHI (input)
See the description of ILO.
A (input) On entry, the N-by-N upper Hessenberg matrix A.
Elements below the subdiagonal must be zero. If
JOB='S', then on exit A and B will have been
simultaneously reduced to generalized Schur form.
If JOB='E', then on exit A will have been des-
troyed. The diagonal blocks will be correct, but
the off-diagonal portion will be meaningless.
LDA (input)
The leading dimension of the array A. LDA >= max(
1, N ).
B (input) On entry, the N-by-N upper triangular matrix B.
Elements below the diagonal must be zero. 2-by-2
blocks in B corresponding to 2-by-2 blocks in A
will be reduced to positive diagonal form. (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='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.
LDB (input)
The leading dimension of the array B. LDB >= max(
1, N ).
ALPHAR (output)
ALPHAR(1:N) will be set to real parts of the diag-
onal 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 (output)
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 uni-
tary 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 (output)
BETA(1:N) will be set to the (real) diagonal ele-
ments 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 general-
ized eigenvalues of the matrix pencil A - wB.
(Note that BETA(1:N) will always be non-negative,
and no BETAI is necessary.)
Q (input/output)
If COMPQ='N', then Q will not be referenced. If
COMPQ='V' or 'I', then the transpose of the
orthogonal transformations which are applied to A
and B on the left will be applied to the array Q
on the right.
LDQ (input)
The leading dimension of the array Q. LDQ >= 1.
If COMPQ='V' or 'I', then LDQ >= N.
Z (input/output)
If COMPZ='N', then Z will not be referenced. If
COMPZ='V' or 'I', then the orthogonal transforma-
tions which are applied to A and B on the right
will be applied to the array Z on the right.
LDZ (input)
The leading dimension of the array Z. LDZ >= 1.
If COMPZ='V' or 'I', then LDZ >= N.
WORK (workspace)
On exit, if INFO >= 0, WORK(1) returns the optimal
LWORK.
LWORK (input)
The dimension of the array WORK. LWORK >=
max(1,N).
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.
INFO (output)
= 0: successful exit
< 0: if INFO = -i, the i-th argument had an ille-
gal 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.