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