LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zheevr.f
Go to the documentation of this file.
1 *> \brief <b> ZHEEVR 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 ZHEEVR + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zheevr.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zheevr.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zheevr.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZHEEVR( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU,
22 * ABSTOL, M, W, Z, LDZ, ISUPPZ, WORK, LWORK,
23 * RWORK, LRWORK, IWORK, LIWORK, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER JOBZ, RANGE, UPLO
27 * INTEGER IL, INFO, IU, LDA, LDZ, LIWORK, LRWORK, LWORK,
28 * $ M, N
29 * DOUBLE PRECISION ABSTOL, VL, VU
30 * ..
31 * .. Array Arguments ..
32 * INTEGER ISUPPZ( * ), IWORK( * )
33 * DOUBLE PRECISION RWORK( * ), W( * )
34 * COMPLEX*16 A( LDA, * ), WORK( * ), Z( LDZ, * )
35 * ..
36 *
37 *
38 *> \par Purpose:
39 * =============
40 *>
41 *> \verbatim
42 *>
43 *> ZHEEVR computes selected eigenvalues and, optionally, eigenvectors
44 *> of a complex Hermitian matrix A. Eigenvalues and eigenvectors can
45 *> be selected by specifying either a range of values or a range of
46 *> indices for the desired eigenvalues.
47 *>
48 *> ZHEEVR first reduces the matrix A to tridiagonal form T with a call
49 *> to ZHETRD. Then, whenever possible, ZHEEVR calls ZSTEMR to compute
50 *> eigenspectrum using Relatively Robust Representations. ZSTEMR
51 *> computes eigenvalues by the dqds algorithm, while orthogonal
52 *> eigenvectors are computed from various "good" L D L^T representations
53 *> (also known as Relatively Robust Representations). Gram-Schmidt
54 *> orthogonalization is avoided as far as possible. More specifically,
55 *> the various steps of the algorithm are as follows.
56 *>
57 *> For each unreduced block (submatrix) of T,
58 *> (a) Compute T - sigma I = L D L^T, so that L and D
59 *> define all the wanted eigenvalues to high relative accuracy.
60 *> This means that small relative changes in the entries of D and L
61 *> cause only small relative changes in the eigenvalues and
62 *> eigenvectors. The standard (unfactored) representation of the
63 *> tridiagonal matrix T does not have this property in general.
64 *> (b) Compute the eigenvalues to suitable accuracy.
65 *> If the eigenvectors are desired, the algorithm attains full
66 *> accuracy of the computed eigenvalues only right before
67 *> the corresponding vectors have to be computed, see steps c) and d).
68 *> (c) For each cluster of close eigenvalues, select a new
69 *> shift close to the cluster, find a new factorization, and refine
70 *> the shifted eigenvalues to suitable accuracy.
71 *> (d) For each eigenvalue with a large enough relative separation compute
72 *> the corresponding eigenvector by forming a rank revealing twisted
73 *> factorization. Go back to (c) for any clusters that remain.
74 *>
75 *> The desired accuracy of the output can be specified by the input
76 *> parameter ABSTOL.
77 *>
78 *> For more details, see DSTEMR's documentation and:
79 *> - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations
80 *> to compute orthogonal eigenvectors of symmetric tridiagonal matrices,"
81 *> Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004.
82 *> - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and
83 *> Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25,
84 *> 2004. Also LAPACK Working Note 154.
85 *> - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric
86 *> tridiagonal eigenvalue/eigenvector problem",
87 *> Computer Science Division Technical Report No. UCB/CSD-97-971,
88 *> UC Berkeley, May 1997.
89 *>
90 *>
91 *> Note 1 : ZHEEVR calls ZSTEMR when the full spectrum is requested
92 *> on machines which conform to the ieee-754 floating point standard.
93 *> ZHEEVR calls DSTEBZ and ZSTEIN on non-ieee machines and
94 *> when partial spectrum requests are made.
95 *>
96 *> Normal execution of ZSTEMR may create NaNs and infinities and
97 *> hence may abort due to a floating point exception in environments
98 *> which do not handle NaNs and infinities in the ieee standard default
99 *> manner.
100 *> \endverbatim
101 *
102 * Arguments:
103 * ==========
104 *
105 *> \param[in] JOBZ
106 *> \verbatim
107 *> JOBZ is CHARACTER*1
108 *> = 'N': Compute eigenvalues only;
109 *> = 'V': Compute eigenvalues and eigenvectors.
110 *> \endverbatim
111 *>
112 *> \param[in] RANGE
113 *> \verbatim
114 *> RANGE is CHARACTER*1
115 *> = 'A': all eigenvalues will be found.
116 *> = 'V': all eigenvalues in the half-open interval (VL,VU]
117 *> will be found.
118 *> = 'I': the IL-th through IU-th eigenvalues will be found.
119 *> For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and
120 *> ZSTEIN are called
121 *> \endverbatim
122 *>
123 *> \param[in] UPLO
124 *> \verbatim
125 *> UPLO is CHARACTER*1
126 *> = 'U': Upper triangle of A is stored;
127 *> = 'L': Lower triangle of A is stored.
128 *> \endverbatim
129 *>
130 *> \param[in] N
131 *> \verbatim
132 *> N is INTEGER
133 *> The order of the matrix A. N >= 0.
134 *> \endverbatim
135 *>
136 *> \param[in,out] A
137 *> \verbatim
138 *> A is COMPLEX*16 array, dimension (LDA, N)
139 *> On entry, the Hermitian matrix A. If UPLO = 'U', the
140 *> leading N-by-N upper triangular part of A contains the
141 *> upper triangular part of the matrix A. If UPLO = 'L',
142 *> the leading N-by-N lower triangular part of A contains
143 *> the lower triangular part of the matrix A.
144 *> On exit, the lower triangle (if UPLO='L') or the upper
145 *> triangle (if UPLO='U') of A, including the diagonal, is
146 *> destroyed.
147 *> \endverbatim
148 *>
149 *> \param[in] LDA
150 *> \verbatim
151 *> LDA is INTEGER
152 *> The leading dimension of the array A. LDA >= max(1,N).
153 *> \endverbatim
154 *>
155 *> \param[in] VL
156 *> \verbatim
157 *> VL is DOUBLE PRECISION
158 *> \endverbatim
159 *>
160 *> \param[in] VU
161 *> \verbatim
162 *> VU is DOUBLE PRECISION
163 *> If RANGE='V', the lower and upper bounds of the interval to
164 *> be searched for eigenvalues. VL < VU.
165 *> Not referenced if RANGE = 'A' or 'I'.
166 *> \endverbatim
167 *>
168 *> \param[in] IL
169 *> \verbatim
170 *> IL is INTEGER
171 *> \endverbatim
172 *>
173 *> \param[in] IU
174 *> \verbatim
175 *> IU is INTEGER
176 *> If RANGE='I', the indices (in ascending order) of the
177 *> smallest and largest eigenvalues to be returned.
178 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
179 *> Not referenced if RANGE = 'A' or 'V'.
180 *> \endverbatim
181 *>
182 *> \param[in] ABSTOL
183 *> \verbatim
184 *> ABSTOL is DOUBLE PRECISION
185 *> The absolute error tolerance for the eigenvalues.
186 *> An approximate eigenvalue is accepted as converged
187 *> when it is determined to lie in an interval [a,b]
188 *> of width less than or equal to
189 *>
190 *> ABSTOL + EPS * max( |a|,|b| ) ,
191 *>
192 *> where EPS is the machine precision. If ABSTOL is less than
193 *> or equal to zero, then EPS*|T| will be used in its place,
194 *> where |T| is the 1-norm of the tridiagonal matrix obtained
195 *> by reducing A to tridiagonal form.
196 *>
197 *> See "Computing Small Singular Values of Bidiagonal Matrices
198 *> with Guaranteed High Relative Accuracy," by Demmel and
199 *> Kahan, LAPACK Working Note #3.
200 *>
201 *> If high relative accuracy is important, set ABSTOL to
202 *> DLAMCH( 'Safe minimum' ). Doing so will guarantee that
203 *> eigenvalues are computed to high relative accuracy when
204 *> possible in future releases. The current code does not
205 *> make any guarantees about high relative accuracy, but
206 *> furutre releases will. See J. Barlow and J. Demmel,
207 *> "Computing Accurate Eigensystems of Scaled Diagonally
208 *> Dominant Matrices", LAPACK Working Note #7, for a discussion
209 *> of which matrices define their eigenvalues to high relative
210 *> accuracy.
211 *> \endverbatim
212 *>
213 *> \param[out] M
214 *> \verbatim
215 *> M is INTEGER
216 *> The total number of eigenvalues found. 0 <= M <= N.
217 *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
218 *> \endverbatim
219 *>
220 *> \param[out] W
221 *> \verbatim
222 *> W is DOUBLE PRECISION array, dimension (N)
223 *> The first M elements contain the selected eigenvalues in
224 *> ascending order.
225 *> \endverbatim
226 *>
227 *> \param[out] Z
228 *> \verbatim
229 *> Z is COMPLEX*16 array, dimension (LDZ, max(1,M))
230 *> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
231 *> contain the orthonormal eigenvectors of the matrix A
232 *> corresponding to the selected eigenvalues, with the i-th
233 *> column of Z holding the eigenvector associated with W(i).
234 *> If JOBZ = 'N', then Z is not referenced.
235 *> Note: the user must ensure that at least max(1,M) columns are
236 *> supplied in the array Z; if RANGE = 'V', the exact value of M
237 *> is not known in advance and an upper bound must be used.
238 *> \endverbatim
239 *>
240 *> \param[in] LDZ
241 *> \verbatim
242 *> LDZ is INTEGER
243 *> The leading dimension of the array Z. LDZ >= 1, and if
244 *> JOBZ = 'V', LDZ >= max(1,N).
245 *> \endverbatim
246 *>
247 *> \param[out] ISUPPZ
248 *> \verbatim
249 *> ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
250 *> The support of the eigenvectors in Z, i.e., the indices
251 *> indicating the nonzero elements in Z. The i-th eigenvector
252 *> is nonzero only in elements ISUPPZ( 2*i-1 ) through
253 *> ISUPPZ( 2*i ).
254 *> Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
255 *> \endverbatim
256 *>
257 *> \param[out] WORK
258 *> \verbatim
259 *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
260 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
261 *> \endverbatim
262 *>
263 *> \param[in] LWORK
264 *> \verbatim
265 *> LWORK is INTEGER
266 *> The length of the array WORK. LWORK >= max(1,2*N).
267 *> For optimal efficiency, LWORK >= (NB+1)*N,
268 *> where NB is the max of the blocksize for ZHETRD and for
269 *> ZUNMTR as returned by ILAENV.
270 *>
271 *> If LWORK = -1, then a workspace query is assumed; the routine
272 *> only calculates the optimal sizes of the WORK, RWORK and
273 *> IWORK arrays, returns these values as the first entries of
274 *> the WORK, RWORK and IWORK arrays, and no error message
275 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
276 *> \endverbatim
277 *>
278 *> \param[out] RWORK
279 *> \verbatim
280 *> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
281 *> On exit, if INFO = 0, RWORK(1) returns the optimal
282 *> (and minimal) LRWORK.
283 *> \endverbatim
284 *>
285 *> \param[in] LRWORK
286 *> \verbatim
287 *> LRWORK is INTEGER
288 *> The length of the array RWORK. LRWORK >= max(1,24*N).
289 *>
290 *> If LRWORK = -1, then a workspace query is assumed; the
291 *> routine only calculates the optimal sizes of the WORK, RWORK
292 *> and IWORK arrays, returns these values as the first entries
293 *> of the WORK, RWORK and IWORK arrays, and no error message
294 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
295 *> \endverbatim
296 *>
297 *> \param[out] IWORK
298 *> \verbatim
299 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
300 *> On exit, if INFO = 0, IWORK(1) returns the optimal
301 *> (and minimal) LIWORK.
302 *> \endverbatim
303 *>
304 *> \param[in] LIWORK
305 *> \verbatim
306 *> LIWORK is INTEGER
307 *> The dimension of the array IWORK. LIWORK >= max(1,10*N).
308 *>
309 *> If LIWORK = -1, then a workspace query is assumed; the
310 *> routine only calculates the optimal sizes of the WORK, RWORK
311 *> and IWORK arrays, returns these values as the first entries
312 *> of the WORK, RWORK and IWORK arrays, and no error message
313 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
314 *> \endverbatim
315 *>
316 *> \param[out] INFO
317 *> \verbatim
318 *> INFO is INTEGER
319 *> = 0: successful exit
320 *> < 0: if INFO = -i, the i-th argument had an illegal value
321 *> > 0: Internal error
322 *> \endverbatim
323 *
324 * Authors:
325 * ========
326 *
327 *> \author Univ. of Tennessee
328 *> \author Univ. of California Berkeley
329 *> \author Univ. of Colorado Denver
330 *> \author NAG Ltd.
331 *
332 *> \date September 2012
333 *
334 *> \ingroup complex16HEeigen
335 *
336 *> \par Contributors:
337 * ==================
338 *>
339 *> Inderjit Dhillon, IBM Almaden, USA \n
340 *> Osni Marques, LBNL/NERSC, USA \n
341 *> Ken Stanley, Computer Science Division, University of
342 *> California at Berkeley, USA \n
343 *> Jason Riedy, Computer Science Division, University of
344 *> California at Berkeley, USA \n
345 *>
346 * =====================================================================
347  SUBROUTINE zheevr( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU,
348  $ abstol, m, w, z, ldz, isuppz, work, lwork,
349  $ rwork, lrwork, iwork, liwork, info )
350 *
351 * -- LAPACK driver routine (version 3.4.2) --
352 * -- LAPACK is a software package provided by Univ. of Tennessee, --
353 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
354 * September 2012
355 *
356 * .. Scalar Arguments ..
357  CHARACTER jobz, range, uplo
358  INTEGER il, info, iu, lda, ldz, liwork, lrwork, lwork,
359  $ m, n
360  DOUBLE PRECISION abstol, vl, vu
361 * ..
362 * .. Array Arguments ..
363  INTEGER isuppz( * ), iwork( * )
364  DOUBLE PRECISION rwork( * ), w( * )
365  COMPLEX*16 a( lda, * ), work( * ), z( ldz, * )
366 * ..
367 *
368 * =====================================================================
369 *
370 * .. Parameters ..
371  DOUBLE PRECISION zero, one, two
372  parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
373 * ..
374 * .. Local Scalars ..
375  LOGICAL alleig, indeig, lower, lquery, test, valeig,
376  $ wantz, tryrac
377  CHARACTER order
378  INTEGER i, ieeeok, iinfo, imax, indibl, indifl, indisp,
379  $ indiwo, indrd, indrdd, indre, indree, indrwk,
380  $ indtau, indwk, indwkn, iscale, itmp1, j, jj,
381  $ liwmin, llwork, llrwork, llwrkn, lrwmin,
382  $ lwkopt, lwmin, nb, nsplit
383  DOUBLE PRECISION abstll, anrm, bignum, eps, rmax, rmin, safmin,
384  $ sigma, smlnum, tmp1, vll, vuu
385 * ..
386 * .. External Functions ..
387  LOGICAL lsame
388  INTEGER ilaenv
389  DOUBLE PRECISION dlamch, zlansy
390  EXTERNAL lsame, ilaenv, dlamch, zlansy
391 * ..
392 * .. External Subroutines ..
393  EXTERNAL dcopy, dscal, dstebz, dsterf, xerbla, zdscal,
395 * ..
396 * .. Intrinsic Functions ..
397  INTRINSIC dble, max, min, sqrt
398 * ..
399 * .. Executable Statements ..
400 *
401 * Test the input parameters.
402 *
403  ieeeok = ilaenv( 10, 'ZHEEVR', 'N', 1, 2, 3, 4 )
404 *
405  lower = lsame( uplo, 'L' )
406  wantz = lsame( jobz, 'V' )
407  alleig = lsame( range, 'A' )
408  valeig = lsame( range, 'V' )
409  indeig = lsame( range, 'I' )
410 *
411  lquery = ( ( lwork.EQ.-1 ) .OR. ( lrwork.EQ.-1 ) .OR.
412  $ ( liwork.EQ.-1 ) )
413 *
414  lrwmin = max( 1, 24*n )
415  liwmin = max( 1, 10*n )
416  lwmin = max( 1, 2*n )
417 *
418  info = 0
419  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
420  info = -1
421  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
422  info = -2
423  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
424  info = -3
425  ELSE IF( n.LT.0 ) THEN
426  info = -4
427  ELSE IF( lda.LT.max( 1, n ) ) THEN
428  info = -6
429  ELSE
430  IF( valeig ) THEN
431  IF( n.GT.0 .AND. vu.LE.vl )
432  $ info = -8
433  ELSE IF( indeig ) THEN
434  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
435  info = -9
436  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
437  info = -10
438  END IF
439  END IF
440  END IF
441  IF( info.EQ.0 ) THEN
442  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
443  info = -15
444  END IF
445  END IF
446 *
447  IF( info.EQ.0 ) THEN
448  nb = ilaenv( 1, 'ZHETRD', uplo, n, -1, -1, -1 )
449  nb = max( nb, ilaenv( 1, 'ZUNMTR', uplo, n, -1, -1, -1 ) )
450  lwkopt = max( ( nb+1 )*n, lwmin )
451  work( 1 ) = lwkopt
452  rwork( 1 ) = lrwmin
453  iwork( 1 ) = liwmin
454 *
455  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
456  info = -18
457  ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
458  info = -20
459  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
460  info = -22
461  END IF
462  END IF
463 *
464  IF( info.NE.0 ) THEN
465  CALL xerbla( 'ZHEEVR', -info )
466  return
467  ELSE IF( lquery ) THEN
468  return
469  END IF
470 *
471 * Quick return if possible
472 *
473  m = 0
474  IF( n.EQ.0 ) THEN
475  work( 1 ) = 1
476  return
477  END IF
478 *
479  IF( n.EQ.1 ) THEN
480  work( 1 ) = 2
481  IF( alleig .OR. indeig ) THEN
482  m = 1
483  w( 1 ) = dble( a( 1, 1 ) )
484  ELSE
485  IF( vl.LT.dble( a( 1, 1 ) ) .AND. vu.GE.dble( a( 1, 1 ) ) )
486  $ THEN
487  m = 1
488  w( 1 ) = dble( a( 1, 1 ) )
489  END IF
490  END IF
491  IF( wantz ) THEN
492  z( 1, 1 ) = one
493  isuppz( 1 ) = 1
494  isuppz( 2 ) = 1
495  END IF
496  return
497  END IF
498 *
499 * Get machine constants.
500 *
501  safmin = dlamch( 'Safe minimum' )
502  eps = dlamch( 'Precision' )
503  smlnum = safmin / eps
504  bignum = one / smlnum
505  rmin = sqrt( smlnum )
506  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
507 *
508 * Scale matrix to allowable range, if necessary.
509 *
510  iscale = 0
511  abstll = abstol
512  IF (valeig) THEN
513  vll = vl
514  vuu = vu
515  END IF
516  anrm = zlansy( 'M', uplo, n, a, lda, rwork )
517  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
518  iscale = 1
519  sigma = rmin / anrm
520  ELSE IF( anrm.GT.rmax ) THEN
521  iscale = 1
522  sigma = rmax / anrm
523  END IF
524  IF( iscale.EQ.1 ) THEN
525  IF( lower ) THEN
526  DO 10 j = 1, n
527  CALL zdscal( n-j+1, sigma, a( j, j ), 1 )
528  10 continue
529  ELSE
530  DO 20 j = 1, n
531  CALL zdscal( j, sigma, a( 1, j ), 1 )
532  20 continue
533  END IF
534  IF( abstol.GT.0 )
535  $ abstll = abstol*sigma
536  IF( valeig ) THEN
537  vll = vl*sigma
538  vuu = vu*sigma
539  END IF
540  END IF
541 
542 * Initialize indices into workspaces. Note: The IWORK indices are
543 * used only if DSTERF or ZSTEMR fail.
544 
545 * WORK(INDTAU:INDTAU+N-1) stores the complex scalar factors of the
546 * elementary reflectors used in ZHETRD.
547  indtau = 1
548 * INDWK is the starting offset of the remaining complex workspace,
549 * and LLWORK is the remaining complex workspace size.
550  indwk = indtau + n
551  llwork = lwork - indwk + 1
552 
553 * RWORK(INDRD:INDRD+N-1) stores the real tridiagonal's diagonal
554 * entries.
555  indrd = 1
556 * RWORK(INDRE:INDRE+N-1) stores the off-diagonal entries of the
557 * tridiagonal matrix from ZHETRD.
558  indre = indrd + n
559 * RWORK(INDRDD:INDRDD+N-1) is a copy of the diagonal entries over
560 * -written by ZSTEMR (the DSTERF path copies the diagonal to W).
561  indrdd = indre + n
562 * RWORK(INDREE:INDREE+N-1) is a copy of the off-diagonal entries over
563 * -written while computing the eigenvalues in DSTERF and ZSTEMR.
564  indree = indrdd + n
565 * INDRWK is the starting offset of the left-over real workspace, and
566 * LLRWORK is the remaining workspace size.
567  indrwk = indree + n
568  llrwork = lrwork - indrwk + 1
569 
570 * IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in DSTEBZ and
571 * stores the block indices of each of the M<=N eigenvalues.
572  indibl = 1
573 * IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in DSTEBZ and
574 * stores the starting and finishing indices of each block.
575  indisp = indibl + n
576 * IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
577 * that corresponding to eigenvectors that fail to converge in
578 * DSTEIN. This information is discarded; if any fail, the driver
579 * returns INFO > 0.
580  indifl = indisp + n
581 * INDIWO is the offset of the remaining integer workspace.
582  indiwo = indifl + n
583 
584 *
585 * Call ZHETRD to reduce Hermitian matrix to tridiagonal form.
586 *
587  CALL zhetrd( uplo, n, a, lda, rwork( indrd ), rwork( indre ),
588  $ work( indtau ), work( indwk ), llwork, iinfo )
589 *
590 * If all eigenvalues are desired
591 * then call DSTERF or ZSTEMR and ZUNMTR.
592 *
593  test = .false.
594  IF( indeig ) THEN
595  IF( il.EQ.1 .AND. iu.EQ.n ) THEN
596  test = .true.
597  END IF
598  END IF
599  IF( ( alleig.OR.test ) .AND. ( ieeeok.EQ.1 ) ) THEN
600  IF( .NOT.wantz ) THEN
601  CALL dcopy( n, rwork( indrd ), 1, w, 1 )
602  CALL dcopy( n-1, rwork( indre ), 1, rwork( indree ), 1 )
603  CALL dsterf( n, w, rwork( indree ), info )
604  ELSE
605  CALL dcopy( n-1, rwork( indre ), 1, rwork( indree ), 1 )
606  CALL dcopy( n, rwork( indrd ), 1, rwork( indrdd ), 1 )
607 *
608  IF (abstol .LE. two*n*eps) THEN
609  tryrac = .true.
610  ELSE
611  tryrac = .false.
612  END IF
613  CALL zstemr( jobz, 'A', n, rwork( indrdd ),
614  $ rwork( indree ), vl, vu, il, iu, m, w,
615  $ z, ldz, n, isuppz, tryrac,
616  $ rwork( indrwk ), llrwork,
617  $ iwork, liwork, info )
618 *
619 * Apply unitary matrix used in reduction to tridiagonal
620 * form to eigenvectors returned by ZSTEIN.
621 *
622  IF( wantz .AND. info.EQ.0 ) THEN
623  indwkn = indwk
624  llwrkn = lwork - indwkn + 1
625  CALL zunmtr( 'L', uplo, 'N', n, m, a, lda,
626  $ work( indtau ), z, ldz, work( indwkn ),
627  $ llwrkn, iinfo )
628  END IF
629  END IF
630 *
631 *
632  IF( info.EQ.0 ) THEN
633  m = n
634  go to 30
635  END IF
636  info = 0
637  END IF
638 *
639 * Otherwise, call DSTEBZ and, if eigenvectors are desired, ZSTEIN.
640 * Also call DSTEBZ and ZSTEIN if ZSTEMR fails.
641 *
642  IF( wantz ) THEN
643  order = 'B'
644  ELSE
645  order = 'E'
646  END IF
647 
648  CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
649  $ rwork( indrd ), rwork( indre ), m, nsplit, w,
650  $ iwork( indibl ), iwork( indisp ), rwork( indrwk ),
651  $ iwork( indiwo ), info )
652 *
653  IF( wantz ) THEN
654  CALL zstein( n, rwork( indrd ), rwork( indre ), m, w,
655  $ iwork( indibl ), iwork( indisp ), z, ldz,
656  $ rwork( indrwk ), iwork( indiwo ), iwork( indifl ),
657  $ info )
658 *
659 * Apply unitary matrix used in reduction to tridiagonal
660 * form to eigenvectors returned by ZSTEIN.
661 *
662  indwkn = indwk
663  llwrkn = lwork - indwkn + 1
664  CALL zunmtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
665  $ ldz, work( indwkn ), llwrkn, iinfo )
666  END IF
667 *
668 * If matrix was scaled, then rescale eigenvalues appropriately.
669 *
670  30 continue
671  IF( iscale.EQ.1 ) THEN
672  IF( info.EQ.0 ) THEN
673  imax = m
674  ELSE
675  imax = info - 1
676  END IF
677  CALL dscal( imax, one / sigma, w, 1 )
678  END IF
679 *
680 * If eigenvalues are not in order, then sort them, along with
681 * eigenvectors.
682 *
683  IF( wantz ) THEN
684  DO 50 j = 1, m - 1
685  i = 0
686  tmp1 = w( j )
687  DO 40 jj = j + 1, m
688  IF( w( jj ).LT.tmp1 ) THEN
689  i = jj
690  tmp1 = w( jj )
691  END IF
692  40 continue
693 *
694  IF( i.NE.0 ) THEN
695  itmp1 = iwork( indibl+i-1 )
696  w( i ) = w( j )
697  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
698  w( j ) = tmp1
699  iwork( indibl+j-1 ) = itmp1
700  CALL zswap( n, z( 1, i ), 1, z( 1, j ), 1 )
701  END IF
702  50 continue
703  END IF
704 *
705 * Set WORK(1) to optimal workspace size.
706 *
707  work( 1 ) = lwkopt
708  rwork( 1 ) = lrwmin
709  iwork( 1 ) = liwmin
710 *
711  return
712 *
713 * End of ZHEEVR
714 *
715  END