dlamch.c

InsightToolkit-1.4.0(有大量的优化算法程序)

/*        Start to find EPS. */

        b = (doublereal) lbeta;
        i__1 = -lt;
        a = pow_di(&b, &i__1);
        leps = a;

/*        Try some tricks to see whether or not this is the correct  EPS. */

        b = 2. / 3;
        sixth = dlamc33_(&b, &half);
        third = dlamc3_(&sixth, &sixth);
        b = dlamc33_(&third, &half);
        b = dlamc3_(&b, &sixth);
        b = abs(b);
        if (b < leps) {
            b = leps;
        }

        leps = one;

        while (leps > b && b > zero) {
            leps = b;
            d__1 = half * leps;
            d__2 = 32. * (leps * leps);
            c = dlamc3_(&d__1, &d__2);
            c = dlamc33_(&half, &c);
            b = dlamc3_(&half, &c);
            c = dlamc33_(&half, &b);
            b = dlamc3_(&half, &c);
        }

        if (a < leps) {
            leps = a;
        }

/*        Computation of EPS complete.

          Now find  EMIN.  Let A = + or - 1, and + or - (1 + BASE**(-3)).
          Keep dividing  A by BETA until (gradual) underflow occurs. This
          is detected when we cannot recover the previous A.
*/

        rbase = one / lbeta;
        small = one;
        for (i = 1; i <= 3; ++i) {
            d__1 = small * rbase;
            small = dlamc3_(&d__1, &zero);
        }
        a = dlamc3_(&one, &small);
        dlamc4_(&ngpmin, &one, &lbeta);
        d__1 = -one;
        dlamc4_(&ngnmin, &d__1, &lbeta);
        dlamc4_(&gpmin, &a, &lbeta);
        d__1 = -a;
        dlamc4_(&gnmin, &d__1, &lbeta);
        ieee = FALSE_;

        if (ngpmin == ngnmin && gpmin == gnmin) {
            if (ngpmin == gpmin) {
                lemin = ngpmin;
/*            ( Non twos-complement machines, no gradual underflow; e.g.,  VAX ) */
            } else if (gpmin - ngpmin == 3) {
                lemin = ngpmin - 1 + lt;
                ieee = TRUE_;
/*            ( Non twos-complement machines, with gradual underflow; e.g., IEEE standard followers ) */
            } else {
                lemin = min(ngpmin,gpmin);
/*            ( A guess; no known machine ) */
                iwarn = TRUE_;
            }

        } else if (ngpmin == gpmin && ngnmin == gnmin) {
            if (abs(ngpmin - ngnmin) == 1) {
                lemin = max(ngpmin,ngnmin);
/*            ( Twos-complement machines, no gradual underflow; e.g., CYBER 205 ) */
            } else {
                lemin = min(ngpmin,ngnmin);
/*            ( A guess; no known machine ) */
                iwarn = TRUE_;
            }
        } else if (abs(ngpmin - ngnmin) == 1 && gpmin == gnmin) {
            if (gpmin - min(ngpmin,ngnmin) == 3) {
                lemin = max(ngpmin,ngnmin) - 1 + lt;
/*            ( Twos-complement machines with gradual underflow; no known machine ) */
            } else {
                lemin = min(ngpmin,ngnmin);
/*            ( A guess; no known machine ) */
                iwarn = TRUE_;
            }
        } else {
            lemin = min(min(min(ngpmin,ngnmin),gpmin),gnmin);
/*         ( A guess; no known machine ) */
            iwarn = TRUE_;
        }
/* ** Comment out this if block if EMIN is ok */
        if (iwarn) {
            first = TRUE_;
            printf("\n\n WARNING. The value EMIN may be incorrect:- ");
            printf("EMIN = %8i\n",lemin);
            printf("If, after inspection, the value EMIN looks acceptable");
            printf("please comment out \n the IF block as marked within the");
            printf("code of routine DLAMC2, \n otherwise supply EMIN");
            printf("explicitly.\n");
        }
