1 SUBROUTINE SLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO ) 2* 3* -- LAPACK routine (version 3.0) -- 4* Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab, 5* Courant Institute, NAG Ltd., and Rice University 6* December 23, 1999 7* 8* .. Scalar Arguments .. 9 INTEGER I, INFO, N 10 REAL DLAM, RHO 11* .. 12* .. Array Arguments .. 13 REAL D( * ), DELTA( * ), Z( * ) 14* .. 15* 16* Purpose 17* ======= 18* 19* This subroutine computes the I-th updated eigenvalue of a symmetric 20* rank-one modification to a diagonal matrix whose elements are 21* given in the array d, and that 22* 23* D(i) < D(j) for i < j 24* 25* and that RHO > 0. This is arranged by the calling routine, and is 26* no loss in generality. The rank-one modified system is thus 27* 28* diag( D ) + RHO * Z * Z_transpose. 29* 30* where we assume the Euclidean norm of Z is 1. 31* 32* The method consists of approximating the rational functions in the 33* secular equation by simpler interpolating rational functions. 34* 35* Arguments 36* ========= 37* 38* N (input) INTEGER 39* The length of all arrays. 40* 41* I (input) INTEGER 42* The index of the eigenvalue to be computed. 1 <= I <= N. 43* 44* D (input) REAL array, dimension (N) 45* The original eigenvalues. It is assumed that they are in 46* order, D(I) < D(J) for I < J. 47* 48* Z (input) REAL array, dimension (N) 49* The components of the updating vector. 50* 51* DELTA (output) REAL array, dimension (N) 52* If N .ne. 1, DELTA contains (D(j) - lambda_I) in its j-th 53* component. If N = 1, then DELTA(1) = 1. The vector DELTA 54* contains the information necessary to construct the 55* eigenvectors. 56* 57* RHO (input) REAL 58* The scalar in the symmetric updating formula. 59* 60* DLAM (output) REAL 61* The computed lambda_I, the I-th updated eigenvalue. 62* 63* INFO (output) INTEGER 64* = 0: successful exit 65* > 0: if INFO = 1, the updating process failed. 66* 67* Internal Parameters 68* =================== 69* 70* Logical variable ORGATI (origin-at-i?) is used for distinguishing 71* whether D(i) or D(i+1) is treated as the origin. 72* 73* ORGATI = .true. origin at i 74* ORGATI = .false. origin at i+1 75* 76* Logical variable SWTCH3 (switch-for-3-poles?) is for noting 77* if we are working with THREE poles! 78* 79* MAXIT is the maximum number of iterations allowed for each 80* eigenvalue. 81* 82* Further Details 83* =============== 84* 85* Based on contributions by 86* Ren-Cang Li, Computer Science Division, University of California 87* at Berkeley, USA 88* 89* ===================================================================== 90* 91* .. Parameters .. 92 INTEGER MAXIT 93 PARAMETER ( MAXIT = 30 ) 94 REAL ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN 95 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0, 96 $ THREE = 3.0E0, FOUR = 4.0E0, EIGHT = 8.0E0, 97 $ TEN = 10.0E0 ) 98* .. 99* .. Local Scalars .. 100 LOGICAL ORGATI, SWTCH, SWTCH3 101 INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER 102 REAL A, B, C, DEL, DLTLB, DLTUB, DPHI, DPSI, DW, 103 $ EPS, ERRETM, ETA, MIDPT, PHI, PREW, PSI, 104 $ RHOINV, TAU, TEMP, TEMP1, W 105* .. 106* .. Local Arrays .. 107 REAL ZZ( 3 ) 108* .. 109* .. External Functions .. 110 REAL SLAMCH 111 EXTERNAL SLAMCH 112* .. 113* .. External Subroutines .. 114 EXTERNAL SLAED5, SLAED6 115* .. 116* .. Intrinsic Functions .. 117 INTRINSIC ABS, MAX, MIN, SQRT 118* .. 119* .. Executable Statements .. 120* 121* Since this routine is called in an inner loop, we do no argument 122* checking. 123* 124* Quick return for N=1 and 2. 125* 126 INFO = 0 127 IF( N.EQ.1 ) THEN 128* 129* Presumably, I=1 upon entry 130* 131 DLAM = D( 1 ) + RHO*Z( 1 )*Z( 1 ) 132 DELTA( 1 ) = ONE 133 RETURN 134 END IF 135 IF( N.EQ.2 ) THEN 136 CALL SLAED5( I, D, Z, DELTA, RHO, DLAM ) 137 RETURN 138 END IF 139* 140* Compute machine epsilon 141* 142 EPS = SLAMCH( 'Epsilon' ) 143 RHOINV = ONE / RHO 144* 145* The case I = N 146* 147 IF( I.EQ.N ) THEN 148* 149* Initialize some basic variables 150* 151 II = N - 1 152 NITER = 1 153* 154* Calculate initial guess 155* 156 MIDPT = RHO / TWO 157* 158* If ||Z||_2 is not one, then TEMP should be set to 159* RHO * ||Z||_2^2 / TWO 160* 161 DO 10 J = 1, N 162 DELTA( J ) = ( D( J )-D( I ) ) - MIDPT 163 10 CONTINUE 164* 165 PSI = ZERO 166 DO 20 J = 1, N - 2 167 PSI = PSI + Z( J )*Z( J ) / DELTA( J ) 168 20 CONTINUE 169* 170 C = RHOINV + PSI 171 W = C + Z( II )*Z( II ) / DELTA( II ) + 172 $ Z( N )*Z( N ) / DELTA( N ) 173* 174 IF( W.LE.ZERO ) THEN 175 TEMP = Z( N-1 )*Z( N-1 ) / ( D( N )-D( N-1 )+RHO ) + 176 $ Z( N )*Z( N ) / RHO 177 IF( C.LE.TEMP ) THEN 178 TAU = RHO 179 ELSE 180 DEL = D( N ) - D( N-1 ) 181 A = -C*DEL + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N ) 182 B = Z( N )*Z( N )*DEL 183 IF( A.LT.ZERO ) THEN 184 TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A ) 185 ELSE 186 TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C ) 187 END IF 188 END IF 189* 190* It can be proved that 191* D(N)+RHO/2 <= LAMBDA(N) < D(N)+TAU <= D(N)+RHO 192* 193 DLTLB = MIDPT 194 DLTUB = RHO 195 ELSE 196 DEL = D( N ) - D( N-1 ) 197 A = -C*DEL + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N ) 198 B = Z( N )*Z( N )*DEL 199 IF( A.LT.ZERO ) THEN 200 TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A ) 201 ELSE 202 TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C ) 203 END IF 204* 205* It can be proved that 206* D(N) < D(N)+TAU < LAMBDA(N) < D(N)+RHO/2 207* 208 DLTLB = ZERO 209 DLTUB = MIDPT 210 END IF 211* 212 DO 30 J = 1, N 213 DELTA( J ) = ( D( J )-D( I ) ) - TAU 214 30 CONTINUE 215* 216* Evaluate PSI and the derivative DPSI 217* 218 DPSI = ZERO 219 PSI = ZERO 220 ERRETM = ZERO 221 DO 40 J = 1, II 222 TEMP = Z( J ) / DELTA( J ) 223 PSI = PSI + Z( J )*TEMP 224 DPSI = DPSI + TEMP*TEMP 225 ERRETM = ERRETM + PSI 226 40 CONTINUE 227 ERRETM = ABS( ERRETM ) 228* 229* Evaluate PHI and the derivative DPHI 230* 231 TEMP = Z( N ) / DELTA( N ) 232 PHI = Z( N )*TEMP 233 DPHI = TEMP*TEMP 234 ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV + 235 $ ABS( TAU )*( DPSI+DPHI ) 236* 237 W = RHOINV + PHI + PSI 238* 239* Test for convergence 240* 241 IF( ABS( W ).LE.EPS*ERRETM ) THEN 242 DLAM = D( I ) + TAU 243 GO TO 250 244 END IF 245* 246 IF( W.LE.ZERO ) THEN 247 DLTLB = MAX( DLTLB, TAU ) 248 ELSE 249 DLTUB = MIN( DLTUB, TAU ) 250 END IF 251* 252* Calculate the new step 253* 254 NITER = NITER + 1 255 C = W - DELTA( N-1 )*DPSI - DELTA( N )*DPHI 256 A = ( DELTA( N-1 )+DELTA( N ) )*W - 257 $ DELTA( N-1 )*DELTA( N )*( DPSI+DPHI ) 258 B = DELTA( N-1 )*DELTA( N )*W 259 IF( C.LT.ZERO ) 260 $ C = ABS( C ) 261 IF( C.EQ.ZERO ) THEN 262* ETA = B/A 263* ETA = RHO - TAU 264 ETA = DLTUB - TAU 265 ELSE IF( A.GE.ZERO ) THEN 266 ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) 267 ELSE 268 ETA = TWO*B / ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) 269 END IF 270* 271* Note, eta should be positive if w is negative, and 272* eta should be negative otherwise. However, 273* if for some reason caused by roundoff, eta*w > 0, 274* we simply use one Newton step instead. This way 275* will guarantee eta*w < 0. 276* 277 IF( W*ETA.GT.ZERO ) 278 $ ETA = -W / ( DPSI+DPHI ) 279 TEMP = TAU + ETA 280 IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN 281 IF( W.LT.ZERO ) THEN 282 ETA = ( DLTUB-TAU ) / TWO 283 ELSE 284 ETA = ( DLTLB-TAU ) / TWO 285 END IF 286 END IF 287 DO 50 J = 1, N 288 DELTA( J ) = DELTA( J ) - ETA 289 50 CONTINUE 290* 291 TAU = TAU + ETA 292* 293* Evaluate PSI and the derivative DPSI 294* 295 DPSI = ZERO 296 PSI = ZERO 297 ERRETM = ZERO 298 DO 60 J = 1, II 299 TEMP = Z( J ) / DELTA( J ) 300 PSI = PSI + Z( J )*TEMP 301 DPSI = DPSI + TEMP*TEMP 302 ERRETM = ERRETM + PSI 303 60 CONTINUE 304 ERRETM = ABS( ERRETM ) 305* 306* Evaluate PHI and the derivative DPHI 307* 308 TEMP = Z( N ) / DELTA( N ) 309 PHI = Z( N )*TEMP 310 DPHI = TEMP*TEMP 311 ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV + 312 $ ABS( TAU )*( DPSI+DPHI ) 313* 314 W = RHOINV + PHI + PSI 315* 316* Main loop to update the values of the array DELTA 317* 318 ITER = NITER + 1 319* 320 DO 90 NITER = ITER, MAXIT 321* 322* Test for convergence 323* 324 IF( ABS( W ).LE.EPS*ERRETM ) THEN 325 DLAM = D( I ) + TAU 326 GO TO 250 327 END IF 328* 329 IF( W.LE.ZERO ) THEN 330 DLTLB = MAX( DLTLB, TAU ) 331 ELSE 332 DLTUB = MIN( DLTUB, TAU ) 333 END IF 334* 335* Calculate the new step 336* 337 C = W - DELTA( N-1 )*DPSI - DELTA( N )*DPHI 338 A = ( DELTA( N-1 )+DELTA( N ) )*W - 339 $ DELTA( N-1 )*DELTA( N )*( DPSI+DPHI ) 340 B = DELTA( N-1 )*DELTA( N )*W 341 IF( A.GE.ZERO ) THEN 342 ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) 343 ELSE 344 ETA = TWO*B / ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) 345 END IF 346* 347* Note, eta should be positive if w is negative, and 348* eta should be negative otherwise. However, 349* if for some reason caused by roundoff, eta*w > 0, 350* we simply use one Newton step instead. This way 351* will guarantee eta*w < 0. 352* 353 IF( W*ETA.GT.ZERO ) 354 $ ETA = -W / ( DPSI+DPHI ) 355 TEMP = TAU + ETA 356 IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN 357 IF( W.LT.ZERO ) THEN 358 ETA = ( DLTUB-TAU ) / TWO 359 ELSE 360 ETA = ( DLTLB-TAU ) / TWO 361 END IF 362 END IF 363 DO 70 J = 1, N 364 DELTA( J ) = DELTA( J ) - ETA 365 70 CONTINUE 366* 367 TAU = TAU + ETA 368* 369* Evaluate PSI and the derivative DPSI 370* 371 DPSI = ZERO 372 PSI = ZERO 373 ERRETM = ZERO 374 DO 80 J = 1, II 375 TEMP = Z( J ) / DELTA( J ) 376 PSI = PSI + Z( J )*TEMP 377 DPSI = DPSI + TEMP*TEMP 378 ERRETM = ERRETM + PSI 379 80 CONTINUE 380 ERRETM = ABS( ERRETM ) 381* 382* Evaluate PHI and the derivative DPHI 383* 384 TEMP = Z( N ) / DELTA( N ) 385 PHI = Z( N )*TEMP 386 DPHI = TEMP*TEMP 387 ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV + 388 $ ABS( TAU )*( DPSI+DPHI ) 389* 390 W = RHOINV + PHI + PSI 391 90 CONTINUE 392* 393* Return with INFO = 1, NITER = MAXIT and not converged 394* 395 INFO = 1 396 DLAM = D( I ) + TAU 397 GO TO 250 398* 399* End for the case I = N 400* 401 ELSE 402* 403* The case for I < N 404* 405 NITER = 1 406 IP1 = I + 1 407* 408* Calculate initial guess 409* 410 DEL = D( IP1 ) - D( I ) 411 MIDPT = DEL / TWO 412 DO 100 J = 1, N 413 DELTA( J ) = ( D( J )-D( I ) ) - MIDPT 414 100 CONTINUE 415* 416 PSI = ZERO 417 DO 110 J = 1, I - 1 418 PSI = PSI + Z( J )*Z( J ) / DELTA( J ) 419 110 CONTINUE 420* 421 PHI = ZERO 422 DO 120 J = N, I + 2, -1 423 PHI = PHI + Z( J )*Z( J ) / DELTA( J ) 424 120 CONTINUE 425 C = RHOINV + PSI + PHI 426 W = C + Z( I )*Z( I ) / DELTA( I ) + 427 $ Z( IP1 )*Z( IP1 ) / DELTA( IP1 ) 428* 429 IF( W.GT.ZERO ) THEN 430* 431* d(i)< the ith eigenvalue < (d(i)+d(i+1))/2 432* 433* We choose d(i) as origin. 434* 435 ORGATI = .TRUE. 436 A = C*DEL + Z( I )*Z( I ) + Z( IP1 )*Z( IP1 ) 437 B = Z( I )*Z( I )*DEL 438 IF( A.GT.ZERO ) THEN 439 TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) 440 ELSE 441 TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) 442 END IF 443 DLTLB = ZERO 444 DLTUB = MIDPT 445 ELSE 446* 447* (d(i)+d(i+1))/2 <= the ith eigenvalue < d(i+1) 448* 449* We choose d(i+1) as origin. 450* 451 ORGATI = .FALSE. 452 A = C*DEL - Z( I )*Z( I ) - Z( IP1 )*Z( IP1 ) 453 B = Z( IP1 )*Z( IP1 )*DEL 454 IF( A.LT.ZERO ) THEN 455 TAU = TWO*B / ( A-SQRT( ABS( A*A+FOUR*B*C ) ) ) 456 ELSE 457 TAU = -( A+SQRT( ABS( A*A+FOUR*B*C ) ) ) / ( TWO*C ) 458 END IF 459 DLTLB = -MIDPT 460 DLTUB = ZERO 461 END IF 462* 463 IF( ORGATI ) THEN 464 DO 130 J = 1, N 465 DELTA( J ) = ( D( J )-D( I ) ) - TAU 466 130 CONTINUE 467 ELSE 468 DO 140 J = 1, N 469 DELTA( J ) = ( D( J )-D( IP1 ) ) - TAU 470 140 CONTINUE 471 END IF 472 IF( ORGATI ) THEN 473 II = I 474 ELSE 475 II = I + 1 476 END IF 477 IIM1 = II - 1 478 IIP1 = II + 1 479* 480* Evaluate PSI and the derivative DPSI 481* 482 DPSI = ZERO 483 PSI = ZERO 484 ERRETM = ZERO 485 DO 150 J = 1, IIM1 486 TEMP = Z( J ) / DELTA( J ) 487 PSI = PSI + Z( J )*TEMP 488 DPSI = DPSI + TEMP*TEMP 489 ERRETM = ERRETM + PSI 490 150 CONTINUE 491 ERRETM = ABS( ERRETM ) 492* 493* Evaluate PHI and the derivative DPHI 494* 495 DPHI = ZERO 496 PHI = ZERO 497 DO 160 J = N, IIP1, -1 498 TEMP = Z( J ) / DELTA( J ) 499 PHI = PHI + Z( J )*TEMP 500 DPHI = DPHI + TEMP*TEMP 501 ERRETM = ERRETM + PHI 502 160 CONTINUE 503* 504 W = RHOINV + PHI + PSI 505* 506* W is the value of the secular function with 507* its ii-th element removed. 508* 509 SWTCH3 = .FALSE. 510 IF( ORGATI ) THEN 511 IF( W.LT.ZERO ) 512 $ SWTCH3 = .TRUE. 513 ELSE 514 IF( W.GT.ZERO ) 515 $ SWTCH3 = .TRUE. 516 END IF 517 IF( II.EQ.1 .OR. II.EQ.N ) 518 $ SWTCH3 = .FALSE. 519* 520 TEMP = Z( II ) / DELTA( II ) 521 DW = DPSI + DPHI + TEMP*TEMP 522 TEMP = Z( II )*TEMP 523 W = W + TEMP 524 ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV + 525 $ THREE*ABS( TEMP ) + ABS( TAU )*DW 526* 527* Test for convergence 528* 529 IF( ABS( W ).LE.EPS*ERRETM ) THEN 530 IF( ORGATI ) THEN 531 DLAM = D( I ) + TAU 532 ELSE 533 DLAM = D( IP1 ) + TAU 534 END IF 535 GO TO 250 536 END IF 537* 538 IF( W.LE.ZERO ) THEN 539 DLTLB = MAX( DLTLB, TAU ) 540 ELSE 541 DLTUB = MIN( DLTUB, TAU ) 542 END IF 543* 544* Calculate the new step 545* 546 NITER = NITER + 1 547 IF( .NOT.SWTCH3 ) THEN 548 IF( ORGATI ) THEN 549 C = W - DELTA( IP1 )*DW - ( D( I )-D( IP1 ) )* 550 $ ( Z( I ) / DELTA( I ) )**2 551 ELSE 552 C = W - DELTA( I )*DW - ( D( IP1 )-D( I ) )* 553 $ ( Z( IP1 ) / DELTA( IP1 ) )**2 554 END IF 555 A = ( DELTA( I )+DELTA( IP1 ) )*W - 556 $ DELTA( I )*DELTA( IP1 )*DW 557 B = DELTA( I )*DELTA( IP1 )*W 558 IF( C.EQ.ZERO ) THEN 559 IF( A.EQ.ZERO ) THEN 560 IF( ORGATI ) THEN 561 A = Z( I )*Z( I ) + DELTA( IP1 )*DELTA( IP1 )* 562 $ ( DPSI+DPHI ) 563 ELSE 564 A = Z( IP1 )*Z( IP1 ) + DELTA( I )*DELTA( I )* 565 $ ( DPSI+DPHI ) 566 END IF 567 END IF 568 ETA = B / A 569 ELSE IF( A.LE.ZERO ) THEN 570 ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) 571 ELSE 572 ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) 573 END IF 574 ELSE 575* 576* Interpolation using THREE most relevant poles 577* 578 TEMP = RHOINV + PSI + PHI 579 IF( ORGATI ) THEN 580 TEMP1 = Z( IIM1 ) / DELTA( IIM1 ) 581 TEMP1 = TEMP1*TEMP1 582 C = TEMP - DELTA( IIP1 )*( DPSI+DPHI ) - 583 $ ( D( IIM1 )-D( IIP1 ) )*TEMP1 584 ZZ( 1 ) = Z( IIM1 )*Z( IIM1 ) 585 ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )* 586 $ ( ( DPSI-TEMP1 )+DPHI ) 587 ELSE 588 TEMP1 = Z( IIP1 ) / DELTA( IIP1 ) 589 TEMP1 = TEMP1*TEMP1 590 C = TEMP - DELTA( IIM1 )*( DPSI+DPHI ) - 591 $ ( D( IIP1 )-D( IIM1 ) )*TEMP1 592 ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )* 593 $ ( DPSI+( DPHI-TEMP1 ) ) 594 ZZ( 3 ) = Z( IIP1 )*Z( IIP1 ) 595 END IF 596 ZZ( 2 ) = Z( II )*Z( II ) 597 CALL SLAED6( NITER, ORGATI, C, DELTA( IIM1 ), ZZ, W, ETA, 598 $ INFO ) 599 IF( INFO.NE.0 ) 600 $ GO TO 250 601 END IF 602* 603* Note, eta should be positive if w is negative, and 604* eta should be negative otherwise. However, 605* if for some reason caused by roundoff, eta*w > 0, 606* we simply use one Newton step instead. This way 607* will guarantee eta*w < 0. 608* 609 IF( W*ETA.GE.ZERO ) 610 $ ETA = -W / DW 611 TEMP = TAU + ETA 612 IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN 613 IF( W.LT.ZERO ) THEN 614 ETA = ( DLTUB-TAU ) / TWO 615 ELSE 616 ETA = ( DLTLB-TAU ) / TWO 617 END IF 618 END IF 619* 620 PREW = W 621* 622 170 CONTINUE 623 DO 180 J = 1, N 624 DELTA( J ) = DELTA( J ) - ETA 625 180 CONTINUE 626* 627* Evaluate PSI and the derivative DPSI 628* 629 DPSI = ZERO 630 PSI = ZERO 631 ERRETM = ZERO 632 DO 190 J = 1, IIM1 633 TEMP = Z( J ) / DELTA( J ) 634 PSI = PSI + Z( J )*TEMP 635 DPSI = DPSI + TEMP*TEMP 636 ERRETM = ERRETM + PSI 637 190 CONTINUE 638 ERRETM = ABS( ERRETM ) 639* 640* Evaluate PHI and the derivative DPHI 641* 642 DPHI = ZERO 643 PHI = ZERO 644 DO 200 J = N, IIP1, -1 645 TEMP = Z( J ) / DELTA( J ) 646 PHI = PHI + Z( J )*TEMP 647 DPHI = DPHI + TEMP*TEMP 648 ERRETM = ERRETM + PHI 649 200 CONTINUE 650* 651 TEMP = Z( II ) / DELTA( II ) 652 DW = DPSI + DPHI + TEMP*TEMP 653 TEMP = Z( II )*TEMP 654 W = RHOINV + PHI + PSI + TEMP 655 ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV + 656 $ THREE*ABS( TEMP ) + ABS( TAU+ETA )*DW 657* 658 SWTCH = .FALSE. 659 IF( ORGATI ) THEN 660 IF( -W.GT.ABS( PREW ) / TEN ) 661 $ SWTCH = .TRUE. 662 ELSE 663 IF( W.GT.ABS( PREW ) / TEN ) 664 $ SWTCH = .TRUE. 665 END IF 666* 667 TAU = TAU + ETA 668* 669* Main loop to update the values of the array DELTA 670* 671 ITER = NITER + 1 672* 673 DO 240 NITER = ITER, MAXIT 674* 675* Test for convergence 676* 677 IF( ABS( W ).LE.EPS*ERRETM ) THEN 678 IF( ORGATI ) THEN 679 DLAM = D( I ) + TAU 680 ELSE 681 DLAM = D( IP1 ) + TAU 682 END IF 683 GO TO 250 684 END IF 685* 686 IF( W.LE.ZERO ) THEN 687 DLTLB = MAX( DLTLB, TAU ) 688 ELSE 689 DLTUB = MIN( DLTUB, TAU ) 690 END IF 691* 692* Calculate the new step 693* 694 IF( .NOT.SWTCH3 ) THEN 695 IF( .NOT.SWTCH ) THEN 696 IF( ORGATI ) THEN 697 C = W - DELTA( IP1 )*DW - 698 $ ( D( I )-D( IP1 ) )*( Z( I ) / DELTA( I ) )**2 699 ELSE 700 C = W - DELTA( I )*DW - ( D( IP1 )-D( I ) )* 701 $ ( Z( IP1 ) / DELTA( IP1 ) )**2 702 END IF 703 ELSE 704 TEMP = Z( II ) / DELTA( II ) 705 IF( ORGATI ) THEN 706 DPSI = DPSI + TEMP*TEMP 707 ELSE 708 DPHI = DPHI + TEMP*TEMP 709 END IF 710 C = W - DELTA( I )*DPSI - DELTA( IP1 )*DPHI 711 END IF 712 A = ( DELTA( I )+DELTA( IP1 ) )*W - 713 $ DELTA( I )*DELTA( IP1 )*DW 714 B = DELTA( I )*DELTA( IP1 )*W 715 IF( C.EQ.ZERO ) THEN 716 IF( A.EQ.ZERO ) THEN 717 IF( .NOT.SWTCH ) THEN 718 IF( ORGATI ) THEN 719 A = Z( I )*Z( I ) + DELTA( IP1 )* 720 $ DELTA( IP1 )*( DPSI+DPHI ) 721 ELSE 722 A = Z( IP1 )*Z( IP1 ) + 723 $ DELTA( I )*DELTA( I )*( DPSI+DPHI ) 724 END IF 725 ELSE 726 A = DELTA( I )*DELTA( I )*DPSI + 727 $ DELTA( IP1 )*DELTA( IP1 )*DPHI 728 END IF 729 END IF 730 ETA = B / A 731 ELSE IF( A.LE.ZERO ) THEN 732 ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) 733 ELSE 734 ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) 735 END IF 736 ELSE 737* 738* Interpolation using THREE most relevant poles 739* 740 TEMP = RHOINV + PSI + PHI 741 IF( SWTCH ) THEN 742 C = TEMP - DELTA( IIM1 )*DPSI - DELTA( IIP1 )*DPHI 743 ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )*DPSI 744 ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )*DPHI 745 ELSE 746 IF( ORGATI ) THEN 747 TEMP1 = Z( IIM1 ) / DELTA( IIM1 ) 748 TEMP1 = TEMP1*TEMP1 749 C = TEMP - DELTA( IIP1 )*( DPSI+DPHI ) - 750 $ ( D( IIM1 )-D( IIP1 ) )*TEMP1 751 ZZ( 1 ) = Z( IIM1 )*Z( IIM1 ) 752 ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )* 753 $ ( ( DPSI-TEMP1 )+DPHI ) 754 ELSE 755 TEMP1 = Z( IIP1 ) / DELTA( IIP1 ) 756 TEMP1 = TEMP1*TEMP1 757 C = TEMP - DELTA( IIM1 )*( DPSI+DPHI ) - 758 $ ( D( IIP1 )-D( IIM1 ) )*TEMP1 759 ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )* 760 $ ( DPSI+( DPHI-TEMP1 ) ) 761 ZZ( 3 ) = Z( IIP1 )*Z( IIP1 ) 762 END IF 763 END IF 764 CALL SLAED6( NITER, ORGATI, C, DELTA( IIM1 ), ZZ, W, ETA, 765 $ INFO ) 766 IF( INFO.NE.0 ) 767 $ GO TO 250 768 END IF 769* 770* Note, eta should be positive if w is negative, and 771* eta should be negative otherwise. However, 772* if for some reason caused by roundoff, eta*w > 0, 773* we simply use one Newton step instead. This way 774* will guarantee eta*w < 0. 775* 776 IF( W*ETA.GE.ZERO ) 777 $ ETA = -W / DW 778 TEMP = TAU + ETA 779 IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN 780 IF( W.LT.ZERO ) THEN 781 ETA = ( DLTUB-TAU ) / TWO 782 ELSE 783 ETA = ( DLTLB-TAU ) / TWO 784 END IF 785 END IF 786* 787 DO 210 J = 1, N 788 DELTA( J ) = DELTA( J ) - ETA 789 210 CONTINUE 790* 791 TAU = TAU + ETA 792 PREW = W 793* 794* Evaluate PSI and the derivative DPSI 795* 796 DPSI = ZERO 797 PSI = ZERO 798 ERRETM = ZERO 799 DO 220 J = 1, IIM1 800 TEMP = Z( J ) / DELTA( J ) 801 PSI = PSI + Z( J )*TEMP 802 DPSI = DPSI + TEMP*TEMP 803 ERRETM = ERRETM + PSI 804 220 CONTINUE 805 ERRETM = ABS( ERRETM ) 806* 807* Evaluate PHI and the derivative DPHI 808* 809 DPHI = ZERO 810 PHI = ZERO 811 DO 230 J = N, IIP1, -1 812 TEMP = Z( J ) / DELTA( J ) 813 PHI = PHI + Z( J )*TEMP 814 DPHI = DPHI + TEMP*TEMP 815 ERRETM = ERRETM + PHI 816 230 CONTINUE 817* 818 TEMP = Z( II ) / DELTA( II ) 819 DW = DPSI + DPHI + TEMP*TEMP 820 TEMP = Z( II )*TEMP 821 W = RHOINV + PHI + PSI + TEMP 822 ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV + 823 $ THREE*ABS( TEMP ) + ABS( TAU )*DW 824 IF( W*PREW.GT.ZERO .AND. ABS( W ).GT.ABS( PREW ) / TEN ) 825 $ SWTCH = .NOT.SWTCH 826* 827 240 CONTINUE 828* 829* Return with INFO = 1, NITER = MAXIT and not converged 830* 831 INFO = 1 832 IF( ORGATI ) THEN 833 DLAM = D( I ) + TAU 834 ELSE 835 DLAM = D( IP1 ) + TAU 836 END IF 837* 838 END IF 839* 840 250 CONTINUE 841* 842 RETURN 843* 844* End of SLAED4 845* 846 END 847