LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
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
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:51
LOGICAL function lsame(CA, CB)
LSAME
Definition: lsame.f:54
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:53
subroutine zstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
ZSTEIN
Definition: zstein.f:182
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:87
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:52
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:61
subroutine zhetrd(UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO)
ZHETRD
Definition: zhetrd.f:193
DOUBLE PRECISION function zlansy(NORM, UPLO, N, A, LDA, WORK)
ZLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex symmetric matrix.
Definition: zlansy.f:124
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:54
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:328
INTEGER function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:81
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
Definition: xerbla-fortran:9
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:64
subroutine dstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
DSTEBZ
Definition: dstebz.f:262
subroutine zunmtr(SIDE, UPLO, TRANS, M, N, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMTR
Definition: zunmtr.f:171
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:347