1 SUBROUTINE DLARRB2( N, D, LLD, IFIRST, ILAST, RTOL1, 2 $ RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK, 3 $ PIVMIN, LGPVMN, LGSPDM, TWIST, 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 IMPLICIT NONE 10* 11* .. Scalar Arguments .. 12 INTEGER IFIRST, ILAST, INFO, N, OFFSET, TWIST 13 DOUBLE PRECISION LGPVMN, LGSPDM, PIVMIN, 14 $ RTOL1, RTOL2 15* .. 16* .. Array Arguments .. 17 INTEGER IWORK( * ) 18 DOUBLE PRECISION D( * ), LLD( * ), W( * ), 19 $ WERR( * ), WGAP( * ), WORK( * ) 20* .. 21* 22* Purpose 23* ======= 24* 25* Given the relatively robust representation(RRR) L D L^T, DLARRB2 26* does "limited" bisection to refine the eigenvalues of L D L^T, 27* W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial 28* guesses for these eigenvalues are input in W, the corresponding estimate 29* of the error in these guesses and their gaps are input in WERR 30* and WGAP, respectively. During bisection, intervals 31* [left, right] are maintained by storing their mid-points and 32* semi-widths in the arrays W and WERR respectively. 33* 34* NOTE: 35* There are very few minor differences between DLARRB from LAPACK 36* and this current subroutine DLARRB2. 37* The most important reason for creating this nearly identical copy 38* is profiling: in the ScaLAPACK MRRR algorithm, eigenvalue computation 39* using DLARRB2 is used for refinement in the construction of 40* the representation tree, as opposed to the initial computation of the 41* eigenvalues for the root RRR which uses DLARRB. When profiling, 42* this allows an easy quantification of refinement work vs. computing 43* eigenvalues of the root. 44* 45* Arguments 46* ========= 47* 48* N (input) INTEGER 49* The order of the matrix. 50* 51* D (input) DOUBLE PRECISION array, dimension (N) 52* The N diagonal elements of the diagonal matrix D. 53* 54* LLD (input) DOUBLE PRECISION array, dimension (N-1) 55* The (N-1) elements L(i)*L(i)*D(i). 56* 57* IFIRST (input) INTEGER 58* The index of the first eigenvalue to be computed. 59* 60* ILAST (input) INTEGER 61* The index of the last eigenvalue to be computed. 62* 63* RTOL1 (input) DOUBLE PRECISION 64* RTOL2 (input) DOUBLE PRECISION 65* Tolerance for the convergence of the bisection intervals. 66* An interval [LEFT,RIGHT] has converged if 67* RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) 68* where GAP is the (estimated) distance to the nearest 69* eigenvalue. 70* 71* OFFSET (input) INTEGER 72* Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET 73* through ILAST-OFFSET elements of these arrays are to be used. 74* 75* W (input/output) DOUBLE PRECISION array, dimension (N) 76* On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are 77* estimates of the eigenvalues of L D L^T indexed IFIRST through ILAST. 78* On output, these estimates are refined. 79* 80* WGAP (input/output) DOUBLE PRECISION array, dimension (N-1) 81* On input, the (estimated) gaps between consecutive 82* eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between 83* eigenvalues I and I+1. Note that if IFIRST.EQ.ILAST 84* then WGAP(IFIRST-OFFSET) must be set to ZERO. 85* On output, these gaps are refined. 86* 87* WERR (input/output) DOUBLE PRECISION array, dimension (N) 88* On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are 89* the errors in the estimates of the corresponding elements in W. 90* On output, these errors are refined. 91* 92* WORK (workspace) DOUBLE PRECISION array, dimension (4*N) 93* Workspace. 94* 95* IWORK (workspace) INTEGER array, dimension (2*N) 96* Workspace. 97* 98* PIVMIN (input) DOUBLE PRECISION 99* The minimum pivot in the sturm sequence. 100* 101* LGPVMN (input) DOUBLE PRECISION 102* Logarithm of PIVMIN, precomputed. 103* 104* LGSPDM (input) DOUBLE PRECISION 105* Logarithm of the spectral diameter, precomputed. 106* 107* TWIST (input) INTEGER 108* The twist index for the twisted factorization that is used 109* for the negcount. 110* TWIST = N: Compute negcount from L D L^T - LAMBDA I = L+ D+ L+^T 111* TWIST = 1: Compute negcount from L D L^T - LAMBDA I = U- D- U-^T 112* TWIST = R: Compute negcount from L D L^T - LAMBDA I = N(r) D(r) N(r) 113* 114* INFO (output) INTEGER 115* Error flag. 116* 117* .. Parameters .. 118 DOUBLE PRECISION ZERO, TWO, HALF 119 PARAMETER ( ZERO = 0.0D0, TWO = 2.0D0, 120 $ HALF = 0.5D0 ) 121 INTEGER MAXITR 122* .. 123* .. Local Scalars .. 124 INTEGER I, I1, II, INDLLD, IP, ITER, J, K, NEGCNT, 125 $ NEXT, NINT, OLNINT, PREV, R 126 DOUBLE PRECISION BACK, CVRGD, GAP, LEFT, LGAP, MID, MNWDTH, 127 $ RGAP, RIGHT, SAVGAP, TMP, WIDTH 128 LOGICAL PARANOID 129* .. 130* .. External Functions .. 131 LOGICAL DISNAN 132 DOUBLE PRECISION DLAMCH 133 INTEGER DLANEG2A 134 EXTERNAL DISNAN, DLAMCH, 135 $ DLANEG2A 136* 137* .. 138* .. Intrinsic Functions .. 139 INTRINSIC ABS, MAX, MIN 140* .. 141* .. Executable Statements .. 142* 143 INFO = 0 144* 145* Turn on paranoid check for rounding errors 146* invalidating uncertainty intervals of eigenvalues 147* 148 PARANOID = .TRUE. 149* 150 MAXITR = INT( ( LGSPDM - LGPVMN ) / LOG( TWO ) ) + 2 151 MNWDTH = TWO * PIVMIN 152* 153 R = TWIST 154* 155 INDLLD = 2*N 156 DO 5 J = 1, N-1 157 I=2*J 158 WORK(INDLLD+I-1) = D(J) 159 WORK(INDLLD+I) = LLD(J) 160 5 CONTINUE 161 WORK(INDLLD+2*N-1) = D(N) 162* 163 IF((R.LT.1).OR.(R.GT.N)) R = N 164* 165* Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ]. 166* The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while 167* Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 ) 168* for an unconverged interval is set to the index of the next unconverged 169* interval, and is -1 or 0 for a converged interval. Thus a linked 170* list of unconverged intervals is set up. 171* 172 I1 = IFIRST 173* The number of unconverged intervals 174 NINT = 0 175* The last unconverged interval found 176 PREV = 0 177 178 RGAP = WGAP( I1-OFFSET ) 179 DO 75 I = I1, ILAST 180 K = 2*I 181 II = I - OFFSET 182 LEFT = W( II ) - WERR( II ) 183 RIGHT = W( II ) + WERR( II ) 184 LGAP = RGAP 185 RGAP = WGAP( II ) 186 GAP = MIN( LGAP, RGAP ) 187 188 IF((ABS(LEFT).LE.16*PIVMIN).OR.(ABS(RIGHT).LE.16*PIVMIN)) 189 $ THEN 190 INFO = -1 191 RETURN 192 ENDIF 193 194 IF( PARANOID ) THEN 195* Make sure that [LEFT,RIGHT] contains the desired eigenvalue 196* Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT 197* 198* Do while( NEGCNT(LEFT).GT.I-1 ) 199* 200 BACK = WERR( II ) 201 20 CONTINUE 202 NEGCNT = DLANEG2A( N, WORK(INDLLD+1), LEFT, PIVMIN, R ) 203 IF( NEGCNT.GT.I-1 ) THEN 204 LEFT = LEFT - BACK 205 BACK = TWO*BACK 206 GO TO 20 207 END IF 208* 209* Do while( NEGCNT(RIGHT).LT.I ) 210* Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT 211* 212 BACK = WERR( II ) 213 50 CONTINUE 214 NEGCNT = DLANEG2A( N, WORK(INDLLD+1),RIGHT, PIVMIN, R ) 215 216 IF( NEGCNT.LT.I ) THEN 217 RIGHT = RIGHT + BACK 218 BACK = TWO*BACK 219 GO TO 50 220 END IF 221 ENDIF 222 223 WIDTH = HALF*ABS( LEFT - RIGHT ) 224 TMP = MAX( ABS( LEFT ), ABS( RIGHT ) ) 225 CVRGD = MAX(RTOL1*GAP,RTOL2*TMP) 226 IF( WIDTH.LE.CVRGD .OR. WIDTH.LE.MNWDTH ) THEN 227* This interval has already converged and does not need refinement. 228* (Note that the gaps might change through refining the 229* eigenvalues, however, they can only get bigger.) 230* Remove it from the list. 231 IWORK( K-1 ) = -1 232* Make sure that I1 always points to the first unconverged interval 233 IF((I.EQ.I1).AND.(I.LT.ILAST)) I1 = I + 1 234 IF((PREV.GE.I1).AND.(I.LE.ILAST)) IWORK( 2*PREV-1 ) = I + 1 235 ELSE 236* unconverged interval found 237 PREV = I 238 NINT = NINT + 1 239 IWORK( K-1 ) = I + 1 240 IWORK( K ) = NEGCNT 241 END IF 242 WORK( K-1 ) = LEFT 243 WORK( K ) = RIGHT 244 75 CONTINUE 245 246* 247* Do while( NINT.GT.0 ), i.e. there are still unconverged intervals 248* and while (ITER.LT.MAXITR) 249* 250 ITER = 0 251 80 CONTINUE 252 PREV = I1 - 1 253 I = I1 254 OLNINT = NINT 255 256 DO 100 IP = 1, OLNINT 257 K = 2*I 258 II = I - OFFSET 259 RGAP = WGAP( II ) 260 LGAP = RGAP 261 IF(II.GT.1) LGAP = WGAP( II-1 ) 262 GAP = MIN( LGAP, RGAP ) 263 NEXT = IWORK( K-1 ) 264 LEFT = WORK( K-1 ) 265 RIGHT = WORK( K ) 266 MID = HALF*( LEFT + RIGHT ) 267* semiwidth of interval 268 WIDTH = RIGHT - MID 269 TMP = MAX( ABS( LEFT ), ABS( RIGHT ) ) 270 CVRGD = MAX(RTOL1*GAP,RTOL2*TMP) 271 IF( ( WIDTH.LE.CVRGD ) .OR. ( WIDTH.LE.MNWDTH ).OR. 272 $ ( ITER.EQ.MAXITR ) )THEN 273* reduce number of unconverged intervals 274 NINT = NINT - 1 275* Mark interval as converged. 276 IWORK( K-1 ) = 0 277 IF( I1.EQ.I ) THEN 278 I1 = NEXT 279 ELSE 280* Prev holds the last unconverged interval previously examined 281 IF(PREV.GE.I1) IWORK( 2*PREV-1 ) = NEXT 282 END IF 283 I = NEXT 284 GO TO 100 285 END IF 286 PREV = I 287* 288* Perform one bisection step 289* 290 NEGCNT = DLANEG2A( N, WORK(INDLLD+1), MID, PIVMIN, R ) 291 IF( NEGCNT.LE.I-1 ) THEN 292 WORK( K-1 ) = MID 293 ELSE 294 WORK( K ) = MID 295 END IF 296 I = NEXT 297 100 CONTINUE 298 ITER = ITER + 1 299* do another loop if there are still unconverged intervals 300* However, in the last iteration, all intervals are accepted 301* since this is the best we can do. 302 IF( ( NINT.GT.0 ).AND.(ITER.LE.MAXITR) ) GO TO 80 303* 304* 305* At this point, all the intervals have converged 306* 307* save this gap to restore it after the loop 308 SAVGAP = WGAP( ILAST-OFFSET ) 309* 310 LEFT = WORK( 2*IFIRST-1 ) 311 DO 110 I = IFIRST, ILAST 312 K = 2*I 313 II = I - OFFSET 314* RIGHT is the right boundary of this current interval 315 RIGHT = WORK( K ) 316* All intervals marked by '0' have been refined. 317 IF( IWORK( K-1 ).EQ.0 ) THEN 318 W( II ) = HALF*( LEFT+RIGHT ) 319 WERR( II ) = RIGHT - W( II ) 320 END IF 321* Left is the boundary of the next interval 322 LEFT = WORK( K +1 ) 323 WGAP( II ) = MAX( ZERO, LEFT - RIGHT ) 324 110 CONTINUE 325* restore the last gap which was overwritten by garbage 326 WGAP( ILAST-OFFSET ) = SAVGAP 327 328 RETURN 329* 330* End of DLARRB2 331* 332 END 333* 334* 335* 336 FUNCTION DLANEG2( N, D, LLD, SIGMA, PIVMIN, R ) 337* 338 IMPLICIT NONE 339* 340 INTEGER DLANEG2 341* 342* .. Scalar Arguments .. 343 INTEGER N, R 344 DOUBLE PRECISION PIVMIN, SIGMA 345* .. 346* .. Array Arguments .. 347 DOUBLE PRECISION D( * ), LLD( * ) 348* 349 DOUBLE PRECISION ZERO 350 PARAMETER ( ZERO = 0.0D0 ) 351 352 INTEGER BLKLEN 353 PARAMETER ( BLKLEN = 2048 ) 354* .. 355* .. Local Scalars .. 356 INTEGER BJ, J, NEG1, NEG2, NEGCNT, TO 357 DOUBLE PRECISION DMINUS, DPLUS, GAMMA, P, S, T, TMP, XSAV 358 LOGICAL SAWNAN 359* .. 360* .. External Functions .. 361 LOGICAL DISNAN 362 EXTERNAL DISNAN 363 364 NEGCNT = 0 365* 366* I) upper part: L D L^T - SIGMA I = L+ D+ L+^T 367* run dstqds block-wise to avoid excessive work when NaNs occur 368* 369 S = ZERO 370 DO 210 BJ = 1, R-1, BLKLEN 371 NEG1 = 0 372 XSAV = S 373 TO = BJ+BLKLEN-1 374 IF ( TO.LE.R-1 ) THEN 375 DO 21 J = BJ, TO 376 T = S - SIGMA 377 DPLUS = D( J ) + T 378 IF( DPLUS.LT.ZERO ) NEG1=NEG1 + 1 379 S = T*LLD( J ) / DPLUS 380 21 CONTINUE 381 ELSE 382 DO 22 J = BJ, R-1 383 T = S - SIGMA 384 DPLUS = D( J ) + T 385 IF( DPLUS.LT.ZERO ) NEG1=NEG1 + 1 386 S = T*LLD( J ) / DPLUS 387 22 CONTINUE 388 ENDIF 389 SAWNAN = DISNAN( S ) 390* 391 IF( SAWNAN ) THEN 392 NEG1 = 0 393 S = XSAV 394 TO = BJ+BLKLEN-1 395 IF ( TO.LE.R-1 ) THEN 396 DO 23 J = BJ, TO 397 T = S - SIGMA 398 DPLUS = D( J ) + T 399 IF(ABS(DPLUS).LT.PIVMIN) 400 $ DPLUS = -PIVMIN 401 TMP = LLD( J ) / DPLUS 402 IF( DPLUS.LT.ZERO ) 403 $ NEG1 = NEG1 + 1 404 S = T*TMP 405 IF( TMP.EQ.ZERO ) S = LLD( J ) 406 23 CONTINUE 407 ELSE 408 DO 24 J = BJ, R-1 409 T = S - SIGMA 410 DPLUS = D( J ) + T 411 IF(ABS(DPLUS).LT.PIVMIN) 412 $ DPLUS = -PIVMIN 413 TMP = LLD( J ) / DPLUS 414 IF( DPLUS.LT.ZERO ) NEG1=NEG1+1 415 S = T*TMP 416 IF( TMP.EQ.ZERO ) S = LLD( J ) 417 24 CONTINUE 418 ENDIF 419 END IF 420 NEGCNT = NEGCNT + NEG1 421 210 CONTINUE 422* 423* II) lower part: L D L^T - SIGMA I = U- D- U-^T 424* 425 P = D( N ) - SIGMA 426 DO 230 BJ = N-1, R, -BLKLEN 427 NEG2 = 0 428 XSAV = P 429 TO = BJ-BLKLEN+1 430 IF ( TO.GE.R ) THEN 431 DO 25 J = BJ, TO, -1 432 DMINUS = LLD( J ) + P 433 IF( DMINUS.LT.ZERO ) NEG2=NEG2+1 434 TMP = P / DMINUS 435 P = TMP * D( J ) - SIGMA 436 25 CONTINUE 437 ELSE 438 DO 26 J = BJ, R, -1 439 DMINUS = LLD( J ) + P 440 IF( DMINUS.LT.ZERO ) NEG2=NEG2+1 441 TMP = P / DMINUS 442 P = TMP * D( J ) - SIGMA 443 26 CONTINUE 444 ENDIF 445 SAWNAN = DISNAN( P ) 446* 447 IF( SAWNAN ) THEN 448 NEG2 = 0 449 P = XSAV 450 TO = BJ-BLKLEN+1 451 IF ( TO.GE.R ) THEN 452 DO 27 J = BJ, TO, -1 453 DMINUS = LLD( J ) + P 454 IF(ABS(DMINUS).LT.PIVMIN) 455 $ DMINUS = -PIVMIN 456 TMP = D( J ) / DMINUS 457 IF( DMINUS.LT.ZERO ) 458 $ NEG2 = NEG2 + 1 459 P = P*TMP - SIGMA 460 IF( TMP.EQ.ZERO ) 461 $ P = D( J ) - SIGMA 462 27 CONTINUE 463 ELSE 464 DO 28 J = BJ, R, -1 465 DMINUS = LLD( J ) + P 466 IF(ABS(DMINUS).LT.PIVMIN) 467 $ DMINUS = -PIVMIN 468 TMP = D( J ) / DMINUS 469 IF( DMINUS.LT.ZERO ) 470 $ NEG2 = NEG2 + 1 471 P = P*TMP - SIGMA 472 IF( TMP.EQ.ZERO ) 473 $ P = D( J ) - SIGMA 474 28 CONTINUE 475 ENDIF 476 END IF 477 NEGCNT = NEGCNT + NEG2 478 230 CONTINUE 479* 480* III) Twist index 481* 482 GAMMA = S + P 483 IF( GAMMA.LT.ZERO ) NEGCNT = NEGCNT+1 484 485 DLANEG2 = NEGCNT 486 END 487* 488* 489* 490 FUNCTION DLANEG2A( N, DLLD, SIGMA, PIVMIN, R ) 491* 492 IMPLICIT NONE 493* 494 INTEGER DLANEG2A 495* 496* .. Scalar Arguments .. 497 INTEGER N, R 498 DOUBLE PRECISION PIVMIN, SIGMA 499* .. 500* .. Array Arguments .. 501 DOUBLE PRECISION DLLD( * ) 502* 503 DOUBLE PRECISION ZERO 504 PARAMETER ( ZERO = 0.0D0 ) 505 506 INTEGER BLKLEN 507 PARAMETER ( BLKLEN = 512 ) 508* 509* .. 510* .. Intrinsic Functions .. 511 INTRINSIC INT 512* .. 513* .. Local Scalars .. 514 INTEGER BJ, I, J, NB, NEG1, NEG2, NEGCNT, NX 515 DOUBLE PRECISION DMINUS, DPLUS, GAMMA, P, S, T, TMP, XSAV 516 LOGICAL SAWNAN 517* .. 518* .. External Functions .. 519 LOGICAL DISNAN 520 EXTERNAL DISNAN 521 522 NEGCNT = 0 523* 524* I) upper part: L D L^T - SIGMA I = L+ D+ L+^T 525* run dstqds block-wise to avoid excessive work when NaNs occur, 526* first in chunks of size BLKLEN and then the remainder 527* 528 NB = INT((R-1)/BLKLEN) 529 NX = NB*BLKLEN 530 S = ZERO 531 DO 210 BJ = 1, NX, BLKLEN 532 NEG1 = 0 533 XSAV = S 534 DO 21 J = BJ, BJ+BLKLEN-1 535 I = 2*J 536 T = S - SIGMA 537 DPLUS = DLLD( I-1 ) + T 538 IF( DPLUS.LT.ZERO ) NEG1=NEG1 + 1 539 S = T*DLLD( I ) / DPLUS 540 21 CONTINUE 541 SAWNAN = DISNAN( S ) 542* 543 IF( SAWNAN ) THEN 544 NEG1 = 0 545 S = XSAV 546 DO 23 J = BJ, BJ+BLKLEN-1 547 I = 2*J 548 T = S - SIGMA 549 DPLUS = DLLD( I-1 ) + T 550 IF(ABS(DPLUS).LT.PIVMIN) 551 $ DPLUS = -PIVMIN 552 TMP = DLLD( I ) / DPLUS 553 IF( DPLUS.LT.ZERO ) 554 $ NEG1 = NEG1 + 1 555 S = T*TMP 556 IF( TMP.EQ.ZERO ) S = DLLD( I ) 557 23 CONTINUE 558 END IF 559 NEGCNT = NEGCNT + NEG1 560 210 CONTINUE 561* 562 NEG1 = 0 563 XSAV = S 564 DO 22 J = NX+1, R-1 565 I = 2*J 566 T = S - SIGMA 567 DPLUS = DLLD( I-1 ) + T 568 IF( DPLUS.LT.ZERO ) NEG1=NEG1 + 1 569 S = T*DLLD( I ) / DPLUS 570 22 CONTINUE 571 SAWNAN = DISNAN( S ) 572* 573 IF( SAWNAN ) THEN 574 NEG1 = 0 575 S = XSAV 576 DO 24 J = NX+1, R-1 577 I = 2*J 578 T = S - SIGMA 579 DPLUS = DLLD( I-1 ) + T 580 IF(ABS(DPLUS).LT.PIVMIN) 581 $ DPLUS = -PIVMIN 582 TMP = DLLD( I ) / DPLUS 583 IF( DPLUS.LT.ZERO ) NEG1=NEG1+1 584 S = T*TMP 585 IF( TMP.EQ.ZERO ) S = DLLD( I ) 586 24 CONTINUE 587 ENDIF 588 NEGCNT = NEGCNT + NEG1 589* 590* II) lower part: L D L^T - SIGMA I = U- D- U-^T 591* 592 NB = INT((N-R)/BLKLEN) 593 NX = N-NB*BLKLEN 594 P = DLLD( 2*N-1 ) - SIGMA 595 DO 230 BJ = N-1, NX, -BLKLEN 596 NEG2 = 0 597 XSAV = P 598 DO 25 J = BJ, BJ-BLKLEN+1, -1 599 I = 2*J 600 DMINUS = DLLD( I ) + P 601 IF( DMINUS.LT.ZERO ) NEG2=NEG2+1 602 TMP = P / DMINUS 603 P = TMP * DLLD( I-1 ) - SIGMA 604 25 CONTINUE 605 SAWNAN = DISNAN( P ) 606* 607 IF( SAWNAN ) THEN 608 NEG2 = 0 609 P = XSAV 610 DO 27 J = BJ, BJ-BLKLEN+1, -1 611 I = 2*J 612 DMINUS = DLLD( I ) + P 613 IF(ABS(DMINUS).LT.PIVMIN) 614 $ DMINUS = -PIVMIN 615 TMP = DLLD( I-1 ) / DMINUS 616 IF( DMINUS.LT.ZERO ) 617 $ NEG2 = NEG2 + 1 618 P = P*TMP - SIGMA 619 IF( TMP.EQ.ZERO ) 620 $ P = DLLD( I-1 ) - SIGMA 621 27 CONTINUE 622 END IF 623 NEGCNT = NEGCNT + NEG2 624 230 CONTINUE 625 626 NEG2 = 0 627 XSAV = P 628 DO 26 J = NX-1, R, -1 629 I = 2*J 630 DMINUS = DLLD( I ) + P 631 IF( DMINUS.LT.ZERO ) NEG2=NEG2+1 632 TMP = P / DMINUS 633 P = TMP * DLLD( I-1 ) - SIGMA 634 26 CONTINUE 635 SAWNAN = DISNAN( P ) 636* 637 IF( SAWNAN ) THEN 638 NEG2 = 0 639 P = XSAV 640 DO 28 J = NX-1, R, -1 641 I = 2*J 642 DMINUS = DLLD( I ) + P 643 IF(ABS(DMINUS).LT.PIVMIN) 644 $ DMINUS = -PIVMIN 645 TMP = DLLD( I-1 ) / DMINUS 646 IF( DMINUS.LT.ZERO ) 647 $ NEG2 = NEG2 + 1 648 P = P*TMP - SIGMA 649 IF( TMP.EQ.ZERO ) 650 $ P = DLLD( I-1 ) - SIGMA 651 28 CONTINUE 652 END IF 653 NEGCNT = NEGCNT + NEG2 654* 655* III) Twist index 656* 657 GAMMA = S + P 658 IF( GAMMA.LT.ZERO ) NEGCNT = NEGCNT+1 659 660 DLANEG2A = NEGCNT 661 END 662 663