LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
zheevx.f
Go to the documentation of this file.
1 *> \brief <b> ZHEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices</b>
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZHEEVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zheevx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zheevx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zheevx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZHEEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU,
22 * ABSTOL, M, W, Z, LDZ, WORK, LWORK, RWORK,
23 * IWORK, IFAIL, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER JOBZ, RANGE, UPLO
27 * INTEGER IL, INFO, IU, LDA, 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, * ), WORK( * ), Z( LDZ, * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> ZHEEVX computes selected eigenvalues and, optionally, eigenvectors
43 *> of a complex Hermitian matrix A. Eigenvalues and eigenvectors can
44 *> be selected by specifying either a range of values or a range of
45 *> indices for the desired eigenvalues.
46 *> \endverbatim
47 *
48 * Arguments:
49 * ==========
50 *
51 *> \param[in] JOBZ
52 *> \verbatim
53 *> JOBZ is CHARACTER*1
54 *> = 'N': Compute eigenvalues only;
55 *> = 'V': Compute eigenvalues and eigenvectors.
56 *> \endverbatim
57 *>
58 *> \param[in] RANGE
59 *> \verbatim
60 *> RANGE is CHARACTER*1
61 *> = 'A': all eigenvalues will be found.
62 *> = 'V': all eigenvalues in the half-open interval (VL,VU]
63 *> will be found.
64 *> = 'I': the IL-th through IU-th eigenvalues will be found.
65 *> \endverbatim
66 *>
67 *> \param[in] UPLO
68 *> \verbatim
69 *> UPLO is CHARACTER*1
70 *> = 'U': Upper triangle of A is stored;
71 *> = 'L': Lower triangle of A is stored.
72 *> \endverbatim
73 *>
74 *> \param[in] N
75 *> \verbatim
76 *> N is INTEGER
77 *> The order of the matrix A. N >= 0.
78 *> \endverbatim
79 *>
80 *> \param[in,out] A
81 *> \verbatim
82 *> A is COMPLEX*16 array, dimension (LDA, N)
83 *> On entry, the Hermitian matrix A. If UPLO = 'U', the
84 *> leading N-by-N upper triangular part of A contains the
85 *> upper triangular part of the matrix A. If UPLO = 'L',
86 *> the leading N-by-N lower triangular part of A contains
87 *> the lower triangular part of the matrix A.
88 *> On exit, the lower triangle (if UPLO='L') or the upper
89 *> triangle (if UPLO='U') of A, including the diagonal, is
90 *> destroyed.
91 *> \endverbatim
92 *>
93 *> \param[in] LDA
94 *> \verbatim
95 *> LDA is INTEGER
96 *> The leading dimension of the array A. LDA >= max(1,N).
97 *> \endverbatim
98 *>
99 *> \param[in] VL
100 *> \verbatim
101 *> VL is DOUBLE PRECISION
102 *> \endverbatim
103 *>
104 *> \param[in] VU
105 *> \verbatim
106 *> VU is DOUBLE PRECISION
107 *> If RANGE='V', the lower and upper bounds of the interval to
108 *> be searched for eigenvalues. VL < VU.
109 *> Not referenced if RANGE = 'A' or 'I'.
110 *> \endverbatim
111 *>
112 *> \param[in] IL
113 *> \verbatim
114 *> IL is INTEGER
115 *> \endverbatim
116 *>
117 *> \param[in] IU
118 *> \verbatim
119 *> IU is INTEGER
120 *> If RANGE='I', the indices (in ascending order) of the
121 *> smallest and largest eigenvalues to be returned.
122 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
123 *> Not referenced if RANGE = 'A' or 'V'.
124 *> \endverbatim
125 *>
126 *> \param[in] ABSTOL
127 *> \verbatim
128 *> ABSTOL is DOUBLE PRECISION
129 *> The absolute error tolerance for the eigenvalues.
130 *> An approximate eigenvalue is accepted as converged
131 *> when it is determined to lie in an interval [a,b]
132 *> of width less than or equal to
133 *>
134 *> ABSTOL + EPS * max( |a|,|b| ) ,
135 *>
136 *> where EPS is the machine precision. If ABSTOL is less than
137 *> or equal to zero, then EPS*|T| will be used in its place,
138 *> where |T| is the 1-norm of the tridiagonal matrix obtained
139 *> by reducing A to tridiagonal form.
140 *>
141 *> Eigenvalues will be computed most accurately when ABSTOL is
142 *> set to twice the underflow threshold 2*DLAMCH('S'), not zero.
143 *> If this routine returns with INFO>0, indicating that some
144 *> eigenvectors did not converge, try setting ABSTOL to
145 *> 2*DLAMCH('S').
146 *>
147 *> See "Computing Small Singular Values of Bidiagonal Matrices
148 *> with Guaranteed High Relative Accuracy," by Demmel and
149 *> Kahan, LAPACK Working Note #3.
150 *> \endverbatim
151 *>
152 *> \param[out] M
153 *> \verbatim
154 *> M is INTEGER
155 *> The total number of eigenvalues found. 0 <= M <= N.
156 *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
157 *> \endverbatim
158 *>
159 *> \param[out] W
160 *> \verbatim
161 *> W is DOUBLE PRECISION array, dimension (N)
162 *> On normal exit, the first M elements contain the selected
163 *> eigenvalues in ascending order.
164 *> \endverbatim
165 *>
166 *> \param[out] Z
167 *> \verbatim
168 *> Z is COMPLEX*16 array, dimension (LDZ, max(1,M))
169 *> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
170 *> contain the orthonormal eigenvectors of the matrix A
171 *> corresponding to the selected eigenvalues, with the i-th
172 *> column of Z holding the eigenvector associated with W(i).
173 *> If an eigenvector fails to converge, then that column of Z
174 *> contains the latest approximation to the eigenvector, and the
175 *> index of the eigenvector is returned in IFAIL.
176 *> If JOBZ = 'N', then Z is not referenced.
177 *> Note: the user must ensure that at least max(1,M) columns are
178 *> supplied in the array Z; if RANGE = 'V', the exact value of M
179 *> is not known in advance and an upper bound must be used.
180 *> \endverbatim
181 *>
182 *> \param[in] LDZ
183 *> \verbatim
184 *> LDZ is INTEGER
185 *> The leading dimension of the array Z. LDZ >= 1, and if
186 *> JOBZ = 'V', LDZ >= max(1,N).
187 *> \endverbatim
188 *>
189 *> \param[out] WORK
190 *> \verbatim
191 *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
192 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
193 *> \endverbatim
194 *>
195 *> \param[in] LWORK
196 *> \verbatim
197 *> LWORK is INTEGER
198 *> The length of the array WORK. LWORK >= 1, when N <= 1;
199 *> otherwise 2*N.
200 *> For optimal efficiency, LWORK >= (NB+1)*N,
201 *> where NB is the max of the blocksize for ZHETRD and for
202 *> ZUNMTR as returned by ILAENV.
203 *>
204 *> If LWORK = -1, then a workspace query is assumed; the routine
205 *> only calculates the optimal size of the WORK array, returns
206 *> this value as the first entry of the WORK array, and no error
207 *> message related to LWORK is issued by XERBLA.
208 *> \endverbatim
209 *>
210 *> \param[out] RWORK
211 *> \verbatim
212 *> RWORK is DOUBLE PRECISION array, dimension (7*N)
213 *> \endverbatim
214 *>
215 *> \param[out] IWORK
216 *> \verbatim
217 *> IWORK is INTEGER array, dimension (5*N)
218 *> \endverbatim
219 *>
220 *> \param[out] IFAIL
221 *> \verbatim
222 *> IFAIL is INTEGER array, dimension (N)
223 *> If JOBZ = 'V', then if INFO = 0, the first M elements of
224 *> IFAIL are zero. If INFO > 0, then IFAIL contains the
225 *> indices of the eigenvectors that failed to converge.
226 *> If JOBZ = 'N', then IFAIL is not referenced.
227 *> \endverbatim
228 *>
229 *> \param[out] INFO
230 *> \verbatim
231 *> INFO is INTEGER
232 *> = 0: successful exit
233 *> < 0: if INFO = -i, the i-th argument had an illegal value
234 *> > 0: if INFO = i, then i eigenvectors failed to converge.
235 *> Their indices are stored in array IFAIL.
236 *> \endverbatim
237 *
238 * Authors:
239 * ========
240 *
241 *> \author Univ. of Tennessee
242 *> \author Univ. of California Berkeley
243 *> \author Univ. of Colorado Denver
244 *> \author NAG Ltd.
245 *
246 *> \date November 2011
247 *
248 *> \ingroup complex16HEeigen
249 *
250 * =====================================================================
251  SUBROUTINE zheevx( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU,
252  $ abstol, m, w, z, ldz, work, lwork, rwork,
253  $ iwork, ifail, info )
254 *
255 * -- LAPACK driver routine (version 3.4.0) --
256 * -- LAPACK is a software package provided by Univ. of Tennessee, --
257 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
258 * November 2011
259 *
260 * .. Scalar Arguments ..
261  CHARACTER JOBZ, RANGE, UPLO
262  INTEGER IL, INFO, IU, LDA, LDZ, LWORK, M, N
263  DOUBLE PRECISION ABSTOL, VL, VU
264 * ..
265 * .. Array Arguments ..
266  INTEGER IFAIL( * ), IWORK( * )
267  DOUBLE PRECISION RWORK( * ), W( * )
268  COMPLEX*16 A( lda, * ), WORK( * ), Z( ldz, * )
269 * ..
270 *
271 * =====================================================================
272 *
273 * .. Parameters ..
274  DOUBLE PRECISION ZERO, ONE
275  parameter( zero = 0.0d+0, one = 1.0d+0 )
276  COMPLEX*16 CONE
277  parameter( cone = ( 1.0d+0, 0.0d+0 ) )
278 * ..
279 * .. Local Scalars ..
280  LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
281  $ wantz
282  CHARACTER ORDER
283  INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
284  $ indisp, indiwk, indrwk, indtau, indwrk, iscale,
285  $ itmp1, j, jj, llwork, lwkmin, lwkopt, nb,
286  $ nsplit
287  DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
288  $ sigma, smlnum, tmp1, vll, vuu
289 * ..
290 * .. External Functions ..
291  LOGICAL LSAME
292  INTEGER ILAENV
293  DOUBLE PRECISION DLAMCH, ZLANHE
294  EXTERNAL lsame, ilaenv, dlamch, zlanhe
295 * ..
296 * .. External Subroutines ..
297  EXTERNAL dcopy, dscal, dstebz, dsterf, xerbla, zdscal,
299  $ zunmtr
300 * ..
301 * .. Intrinsic Functions ..
302  INTRINSIC dble, max, min, sqrt
303 * ..
304 * .. Executable Statements ..
305 *
306 * Test the input parameters.
307 *
308  lower = lsame( uplo, 'L' )
309  wantz = lsame( jobz, 'V' )
310  alleig = lsame( range, 'A' )
311  valeig = lsame( range, 'V' )
312  indeig = lsame( range, 'I' )
313  lquery = ( lwork.EQ.-1 )
314 *
315  info = 0
316  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
317  info = -1
318  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
319  info = -2
320  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
321  info = -3
322  ELSE IF( n.LT.0 ) THEN
323  info = -4
324  ELSE IF( lda.LT.max( 1, n ) ) THEN
325  info = -6
326  ELSE
327  IF( valeig ) THEN
328  IF( n.GT.0 .AND. vu.LE.vl )
329  $ info = -8
330  ELSE IF( indeig ) THEN
331  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
332  info = -9
333  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
334  info = -10
335  END IF
336  END IF
337  END IF
338  IF( info.EQ.0 ) THEN
339  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
340  info = -15
341  END IF
342  END IF
343 *
344  IF( info.EQ.0 ) THEN
345  IF( n.LE.1 ) THEN
346  lwkmin = 1
347  work( 1 ) = lwkmin
348  ELSE
349  lwkmin = 2*n
350  nb = ilaenv( 1, 'ZHETRD', uplo, n, -1, -1, -1 )
351  nb = max( nb, ilaenv( 1, 'ZUNMTR', uplo, n, -1, -1, -1 ) )
352  lwkopt = max( 1, ( nb + 1 )*n )
353  work( 1 ) = lwkopt
354  END IF
355 *
356  IF( lwork.LT.lwkmin .AND. .NOT.lquery )
357  $ info = -17
358  END IF
359 *
360  IF( info.NE.0 ) THEN
361  CALL xerbla( 'ZHEEVX', -info )
362  RETURN
363  ELSE IF( lquery ) THEN
364  RETURN
365  END IF
366 *
367 * Quick return if possible
368 *
369  m = 0
370  IF( n.EQ.0 ) THEN
371  RETURN
372  END IF
373 *
374  IF( n.EQ.1 ) THEN
375  IF( alleig .OR. indeig ) THEN
376  m = 1
377  w( 1 ) = a( 1, 1 )
378  ELSE IF( valeig ) THEN
379  IF( vl.LT.dble( a( 1, 1 ) ) .AND. vu.GE.dble( a( 1, 1 ) ) )
380  $ THEN
381  m = 1
382  w( 1 ) = a( 1, 1 )
383  END IF
384  END IF
385  IF( wantz )
386  $ z( 1, 1 ) = cone
387  RETURN
388  END IF
389 *
390 * Get machine constants.
391 *
392  safmin = dlamch( 'Safe minimum' )
393  eps = dlamch( 'Precision' )
394  smlnum = safmin / eps
395  bignum = one / smlnum
396  rmin = sqrt( smlnum )
397  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
398 *
399 * Scale matrix to allowable range, if necessary.
400 *
401  iscale = 0
402  abstll = abstol
403  IF( valeig ) THEN
404  vll = vl
405  vuu = vu
406  END IF
407  anrm = zlanhe( 'M', uplo, n, a, lda, rwork )
408  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
409  iscale = 1
410  sigma = rmin / anrm
411  ELSE IF( anrm.GT.rmax ) THEN
412  iscale = 1
413  sigma = rmax / anrm
414  END IF
415  IF( iscale.EQ.1 ) THEN
416  IF( lower ) THEN
417  DO 10 j = 1, n
418  CALL zdscal( n-j+1, sigma, a( j, j ), 1 )
419  10 CONTINUE
420  ELSE
421  DO 20 j = 1, n
422  CALL zdscal( j, sigma, a( 1, j ), 1 )
423  20 CONTINUE
424  END IF
425  IF( abstol.GT.0 )
426  $ abstll = abstol*sigma
427  IF( valeig ) THEN
428  vll = vl*sigma
429  vuu = vu*sigma
430  END IF
431  END IF
432 *
433 * Call ZHETRD to reduce Hermitian matrix to tridiagonal form.
434 *
435  indd = 1
436  inde = indd + n
437  indrwk = inde + n
438  indtau = 1
439  indwrk = indtau + n
440  llwork = lwork - indwrk + 1
441  CALL zhetrd( uplo, n, a, lda, rwork( indd ), rwork( inde ),
442  $ work( indtau ), work( indwrk ), llwork, iinfo )
443 *
444 * If all eigenvalues are desired and ABSTOL is less than or equal to
445 * zero, then call DSTERF or ZUNGTR and ZSTEQR. If this fails for
446 * some eigenvalue, then try DSTEBZ.
447 *
448  test = .false.
449  IF( indeig ) THEN
450  IF( il.EQ.1 .AND. iu.EQ.n ) THEN
451  test = .true.
452  END IF
453  END IF
454  IF( ( alleig .OR. test ) .AND. ( abstol.LE.zero ) ) THEN
455  CALL dcopy( n, rwork( indd ), 1, w, 1 )
456  indee = indrwk + 2*n
457  IF( .NOT.wantz ) THEN
458  CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
459  CALL dsterf( n, w, rwork( indee ), info )
460  ELSE
461  CALL zlacpy( 'A', n, n, a, lda, z, ldz )
462  CALL zungtr( uplo, n, z, ldz, work( indtau ),
463  $ work( indwrk ), llwork, iinfo )
464  CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
465  CALL zsteqr( jobz, n, w, rwork( indee ), z, ldz,
466  $ rwork( indrwk ), info )
467  IF( info.EQ.0 ) THEN
468  DO 30 i = 1, n
469  ifail( i ) = 0
470  30 CONTINUE
471  END IF
472  END IF
473  IF( info.EQ.0 ) THEN
474  m = n
475  GO TO 40
476  END IF
477  info = 0
478  END IF
479 *
480 * Otherwise, call DSTEBZ and, if eigenvectors are desired, ZSTEIN.
481 *
482  IF( wantz ) THEN
483  order = 'B'
484  ELSE
485  order = 'E'
486  END IF
487  indibl = 1
488  indisp = indibl + n
489  indiwk = indisp + n
490  CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
491  $ rwork( indd ), rwork( inde ), m, nsplit, w,
492  $ iwork( indibl ), iwork( indisp ), rwork( indrwk ),
493  $ iwork( indiwk ), info )
494 *
495  IF( wantz ) THEN
496  CALL zstein( n, rwork( indd ), rwork( inde ), m, w,
497  $ iwork( indibl ), iwork( indisp ), z, ldz,
498  $ rwork( indrwk ), iwork( indiwk ), ifail, info )
499 *
500 * Apply unitary matrix used in reduction to tridiagonal
501 * form to eigenvectors returned by ZSTEIN.
502 *
503  CALL zunmtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
504  $ ldz, work( indwrk ), llwork, iinfo )
505  END IF
506 *
507 * If matrix was scaled, then rescale eigenvalues appropriately.
508 *
509  40 CONTINUE
510  IF( iscale.EQ.1 ) THEN
511  IF( info.EQ.0 ) THEN
512  imax = m
513  ELSE
514  imax = info - 1
515  END IF
516  CALL dscal( imax, one / sigma, w, 1 )
517  END IF
518 *
519 * If eigenvalues are not in order, then sort them, along with
520 * eigenvectors.
521 *
522  IF( wantz ) THEN
523  DO 60 j = 1, m - 1
524  i = 0
525  tmp1 = w( j )
526  DO 50 jj = j + 1, m
527  IF( w( jj ).LT.tmp1 ) THEN
528  i = jj
529  tmp1 = w( jj )
530  END IF
531  50 CONTINUE
532 *
533  IF( i.NE.0 ) THEN
534  itmp1 = iwork( indibl+i-1 )
535  w( i ) = w( j )
536  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
537  w( j ) = tmp1
538  iwork( indibl+j-1 ) = itmp1
539  CALL zswap( n, z( 1, i ), 1, z( 1, j ), 1 )
540  IF( info.NE.0 ) THEN
541  itmp1 = ifail( i )
542  ifail( i ) = ifail( j )
543  ifail( j ) = itmp1
544  END IF
545  END IF
546  60 CONTINUE
547  END IF
548 *
549 * Set WORK(1) to optimal complex workspace size.
550 *
551  work( 1 ) = lwkopt
552 *
553  RETURN
554 *
555 * End of ZHEEVX
556 *
557  END
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:54
subroutine zstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
ZSTEIN
Definition: zstein.f:184
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:88
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zhetrd(UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO)
ZHETRD
Definition: zhetrd.f:194
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
subroutine zsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
ZSTEQR
Definition: zsteqr.f:134
subroutine dstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
DSTEBZ
Definition: dstebz.f:265
subroutine zunmtr(SIDE, UPLO, TRANS, M, N, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMTR
Definition: zunmtr.f:173
subroutine zungtr(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGTR
Definition: zungtr.f:125
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:254