LAPACK: SRC/dlaed6.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     February 2007
00007 *
00008 *     .. Scalar Arguments ..
00009       LOGICAL            ORGATI
00010       INTEGER            INFO, KNITER
00011       DOUBLE PRECISION   FINIT, RHO, TAU
00012 *     ..
00013 *     .. Array Arguments ..
00014       DOUBLE PRECISION   D( 3 ), Z( 3 )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  DLAED6 computes the positive or negative root (closest to the origin)
00021 *  of
00022 *                   z(1)        z(2)        z(3)
00023 *  f(x) =   rho + --------- + ---------- + ---------
00024 *                  d(1)-x      d(2)-x      d(3)-x
00025 *
00026 *  It is assumed that
00027 *
00028 *        if ORGATI = .true. the root is between d(2) and d(3);
00029 *        otherwise it is between d(1) and d(2)
00030 *
00031 *  This routine will be called by DLAED4 when necessary. In most cases,
00032 *  the root sought is the smallest in magnitude, though it might not be
00033 *  in some extremely rare situations.
00034 *
00035 *  Arguments
00036 *  =========
00037 *
00038 *  KNITER       (input) INTEGER
00039 *               Refer to DLAED4 for its significance.
00040 *
00041 *  ORGATI       (input) LOGICAL
00042 *               If ORGATI is true, the needed root is between d(2) and
00043 *               d(3); otherwise it is between d(1) and d(2).  See
00044 *               DLAED4 for further details.
00045 *
00046 *  RHO          (input) DOUBLE PRECISION
00047 *               Refer to the equation f(x) above.
00048 *
00049 *  D            (input) DOUBLE PRECISION array, dimension (3)
00050 *               D satisfies d(1) < d(2) < d(3).
00051 *
00052 *  Z            (input) DOUBLE PRECISION array, dimension (3)
00053 *               Each of the elements in z must be positive.
00054 *
00055 *  FINIT        (input) DOUBLE PRECISION
00056 *               The value of f at 0. It is more accurate than the one
00057 *               evaluated inside this routine (if someone wants to do
00058 *               so).
00059 *
00060 *  TAU          (output) DOUBLE PRECISION
00061 *               The root of the equation f(x).
00062 *
00063 *  INFO         (output) INTEGER
00064 *               = 0: successful exit
00065 *               > 0: if INFO = 1, failure to converge
00066 *
00067 *  Further Details
00068 *  ===============
00069 *
00070 *  30/06/99: Based on contributions by
00071 *     Ren-Cang Li, Computer Science Division, University of California
00072 *     at Berkeley, USA
00073 *
00074 *  10/02/03: This version has a few statements commented out for thread
00075 *  safety (machine parameters are computed on each entry). SJH.
00076 *
00077 *  05/10/06: Modified from a new version of Ren-Cang Li, use
00078 *     Gragg-Thornton-Warner cubic convergent scheme for better stability.
00079 *
00080 *  =====================================================================
00081 *
00082 *     .. Parameters ..
00083       INTEGER            MAXIT
00084       PARAMETER          ( MAXIT = 40 )
00085       DOUBLE PRECISION   ZERO, ONE, TWO, THREE, FOUR, EIGHT
00086       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
00087      $                   THREE = 3.0D0, FOUR = 4.0D0, EIGHT = 8.0D0 )
00088 *     ..
00089 *     .. External Functions ..
00090       DOUBLE PRECISION   DLAMCH
00091       EXTERNAL           DLAMCH
00092 *     ..
00093 *     .. Local Arrays ..
00094       DOUBLE PRECISION   DSCALE( 3 ), ZSCALE( 3 )
00095 *     ..
00096 *     .. Local Scalars ..
00097       LOGICAL            SCALE
00098       INTEGER            I, ITER, NITER
00099       DOUBLE PRECISION   A, B, BASE, C, DDF, DF, EPS, ERRETM, ETA, F,
00100      $                   FC, SCLFAC, SCLINV, SMALL1, SMALL2, SMINV1,
00101      $                   SMINV2, TEMP, TEMP1, TEMP2, TEMP3, TEMP4, 
00102      $                   LBD, UBD
00103 *     ..
00104 *     .. Intrinsic Functions ..
00105       INTRINSIC          ABS, INT, LOG, MAX, MIN, SQRT
00106 *     ..
00107 *     .. Executable Statements ..
00108 *
00109       INFO = 0
00110 *
00111       IF( ORGATI ) THEN
00112          LBD = D(2)
00113          UBD = D(3)
00114       ELSE
00115          LBD = D(1)
00116          UBD = D(2)
00117       END IF
00118       IF( FINIT .LT. ZERO )THEN
00119          LBD = ZERO
00120       ELSE
00121          UBD = ZERO 
00122       END IF
00123 *
00124       NITER = 1
00125       TAU = ZERO
00126       IF( KNITER.EQ.2 ) THEN
00127          IF( ORGATI ) THEN
00128             TEMP = ( D( 3 )-D( 2 ) ) / TWO
00129             C = RHO + Z( 1 ) / ( ( D( 1 )-D( 2 ) )-TEMP )
00130             A = C*( D( 2 )+D( 3 ) ) + Z( 2 ) + Z( 3 )
00131             B = C*D( 2 )*D( 3 ) + Z( 2 )*D( 3 ) + Z( 3 )*D( 2 )
00132          ELSE
00133             TEMP = ( D( 1 )-D( 2 ) ) / TWO
00134             C = RHO + Z( 3 ) / ( ( D( 3 )-D( 2 ) )-TEMP )
00135             A = C*( D( 1 )+D( 2 ) ) + Z( 1 ) + Z( 2 )
00136             B = C*D( 1 )*D( 2 ) + Z( 1 )*D( 2 ) + Z( 2 )*D( 1 )
00137          END IF
00138          TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) )
00139          A = A / TEMP
00140          B = B / TEMP
00141          C = C / TEMP
00142          IF( C.EQ.ZERO ) THEN
00143             TAU = B / A
00144          ELSE IF( A.LE.ZERO ) THEN
00145             TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
00146          ELSE
00147             TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
00148          END IF
00149          IF( TAU .LT. LBD .OR. TAU .GT. UBD )
00150      $      TAU = ( LBD+UBD )/TWO
00151          IF( D(1).EQ.TAU .OR. D(2).EQ.TAU .OR. D(3).EQ.TAU ) THEN
00152             TAU = ZERO
00153          ELSE
00154             TEMP = FINIT + TAU*Z(1)/( D(1)*( D( 1 )-TAU ) ) +
00155      $                     TAU*Z(2)/( D(2)*( D( 2 )-TAU ) ) +
00156      $                     TAU*Z(3)/( D(3)*( D( 3 )-TAU ) )
00157             IF( TEMP .LE. ZERO )THEN
00158                LBD = TAU
00159             ELSE
00160                UBD = TAU
00161             END IF
00162             IF( ABS( FINIT ).LE.ABS( TEMP ) )
00163      $         TAU = ZERO
00164          END IF
00165       END IF
00166 *
00167 *     get machine parameters for possible scaling to avoid overflow
00168 *
00169 *     modified by Sven: parameters SMALL1, SMINV1, SMALL2,
00170 *     SMINV2, EPS are not SAVEd anymore between one call to the
00171 *     others but recomputed at each call
00172 *
00173       EPS = DLAMCH( 'Epsilon' )
00174       BASE = DLAMCH( 'Base' )
00175       SMALL1 = BASE**( INT( LOG( DLAMCH( 'SafMin' ) ) / LOG( BASE ) /
00176      $         THREE ) )
00177       SMINV1 = ONE / SMALL1
00178       SMALL2 = SMALL1*SMALL1
00179       SMINV2 = SMINV1*SMINV1
00180 *
00181 *     Determine if scaling of inputs necessary to avoid overflow
00182 *     when computing 1/TEMP**3
00183 *
00184       IF( ORGATI ) THEN
00185          TEMP = MIN( ABS( D( 2 )-TAU ), ABS( D( 3 )-TAU ) )
00186       ELSE
00187          TEMP = MIN( ABS( D( 1 )-TAU ), ABS( D( 2 )-TAU ) )
00188       END IF
00189       SCALE = .FALSE.
00190       IF( TEMP.LE.SMALL1 ) THEN
00191          SCALE = .TRUE.
00192          IF( TEMP.LE.SMALL2 ) THEN
00193 *
00194 *        Scale up by power of radix nearest 1/SAFMIN**(2/3)
00195 *
00196             SCLFAC = SMINV2
00197             SCLINV = SMALL2
00198          ELSE
00199 *
00200 *        Scale up by power of radix nearest 1/SAFMIN**(1/3)
00201 *
00202             SCLFAC = SMINV1
00203             SCLINV = SMALL1
00204          END IF
00205 *
00206 *        Scaling up safe because D, Z, TAU scaled elsewhere to be O(1)
00207 *
00208          DO 10 I = 1, 3
00209             DSCALE( I ) = D( I )*SCLFAC
00210             ZSCALE( I ) = Z( I )*SCLFAC
00211    10    CONTINUE
00212          TAU = TAU*SCLFAC
00213          LBD = LBD*SCLFAC
00214          UBD = UBD*SCLFAC
00215       ELSE
00216 *
00217 *        Copy D and Z to DSCALE and ZSCALE
00218 *
00219          DO 20 I = 1, 3
00220             DSCALE( I ) = D( I )
00221             ZSCALE( I ) = Z( I )
00222    20    CONTINUE
00223       END IF
00224 *
00225       FC = ZERO
00226       DF = ZERO
00227       DDF = ZERO
00228       DO 30 I = 1, 3
00229          TEMP = ONE / ( DSCALE( I )-TAU )
00230          TEMP1 = ZSCALE( I )*TEMP
00231          TEMP2 = TEMP1*TEMP
00232          TEMP3 = TEMP2*TEMP
00233          FC = FC + TEMP1 / DSCALE( I )
00234          DF = DF + TEMP2
00235          DDF = DDF + TEMP3
00236    30 CONTINUE
00237       F = FINIT + TAU*FC
00238 *
00239       IF( ABS( F ).LE.ZERO )
00240      $   GO TO 60
00241       IF( F .LE. ZERO )THEN
00242          LBD = TAU
00243       ELSE
00244          UBD = TAU
00245       END IF
00246 *
00247 *        Iteration begins -- Use Gragg-Thornton-Warner cubic convergent
00248 *                            scheme
00249 *
00250 *     It is not hard to see that
00251 *
00252 *           1) Iterations will go up monotonically
00253 *              if FINIT < 0;
00254 *
00255 *           2) Iterations will go down monotonically
00256 *              if FINIT > 0.
00257 *
00258       ITER = NITER + 1
00259 *
00260       DO 50 NITER = ITER, MAXIT
00261 *
00262          IF( ORGATI ) THEN
00263             TEMP1 = DSCALE( 2 ) - TAU
00264             TEMP2 = DSCALE( 3 ) - TAU
00265          ELSE
00266             TEMP1 = DSCALE( 1 ) - TAU
00267             TEMP2 = DSCALE( 2 ) - TAU
00268          END IF
00269          A = ( TEMP1+TEMP2 )*F - TEMP1*TEMP2*DF
00270          B = TEMP1*TEMP2*F
00271          C = F - ( TEMP1+TEMP2 )*DF + TEMP1*TEMP2*DDF
00272          TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) )
00273          A = A / TEMP
00274          B = B / TEMP
00275          C = C / TEMP
00276          IF( C.EQ.ZERO ) THEN
00277             ETA = B / A
00278          ELSE IF( A.LE.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          IF( F*ETA.GE.ZERO ) THEN
00284             ETA = -F / DF
00285          END IF
00286 *
00287          TAU = TAU + ETA
00288          IF( TAU .LT. LBD .OR. TAU .GT. UBD )
00289      $      TAU = ( LBD + UBD )/TWO 
00290 *
00291          FC = ZERO
00292          ERRETM = ZERO
00293          DF = ZERO
00294          DDF = ZERO
00295          DO 40 I = 1, 3
00296             TEMP = ONE / ( DSCALE( I )-TAU )
00297             TEMP1 = ZSCALE( I )*TEMP
00298             TEMP2 = TEMP1*TEMP
00299             TEMP3 = TEMP2*TEMP
00300             TEMP4 = TEMP1 / DSCALE( I )
00301             FC = FC + TEMP4
00302             ERRETM = ERRETM + ABS( TEMP4 )
00303             DF = DF + TEMP2
00304             DDF = DDF + TEMP3
00305    40    CONTINUE
00306          F = FINIT + TAU*FC
00307          ERRETM = EIGHT*( ABS( FINIT )+ABS( TAU )*ERRETM ) +
00308      $            ABS( TAU )*DF
00309          IF( ABS( F ).LE.EPS*ERRETM )
00310      $      GO TO 60
00311          IF( F .LE. ZERO )THEN
00312             LBD = TAU
00313          ELSE
00314             UBD = TAU
00315          END IF
00316    50 CONTINUE
00317       INFO = 1
00318    60 CONTINUE
00319 *
00320 *     Undo scaling
00321 *
00322       IF( SCALE )
00323      $   TAU = TAU*SCLINV
00324       RETURN
00325 *
00326 *     End of DLAED6
00327 *
00328       END