slamch.c

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

          that relevant values are stored  and not held in registers, or
          are not affected by optimizers.

          SLAMC1 returns the parameters  LBETA, LT, LRND and LIEEE1.
*/

        slamc1_(&lbeta, &lt, &lrnd, &lieee1);

/*        Start to find EPS. */

        b = (real) lbeta;
        i__1 = -lt;
        a = pow_ri(&b, &i__1);
        leps = a;

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

        b = 2.f / 3;
        r__1 = -half;
        sixth = slamc3_(&b, &r__1);
        third = slamc3_(&sixth, &sixth);
        b = slamc3_(&third, &r__1);
        b = slamc3_(&b, &sixth);
        b = abs(b);
        if (b < leps) {
            b = leps;
        }

        leps = one;

        while (leps > b && b > zero) {
            leps = b;
            r__1 = half * leps;
            r__2 = 32.0f * leps * leps;
            c = slamc3_(&r__1, &r__2);
            r__1 = -c;
            c = slamc3_(&half, &r__1);
            b = slamc3_(&half, &c);
            r__1 = -b;
            c = slamc3_(&half, &r__1);
            b = slamc3_(&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) {
            r__1 = small * rbase;
            small = slamc3_(&r__1, &zero);
        }
        a = slamc3_(&one, &small);
        slamc4_(&ngpmin, &one, &lbeta);
        r__1 = -one;
        slamc4_(&ngnmin, &r__1, &lbeta);
        slamc4_(&gpmin, &a, &lbeta);
        r__1 = -a;
        slamc4_(&gnmin, &r__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 SLAMC2, \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 SLAMC1. 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) {
            r__1 = lrmin * rbase;
            lrmin = slamc3_(&r__1, &zero);
        }

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

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

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


// Microsoft Visual C++ 2003 produces bad code when the following
// routine is optimized.  Turn off the optimization for this one
// routine and turn back on any optimizations after this routine.
#if (_MSC_VER >= 1310)
#pragma optimize("", off)
#endif

real slamc3_(real *a, real *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
    =======

    SLAMC3  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) REAL
            The values A and B.

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

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

// Turn the optimizations back on for Visual Studio .NET 2003
#if (_MSC_VER >= 1310)
#pragma optimize("", on)
#endif


/* Subroutine */ void slamc4_(integer *emin, real *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
    =======

    SLAMC4 is a service routine for SLAMC2.

    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) REAL
            The starting point for determining EMIN.

    BASE    (input) INTEGER
            The base of the machine.

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

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

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

/* Subroutine */ void slamc5_(integer *beta, integer *p, integer *emin,
        logical *ieee, integer *emax, real *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
    =======

    SLAMC5 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) REAL
            The largest machine floating-point number.

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

    /* Table of constant values */
    static real c_b5 = 0.f;
    /* System generated locals */
    real r__1;
    /* Local variables */
    static integer lexp;
    static real oldy;
    static integer uexp, i;
    static real y, z;
    static integer nbits;
    static real 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.f / *beta;
    z = *beta - 1.f;
    y = 0.f;
    for (i = 1; i <= *p; ++i) {
        z *= recbas;
        if (y < 1.f) {
            oldy = y;
        }
        y = slamc3_(&y, &z);
    }
    if (y >= 1.f) {
        y = oldy;
    }

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

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

    *rmax = y;
} /* slamc5_ */