/* **     Assume IEEE arithmetic if we found denormalised  numbers above,
          or if arithmetic seems to round in the  IEEE style,  determined
          in routine DLAMC1. A true IEEE machine should have both  things
          true; however, faulty machines may have one or the other.
*/

        ieee = ieee || lieee1;

/*        Compute  RMIN by successive division by  BETA. We could compute
          RMIN as BASE**( EMIN - 1 ),  but some machines underflow during
          this computation.
*/

        lrmin = one;
        for (i = 1; i <= 1-lemin; ++i) {
            d__1 = lrmin * rbase;
            lrmin = dlamc3_(&d__1, &zero);
        }

/*        Finally, call DLAMC5 to compute EMAX and RMAX. */

        dlamc5_(&lbeta, &lt, &lemin, &ieee, &lemax, &lrmax);
    }

    *beta = lbeta;
    *t = lt;
    *rnd = lrnd;
    *eps = leps;
    *emin = lemin;
    *rmin = lrmin;
    *emax = lemax;
    *rmax = lrmax;
} /* dlamc2_ */

doublereal dlamc3_(doublereal *a, doublereal *b)
{
/*  -- LAPACK auxiliary routine (version 2.0) --
       Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
       Courant Institute, Argonne National Lab, and Rice University
       October 31, 1992

    Purpose
    =======

    DLAMC3  is intended to force  A  and  B  to be stored prior to doing
    the addition of  A  and  B ,  for use in situations where optimizers
    might hold one of these in a register.

    Arguments
    =========

    A, B    (input) DOUBLE PRECISION
            The values A and B.

   =====================================================================
*/

    return *a + *b;
} /* dlamc3_ */

doublereal dlamc33_(doublereal *a, doublereal *b)
{
/*  Purpose
    =======

    As DLAMC3, but subtract A and B instead of adding them.

    Arguments
    =========

    A, B    (input) DOUBLE PRECISION
            The values A and B.

   =====================================================================
*/
    return *a - *b;

} /* dlamc33_ */

/* Subroutine */ void dlamc4_(integer *emin, doublereal *start, integer *base)
{
/*  -- LAPACK auxiliary routine (version 2.0) --
       Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
       Courant Institute, Argonne National Lab, and Rice University
       October 31, 1992

    Purpose
    =======

    DLAMC4 is a service routine for DLAMC2.

    Arguments
    =========

    EMIN    (output) EMIN
            The minimum exponent before (gradual) underflow, computed by
            setting A = START and dividing by BASE until the previous A
            can not be recovered.

    START   (input) DOUBLE PRECISION
            The starting point for determining EMIN.

    BASE    (input) INTEGER
            The base of the machine.

   =====================================================================
*/
    /* System generated locals */
    doublereal d__1;
    /* Local variables */
    static doublereal zero = 0., a;
    static integer i;
    static doublereal rbase, b1, b2, c1, c2, d1, d2;
    static doublereal one = 1.;

    a = *start;
    rbase = one / *base;
    *emin = 1;
    d__1 = a * rbase;
    b1 = dlamc3_(&d__1, &zero);
    c1 = c2 = d1 = d2 = a;
    while (c1 == a && c2 == a && d1 == a && d2 == a) {
        --(*emin);
        a = b1;
        d__1 = a / *base;
        b1 = dlamc3_(&d__1, &zero);
        d__1 = b1 * *base;
        c1 = dlamc3_(&d__1, &zero);
        d1 = zero;
        for (i = 1; i <= *base; ++i) {
            d1 += b1;
        }
        d__1 = a * rbase;
        b2 = dlamc3_(&d__1, &zero);
        d__1 = b2 / rbase;
        c2 = dlamc3_(&d__1, &zero);
        d2 = zero;
        for (i = 1; i <= *base; ++i) {
            d2 += b2;
        }
    }
} /* dlamc4_ */

