Generalized Eigenvalues and Eigenvectors for 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 DGEGV	(JOBVL, JOBVR, N, DA, LDA, DB, LDB, DALPHAR, DALPHAI, DBETA, DVL, LDVL, DVR, LDVR, DWORK, LDWORK, INFO)
CALL SGEGV	(JOBVL, JOBVR, N, SA, LDA, SB, LDB, SALPHAR, SALPHAI, SBETA, SVL, LDVL, SVR, LDVR, SWORK, LDWORK, INFO)
CALL ZGEGV	(JOBVL, JOBVR, N, ZA, LDA, ZB, LDB, ZALPHA, ZBETA, ZVL, LDVL, ZVR, LDVR, ZWORK, LDWORK, DWORK2, INFO)
CALL CGEGV	(JOBVL, JOBVR, N, CA, LDA, CB, LDB, CALPHA, CBETA, CVL, LDVL, CVR, LDVR, CWORK, LDWORK, SWORK2, INFO)
void dgegv	(char jobvl, char jobvr, int n, double *da, int lda, double *db, int ldb, double *dalphar, double *dalphai, double *dbeta, double *dvl, int ldvl, double *dvr, int ldvr, int *info)
void sgegv	(char jobvl, char jobvr, int n, float *sa, int lda, float *sb, int ldb, float *salphar, float *dalphai, float *sbeta, float *svl, int ldvl, float *svr, int ldvr, int *info)
void zgegv	(char jobvl, char jobvr, int n, doublecomplex *za, int lda, doublecomplex *zb, int ldb, doublecomplex *zalpha, doublecomplex *zbeta, doublecomplex *zvl, int ldvl, doublecomplex *zvr, int ldvr, int *info)
void cgegv	(char jobvl, char jobvr, int n, complex *ca, int lda, complex *cb, int ldb, complex *calpha, complex *cbeta, complex *cvl, int ldvl, complex *cvr, int ldvr, int *info)

Arguments

JOBVL	Indicates whether the left generalized eigenvectors will be computed. Legal values for JOBVL are shown below. Any values not shown below are illegal.
	'N' or 'n'	Do not compute the left generalized eigenvectors.
	'V' or 'v'	Compute the left generalized eigenvectors.
JOBVR	Indicates whether the right generalized eigenvectors will be computed. Legal values for JOBVR are shown below. Any values not shown below are illegal.
	'N' or 'n'	Do not compute the right generalized eigenvectors.
	'V' or 'v'	Compute the right generalized eigenvectors.
N	Number of rows and columns in the matrices A, B, VL, and VR. N 0.
xA	On entry, the first of the pair of matrices whose generalized eigenvalues are to be computed. On exit, if JOBVL or JOBVR = 'V' or 'v', then A contains the real Schur form of the balanced form of A. Otherwise A is overwritten.
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 are to be computed. On exit, if JOBVL or JOBVR = 'V' or 'v', then B contains the real Schur form of the balanced form of B. Otherwise B is overwritten.
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 will contain the generalized eigenvalues. 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 be always 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) as 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||.
xVL	On exit, if JOBVL = 'V' or 'v' then VL contains the left generalized eigenvectors. VL is not used if JOBVL = 'N' or 'n'. Real eigenvectors occupy one column. Complex eigenvectors occupy two columns, the first for the real part and the second for the imaginary part. Complex eigenvectors correspond to an eigenvalue with positive imaginary part. With the exception of those eigenvalues with = = 0, each eigenvector will be scaled so the largest component will have abs(real part) + abs(imag. part) = 1. For eigenvalues with = = 0, a zero vector will be returned as the corresponding eigenvector.
LDVL	Leading dimension of the array VL as specified in a dimension or type statement. LDVL 1. If JOBVL = 'V' or 'v' then LDVL max(1, N).
xVR	On exit, if JOBVR = 'V' or 'v' then VR contains the right generalized eigenvectors. VR is not used if JOBVR = 'N' or 'n'. Real eigenvectors occupy one column. Complex eigenvectors occupy two columns, the first for the real part and the second for the imaginary part. Complex eigenvectors correspond to an eigenvalue with positive imaginary part. With the exception of those eigenvalues with = = 0, each eigenvector will be scaled so the largest component will have abs(real part) + abs(imag. part) = 1. For eigenvalues with = = 0, a zero vector will be returned as the corresponding eigenvector.
LDVR	Leading dimension of the array VR as specified in a dimension or type statement. LDVR 1. If JOBVR = 'V' or 'v' then LDVR 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, 8 × N) for real subroutines or max(1, 2 × N) for complex subroutines.
xWORK2	Scratch array with a dimension of 8 × N for complex subroutines.
INFO	On exit:
	INFO = 0	Subroutine completed normally.
	INFO < 0	The ith argument, where i = |INFO|, had an illegal value. The QZ iteration failed. No eigenvectors have been calculated, but ALPHAR(j), ALPHAI(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 = 8 * N)

C

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

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

INTEGER ICOL, INFO, IROW

C

EXTERNAL DGEGV

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

C Print the initial value of the arrays.

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

C Compute and print the eigenvalues.

C

CALL DGEGV ('NO LEFT EIGENVECTORS',

$ 'NO RIGHT EIGENVECTORS', N, A, LDA, B, LDB,

$ ALPHAR, ALPHAI, BETA, SCRATCH, 1, SCRATCH, 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 DGEGV, argument #', I2)

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

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

1060 FORMAT (1X, '(', F5.2, ',', F5.2, ')', T19, F5.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)