stgevc - compute some or all of the right and/or left generalized eigenvectors of a pair of real upper triangular matrices (A,B)
SUBROUTINE STGEVC( SIDE, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, * LDVL, VR, LDVR, MM, M, WORK, INFO) CHARACTER * 1 SIDE, HOWMNY INTEGER N, LDA, LDB, LDVL, LDVR, MM, M, INFO LOGICAL SELECT(*) REAL A(LDA,*), B(LDB,*), VL(LDVL,*), VR(LDVR,*), WORK(*)
SUBROUTINE STGEVC_64( SIDE, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, * LDVL, VR, LDVR, MM, M, WORK, INFO) CHARACTER * 1 SIDE, HOWMNY INTEGER*8 N, LDA, LDB, LDVL, LDVR, MM, M, INFO LOGICAL*8 SELECT(*) REAL A(LDA,*), B(LDB,*), VL(LDVL,*), VR(LDVR,*), WORK(*)
SUBROUTINE TGEVC( SIDE, HOWMNY, SELECT, [N], A, [LDA], B, [LDB], VL, * [LDVL], VR, [LDVR], MM, M, [WORK], [INFO]) CHARACTER(LEN=1) :: SIDE, HOWMNY INTEGER :: N, LDA, LDB, LDVL, LDVR, MM, M, INFO LOGICAL, DIMENSION(:) :: SELECT REAL, DIMENSION(:) :: WORK REAL, DIMENSION(:,:) :: A, B, VL, VR
SUBROUTINE TGEVC_64( SIDE, HOWMNY, SELECT, [N], A, [LDA], B, [LDB], * VL, [LDVL], VR, [LDVR], MM, M, [WORK], [INFO]) CHARACTER(LEN=1) :: SIDE, HOWMNY INTEGER(8) :: N, LDA, LDB, LDVL, LDVR, MM, M, INFO LOGICAL(8), DIMENSION(:) :: SELECT REAL, DIMENSION(:) :: WORK REAL, DIMENSION(:,:) :: A, B, VL, VR
#include <sunperf.h>
void stgevc(char side, char howmny, logical *select, int n, float *a, int lda, float *b, int ldb, float *vl, int ldvl, float *vr, int ldvr, int mm, int *m, int *info);
void stgevc_64(char side, char howmny, logical *select, long n, float *a, long lda, float *b, long ldb, float *vl, long ldvl, float *vr, long ldvr, long mm, long *m, long *info);
stgevc computes some or all of the right and/or left generalized eigenvectors of a pair of real upper triangular matrices (A,B).
The right generalized eigenvector x and the left generalized eigenvector y of (A,B) corresponding to a generalized eigenvalue w are defined by:
(A - wB) * x = 0 and y**H * (A - wB) = 0
where y**H denotes the conjugate tranpose of y.
If an eigenvalue w is determined by zero diagonal elements of both A and B, a unit vector is returned as the corresponding eigenvector.
If all eigenvectors are requested, the routine may either return the matrices X and/or Y of right or left eigenvectors of (A,B), or the products Z*X and/or Q*Y, where Z and Q are input orthogonal matrices. If (A,B) was obtained from the generalized real-Schur factorization of an original pair of matrices
(A0,B0) = (Q*A*Z**H,Q*B*Z**H),
then Z*X and Q*Y are the matrices of right or left eigenvectors of A.
A must be block upper triangular, with 1-by-1 and 2-by-2 diagonal blocks. Corresponding to each 2-by-2 diagonal block is a complex conjugate pair of eigenvalues and eigenvectors; only one
eigenvector of the pair is computed, namely the one corresponding to the eigenvalue with positive imaginary part.
= 'R': compute right eigenvectors only;
= 'L': compute left eigenvectors only;
= 'B': compute both right and left eigenvectors.
= 'A': compute all right and/or left eigenvectors;
= 'B': compute all right and/or left eigenvectors, and backtransform them using the input matrices supplied in VR and/or VL; = 'S': compute selected right and/or left eigenvectors, specified by the logical array SELECT.
SELECT(j) must be set to .TRUE. To select
the complex eigenvector corresponding to a complex conjugate
pair w(j) and w(j+1), either SELECT(j) or SELECT(j+1) must
be set to .TRUE..
A complex eigenvector corresponding to a complex eigenvalue is stored in two consecutive columns, the first holding the real part, and the second the imaginary part.
max(1,N) if SIDE = 'L' or 'B'; LDVL > = 1 otherwise.
A complex eigenvector corresponding to a complex eigenvalue is stored in two consecutive columns, the first holding the real part and the second the imaginary part.
max(1,N) if SIDE = 'R' or 'B'; LDVR > = 1 otherwise.
dimension(6*N)
= 0: successful exit.
< 0: if INFO = -i, the i-th argument had an illegal value.
> 0: the 2-by-2 block (INFO:INFO+1) does not have a complex eigenvalue.
Allocation of workspace:
---------- -- ---------
WORK( j ) = 1-norm of j-th column of A, above the diagonal WORK( N+j ) = 1-norm of j-th column of B, above the diagonal WORK( 2*N+1:3*N ) = real part of eigenvector
WORK( 3*N+1:4*N ) = imaginary part of eigenvector
WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector
Rowwise vs. columnwise solution methods:
------- -- ---------- -------- -------
Finding a generalized eigenvector consists basically of solving the singular triangular system
(A - w B) x = 0 (for right) or: (A - w B)**H y = 0 (for left)
Consider finding the i-th right eigenvector (assume all eigenvalues are real). The equation to be solved is: i
0 = sum C(j,k) v(k) = sum C(j,k) v(k) for j = i,. . .,1
k =j k =j
where C = (A - w B) (The components v(i+1:n) are 0.)
The ``rowwise'' method is:
(1) v(i) : = 1
for j = i-1,. . .,1:
i
(2) compute s = - sum C(j,k) v(k) and
k =j+1
(3) v(j) : = s / C(j,j)
Step 2 is sometimes called the ``dot product'' step, since it is an inner product between the j-th row and the portion of the eigenvector that has been computed so far.
The ``columnwise'' method consists basically in doing the sums
for all the rows in parallel. As each v(j) is computed, the
contribution of v(j) times the j-th column of C is added to the
partial sums. Since FORTRAN arrays are stored columnwise, this has
the advantage that at each step, the elements of C that are accessed
are adjacent to one another, whereas with the rowwise method, the
elements accessed at a step are spaced LDA (and LDB) words apart.
When finding left eigenvectors, the matrix in question is the transpose of the one in storage, so the rowwise method then actually accesses columns of A and B at each step, and so is the preferred method.