slamch.c

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

/* L30: */
        }

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

/*<          CALL SLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX ) >*/
        slamc5_(&lbeta, &lt, &lemin, &ieee, &lemax, &lrmax);
/*<       END IF >*/
    }

/*<       BETA = LBETA >*/
    *beta = lbeta;
/*<       T = LT >*/
    *t = lt;
/*<       RND = LRND >*/
    *rnd = lrnd;
/*<       EPS = LEPS >*/
    *eps = leps;
/*<       EMIN = LEMIN >*/
    *emin = lemin;
/*<       RMIN = LRMIN >*/
    *rmin = lrmin;
/*<       EMAX = LEMAX >*/
    *emax = lemax;
/*<       RMAX = LRMAX >*/
    *rmax = lrmax;

/*<       RETURN >*/
    return 0;

/*<  9 >*/

/*     End of SLAMC2 */

/*<       END >*/
} /* slamc2_ */


/* *********************************************************************** */

/*<       REAL             FUNCTION SLAMC3( A, B ) >*/
doublereal slamc3_(real *a, real *b)
{
    /* System generated locals */
    real ret_val;


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

/*     .. Scalar Arguments .. */
/*<       REAL               A, B >*/
/*     .. */

/*  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. */

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

/*     .. Executable Statements .. */

/*<       SLAMC3 = A + B >*/
    ret_val = *a + *b;

/*<       RETURN >*/
    return ret_val;

/*     End of SLAMC3 */

/*<       END >*/
} /* slamc3_ */


/* *********************************************************************** */

/*<       SUBROUTINE SLAMC4( EMIN, START, BASE ) >*/
/* Subroutine */ int slamc4_(integer *emin, real *start, integer *base)
{
    /* System generated locals */
    integer i__1;
    real r__1;

    /* Local variables */
    real a;
    integer i__;
    real b1, b2, c1, c2, d1, d2, one, zero, rbase;
    extern doublereal slamc3_(real *, real *);


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

/*     .. Scalar Arguments .. */
/*<       INTEGER            BASE, EMIN >*/
/*<       REAL               START >*/
/*     .. */

/*  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. */

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

/*     .. Local Scalars .. */
/*<       INTEGER            I >*/
/*<       REAL               A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO >*/
/*     .. */
/*     .. External Functions .. */
/*<       REAL               SLAMC3 >*/
/*<       EXTERNAL           SLAMC3 >*/
/*     .. */
/*     .. Executable Statements .. */

/*<       A = START >*/
    a = *start;
/*<       ONE = 1 >*/
    one = (float)1.;
/*<       RBASE = ONE / BASE >*/
    rbase = one / *base;
/*<       ZERO = 0 >*/
    zero = (float)0.;
/*<       EMIN = 1 >*/
    *emin = 1;
/*<       B1 = SLAMC3( A*RBASE, ZERO ) >*/
    r__1 = a * rbase;
    b1 = slamc3_(&r__1, &zero);
/*<       C1 = A >*/
    c1 = a;
/*<       C2 = A >*/
    c2 = a;
/*<       D1 = A >*/
    d1 = a;
/*<       D2 = A >*/
    d2 = a;
/* +    WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. */
/*    $       ( D1.EQ.A ).AND.( D2.EQ.A )      )LOOP */
/*<    10 CONTINUE >*/
L10:
/*<    >*/
    if (c1 == a && c2 == a && d1 == a && d2 == a) {
/*<          EMIN = EMIN - 1 >*/
        --(*emin);
/*<          A = B1 >*/
        a = b1;
/*<          B1 = SLAMC3( A / BASE, ZERO ) >*/
        r__1 = a / *base;
        b1 = slamc3_(&r__1, &zero);
/*<          C1 = SLAMC3( B1*BASE, ZERO ) >*/
        r__1 = b1 * *base;
        c1 = slamc3_(&r__1, &zero);
/*<          D1 = ZERO >*/
        d1 = zero;
/*<          DO 20 I = 1, BASE >*/
        i__1 = *base;
        for (i__ = 1; i__ <= i__1; ++i__) {
/*<             D1 = D1 + B1 >*/
            d1 += b1;
/*<    20    CONTINUE >*/
/* L20: */
        }
/*<          B2 = SLAMC3( A*RBASE, ZERO ) >*/
        r__1 = a * rbase;
        b2 = slamc3_(&r__1, &zero);
/*<          C2 = SLAMC3( B2 / RBASE, ZERO ) >*/
        r__1 = b2 / rbase;
        c2 = slamc3_(&r__1, &zero);
/*<          D2 = ZERO >*/
        d2 = zero;
/*<          DO 30 I = 1, BASE >*/
        i__1 = *base;
        for (i__ = 1; i__ <= i__1; ++i__) {
/*<             D2 = D2 + B2 >*/
            d2 += b2;
/*<    30    CONTINUE >*/
/* L30: */
        }
/*<          GO TO 10 >*/
        goto L10;
/*<       END IF >*/
    }
/* +    END WHILE */

/*<       RETURN >*/
    return 0;

/*     End of SLAMC4 */

/*<       END >*/
} /* slamc4_ */


