1 SUBROUTINE DLAEBZ( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL, 2 $ RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT, 3 $ NAB, WORK, IWORK, INFO ) 4* 5* -- LAPACK auxiliary routine (version 3.0) -- 6* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 7* Courant Institute, Argonne National Lab, and Rice University 8* June 30, 1999 9* 10* .. Scalar Arguments .. 11 INTEGER IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX 12 DOUBLE PRECISION ABSTOL, PIVMIN, RELTOL 13* .. 14* .. Array Arguments .. 15 INTEGER IWORK( * ), NAB( MMAX, * ), NVAL( * ) 16 DOUBLE PRECISION AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ), 17 $ WORK( * ) 18* .. 19* 20* Purpose 21* ======= 22* 23* DLAEBZ contains the iteration loops which compute and use the 24* function N(w), which is the count of eigenvalues of a symmetric 25* tridiagonal matrix T less than or equal to its argument w. It 26* performs a choice of two types of loops: 27* 28* IJOB=1, followed by 29* IJOB=2: It takes as input a list of intervals and returns a list of 30* sufficiently small intervals whose union contains the same 31* eigenvalues as the union of the original intervals. 32* The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP. 33* The output interval (AB(j,1),AB(j,2)] will contain 34* eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT. 35* 36* IJOB=3: It performs a binary search in each input interval 37* (AB(j,1),AB(j,2)] for a point w(j) such that 38* N(w(j))=NVAL(j), and uses C(j) as the starting point of 39* the search. If such a w(j) is found, then on output 40* AB(j,1)=AB(j,2)=w. If no such w(j) is found, then on output 41* (AB(j,1),AB(j,2)] will be a small interval containing the 42* point where N(w) jumps through NVAL(j), unless that point 43* lies outside the initial interval. 44* 45* Note that the intervals are in all cases half-open intervals, 46* i.e., of the form (a,b] , which includes b but not a . 47* 48* To avoid underflow, the matrix should be scaled so that its largest 49* element is no greater than overflow**(1/2) * underflow**(1/4) 50* in absolute value. To assure the most accurate computation 51* of small eigenvalues, the matrix should be scaled to be 52* not much smaller than that, either. 53* 54* See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal 55* Matrix", Report CS41, Computer Science Dept., Stanford 56* University, July 21, 1966 57* 58* Note: the arguments are, in general, *not* checked for unreasonable 59* values. 60* 61* Arguments 62* ========= 63* 64* IJOB (input) INTEGER 65* Specifies what is to be done: 66* = 1: Compute NAB for the initial intervals. 67* = 2: Perform bisection iteration to find eigenvalues of T. 68* = 3: Perform bisection iteration to invert N(w), i.e., 69* to find a point which has a specified number of 70* eigenvalues of T to its left. 71* Other values will cause DLAEBZ to return with INFO=-1. 72* 73* NITMAX (input) INTEGER 74* The maximum number of "levels" of bisection to be 75* performed, i.e., an interval of width W will not be made 76* smaller than 2^(-NITMAX) * W. If not all intervals 77* have converged after NITMAX iterations, then INFO is set 78* to the number of non-converged intervals. 79* 80* N (input) INTEGER 81* The dimension n of the tridiagonal matrix T. It must be at 82* least 1. 83* 84* MMAX (input) INTEGER 85* The maximum number of intervals. If more than MMAX intervals 86* are generated, then DLAEBZ will quit with INFO=MMAX+1. 87* 88* MINP (input) INTEGER 89* The initial number of intervals. It may not be greater than 90* MMAX. 91* 92* NBMIN (input) INTEGER 93* The smallest number of intervals that should be processed 94* using a vector loop. If zero, then only the scalar loop 95* will be used. 96* 97* ABSTOL (input) DOUBLE PRECISION 98* The minimum (absolute) width of an interval. When an 99* interval is narrower than ABSTOL, or than RELTOL times the 100* larger (in magnitude) endpoint, then it is considered to be 101* sufficiently small, i.e., converged. This must be at least 102* zero. 103* 104* RELTOL (input) DOUBLE PRECISION 105* The minimum relative width of an interval. When an interval 106* is narrower than ABSTOL, or than RELTOL times the larger (in 107* magnitude) endpoint, then it is considered to be 108* sufficiently small, i.e., converged. Note: this should 109* always be at least radix*machine epsilon. 110* 111* PIVMIN (input) DOUBLE PRECISION 112* The minimum absolute value of a "pivot" in the Sturm 113* sequence loop. This *must* be at least max |e(j)**2| * 114* safe_min and at least safe_min, where safe_min is at least 115* the smallest number that can divide one without overflow. 116* 117* D (input) DOUBLE PRECISION array, dimension (N) 118* The diagonal elements of the tridiagonal matrix T. 119* 120* E (input) DOUBLE PRECISION array, dimension (N) 121* The offdiagonal elements of the tridiagonal matrix T in 122* positions 1 through N-1. E(N) is arbitrary. 123* 124* E2 (input) DOUBLE PRECISION array, dimension (N) 125* The squares of the offdiagonal elements of the tridiagonal 126* matrix T. E2(N) is ignored. 127* 128* NVAL (input/output) INTEGER array, dimension (MINP) 129* If IJOB=1 or 2, not referenced. 130* If IJOB=3, the desired values of N(w). The elements of NVAL 131* will be reordered to correspond with the intervals in AB. 132* Thus, NVAL(j) on output will not, in general be the same as 133* NVAL(j) on input, but it will correspond with the interval 134* (AB(j,1),AB(j,2)] on output. 135* 136* AB (input/output) DOUBLE PRECISION array, dimension (MMAX,2) 137* The endpoints of the intervals. AB(j,1) is a(j), the left 138* endpoint of the j-th interval, and AB(j,2) is b(j), the 139* right endpoint of the j-th interval. The input intervals 140* will, in general, be modified, split, and reordered by the 141* calculation. 142* 143* C (input/output) DOUBLE PRECISION array, dimension (MMAX) 144* If IJOB=1, ignored. 145* If IJOB=2, workspace. 146* If IJOB=3, then on input C(j) should be initialized to the 147* first search point in the binary search. 148* 149* MOUT (output) INTEGER 150* If IJOB=1, the number of eigenvalues in the intervals. 151* If IJOB=2 or 3, the number of intervals output. 152* If IJOB=3, MOUT will equal MINP. 153* 154* NAB (input/output) INTEGER array, dimension (MMAX,2) 155* If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)). 156* If IJOB=2, then on input, NAB(i,j) should be set. It must 157* satisfy the condition: 158* N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)), 159* which means that in interval i only eigenvalues 160* NAB(i,1)+1,...,NAB(i,2) will be considered. Usually, 161* NAB(i,j)=N(AB(i,j)), from a previous call to DLAEBZ with 162* IJOB=1. 163* On output, NAB(i,j) will contain 164* max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of 165* the input interval that the output interval 166* (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the 167* the input values of NAB(k,1) and NAB(k,2). 168* If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)), 169* unless N(w) > NVAL(i) for all search points w , in which 170* case NAB(i,1) will not be modified, i.e., the output 171* value will be the same as the input value (modulo 172* reorderings -- see NVAL and AB), or unless N(w) < NVAL(i) 173* for all search points w , in which case NAB(i,2) will 174* not be modified. Normally, NAB should be set to some 175* distinctive value(s) before DLAEBZ is called. 176* 177* WORK (workspace) DOUBLE PRECISION array, dimension (MMAX) 178* Workspace. 179* 180* IWORK (workspace) INTEGER array, dimension (MMAX) 181* Workspace. 182* 183* INFO (output) INTEGER 184* = 0: All intervals converged. 185* = 1--MMAX: The last INFO intervals did not converge. 186* = MMAX+1: More than MMAX intervals were generated. 187* 188* Further Details 189* =============== 190* 191* This routine is intended to be called only by other LAPACK 192* routines, thus the interface is less user-friendly. It is intended 193* for two purposes: 194* 195* (a) finding eigenvalues. In this case, DLAEBZ should have one or 196* more initial intervals set up in AB, and DLAEBZ should be called 197* with IJOB=1. This sets up NAB, and also counts the eigenvalues. 198* Intervals with no eigenvalues would usually be thrown out at 199* this point. Also, if not all the eigenvalues in an interval i 200* are desired, NAB(i,1) can be increased or NAB(i,2) decreased. 201* For example, set NAB(i,1)=NAB(i,2)-1 to get the largest 202* eigenvalue. DLAEBZ is then called with IJOB=2 and MMAX 203* no smaller than the value of MOUT returned by the call with 204* IJOB=1. After this (IJOB=2) call, eigenvalues NAB(i,1)+1 205* through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the 206* tolerance specified by ABSTOL and RELTOL. 207* 208* (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l). 209* In this case, start with a Gershgorin interval (a,b). Set up 210* AB to contain 2 search intervals, both initially (a,b). One 211* NVAL element should contain f-1 and the other should contain l 212* , while C should contain a and b, resp. NAB(i,1) should be -1 213* and NAB(i,2) should be N+1, to flag an error if the desired 214* interval does not lie in (a,b). DLAEBZ is then called with 215* IJOB=3. On exit, if w(f-1) < w(f), then one of the intervals -- 216* j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while 217* if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r 218* >= 0, then the interval will have N(AB(j,1))=NAB(j,1)=f-k and 219* N(AB(j,2))=NAB(j,2)=f+r. The cases w(l) < w(l+1) and 220* w(l-r)=...=w(l+k) are handled similarly. 221* 222* ===================================================================== 223* 224* .. Parameters .. 225 DOUBLE PRECISION ZERO, TWO, HALF 226 PARAMETER ( ZERO = 0.0D0, TWO = 2.0D0, 227 $ HALF = 1.0D0 / TWO ) 228* .. 229* .. Local Scalars .. 230 INTEGER ITMP1, ITMP2, J, JI, JIT, JP, KF, KFNEW, KL, 231 $ KLNEW 232 DOUBLE PRECISION TMP1, TMP2 233* .. 234* .. Intrinsic Functions .. 235 INTRINSIC ABS, MAX, MIN 236* .. 237* .. Executable Statements .. 238* 239* Check for Errors 240* 241 INFO = 0 242 IF( IJOB.LT.1 .OR. IJOB.GT.3 ) THEN 243 INFO = -1 244 RETURN 245 END IF 246* 247* Initialize NAB 248* 249 IF( IJOB.EQ.1 ) THEN 250* 251* Compute the number of eigenvalues in the initial intervals. 252* 253 MOUT = 0 254*DIR$ NOVECTOR 255 DO 30 JI = 1, MINP 256 DO 20 JP = 1, 2 257 TMP1 = D( 1 ) - AB( JI, JP ) 258 IF( ABS( TMP1 ).LT.PIVMIN ) 259 $ TMP1 = -PIVMIN 260 NAB( JI, JP ) = 0 261 IF( TMP1.LE.ZERO ) 262 $ NAB( JI, JP ) = 1 263* 264 DO 10 J = 2, N 265 TMP1 = D( J ) - E2( J-1 ) / TMP1 - AB( JI, JP ) 266 IF( ABS( TMP1 ).LT.PIVMIN ) 267 $ TMP1 = -PIVMIN 268 IF( TMP1.LE.ZERO ) 269 $ NAB( JI, JP ) = NAB( JI, JP ) + 1 270 10 CONTINUE 271 20 CONTINUE 272 MOUT = MOUT + NAB( JI, 2 ) - NAB( JI, 1 ) 273 30 CONTINUE 274 RETURN 275 END IF 276* 277* Initialize for loop 278* 279* KF and KL have the following meaning: 280* Intervals 1,...,KF-1 have converged. 281* Intervals KF,...,KL still need to be refined. 282* 283 KF = 1 284 KL = MINP 285* 286* If IJOB=2, initialize C. 287* If IJOB=3, use the user-supplied starting point. 288* 289 IF( IJOB.EQ.2 ) THEN 290 DO 40 JI = 1, MINP 291 C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) ) 292 40 CONTINUE 293 END IF 294* 295* Iteration loop 296* 297 DO 130 JIT = 1, NITMAX 298* 299* Loop over intervals 300* 301 IF( KL-KF+1.GE.NBMIN .AND. NBMIN.GT.0 ) THEN 302* 303* Begin of Parallel Version of the loop 304* 305 DO 60 JI = KF, KL 306* 307* Compute N(c), the number of eigenvalues less than c 308* 309 WORK( JI ) = D( 1 ) - C( JI ) 310 IWORK( JI ) = 0 311 IF( WORK( JI ).LE.PIVMIN ) THEN 312 IWORK( JI ) = 1 313 WORK( JI ) = MIN( WORK( JI ), -PIVMIN ) 314 END IF 315* 316 DO 50 J = 2, N 317 WORK( JI ) = D( J ) - E2( J-1 ) / WORK( JI ) - C( JI ) 318 IF( WORK( JI ).LE.PIVMIN ) THEN 319 IWORK( JI ) = IWORK( JI ) + 1 320 WORK( JI ) = MIN( WORK( JI ), -PIVMIN ) 321 END IF 322 50 CONTINUE 323 60 CONTINUE 324* 325 IF( IJOB.LE.2 ) THEN 326* 327* IJOB=2: Choose all intervals containing eigenvalues. 328* 329 KLNEW = KL 330 DO 70 JI = KF, KL 331* 332* Insure that N(w) is monotone 333* 334 IWORK( JI ) = MIN( NAB( JI, 2 ), 335 $ MAX( NAB( JI, 1 ), IWORK( JI ) ) ) 336* 337* Update the Queue -- add intervals if both halves 338* contain eigenvalues. 339* 340 IF( IWORK( JI ).EQ.NAB( JI, 2 ) ) THEN 341* 342* No eigenvalue in the upper interval: 343* just use the lower interval. 344* 345 AB( JI, 2 ) = C( JI ) 346* 347 ELSE IF( IWORK( JI ).EQ.NAB( JI, 1 ) ) THEN 348* 349* No eigenvalue in the lower interval: 350* just use the upper interval. 351* 352 AB( JI, 1 ) = C( JI ) 353 ELSE 354 KLNEW = KLNEW + 1 355 IF( KLNEW.LE.MMAX ) THEN 356* 357* Eigenvalue in both intervals -- add upper to 358* queue. 359* 360 AB( KLNEW, 2 ) = AB( JI, 2 ) 361 NAB( KLNEW, 2 ) = NAB( JI, 2 ) 362 AB( KLNEW, 1 ) = C( JI ) 363 NAB( KLNEW, 1 ) = IWORK( JI ) 364 AB( JI, 2 ) = C( JI ) 365 NAB( JI, 2 ) = IWORK( JI ) 366 ELSE 367 INFO = MMAX + 1 368 END IF 369 END IF 370 70 CONTINUE 371 IF( INFO.NE.0 ) 372 $ RETURN 373 KL = KLNEW 374 ELSE 375* 376* IJOB=3: Binary search. Keep only the interval containing 377* w s.t. N(w) = NVAL 378* 379 DO 80 JI = KF, KL 380 IF( IWORK( JI ).LE.NVAL( JI ) ) THEN 381 AB( JI, 1 ) = C( JI ) 382 NAB( JI, 1 ) = IWORK( JI ) 383 END IF 384 IF( IWORK( JI ).GE.NVAL( JI ) ) THEN 385 AB( JI, 2 ) = C( JI ) 386 NAB( JI, 2 ) = IWORK( JI ) 387 END IF 388 80 CONTINUE 389 END IF 390* 391 ELSE 392* 393* End of Parallel Version of the loop 394* 395* Begin of Serial Version of the loop 396* 397 KLNEW = KL 398 DO 100 JI = KF, KL 399* 400* Compute N(w), the number of eigenvalues less than w 401* 402 TMP1 = C( JI ) 403 TMP2 = D( 1 ) - TMP1 404 ITMP1 = 0 405 IF( TMP2.LE.PIVMIN ) THEN 406 ITMP1 = 1 407 TMP2 = MIN( TMP2, -PIVMIN ) 408 END IF 409* 410* A series of compiler directives to defeat vectorization 411* for the next loop 412* 413*$PL$ CMCHAR=' ' 414CDIR$ NEXTSCALAR 415C$DIR SCALAR 416CDIR$ NEXT SCALAR 417CVD$L NOVECTOR 418CDEC$ NOVECTOR 419CVD$ NOVECTOR 420*VDIR NOVECTOR 421*VOCL LOOP,SCALAR 422CIBM PREFER SCALAR 423*$PL$ CMCHAR='*' 424* 425 DO 90 J = 2, N 426 TMP2 = D( J ) - E2( J-1 ) / TMP2 - TMP1 427 IF( TMP2.LE.PIVMIN ) THEN 428 ITMP1 = ITMP1 + 1 429 TMP2 = MIN( TMP2, -PIVMIN ) 430 END IF 431 90 CONTINUE 432* 433 IF( IJOB.LE.2 ) THEN 434* 435* IJOB=2: Choose all intervals containing eigenvalues. 436* 437* Insure that N(w) is monotone 438* 439 ITMP1 = MIN( NAB( JI, 2 ), 440 $ MAX( NAB( JI, 1 ), ITMP1 ) ) 441* 442* Update the Queue -- add intervals if both halves 443* contain eigenvalues. 444* 445 IF( ITMP1.EQ.NAB( JI, 2 ) ) THEN 446* 447* No eigenvalue in the upper interval: 448* just use the lower interval. 449* 450 AB( JI, 2 ) = TMP1 451* 452 ELSE IF( ITMP1.EQ.NAB( JI, 1 ) ) THEN 453* 454* No eigenvalue in the lower interval: 455* just use the upper interval. 456* 457 AB( JI, 1 ) = TMP1 458 ELSE IF( KLNEW.LT.MMAX ) THEN 459* 460* Eigenvalue in both intervals -- add upper to queue. 461* 462 KLNEW = KLNEW + 1 463 AB( KLNEW, 2 ) = AB( JI, 2 ) 464 NAB( KLNEW, 2 ) = NAB( JI, 2 ) 465 AB( KLNEW, 1 ) = TMP1 466 NAB( KLNEW, 1 ) = ITMP1 467 AB( JI, 2 ) = TMP1 468 NAB( JI, 2 ) = ITMP1 469 ELSE 470 INFO = MMAX + 1 471 RETURN 472 END IF 473 ELSE 474* 475* IJOB=3: Binary search. Keep only the interval 476* containing w s.t. N(w) = NVAL 477* 478 IF( ITMP1.LE.NVAL( JI ) ) THEN 479 AB( JI, 1 ) = TMP1 480 NAB( JI, 1 ) = ITMP1 481 END IF 482 IF( ITMP1.GE.NVAL( JI ) ) THEN 483 AB( JI, 2 ) = TMP1 484 NAB( JI, 2 ) = ITMP1 485 END IF 486 END IF 487 100 CONTINUE 488 KL = KLNEW 489* 490* End of Serial Version of the loop 491* 492 END IF 493* 494* Check for convergence 495* 496 KFNEW = KF 497 DO 110 JI = KF, KL 498 TMP1 = ABS( AB( JI, 2 )-AB( JI, 1 ) ) 499 TMP2 = MAX( ABS( AB( JI, 2 ) ), ABS( AB( JI, 1 ) ) ) 500 IF( TMP1.LT.MAX( ABSTOL, PIVMIN, RELTOL*TMP2 ) .OR. 501 $ NAB( JI, 1 ).GE.NAB( JI, 2 ) ) THEN 502* 503* Converged -- Swap with position KFNEW, 504* then increment KFNEW 505* 506 IF( JI.GT.KFNEW ) THEN 507 TMP1 = AB( JI, 1 ) 508 TMP2 = AB( JI, 2 ) 509 ITMP1 = NAB( JI, 1 ) 510 ITMP2 = NAB( JI, 2 ) 511 AB( JI, 1 ) = AB( KFNEW, 1 ) 512 AB( JI, 2 ) = AB( KFNEW, 2 ) 513 NAB( JI, 1 ) = NAB( KFNEW, 1 ) 514 NAB( JI, 2 ) = NAB( KFNEW, 2 ) 515 AB( KFNEW, 1 ) = TMP1 516 AB( KFNEW, 2 ) = TMP2 517 NAB( KFNEW, 1 ) = ITMP1 518 NAB( KFNEW, 2 ) = ITMP2 519 IF( IJOB.EQ.3 ) THEN 520 ITMP1 = NVAL( JI ) 521 NVAL( JI ) = NVAL( KFNEW ) 522 NVAL( KFNEW ) = ITMP1 523 END IF 524 END IF 525 KFNEW = KFNEW + 1 526 END IF 527 110 CONTINUE 528 KF = KFNEW 529* 530* Choose Midpoints 531* 532 DO 120 JI = KF, KL 533 C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) ) 534 120 CONTINUE 535* 536* If no more intervals to refine, quit. 537* 538 IF( KF.GT.KL ) 539 $ GO TO 140 540 130 CONTINUE 541* 542* Converged 543* 544 140 CONTINUE 545 INFO = MAX( KL+1-KF, 0 ) 546 MOUT = KL 547* 548 RETURN 549* 550* End of DLAEBZ 551* 552 END 553