LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
chegvd.f
Go to the documentation of this file.
1 *> \brief \b CHEGVD
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CHEGVD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chegvd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chegvd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chegvd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CHEGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
22 * LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER JOBZ, UPLO
26 * INTEGER INFO, ITYPE, LDA, LDB, LIWORK, LRWORK, LWORK, N
27 * ..
28 * .. Array Arguments ..
29 * INTEGER IWORK( * )
30 * REAL RWORK( * ), W( * )
31 * COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> CHEGVD computes all the eigenvalues, and optionally, the eigenvectors
41 *> of a complex generalized Hermitian-definite eigenproblem, of the form
42 *> A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. Here A and
43 *> B are assumed to be Hermitian and B is also positive definite.
44 *> If eigenvectors are desired, it uses a divide and conquer algorithm.
45 *>
46 *> The divide and conquer algorithm makes very mild assumptions about
47 *> floating point arithmetic. It will work on machines with a guard
48 *> digit in add/subtract, or on those binary machines without guard
49 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
50 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
51 *> without guard digits, but we know of none.
52 *> \endverbatim
53 *
54 * Arguments:
55 * ==========
56 *
57 *> \param[in] ITYPE
58 *> \verbatim
59 *> ITYPE is INTEGER
60 *> Specifies the problem type to be solved:
61 *> = 1: A*x = (lambda)*B*x
62 *> = 2: A*B*x = (lambda)*x
63 *> = 3: B*A*x = (lambda)*x
64 *> \endverbatim
65 *>
66 *> \param[in] JOBZ
67 *> \verbatim
68 *> JOBZ is CHARACTER*1
69 *> = 'N': Compute eigenvalues only;
70 *> = 'V': Compute eigenvalues and eigenvectors.
71 *> \endverbatim
72 *>
73 *> \param[in] UPLO
74 *> \verbatim
75 *> UPLO is CHARACTER*1
76 *> = 'U': Upper triangles of A and B are stored;
77 *> = 'L': Lower triangles of A and B are stored.
78 *> \endverbatim
79 *>
80 *> \param[in] N
81 *> \verbatim
82 *> N is INTEGER
83 *> The order of the matrices A and B. N >= 0.
84 *> \endverbatim
85 *>
86 *> \param[in,out] A
87 *> \verbatim
88 *> A is COMPLEX array, dimension (LDA, N)
89 *> On entry, the Hermitian matrix A. If UPLO = 'U', the
90 *> leading N-by-N upper triangular part of A contains the
91 *> upper triangular part of the matrix A. If UPLO = 'L',
92 *> the leading N-by-N lower triangular part of A contains
93 *> the lower triangular part of the matrix A.
94 *>
95 *> On exit, if JOBZ = 'V', then if INFO = 0, A contains the
96 *> matrix Z of eigenvectors. The eigenvectors are normalized
97 *> as follows:
98 *> if ITYPE = 1 or 2, Z**H*B*Z = I;
99 *> if ITYPE = 3, Z**H*inv(B)*Z = I.
100 *> If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
101 *> or the lower triangle (if UPLO='L') of A, including the
102 *> diagonal, is destroyed.
103 *> \endverbatim
104 *>
105 *> \param[in] LDA
106 *> \verbatim
107 *> LDA is INTEGER
108 *> The leading dimension of the array A. LDA >= max(1,N).
109 *> \endverbatim
110 *>
111 *> \param[in,out] B
112 *> \verbatim
113 *> B is COMPLEX array, dimension (LDB, N)
114 *> On entry, the Hermitian matrix B. If UPLO = 'U', the
115 *> leading N-by-N upper triangular part of B contains the
116 *> upper triangular part of the matrix B. If UPLO = 'L',
117 *> the leading N-by-N lower triangular part of B contains
118 *> the lower triangular part of the matrix B.
119 *>
120 *> On exit, if INFO <= N, the part of B containing the matrix is
121 *> overwritten by the triangular factor U or L from the Cholesky
122 *> factorization B = U**H*U or B = L*L**H.
123 *> \endverbatim
124 *>
125 *> \param[in] LDB
126 *> \verbatim
127 *> LDB is INTEGER
128 *> The leading dimension of the array B. LDB >= max(1,N).
129 *> \endverbatim
130 *>
131 *> \param[out] W
132 *> \verbatim
133 *> W is REAL array, dimension (N)
134 *> If INFO = 0, the eigenvalues in ascending order.
135 *> \endverbatim
136 *>
137 *> \param[out] WORK
138 *> \verbatim
139 *> WORK is COMPLEX array, dimension (MAX(1,LWORK))
140 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
141 *> \endverbatim
142 *>
143 *> \param[in] LWORK
144 *> \verbatim
145 *> LWORK is INTEGER
146 *> The length of the array WORK.
147 *> If N <= 1, LWORK >= 1.
148 *> If JOBZ = 'N' and N > 1, LWORK >= N + 1.
149 *> If JOBZ = 'V' and N > 1, LWORK >= 2*N + N**2.
150 *>
151 *> If LWORK = -1, then a workspace query is assumed; the routine
152 *> only calculates the optimal sizes of the WORK, RWORK and
153 *> IWORK arrays, returns these values as the first entries of
154 *> the WORK, RWORK and IWORK arrays, and no error message
155 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
156 *> \endverbatim
157 *>
158 *> \param[out] RWORK
159 *> \verbatim
160 *> RWORK is REAL array, dimension (MAX(1,LRWORK))
161 *> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
162 *> \endverbatim
163 *>
164 *> \param[in] LRWORK
165 *> \verbatim
166 *> LRWORK is INTEGER
167 *> The dimension of the array RWORK.
168 *> If N <= 1, LRWORK >= 1.
169 *> If JOBZ = 'N' and N > 1, LRWORK >= N.
170 *> If JOBZ = 'V' and N > 1, LRWORK >= 1 + 5*N + 2*N**2.
171 *>
172 *> If LRWORK = -1, then a workspace query is assumed; the
173 *> routine only calculates the optimal sizes of the WORK, RWORK
174 *> and IWORK arrays, returns these values as the first entries
175 *> of the WORK, RWORK and IWORK arrays, and no error message
176 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
177 *> \endverbatim
178 *>
179 *> \param[out] IWORK
180 *> \verbatim
181 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
182 *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
183 *> \endverbatim
184 *>
185 *> \param[in] LIWORK
186 *> \verbatim
187 *> LIWORK is INTEGER
188 *> The dimension of the array IWORK.
189 *> If N <= 1, LIWORK >= 1.
190 *> If JOBZ = 'N' and N > 1, LIWORK >= 1.
191 *> If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N.
192 *>
193 *> If LIWORK = -1, then a workspace query is assumed; the
194 *> routine only calculates the optimal sizes of the WORK, RWORK
195 *> and IWORK arrays, returns these values as the first entries
196 *> of the WORK, RWORK and IWORK arrays, and no error message
197 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
198 *> \endverbatim
199 *>
200 *> \param[out] INFO
201 *> \verbatim
202 *> INFO is INTEGER
203 *> = 0: successful exit
204 *> < 0: if INFO = -i, the i-th argument had an illegal value
205 *> > 0: CPOTRF or CHEEVD returned an error code:
206 *> <= N: if INFO = i and JOBZ = 'N', then the algorithm
207 *> failed to converge; i off-diagonal elements of an
208 *> intermediate tridiagonal form did not converge to
209 *> zero;
210 *> if INFO = i and JOBZ = 'V', then the algorithm
211 *> failed to compute an eigenvalue while working on
212 *> the submatrix lying in rows and columns INFO/(N+1)
213 *> through mod(INFO,N+1);
214 *> > N: if INFO = N + i, for 1 <= i <= N, then the leading
215 *> minor of order i of B is not positive definite.
216 *> The factorization of B could not be completed and
217 *> no eigenvalues or eigenvectors were computed.
218 *> \endverbatim
219 *
220 * Authors:
221 * ========
222 *
223 *> \author Univ. of Tennessee
224 *> \author Univ. of California Berkeley
225 *> \author Univ. of Colorado Denver
226 *> \author NAG Ltd.
227 *
228 *> \ingroup complexHEeigen
229 *
230 *> \par Further Details:
231 * =====================
232 *>
233 *> \verbatim
234 *>
235 *> Modified so that no backsubstitution is performed if CHEEVD fails to
236 *> converge (NEIG in old code could be greater than N causing out of
237 *> bounds reference to A - reported by Ralf Meyer). Also corrected the
238 *> description of INFO and the test on ITYPE. Sven, 16 Feb 05.
239 *> \endverbatim
240 *
241 *> \par Contributors:
242 * ==================
243 *>
244 *> Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
245 *>
246 * =====================================================================
247  SUBROUTINE chegvd( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
248  $ LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO )
249 *
250 * -- LAPACK driver routine --
251 * -- LAPACK is a software package provided by Univ. of Tennessee, --
252 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
253 *
254 * .. Scalar Arguments ..
255  CHARACTER JOBZ, UPLO
256  INTEGER INFO, ITYPE, LDA, LDB, LIWORK, LRWORK, LWORK, N
257 * ..
258 * .. Array Arguments ..
259  INTEGER IWORK( * )
260  REAL RWORK( * ), W( * )
261  COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
262 * ..
263 *
264 * =====================================================================
265 *
266 * .. Parameters ..
267  COMPLEX CONE
268  parameter( cone = ( 1.0e+0, 0.0e+0 ) )
269 * ..
270 * .. Local Scalars ..
271  LOGICAL LQUERY, UPPER, WANTZ
272  CHARACTER TRANS
273  INTEGER LIOPT, LIWMIN, LOPT, LROPT, LRWMIN, LWMIN
274 * ..
275 * .. External Functions ..
276  LOGICAL LSAME
277  EXTERNAL lsame
278 * ..
279 * .. External Subroutines ..
280  EXTERNAL cheevd, chegst, cpotrf, ctrmm, ctrsm, xerbla
281 * ..
282 * .. Intrinsic Functions ..
283  INTRINSIC max, real
284 * ..
285 * .. Executable Statements ..
286 *
287 * Test the input parameters.
288 *
289  wantz = lsame( jobz, 'V' )
290  upper = lsame( uplo, 'U' )
291  lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 .OR. liwork.EQ.-1 )
292 *
293  info = 0
294  IF( n.LE.1 ) THEN
295  lwmin = 1
296  lrwmin = 1
297  liwmin = 1
298  ELSE IF( wantz ) THEN
299  lwmin = 2*n + n*n
300  lrwmin = 1 + 5*n + 2*n*n
301  liwmin = 3 + 5*n
302  ELSE
303  lwmin = n + 1
304  lrwmin = n
305  liwmin = 1
306  END IF
307  lopt = lwmin
308  lropt = lrwmin
309  liopt = liwmin
310  IF( itype.LT.1 .OR. itype.GT.3 ) THEN
311  info = -1
312  ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
313  info = -2
314  ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
315  info = -3
316  ELSE IF( n.LT.0 ) THEN
317  info = -4
318  ELSE IF( lda.LT.max( 1, n ) ) THEN
319  info = -6
320  ELSE IF( ldb.LT.max( 1, n ) ) THEN
321  info = -8
322  END IF
323 *
324  IF( info.EQ.0 ) THEN
325  work( 1 ) = lopt
326  rwork( 1 ) = lropt
327  iwork( 1 ) = liopt
328 *
329  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
330  info = -11
331  ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
332  info = -13
333  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
334  info = -15
335  END IF
336  END IF
337 *
338  IF( info.NE.0 ) THEN
339  CALL xerbla( 'CHEGVD', -info )
340  RETURN
341  ELSE IF( lquery ) THEN
342  RETURN
343  END IF
344 *
345 * Quick return if possible
346 *
347  IF( n.EQ.0 )
348  $ RETURN
349 *
350 * Form a Cholesky factorization of B.
351 *
352  CALL cpotrf( uplo, n, b, ldb, info )
353  IF( info.NE.0 ) THEN
354  info = n + info
355  RETURN
356  END IF
357 *
358 * Transform problem to standard eigenvalue problem and solve.
359 *
360  CALL chegst( itype, uplo, n, a, lda, b, ldb, info )
361  CALL cheevd( jobz, uplo, n, a, lda, w, work, lwork, rwork, lrwork,
362  $ iwork, liwork, info )
363  lopt = max( real( lopt ), real( work( 1 ) ) )
364  lropt = max( real( lropt ), real( rwork( 1 ) ) )
365  liopt = max( real( liopt ), real( iwork( 1 ) ) )
366 *
367  IF( wantz .AND. info.EQ.0 ) THEN
368 *
369 * Backtransform eigenvectors to the original problem.
370 *
371  IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
372 *
373 * For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
374 * backtransform eigenvectors: x = inv(L)**H *y or inv(U)*y
375 *
376  IF( upper ) THEN
377  trans = 'N'
378  ELSE
379  trans = 'C'
380  END IF
381 *
382  CALL ctrsm( 'Left', uplo, trans, 'Non-unit', n, n, cone,
383  $ b, ldb, a, lda )
384 *
385  ELSE IF( itype.EQ.3 ) THEN
386 *
387 * For B*A*x=(lambda)*x;
388 * backtransform eigenvectors: x = L*y or U**H *y
389 *
390  IF( upper ) THEN
391  trans = 'C'
392  ELSE
393  trans = 'N'
394  END IF
395 *
396  CALL ctrmm( 'Left', uplo, trans, 'Non-unit', n, n, cone,
397  $ b, ldb, a, lda )
398  END IF
399  END IF
400 *
401  work( 1 ) = lopt
402  rwork( 1 ) = lropt
403  iwork( 1 ) = liopt
404 *
405  RETURN
406 *
407 * End of CHEGVD
408 *
409  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ctrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRMM
Definition: ctrmm.f:177
subroutine ctrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRSM
Definition: ctrsm.f:180
subroutine chegst(ITYPE, UPLO, N, A, LDA, B, LDB, INFO)
CHEGST
Definition: chegst.f:128
subroutine chegvd(ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
CHEGVD
Definition: chegvd.f:249
subroutine cheevd(JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
CHEEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices
Definition: cheevd.f:205
subroutine cpotrf(UPLO, N, A, LDA, INFO)
CPOTRF
Definition: cpotrf.f:107