Annotation of rpl/lapack/lapack/dlaed6.f, revision 1.21

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

CVSweb interface <joel.bertrand@systella.fr>