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