LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
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
subroutine zpotrf(UPLO, N, A, LDA, INFO)
ZPOTRF VARIANT: right looking block version of the algorithm, calling Level 3 BLAS.
Definition: zpotrf.f:101
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
Definition: ztrsm.f:181
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:61
subroutine zhegvx(ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, LWORK, RWORK, IWORK, IFAIL, INFO)
ZHEGST
Definition: zhegvx.f:297
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:54
subroutine zhegst(ITYPE, UPLO, N, A, LDA, B, LDB, INFO)
ZHEGST
Definition: zhegst.f:128
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:81
subroutine zheevx(JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, LWORK, RWORK, IWORK, IFAIL, INFO)
ZHEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices ...
Definition: zheevx.f:251
subroutine ztrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRMM
Definition: ztrmm.f:178