LAPACK: SRC/dlaed4.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            I, INFO, N
00010       DOUBLE PRECISION   DLAM, RHO
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   D( * ), DELTA( * ), Z( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  This subroutine computes the I-th updated eigenvalue of a symmetric
00020 *  rank-one modification to a diagonal matrix whose elements are
00021 *  given in the array d, and that
00022 *
00023 *             D(i) < D(j)  for  i < j
00024 *
00025 *  and that RHO > 0.  This is arranged by the calling routine, and is
00026 *  no loss in generality.  The rank-one modified system is thus
00027 *
00028 *             diag( D )  +  RHO *  Z * Z_transpose.
00029 *
00030 *  where we assume the Euclidean norm of Z is 1.
00031 *
00032 *  The method consists of approximating the rational functions in the
00033 *  secular equation by simpler interpolating rational functions.
00034 *
00035 *  Arguments
00036 *  =========
00037 *
00038 *  N      (input) INTEGER
00039 *         The length of all arrays.
00040 *
00041 *  I      (input) INTEGER
00042 *         The index of the eigenvalue to be computed.  1 <= I <= N.
00043 *
00044 *  D      (input) DOUBLE PRECISION array, dimension (N)
00045 *         The original eigenvalues.  It is assumed that they are in
00046 *         order, D(I) < D(J)  for I < J.
00047 *
00048 *  Z      (input) DOUBLE PRECISION array, dimension (N)
00049 *         The components of the updating vector.
00050 *
00051 *  DELTA  (output) DOUBLE PRECISION array, dimension (N)
00052 *         If N .GT. 2, DELTA contains (D(j) - lambda_I) in its  j-th
00053 *         component.  If N = 1, then DELTA(1) = 1. If N = 2, see DLAED5
00054 *         for detail. The vector DELTA contains the information necessary
00055 *         to construct the eigenvectors by DLAED3 and DLAED9.
00056 *
00057 *  RHO    (input) DOUBLE PRECISION
00058 *         The scalar in the symmetric updating formula.
00059 *
00060 *  DLAM   (output) DOUBLE PRECISION
00061 *         The computed lambda_I, the I-th updated eigenvalue.
00062 *
00063 *  INFO   (output) INTEGER
00064 *         = 0:  successful exit
00065 *         > 0:  if INFO = 1, the updating process failed.
00066 *
00067 *  Internal Parameters
00068 *  ===================
00069 *
00070 *  Logical variable ORGATI (origin-at-i?) is used for distinguishing
00071 *  whether D(i) or D(i+1) is treated as the origin.
00072 *
00073 *            ORGATI = .true.    origin at i
00074 *            ORGATI = .false.   origin at i+1
00075 *
00076 *   Logical variable SWTCH3 (switch-for-3-poles?) is for noting
00077 *   if we are working with THREE poles!
00078 *
00079 *   MAXIT is the maximum number of iterations allowed for each
00080 *   eigenvalue.
00081 *
00082 *  Further Details
00083 *  ===============
00084 *
00085 *  Based on contributions by
00086 *     Ren-Cang Li, Computer Science Division, University of California
00087 *     at Berkeley, USA
00088 *
00089 *  =====================================================================
00090 *
00091 *     .. Parameters ..
00092       INTEGER            MAXIT
00093       PARAMETER          ( MAXIT = 30 )
00094       DOUBLE PRECISION   ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN
00095       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
00096      $                   THREE = 3.0D0, FOUR = 4.0D0, EIGHT = 8.0D0,
00097      $                   TEN = 10.0D0 )
00098 *     ..
00099 *     .. Local Scalars ..
00100       LOGICAL            ORGATI, SWTCH, SWTCH3
00101       INTEGER            II, IIM1, IIP1, IP1, ITER, J, NITER
00102       DOUBLE PRECISION   A, B, C, DEL, DLTLB, DLTUB, DPHI, DPSI, DW,
00103      $                   EPS, ERRETM, ETA, MIDPT, PHI, PREW, PSI,
00104      $                   RHOINV, TAU, TEMP, TEMP1, W
00105 *     ..
00106 *     .. Local Arrays ..
00107       DOUBLE PRECISION   ZZ( 3 )
00108 *     ..
00109 *     .. External Functions ..
00110       DOUBLE PRECISION   DLAMCH
00111       EXTERNAL           DLAMCH
00112 *     ..
00113 *     .. External Subroutines ..
00114       EXTERNAL           DLAED5, DLAED6
00115 *     ..
00116 *     .. Intrinsic Functions ..
00117       INTRINSIC          ABS, MAX, MIN, SQRT
00118 *     ..
00119 *     .. Executable Statements ..
00120 *
00121 *     Since this routine is called in an inner loop, we do no argument
00122 *     checking.
00123 *
00124 *     Quick return for N=1 and 2.
00125 *
00126       INFO = 0
00127       IF( N.EQ.1 ) THEN
00128 *
00129 *         Presumably, I=1 upon entry
00130 *
00131          DLAM = D( 1 ) + RHO*Z( 1 )*Z( 1 )
00132          DELTA( 1 ) = ONE
00133          RETURN
00134       END IF
00135       IF( N.EQ.2 ) THEN
00136          CALL DLAED5( I, D, Z, DELTA, RHO, DLAM )
00137          RETURN
00138       END IF
00139 *
00140 *     Compute machine epsilon
00141 *
00142       EPS = DLAMCH( 'Epsilon' )
00143       RHOINV = ONE / RHO
00144 *
00145 *     The case I = N
00146 *
00147       IF( I.EQ.N ) THEN
00148 *
00149 *        Initialize some basic variables
00150 *
00151          II = N - 1
00152          NITER = 1
00153 *
00154 *        Calculate initial guess
00155 *
00156          MIDPT = RHO / TWO
00157 *
00158 *        If ||Z||_2 is not one, then TEMP should be set to
00159 *        RHO * ||Z||_2^2 / TWO
00160 *
00161          DO 10 J = 1, N
00162             DELTA( J ) = ( D( J )-D( I ) ) - MIDPT
00163    10    CONTINUE
00164 *
00165          PSI = ZERO
00166          DO 20 J = 1, N - 2
00167             PSI = PSI + Z( J )*Z( J ) / DELTA( J )
00168    20    CONTINUE
00169 *
00170          C = RHOINV + PSI
00171          W = C + Z( II )*Z( II ) / DELTA( II ) +
00172      $       Z( N )*Z( N ) / DELTA( N )
00173 *
00174          IF( W.LE.ZERO ) THEN
00175             TEMP = Z( N-1 )*Z( N-1 ) / ( D( N )-D( N-1 )+RHO ) +
00176      $             Z( N )*Z( N ) / RHO
00177             IF( C.LE.TEMP ) THEN
00178                TAU = RHO
00179             ELSE
00180                DEL = D( N ) - D( N-1 )
00181                A = -C*DEL + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )
00182                B = Z( N )*Z( N )*DEL
00183                IF( A.LT.ZERO ) THEN
00184                   TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )
00185                ELSE
00186                   TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )
00187                END IF
00188             END IF
00189 *
00190 *           It can be proved that
00191 *               D(N)+RHO/2 <= LAMBDA(N) < D(N)+TAU <= D(N)+RHO
00192 *
00193             DLTLB = MIDPT
00194             DLTUB = RHO
00195          ELSE
00196             DEL = D( N ) - D( N-1 )
00197             A = -C*DEL + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )
00198             B = Z( N )*Z( N )*DEL
00199             IF( A.LT.ZERO ) THEN
00200                TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )
00201             ELSE
00202                TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )
00203             END IF
00204 *
00205 *           It can be proved that
00206 *               D(N) < D(N)+TAU < LAMBDA(N) < D(N)+RHO/2
00207 *
00208             DLTLB = ZERO
00209             DLTUB = MIDPT
00210          END IF
00211 *
00212          DO 30 J = 1, N
00213             DELTA( J ) = ( D( J )-D( I ) ) - TAU
00214    30    CONTINUE
00215 *
00216 *        Evaluate PSI and the derivative DPSI
00217 *
00218          DPSI = ZERO
00219          PSI = ZERO
00220          ERRETM = ZERO
00221          DO 40 J = 1, II
00222             TEMP = Z( J ) / DELTA( J )
00223             PSI = PSI + Z( J )*TEMP
00224             DPSI = DPSI + TEMP*TEMP
00225             ERRETM = ERRETM + PSI
00226    40    CONTINUE
00227          ERRETM = ABS( ERRETM )
00228 *
00229 *        Evaluate PHI and the derivative DPHI
00230 *
00231          TEMP = Z( N ) / DELTA( N )
00232          PHI = Z( N )*TEMP
00233          DPHI = TEMP*TEMP
00234          ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +
00235      $            ABS( TAU )*( DPSI+DPHI )
00236 *
00237          W = RHOINV + PHI + PSI
00238 *
00239 *        Test for convergence
00240 *
00241          IF( ABS( W ).LE.EPS*ERRETM ) THEN
00242             DLAM = D( I ) + TAU
00243             GO TO 250
00244          END IF
00245 *
00246          IF( W.LE.ZERO ) THEN
00247             DLTLB = MAX( DLTLB, TAU )
00248          ELSE
00249             DLTUB = MIN( DLTUB, TAU )
00250          END IF
00251 *
00252 *        Calculate the new step
00253 *
00254          NITER = NITER + 1
00255          C = W - DELTA( N-1 )*DPSI - DELTA( N )*DPHI
00256          A = ( DELTA( N-1 )+DELTA( N ) )*W -
00257      $       DELTA( N-1 )*DELTA( N )*( DPSI+DPHI )
00258          B = DELTA( N-1 )*DELTA( N )*W
00259          IF( C.LT.ZERO )
00260      $      C = ABS( C )
00261          IF( C.EQ.ZERO ) THEN
00262 *          ETA = B/A
00263 *           ETA = RHO - TAU
00264             ETA = DLTUB - TAU
00265          ELSE IF( A.GE.ZERO ) THEN
00266             ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
00267          ELSE
00268             ETA = TWO*B / ( A-SQRT( ABS( A*A-FOUR*B*C ) ) )
00269          END IF
00270 *
00271 *        Note, eta should be positive if w is negative, and
00272 *        eta should be negative otherwise. However,
00273 *        if for some reason caused by roundoff, eta*w > 0,
00274 *        we simply use one Newton step instead. This way
00275 *        will guarantee eta*w < 0.
00276 *
00277          IF( W*ETA.GT.ZERO )
00278      $      ETA = -W / ( DPSI+DPHI )
00279          TEMP = TAU + ETA
00280          IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
00281             IF( W.LT.ZERO ) THEN
00282                ETA = ( DLTUB-TAU ) / TWO
00283             ELSE
00284                ETA = ( DLTLB-TAU ) / TWO
00285             END IF
00286          END IF
00287          DO 50 J = 1, N
00288             DELTA( J ) = DELTA( J ) - ETA
00289    50    CONTINUE
00290 *
00291          TAU = TAU + ETA
00292 *
00293 *        Evaluate PSI and the derivative DPSI
00294 *
00295          DPSI = ZERO
00296          PSI = ZERO
00297          ERRETM = ZERO
00298          DO 60 J = 1, II
00299             TEMP = Z( J ) / DELTA( J )
00300             PSI = PSI + Z( J )*TEMP
00301             DPSI = DPSI + TEMP*TEMP
00302             ERRETM = ERRETM + PSI
00303    60    CONTINUE
00304          ERRETM = ABS( ERRETM )
00305 *
00306 *        Evaluate PHI and the derivative DPHI
00307 *
00308          TEMP = Z( N ) / DELTA( N )
00309          PHI = Z( N )*TEMP
00310          DPHI = TEMP*TEMP
00311          ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +
00312      $            ABS( TAU )*( DPSI+DPHI )
00313 *
00314          W = RHOINV + PHI + PSI
00315 *
00316 *        Main loop to update the values of the array   DELTA
00317 *
00318          ITER = NITER + 1
00319 *
00320          DO 90 NITER = ITER, MAXIT
00321 *
00322 *           Test for convergence
00323 *
00324             IF( ABS( W ).LE.EPS*ERRETM ) THEN
00325                DLAM = D( I ) + TAU
00326                GO TO 250
00327             END IF
00328 *
00329             IF( W.LE.ZERO ) THEN
00330                DLTLB = MAX( DLTLB, TAU )
00331             ELSE
00332                DLTUB = MIN( DLTUB, TAU )
00333             END IF
00334 *
00335 *           Calculate the new step
00336 *
00337             C = W - DELTA( N-1 )*DPSI - DELTA( N )*DPHI
00338             A = ( DELTA( N-1 )+DELTA( N ) )*W -
00339      $          DELTA( N-1 )*DELTA( N )*( DPSI+DPHI )
00340             B = DELTA( N-1 )*DELTA( N )*W
00341             IF( A.GE.ZERO ) THEN
00342                ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
00343             ELSE
00344                ETA = TWO*B / ( A-SQRT( ABS( A*A-FOUR*B*C ) ) )
00345             END IF
00346 *
00347 *           Note, eta should be positive if w is negative, and
00348 *           eta should be negative otherwise. However,
00349 *           if for some reason caused by roundoff, eta*w > 0,
00350 *           we simply use one Newton step instead. This way
00351 *           will guarantee eta*w < 0.
00352 *
00353             IF( W*ETA.GT.ZERO )
00354      $         ETA = -W / ( DPSI+DPHI )
00355             TEMP = TAU + ETA
00356             IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
00357                IF( W.LT.ZERO ) THEN
00358                   ETA = ( DLTUB-TAU ) / TWO
00359                ELSE
00360                   ETA = ( DLTLB-TAU ) / TWO
00361                END IF
00362             END IF
00363             DO 70 J = 1, N
00364                DELTA( J ) = DELTA( J ) - ETA
00365    70       CONTINUE
00366 *
00367             TAU = TAU + ETA
00368 *
00369 *           Evaluate PSI and the derivative DPSI
00370 *
00371             DPSI = ZERO
00372             PSI = ZERO
00373             ERRETM = ZERO
00374             DO 80 J = 1, II
00375                TEMP = Z( J ) / DELTA( J )
00376                PSI = PSI + Z( J )*TEMP
00377                DPSI = DPSI + TEMP*TEMP
00378                ERRETM = ERRETM + PSI
00379    80       CONTINUE
00380             ERRETM = ABS( ERRETM )
00381 *
00382 *           Evaluate PHI and the derivative DPHI
00383 *
00384             TEMP = Z( N ) / DELTA( N )
00385             PHI = Z( N )*TEMP
00386             DPHI = TEMP*TEMP
00387             ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +
00388      $               ABS( TAU )*( DPSI+DPHI )
00389 *
00390             W = RHOINV + PHI + PSI
00391    90    CONTINUE
00392 *
00393 *        Return with INFO = 1, NITER = MAXIT and not converged
00394 *
00395          INFO = 1
00396          DLAM = D( I ) + TAU
00397          GO TO 250
00398 *
00399 *        End for the case I = N
00400 *
00401       ELSE
00402 *
00403 *        The case for I < N
00404 *
00405          NITER = 1
00406          IP1 = I + 1
00407 *
00408 *        Calculate initial guess
00409 *
00410          DEL = D( IP1 ) - D( I )
00411          MIDPT = DEL / TWO
00412          DO 100 J = 1, N
00413             DELTA( J ) = ( D( J )-D( I ) ) - MIDPT
00414   100    CONTINUE
00415 *
00416          PSI = ZERO
00417          DO 110 J = 1, I - 1
00418             PSI = PSI + Z( J )*Z( J ) / DELTA( J )
00419   110    CONTINUE
00420 *
00421          PHI = ZERO
00422          DO 120 J = N, I + 2, -1
00423             PHI = PHI + Z( J )*Z( J ) / DELTA( J )
00424   120    CONTINUE
00425          C = RHOINV + PSI + PHI
00426          W = C + Z( I )*Z( I ) / DELTA( I ) +
00427      $       Z( IP1 )*Z( IP1 ) / DELTA( IP1 )
00428 *
00429          IF( W.GT.ZERO ) THEN
00430 *
00431 *           d(i)< the ith eigenvalue < (d(i)+d(i+1))/2
00432 *
00433 *           We choose d(i) as origin.
00434 *
00435             ORGATI = .TRUE.
00436             A = C*DEL + Z( I )*Z( I ) + Z( IP1 )*Z( IP1 )
00437             B = Z( I )*Z( I )*DEL
00438             IF( A.GT.ZERO ) THEN
00439                TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
00440             ELSE
00441                TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
00442             END IF
00443             DLTLB = ZERO
00444             DLTUB = MIDPT
00445          ELSE
00446 *
00447 *           (d(i)+d(i+1))/2 <= the ith eigenvalue < d(i+1)
00448 *
00449 *           We choose d(i+1) as origin.
00450 *
00451             ORGATI = .FALSE.
00452             A = C*DEL - Z( I )*Z( I ) - Z( IP1 )*Z( IP1 )
00453             B = Z( IP1 )*Z( IP1 )*DEL
00454             IF( A.LT.ZERO ) THEN
00455                TAU = TWO*B / ( A-SQRT( ABS( A*A+FOUR*B*C ) ) )
00456             ELSE
00457                TAU = -( A+SQRT( ABS( A*A+FOUR*B*C ) ) ) / ( TWO*C )
00458             END IF
00459             DLTLB = -MIDPT
00460             DLTUB = ZERO
00461          END IF
00462 *
00463          IF( ORGATI ) THEN
00464             DO 130 J = 1, N
00465                DELTA( J ) = ( D( J )-D( I ) ) - TAU
00466   130       CONTINUE
00467          ELSE
00468             DO 140 J = 1, N
00469                DELTA( J ) = ( D( J )-D( IP1 ) ) - TAU
00470   140       CONTINUE
00471          END IF
00472          IF( ORGATI ) THEN
00473             II = I
00474          ELSE
00475             II = I + 1
00476          END IF
00477          IIM1 = II - 1
00478          IIP1 = II + 1
00479 *
00480 *        Evaluate PSI and the derivative DPSI
00481 *
00482          DPSI = ZERO
00483          PSI = ZERO
00484          ERRETM = ZERO
00485          DO 150 J = 1, IIM1
00486             TEMP = Z( J ) / DELTA( J )
00487             PSI = PSI + Z( J )*TEMP
00488             DPSI = DPSI + TEMP*TEMP
00489             ERRETM = ERRETM + PSI
00490   150    CONTINUE
00491          ERRETM = ABS( ERRETM )
00492 *
00493 *        Evaluate PHI and the derivative DPHI
00494 *
00495          DPHI = ZERO
00496          PHI = ZERO
00497          DO 160 J = N, IIP1, -1
00498             TEMP = Z( J ) / DELTA( J )
00499             PHI = PHI + Z( J )*TEMP
00500             DPHI = DPHI + TEMP*TEMP
00501             ERRETM = ERRETM + PHI
00502   160    CONTINUE
00503 *
00504          W = RHOINV + PHI + PSI
00505 *
00506 *        W is the value of the secular function with
00507 *        its ii-th element removed.
00508 *
00509          SWTCH3 = .FALSE.
00510          IF( ORGATI ) THEN
00511             IF( W.LT.ZERO )
00512      $         SWTCH3 = .TRUE.
00513          ELSE
00514             IF( W.GT.ZERO )
00515      $         SWTCH3 = .TRUE.
00516          END IF
00517          IF( II.EQ.1 .OR. II.EQ.N )
00518      $      SWTCH3 = .FALSE.
00519 *
00520          TEMP = Z( II ) / DELTA( II )
00521          DW = DPSI + DPHI + TEMP*TEMP
00522          TEMP = Z( II )*TEMP
00523          W = W + TEMP
00524          ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV +
00525      $            THREE*ABS( TEMP ) + ABS( TAU )*DW
00526 *
00527 *        Test for convergence
00528 *
00529          IF( ABS( W ).LE.EPS*ERRETM ) THEN
00530             IF( ORGATI ) THEN
00531                DLAM = D( I ) + TAU
00532             ELSE
00533                DLAM = D( IP1 ) + TAU
00534             END IF
00535             GO TO 250
00536          END IF
00537 *
00538          IF( W.LE.ZERO ) THEN
00539             DLTLB = MAX( DLTLB, TAU )
00540          ELSE
00541             DLTUB = MIN( DLTUB, TAU )
00542          END IF
00543 *
00544 *        Calculate the new step
00545 *
00546          NITER = NITER + 1
00547          IF( .NOT.SWTCH3 ) THEN
00548             IF( ORGATI ) THEN
00549                C = W - DELTA( IP1 )*DW - ( D( I )-D( IP1 ) )*
00550      $             ( Z( I ) / DELTA( I ) )**2
00551             ELSE
00552                C = W - DELTA( I )*DW - ( D( IP1 )-D( I ) )*
00553      $             ( Z( IP1 ) / DELTA( IP1 ) )**2
00554             END IF
00555             A = ( DELTA( I )+DELTA( IP1 ) )*W -
00556      $          DELTA( I )*DELTA( IP1 )*DW
00557             B = DELTA( I )*DELTA( IP1 )*W
00558             IF( C.EQ.ZERO ) THEN
00559                IF( A.EQ.ZERO ) THEN
00560                   IF( ORGATI ) THEN
00561                      A = Z( I )*Z( I ) + DELTA( IP1 )*DELTA( IP1 )*
00562      $                   ( DPSI+DPHI )
00563                   ELSE
00564                      A = Z( IP1 )*Z( IP1 ) + DELTA( I )*DELTA( I )*
00565      $                   ( DPSI+DPHI )
00566                   END IF
00567                END IF
00568                ETA = B / A
00569             ELSE IF( A.LE.ZERO ) THEN
00570                ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
00571             ELSE
00572                ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
00573             END IF
00574          ELSE
00575 *
00576 *           Interpolation using THREE most relevant poles
00577 *
00578             TEMP = RHOINV + PSI + PHI
00579             IF( ORGATI ) THEN
00580                TEMP1 = Z( IIM1 ) / DELTA( IIM1 )
00581                TEMP1 = TEMP1*TEMP1
00582                C = TEMP - DELTA( IIP1 )*( DPSI+DPHI ) -
00583      $             ( D( IIM1 )-D( IIP1 ) )*TEMP1
00584                ZZ( 1 ) = Z( IIM1 )*Z( IIM1 )
00585                ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )*
00586      $                   ( ( DPSI-TEMP1 )+DPHI )
00587             ELSE
00588                TEMP1 = Z( IIP1 ) / DELTA( IIP1 )
00589                TEMP1 = TEMP1*TEMP1
00590                C = TEMP - DELTA( IIM1 )*( DPSI+DPHI ) -
00591      $             ( D( IIP1 )-D( IIM1 ) )*TEMP1
00592                ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )*
00593      $                   ( DPSI+( DPHI-TEMP1 ) )
00594                ZZ( 3 ) = Z( IIP1 )*Z( IIP1 )
00595             END IF
00596             ZZ( 2 ) = Z( II )*Z( II )
00597             CALL DLAED6( NITER, ORGATI, C, DELTA( IIM1 ), ZZ, W, ETA,
00598      $                   INFO )
00599             IF( INFO.NE.0 )
00600      $         GO TO 250
00601          END IF
00602 *
00603 *        Note, eta should be positive if w is negative, and
00604 *        eta should be negative otherwise. However,
00605 *        if for some reason caused by roundoff, eta*w > 0,
00606 *        we simply use one Newton step instead. This way
00607 *        will guarantee eta*w < 0.
00608 *
00609          IF( W*ETA.GE.ZERO )
00610      $      ETA = -W / DW
00611          TEMP = TAU + ETA
00612          IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
00613             IF( W.LT.ZERO ) THEN
00614                ETA = ( DLTUB-TAU ) / TWO
00615             ELSE
00616                ETA = ( DLTLB-TAU ) / TWO
00617             END IF
00618          END IF
00619 *
00620          PREW = W
00621 *
00622          DO 180 J = 1, N
00623             DELTA( J ) = DELTA( J ) - ETA
00624   180    CONTINUE
00625 *
00626 *        Evaluate PSI and the derivative DPSI
00627 *
00628          DPSI = ZERO
00629          PSI = ZERO
00630          ERRETM = ZERO
00631          DO 190 J = 1, IIM1
00632             TEMP = Z( J ) / DELTA( J )
00633             PSI = PSI + Z( J )*TEMP
00634             DPSI = DPSI + TEMP*TEMP
00635             ERRETM = ERRETM + PSI
00636   190    CONTINUE
00637          ERRETM = ABS( ERRETM )
00638 *
00639 *        Evaluate PHI and the derivative DPHI
00640 *
00641          DPHI = ZERO
00642          PHI = ZERO
00643          DO 200 J = N, IIP1, -1
00644             TEMP = Z( J ) / DELTA( J )
00645             PHI = PHI + Z( J )*TEMP
00646             DPHI = DPHI + TEMP*TEMP
00647             ERRETM = ERRETM + PHI
00648   200    CONTINUE
00649 *
00650          TEMP = Z( II ) / DELTA( II )
00651          DW = DPSI + DPHI + TEMP*TEMP
00652          TEMP = Z( II )*TEMP
00653          W = RHOINV + PHI + PSI + TEMP
00654          ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV +
00655      $            THREE*ABS( TEMP ) + ABS( TAU+ETA )*DW
00656 *
00657          SWTCH = .FALSE.
00658          IF( ORGATI ) THEN
00659             IF( -W.GT.ABS( PREW ) / TEN )
00660      $         SWTCH = .TRUE.
00661          ELSE
00662             IF( W.GT.ABS( PREW ) / TEN )
00663      $         SWTCH = .TRUE.
00664          END IF
00665 *
00666          TAU = TAU + ETA
00667 *
00668 *        Main loop to update the values of the array   DELTA
00669 *
00670          ITER = NITER + 1
00671 *
00672          DO 240 NITER = ITER, MAXIT
00673 *
00674 *           Test for convergence
00675 *
00676             IF( ABS( W ).LE.EPS*ERRETM ) THEN
00677                IF( ORGATI ) THEN
00678                   DLAM = D( I ) + TAU
00679                ELSE
00680                   DLAM = D( IP1 ) + TAU
00681                END IF
00682                GO TO 250
00683             END IF
00684 *
00685             IF( W.LE.ZERO ) THEN
00686                DLTLB = MAX( DLTLB, TAU )
00687             ELSE
00688                DLTUB = MIN( DLTUB, TAU )
00689             END IF
00690 *
00691 *           Calculate the new step
00692 *
00693             IF( .NOT.SWTCH3 ) THEN
00694                IF( .NOT.SWTCH ) THEN
00695                   IF( ORGATI ) THEN
00696                      C = W - DELTA( IP1 )*DW -
00697      $                   ( D( I )-D( IP1 ) )*( Z( I ) / DELTA( I ) )**2
00698                   ELSE
00699                      C = W - DELTA( I )*DW - ( D( IP1 )-D( I ) )*
00700      $                   ( Z( IP1 ) / DELTA( IP1 ) )**2
00701                   END IF
00702                ELSE
00703                   TEMP = Z( II ) / DELTA( II )
00704                   IF( ORGATI ) THEN
00705                      DPSI = DPSI + TEMP*TEMP
00706                   ELSE
00707                      DPHI = DPHI + TEMP*TEMP
00708                   END IF
00709                   C = W - DELTA( I )*DPSI - DELTA( IP1 )*DPHI
00710                END IF
00711                A = ( DELTA( I )+DELTA( IP1 ) )*W -
00712      $             DELTA( I )*DELTA( IP1 )*DW
00713                B = DELTA( I )*DELTA( IP1 )*W
00714                IF( C.EQ.ZERO ) THEN
00715                   IF( A.EQ.ZERO ) THEN
00716                      IF( .NOT.SWTCH ) THEN
00717                         IF( ORGATI ) THEN
00718                            A = Z( I )*Z( I ) + DELTA( IP1 )*
00719      $                         DELTA( IP1 )*( DPSI+DPHI )
00720                         ELSE
00721                            A = Z( IP1 )*Z( IP1 ) +
00722      $                         DELTA( I )*DELTA( I )*( DPSI+DPHI )
00723                         END IF
00724                      ELSE
00725                         A = DELTA( I )*DELTA( I )*DPSI +
00726      $                      DELTA( IP1 )*DELTA( IP1 )*DPHI
00727                      END IF
00728                   END IF
00729                   ETA = B / A
00730                ELSE IF( A.LE.ZERO ) THEN
00731                   ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
00732                ELSE
00733                   ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
00734                END IF
00735             ELSE
00736 *
00737 *              Interpolation using THREE most relevant poles
00738 *
00739                TEMP = RHOINV + PSI + PHI
00740                IF( SWTCH ) THEN
00741                   C = TEMP - DELTA( IIM1 )*DPSI - DELTA( IIP1 )*DPHI
00742                   ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )*DPSI
00743                   ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )*DPHI
00744                ELSE
00745                   IF( ORGATI ) THEN
00746                      TEMP1 = Z( IIM1 ) / DELTA( IIM1 )
00747                      TEMP1 = TEMP1*TEMP1
00748                      C = TEMP - DELTA( IIP1 )*( DPSI+DPHI ) -
00749      $                   ( D( IIM1 )-D( IIP1 ) )*TEMP1
00750                      ZZ( 1 ) = Z( IIM1 )*Z( IIM1 )
00751                      ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )*
00752      $                         ( ( DPSI-TEMP1 )+DPHI )
00753                   ELSE
00754                      TEMP1 = Z( IIP1 ) / DELTA( IIP1 )
00755                      TEMP1 = TEMP1*TEMP1
00756                      C = TEMP - DELTA( IIM1 )*( DPSI+DPHI ) -
00757      $                   ( D( IIP1 )-D( IIM1 ) )*TEMP1
00758                      ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )*
00759      $                         ( DPSI+( DPHI-TEMP1 ) )
00760                      ZZ( 3 ) = Z( IIP1 )*Z( IIP1 )
00761                   END IF
00762                END IF
00763                CALL DLAED6( NITER, ORGATI, C, DELTA( IIM1 ), ZZ, W, ETA,
00764      $                      INFO )
00765                IF( INFO.NE.0 )
00766      $            GO TO 250
00767             END IF
00768 *
00769 *           Note, eta should be positive if w is negative, and
00770 *           eta should be negative otherwise. However,
00771 *           if for some reason caused by roundoff, eta*w > 0,
00772 *           we simply use one Newton step instead. This way
00773 *           will guarantee eta*w < 0.
00774 *
00775             IF( W*ETA.GE.ZERO )
00776      $         ETA = -W / DW
00777             TEMP = TAU + ETA
00778             IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
00779                IF( W.LT.ZERO ) THEN
00780                   ETA = ( DLTUB-TAU ) / TWO
00781                ELSE
00782                   ETA = ( DLTLB-TAU ) / TWO
00783                END IF
00784             END IF
00785 *
00786             DO 210 J = 1, N
00787                DELTA( J ) = DELTA( J ) - ETA
00788   210       CONTINUE
00789 *
00790             TAU = TAU + ETA
00791             PREW = W
00792 *
00793 *           Evaluate PSI and the derivative DPSI
00794 *
00795             DPSI = ZERO
00796             PSI = ZERO
00797             ERRETM = ZERO
00798             DO 220 J = 1, IIM1
00799                TEMP = Z( J ) / DELTA( J )
00800                PSI = PSI + Z( J )*TEMP
00801                DPSI = DPSI + TEMP*TEMP
00802                ERRETM = ERRETM + PSI
00803   220       CONTINUE
00804             ERRETM = ABS( ERRETM )
00805 *
00806 *           Evaluate PHI and the derivative DPHI
00807 *
00808             DPHI = ZERO
00809             PHI = ZERO
00810             DO 230 J = N, IIP1, -1
00811                TEMP = Z( J ) / DELTA( J )
00812                PHI = PHI + Z( J )*TEMP
00813                DPHI = DPHI + TEMP*TEMP
00814                ERRETM = ERRETM + PHI
00815   230       CONTINUE
00816 *
00817             TEMP = Z( II ) / DELTA( II )
00818             DW = DPSI + DPHI + TEMP*TEMP
00819             TEMP = Z( II )*TEMP
00820             W = RHOINV + PHI + PSI + TEMP
00821             ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV +
00822      $               THREE*ABS( TEMP ) + ABS( TAU )*DW
00823             IF( W*PREW.GT.ZERO .AND. ABS( W ).GT.ABS( PREW ) / TEN )
00824      $         SWTCH = .NOT.SWTCH
00825 *
00826   240    CONTINUE
00827 *
00828 *        Return with INFO = 1, NITER = MAXIT and not converged
00829 *
00830          INFO = 1
00831          IF( ORGATI ) THEN
00832             DLAM = D( I ) + TAU
00833          ELSE
00834             DLAM = D( IP1 ) + TAU
00835          END IF
00836 *
00837       END IF
00838 *
00839   250 CONTINUE
00840 *
00841       RETURN
00842 *
00843 *     End of DLAED4
00844 *
00845       END