LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 ZSTEMR'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 *> If RANGE='V', the lower bound of the interval to
159 *> be searched for eigenvalues. VL < VU.
160 *> Not referenced if RANGE = 'A' or 'I'.
161 *> \endverbatim
162 *>
163 *> \param[in] VU
164 *> \verbatim
165 *> VU is DOUBLE PRECISION
166 *> If RANGE='V', the upper bound of the interval to
167 *> be searched for eigenvalues. VL < VU.
168 *> Not referenced if RANGE = 'A' or 'I'.
169 *> \endverbatim
170 *>
171 *> \param[in] IL
172 *> \verbatim
173 *> IL is INTEGER
174 *> If RANGE='I', the index of the
175 *> smallest eigenvalue to be returned.
176 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
177 *> Not referenced if RANGE = 'A' or 'V'.
178 *> \endverbatim
179 *>
180 *> \param[in] IU
181 *> \verbatim
182 *> IU is INTEGER
183 *> If RANGE='I', the index of the
184 *> largest eigenvalue to be returned.
185 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
186 *> Not referenced if RANGE = 'A' or 'V'.
187 *> \endverbatim
188 *>
189 *> \param[in] ABSTOL
190 *> \verbatim
191 *> ABSTOL is DOUBLE PRECISION
192 *> The absolute error tolerance for the eigenvalues.
193 *> An approximate eigenvalue is accepted as converged
194 *> when it is determined to lie in an interval [a,b]
195 *> of width less than or equal to
196 *>
197 *> ABSTOL + EPS * max( |a|,|b| ) ,
198 *>
199 *> where EPS is the machine precision. If ABSTOL is less than
200 *> or equal to zero, then EPS*|T| will be used in its place,
201 *> where |T| is the 1-norm of the tridiagonal matrix obtained
202 *> by reducing A to tridiagonal form.
203 *>
204 *> See "Computing Small Singular Values of Bidiagonal Matrices
205 *> with Guaranteed High Relative Accuracy," by Demmel and
206 *> Kahan, LAPACK Working Note #3.
207 *>
208 *> If high relative accuracy is important, set ABSTOL to
209 *> DLAMCH( 'Safe minimum' ). Doing so will guarantee that
210 *> eigenvalues are computed to high relative accuracy when
211 *> possible in future releases. The current code does not
212 *> make any guarantees about high relative accuracy, but
213 *> future releases will. See J. Barlow and J. Demmel,
214 *> "Computing Accurate Eigensystems of Scaled Diagonally
215 *> Dominant Matrices", LAPACK Working Note #7, for a discussion
216 *> of which matrices define their eigenvalues to high relative
217 *> accuracy.
218 *> \endverbatim
219 *>
220 *> \param[out] M
221 *> \verbatim
222 *> M is INTEGER
223 *> The total number of eigenvalues found. 0 <= M <= N.
224 *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
225 *> \endverbatim
226 *>
227 *> \param[out] W
228 *> \verbatim
229 *> W is DOUBLE PRECISION array, dimension (N)
230 *> The first M elements contain the selected eigenvalues in
231 *> ascending order.
232 *> \endverbatim
233 *>
234 *> \param[out] Z
235 *> \verbatim
236 *> Z is COMPLEX*16 array, dimension (LDZ, max(1,M))
237 *> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
238 *> contain the orthonormal eigenvectors of the matrix A
239 *> corresponding to the selected eigenvalues, with the i-th
240 *> column of Z holding the eigenvector associated with W(i).
241 *> If JOBZ = 'N', then Z is not referenced.
242 *> Note: the user must ensure that at least max(1,M) columns are
243 *> supplied in the array Z; if RANGE = 'V', the exact value of M
244 *> is not known in advance and an upper bound must be used.
245 *> \endverbatim
246 *>
247 *> \param[in] LDZ
248 *> \verbatim
249 *> LDZ is INTEGER
250 *> The leading dimension of the array Z. LDZ >= 1, and if
251 *> JOBZ = 'V', LDZ >= max(1,N).
252 *> \endverbatim
253 *>
254 *> \param[out] ISUPPZ
255 *> \verbatim
256 *> ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
257 *> The support of the eigenvectors in Z, i.e., the indices
258 *> indicating the nonzero elements in Z. The i-th eigenvector
259 *> is nonzero only in elements ISUPPZ( 2*i-1 ) through
260 *> ISUPPZ( 2*i ). This is an output of ZSTEMR (tridiagonal
261 *> matrix). The support of the eigenvectors of A is typically
262 *> 1:N because of the unitary transformations applied by ZUNMTR.
263 *> Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
264 *> \endverbatim
265 *>
266 *> \param[out] WORK
267 *> \verbatim
268 *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
269 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
270 *> \endverbatim
271 *>
272 *> \param[in] LWORK
273 *> \verbatim
274 *> LWORK is INTEGER
275 *> The length of the array WORK. LWORK >= max(1,2*N).
276 *> For optimal efficiency, LWORK >= (NB+1)*N,
277 *> where NB is the max of the blocksize for ZHETRD and for
278 *> ZUNMTR as returned by ILAENV.
279 *>
280 *> If LWORK = -1, then a workspace query is assumed; the routine
281 *> only calculates the optimal sizes of the WORK, RWORK and
282 *> IWORK arrays, returns these values as the first entries of
283 *> the WORK, RWORK and IWORK arrays, and no error message
284 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
285 *> \endverbatim
286 *>
287 *> \param[out] RWORK
288 *> \verbatim
289 *> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
290 *> On exit, if INFO = 0, RWORK(1) returns the optimal
291 *> (and minimal) LRWORK.
292 *> \endverbatim
293 *>
294 *> \param[in] LRWORK
295 *> \verbatim
296 *> LRWORK is INTEGER
297 *> The length of the array RWORK. LRWORK >= max(1,24*N).
298 *>
299 *> If LRWORK = -1, then a workspace query is assumed; the
300 *> routine only calculates the optimal sizes of the WORK, RWORK
301 *> and IWORK arrays, returns these values as the first entries
302 *> of the WORK, RWORK and IWORK arrays, and no error message
303 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
304 *> \endverbatim
305 *>
306 *> \param[out] IWORK
307 *> \verbatim
308 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
309 *> On exit, if INFO = 0, IWORK(1) returns the optimal
310 *> (and minimal) LIWORK.
311 *> \endverbatim
312 *>
313 *> \param[in] LIWORK
314 *> \verbatim
315 *> LIWORK is INTEGER
316 *> The dimension of the array IWORK. LIWORK >= max(1,10*N).
317 *>
318 *> If LIWORK = -1, then a workspace query is assumed; the
319 *> routine only calculates the optimal sizes of the WORK, RWORK
320 *> and IWORK arrays, returns these values as the first entries
321 *> of the WORK, RWORK and IWORK arrays, and no error message
322 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
323 *> \endverbatim
324 *>
325 *> \param[out] INFO
326 *> \verbatim
327 *> INFO is INTEGER
328 *> = 0: successful exit
329 *> < 0: if INFO = -i, the i-th argument had an illegal value
330 *> > 0: Internal error
331 *> \endverbatim
332 *
333 * Authors:
334 * ========
335 *
336 *> \author Univ. of Tennessee
337 *> \author Univ. of California Berkeley
338 *> \author Univ. of Colorado Denver
339 *> \author NAG Ltd.
340 *
341 *> \ingroup complex16HEeigen
342 *
343 *> \par Contributors:
344 * ==================
345 *>
346 *> Inderjit Dhillon, IBM Almaden, USA \n
347 *> Osni Marques, LBNL/NERSC, USA \n
348 *> Ken Stanley, Computer Science Division, University of
349 *> California at Berkeley, USA \n
350 *> Jason Riedy, Computer Science Division, University of
351 *> California at Berkeley, USA \n
352 *>
353 * =====================================================================
354  SUBROUTINE zheevr( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU,
355  $ ABSTOL, M, W, Z, LDZ, ISUPPZ, WORK, LWORK,
356  $ RWORK, LRWORK, IWORK, LIWORK, INFO )
357 *
358 * -- LAPACK driver routine --
359 * -- LAPACK is a software package provided by Univ. of Tennessee, --
360 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
361 *
362 * .. Scalar Arguments ..
363  CHARACTER JOBZ, RANGE, UPLO
364  INTEGER IL, INFO, IU, LDA, LDZ, LIWORK, LRWORK, LWORK,
365  $ m, n
366  DOUBLE PRECISION ABSTOL, VL, VU
367 * ..
368 * .. Array Arguments ..
369  INTEGER ISUPPZ( * ), IWORK( * )
370  DOUBLE PRECISION RWORK( * ), W( * )
371  COMPLEX*16 A( LDA, * ), WORK( * ), Z( LDZ, * )
372 * ..
373 *
374 * =====================================================================
375 *
376 * .. Parameters ..
377  DOUBLE PRECISION ZERO, ONE, TWO
378  PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
379 * ..
380 * .. Local Scalars ..
381  LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
382  $ WANTZ, TRYRAC
383  CHARACTER ORDER
384  INTEGER I, IEEEOK, IINFO, IMAX, INDIBL, INDIFL, INDISP,
385  $ indiwo, indrd, indrdd, indre, indree, indrwk,
386  $ indtau, indwk, indwkn, iscale, itmp1, j, jj,
387  $ liwmin, llwork, llrwork, llwrkn, lrwmin,
388  $ lwkopt, lwmin, nb, nsplit
389  DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
390  $ SIGMA, SMLNUM, TMP1, VLL, VUU
391 * ..
392 * .. External Functions ..
393  LOGICAL LSAME
394  INTEGER ILAENV
395  DOUBLE PRECISION DLAMCH, ZLANSY
396  EXTERNAL lsame, ilaenv, dlamch, zlansy
397 * ..
398 * .. External Subroutines ..
399  EXTERNAL dcopy, dscal, dstebz, dsterf, xerbla, zdscal,
401 * ..
402 * .. Intrinsic Functions ..
403  INTRINSIC dble, max, min, sqrt
404 * ..
405 * .. Executable Statements ..
406 *
407 * Test the input parameters.
408 *
409  ieeeok = ilaenv( 10, 'ZHEEVR', 'N', 1, 2, 3, 4 )
410 *
411  lower = lsame( uplo, 'L' )
412  wantz = lsame( jobz, 'V' )
413  alleig = lsame( range, 'A' )
414  valeig = lsame( range, 'V' )
415  indeig = lsame( range, 'I' )
416 *
417  lquery = ( ( lwork.EQ.-1 ) .OR. ( lrwork.EQ.-1 ) .OR.
418  $ ( liwork.EQ.-1 ) )
419 *
420  lrwmin = max( 1, 24*n )
421  liwmin = max( 1, 10*n )
422  lwmin = max( 1, 2*n )
423 *
424  info = 0
425  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
426  info = -1
427  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
428  info = -2
429  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
430  info = -3
431  ELSE IF( n.LT.0 ) THEN
432  info = -4
433  ELSE IF( lda.LT.max( 1, n ) ) THEN
434  info = -6
435  ELSE
436  IF( valeig ) THEN
437  IF( n.GT.0 .AND. vu.LE.vl )
438  $ info = -8
439  ELSE IF( indeig ) THEN
440  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
441  info = -9
442  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
443  info = -10
444  END IF
445  END IF
446  END IF
447  IF( info.EQ.0 ) THEN
448  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
449  info = -15
450  END IF
451  END IF
452 *
453  IF( info.EQ.0 ) THEN
454  nb = ilaenv( 1, 'ZHETRD', uplo, n, -1, -1, -1 )
455  nb = max( nb, ilaenv( 1, 'ZUNMTR', uplo, n, -1, -1, -1 ) )
456  lwkopt = max( ( nb+1 )*n, lwmin )
457  work( 1 ) = lwkopt
458  rwork( 1 ) = lrwmin
459  iwork( 1 ) = liwmin
460 *
461  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
462  info = -18
463  ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
464  info = -20
465  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
466  info = -22
467  END IF
468  END IF
469 *
470  IF( info.NE.0 ) THEN
471  CALL xerbla( 'ZHEEVR', -info )
472  RETURN
473  ELSE IF( lquery ) THEN
474  RETURN
475  END IF
476 *
477 * Quick return if possible
478 *
479  m = 0
480  IF( n.EQ.0 ) THEN
481  work( 1 ) = 1
482  RETURN
483  END IF
484 *
485  IF( n.EQ.1 ) THEN
486  work( 1 ) = 2
487  IF( alleig .OR. indeig ) THEN
488  m = 1
489  w( 1 ) = dble( a( 1, 1 ) )
490  ELSE
491  IF( vl.LT.dble( a( 1, 1 ) ) .AND. vu.GE.dble( a( 1, 1 ) ) )
492  $ THEN
493  m = 1
494  w( 1 ) = dble( a( 1, 1 ) )
495  END IF
496  END IF
497  IF( wantz ) THEN
498  z( 1, 1 ) = one
499  isuppz( 1 ) = 1
500  isuppz( 2 ) = 1
501  END IF
502  RETURN
503  END IF
504 *
505 * Get machine constants.
506 *
507  safmin = dlamch( 'Safe minimum' )
508  eps = dlamch( 'Precision' )
509  smlnum = safmin / eps
510  bignum = one / smlnum
511  rmin = sqrt( smlnum )
512  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
513 *
514 * Scale matrix to allowable range, if necessary.
515 *
516  iscale = 0
517  abstll = abstol
518  IF (valeig) THEN
519  vll = vl
520  vuu = vu
521  END IF
522  anrm = zlansy( 'M', uplo, n, a, lda, rwork )
523  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
524  iscale = 1
525  sigma = rmin / anrm
526  ELSE IF( anrm.GT.rmax ) THEN
527  iscale = 1
528  sigma = rmax / anrm
529  END IF
530  IF( iscale.EQ.1 ) THEN
531  IF( lower ) THEN
532  DO 10 j = 1, n
533  CALL zdscal( n-j+1, sigma, a( j, j ), 1 )
534  10 CONTINUE
535  ELSE
536  DO 20 j = 1, n
537  CALL zdscal( j, sigma, a( 1, j ), 1 )
538  20 CONTINUE
539  END IF
540  IF( abstol.GT.0 )
541  $ abstll = abstol*sigma
542  IF( valeig ) THEN
543  vll = vl*sigma
544  vuu = vu*sigma
545  END IF
546  END IF
547 
548 * Initialize indices into workspaces. Note: The IWORK indices are
549 * used only if DSTERF or ZSTEMR fail.
550 
551 * WORK(INDTAU:INDTAU+N-1) stores the complex scalar factors of the
552 * elementary reflectors used in ZHETRD.
553  indtau = 1
554 * INDWK is the starting offset of the remaining complex workspace,
555 * and LLWORK is the remaining complex workspace size.
556  indwk = indtau + n
557  llwork = lwork - indwk + 1
558 
559 * RWORK(INDRD:INDRD+N-1) stores the real tridiagonal's diagonal
560 * entries.
561  indrd = 1
562 * RWORK(INDRE:INDRE+N-1) stores the off-diagonal entries of the
563 * tridiagonal matrix from ZHETRD.
564  indre = indrd + n
565 * RWORK(INDRDD:INDRDD+N-1) is a copy of the diagonal entries over
566 * -written by ZSTEMR (the DSTERF path copies the diagonal to W).
567  indrdd = indre + n
568 * RWORK(INDREE:INDREE+N-1) is a copy of the off-diagonal entries over
569 * -written while computing the eigenvalues in DSTERF and ZSTEMR.
570  indree = indrdd + n
571 * INDRWK is the starting offset of the left-over real workspace, and
572 * LLRWORK is the remaining workspace size.
573  indrwk = indree + n
574  llrwork = lrwork - indrwk + 1
575 
576 * IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in DSTEBZ and
577 * stores the block indices of each of the M<=N eigenvalues.
578  indibl = 1
579 * IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in DSTEBZ and
580 * stores the starting and finishing indices of each block.
581  indisp = indibl + n
582 * IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
583 * that corresponding to eigenvectors that fail to converge in
584 * DSTEIN. This information is discarded; if any fail, the driver
585 * returns INFO > 0.
586  indifl = indisp + n
587 * INDIWO is the offset of the remaining integer workspace.
588  indiwo = indifl + n
589 
590 *
591 * Call ZHETRD to reduce Hermitian matrix to tridiagonal form.
592 *
593  CALL zhetrd( uplo, n, a, lda, rwork( indrd ), rwork( indre ),
594  $ work( indtau ), work( indwk ), llwork, iinfo )
595 *
596 * If all eigenvalues are desired
597 * then call DSTERF or ZSTEMR and ZUNMTR.
598 *
599  test = .false.
600  IF( indeig ) THEN
601  IF( il.EQ.1 .AND. iu.EQ.n ) THEN
602  test = .true.
603  END IF
604  END IF
605  IF( ( alleig.OR.test ) .AND. ( ieeeok.EQ.1 ) ) THEN
606  IF( .NOT.wantz ) THEN
607  CALL dcopy( n, rwork( indrd ), 1, w, 1 )
608  CALL dcopy( n-1, rwork( indre ), 1, rwork( indree ), 1 )
609  CALL dsterf( n, w, rwork( indree ), info )
610  ELSE
611  CALL dcopy( n-1, rwork( indre ), 1, rwork( indree ), 1 )
612  CALL dcopy( n, rwork( indrd ), 1, rwork( indrdd ), 1 )
613 *
614  IF (abstol .LE. two*n*eps) THEN
615  tryrac = .true.
616  ELSE
617  tryrac = .false.
618  END IF
619  CALL zstemr( jobz, 'A', n, rwork( indrdd ),
620  $ rwork( indree ), vl, vu, il, iu, m, w,
621  $ z, ldz, n, isuppz, tryrac,
622  $ rwork( indrwk ), llrwork,
623  $ iwork, liwork, info )
624 *
625 * Apply unitary matrix used in reduction to tridiagonal
626 * form to eigenvectors returned by ZSTEMR.
627 *
628  IF( wantz .AND. info.EQ.0 ) THEN
629  indwkn = indwk
630  llwrkn = lwork - indwkn + 1
631  CALL zunmtr( 'L', uplo, 'N', n, m, a, lda,
632  $ work( indtau ), z, ldz, work( indwkn ),
633  $ llwrkn, iinfo )
634  END IF
635  END IF
636 *
637 *
638  IF( info.EQ.0 ) THEN
639  m = n
640  GO TO 30
641  END IF
642  info = 0
643  END IF
644 *
645 * Otherwise, call DSTEBZ and, if eigenvectors are desired, ZSTEIN.
646 * Also call DSTEBZ and ZSTEIN if ZSTEMR fails.
647 *
648  IF( wantz ) THEN
649  order = 'B'
650  ELSE
651  order = 'E'
652  END IF
653 
654  CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
655  $ rwork( indrd ), rwork( indre ), m, nsplit, w,
656  $ iwork( indibl ), iwork( indisp ), rwork( indrwk ),
657  $ iwork( indiwo ), info )
658 *
659  IF( wantz ) THEN
660  CALL zstein( n, rwork( indrd ), rwork( indre ), m, w,
661  $ iwork( indibl ), iwork( indisp ), z, ldz,
662  $ rwork( indrwk ), iwork( indiwo ), iwork( indifl ),
663  $ info )
664 *
665 * Apply unitary matrix used in reduction to tridiagonal
666 * form to eigenvectors returned by ZSTEIN.
667 *
668  indwkn = indwk
669  llwrkn = lwork - indwkn + 1
670  CALL zunmtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
671  $ ldz, work( indwkn ), llwrkn, iinfo )
672  END IF
673 *
674 * If matrix was scaled, then rescale eigenvalues appropriately.
675 *
676  30 CONTINUE
677  IF( iscale.EQ.1 ) THEN
678  IF( info.EQ.0 ) THEN
679  imax = m
680  ELSE
681  imax = info - 1
682  END IF
683  CALL dscal( imax, one / sigma, w, 1 )
684  END IF
685 *
686 * If eigenvalues are not in order, then sort them, along with
687 * eigenvectors.
688 *
689  IF( wantz ) THEN
690  DO 50 j = 1, m - 1
691  i = 0
692  tmp1 = w( j )
693  DO 40 jj = j + 1, m
694  IF( w( jj ).LT.tmp1 ) THEN
695  i = jj
696  tmp1 = w( jj )
697  END IF
698  40 CONTINUE
699 *
700  IF( i.NE.0 ) THEN
701  itmp1 = iwork( indibl+i-1 )
702  w( i ) = w( j )
703  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
704  w( j ) = tmp1
705  iwork( indibl+j-1 ) = itmp1
706  CALL zswap( n, z( 1, i ), 1, z( 1, j ), 1 )
707  END IF
708  50 CONTINUE
709  END IF
710 *
711 * Set WORK(1) to optimal workspace size.
712 *
713  work( 1 ) = lwkopt
714  rwork( 1 ) = lrwmin
715  iwork( 1 ) = liwmin
716 *
717  RETURN
718 *
719 * End of ZHEEVR
720 *
721  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
DSTEBZ
Definition: dstebz.f:273
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:86
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:78
subroutine zhetrd(UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO)
ZHETRD
Definition: zhetrd.f:192
subroutine zheevr(JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, ISUPPZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
ZHEEVR computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices
Definition: zheevr.f:357
subroutine zstemr(JOBZ, RANGE, N, D, E, VL, VU, IL, IU, M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK, IWORK, LIWORK, INFO)
ZSTEMR
Definition: zstemr.f:338
subroutine zunmtr(SIDE, UPLO, TRANS, M, N, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMTR
Definition: zunmtr.f:171
subroutine zstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
ZSTEIN
Definition: zstein.f:182
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79