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*> \date June 2016 254* 255*> \ingroup complexHEeigen 256* 257* ===================================================================== 258 SUBROUTINE CHEEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, 259 $ ABSTOL, M, W, Z, LDZ, WORK, LWORK, RWORK, 260 $ IWORK, IFAIL, INFO ) 261* 262* -- LAPACK driver routine (version 3.7.0) -- 263* -- LAPACK is a software package provided by Univ. of Tennessee, -- 264* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 265* June 2016 266* 267* .. Scalar Arguments .. 268 CHARACTER JOBZ, RANGE, UPLO 269 INTEGER IL, INFO, IU, LDA, LDZ, LWORK, M, N 270 REAL ABSTOL, VL, VU 271* .. 272* .. Array Arguments .. 273 INTEGER IFAIL( * ), IWORK( * ) 274 REAL RWORK( * ), W( * ) 275 COMPLEX A( LDA, * ), WORK( * ), Z( LDZ, * ) 276* .. 277* 278* ===================================================================== 279* 280* .. Parameters .. 281 REAL ZERO, ONE 282 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 283 COMPLEX CONE 284 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) 285* .. 286* .. Local Scalars .. 287 LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG, 288 $ WANTZ 289 CHARACTER ORDER 290 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL, 291 $ INDISP, INDIWK, INDRWK, INDTAU, INDWRK, ISCALE, 292 $ ITMP1, J, JJ, LLWORK, LWKMIN, LWKOPT, NB, 293 $ NSPLIT 294 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, 295 $ SIGMA, SMLNUM, TMP1, VLL, VUU 296* .. 297* .. External Functions .. 298 LOGICAL LSAME 299 INTEGER ILAENV 300 REAL SLAMCH, CLANHE 301 EXTERNAL LSAME, ILAENV, SLAMCH, CLANHE 302* .. 303* .. External Subroutines .. 304 EXTERNAL SCOPY, SSCAL, SSTEBZ, SSTERF, XERBLA, CSSCAL, 305 $ CHETRD, CLACPY, CSTEIN, CSTEQR, CSWAP, CUNGTR, 306 $ CUNMTR 307* .. 308* .. Intrinsic Functions .. 309 INTRINSIC REAL, MAX, MIN, SQRT 310* .. 311* .. Executable Statements .. 312* 313* Test the input parameters. 314* 315 LOWER = LSAME( UPLO, 'L' ) 316 WANTZ = LSAME( JOBZ, 'V' ) 317 ALLEIG = LSAME( RANGE, 'A' ) 318 VALEIG = LSAME( RANGE, 'V' ) 319 INDEIG = LSAME( RANGE, 'I' ) 320 LQUERY = ( LWORK.EQ.-1 ) 321* 322 INFO = 0 323 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 324 INFO = -1 325 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN 326 INFO = -2 327 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 328 INFO = -3 329 ELSE IF( N.LT.0 ) THEN 330 INFO = -4 331 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 332 INFO = -6 333 ELSE 334 IF( VALEIG ) THEN 335 IF( N.GT.0 .AND. VU.LE.VL ) 336 $ INFO = -8 337 ELSE IF( INDEIG ) THEN 338 IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN 339 INFO = -9 340 ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN 341 INFO = -10 342 END IF 343 END IF 344 END IF 345 IF( INFO.EQ.0 ) THEN 346 IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 347 INFO = -15 348 END IF 349 END IF 350* 351 IF( INFO.EQ.0 ) THEN 352 IF( N.LE.1 ) THEN 353 LWKMIN = 1 354 WORK( 1 ) = LWKMIN 355 ELSE 356 LWKMIN = 2*N 357 NB = ILAENV( 1, 'CHETRD', UPLO, N, -1, -1, -1 ) 358 NB = MAX( NB, ILAENV( 1, 'CUNMTR', UPLO, N, -1, -1, -1 ) ) 359 LWKOPT = MAX( 1, ( NB + 1 )*N ) 360 WORK( 1 ) = LWKOPT 361 END IF 362* 363 IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) 364 $ INFO = -17 365 END IF 366* 367 IF( INFO.NE.0 ) THEN 368 CALL XERBLA( 'CHEEVX', -INFO ) 369 RETURN 370 ELSE IF( LQUERY ) THEN 371 RETURN 372 END IF 373* 374* Quick return if possible 375* 376 M = 0 377 IF( N.EQ.0 ) THEN 378 RETURN 379 END IF 380* 381 IF( N.EQ.1 ) THEN 382 IF( ALLEIG .OR. INDEIG ) THEN 383 M = 1 384 W( 1 ) = A( 1, 1 ) 385 ELSE IF( VALEIG ) THEN 386 IF( VL.LT.REAL( A( 1, 1 ) ) .AND. VU.GE.REAL( A( 1, 1 ) ) ) 387 $ THEN 388 M = 1 389 W( 1 ) = A( 1, 1 ) 390 END IF 391 END IF 392 IF( WANTZ ) 393 $ Z( 1, 1 ) = CONE 394 RETURN 395 END IF 396* 397* Get machine constants. 398* 399 SAFMIN = SLAMCH( 'Safe minimum' ) 400 EPS = SLAMCH( 'Precision' ) 401 SMLNUM = SAFMIN / EPS 402 BIGNUM = ONE / SMLNUM 403 RMIN = SQRT( SMLNUM ) 404 RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) ) 405* 406* Scale matrix to allowable range, if necessary. 407* 408 ISCALE = 0 409 ABSTLL = ABSTOL 410 IF( VALEIG ) THEN 411 VLL = VL 412 VUU = VU 413 END IF 414 ANRM = CLANHE( 'M', UPLO, N, A, LDA, RWORK ) 415 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 416 ISCALE = 1 417 SIGMA = RMIN / ANRM 418 ELSE IF( ANRM.GT.RMAX ) THEN 419 ISCALE = 1 420 SIGMA = RMAX / ANRM 421 END IF 422 IF( ISCALE.EQ.1 ) THEN 423 IF( LOWER ) THEN 424 DO 10 J = 1, N 425 CALL CSSCAL( N-J+1, SIGMA, A( J, J ), 1 ) 426 10 CONTINUE 427 ELSE 428 DO 20 J = 1, N 429 CALL CSSCAL( J, SIGMA, A( 1, J ), 1 ) 430 20 CONTINUE 431 END IF 432 IF( ABSTOL.GT.0 ) 433 $ ABSTLL = ABSTOL*SIGMA 434 IF( VALEIG ) THEN 435 VLL = VL*SIGMA 436 VUU = VU*SIGMA 437 END IF 438 END IF 439* 440* Call CHETRD to reduce Hermitian matrix to tridiagonal form. 441* 442 INDD = 1 443 INDE = INDD + N 444 INDRWK = INDE + N 445 INDTAU = 1 446 INDWRK = INDTAU + N 447 LLWORK = LWORK - INDWRK + 1 448 CALL CHETRD( UPLO, N, A, LDA, RWORK( INDD ), RWORK( INDE ), 449 $ WORK( INDTAU ), WORK( INDWRK ), LLWORK, IINFO ) 450* 451* If all eigenvalues are desired and ABSTOL is less than or equal to 452* zero, then call SSTERF or CUNGTR and CSTEQR. If this fails for 453* some eigenvalue, then try SSTEBZ. 454* 455 TEST = .FALSE. 456 IF( INDEIG ) THEN 457 IF( IL.EQ.1 .AND. IU.EQ.N ) THEN 458 TEST = .TRUE. 459 END IF 460 END IF 461 IF( ( ALLEIG .OR. TEST ) .AND. ( ABSTOL.LE.ZERO ) ) THEN 462 CALL SCOPY( N, RWORK( INDD ), 1, W, 1 ) 463 INDEE = INDRWK + 2*N 464 IF( .NOT.WANTZ ) THEN 465 CALL SCOPY( N-1, RWORK( INDE ), 1, RWORK( INDEE ), 1 ) 466 CALL SSTERF( N, W, RWORK( INDEE ), INFO ) 467 ELSE 468 CALL CLACPY( 'A', N, N, A, LDA, Z, LDZ ) 469 CALL CUNGTR( UPLO, N, Z, LDZ, WORK( INDTAU ), 470 $ WORK( INDWRK ), LLWORK, IINFO ) 471 CALL SCOPY( N-1, RWORK( INDE ), 1, RWORK( INDEE ), 1 ) 472 CALL CSTEQR( JOBZ, N, W, RWORK( INDEE ), Z, LDZ, 473 $ RWORK( INDRWK ), INFO ) 474 IF( INFO.EQ.0 ) THEN 475 DO 30 I = 1, N 476 IFAIL( I ) = 0 477 30 CONTINUE 478 END IF 479 END IF 480 IF( INFO.EQ.0 ) THEN 481 M = N 482 GO TO 40 483 END IF 484 INFO = 0 485 END IF 486* 487* Otherwise, call SSTEBZ and, if eigenvectors are desired, CSTEIN. 488* 489 IF( WANTZ ) THEN 490 ORDER = 'B' 491 ELSE 492 ORDER = 'E' 493 END IF 494 INDIBL = 1 495 INDISP = INDIBL + N 496 INDIWK = INDISP + N 497 CALL SSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL, 498 $ RWORK( INDD ), RWORK( INDE ), M, NSPLIT, W, 499 $ IWORK( INDIBL ), IWORK( INDISP ), RWORK( INDRWK ), 500 $ IWORK( INDIWK ), INFO ) 501* 502 IF( WANTZ ) THEN 503 CALL CSTEIN( N, RWORK( INDD ), RWORK( INDE ), M, W, 504 $ IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ, 505 $ RWORK( INDRWK ), IWORK( INDIWK ), IFAIL, INFO ) 506* 507* Apply unitary matrix used in reduction to tridiagonal 508* form to eigenvectors returned by CSTEIN. 509* 510 CALL CUNMTR( 'L', UPLO, 'N', N, M, A, LDA, WORK( INDTAU ), Z, 511 $ LDZ, WORK( INDWRK ), LLWORK, IINFO ) 512 END IF 513* 514* If matrix was scaled, then rescale eigenvalues appropriately. 515* 516 40 CONTINUE 517 IF( ISCALE.EQ.1 ) THEN 518 IF( INFO.EQ.0 ) THEN 519 IMAX = M 520 ELSE 521 IMAX = INFO - 1 522 END IF 523 CALL SSCAL( IMAX, ONE / SIGMA, W, 1 ) 524 END IF 525* 526* If eigenvalues are not in order, then sort them, along with 527* eigenvectors. 528* 529 IF( WANTZ ) THEN 530 DO 60 J = 1, M - 1 531 I = 0 532 TMP1 = W( J ) 533 DO 50 JJ = J + 1, M 534 IF( W( JJ ).LT.TMP1 ) THEN 535 I = JJ 536 TMP1 = W( JJ ) 537 END IF 538 50 CONTINUE 539* 540 IF( I.NE.0 ) THEN 541 ITMP1 = IWORK( INDIBL+I-1 ) 542 W( I ) = W( J ) 543 IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 ) 544 W( J ) = TMP1 545 IWORK( INDIBL+J-1 ) = ITMP1 546 CALL CSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 ) 547 IF( INFO.NE.0 ) THEN 548 ITMP1 = IFAIL( I ) 549 IFAIL( I ) = IFAIL( J ) 550 IFAIL( J ) = ITMP1 551 END IF 552 END IF 553 60 CONTINUE 554 END IF 555* 556* Set WORK(1) to optimal complex workspace size. 557* 558 WORK( 1 ) = LWKOPT 559* 560 RETURN 561* 562* End of CHEEVX 563* 564 END 565