1 SUBROUTINE ZHEGV( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
2 $ LWORK, RWORK, INFO )
3 *
4 * -- LAPACK driver routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 CHARACTER JOBZ, UPLO
11 INTEGER INFO, ITYPE, LDA, LDB, LWORK, N
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION RWORK( * ), W( * )
15 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZHEGV computes all the eigenvalues, and optionally, the eigenvectors
22 * of a complex generalized Hermitian-definite eigenproblem, of the form
23 * A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x.
24 * Here A and B are assumed to be Hermitian and B is also
25 * positive definite.
26 *
27 * Arguments
28 * =========
29 *
30 * ITYPE (input) INTEGER
31 * Specifies the problem type to be solved:
32 * = 1: A*x = (lambda)*B*x
33 * = 2: A*B*x = (lambda)*x
34 * = 3: B*A*x = (lambda)*x
35 *
36 * JOBZ (input) CHARACTER*1
37 * = 'N': Compute eigenvalues only;
38 * = 'V': Compute eigenvalues and eigenvectors.
39 *
40 * UPLO (input) CHARACTER*1
41 * = 'U': Upper triangles of A and B are stored;
42 * = 'L': Lower triangles of A and B are stored.
43 *
44 * N (input) INTEGER
45 * The order of the matrices A and B. N >= 0.
46 *
47 * A (input/output) COMPLEX*16 array, dimension (LDA, N)
48 * On entry, the Hermitian matrix A. If UPLO = 'U', the
49 * leading N-by-N upper triangular part of A contains the
50 * upper triangular part of the matrix A. If UPLO = 'L',
51 * the leading N-by-N lower triangular part of A contains
52 * the lower triangular part of the matrix A.
53 *
54 * On exit, if JOBZ = 'V', then if INFO = 0, A contains the
55 * matrix Z of eigenvectors. The eigenvectors are normalized
56 * as follows:
57 * if ITYPE = 1 or 2, Z**H*B*Z = I;
58 * if ITYPE = 3, Z**H*inv(B)*Z = I.
59 * If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
60 * or the lower triangle (if UPLO='L') of A, including the
61 * diagonal, is destroyed.
62 *
63 * LDA (input) INTEGER
64 * The leading dimension of the array A. LDA >= max(1,N).
65 *
66 * B (input/output) COMPLEX*16 array, dimension (LDB, N)
67 * On entry, the Hermitian positive definite matrix B.
68 * If UPLO = 'U', the leading N-by-N upper triangular part of B
69 * contains the upper triangular part of the matrix B.
70 * If UPLO = 'L', the leading N-by-N lower triangular part of B
71 * contains the lower triangular part of the matrix B.
72 *
73 * On exit, if INFO <= N, the part of B containing the matrix is
74 * overwritten by the triangular factor U or L from the Cholesky
75 * factorization B = U**H*U or B = L*L**H.
76 *
77 * LDB (input) INTEGER
78 * The leading dimension of the array B. LDB >= max(1,N).
79 *
80 * W (output) DOUBLE PRECISION array, dimension (N)
81 * If INFO = 0, the eigenvalues in ascending order.
82 *
83 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
84 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
85 *
86 * LWORK (input) INTEGER
87 * The length of the array WORK. LWORK >= max(1,2*N-1).
88 * For optimal efficiency, LWORK >= (NB+1)*N,
89 * where NB is the blocksize for ZHETRD returned by ILAENV.
90 *
91 * If LWORK = -1, then a workspace query is assumed; the routine
92 * only calculates the optimal size of the WORK array, returns
93 * this value as the first entry of the WORK array, and no error
94 * message related to LWORK is issued by XERBLA.
95 *
96 * RWORK (workspace) DOUBLE PRECISION array, dimension (max(1, 3*N-2))
97 *
98 * INFO (output) INTEGER
99 * = 0: successful exit
100 * < 0: if INFO = -i, the i-th argument had an illegal value
101 * > 0: ZPOTRF or ZHEEV returned an error code:
102 * <= N: if INFO = i, ZHEEV failed to converge;
103 * i off-diagonal elements of an intermediate
104 * tridiagonal form did not converge to zero;
105 * > N: if INFO = N + i, for 1 <= i <= N, then the leading
106 * minor of order i of B is not positive definite.
107 * The factorization of B could not be completed and
108 * no eigenvalues or eigenvectors were computed.
109 *
110 * =====================================================================
111 *
112 * .. Parameters ..
113 COMPLEX*16 ONE
114 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
115 * ..
116 * .. Local Scalars ..
117 LOGICAL LQUERY, UPPER, WANTZ
118 CHARACTER TRANS
119 INTEGER LWKOPT, NB, NEIG
120 * ..
121 * .. External Functions ..
122 LOGICAL LSAME
123 INTEGER ILAENV
124 EXTERNAL LSAME, ILAENV
125 * ..
126 * .. External Subroutines ..
127 EXTERNAL XERBLA, ZHEEV, ZHEGST, ZPOTRF, ZTRMM, ZTRSM
128 * ..
129 * .. Intrinsic Functions ..
130 INTRINSIC MAX
131 * ..
132 * .. Executable Statements ..
133 *
134 * Test the input parameters.
135 *
136 WANTZ = LSAME( JOBZ, 'V' )
137 UPPER = LSAME( UPLO, 'U' )
138 LQUERY = ( LWORK.EQ.-1 )
139 *
140 INFO = 0
141 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
142 INFO = -1
143 ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
144 INFO = -2
145 ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
146 INFO = -3
147 ELSE IF( N.LT.0 ) THEN
148 INFO = -4
149 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
150 INFO = -6
151 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
152 INFO = -8
153 END IF
154 *
155 IF( INFO.EQ.0 ) THEN
156 NB = ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 )
157 LWKOPT = MAX( 1, ( NB + 1 )*N )
158 WORK( 1 ) = LWKOPT
159 *
160 IF( LWORK.LT.MAX( 1, 2*N - 1 ) .AND. .NOT.LQUERY ) THEN
161 INFO = -11
162 END IF
163 END IF
164 *
165 IF( INFO.NE.0 ) THEN
166 CALL XERBLA( 'ZHEGV ', -INFO )
167 RETURN
168 ELSE IF( LQUERY ) THEN
169 RETURN
170 END IF
171 *
172 * Quick return if possible
173 *
174 IF( N.EQ.0 )
175 $ RETURN
176 *
177 * Form a Cholesky factorization of B.
178 *
179 CALL ZPOTRF( UPLO, N, B, LDB, INFO )
180 IF( INFO.NE.0 ) THEN
181 INFO = N + INFO
182 RETURN
183 END IF
184 *
185 * Transform problem to standard eigenvalue problem and solve.
186 *
187 CALL ZHEGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
188 CALL ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, INFO )
189 *
190 IF( WANTZ ) THEN
191 *
192 * Backtransform eigenvectors to the original problem.
193 *
194 NEIG = N
195 IF( INFO.GT.0 )
196 $ NEIG = INFO - 1
197 IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
198 *
199 * For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
200 * backtransform eigenvectors: x = inv(L)**H *y or inv(U)*y
201 *
202 IF( UPPER ) THEN
203 TRANS = 'N'
204 ELSE
205 TRANS = 'C'
206 END IF
207 *
208 CALL ZTRSM( 'Left', UPLO, TRANS, 'Non-unit', N, NEIG, ONE,
209 $ B, LDB, A, LDA )
210 *
211 ELSE IF( ITYPE.EQ.3 ) THEN
212 *
213 * For B*A*x=(lambda)*x;
214 * backtransform eigenvectors: x = L*y or U**H *y
215 *
216 IF( UPPER ) THEN
217 TRANS = 'C'
218 ELSE
219 TRANS = 'N'
220 END IF
221 *
222 CALL ZTRMM( 'Left', UPLO, TRANS, 'Non-unit', N, NEIG, ONE,
223 $ B, LDB, A, LDA )
224 END IF
225 END IF
226 *
227 WORK( 1 ) = LWKOPT
228 *
229 RETURN
230 *
231 * End of ZHEGV
232 *
233 END
2 $ LWORK, RWORK, INFO )
3 *
4 * -- LAPACK driver routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 CHARACTER JOBZ, UPLO
11 INTEGER INFO, ITYPE, LDA, LDB, LWORK, N
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION RWORK( * ), W( * )
15 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZHEGV computes all the eigenvalues, and optionally, the eigenvectors
22 * of a complex generalized Hermitian-definite eigenproblem, of the form
23 * A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x.
24 * Here A and B are assumed to be Hermitian and B is also
25 * positive definite.
26 *
27 * Arguments
28 * =========
29 *
30 * ITYPE (input) INTEGER
31 * Specifies the problem type to be solved:
32 * = 1: A*x = (lambda)*B*x
33 * = 2: A*B*x = (lambda)*x
34 * = 3: B*A*x = (lambda)*x
35 *
36 * JOBZ (input) CHARACTER*1
37 * = 'N': Compute eigenvalues only;
38 * = 'V': Compute eigenvalues and eigenvectors.
39 *
40 * UPLO (input) CHARACTER*1
41 * = 'U': Upper triangles of A and B are stored;
42 * = 'L': Lower triangles of A and B are stored.
43 *
44 * N (input) INTEGER
45 * The order of the matrices A and B. N >= 0.
46 *
47 * A (input/output) COMPLEX*16 array, dimension (LDA, N)
48 * On entry, the Hermitian matrix A. If UPLO = 'U', the
49 * leading N-by-N upper triangular part of A contains the
50 * upper triangular part of the matrix A. If UPLO = 'L',
51 * the leading N-by-N lower triangular part of A contains
52 * the lower triangular part of the matrix A.
53 *
54 * On exit, if JOBZ = 'V', then if INFO = 0, A contains the
55 * matrix Z of eigenvectors. The eigenvectors are normalized
56 * as follows:
57 * if ITYPE = 1 or 2, Z**H*B*Z = I;
58 * if ITYPE = 3, Z**H*inv(B)*Z = I.
59 * If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
60 * or the lower triangle (if UPLO='L') of A, including the
61 * diagonal, is destroyed.
62 *
63 * LDA (input) INTEGER
64 * The leading dimension of the array A. LDA >= max(1,N).
65 *
66 * B (input/output) COMPLEX*16 array, dimension (LDB, N)
67 * On entry, the Hermitian positive definite matrix B.
68 * If UPLO = 'U', the leading N-by-N upper triangular part of B
69 * contains the upper triangular part of the matrix B.
70 * If UPLO = 'L', the leading N-by-N lower triangular part of B
71 * contains the lower triangular part of the matrix B.
72 *
73 * On exit, if INFO <= N, the part of B containing the matrix is
74 * overwritten by the triangular factor U or L from the Cholesky
75 * factorization B = U**H*U or B = L*L**H.
76 *
77 * LDB (input) INTEGER
78 * The leading dimension of the array B. LDB >= max(1,N).
79 *
80 * W (output) DOUBLE PRECISION array, dimension (N)
81 * If INFO = 0, the eigenvalues in ascending order.
82 *
83 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
84 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
85 *
86 * LWORK (input) INTEGER
87 * The length of the array WORK. LWORK >= max(1,2*N-1).
88 * For optimal efficiency, LWORK >= (NB+1)*N,
89 * where NB is the blocksize for ZHETRD returned by ILAENV.
90 *
91 * If LWORK = -1, then a workspace query is assumed; the routine
92 * only calculates the optimal size of the WORK array, returns
93 * this value as the first entry of the WORK array, and no error
94 * message related to LWORK is issued by XERBLA.
95 *
96 * RWORK (workspace) DOUBLE PRECISION array, dimension (max(1, 3*N-2))
97 *
98 * INFO (output) INTEGER
99 * = 0: successful exit
100 * < 0: if INFO = -i, the i-th argument had an illegal value
101 * > 0: ZPOTRF or ZHEEV returned an error code:
102 * <= N: if INFO = i, ZHEEV failed to converge;
103 * i off-diagonal elements of an intermediate
104 * tridiagonal form did not converge to zero;
105 * > N: if INFO = N + i, for 1 <= i <= N, then the leading
106 * minor of order i of B is not positive definite.
107 * The factorization of B could not be completed and
108 * no eigenvalues or eigenvectors were computed.
109 *
110 * =====================================================================
111 *
112 * .. Parameters ..
113 COMPLEX*16 ONE
114 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
115 * ..
116 * .. Local Scalars ..
117 LOGICAL LQUERY, UPPER, WANTZ
118 CHARACTER TRANS
119 INTEGER LWKOPT, NB, NEIG
120 * ..
121 * .. External Functions ..
122 LOGICAL LSAME
123 INTEGER ILAENV
124 EXTERNAL LSAME, ILAENV
125 * ..
126 * .. External Subroutines ..
127 EXTERNAL XERBLA, ZHEEV, ZHEGST, ZPOTRF, ZTRMM, ZTRSM
128 * ..
129 * .. Intrinsic Functions ..
130 INTRINSIC MAX
131 * ..
132 * .. Executable Statements ..
133 *
134 * Test the input parameters.
135 *
136 WANTZ = LSAME( JOBZ, 'V' )
137 UPPER = LSAME( UPLO, 'U' )
138 LQUERY = ( LWORK.EQ.-1 )
139 *
140 INFO = 0
141 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
142 INFO = -1
143 ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
144 INFO = -2
145 ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
146 INFO = -3
147 ELSE IF( N.LT.0 ) THEN
148 INFO = -4
149 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
150 INFO = -6
151 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
152 INFO = -8
153 END IF
154 *
155 IF( INFO.EQ.0 ) THEN
156 NB = ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 )
157 LWKOPT = MAX( 1, ( NB + 1 )*N )
158 WORK( 1 ) = LWKOPT
159 *
160 IF( LWORK.LT.MAX( 1, 2*N - 1 ) .AND. .NOT.LQUERY ) THEN
161 INFO = -11
162 END IF
163 END IF
164 *
165 IF( INFO.NE.0 ) THEN
166 CALL XERBLA( 'ZHEGV ', -INFO )
167 RETURN
168 ELSE IF( LQUERY ) THEN
169 RETURN
170 END IF
171 *
172 * Quick return if possible
173 *
174 IF( N.EQ.0 )
175 $ RETURN
176 *
177 * Form a Cholesky factorization of B.
178 *
179 CALL ZPOTRF( UPLO, N, B, LDB, INFO )
180 IF( INFO.NE.0 ) THEN
181 INFO = N + INFO
182 RETURN
183 END IF
184 *
185 * Transform problem to standard eigenvalue problem and solve.
186 *
187 CALL ZHEGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
188 CALL ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, INFO )
189 *
190 IF( WANTZ ) THEN
191 *
192 * Backtransform eigenvectors to the original problem.
193 *
194 NEIG = N
195 IF( INFO.GT.0 )
196 $ NEIG = INFO - 1
197 IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
198 *
199 * For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
200 * backtransform eigenvectors: x = inv(L)**H *y or inv(U)*y
201 *
202 IF( UPPER ) THEN
203 TRANS = 'N'
204 ELSE
205 TRANS = 'C'
206 END IF
207 *
208 CALL ZTRSM( 'Left', UPLO, TRANS, 'Non-unit', N, NEIG, ONE,
209 $ B, LDB, A, LDA )
210 *
211 ELSE IF( ITYPE.EQ.3 ) THEN
212 *
213 * For B*A*x=(lambda)*x;
214 * backtransform eigenvectors: x = L*y or U**H *y
215 *
216 IF( UPPER ) THEN
217 TRANS = 'C'
218 ELSE
219 TRANS = 'N'
220 END IF
221 *
222 CALL ZTRMM( 'Left', UPLO, TRANS, 'Non-unit', N, NEIG, ONE,
223 $ B, LDB, A, LDA )
224 END IF
225 END IF
226 *
227 WORK( 1 ) = LWKOPT
228 *
229 RETURN
230 *
231 * End of ZHEGV
232 *
233 END