1*> \brief <b> CHEEVX 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 CHEEVX + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cheevx.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cheevx.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cheevx.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE CHEEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, 22* ABSTOL, M, W, Z, LDZ, WORK, LWORK, RWORK, 23* IWORK, IFAIL, INFO ) 24* 25* .. Scalar Arguments .. 26* CHARACTER JOBZ, RANGE, UPLO 27* INTEGER IL, INFO, IU, LDA, LDZ, LWORK, M, N 28* REAL ABSTOL, VL, VU 29* .. 30* .. Array Arguments .. 31* INTEGER IFAIL( * ), IWORK( * ) 32* REAL RWORK( * ), W( * ) 33* COMPLEX A( LDA, * ), WORK( * ), Z( LDZ, * ) 34* .. 35* 36* 37*> \par Purpose: 38* ============= 39*> 40*> \verbatim 41*> 42*> CHEEVX computes selected eigenvalues and, optionally, eigenvectors 43*> of a complex Hermitian matrix A. Eigenvalues and eigenvectors can 44*> be selected by specifying either a range of values or a range of 45*> indices for the desired eigenvalues. 46*> \endverbatim 47* 48* Arguments: 49* ========== 50* 51*> \param[in] JOBZ 52*> \verbatim 53*> JOBZ is CHARACTER*1 54*> = 'N': Compute eigenvalues only; 55*> = 'V': Compute eigenvalues and eigenvectors. 56*> \endverbatim 57*> 58*> \param[in] RANGE 59*> \verbatim 60*> RANGE is CHARACTER*1 61*> = 'A': all eigenvalues will be found. 62*> = 'V': all eigenvalues in the half-open interval (VL,VU] 63*> will be found. 64*> = 'I': the IL-th through IU-th eigenvalues will be found. 65*> \endverbatim 66*> 67*> \param[in] UPLO 68*> \verbatim 69*> UPLO is CHARACTER*1 70*> = 'U': Upper triangle of A is stored; 71*> = 'L': Lower triangle of A is stored. 72*> \endverbatim 73*> 74*> \param[in] N 75*> \verbatim 76*> N is INTEGER 77*> The order of the matrix A. N >= 0. 78*> \endverbatim 79*> 80*> \param[in,out] A 81*> \verbatim 82*> A is COMPLEX array, dimension (LDA, N) 83*> On entry, the Hermitian matrix A. If UPLO = 'U', the 84*> leading N-by-N upper triangular part of A contains the 85*> upper triangular part of the matrix A. If UPLO = 'L', 86*> the leading N-by-N lower triangular part of A contains 87*> the lower triangular part of the matrix A. 88*> On exit, the lower triangle (if UPLO='L') or the upper 89*> triangle (if UPLO='U') of A, including the diagonal, is 90*> destroyed. 91*> \endverbatim 92*> 93*> \param[in] LDA 94*> \verbatim 95*> LDA is INTEGER 96*> The leading dimension of the array A. LDA >= max(1,N). 97*> \endverbatim 98*> 99*> \param[in] VL 100*> \verbatim 101*> VL is REAL 102*> If RANGE='V', the lower bound of the interval to 103*> be searched for eigenvalues. VL < VU. 104*> Not referenced if RANGE = 'A' or 'I'. 105*> \endverbatim 106*> 107*> \param[in] VU 108*> \verbatim 109*> VU is REAL 110*> If RANGE='V', the upper bound of the interval to 111*> be searched for eigenvalues. VL < VU. 112*> Not referenced if RANGE = 'A' or 'I'. 113*> \endverbatim 114*> 115*> \param[in] IL 116*> \verbatim 117*> IL is INTEGER 118*> If RANGE='I', the index of the 119*> smallest eigenvalue to be returned. 120*> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. 121*> Not referenced if RANGE = 'A' or 'V'. 122*> \endverbatim 123*> 124*> \param[in] IU 125*> \verbatim 126*> IU is INTEGER 127*> If RANGE='I', the index of the 128*> largest eigenvalue to be returned. 129*> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. 130*> Not referenced if RANGE = 'A' or 'V'. 131*> \endverbatim 132*> 133*> \param[in] ABSTOL 134*> \verbatim 135*> ABSTOL is REAL 136*> The absolute error tolerance for the eigenvalues. 137*> An approximate eigenvalue is accepted as converged 138*> when it is determined to lie in an interval [a,b] 139*> of width less than or equal to 140*> 141*> ABSTOL + EPS * max( |a|,|b| ) , 142*> 143*> where EPS is the machine precision. If ABSTOL is less than 144*> or equal to zero, then EPS*|T| will be used in its place, 145*> where |T| is the 1-norm of the tridiagonal matrix obtained 146*> by reducing A to tridiagonal form. 147*> 148*> Eigenvalues will be computed most accurately when ABSTOL is 149*> set to twice the underflow threshold 2*SLAMCH('S'), not zero. 150*> If this routine returns with INFO>0, indicating that some 151*> eigenvectors did not converge, try setting ABSTOL to 152*> 2*SLAMCH('S'). 153*> 154*> See "Computing Small Singular Values of Bidiagonal Matrices 155*> with Guaranteed High Relative Accuracy," by Demmel and 156*> Kahan, LAPACK Working Note #3. 157*> \endverbatim 158*> 159*> \param[out] M 160*> \verbatim 161*> M is INTEGER 162*> The total number of eigenvalues found. 0 <= M <= N. 163*> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. 164*> \endverbatim 165*> 166*> \param[out] W 167*> \verbatim 168*> W is REAL array, dimension (N) 169*> On normal exit, the first M elements contain the selected 170*> eigenvalues in ascending order. 171*> \endverbatim 172*> 173*> \param[out] Z 174*> \verbatim 175*> Z is COMPLEX array, dimension (LDZ, max(1,M)) 176*> If JOBZ = 'V', then if INFO = 0, the first M columns of Z 177*> contain the orthonormal eigenvectors of the matrix A 178*> corresponding to the selected eigenvalues, with the i-th 179*> column of Z holding the eigenvector associated with W(i). 180*> If an eigenvector fails to converge, then that column of Z 181*> contains the latest approximation to the eigenvector, and the 182*> index of the eigenvector is returned in IFAIL. 183*> If JOBZ = 'N', then Z is not referenced. 184*> Note: the user must ensure that at least max(1,M) columns are 185*> supplied in the array Z; if RANGE = 'V', the exact value of M 186*> is not known in advance and an upper bound must be used. 187*> \endverbatim 188*> 189*> \param[in] LDZ 190*> \verbatim 191*> LDZ is INTEGER 192*> The leading dimension of the array Z. LDZ >= 1, and if 193*> JOBZ = 'V', LDZ >= max(1,N). 194*> \endverbatim 195*> 196*> \param[out] WORK 197*> \verbatim 198*> WORK is COMPLEX array, dimension (MAX(1,LWORK)) 199*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 200*> \endverbatim 201*> 202*> \param[in] LWORK 203*> \verbatim 204*> LWORK is INTEGER 205*> The length of the array WORK. LWORK >= 1, when N <= 1; 206*> otherwise 2*N. 207*> For optimal efficiency, LWORK >= (NB+1)*N, 208*> where NB is the max of the blocksize for CHETRD and for 209*> CUNMTR as returned by ILAENV. 210*> 211*> If LWORK = -1, then a workspace query is assumed; the routine 212*> only calculates the optimal size of the WORK array, returns 213*> this value as the first entry of the WORK array, and no error 214*> message related to LWORK is issued by XERBLA. 215*> \endverbatim 216*> 217*> \param[out] RWORK 218*> \verbatim 219*> RWORK is REAL array, dimension (7*N) 220*> \endverbatim 221*> 222*> \param[out] IWORK 223*> \verbatim 224*> IWORK is INTEGER array, dimension (5*N) 225*> \endverbatim 226*> 227*> \param[out] IFAIL 228*> \verbatim 229*> IFAIL is INTEGER array, dimension (N) 230*> If JOBZ = 'V', then if INFO = 0, the first M elements of 231*> IFAIL are zero. If INFO > 0, then IFAIL contains the 232*> indices of the eigenvectors that failed to converge. 233*> If JOBZ = 'N', then IFAIL is not referenced. 234*> \endverbatim 235*> 236*> \param[out] INFO 237*> \verbatim 238*> INFO is INTEGER 239*> = 0: successful exit 240*> < 0: if INFO = -i, the i-th argument had an illegal value 241*> > 0: if INFO = i, then i eigenvectors failed to converge. 242*> Their indices are stored in array IFAIL. 243*> \endverbatim 244* 245* Authors: 246* ======== 247* 248*> \author Univ. of Tennessee 249*> \author Univ. of California Berkeley 250*> \author Univ. of Colorado Denver 251*> \author NAG Ltd. 252* 253*> \ingroup complexHEeigen 254* 255* ===================================================================== 256 SUBROUTINE CHEEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, 257 $ ABSTOL, M, W, Z, LDZ, WORK, LWORK, RWORK, 258 $ IWORK, IFAIL, INFO ) 259* 260* -- LAPACK driver routine -- 261* -- LAPACK is a software package provided by Univ. of Tennessee, -- 262* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 263* 264* .. Scalar Arguments .. 265 CHARACTER JOBZ, RANGE, UPLO 266 INTEGER IL, INFO, IU, LDA, LDZ, LWORK, M, N 267 REAL ABSTOL, VL, VU 268* .. 269* .. Array Arguments .. 270 INTEGER IFAIL( * ), IWORK( * ) 271 REAL RWORK( * ), W( * ) 272 COMPLEX A( LDA, * ), WORK( * ), Z( LDZ, * ) 273* .. 274* 275* ===================================================================== 276* 277* .. Parameters .. 278 REAL ZERO, ONE 279 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 280 COMPLEX CONE 281 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) 282* .. 283* .. Local Scalars .. 284 LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG, 285 $ WANTZ 286 CHARACTER ORDER 287 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL, 288 $ INDISP, INDIWK, INDRWK, INDTAU, INDWRK, ISCALE, 289 $ ITMP1, J, JJ, LLWORK, LWKMIN, LWKOPT, NB, 290 $ NSPLIT 291 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, 292 $ SIGMA, SMLNUM, TMP1, VLL, VUU 293* .. 294* .. External Functions .. 295 LOGICAL LSAME 296 INTEGER ILAENV 297 REAL SLAMCH, CLANHE 298 EXTERNAL LSAME, ILAENV, SLAMCH, CLANHE 299* .. 300* .. External Subroutines .. 301 EXTERNAL SCOPY, SSCAL, SSTEBZ, SSTERF, XERBLA, CSSCAL, 302 $ CHETRD, CLACPY, CSTEIN, CSTEQR, CSWAP, CUNGTR, 303 $ CUNMTR 304* .. 305* .. Intrinsic Functions .. 306 INTRINSIC REAL, MAX, MIN, SQRT 307* .. 308* .. Executable Statements .. 309* 310* Test the input parameters. 311* 312 LOWER = LSAME( UPLO, 'L' ) 313 WANTZ = LSAME( JOBZ, 'V' ) 314 ALLEIG = LSAME( RANGE, 'A' ) 315 VALEIG = LSAME( RANGE, 'V' ) 316 INDEIG = LSAME( RANGE, 'I' ) 317 LQUERY = ( LWORK.EQ.-1 ) 318* 319 INFO = 0 320 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 321 INFO = -1 322 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN 323 INFO = -2 324 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 325 INFO = -3 326 ELSE IF( N.LT.0 ) THEN 327 INFO = -4 328 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 329 INFO = -6 330 ELSE 331 IF( VALEIG ) THEN 332 IF( N.GT.0 .AND. VU.LE.VL ) 333 $ INFO = -8 334 ELSE IF( INDEIG ) THEN 335 IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN 336 INFO = -9 337 ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN 338 INFO = -10 339 END IF 340 END IF 341 END IF 342 IF( INFO.EQ.0 ) THEN 343 IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 344 INFO = -15 345 END IF 346 END IF 347* 348 IF( INFO.EQ.0 ) THEN 349 IF( N.LE.1 ) THEN 350 LWKMIN = 1 351 WORK( 1 ) = LWKMIN 352 ELSE 353 LWKMIN = 2*N 354 NB = ILAENV( 1, 'CHETRD', UPLO, N, -1, -1, -1 ) 355 NB = MAX( NB, ILAENV( 1, 'CUNMTR', UPLO, N, -1, -1, -1 ) ) 356 LWKOPT = MAX( 1, ( NB + 1 )*N ) 357 WORK( 1 ) = LWKOPT 358 END IF 359* 360 IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) 361 $ INFO = -17 362 END IF 363* 364 IF( INFO.NE.0 ) THEN 365 CALL XERBLA( 'CHEEVX', -INFO ) 366 RETURN 367 ELSE IF( LQUERY ) THEN 368 RETURN 369 END IF 370* 371* Quick return if possible 372* 373 M = 0 374 IF( N.EQ.0 ) THEN 375 RETURN 376 END IF 377* 378 IF( N.EQ.1 ) THEN 379 IF( ALLEIG .OR. INDEIG ) THEN 380 M = 1 381 W( 1 ) = REAL( A( 1, 1 ) ) 382 ELSE IF( VALEIG ) THEN 383 IF( VL.LT.REAL( A( 1, 1 ) ) .AND. VU.GE.REAL( A( 1, 1 ) ) ) 384 $ THEN 385 M = 1 386 W( 1 ) = REAL( A( 1, 1 ) ) 387 END IF 388 END IF 389 IF( WANTZ ) 390 $ Z( 1, 1 ) = CONE 391 RETURN 392 END IF 393* 394* Get machine constants. 395* 396 SAFMIN = SLAMCH( 'Safe minimum' ) 397 EPS = SLAMCH( 'Precision' ) 398 SMLNUM = SAFMIN / EPS 399 BIGNUM = ONE / SMLNUM 400 RMIN = SQRT( SMLNUM ) 401 RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) ) 402* 403* Scale matrix to allowable range, if necessary. 404* 405 ISCALE = 0 406 ABSTLL = ABSTOL 407 IF( VALEIG ) THEN 408 VLL = VL 409 VUU = VU 410 END IF 411 ANRM = CLANHE( 'M', UPLO, N, A, LDA, RWORK ) 412 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 413 ISCALE = 1 414 SIGMA = RMIN / ANRM 415 ELSE IF( ANRM.GT.RMAX ) THEN 416 ISCALE = 1 417 SIGMA = RMAX / ANRM 418 END IF 419 IF( ISCALE.EQ.1 ) THEN 420 IF( LOWER ) THEN 421 DO 10 J = 1, N 422 CALL CSSCAL( N-J+1, SIGMA, A( J, J ), 1 ) 423 10 CONTINUE 424 ELSE 425 DO 20 J = 1, N 426 CALL CSSCAL( J, SIGMA, A( 1, J ), 1 ) 427 20 CONTINUE 428 END IF 429 IF( ABSTOL.GT.0 ) 430 $ ABSTLL = ABSTOL*SIGMA 431 IF( VALEIG ) THEN 432 VLL = VL*SIGMA 433 VUU = VU*SIGMA 434 END IF 435 END IF 436* 437* Call CHETRD to reduce Hermitian matrix to tridiagonal form. 438* 439 INDD = 1 440 INDE = INDD + N 441 INDRWK = INDE + N 442 INDTAU = 1 443 INDWRK = INDTAU + N 444 LLWORK = LWORK - INDWRK + 1 445 CALL CHETRD( UPLO, N, A, LDA, RWORK( INDD ), RWORK( INDE ), 446 $ WORK( INDTAU ), WORK( INDWRK ), LLWORK, IINFO ) 447* 448* If all eigenvalues are desired and ABSTOL is less than or equal to 449* zero, then call SSTERF or CUNGTR and CSTEQR. If this fails for 450* some eigenvalue, then try SSTEBZ. 451* 452 TEST = .FALSE. 453 IF( INDEIG ) THEN 454 IF( IL.EQ.1 .AND. IU.EQ.N ) THEN 455 TEST = .TRUE. 456 END IF 457 END IF 458 IF( ( ALLEIG .OR. TEST ) .AND. ( ABSTOL.LE.ZERO ) ) THEN 459 CALL SCOPY( N, RWORK( INDD ), 1, W, 1 ) 460 INDEE = INDRWK + 2*N 461 IF( .NOT.WANTZ ) THEN 462 CALL SCOPY( N-1, RWORK( INDE ), 1, RWORK( INDEE ), 1 ) 463 CALL SSTERF( N, W, RWORK( INDEE ), INFO ) 464 ELSE 465 CALL CLACPY( 'A', N, N, A, LDA, Z, LDZ ) 466 CALL CUNGTR( UPLO, N, Z, LDZ, WORK( INDTAU ), 467 $ WORK( INDWRK ), LLWORK, IINFO ) 468 CALL SCOPY( N-1, RWORK( INDE ), 1, RWORK( INDEE ), 1 ) 469 CALL CSTEQR( JOBZ, N, W, RWORK( INDEE ), Z, LDZ, 470 $ RWORK( INDRWK ), INFO ) 471 IF( INFO.EQ.0 ) THEN 472 DO 30 I = 1, N 473 IFAIL( I ) = 0 474 30 CONTINUE 475 END IF 476 END IF 477 IF( INFO.EQ.0 ) THEN 478 M = N 479 GO TO 40 480 END IF 481 INFO = 0 482 END IF 483* 484* Otherwise, call SSTEBZ and, if eigenvectors are desired, CSTEIN. 485* 486 IF( WANTZ ) THEN 487 ORDER = 'B' 488 ELSE 489 ORDER = 'E' 490 END IF 491 INDIBL = 1 492 INDISP = INDIBL + N 493 INDIWK = INDISP + N 494 CALL SSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL, 495 $ RWORK( INDD ), RWORK( INDE ), M, NSPLIT, W, 496 $ IWORK( INDIBL ), IWORK( INDISP ), RWORK( INDRWK ), 497 $ IWORK( INDIWK ), INFO ) 498* 499 IF( WANTZ ) THEN 500 CALL CSTEIN( N, RWORK( INDD ), RWORK( INDE ), M, W, 501 $ IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ, 502 $ RWORK( INDRWK ), IWORK( INDIWK ), IFAIL, INFO ) 503* 504* Apply unitary matrix used in reduction to tridiagonal 505* form to eigenvectors returned by CSTEIN. 506* 507 CALL CUNMTR( 'L', UPLO, 'N', N, M, A, LDA, WORK( INDTAU ), Z, 508 $ LDZ, WORK( INDWRK ), LLWORK, IINFO ) 509 END IF 510* 511* If matrix was scaled, then rescale eigenvalues appropriately. 512* 513 40 CONTINUE 514 IF( ISCALE.EQ.1 ) THEN 515 IF( INFO.EQ.0 ) THEN 516 IMAX = M 517 ELSE 518 IMAX = INFO - 1 519 END IF 520 CALL SSCAL( IMAX, ONE / SIGMA, W, 1 ) 521 END IF 522* 523* If eigenvalues are not in order, then sort them, along with 524* eigenvectors. 525* 526 IF( WANTZ ) THEN 527 DO 60 J = 1, M - 1 528 I = 0 529 TMP1 = W( J ) 530 DO 50 JJ = J + 1, M 531 IF( W( JJ ).LT.TMP1 ) THEN 532 I = JJ 533 TMP1 = W( JJ ) 534 END IF 535 50 CONTINUE 536* 537 IF( I.NE.0 ) THEN 538 ITMP1 = IWORK( INDIBL+I-1 ) 539 W( I ) = W( J ) 540 IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 ) 541 W( J ) = TMP1 542 IWORK( INDIBL+J-1 ) = ITMP1 543 CALL CSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 ) 544 IF( INFO.NE.0 ) THEN 545 ITMP1 = IFAIL( I ) 546 IFAIL( I ) = IFAIL( J ) 547 IFAIL( J ) = ITMP1 548 END IF 549 END IF 550 60 CONTINUE 551 END IF 552* 553* Set WORK(1) to optimal complex workspace size. 554* 555 WORK( 1 ) = LWKOPT 556* 557 RETURN 558* 559* End of CHEEVX 560* 561 END 562