LAPACK: SRC/slasd4.f Source File

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