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