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