/* Subroutine */ void dlamc5_(integer *beta, integer *p, integer *emin,
        logical *ieee, integer *emax, doublereal *rmax)
{
/*  -- LAPACK auxiliary routine (version 2.0) --
       Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
       Courant Institute, Argonne National Lab, and Rice University
       October 31, 1992

    Purpose
    =======

    DLAMC5 attempts to compute RMAX, the largest machine floating-point
    number, without overflow.  It assumes that EMAX + abs(EMIN) sum
    approximately to a power of 2.  It will fail on machines where this
    assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
    EMAX = 28718).  It will also fail if the value supplied for EMIN is
    too large (i.e. too close to zero), probably with overflow.

    Arguments
    =========

    BETA    (input) INTEGER
            The base of floating-point arithmetic.

    P       (input) INTEGER
            The number of base BETA digits in the mantissa of a
            floating-point value.

    EMIN    (input) INTEGER
            The minimum exponent before (gradual) underflow.

    IEEE    (input) LOGICAL
            A logical flag specifying whether or not the arithmetic
            system is thought to comply with the IEEE standard.

    EMAX    (output) INTEGER
            The largest exponent before overflow

    RMAX    (output) DOUBLE PRECISION
            The largest machine floating-point number.

   =====================================================================
*/

    /* Table of constant values */
    static doublereal c_b5 = 0.;
    /* System generated locals */
    doublereal d__1;
    /* Local variables */
    static integer lexp;
    static doublereal oldy;
    static integer uexp, i;
    static doublereal y, z;
    static integer nbits;
    static doublereal recbas;
    static integer exbits, expsum, try;

/*     First compute LEXP and UEXP, two powers of 2 that bound
       abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
       approximately to the bound that is closest to abs(EMIN).
       (EMAX is the exponent of the required number RMAX).
*/

    lexp = 1;
    exbits = 1;
    while ((try = lexp << 1) <= -(*emin)) {
        lexp = try;
        ++exbits;
    }
    if (lexp == -(*emin)) {
        uexp = lexp;
    } else {
        uexp = try;
        ++exbits;
    }

/*     Now -LEXP is less than or equal to EMIN, and -UEXP is greater
       than or equal to EMIN. EXBITS is the number of bits needed to
       store the exponent.
*/

    if (uexp + *emin > -lexp - *emin) {
        expsum = lexp << 1;
    } else {
        expsum = uexp << 1;
    }

/*     EXPSUM is the exponent range, approximately equal to EMAX - EMIN + 1 . */

    *emax = expsum + *emin - 1;
    nbits = exbits + 1 + *p;

/*     NBITS is the total number of bits needed to store a floating-point number. */

    if (nbits % 2 == 1 && *beta == 2) {

/*        Either there are an odd number of bits used to store a
          floating-point number, which is unlikely, or some bits are
          not used in the representation of numbers, which is possible,
          (e.g. Cray machines) or the mantissa has an implicit bit,
          (e.g. IEEE machines, Dec Vax machines), which is perhaps the
          most likely. We have to assume the last alternative.
          If this is true, then we need to reduce EMAX by one because
          there must be some way of representing zero in an implicit-bit
          system. On machines like Cray, we are reducing EMAX by one
          unnecessarily.
*/

        --(*emax);
    }

    if (*ieee) {

/*        Assume we are on an IEEE machine which reserves one exponent for infinity and NaN.  */

        --(*emax);
    }

/*     Now create RMAX, the largest machine number, which should
       be equal to (1.0 - BETA**(-P)) * BETA**EMAX .

       First compute 1.0 - BETA**(-P), being careful that the result is less than 1.0 . */

    recbas = 1. / *beta;
    z = *beta - 1.;
    y = 0.;
    for (i = 1; i <= *p; ++i) {
        z *= recbas;
        if (y < 1.) {
            oldy = y;
        }
        y = dlamc3_(&y, &z);
    }
    if (y >= 1.) {
        y = oldy;
    }

/*     Now multiply by BETA**EMAX to get RMAX. */

    for (i = 1; i <= *emax; ++i) {
        d__1 = y * *beta;
        y = dlamc3_(&d__1, &c_b5);
    }

    *rmax = y;
} /* dlamc5_ */