/* *********************************************************************** */

/*<       SUBROUTINE SLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX ) >*/
/* Subroutine */ int slamc5_(integer *beta, integer *p, integer *emin, 
        logical *ieee, integer *emax, real *rmax)
{
    /* System generated locals */
    integer i__1;
    real r__1;

    /* Local variables */
    integer i__;
    real y, z__;
    integer try__, lexp;
    real oldy=0;
    integer uexp, nbits;
    extern doublereal slamc3_(real *, real *);
    real recbas;
    integer exbits, expsum;


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

/*     .. Scalar Arguments .. */
/*<       LOGICAL            IEEE >*/
/*<       INTEGER            BETA, EMAX, EMIN, P >*/
/*<       REAL               RMAX >*/
/*     .. */

/*  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. */

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

/*     .. Parameters .. */
/*<       REAL               ZERO, ONE >*/
/*<       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 ) >*/
/*     .. */
/*     .. Local Scalars .. */
/*<       INTEGER            EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP >*/
/*<       REAL               OLDY, RECBAS, Y, Z >*/
/*     .. */
/*     .. External Functions .. */
/*<       REAL               SLAMC3 >*/
/*<       EXTERNAL           SLAMC3 >*/
/*     .. */
/*     .. Intrinsic Functions .. */
/*<       INTRINSIC          MOD >*/
/*     .. */
/*     .. Executable Statements .. */

/*     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 >*/
    lexp = 1;
/*<       EXBITS = 1 >*/
    exbits = 1;
/*<    10 CONTINUE >*/
L10:
/*<       TRY = LEXP*2 >*/
    try__ = lexp << 1;
/*<       IF( TRY.LE.( -EMIN ) ) THEN >*/
    if (try__ <= -(*emin)) {
/*<          LEXP = TRY >*/
        lexp = try__;
/*<          EXBITS = EXBITS + 1 >*/
        ++exbits;
/*<          GO TO 10 >*/
        goto L10;
/*<       END IF >*/
    }
/*<       IF( LEXP.EQ.-EMIN ) THEN >*/
    if (lexp == -(*emin)) {
/*<          UEXP = LEXP >*/
        uexp = lexp;
/*<       ELSE >*/
    } else {
/*<          UEXP = TRY >*/
        uexp = try__;
/*<          EXBITS = EXBITS + 1 >*/
        ++exbits;
/*<       END IF >*/
    }

/*     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 ).GT.( -LEXP-EMIN ) ) THEN >*/
    if (uexp + *emin > -lexp - *emin) {
/*<          EXPSUM = 2*LEXP >*/
        expsum = lexp << 1;
/*<       ELSE >*/
    } else {
/*<          EXPSUM = 2*UEXP >*/
        expsum = uexp << 1;
/*<       END IF >*/
    }

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

/*<       EMAX = EXPSUM + EMIN - 1 >*/
    *emax = expsum + *emin - 1;
/*<       NBITS = 1 + EXBITS + P >*/
    nbits = exbits + 1 + *p;

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

/*<       IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN >*/
    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 = EMAX - 1 >*/
        --(*emax);
/*<       END IF >*/
    }

/*<       IF( IEEE ) THEN >*/
    if (*ieee) {

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

/*<          EMAX = EMAX - 1 >*/
        --(*emax);
/*<       END IF >*/
    }

/*     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 = ONE / BETA >*/
    recbas = (float)1. / *beta;
/*<       Z = BETA - ONE >*/
    z__ = *beta - (float)1.;
/*<       Y = ZERO >*/
    y = (float)0.;
/*<       DO 20 I = 1, P >*/
    i__1 = *p;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          Z = Z*RECBAS >*/
        z__ *= recbas;
/*<    >*/
        if (y < (float)1.) {
            oldy = y;
        }
/*<          Y = SLAMC3( Y, Z ) >*/
        y = slamc3_(&y, &z__);
/*<    20 CONTINUE >*/
/* L20: */
    }
/*<    >*/
    if (y >= (float)1.) {
        y = oldy;
    }

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

/*<       DO 30 I = 1, EMAX >*/
    i__1 = *emax;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          Y = SLAMC3( Y*BETA, ZERO ) >*/
        r__1 = y * *beta;
        y = slamc3_(&r__1, &c_b32);
/*<    30 CONTINUE >*/
/* L30: */
    }

/*<       RMAX = Y >*/
    *rmax = y;
/*<       RETURN >*/
    return 0;

/*     End of SLAMC5 */

/*<       END >*/
} /* slamc5_ */

#ifdef __cplusplus
        }
#endif