1 SUBROUTINE DSTEGR2( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, 2 $ M, W, Z, LDZ, NZC, ISUPPZ, WORK, LWORK, IWORK, 3 $ LIWORK, DOL, DOU, ZOFFSET, INFO ) 4* 5* -- ScaLAPACK auxiliary routine (version 2.0) -- 6* Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver 7* July 4, 2010 8* 9* .. Scalar Arguments .. 10 CHARACTER JOBZ, RANGE 11 INTEGER DOL, DOU, IL, INFO, IU, 12 $ LDZ, NZC, LIWORK, LWORK, M, N, ZOFFSET 13 DOUBLE PRECISION VL, VU 14 15* .. 16* .. Array Arguments .. 17 INTEGER ISUPPZ( * ), IWORK( * ) 18 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ) 19 DOUBLE PRECISION Z( LDZ, * ) 20* .. 21* 22* Purpose 23* ======= 24* 25* DSTEGR2 computes selected eigenvalues and, optionally, eigenvectors 26* of a real symmetric tridiagonal matrix T. It is invoked in the 27* ScaLAPACK MRRR driver PDSYEVR and the corresponding Hermitian 28* version either when only eigenvalues are to be computed, or when only 29* a single processor is used (the sequential-like case). 30* 31* DSTEGR2 has been adapted from LAPACK's DSTEGR. Please note the 32* following crucial changes. 33* 34* 1. The calling sequence has two additional INTEGER parameters, 35* DOL and DOU, that should satisfy M>=DOU>=DOL>=1. 36* DSTEGR2 ONLY computes the eigenpairs 37* corresponding to eigenvalues DOL through DOU in W. (That is, 38* instead of computing the eigenpairs belonging to W(1) 39* through W(M), only the eigenvectors belonging to eigenvalues 40* W(DOL) through W(DOU) are computed. In this case, only the 41* eigenvalues DOL:DOU are guaranteed to be fully accurate. 42* 43* 2. M is NOT the number of eigenvalues specified by RANGE, but is 44* M = DOU - DOL + 1. This concerns the case where only eigenvalues 45* are computed, but on more than one processor. Thus, in this case 46* M refers to the number of eigenvalues computed on this processor. 47* 48* 3. The arrays W and Z might not contain all the wanted eigenpairs 49* locally, instead this information is distributed over other 50* processors. 51* 52* Arguments 53* ========= 54* 55* JOBZ (input) CHARACTER*1 56* = 'N': Compute eigenvalues only; 57* = 'V': Compute eigenvalues and eigenvectors. 58* 59* RANGE (input) CHARACTER*1 60* = 'A': all eigenvalues will be found. 61* = 'V': all eigenvalues in the half-open interval (VL,VU] 62* will be found. 63* = 'I': the IL-th through IU-th eigenvalues will be found. 64* 65* N (input) INTEGER 66* The order of the matrix. N >= 0. 67* 68* D (input/output) DOUBLE PRECISION array, dimension (N) 69* On entry, the N diagonal elements of the tridiagonal matrix 70* T. On exit, D is overwritten. 71* 72* E (input/output) DOUBLE PRECISION array, dimension (N) 73* On entry, the (N-1) subdiagonal elements of the tridiagonal 74* matrix T in elements 1 to N-1 of E. E(N) need not be set on 75* input, but is used internally as workspace. 76* On exit, E is overwritten. 77* 78* VL (input) DOUBLE PRECISION 79* VU (input) DOUBLE PRECISION 80* If RANGE='V', the lower and upper bounds of the interval to 81* be searched for eigenvalues. VL < VU. 82* Not referenced if RANGE = 'A' or 'I'. 83* 84* IL (input) INTEGER 85* IU (input) INTEGER 86* If RANGE='I', the indices (in ascending order) of the 87* smallest and largest eigenvalues to be returned. 88* 1 <= IL <= IU <= N, if N > 0. 89* Not referenced if RANGE = 'A' or 'V'. 90* 91* M (output) INTEGER 92* Globally summed over all processors, M equals 93* the total number of eigenvalues found. 0 <= M <= N. 94* If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. 95* The local output equals M = DOU - DOL + 1. 96* 97* W (output) DOUBLE PRECISION array, dimension (N) 98* The first M elements contain the selected eigenvalues in 99* ascending order. Note that immediately after exiting this 100* routine, only the eigenvalues from 101* position DOL:DOU are to reliable on this processor 102* because the eigenvalue computation is done in parallel. 103* Other processors will hold reliable information on other 104* parts of the W array. This information is communicated in 105* the ScaLAPACK driver. 106* 107* Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) ) 108* If JOBZ = 'V', and if INFO = 0, then the first M columns of Z 109* contain some of the orthonormal eigenvectors of the matrix T 110* corresponding to the selected eigenvalues, with the i-th 111* column of Z holding the eigenvector associated with W(i). 112* If JOBZ = 'N', then Z is not referenced. 113* Note: the user must ensure that at least max(1,M) columns are 114* supplied in the array Z; if RANGE = 'V', the exact value of M 115* is not known in advance and can be computed with a workspace 116* query by setting NZC = -1, see below. 117* 118* LDZ (input) INTEGER 119* The leading dimension of the array Z. LDZ >= 1, and if 120* JOBZ = 'V', then LDZ >= max(1,N). 121* 122* NZC (input) INTEGER 123* The number of eigenvectors to be held in the array Z. 124* If RANGE = 'A', then NZC >= max(1,N). 125* If RANGE = 'V', then NZC >= the number of eigenvalues in (VL,VU]. 126* If RANGE = 'I', then NZC >= IU-IL+1. 127* If NZC = -1, then a workspace query is assumed; the 128* routine calculates the number of columns of the array Z that 129* are needed to hold the eigenvectors. 130* This value is returned as the first entry of the Z array, and 131* no error message related to NZC is issued. 132* 133* ISUPPZ (output) INTEGER ARRAY, dimension ( 2*max(1,M) ) 134* The support of the eigenvectors in Z, i.e., the indices 135* indicating the nonzero elements in Z. The i-th computed eigenvector 136* is nonzero only in elements ISUPPZ( 2*i-1 ) through 137* ISUPPZ( 2*i ). This is relevant in the case when the matrix 138* is split. ISUPPZ is only set if N>2. 139* 140* WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) 141* On exit, if INFO = 0, WORK(1) returns the optimal 142* (and minimal) LWORK. 143* 144* LWORK (input) INTEGER 145* The dimension of the array WORK. LWORK >= max(1,18*N) 146* if JOBZ = 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'. 147* If LWORK = -1, then a workspace query is assumed; the routine 148* only calculates the optimal size of the WORK array, returns 149* this value as the first entry of the WORK array, and no error 150* message related to LWORK is issued. 151* 152* IWORK (workspace/output) INTEGER array, dimension (LIWORK) 153* On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 154* 155* LIWORK (input) INTEGER 156* The dimension of the array IWORK. LIWORK >= max(1,10*N) 157* if the eigenvectors are desired, and LIWORK >= max(1,8*N) 158* if only the eigenvalues are to be computed. 159* If LIWORK = -1, then a workspace query is assumed; the 160* routine only calculates the optimal size of the IWORK array, 161* returns this value as the first entry of the IWORK array, and 162* no error message related to LIWORK is issued. 163* 164* DOL (input) INTEGER 165* DOU (input) INTEGER 166* From the eigenvalues W(1:M), only eigenvectors 167* Z(:,DOL) to Z(:,DOU) are computed. 168* If DOL > 1, then Z(:,DOL-1-ZOFFSET) is used and overwritten. 169* If DOU < M, then Z(:,DOU+1-ZOFFSET) is used and overwritten. 170* 171* ZOFFSET (input) INTEGER 172* Offset for storing the eigenpairs when Z is distributed 173* in 1D-cyclic fashion 174* 175* INFO (output) INTEGER 176* On exit, INFO 177* = 0: successful exit 178* other:if INFO = -i, the i-th argument had an illegal value 179* if INFO = 10X, internal error in DLARRE2, 180* if INFO = 20X, internal error in DLARRV. 181* Here, the digit X = ABS( IINFO ) < 10, where IINFO is 182* the nonzero error code returned by DLARRE2 or 183* DLARRV, respectively. 184* 185* ===================================================================== 186* 187* .. Parameters .. 188 DOUBLE PRECISION ZERO, ONE, FOUR, MINRGP 189 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, 190 $ FOUR = 4.0D0, 191 $ MINRGP = 1.0D-3 ) 192* .. 193* .. Local Scalars .. 194 LOGICAL ALLEIG, INDEIG, LQUERY, VALEIG, WANTZ, ZQUERY 195 INTEGER I, IIL, IINDBL, IINDW, IINDWK, IINFO, IINSPL, 196 $ IIU, INDE2, INDERR, INDGP, INDGRS, INDWRK, 197 $ ITMP, ITMP2, J, JJ, LIWMIN, LWMIN, NSPLIT, 198 $ NZCMIN 199 DOUBLE PRECISION BIGNUM, EPS, PIVMIN, RMAX, RMIN, RTOL1, RTOL2, 200 $ SAFMIN, SCALE, SMLNUM, THRESH, TMP, TNRM, WL, 201 $ WU 202* .. 203* .. External Functions .. 204 LOGICAL LSAME 205 DOUBLE PRECISION DLAMCH, DLANST 206 EXTERNAL LSAME, DLAMCH, DLANST 207* .. 208* .. External Subroutines .. 209 EXTERNAL DCOPY, DLAE2, DLAEV2, DLARRC, DLARRE2, 210 $ DLARRV, DLASRT, DSCAL, DSWAP 211* .. 212* .. Intrinsic Functions .. 213 INTRINSIC DBLE, MAX, MIN, SQRT 214* .. 215* .. Executable Statements .. 216* 217* Test the input parameters. 218* 219 WANTZ = LSAME( JOBZ, 'V' ) 220 ALLEIG = LSAME( RANGE, 'A' ) 221 VALEIG = LSAME( RANGE, 'V' ) 222 INDEIG = LSAME( RANGE, 'I' ) 223* 224 LQUERY = ( ( LWORK.EQ.-1 ).OR.( LIWORK.EQ.-1 ) ) 225 ZQUERY = ( NZC.EQ.-1 ) 226 227* DSTEGR2 needs WORK of size 6*N, IWORK of size 3*N. 228* In addition, DLARRE2 needs WORK of size 6*N, IWORK of size 5*N. 229* Furthermore, DLARRV needs WORK of size 12*N, IWORK of size 7*N. 230 IF( WANTZ ) THEN 231 LWMIN = 18*N 232 LIWMIN = 10*N 233 ELSE 234* need less workspace if only the eigenvalues are wanted 235 LWMIN = 12*N 236 LIWMIN = 8*N 237 ENDIF 238 239 WL = ZERO 240 WU = ZERO 241 IIL = 0 242 IIU = 0 243 244 IF( VALEIG ) THEN 245* We do not reference VL, VU in the cases RANGE = 'I','A' 246* The interval (WL, WU] contains all the wanted eigenvalues. 247* It is either given by the user or computed in DLARRE2. 248 WL = VL 249 WU = VU 250 ELSEIF( INDEIG ) THEN 251* We do not reference IL, IU in the cases RANGE = 'V','A' 252 IIL = IL 253 IIU = IU 254 ENDIF 255* 256 INFO = 0 257 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 258 INFO = -1 259 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN 260 INFO = -2 261 ELSE IF( N.LT.0 ) THEN 262 INFO = -3 263 ELSE IF( VALEIG .AND. N.GT.0 .AND. WU.LE.WL ) THEN 264 INFO = -7 265 ELSE IF( INDEIG .AND. ( IIL.LT.1 .OR. IIL.GT.N ) ) THEN 266 INFO = -8 267 ELSE IF( INDEIG .AND. ( IIU.LT.IIL .OR. IIU.GT.N ) ) THEN 268 INFO = -9 269 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 270 INFO = -13 271 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 272 INFO = -17 273 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 274 INFO = -19 275 END IF 276* 277* Get machine constants. 278* 279 SAFMIN = DLAMCH( 'Safe minimum' ) 280 EPS = DLAMCH( 'Precision' ) 281 SMLNUM = SAFMIN / EPS 282 BIGNUM = ONE / SMLNUM 283 RMIN = SQRT( SMLNUM ) 284 RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) ) 285* 286 IF( INFO.EQ.0 ) THEN 287 WORK( 1 ) = LWMIN 288 IWORK( 1 ) = LIWMIN 289* 290 IF( WANTZ .AND. ALLEIG ) THEN 291 NZCMIN = N 292 IIL = 1 293 IIU = N 294 ELSE IF( WANTZ .AND. VALEIG ) THEN 295 CALL DLARRC( 'T', N, VL, VU, D, E, SAFMIN, 296 $ NZCMIN, ITMP, ITMP2, INFO ) 297 IIL = ITMP+1 298 IIU = ITMP2 299 ELSE IF( WANTZ .AND. INDEIG ) THEN 300 NZCMIN = IIU-IIL+1 301 ELSE 302* WANTZ .EQ. FALSE. 303 NZCMIN = 0 304 ENDIF 305 IF( ZQUERY .AND. INFO.EQ.0 ) THEN 306 Z( 1,1 ) = NZCMIN 307 ELSE IF( NZC.LT.NZCMIN .AND. .NOT.ZQUERY ) THEN 308 INFO = -14 309 END IF 310 END IF 311 312 IF ( WANTZ ) THEN 313 IF ( DOL.LT.1 .OR. DOL.GT.NZCMIN ) THEN 314 INFO = -20 315 ENDIF 316 IF ( DOU.LT.1 .OR. DOU.GT.NZCMIN .OR. DOU.LT.DOL) THEN 317 INFO = -21 318 ENDIF 319 ENDIF 320 321 IF( INFO.NE.0 ) THEN 322* 323C Disable sequential error handler 324C for parallel case 325C CALL XERBLA( 'DSTEGR2', -INFO ) 326* 327 RETURN 328 ELSE IF( LQUERY .OR. ZQUERY ) THEN 329 RETURN 330 END IF 331* 332* Quick return if possible 333* 334 M = 0 335 IF( N.EQ.0 ) 336 $ RETURN 337* 338 IF( N.EQ.1 ) THEN 339 IF( ALLEIG .OR. INDEIG ) THEN 340 M = 1 341 W( 1 ) = D( 1 ) 342 ELSE 343 IF( WL.LT.D( 1 ) .AND. WU.GE.D( 1 ) ) THEN 344 M = 1 345 W( 1 ) = D( 1 ) 346 END IF 347 END IF 348 IF( WANTZ ) 349 $ Z( 1, 1 ) = ONE 350 RETURN 351 END IF 352* 353 INDGRS = 1 354 INDERR = 2*N + 1 355 INDGP = 3*N + 1 356 INDE2 = 5*N + 1 357 INDWRK = 6*N + 1 358* 359 IINSPL = 1 360 IINDBL = N + 1 361 IINDW = 2*N + 1 362 IINDWK = 3*N + 1 363* 364* Scale matrix to allowable range, if necessary. 365* 366 SCALE = ONE 367 TNRM = DLANST( 'M', N, D, E ) 368 IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN 369 SCALE = RMIN / TNRM 370 ELSE IF( TNRM.GT.RMAX ) THEN 371 SCALE = RMAX / TNRM 372 END IF 373 IF( SCALE.NE.ONE ) THEN 374 CALL DSCAL( N, SCALE, D, 1 ) 375 CALL DSCAL( N-1, SCALE, E, 1 ) 376 TNRM = TNRM*SCALE 377 IF( VALEIG ) THEN 378* If eigenvalues in interval have to be found, 379* scale (WL, WU] accordingly 380 WL = WL*SCALE 381 WU = WU*SCALE 382 ENDIF 383 END IF 384* 385* Compute the desired eigenvalues of the tridiagonal after splitting 386* into smaller subblocks if the corresponding off-diagonal elements 387* are small 388* THRESH is the splitting parameter for DLARRE2 389* A negative THRESH forces the old splitting criterion based on the 390* size of the off-diagonal. A positive THRESH switches to splitting 391* which preserves relative accuracy. 392* 393 IINFO = -1 394* Set the splitting criterion 395 IF (IINFO.EQ.0) THEN 396 THRESH = EPS 397 ELSE 398 THRESH = -EPS 399 ENDIF 400* 401* Store the squares of the offdiagonal values of T 402 DO 5 J = 1, N-1 403 WORK( INDE2+J-1 ) = E(J)**2 404 5 CONTINUE 405 406* Set the tolerance parameters for bisection 407 IF( .NOT.WANTZ ) THEN 408* DLARRE2 computes the eigenvalues to full precision. 409 RTOL1 = FOUR * EPS 410 RTOL2 = FOUR * EPS 411 ELSE 412* DLARRE2 computes the eigenvalues to less than full precision. 413* DLARRV will refine the eigenvalue approximations, and we can 414* need less accurate initial bisection in DLARRE2. 415* Note: these settings do only affect the subset case and DLARRE2 416 RTOL1 = SQRT(EPS) 417 RTOL2 = MAX( SQRT(EPS)*5.0D-3, FOUR * EPS ) 418 ENDIF 419 CALL DLARRE2( RANGE, N, WL, WU, IIL, IIU, D, E, 420 $ WORK(INDE2), RTOL1, RTOL2, THRESH, NSPLIT, 421 $ IWORK( IINSPL ), M, DOL, DOU, 422 $ W, WORK( INDERR ), 423 $ WORK( INDGP ), IWORK( IINDBL ), 424 $ IWORK( IINDW ), WORK( INDGRS ), PIVMIN, 425 $ WORK( INDWRK ), IWORK( IINDWK ), IINFO ) 426 IF( IINFO.NE.0 ) THEN 427 INFO = 100 + ABS( IINFO ) 428 RETURN 429 END IF 430* Note that if RANGE .NE. 'V', DLARRE2 computes bounds on the desired 431* part of the spectrum. All desired eigenvalues are contained in 432* (WL,WU] 433 434 435 IF( WANTZ ) THEN 436* 437* Compute the desired eigenvectors corresponding to the computed 438* eigenvalues 439* 440 CALL DLARRV( N, WL, WU, D, E, 441 $ PIVMIN, IWORK( IINSPL ), M, 442 $ DOL, DOU, MINRGP, RTOL1, RTOL2, 443 $ W, WORK( INDERR ), WORK( INDGP ), IWORK( IINDBL ), 444 $ IWORK( IINDW ), WORK( INDGRS ), Z, LDZ, 445 $ ISUPPZ, WORK( INDWRK ), IWORK( IINDWK ), IINFO ) 446 IF( IINFO.NE.0 ) THEN 447 INFO = 200 + ABS( IINFO ) 448 RETURN 449 END IF 450 ELSE 451* DLARRE2 computes eigenvalues of the (shifted) root representation 452* DLARRV returns the eigenvalues of the unshifted matrix. 453* However, if the eigenvectors are not desired by the user, we need 454* to apply the corresponding shifts from DLARRE2 to obtain the 455* eigenvalues of the original matrix. 456 DO 20 J = 1, M 457 ITMP = IWORK( IINDBL+J-1 ) 458 W( J ) = W( J ) + E( IWORK( IINSPL+ITMP-1 ) ) 459 20 CONTINUE 460 END IF 461* 462 463* 464* If matrix was scaled, then rescale eigenvalues appropriately. 465* 466 IF( SCALE.NE.ONE ) THEN 467 CALL DSCAL( M, ONE / SCALE, W, 1 ) 468 END IF 469* 470* Correct M if needed 471* 472 IF ( WANTZ ) THEN 473 IF( DOL.NE.1 .OR. DOU.NE.M ) THEN 474 M = DOU - DOL +1 475 ENDIF 476 ENDIF 477* 478* If eigenvalues are not in increasing order, then sort them, 479* possibly along with eigenvectors. 480* 481 IF( NSPLIT.GT.1 ) THEN 482 IF( .NOT. WANTZ ) THEN 483 CALL DLASRT( 'I', DOU - DOL +1, W(DOL), IINFO ) 484 IF( IINFO.NE.0 ) THEN 485 INFO = 3 486 RETURN 487 END IF 488 ELSE 489 DO 60 J = DOL, DOU - 1 490 I = 0 491 TMP = W( J ) 492 DO 50 JJ = J + 1, M 493 IF( W( JJ ).LT.TMP ) THEN 494 I = JJ 495 TMP = W( JJ ) 496 END IF 497 50 CONTINUE 498 IF( I.NE.0 ) THEN 499 W( I ) = W( J ) 500 W( J ) = TMP 501 IF( WANTZ ) THEN 502 CALL DSWAP( N, Z( 1, I-ZOFFSET ), 503 $ 1, Z( 1, J-ZOFFSET ), 1 ) 504 ITMP = ISUPPZ( 2*I-1 ) 505 ISUPPZ( 2*I-1 ) = ISUPPZ( 2*J-1 ) 506 ISUPPZ( 2*J-1 ) = ITMP 507 ITMP = ISUPPZ( 2*I ) 508 ISUPPZ( 2*I ) = ISUPPZ( 2*J ) 509 ISUPPZ( 2*J ) = ITMP 510 END IF 511 END IF 512 60 CONTINUE 513 END IF 514 ENDIF 515* 516 WORK( 1 ) = LWMIN 517 IWORK( 1 ) = LIWMIN 518 RETURN 519* 520* End of DSTEGR2 521* 522 END 523