LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
cheevr.f
Go to the documentation of this file.
1 *> \brief <b> CHEEVR 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 CHEEVR + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cheevr.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cheevr.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cheevr.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CHEEVR( 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 * REAL ABSTOL, VL, VU
30 * ..
31 * .. Array Arguments ..
32 * INTEGER ISUPPZ( * ), IWORK( * )
33 * REAL RWORK( * ), W( * )
34 * COMPLEX A( LDA, * ), WORK( * ), Z( LDZ, * )
35 * ..
36 *
37 *
38 *> \par Purpose:
39 * =============
40 *>
41 *> \verbatim
42 *>
43 *> CHEEVR 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 *> CHEEVR first reduces the matrix A to tridiagonal form T with a call
49 *> to CHETRD. Then, whenever possible, CHEEVR calls CSTEMR to compute
50 *> the eigenspectrum using Relatively Robust Representations. CSTEMR
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 : CHEEVR calls CSTEMR when the full spectrum is requested
92 *> on machines which conform to the ieee-754 floating point standard.
93 *> CHEEVR calls SSTEBZ and CSTEIN on non-ieee machines and
94 *> when partial spectrum requests are made.
95 *>
96 *> Normal execution of CSTEMR 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, SSTEBZ and
120 *> CSTEIN 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 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 REAL
158 *> \endverbatim
159 *>
160 *> \param[in] VU
161 *> \verbatim
162 *> VU is REAL
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 REAL
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 *> SLAMCH( '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 REAL 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 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 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 CHETRD and for
269 *> CUNMTR 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 REAL 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 complexHEeigen
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 cheevr( 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  REAL abstol, vl, vu
361 * ..
362 * .. Array Arguments ..
363  INTEGER isuppz( * ), iwork( * )
364  REAL rwork( * ), w( * )
365  COMPLEX a( lda, * ), work( * ), z( ldz, * )
366 * ..
367 *
368 * =====================================================================
369 *
370 * .. Parameters ..
371  REAL zero, one, two
372  parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+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  REAL abstll, anrm, bignum, eps, rmax, rmin, safmin,
384  $ sigma, smlnum, tmp1, vll, vuu
385 * ..
386 * .. External Functions ..
387  LOGICAL lsame
388  INTEGER ilaenv
389  REAL clansy, slamch
390  EXTERNAL lsame, ilaenv, clansy, slamch
391 * ..
392 * .. External Subroutines ..
393  EXTERNAL chetrd, csscal, cstemr, cstein, cswap, cunmtr,
395 * ..
396 * .. Intrinsic Functions ..
397  INTRINSIC max, min, REAL, sqrt
398 * ..
399 * .. Executable Statements ..
400 *
401 * Test the input parameters.
402 *
403  ieeeok = ilaenv( 10, 'CHEEVR', '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, 'CHETRD', uplo, n, -1, -1, -1 )
449  nb = max( nb, ilaenv( 1, 'CUNMTR', 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( 'CHEEVR', -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 ) = REAL( A( 1, 1 ) )
484  ELSE
485  IF( vl.LT.REAL( A( 1, 1 ) ) .AND. vu.GE.REAL( A( 1, 1 ) ) )
486  $ THEN
487  m = 1
488  w( 1 ) = REAL( 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 = slamch( 'Safe minimum' )
502  eps = slamch( '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 = clansy( '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 csscal( n-j+1, sigma, a( j, j ), 1 )
528  10 CONTINUE
529  ELSE
530  DO 20 j = 1, n
531  CALL csscal( 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 SSTERF or CSTEMR fail.
544 
545 * WORK(INDTAU:INDTAU+N-1) stores the complex scalar factors of the
546 * elementary reflectors used in CHETRD.
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 CHETRD.
558  indre = indrd + n
559 * RWORK(INDRDD:INDRDD+N-1) is a copy of the diagonal entries over
560 * -written by CSTEMR (the SSTERF 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 SSTERF and CSTEMR.
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 SSTEBZ 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 SSTEBZ 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 * SSTEIN. 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 CHETRD to reduce Hermitian matrix to tridiagonal form.
586 *
587  CALL chetrd( 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 SSTERF or CSTEMR and CUNMTR.
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 scopy( n, rwork( indrd ), 1, w, 1 )
602  CALL scopy( n-1, rwork( indre ), 1, rwork( indree ), 1 )
603  CALL ssterf( n, w, rwork( indree ), info )
604  ELSE
605  CALL scopy( n-1, rwork( indre ), 1, rwork( indree ), 1 )
606  CALL scopy( 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 cstemr( 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 CSTEIN.
621 *
622  IF( wantz .AND. info.EQ.0 ) THEN
623  indwkn = indwk
624  llwrkn = lwork - indwkn + 1
625  CALL cunmtr( '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 SSTEBZ and, if eigenvectors are desired, CSTEIN.
640 * Also call SSTEBZ and CSTEIN if CSTEMR fails.
641 *
642  IF( wantz ) THEN
643  order = 'B'
644  ELSE
645  order = 'E'
646  END IF
647 
648  CALL sstebz( 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 cstein( 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 CSTEIN.
661 *
662  indwkn = indwk
663  llwrkn = lwork - indwkn + 1
664  CALL cunmtr( '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 sscal( 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 cswap( 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 CHEEVR
714 *
715  END
LOGICAL function lsame(CA, CB)
LSAME
Definition: lsame.f:54
REAL function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
subroutine cstemr(JOBZ, RANGE, N, D, E, VL, VU, IL, IU, M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK, IWORK, LIWORK, INFO)
CSTEMR
Definition: cstemr.f:328
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:61
subroutine cunmtr(SIDE, UPLO, TRANS, M, N, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMTR
Definition: cunmtr.f:172
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:51
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:52
subroutine cheevr(JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, ISUPPZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
CHEEVR computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices ...
Definition: cheevr.f:347
INTEGER function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:81
subroutine cstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
CSTEIN
Definition: cstein.f:182
subroutine sstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
SSTEBZ
Definition: sstebz.f:262
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
Definition: xerbla-fortran:9
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:53
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:87
REAL function clansy(NORM, UPLO, N, A, LDA, WORK)
CLANSY 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: clansy.f:124
subroutine chetrd(UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO)
CHETRD
Definition: chetrd.f:193
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:54