LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zhegvx.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 ZHEGVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhegvx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhegvx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhegvx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZHEGVX( ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB,
22 * VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK,
23 * LWORK, RWORK, IWORK, IFAIL, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER JOBZ, RANGE, UPLO
27 * INTEGER IL, INFO, ITYPE, IU, LDA, LDB, LDZ, LWORK, M, N
28 * DOUBLE PRECISION ABSTOL, VL, VU
29 * ..
30 * .. Array Arguments ..
31 * INTEGER IFAIL( * ), IWORK( * )
32 * DOUBLE PRECISION RWORK( * ), W( * )
33 * COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ),
34 * $ Z( LDZ, * )
35 * ..
36 *
37 *
38 *> \par Purpose:
39 * =============
40 *>
41 *> \verbatim
42 *>
43 *> ZHEGVX computes selected eigenvalues, and optionally, eigenvectors
44 *> of a complex generalized Hermitian-definite eigenproblem, of the form
45 *> A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. Here A and
46 *> B are assumed to be Hermitian and B is also positive definite.
47 *> Eigenvalues and eigenvectors can be selected by specifying either a
48 *> range of values or a range of indices for the desired eigenvalues.
49 *> \endverbatim
50 *
51 * Arguments:
52 * ==========
53 *
54 *> \param[in] ITYPE
55 *> \verbatim
56 *> ITYPE is INTEGER
57 *> Specifies the problem type to be solved:
58 *> = 1: A*x = (lambda)*B*x
59 *> = 2: A*B*x = (lambda)*x
60 *> = 3: B*A*x = (lambda)*x
61 *> \endverbatim
62 *>
63 *> \param[in] JOBZ
64 *> \verbatim
65 *> JOBZ is CHARACTER*1
66 *> = 'N': Compute eigenvalues only;
67 *> = 'V': Compute eigenvalues and eigenvectors.
68 *> \endverbatim
69 *>
70 *> \param[in] RANGE
71 *> \verbatim
72 *> RANGE is CHARACTER*1
73 *> = 'A': all eigenvalues will be found.
74 *> = 'V': all eigenvalues in the half-open interval (VL,VU]
75 *> will be found.
76 *> = 'I': the IL-th through IU-th eigenvalues will be found.
77 *> \endverbatim
78 *>
79 *> \param[in] UPLO
80 *> \verbatim
81 *> UPLO is CHARACTER*1
82 *> = 'U': Upper triangles of A and B are stored;
83 *> = 'L': Lower triangles of A and B are stored.
84 *> \endverbatim
85 *>
86 *> \param[in] N
87 *> \verbatim
88 *> N is INTEGER
89 *> The order of the matrices A and B. N >= 0.
90 *> \endverbatim
91 *>
92 *> \param[in,out] A
93 *> \verbatim
94 *> A is COMPLEX*16 array, dimension (LDA, N)
95 *> On entry, the Hermitian matrix A. If UPLO = 'U', the
96 *> leading N-by-N upper triangular part of A contains the
97 *> upper triangular part of the matrix A. If UPLO = 'L',
98 *> the leading N-by-N lower triangular part of A contains
99 *> the lower triangular part of the matrix A.
100 *>
101 *> On exit, the lower triangle (if UPLO='L') or the upper
102 *> triangle (if UPLO='U') of A, including the diagonal, is
103 *> destroyed.
104 *> \endverbatim
105 *>
106 *> \param[in] LDA
107 *> \verbatim
108 *> LDA is INTEGER
109 *> The leading dimension of the array A. LDA >= max(1,N).
110 *> \endverbatim
111 *>
112 *> \param[in,out] B
113 *> \verbatim
114 *> B is COMPLEX*16 array, dimension (LDB, N)
115 *> On entry, the Hermitian matrix B. If UPLO = 'U', the
116 *> leading N-by-N upper triangular part of B contains the
117 *> upper triangular part of the matrix B. If UPLO = 'L',
118 *> the leading N-by-N lower triangular part of B contains
119 *> the lower triangular part of the matrix B.
120 *>
121 *> On exit, if INFO <= N, the part of B containing the matrix is
122 *> overwritten by the triangular factor U or L from the Cholesky
123 *> factorization B = U**H*U or B = L*L**H.
124 *> \endverbatim
125 *>
126 *> \param[in] LDB
127 *> \verbatim
128 *> LDB is INTEGER
129 *> The leading dimension of the array B. LDB >= max(1,N).
130 *> \endverbatim
131 *>
132 *> \param[in] VL
133 *> \verbatim
134 *> VL is DOUBLE PRECISION
135 *> \endverbatim
136 *>
137 *> \param[in] VU
138 *> \verbatim
139 *> VU is DOUBLE PRECISION
140 *>
141 *> If RANGE='V', the lower and upper bounds of the interval to
142 *> be searched for eigenvalues. VL < VU.
143 *> Not referenced if RANGE = 'A' or 'I'.
144 *> \endverbatim
145 *>
146 *> \param[in] IL
147 *> \verbatim
148 *> IL is INTEGER
149 *> \endverbatim
150 *>
151 *> \param[in] IU
152 *> \verbatim
153 *> IU is INTEGER
154 *>
155 *> If RANGE='I', the indices (in ascending order) of the
156 *> smallest and largest eigenvalues to be returned.
157 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
158 *> Not referenced if RANGE = 'A' or 'V'.
159 *> \endverbatim
160 *>
161 *> \param[in] ABSTOL
162 *> \verbatim
163 *> ABSTOL is DOUBLE PRECISION
164 *> The absolute error tolerance for the eigenvalues.
165 *> An approximate eigenvalue is accepted as converged
166 *> when it is determined to lie in an interval [a,b]
167 *> of width less than or equal to
168 *>
169 *> ABSTOL + EPS * max( |a|,|b| ) ,
170 *>
171 *> where EPS is the machine precision. If ABSTOL is less than
172 *> or equal to zero, then EPS*|T| will be used in its place,
173 *> where |T| is the 1-norm of the tridiagonal matrix obtained
174 *> by reducing C to tridiagonal form, where C is the symmetric
175 *> matrix of the standard symmetric problem to which the
176 *> generalized problem is transformed.
177 *>
178 *> Eigenvalues will be computed most accurately when ABSTOL is
179 *> set to twice the underflow threshold 2*DLAMCH('S'), not zero.
180 *> If this routine returns with INFO>0, indicating that some
181 *> eigenvectors did not converge, try setting ABSTOL to
182 *> 2*DLAMCH('S').
183 *> \endverbatim
184 *>
185 *> \param[out] M
186 *> \verbatim
187 *> M is INTEGER
188 *> The total number of eigenvalues found. 0 <= M <= N.
189 *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
190 *> \endverbatim
191 *>
192 *> \param[out] W
193 *> \verbatim
194 *> W is DOUBLE PRECISION array, dimension (N)
195 *> The first M elements contain the selected
196 *> eigenvalues in ascending order.
197 *> \endverbatim
198 *>
199 *> \param[out] Z
200 *> \verbatim
201 *> Z is COMPLEX*16 array, dimension (LDZ, max(1,M))
202 *> If JOBZ = 'N', then Z is not referenced.
203 *> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
204 *> contain the orthonormal eigenvectors of the matrix A
205 *> corresponding to the selected eigenvalues, with the i-th
206 *> column of Z holding the eigenvector associated with W(i).
207 *> The eigenvectors are normalized as follows:
208 *> if ITYPE = 1 or 2, Z**T*B*Z = I;
209 *> if ITYPE = 3, Z**T*inv(B)*Z = I.
210 *>
211 *> If an eigenvector fails to converge, then that column of Z
212 *> contains the latest approximation to the eigenvector, and the
213 *> index of the eigenvector is returned in IFAIL.
214 *> Note: the user must ensure that at least max(1,M) columns are
215 *> supplied in the array Z; if RANGE = 'V', the exact value of M
216 *> is not known in advance and an upper bound must be used.
217 *> \endverbatim
218 *>
219 *> \param[in] LDZ
220 *> \verbatim
221 *> LDZ is INTEGER
222 *> The leading dimension of the array Z. LDZ >= 1, and if
223 *> JOBZ = 'V', LDZ >= max(1,N).
224 *> \endverbatim
225 *>
226 *> \param[out] WORK
227 *> \verbatim
228 *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
229 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
230 *> \endverbatim
231 *>
232 *> \param[in] LWORK
233 *> \verbatim
234 *> LWORK is INTEGER
235 *> The length of the array WORK. LWORK >= max(1,2*N).
236 *> For optimal efficiency, LWORK >= (NB+1)*N,
237 *> where NB is the blocksize for ZHETRD returned by ILAENV.
238 *>
239 *> If LWORK = -1, then a workspace query is assumed; the routine
240 *> only calculates the optimal size of the WORK array, returns
241 *> this value as the first entry of the WORK array, and no error
242 *> message related to LWORK is issued by XERBLA.
243 *> \endverbatim
244 *>
245 *> \param[out] RWORK
246 *> \verbatim
247 *> RWORK is DOUBLE PRECISION array, dimension (7*N)
248 *> \endverbatim
249 *>
250 *> \param[out] IWORK
251 *> \verbatim
252 *> IWORK is INTEGER array, dimension (5*N)
253 *> \endverbatim
254 *>
255 *> \param[out] IFAIL
256 *> \verbatim
257 *> IFAIL is INTEGER array, dimension (N)
258 *> If JOBZ = 'V', then if INFO = 0, the first M elements of
259 *> IFAIL are zero. If INFO > 0, then IFAIL contains the
260 *> indices of the eigenvectors that failed to converge.
261 *> If JOBZ = 'N', then IFAIL is not referenced.
262 *> \endverbatim
263 *>
264 *> \param[out] INFO
265 *> \verbatim
266 *> INFO is INTEGER
267 *> = 0: successful exit
268 *> < 0: if INFO = -i, the i-th argument had an illegal value
269 *> > 0: ZPOTRF or ZHEEVX returned an error code:
270 *> <= N: if INFO = i, ZHEEVX failed to converge;
271 *> i eigenvectors failed to converge. Their indices
272 *> are stored in array IFAIL.
273 *> > N: if INFO = N + i, for 1 <= i <= N, then the leading
274 *> minor of order i of B is not positive definite.
275 *> The factorization of B could not be completed and
276 *> no eigenvalues or eigenvectors were computed.
277 *> \endverbatim
278 *
279 * Authors:
280 * ========
281 *
282 *> \author Univ. of Tennessee
283 *> \author Univ. of California Berkeley
284 *> \author Univ. of Colorado Denver
285 *> \author NAG Ltd.
286 *
287 *> \date November 2011
288 *
289 *> \ingroup complex16HEeigen
290 *
291 *> \par Contributors:
292 * ==================
293 *>
294 *> Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
295 *
296 * =====================================================================
297  SUBROUTINE zhegvx( ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB,
298  $ vl, vu, il, iu, abstol, m, w, z, ldz, work,
299  $ lwork, rwork, iwork, ifail, info )
300 *
301 * -- LAPACK driver routine (version 3.4.0) --
302 * -- LAPACK is a software package provided by Univ. of Tennessee, --
303 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
304 * November 2011
305 *
306 * .. Scalar Arguments ..
307  CHARACTER jobz, range, uplo
308  INTEGER il, info, itype, iu, lda, ldb, ldz, lwork, m, n
309  DOUBLE PRECISION abstol, vl, vu
310 * ..
311 * .. Array Arguments ..
312  INTEGER ifail( * ), iwork( * )
313  DOUBLE PRECISION rwork( * ), w( * )
314  COMPLEX*16 a( lda, * ), b( ldb, * ), work( * ),
315  $ z( ldz, * )
316 * ..
317 *
318 * =====================================================================
319 *
320 * .. Parameters ..
321  COMPLEX*16 cone
322  parameter( cone = ( 1.0d+0, 0.0d+0 ) )
323 * ..
324 * .. Local Scalars ..
325  LOGICAL alleig, indeig, lquery, upper, valeig, wantz
326  CHARACTER trans
327  INTEGER lwkopt, nb
328 * ..
329 * .. External Functions ..
330  LOGICAL lsame
331  INTEGER ilaenv
332  EXTERNAL lsame, ilaenv
333 * ..
334 * .. External Subroutines ..
335  EXTERNAL xerbla, zheevx, zhegst, zpotrf, ztrmm, ztrsm
336 * ..
337 * .. Intrinsic Functions ..
338  INTRINSIC max, min
339 * ..
340 * .. Executable Statements ..
341 *
342 * Test the input parameters.
343 *
344  wantz = lsame( jobz, 'V' )
345  upper = lsame( uplo, 'U' )
346  alleig = lsame( range, 'A' )
347  valeig = lsame( range, 'V' )
348  indeig = lsame( range, 'I' )
349  lquery = ( lwork.EQ.-1 )
350 *
351  info = 0
352  IF( itype.LT.1 .OR. itype.GT.3 ) THEN
353  info = -1
354  ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
355  info = -2
356  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
357  info = -3
358  ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
359  info = -4
360  ELSE IF( n.LT.0 ) THEN
361  info = -5
362  ELSE IF( lda.LT.max( 1, n ) ) THEN
363  info = -7
364  ELSE IF( ldb.LT.max( 1, n ) ) THEN
365  info = -9
366  ELSE
367  IF( valeig ) THEN
368  IF( n.GT.0 .AND. vu.LE.vl )
369  $ info = -11
370  ELSE IF( indeig ) THEN
371  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
372  info = -12
373  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
374  info = -13
375  END IF
376  END IF
377  END IF
378  IF (info.EQ.0) THEN
379  IF (ldz.LT.1 .OR. (wantz .AND. ldz.LT.n)) THEN
380  info = -18
381  END IF
382  END IF
383 *
384  IF( info.EQ.0 ) THEN
385  nb = ilaenv( 1, 'ZHETRD', uplo, n, -1, -1, -1 )
386  lwkopt = max( 1, ( nb + 1 )*n )
387  work( 1 ) = lwkopt
388 *
389  IF( lwork.LT.max( 1, 2*n ) .AND. .NOT.lquery ) THEN
390  info = -20
391  END IF
392  END IF
393 *
394  IF( info.NE.0 ) THEN
395  CALL xerbla( 'ZHEGVX', -info )
396  return
397  ELSE IF( lquery ) THEN
398  return
399  END IF
400 *
401 * Quick return if possible
402 *
403  m = 0
404  IF( n.EQ.0 ) THEN
405  return
406  END IF
407 *
408 * Form a Cholesky factorization of B.
409 *
410  CALL zpotrf( uplo, n, b, ldb, info )
411  IF( info.NE.0 ) THEN
412  info = n + info
413  return
414  END IF
415 *
416 * Transform problem to standard eigenvalue problem and solve.
417 *
418  CALL zhegst( itype, uplo, n, a, lda, b, ldb, info )
419  CALL zheevx( jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol,
420  $ m, w, z, ldz, work, lwork, rwork, iwork, ifail,
421  $ info )
422 *
423  IF( wantz ) THEN
424 *
425 * Backtransform eigenvectors to the original problem.
426 *
427  IF( info.GT.0 )
428  $ m = info - 1
429  IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
430 *
431 * For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
432 * backtransform eigenvectors: x = inv(L)**H *y or inv(U)*y
433 *
434  IF( upper ) THEN
435  trans = 'N'
436  ELSE
437  trans = 'C'
438  END IF
439 *
440  CALL ztrsm( 'Left', uplo, trans, 'Non-unit', n, m, cone, b,
441  $ ldb, z, ldz )
442 *
443  ELSE IF( itype.EQ.3 ) THEN
444 *
445 * For B*A*x=(lambda)*x;
446 * backtransform eigenvectors: x = L*y or U**H *y
447 *
448  IF( upper ) THEN
449  trans = 'C'
450  ELSE
451  trans = 'N'
452  END IF
453 *
454  CALL ztrmm( 'Left', uplo, trans, 'Non-unit', n, m, cone, b,
455  $ ldb, z, ldz )
456  END IF
457  END IF
458 *
459 * Set WORK(1) to optimal complex workspace size.
460 *
461  work( 1 ) = lwkopt
462 *
463  return
464 *
465 * End of ZHEGVX
466 *
467  END