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