Return to dlaed6.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence.

    1:       SUBROUTINE DLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO )
    2: *
    3: *  -- LAPACK routine (version 3.2) --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *     February 2007
    7: *
    8: *     .. Scalar Arguments ..
    9:       LOGICAL            ORGATI
   10:       INTEGER            INFO, KNITER
   11:       DOUBLE PRECISION   FINIT, RHO, TAU
   12: *     ..
   13: *     .. Array Arguments ..
   14:       DOUBLE PRECISION   D( 3 ), Z( 3 )
   15: *     ..
   16: *
   17: *  Purpose
   18: *  =======
   19: *
   20: *  DLAED6 computes the positive or negative root (closest to the origin)
   21: *  of
   22: *                   z(1)        z(2)        z(3)
   23: *  f(x) =   rho + --------- + ---------- + ---------
   24: *                  d(1)-x      d(2)-x      d(3)-x
   25: *
   26: *  It is assumed that
   27: *
   28: *        if ORGATI = .true. the root is between d(2) and d(3);
   29: *        otherwise it is between d(1) and d(2)
   30: *
   31: *  This routine will be called by DLAED4 when necessary. In most cases,
   32: *  the root sought is the smallest in magnitude, though it might not be
   33: *  in some extremely rare situations.
   34: *
   35: *  Arguments
   36: *  =========
   37: *
   38: *  KNITER       (input) INTEGER
   39: *               Refer to DLAED4 for its significance.
   40: *
   41: *  ORGATI       (input) LOGICAL
   42: *               If ORGATI is true, the needed root is between d(2) and
   43: *               d(3); otherwise it is between d(1) and d(2).  See
   44: *               DLAED4 for further details.
   45: *
   46: *  RHO          (input) DOUBLE PRECISION
   47: *               Refer to the equation f(x) above.
   48: *
   49: *  D            (input) DOUBLE PRECISION array, dimension (3)
   50: *               D satisfies d(1) < d(2) < d(3).
   51: *
   52: *  Z            (input) DOUBLE PRECISION array, dimension (3)
   53: *               Each of the elements in z must be positive.
   54: *
   55: *  FINIT        (input) DOUBLE PRECISION
   56: *               The value of f at 0. It is more accurate than the one
   57: *               evaluated inside this routine (if someone wants to do
   58: *               so).
   59: *
   60: *  TAU          (output) DOUBLE PRECISION
   61: *               The root of the equation f(x).
   62: *
   63: *  INFO         (output) INTEGER
   64: *               = 0: successful exit
   65: *               > 0: if INFO = 1, failure to converge
   66: *
   67: *  Further Details
   68: *  ===============
   69: *
   70: *  30/06/99: Based on contributions by
   71: *     Ren-Cang Li, Computer Science Division, University of California
   72: *     at Berkeley, USA
   73: *
   74: *  10/02/03: This version has a few statements commented out for thread
   75: *  safety (machine parameters are computed on each entry). SJH.
   76: *
   77: *  05/10/06: Modified from a new version of Ren-Cang Li, use
   78: *     Gragg-Thornton-Warner cubic convergent scheme for better stability.
   79: *
   80: *  =====================================================================
   81: *
   82: *     .. Parameters ..
   83:       INTEGER            MAXIT
   84:       PARAMETER          ( MAXIT = 40 )
   85:       DOUBLE PRECISION   ZERO, ONE, TWO, THREE, FOUR, EIGHT
   86:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
   87:      $                   THREE = 3.0D0, FOUR = 4.0D0, EIGHT = 8.0D0 )
   88: *     ..
   89: *     .. External Functions ..
   90:       DOUBLE PRECISION   DLAMCH
   91:       EXTERNAL           DLAMCH
   92: *     ..
   93: *     .. Local Arrays ..
   94:       DOUBLE PRECISION   DSCALE( 3 ), ZSCALE( 3 )
   95: *     ..
   96: *     .. Local Scalars ..
   97:       LOGICAL            SCALE
   98:       INTEGER            I, ITER, NITER
   99:       DOUBLE PRECISION   A, B, BASE, C, DDF, DF, EPS, ERRETM, ETA, F,
  100:      $                   FC, SCLFAC, SCLINV, SMALL1, SMALL2, SMINV1,
  101:      $                   SMINV2, TEMP, TEMP1, TEMP2, TEMP3, TEMP4, 
  102:      $                   LBD, UBD
  103: *     ..
  104: *     .. Intrinsic Functions ..
  105:       INTRINSIC          ABS, INT, LOG, MAX, MIN, SQRT
  106: *     ..
  107: *     .. Executable Statements ..
  108: *
  109:       INFO = 0
  110: *
  111:       IF( ORGATI ) THEN
  112:          LBD = D(2)
  113:          UBD = D(3)
  114:       ELSE
  115:          LBD = D(1)
  116:          UBD = D(2)
  117:       END IF
  118:       IF( FINIT .LT. ZERO )THEN
  119:          LBD = ZERO
  120:       ELSE
  121:          UBD = ZERO 
  122:       END IF
  123: *
  124:       NITER = 1
  125:       TAU = ZERO
  126:       IF( KNITER.EQ.2 ) THEN
  127:          IF( ORGATI ) THEN
  128:             TEMP = ( D( 3 )-D( 2 ) ) / TWO
  129:             C = RHO + Z( 1 ) / ( ( D( 1 )-D( 2 ) )-TEMP )
  130:             A = C*( D( 2 )+D( 3 ) ) + Z( 2 ) + Z( 3 )
  131:             B = C*D( 2 )*D( 3 ) + Z( 2 )*D( 3 ) + Z( 3 )*D( 2 )
  132:          ELSE
  133:             TEMP = ( D( 1 )-D( 2 ) ) / TWO
  134:             C = RHO + Z( 3 ) / ( ( D( 3 )-D( 2 ) )-TEMP )
  135:             A = C*( D( 1 )+D( 2 ) ) + Z( 1 ) + Z( 2 )
  136:             B = C*D( 1 )*D( 2 ) + Z( 1 )*D( 2 ) + Z( 2 )*D( 1 )
  137:          END IF
  138:          TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) )
  139:          A = A / TEMP
  140:          B = B / TEMP
  141:          C = C / TEMP
  142:          IF( C.EQ.ZERO ) THEN
  143:             TAU = B / A
  144:          ELSE IF( A.LE.ZERO ) THEN
  145:             TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
  146:          ELSE
  147:             TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
  148:          END IF
  149:          IF( TAU .LT. LBD .OR. TAU .GT. UBD )
  150:      $      TAU = ( LBD+UBD )/TWO
  151:          IF( D(1).EQ.TAU .OR. D(2).EQ.TAU .OR. D(3).EQ.TAU ) THEN
  152:             TAU = ZERO
  153:          ELSE
  154:             TEMP = FINIT + TAU*Z(1)/( D(1)*( D( 1 )-TAU ) ) +
  155:      $                     TAU*Z(2)/( D(2)*( D( 2 )-TAU ) ) +
  156:      $                     TAU*Z(3)/( D(3)*( D( 3 )-TAU ) )
  157:             IF( TEMP .LE. ZERO )THEN
  158:                LBD = TAU
  159:             ELSE
  160:                UBD = TAU
  161:             END IF
  162:             IF( ABS( FINIT ).LE.ABS( TEMP ) )
  163:      $         TAU = ZERO
  164:          END IF
  165:       END IF
  166: *
  167: *     get machine parameters for possible scaling to avoid overflow
  168: *
  169: *     modified by Sven: parameters SMALL1, SMINV1, SMALL2,
  170: *     SMINV2, EPS are not SAVEd anymore between one call to the
  171: *     others but recomputed at each call
  172: *
  173:       EPS = DLAMCH( 'Epsilon' )
  174:       BASE = DLAMCH( 'Base' )
  175:       SMALL1 = BASE**( INT( LOG( DLAMCH( 'SafMin' ) ) / LOG( BASE ) /
  176:      $         THREE ) )
  177:       SMINV1 = ONE / SMALL1
  178:       SMALL2 = SMALL1*SMALL1
  179:       SMINV2 = SMINV1*SMINV1
  180: *
  181: *     Determine if scaling of inputs necessary to avoid overflow
  182: *     when computing 1/TEMP**3
  183: *
  184:       IF( ORGATI ) THEN
  185:          TEMP = MIN( ABS( D( 2 )-TAU ), ABS( D( 3 )-TAU ) )
  186:       ELSE
  187:          TEMP = MIN( ABS( D( 1 )-TAU ), ABS( D( 2 )-TAU ) )
  188:       END IF
  189:       SCALE = .FALSE.
  190:       IF( TEMP.LE.SMALL1 ) THEN
  191:          SCALE = .TRUE.
  192:          IF( TEMP.LE.SMALL2 ) THEN
  193: *
  194: *        Scale up by power of radix nearest 1/SAFMIN**(2/3)
  195: *
  196:             SCLFAC = SMINV2
  197:             SCLINV = SMALL2
  198:          ELSE
  199: *
  200: *        Scale up by power of radix nearest 1/SAFMIN**(1/3)
  201: *
  202:             SCLFAC = SMINV1
  203:             SCLINV = SMALL1
  204:          END IF
  205: *
  206: *        Scaling up safe because D, Z, TAU scaled elsewhere to be O(1)
  207: *
  208:          DO 10 I = 1, 3
  209:             DSCALE( I ) = D( I )*SCLFAC
  210:             ZSCALE( I ) = Z( I )*SCLFAC
  211:    10    CONTINUE
  212:          TAU = TAU*SCLFAC
  213:          LBD = LBD*SCLFAC
  214:          UBD = UBD*SCLFAC
  215:       ELSE
  216: *
  217: *        Copy D and Z to DSCALE and ZSCALE
  218: *
  219:          DO 20 I = 1, 3
  220:             DSCALE( I ) = D( I )
  221:             ZSCALE( I ) = Z( I )
  222:    20    CONTINUE
  223:       END IF
  224: *
  225:       FC = ZERO
  226:       DF = ZERO
  227:       DDF = ZERO
  228:       DO 30 I = 1, 3
  229:          TEMP = ONE / ( DSCALE( I )-TAU )
  230:          TEMP1 = ZSCALE( I )*TEMP
  231:          TEMP2 = TEMP1*TEMP
  232:          TEMP3 = TEMP2*TEMP
  233:          FC = FC + TEMP1 / DSCALE( I )
  234:          DF = DF + TEMP2
  235:          DDF = DDF + TEMP3
  236:    30 CONTINUE
  237:       F = FINIT + TAU*FC
  238: *
  239:       IF( ABS( F ).LE.ZERO )
  240:      $   GO TO 60
  241:       IF( F .LE. ZERO )THEN
  242:          LBD = TAU
  243:       ELSE
  244:          UBD = TAU
  245:       END IF
  246: *
  247: *        Iteration begins -- Use Gragg-Thornton-Warner cubic convergent
  248: *                            scheme
  249: *
  250: *     It is not hard to see that
  251: *
  252: *           1) Iterations will go up monotonically
  253: *              if FINIT < 0;
  254: *
  255: *           2) Iterations will go down monotonically
  256: *              if FINIT > 0.
  257: *
  258:       ITER = NITER + 1
  259: *
  260:       DO 50 NITER = ITER, MAXIT
  261: *
  262:          IF( ORGATI ) THEN
  263:             TEMP1 = DSCALE( 2 ) - TAU
  264:             TEMP2 = DSCALE( 3 ) - TAU
  265:          ELSE
  266:             TEMP1 = DSCALE( 1 ) - TAU
  267:             TEMP2 = DSCALE( 2 ) - TAU
  268:          END IF
  269:          A = ( TEMP1+TEMP2 )*F - TEMP1*TEMP2*DF
  270:          B = TEMP1*TEMP2*F
  271:          C = F - ( TEMP1+TEMP2 )*DF + TEMP1*TEMP2*DDF
  272:          TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) )
  273:          A = A / TEMP
  274:          B = B / TEMP
  275:          C = C / TEMP
  276:          IF( C.EQ.ZERO ) THEN
  277:             ETA = B / A
  278:          ELSE IF( A.LE.ZERO ) THEN
  279:             ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
  280:          ELSE
  281:             ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
  282:          END IF
  283:          IF( F*ETA.GE.ZERO ) THEN
  284:             ETA = -F / DF
  285:          END IF
  286: *
  287:          TAU = TAU + ETA
  288:          IF( TAU .LT. LBD .OR. TAU .GT. UBD )
  289:      $      TAU = ( LBD + UBD )/TWO 
  290: *
  291:          FC = ZERO
  292:          ERRETM = ZERO
  293:          DF = ZERO
  294:          DDF = ZERO
  295:          DO 40 I = 1, 3
  296:             TEMP = ONE / ( DSCALE( I )-TAU )
  297:             TEMP1 = ZSCALE( I )*TEMP
  298:             TEMP2 = TEMP1*TEMP
  299:             TEMP3 = TEMP2*TEMP
  300:             TEMP4 = TEMP1 / DSCALE( I )
  301:             FC = FC + TEMP4
  302:             ERRETM = ERRETM + ABS( TEMP4 )
  303:             DF = DF + TEMP2
  304:             DDF = DDF + TEMP3
  305:    40    CONTINUE
  306:          F = FINIT + TAU*FC
  307:          ERRETM = EIGHT*( ABS( FINIT )+ABS( TAU )*ERRETM ) +
  308:      $            ABS( TAU )*DF
  309:          IF( ABS( F ).LE.EPS*ERRETM )
  310:      $      GO TO 60
  311:          IF( F .LE. ZERO )THEN
  312:             LBD = TAU
  313:          ELSE
  314:             UBD = TAU
  315:          END IF
  316:    50 CONTINUE
  317:       INFO = 1
  318:    60 CONTINUE
  319: *
  320: *     Undo scaling
  321: *
  322:       IF( SCALE )
  323:      $   TAU = TAU*SCLINV
  324:       RETURN
  325: *
  326: *     End of DLAED6
  327: *
  328:       END