Generalized Eigenvalues and Eigenvectors for Hermitian 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 ZHEGV	(ITYPE, JOBZ, UPLO, N, ZA, LDA, ZB, LDB, DW, ZWORK, LDWORK, DWORK2, INFO)
CALL CHEGV	(ITYPE, JOBZ, UPLO, N, CA, LDA, CB, LDB, SW, CWORK, LDWORK, SWORK2, INFO)
void zhegv	(int itype, char jobz, char uplo, int n, doublecomplex *za, int lda, doublecomplex *zb, int ldb, double *dw, int *info)
void chegv	(int itype, char jobz, char uplo, int n, complex *ca, int lda, complex *cb, int ldb, float *sw, int *info)

Arguments

ITYPE	Indicates the problem to be solved. Legal values for ITYPE are shown below. Any values not shown below are illegal.
	1	Ax = Bx
	2	ABx = x
	3	BAx = x
JOBZ	Indicates whether to compute eigenvalues only or to compute both eigenvalues and eigenvectors. The legal values for JOBZ are listed below. Any values not listed below are illegal.
	'N' or 'n'	Compute eigenvalues only.
	'V' or 'v'	Compute both eigenvalues and eigenvectors.
UPLO	Indicates whether xA and xB contain the upper or lower triangle of the matrix. The legal values for UPLO are listed below. Any values not listed below are illegal.
	'U' or 'u'	xA and xB both contain the upper triangle.
	'L' or 'l'	xA and xB both contain the lower triangle.
N	Order of the matrix A. N 0.
xA	On entry, the upper or lower triangle of the matrix A. On exit, if JOBZ = 'V' or 'v' then A contains the matrix Z of eigenvectors normalized as shown below:
	ITYPE	Eigenvectors stored in A
	1 or 2	ZH × B × Z = I
	3	ZH × B-1 × Z = I
	If JOBZ = 'N' or 'n' and UPLO = 'U' or 'u' then the upper triangle of A including the diagonal is overwritten. If JOBZ = 'N' or 'n' and UPLO = 'L' or 'l' then the lower triangle of A including the diagonal 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 upper or lower triangle of the matrix B. On exit, if UPLO = 'U' or 'u' then B contains the matrix U from the Cholesky factorization UHU of B. If UPLO = 'L' or 'l' then B contains the matrix L from the Cholesky factorization LLH of B.
LDB	Leading dimension of the array B as specified in a dimension or type statement. LDB max(1, N).
xW	On exit, W(1:NFOUND) contains the computed eigenvalues in ascending order.
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, 2 × N - 1).
xWORK2	Scratch array with a dimension of max(1, 3 × N - 2).
INFO	On exit:
	INFO = 0	Subroutine completed normally.
	INFO < 0	The ith argument, where i = |INFO|, had an illegal value.
	0 < INFO N	Convergence failure
	INFO > N	The leading minor of order i of B, where i = INFO-N, is not positive definite.

Sample Program

PROGRAM TEST

IMPLICIT NONE

C

DOUBLE PRECISION ZERO

INTEGER ITYPE, LDA, LDB, LDWORK, LDWRK2, N

PARAMETER (ITYPE = 1)

PARAMETER (N = 3)

PARAMETER (LDA = N)

PARAMETER (LDB = N)

PARAMETER (LDWORK = 2 * N - 1)

PARAMETER (LDWRK2 = 3 * N - 2)

PARAMETER (ZERO = 0.0D0)

C

DOUBLE PRECISION EVALS(N), WORK2(LDWRK2)

COMPLEX*16 A(LDA,N), B(LDB,N), WORK(LDWORK)

INTEGER ICOL, INFO, IROW

C

EXTERNAL ZHEGV

INTRINSIC ABS, CONJG

C

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

C shown below.

C

C 8 2-2i 2-2i 4 1-i 1-i

C A = 2+2i 16 4-4i B = 1+i 8 2-2i

C 2+2i 4+4i 32 1+i 2+2i 16

C

DATA A / (8.0D0,8D8), (8D8,8D8), (8D8,8D8),

$ (2.0D0,-2.0D0), (1.6D1,8D8), (8D8,8D8),

