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