LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zheevr_2stage.f
Go to the documentation of this file.
1*> \brief <b> ZHEEVR_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices</b>
2*
3* @precisions fortran z -> s d c
4*
5* =========== DOCUMENTATION ===========
6*
7* Online html documentation available at
8* http://www.netlib.org/lapack/explore-html/
9*
10*> \htmlonly
11*> Download ZHEEVR_2STAGE + dependencies
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zheevr_2stage.f">
13*> [TGZ]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zheevr_2stage.f">
15*> [ZIP]</a>
16*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zheevr_2stage.f">
17*> [TXT]</a>
18*> \endhtmlonly
19*
20* Definition:
21* ===========
22*
23* SUBROUTINE ZHEEVR_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* DOUBLE PRECISION ABSTOL, VL, VU
35* ..
36* .. Array Arguments ..
37* INTEGER ISUPPZ( * ), IWORK( * )
38* DOUBLE PRECISION RWORK( * ), W( * )
39* COMPLEX*16 A( LDA, * ), WORK( * ), Z( LDZ, * )
40* ..
41*
42*
43*> \par Purpose:
44* =============
45*>
46*> \verbatim
47*>
48*> ZHEEVR_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*> ZHEEVR_2STAGE first reduces the matrix A to tridiagonal form T with a call
55*> to ZHETRD. Then, whenever possible, ZHEEVR_2STAGE calls ZSTEMR to compute
56*> eigenspectrum using Relatively Robust Representations. ZSTEMR
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 ZSTEMR'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 : ZHEEVR_2STAGE calls ZSTEMR when the full spectrum is requested
98*> on machines which conform to the ieee-754 floating point standard.
99*> ZHEEVR_2STAGE calls DSTEBZ and ZSTEIN on non-ieee machines and
100*> when partial spectrum requests are made.
101*>
102*> Normal execution of ZSTEMR 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, DSTEBZ and
127*> ZSTEIN 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*16 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 DOUBLE PRECISION
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 DOUBLE PRECISION
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 DOUBLE PRECISION
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*> DLAMCH( '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 DOUBLE PRECISION 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*16 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 ZSTEMR (tridiagonal
268*> matrix). The support of the eigenvectors of A is typically
269*> 1:N because of the unitary transformations applied by ZUNMTR.
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*16 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 DOUBLE PRECISION 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 heevr_2stage
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 zheevr_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 DOUBLE PRECISION ABSTOL, VL, VU
418* ..
419* .. Array Arguments ..
420 INTEGER ISUPPZ( * ), IWORK( * )
421 DOUBLE PRECISION RWORK( * ), W( * )
422 COMPLEX*16 A( LDA, * ), WORK( * ), Z( LDZ, * )
423* ..
424*
425* =====================================================================
426*
427* .. Parameters ..
428 DOUBLE PRECISION ZERO, ONE, TWO
429 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0, two = 2.0d+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 DOUBLE PRECISION 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 DOUBLE PRECISION DLAMCH, ZLANSY
447 EXTERNAL lsame, dlamch, zlansy, ilaenv, ilaenv2stage
448* ..
449* .. External Subroutines ..
450 EXTERNAL dcopy, dscal, dstebz, dsterf, xerbla, zdscal,
452* ..
453* .. Intrinsic Functions ..
454 INTRINSIC dble, max, min, sqrt
455* ..
456* .. Executable Statements ..
457*
458* Test the input parameters.
459*
460 ieeeok = ilaenv( 10, 'ZHEEVR', '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, 'ZHETRD_2STAGE', jobz, n, -1, -1, -1 )
472 ib = ilaenv2stage( 2, 'ZHETRD_2STAGE', jobz, n, kd, -1, -1 )
473 lhtrd = ilaenv2stage( 3, 'ZHETRD_2STAGE', jobz, n, kd, ib, -1 )
474 lwtrd = ilaenv2stage( 4, 'ZHETRD_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( 'ZHEEVR_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 ) = dble( a( 1, 1 ) )
542 ELSE
543 IF( vl.LT.dble( a( 1, 1 ) ) .AND. vu.GE.dble( a( 1, 1 ) ) )
544 $ THEN
545 m = 1
546 w( 1 ) = dble( 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 = dlamch( 'Safe minimum' )
560 eps = dlamch( '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 = zlansy( '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 zdscal( n-j+1, sigma, a( j, j ), 1 )
586 10 CONTINUE
587 ELSE
588 DO 20 j = 1, n
589 CALL zdscal( 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 DSTERF or ZSTEMR fail.
602
603* WORK(INDTAU:INDTAU+N-1) stores the complex scalar factors of the
604* elementary reflectors used in ZHETRD.
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 ZHETRD.
617 indre = indrd + n
618* RWORK(INDRDD:INDRDD+N-1) is a copy of the diagonal entries over
619* -written by ZSTEMR (the DSTERF 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 DSTERF and ZSTEMR.
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 DSTEBZ 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 DSTEBZ 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* ZSTEIN. 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 ZHETRD_2STAGE to reduce Hermitian matrix to tridiagonal form.
645*
646 CALL zhetrd_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 DSTERF or ZSTEMR and ZUNMTR.
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 dcopy( n, rwork( indrd ), 1, w, 1 )
663 CALL dcopy( n-1, rwork( indre ), 1, rwork( indree ), 1 )
664 CALL dsterf( n, w, rwork( indree ), info )
665 ELSE
666 CALL dcopy( n-1, rwork( indre ), 1, rwork( indree ), 1 )
667 CALL dcopy( 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 zstemr( 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 ZSTEMR.
682*
683 IF( wantz .AND. info.EQ.0 ) THEN
684 indwkn = indwk
685 llwrkn = lwork - indwkn + 1
686 CALL zunmtr( '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 DSTEBZ and, if eigenvectors are desired, ZSTEIN.
701* Also call DSTEBZ and ZSTEIN if ZSTEMR fails.
702*
703 IF( wantz ) THEN
704 order = 'B'
705 ELSE
706 order = 'E'
707 END IF
708
709 CALL dstebz( 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 zstein( 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 ZSTEIN.
722*
723 indwkn = indwk
724 llwrkn = lwork - indwkn + 1
725 CALL zunmtr( '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 dscal( 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 zswap( 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 ZHEEVR_2STAGE
775*
776 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine zheevr_2stage(jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w, z, ldz, isuppz, work, lwork, rwork, lrwork, iwork, liwork, info)
ZHEEVR_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE mat...
subroutine zhetrd_2stage(vect, uplo, n, a, lda, d, e, tau, hous2, lhous2, work, lwork, info)
ZHETRD_2STAGE
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine zdscal(n, da, zx, incx)
ZDSCAL
Definition zdscal.f:78
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 zstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
ZSTEIN
Definition zstein.f:182
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:339
subroutine dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:86
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81
subroutine zunmtr(side, uplo, trans, m, n, a, lda, tau, c, ldc, work, lwork, info)
ZUNMTR
Definition zunmtr.f:171