Generalized Eigenvalues, Schur Form, and Schur Vectors of General Matrices

A generalized eigenvalue is a ratio = / such that AX = BX or
det(A-B) = 0 for general matrices A, B, and X. is typically represented as a ratio rather than as a scalar because there are reasonable interpretations for = 0, = 0, and for = = 0. A right generalized eigenvector x corresponding to a generalized eigenvalue is defined by (A - B)x = 0. A left generalized eigenvector x corresponding to a generalized eigenvalue is defined by
(A - B)^Hx = 0. A good reference for generalized eigenproblems is the book Matrix Computations, 2nd. ed. by Golub and van Loan (1989, The Johns Hopkins University Press).

Calling Sequence

CALL DGEGS	(JOBVSL, JOBVSR, N, DA, LDA, DB, LDB, DALPHAR, DALPHAI, DBETA, DVSL, LDVSL, DVSR, LDVSR, DWORK, LDWORK, INFO)
CALL SGEGS	(JOBVSL, JOBVSR, N, SA, LDA, SB, LDB, SALPHAR, SALPHAI, SBETA, SVSL, LDVSL, SVSR, LDVSR, SWORK, LDWORK, INFO)
CALL ZGEGS	(JOBVSL, JOBVSR, N, ZA, LDA, ZB, LDB, ZALPHA, ZBETA, ZVSL, LDVSL, ZVSR, LDVSR, ZWORK, LDWORK, DWORK2, INFO)
CALL CGEGS	(JOBVSL, JOBVSR, N, CA, LDA, CB, LDB, CALPHA, CBETA, CVSL, LDVSL, CVSR, LDVSR, CWORK, LDWORK, SWORK2, INFO)
void dgegs	(char jobvsl, char jobvsr, int n, double *da, int lda, double *db, int ldb, double *dalphar, double *dalphai, double *dbeta, double *dvsl, int ldvsl, double *dvsr, int ldvsr, int *info)
void sgegs	(char jobvsl, char jobvsr, int n, float *sa, int lda, float *sb, int ldb, float *salphar, float *dalphai, float *sbeta, float *svsl, int ldvsl, float *svsr, int ldvsr, int *info)
void zgegs	(char jobvsl, char jobvsr, int n, doublecomplex *za, int lda, doublecomplex *zb, int ldb, doublecomplex *zalpha, doublecomplex *zbeta, doublecomplex *zvsl, int ldvsl, doublecomplex *zvsr, int ldvsr, int *info)
void cgegs	(char jobvsl, char jobvsr, int n, complex *ca, int lda, complex *cb, int ldb, complex *calpha, complex *cbeta, complex *cvsl, int ldvsl, complex *cvsr, int ldvsr, int *info)

Arguments

JOBVSL	Indicates whether or not to compute the left Schur vectors. The legal values of JOBVSL are shown below. Any values not shown below are illegal.
	'N' or 'n'	Do not compute the left Schur vectors.
	'V' or 'v'	Compute the left Schur vectors.
JOBVSR	Indicates whether or not to compute the right Schur vectors. The legal values of JOBVSR are shown below. Any values not shown below are illegal.
	'N' or 'n'	Do not compute the right Schur vectors.
	'V' or 'v'	Compute the right Schur vectors.
N	Number of rows and columns in the matrices A, B, VSL, and VSR. N 0.
xA	On entry, the first of the pair of matrices whose generalized eigenvalues (and Schur vectors) are to be computed. On exit, the generalized Schur form of A. Note: to avoid overflow, the Frobenius norm of A should be less than the overflow threshold.
LDA	Leading dimension of the array A as specified in a dimension or type statement. LDA max(1, N).
xB	On entry, the second of the pair of matrices whose generalized eigenvalues to be computed. On exit, the generalized Schur form of B. Note: to avoid overflow, the Frobenius norm of B should be less than the overflow threshold.
LDB	Leading dimension of the array B as specified in a dimension or type statement. LDB max(1, N).
xALPHAR, xALPHAI, xBETA
	(For real subroutines) On exit, (ALPHAR(j) + i × ALPHAI(j)) / BETA(j) for 1 j N contain the generalized eigenvalues. ALPHAR(j) + i × ALPHAI(j) and BETA(j) for 1 j N are the diagonals of the complex Schur form (A,B) that would result if the 2×2 diagonal blocks of the real Schur form of (A,B) were further reduced to triangular form using 2×2 complex unitary transformations. If ALPHAI(j) = 0, then the jth eigenvalue is real; if positive, then the jth and (j+1)th eigenvalues are a complex conjugate pair, with ALPHAI(j+1) negative.
	Note: The quotients ALPHAR(j) / BETA(j) and ALPHAI(j) / BETA(j) may easily over- or underflow, and any or all of ALPHAR(j), ALPHAI(j), BETA(j) may be zero. Thus, the user should avoid naively computing the ratio / . However, ALPHAR and ALPHAI will always be less than and usually comparable with ||A|| in magnitude, and BETA always less than and usually comparable with ||B||.
