SRC/slaed4.f

   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