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