xALPHA, xBETA (For complex subroutines)
	On exit, ALPHA(j) / BETA(j) for 1 j N are the generalized eigenvalues. ALPHA(j) and BETA(j) for 1 j N are the diagonals of the complex Schur form (A,B) computed by xGEGS. BETA(j) is non-negative and real.
	Note: The quotients ALPHA(j) / BETA(j) may easily over- or underflow, and either or both ALPHA(j) or BETA(j) may be zero. Thus, the user should avoid naively computing the ratio / . However, ALPHA will always be less than and usually comparable with ||A|| in magnitude, and BETA always less than and usually comparable with ||B||.
xVSL	On exit, if JOBVSL = 'V' or 'v' then VSL will contain the left Schur vectors. If JOBVSL = 'N' or 'n' then VSL is not used.
LDVSL	Leading dimension of the array VSL as specified in a dimension or type statement. LDVSL 1. If JOBVSL = 'V' or 'v' then LDVSL max(1, N).
xVSR	On exit, if JOBVSR = 'V' or 'v' then VSR will contain the right Schur vectors. If JOBVSR = 'N' or 'n' then VSR is not used.
LDVSR	Leading dimension of the array VSR as specified in a dimension or type statement. LDVSR 1. If JOBVSR = 'V' or 'v', LDVSR max(1, N).
xWORK	Scratch array with a dimension of LDWORK.
LDWORK	Leading dimension of the array WORK as specified in a dimension or type statement. LDWORK max(1, 4 × N) for real subroutines or max(1, 2 × N) for complex subroutines.
xWORK2	Scratch array with a dimension of 3 × N for complex subroutines.
INFO	On exit:
	INFO = 0	Subroutine completed normally.
	INFO < 0	The ith argument, where i = |INFO|, had an illegal value.
	1 INFO N	The QZ iteration failed. (A,B) are not in Schur form, but ALPHAR(j), ALPHAI(j), and BETA(j) -- or ALPHA(j) and BETA(j) -- will be correct for INFO+1 j N.
	INFO > N	May indicate an internal LAPACK error.

Sample Program

PROGRAM TEST

IMPLICIT NONE

C

INTEGER LDA, LDB, LDWORK, M, N

PARAMETER (M = 3)

PARAMETER (N = 3)

PARAMETER (LDA = M)

PARAMETER (LDB = M)

PARAMETER (LDWORK = 4 * N)

C

DOUBLE PRECISION A(LDA,N), ALPHAI(N), ALPHAR(N), B(LDB,N)

DOUBLE PRECISION BETA(N), TEMP, WORK(LDWORK)

INTEGER ICOL, INFO, IROW

C

EXTERNAL DGEGS

INTRINSIC ABS, DCMPLX

C

C Initialize the arrays A and B to store the matrices A and B

C shown below.

C

C 1 1 1 1 1 1

C A = 2 1 1 B = 1 0 0

C 1 3 1 0 1 0

C

DATA A / 1.0D0, 2.0D0, 1.0D0,

$ 1.0D0, 1.0D0, 3.0D0,

$ 1.0D0, 1.0D0, 1.0D0 /

DATA B / 1.0D0, 1.0D0, 0.0D0,

$ 1.0D0, 0.0D0, 1.0D0,

$ 1.0D0, 0.0D0, 0.0D0 /

C

PRINT 1000

PRINT 1010, ((A(IROW,ICOL), ICOL = 1, N), IROW = 1, N)

PRINT 1020

PRINT 1010, ((B(IROW,ICOL), ICOL = 1, N), IROW = 1, N)

C

CALL DGEGS ('NO LEFT SCHUR VECTORS',

$ 'NO RIGHT SCHUR VECTORS', N, A, LDA, B, LDB,

$ ALPHAR, ALPHAI, BETA, TEMP, 1, TEMP, 1,

$ WORK, LDWORK, INFO)

IF (INFO .NE. 0) THEN

IF (INFO .LT. 0) THEN

PRINT 1030, ABS(INFO)

STOP 1

ELSE

PRINT 1040, INFO

STOP 2

END IF

END IF

PRINT 1050

PRINT 1060, (ALPHAR(IROW), ALPHAI(IROW), BETA(IROW),

$ IROW = 1, M)

PRINT 1070

PRINT 1080, (DCMPLX(ALPHAR(IROW), ALPHAI(IROW)) /

$ BETA(IROW), IROW = 1, M)

C

1000 FORMAT (1X, 'A:')

1010 FORMAT (3(3X, F4.1))

1020 FORMAT (/1X, 'B:')

1030 FORMAT (1X, 'Illegal argument to DGEGS, argument #', I2)

1040 FORMAT (1X, 'Internal failure in DGEGS, INFO = ', I4)

1050 FORMAT (/1X, T4, 'ALPHA', T20, 'BETA')

1060 FORMAT (1X, '(', F5.2, ',', F5.2, ')', T19, F6.2)

1070 FORMAT (/1X, 'ALPHA / BETA')

1080 FORMAT (1X, '(', F5.2, ',', F5.2, ')')

C

END

Sample Output

A:
1.0 1.0 1.0
2.0 1.0 1.0
1.0 3.0 1.0
B:
1.0 1.0 1.0
1.0 0.0 0.0
0.0 1.0 0.0
ALPHA BETA
( 0.71, 0.00) 0.71
( 0.82, 0.00) 0.82
( 3.46, 0.00) 1.73
ALPHA / BETA
( 1.00, 0.00)
( 1.00, 0.00)
( 2.00, 0.00)