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