LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zhegvd.f
Go to the documentation of this file.
1 *> \brief \b ZHEGST
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZHEGVD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhegvd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhegvd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhegvd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZHEGVD( 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 * DOUBLE PRECISION RWORK( * ), W( * )
31 * COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> ZHEGVD 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*16 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*16 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 DOUBLE PRECISION 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*16 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 DOUBLE PRECISION 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: ZPOTRF or ZHEEVD 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 *> \date November 2011
229 *
230 *> \ingroup complex16HEeigen
231 *
232 *> \par Further Details:
233 * =====================
234 *>
235 *> \verbatim
236 *>
237 *> Modified so that no backsubstitution is performed if ZHEEVD fails to
238 *> converge (NEIG in old code could be greater than N causing out of
239 *> bounds reference to A - reported by Ralf Meyer). Also corrected the
240 *> description of INFO and the test on ITYPE. Sven, 16 Feb 05.
241 *> \endverbatim
242 *
243 *> \par Contributors:
244 * ==================
245 *>
246 *> Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
247 *>
248 * =====================================================================
249  SUBROUTINE zhegvd( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
250  $ lwork, rwork, lrwork, iwork, liwork, info )
251 *
252 * -- LAPACK driver routine (version 3.4.0) --
253 * -- LAPACK is a software package provided by Univ. of Tennessee, --
254 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
255 * November 2011
256 *
257 * .. Scalar Arguments ..
258  CHARACTER jobz, uplo
259  INTEGER info, itype, lda, ldb, liwork, lrwork, lwork, n
260 * ..
261 * .. Array Arguments ..
262  INTEGER iwork( * )
263  DOUBLE PRECISION rwork( * ), w( * )
264  COMPLEX*16 a( lda, * ), b( ldb, * ), work( * )
265 * ..
266 *
267 * =====================================================================
268 *
269 * .. Parameters ..
270  COMPLEX*16 cone
271  parameter( cone = ( 1.0d+0, 0.0d+0 ) )
272 * ..
273 * .. Local Scalars ..
274  LOGICAL lquery, upper, wantz
275  CHARACTER trans
276  INTEGER liopt, liwmin, lopt, lropt, lrwmin, lwmin
277 * ..
278 * .. External Functions ..
279  LOGICAL lsame
280  EXTERNAL lsame
281 * ..
282 * .. External Subroutines ..
283  EXTERNAL xerbla, zheevd, zhegst, zpotrf, ztrmm, ztrsm
284 * ..
285 * .. Intrinsic Functions ..
286  INTRINSIC dble, max
287 * ..
288 * .. Executable Statements ..
289 *
290 * Test the input parameters.
291 *
292  wantz = lsame( jobz, 'V' )
293  upper = lsame( uplo, 'U' )
294  lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 .OR. liwork.EQ.-1 )
295 *
296  info = 0
297  IF( n.LE.1 ) THEN
298  lwmin = 1
299  lrwmin = 1
300  liwmin = 1
301  ELSE IF( wantz ) THEN
302  lwmin = 2*n + n*n
303  lrwmin = 1 + 5*n + 2*n*n
304  liwmin = 3 + 5*n
305  ELSE
306  lwmin = n + 1
307  lrwmin = n
308  liwmin = 1
309  END IF
310  lopt = lwmin
311  lropt = lrwmin
312  liopt = liwmin
313  IF( itype.LT.1 .OR. itype.GT.3 ) THEN
314  info = -1
315  ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
316  info = -2
317  ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
318  info = -3
319  ELSE IF( n.LT.0 ) THEN
320  info = -4
321  ELSE IF( lda.LT.max( 1, n ) ) THEN
322  info = -6
323  ELSE IF( ldb.LT.max( 1, n ) ) THEN
324  info = -8
325  END IF
326 *
327  IF( info.EQ.0 ) THEN
328  work( 1 ) = lopt
329  rwork( 1 ) = lropt
330  iwork( 1 ) = liopt
331 *
332  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
333  info = -11
334  ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
335  info = -13
336  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
337  info = -15
338  END IF
339  END IF
340 *
341  IF( info.NE.0 ) THEN
342  CALL xerbla( 'ZHEGVD', -info )
343  return
344  ELSE IF( lquery ) THEN
345  return
346  END IF
347 *
348 * Quick return if possible
349 *
350  IF( n.EQ.0 )
351  $ return
352 *
353 * Form a Cholesky factorization of B.
354 *
355  CALL zpotrf( uplo, n, b, ldb, info )
356  IF( info.NE.0 ) THEN
357  info = n + info
358  return
359  END IF
360 *
361 * Transform problem to standard eigenvalue problem and solve.
362 *
363  CALL zhegst( itype, uplo, n, a, lda, b, ldb, info )
364  CALL zheevd( jobz, uplo, n, a, lda, w, work, lwork, rwork, lrwork,
365  $ iwork, liwork, info )
366  lopt = max( dble( lopt ), dble( work( 1 ) ) )
367  lropt = max( dble( lropt ), dble( rwork( 1 ) ) )
368  liopt = max( dble( liopt ), dble( iwork( 1 ) ) )
369 *
370  IF( wantz .AND. info.EQ.0 ) THEN
371 *
372 * Backtransform eigenvectors to the original problem.
373 *
374  IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
375 *
376 * For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
377 * backtransform eigenvectors: x = inv(L)**H *y or inv(U)*y
378 *
379  IF( upper ) THEN
380  trans = 'N'
381  ELSE
382  trans = 'C'
383  END IF
384 *
385  CALL ztrsm( 'Left', uplo, trans, 'Non-unit', n, n, cone,
386  $ b, ldb, a, lda )
387 *
388  ELSE IF( itype.EQ.3 ) THEN
389 *
390 * For B*A*x=(lambda)*x;
391 * backtransform eigenvectors: x = L*y or U**H *y
392 *
393  IF( upper ) THEN
394  trans = 'C'
395  ELSE
396  trans = 'N'
397  END IF
398 *
399  CALL ztrmm( 'Left', uplo, trans, 'Non-unit', n, n, cone,
400  $ b, ldb, a, lda )
401  END IF
402  END IF
403 *
404  work( 1 ) = lopt
405  rwork( 1 ) = lropt
406  iwork( 1 ) = liopt
407 *
408  return
409 *
410 * End of ZHEGVD
411 *
412  END