$ (2.0D0,-2.0D0), (4.0D0,-4.0D0), (3.2D1,8D8) /

DATA B / (4.0D0,8D8), (8D8,8D8), (8D8,8D8),

$ (1.0D0,-1.0D0), (8.0D0,8D8), (8D8,8D8),

$ (1.0D0,-1.0D0), (2.0D0,-2.0D0), (1.6D1,8D8) /

C

PRINT 1000

DO 100, IROW = 1, N

PRINT 1010, (CONJG(A(ICOL,IROW)), ICOL = 1, IROW - 1),

$ CMPLX(DBLE(A(IROW,IROW)), ZERO),

$ (A(IROW,ICOL), ICOL = IROW + 1, N)

100 CONTINUE

PRINT 1020

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

PRINT 1030

DO 110, IROW = 1, N

PRINT 1010, (CONJG(B(ICOL,IROW)), ICOL = 1, IROW - 1),

$ CMPLX(DBLE(A(IROW,IROW)), ZERO),

$ (B(IROW,ICOL), ICOL = IROW + 1, N)

110 CONTINUE

PRINT 1040

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

C

CALL ZHEGV (ITYPE, 'VALUES AND EIGENVECTORS',

$ 'UPPER TRIANGLE OF A STORED', N, A, LDA, B, LDB,

$ EVALS, WORK, LDWORK, WORK2, INFO)

IF (INFO .LT. 0) THEN

PRINT 1050, ABS(INFO)

STOP 1

ELSE IF (INFO .GT. N) THEN

PRINT 1060, INFO

STOP 2

ELSE IF (INFO .GT. 0) THEN

PRINT 1070, INFO

STOP 3

END IF

PRINT 1080

DO 120, IROW = 1, N

PRINT 1090, EVALS(IROW), (A(IROW,ICOL), ICOL = 1, N)

120 CONTINUE

C

1000 FORMAT (1X, 'A in full form:')

1010 FORMAT (3(3X, '(', F5.2, ',', F5.2, ')'))

1020 FORMAT (/1X, 'A in Hermitian form: (* in unused elements)')

1030 FORMAT (/1X, 'B in full form:')

1040 FORMAT (/1X, 'B in Hermitian form: (* in unused elements)')

1050 FORMAT (1X, 'Illegal argument to ZHEGV, argument #', I2)

1060 FORMAT (1X, 'B is not positive definite, INFO = ', I4)

1070 FORMAT (1X, 'Convergence failure, INFO = ', I4)

1080 FORMAT (/1X, 'Eigenvalue', 20X, 'Eigenvector**T')

1090 FORMAT (1X, F8.3, 5X, '[', 3('(', F5.3, ',', F5.3, ')', 2X),

$ ']')

C

END

Sample Output

A in full form:

( 8.00, 0.00) ( 2.00,-2.00) ( 2.00,-2.00)

( 2.00, 2.00) (16.00, 0.00) ( 4.00,-4.00)

( 2.00, 2.00) ( 4.00, 4.00) (32.00, 0.00)

A in Hermitian form: (* in unused elements)

( 8.00,*****) ( 2.00,-2.00) ( 2.00,-2.00)

(*****,*****) (16.00,*****) ( 4.00,-4.00)

(*****,*****) (*****,*****) (32.00,*****)

B in full form:

( 8.00, 0.00) ( 1.00,-1.00) ( 1.00,-1.00)

( 1.00, 1.00) (16.00, 0.00) ( 2.00,-2.00)

( 1.00, 1.00) ( 2.00, 2.00) (32.00, 0.00)

B in Hermitian form: (* in unused elements)

( 4.00,*****) ( 1.00,-1.00) ( 1.00,-1.00)

(*****,*****) ( 8.00,*****) ( 2.00,-2.00)

(*****,*****) (*****,*****) (16.00,*****)

Eigenvalue Eigenvector**T

2.000 [(0.500,0.000) (-.091,0.091) (-.070,0.035) ]

2.000 [(0.000,0.000) (0.365,0.000) (-.052,0.070) ]

2.000 [(0.000,0.000) (0.000,0.000) (0.261,0.000) ]