shgeqz - 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 addi- tion, 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);
Oracle Solaris Studio Performance Library shgeqz(3P)
NAME
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 addi-
tion, the pair A,B may be reduced to generalized Schur form
SYNOPSIS
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);
PURPOSE
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 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-triangu-
lar 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 com-
puted, 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)*
Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrixi-
genvalue Problems", SIAM J. Numer. Anal., 10(1973),p. 241--256.
ARGUMENTS
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
destroyed. 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 corre-
sponding 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 posi-
tive.) 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 diagonal ele-
ments 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 pen-
cil 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 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 pen-
cil A - wB.
BETA (output)
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 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 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 pen-
cil A - wB. (Note that BETA(1:N) will always be non-nega-
tive, 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 transformations 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 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.
FURTHER DETAILS
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.
7 Nov 2015 shgeqz(3P)