Chris@202: /* dlamch.f -- translated by f2c (version 20061008). Chris@202: You must link the resulting object file with libf2c: Chris@202: on Microsoft Windows system, link with libf2c.lib; Chris@202: on Linux or Unix systems, link with .../path/to/libf2c.a -lm Chris@202: or, if you install libf2c.a in a standard place, with -lf2c -lm Chris@202: -- in that order, at the end of the command line, as in Chris@202: cc *.o -lf2c -lm Chris@202: Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., Chris@202: Chris@202: http://www.netlib.org/f2c/libf2c.zip Chris@202: */ Chris@202: Chris@202: #include "f2c.h" Chris@202: #include "blaswrap.h" Chris@202: Chris@202: /* Table of constant values */ Chris@202: Chris@202: static integer c__1 = 1; Chris@202: static doublereal c_b32 = 0.; Chris@202: Chris@202: doublereal dlamch_(char *cmach) Chris@202: { Chris@202: /* Initialized data */ Chris@202: Chris@202: static logical first = TRUE_; Chris@202: Chris@202: /* System generated locals */ Chris@202: integer i__1; Chris@202: doublereal ret_val; Chris@202: Chris@202: /* Builtin functions */ Chris@202: double pow_di(doublereal *, integer *); Chris@202: Chris@202: /* Local variables */ Chris@202: static doublereal t; Chris@202: integer it; Chris@202: static doublereal rnd, eps, base; Chris@202: integer beta; Chris@202: static doublereal emin, prec, emax; Chris@202: integer imin, imax; Chris@202: logical lrnd; Chris@202: static doublereal rmin, rmax; Chris@202: doublereal rmach; Chris@202: extern logical lsame_(char *, char *); Chris@202: doublereal small; Chris@202: static doublereal sfmin; Chris@202: extern /* Subroutine */ int dlamc2_(integer *, integer *, logical *, Chris@202: doublereal *, integer *, doublereal *, integer *, doublereal *); Chris@202: Chris@202: Chris@202: /* -- LAPACK auxiliary routine (version 3.2) -- */ Chris@202: /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ Chris@202: /* November 2006 */ Chris@202: Chris@202: /* .. Scalar Arguments .. */ Chris@202: /* .. */ Chris@202: Chris@202: /* Purpose */ Chris@202: /* ======= */ Chris@202: Chris@202: /* DLAMCH determines double precision machine parameters. */ Chris@202: Chris@202: /* Arguments */ Chris@202: /* ========= */ Chris@202: Chris@202: /* CMACH (input) CHARACTER*1 */ Chris@202: /* Specifies the value to be returned by DLAMCH: */ Chris@202: /* = 'E' or 'e', DLAMCH := eps */ Chris@202: /* = 'S' or 's , DLAMCH := sfmin */ Chris@202: /* = 'B' or 'b', DLAMCH := base */ Chris@202: /* = 'P' or 'p', DLAMCH := eps*base */ Chris@202: /* = 'N' or 'n', DLAMCH := t */ Chris@202: /* = 'R' or 'r', DLAMCH := rnd */ Chris@202: /* = 'M' or 'm', DLAMCH := emin */ Chris@202: /* = 'U' or 'u', DLAMCH := rmin */ Chris@202: /* = 'L' or 'l', DLAMCH := emax */ Chris@202: /* = 'O' or 'o', DLAMCH := rmax */ Chris@202: Chris@202: /* where */ Chris@202: Chris@202: /* eps = relative machine precision */ Chris@202: /* sfmin = safe minimum, such that 1/sfmin does not overflow */ Chris@202: /* base = base of the machine */ Chris@202: /* prec = eps*base */ Chris@202: /* t = number of (base) digits in the mantissa */ Chris@202: /* rnd = 1.0 when rounding occurs in addition, 0.0 otherwise */ Chris@202: /* emin = minimum exponent before (gradual) underflow */ Chris@202: /* rmin = underflow threshold - base**(emin-1) */ Chris@202: /* emax = largest exponent before overflow */ Chris@202: /* rmax = overflow threshold - (base**emax)*(1-eps) */ Chris@202: Chris@202: /* ===================================================================== */ Chris@202: Chris@202: /* .. Parameters .. */ Chris@202: /* .. */ Chris@202: /* .. Local Scalars .. */ Chris@202: /* .. */ Chris@202: /* .. External Functions .. */ Chris@202: /* .. */ Chris@202: /* .. External Subroutines .. */ Chris@202: /* .. */ Chris@202: /* .. Save statement .. */ Chris@202: /* .. */ Chris@202: /* .. Data statements .. */ Chris@202: /* .. */ Chris@202: /* .. Executable Statements .. */ Chris@202: Chris@202: if (first) { Chris@202: dlamc2_(&beta, &it, &lrnd, &eps, &imin, &rmin, &imax, &rmax); Chris@202: base = (doublereal) beta; Chris@202: t = (doublereal) it; Chris@202: if (lrnd) { Chris@202: rnd = 1.; Chris@202: i__1 = 1 - it; Chris@202: eps = pow_di(&base, &i__1) / 2; Chris@202: } else { Chris@202: rnd = 0.; Chris@202: i__1 = 1 - it; Chris@202: eps = pow_di(&base, &i__1); Chris@202: } Chris@202: prec = eps * base; Chris@202: emin = (doublereal) imin; Chris@202: emax = (doublereal) imax; Chris@202: sfmin = rmin; Chris@202: small = 1. / rmax; Chris@202: if (small >= sfmin) { Chris@202: Chris@202: /* Use SMALL plus a bit, to avoid the possibility of rounding */ Chris@202: /* causing overflow when computing 1/sfmin. */ Chris@202: Chris@202: sfmin = small * (eps + 1.); Chris@202: } Chris@202: } Chris@202: Chris@202: if (lsame_(cmach, "E")) { Chris@202: rmach = eps; Chris@202: } else if (lsame_(cmach, "S")) { Chris@202: rmach = sfmin; Chris@202: } else if (lsame_(cmach, "B")) { Chris@202: rmach = base; Chris@202: } else if (lsame_(cmach, "P")) { Chris@202: rmach = prec; Chris@202: } else if (lsame_(cmach, "N")) { Chris@202: rmach = t; Chris@202: } else if (lsame_(cmach, "R")) { Chris@202: rmach = rnd; Chris@202: } else if (lsame_(cmach, "M")) { Chris@202: rmach = emin; Chris@202: } else if (lsame_(cmach, "U")) { Chris@202: rmach = rmin; Chris@202: } else if (lsame_(cmach, "L")) { Chris@202: rmach = emax; Chris@202: } else if (lsame_(cmach, "O")) { Chris@202: rmach = rmax; Chris@202: } Chris@202: Chris@202: ret_val = rmach; Chris@202: first = FALSE_; Chris@202: return ret_val; Chris@202: Chris@202: /* End of DLAMCH */ Chris@202: Chris@202: } /* dlamch_ */ Chris@202: Chris@202: Chris@202: /* *********************************************************************** */ Chris@202: Chris@202: /* Subroutine */ int dlamc1_(integer *beta, integer *t, logical *rnd, logical Chris@202: *ieee1) Chris@202: { Chris@202: /* Initialized data */ Chris@202: Chris@202: static logical first = TRUE_; Chris@202: Chris@202: /* System generated locals */ Chris@202: doublereal d__1, d__2; Chris@202: Chris@202: /* Local variables */ Chris@202: doublereal a, b, c__, f, t1, t2; Chris@202: static integer lt; Chris@202: doublereal one, qtr; Chris@202: static logical lrnd; Chris@202: static integer lbeta; Chris@202: doublereal savec; Chris@202: extern doublereal dlamc3_(doublereal *, doublereal *); Chris@202: static logical lieee1; Chris@202: Chris@202: Chris@202: /* -- LAPACK auxiliary routine (version 3.2) -- */ Chris@202: /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ Chris@202: /* November 2006 */ Chris@202: Chris@202: /* .. Scalar Arguments .. */ Chris@202: /* .. */ Chris@202: Chris@202: /* Purpose */ Chris@202: /* ======= */ Chris@202: Chris@202: /* DLAMC1 determines the machine parameters given by BETA, T, RND, and */ Chris@202: /* IEEE1. */ Chris@202: Chris@202: /* Arguments */ Chris@202: /* ========= */ Chris@202: Chris@202: /* BETA (output) INTEGER */ Chris@202: /* The base of the machine. */ Chris@202: Chris@202: /* T (output) INTEGER */ Chris@202: /* The number of ( BETA ) digits in the mantissa. */ Chris@202: Chris@202: /* RND (output) LOGICAL */ Chris@202: /* Specifies whether proper rounding ( RND = .TRUE. ) or */ Chris@202: /* chopping ( RND = .FALSE. ) occurs in addition. This may not */ Chris@202: /* be a reliable guide to the way in which the machine performs */ Chris@202: /* its arithmetic. */ Chris@202: Chris@202: /* IEEE1 (output) LOGICAL */ Chris@202: /* Specifies whether rounding appears to be done in the IEEE */ Chris@202: /* 'round to nearest' style. */ Chris@202: Chris@202: /* Further Details */ Chris@202: /* =============== */ Chris@202: Chris@202: /* The routine is based on the routine ENVRON by Malcolm and */ Chris@202: /* incorporates suggestions by Gentleman and Marovich. See */ Chris@202: Chris@202: /* Malcolm M. A. (1972) Algorithms to reveal properties of */ Chris@202: /* floating-point arithmetic. Comms. of the ACM, 15, 949-951. */ Chris@202: Chris@202: /* Gentleman W. M. and Marovich S. B. (1974) More on algorithms */ Chris@202: /* that reveal properties of floating point arithmetic units. */ Chris@202: /* Comms. of the ACM, 17, 276-277. */ Chris@202: Chris@202: /* ===================================================================== */ Chris@202: Chris@202: /* .. Local Scalars .. */ Chris@202: /* .. */ Chris@202: /* .. External Functions .. */ Chris@202: /* .. */ Chris@202: /* .. Save statement .. */ Chris@202: /* .. */ Chris@202: /* .. Data statements .. */ Chris@202: /* .. */ Chris@202: /* .. Executable Statements .. */ Chris@202: Chris@202: if (first) { Chris@202: one = 1.; Chris@202: Chris@202: /* LBETA, LIEEE1, LT and LRND are the local values of BETA, */ Chris@202: /* IEEE1, T and RND. */ Chris@202: Chris@202: /* Throughout this routine we use the function DLAMC3 to ensure */ Chris@202: /* that relevant values are stored and not held in registers, or */ Chris@202: /* are not affected by optimizers. */ Chris@202: Chris@202: /* Compute a = 2.0**m with the smallest positive integer m such */ Chris@202: /* that */ Chris@202: Chris@202: /* fl( a + 1.0 ) = a. */ Chris@202: Chris@202: a = 1.; Chris@202: c__ = 1.; Chris@202: Chris@202: /* + WHILE( C.EQ.ONE )LOOP */ Chris@202: L10: Chris@202: if (c__ == one) { Chris@202: a *= 2; Chris@202: c__ = dlamc3_(&a, &one); Chris@202: d__1 = -a; Chris@202: c__ = dlamc3_(&c__, &d__1); Chris@202: goto L10; Chris@202: } Chris@202: /* + END WHILE */ Chris@202: Chris@202: /* Now compute b = 2.0**m with the smallest positive integer m */ Chris@202: /* such that */ Chris@202: Chris@202: /* fl( a + b ) .gt. a. */ Chris@202: Chris@202: b = 1.; Chris@202: c__ = dlamc3_(&a, &b); Chris@202: Chris@202: /* + WHILE( C.EQ.A )LOOP */ Chris@202: L20: Chris@202: if (c__ == a) { Chris@202: b *= 2; Chris@202: c__ = dlamc3_(&a, &b); Chris@202: goto L20; Chris@202: } Chris@202: /* + END WHILE */ Chris@202: Chris@202: /* Now compute the base. a and c are neighbouring floating point */ Chris@202: /* numbers in the interval ( beta**t, beta**( t + 1 ) ) and so */ Chris@202: /* their difference is beta. Adding 0.25 to c is to ensure that it */ Chris@202: /* is truncated to beta and not ( beta - 1 ). */ Chris@202: Chris@202: qtr = one / 4; Chris@202: savec = c__; Chris@202: d__1 = -a; Chris@202: c__ = dlamc3_(&c__, &d__1); Chris@202: lbeta = (integer) (c__ + qtr); Chris@202: Chris@202: /* Now determine whether rounding or chopping occurs, by adding a */ Chris@202: /* bit less than beta/2 and a bit more than beta/2 to a. */ Chris@202: Chris@202: b = (doublereal) lbeta; Chris@202: d__1 = b / 2; Chris@202: d__2 = -b / 100; Chris@202: f = dlamc3_(&d__1, &d__2); Chris@202: c__ = dlamc3_(&f, &a); Chris@202: if (c__ == a) { Chris@202: lrnd = TRUE_; Chris@202: } else { Chris@202: lrnd = FALSE_; Chris@202: } Chris@202: d__1 = b / 2; Chris@202: d__2 = b / 100; Chris@202: f = dlamc3_(&d__1, &d__2); Chris@202: c__ = dlamc3_(&f, &a); Chris@202: if (lrnd && c__ == a) { Chris@202: lrnd = FALSE_; Chris@202: } Chris@202: Chris@202: /* Try and decide whether rounding is done in the IEEE 'round to */ Chris@202: /* nearest' style. B/2 is half a unit in the last place of the two */ Chris@202: /* numbers A and SAVEC. Furthermore, A is even, i.e. has last bit */ Chris@202: /* zero, and SAVEC is odd. Thus adding B/2 to A should not change */ Chris@202: /* A, but adding B/2 to SAVEC should change SAVEC. */ Chris@202: Chris@202: d__1 = b / 2; Chris@202: t1 = dlamc3_(&d__1, &a); Chris@202: d__1 = b / 2; Chris@202: t2 = dlamc3_(&d__1, &savec); Chris@202: lieee1 = t1 == a && t2 > savec && lrnd; Chris@202: Chris@202: /* Now find the mantissa, t. It should be the integer part of */ Chris@202: /* log to the base beta of a, however it is safer to determine t */ Chris@202: /* by powering. So we find t as the smallest positive integer for */ Chris@202: /* which */ Chris@202: Chris@202: /* fl( beta**t + 1.0 ) = 1.0. */ Chris@202: Chris@202: lt = 0; Chris@202: a = 1.; Chris@202: c__ = 1.; Chris@202: Chris@202: /* + WHILE( C.EQ.ONE )LOOP */ Chris@202: L30: Chris@202: if (c__ == one) { Chris@202: ++lt; Chris@202: a *= lbeta; Chris@202: c__ = dlamc3_(&a, &one); Chris@202: d__1 = -a; Chris@202: c__ = dlamc3_(&c__, &d__1); Chris@202: goto L30; Chris@202: } Chris@202: /* + END WHILE */ Chris@202: Chris@202: } Chris@202: Chris@202: *beta = lbeta; Chris@202: *t = lt; Chris@202: *rnd = lrnd; Chris@202: *ieee1 = lieee1; Chris@202: first = FALSE_; Chris@202: return 0; Chris@202: Chris@202: /* End of DLAMC1 */ Chris@202: Chris@202: } /* dlamc1_ */ Chris@202: Chris@202: Chris@202: /* *********************************************************************** */ Chris@202: Chris@202: /* Subroutine */ int dlamc2_(integer *beta, integer *t, logical *rnd, Chris@202: doublereal *eps, integer *emin, doublereal *rmin, integer *emax, Chris@202: doublereal *rmax) Chris@202: { Chris@202: /* Initialized data */ Chris@202: Chris@202: static logical first = TRUE_; Chris@202: static logical iwarn = FALSE_; Chris@202: Chris@202: /* Format strings */ Chris@202: static char fmt_9999[] = "(//\002 WARNING. The value EMIN may be incorre" Chris@202: "ct:-\002,\002 EMIN = \002,i8,/\002 If, after inspection, the va" Chris@202: "lue EMIN looks\002,\002 acceptable please comment out \002,/\002" Chris@202: " the IF block as marked within the code of routine\002,\002 DLAM" Chris@202: "C2,\002,/\002 otherwise supply EMIN explicitly.\002,/)"; Chris@202: Chris@202: /* System generated locals */ Chris@202: integer i__1; Chris@202: doublereal d__1, d__2, d__3, d__4, d__5; Chris@202: Chris@202: /* Builtin functions */ Chris@202: double pow_di(doublereal *, integer *); Chris@202: integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void); Chris@202: Chris@202: /* Local variables */ Chris@202: doublereal a, b, c__; Chris@202: integer i__; Chris@202: static integer lt; Chris@202: doublereal one, two; Chris@202: logical ieee; Chris@202: doublereal half; Chris@202: logical lrnd; Chris@202: static doublereal leps; Chris@202: doublereal zero; Chris@202: static integer lbeta; Chris@202: doublereal rbase; Chris@202: static integer lemin, lemax; Chris@202: integer gnmin; Chris@202: doublereal small; Chris@202: integer gpmin; Chris@202: doublereal third; Chris@202: static doublereal lrmin, lrmax; Chris@202: doublereal sixth; Chris@202: extern /* Subroutine */ int dlamc1_(integer *, integer *, logical *, Chris@202: logical *); Chris@202: extern doublereal dlamc3_(doublereal *, doublereal *); Chris@202: logical lieee1; Chris@202: extern /* Subroutine */ int dlamc4_(integer *, doublereal *, integer *), Chris@202: dlamc5_(integer *, integer *, integer *, logical *, integer *, Chris@202: doublereal *); Chris@202: integer ngnmin, ngpmin; Chris@202: Chris@202: /* Fortran I/O blocks */ Chris@202: static cilist io___58 = { 0, 6, 0, fmt_9999, 0 }; Chris@202: Chris@202: Chris@202: Chris@202: /* -- LAPACK auxiliary routine (version 3.2) -- */ Chris@202: /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ Chris@202: /* November 2006 */ Chris@202: Chris@202: /* .. Scalar Arguments .. */ Chris@202: /* .. */ Chris@202: Chris@202: /* Purpose */ Chris@202: /* ======= */ Chris@202: Chris@202: /* DLAMC2 determines the machine parameters specified in its argument */ Chris@202: /* list. */ Chris@202: Chris@202: /* Arguments */ Chris@202: /* ========= */ Chris@202: Chris@202: /* BETA (output) INTEGER */ Chris@202: /* The base of the machine. */ Chris@202: Chris@202: /* T (output) INTEGER */ Chris@202: /* The number of ( BETA ) digits in the mantissa. */ Chris@202: Chris@202: /* RND (output) LOGICAL */ Chris@202: /* Specifies whether proper rounding ( RND = .TRUE. ) or */ Chris@202: /* chopping ( RND = .FALSE. ) occurs in addition. This may not */ Chris@202: /* be a reliable guide to the way in which the machine performs */ Chris@202: /* its arithmetic. */ Chris@202: Chris@202: /* EPS (output) DOUBLE PRECISION */ Chris@202: /* The smallest positive number such that */ Chris@202: Chris@202: /* fl( 1.0 - EPS ) .LT. 1.0, */ Chris@202: Chris@202: /* where fl denotes the computed value. */ Chris@202: Chris@202: /* EMIN (output) INTEGER */ Chris@202: /* The minimum exponent before (gradual) underflow occurs. */ Chris@202: Chris@202: /* RMIN (output) DOUBLE PRECISION */ Chris@202: /* The smallest normalized number for the machine, given by */ Chris@202: /* BASE**( EMIN - 1 ), where BASE is the floating point value */ Chris@202: /* of BETA. */ Chris@202: Chris@202: /* EMAX (output) INTEGER */ Chris@202: /* The maximum exponent before overflow occurs. */ Chris@202: Chris@202: /* RMAX (output) DOUBLE PRECISION */ Chris@202: /* The largest positive number for the machine, given by */ Chris@202: /* BASE**EMAX * ( 1 - EPS ), where BASE is the floating point */ Chris@202: /* value of BETA. */ Chris@202: Chris@202: /* Further Details */ Chris@202: /* =============== */ Chris@202: Chris@202: /* The computation of EPS is based on a routine PARANOIA by */ Chris@202: /* W. Kahan of the University of California at Berkeley. */ Chris@202: Chris@202: /* ===================================================================== */ Chris@202: Chris@202: /* .. Local Scalars .. */ Chris@202: /* .. */ Chris@202: /* .. External Functions .. */ Chris@202: /* .. */ Chris@202: /* .. External Subroutines .. */ Chris@202: /* .. */ Chris@202: /* .. Intrinsic Functions .. */ Chris@202: /* .. */ Chris@202: /* .. Save statement .. */ Chris@202: /* .. */ Chris@202: /* .. Data statements .. */ Chris@202: /* .. */ Chris@202: /* .. Executable Statements .. */ Chris@202: Chris@202: if (first) { Chris@202: zero = 0.; Chris@202: one = 1.; Chris@202: two = 2.; Chris@202: Chris@202: /* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of */ Chris@202: /* BETA, T, RND, EPS, EMIN and RMIN. */ Chris@202: Chris@202: /* Throughout this routine we use the function DLAMC3 to ensure */ Chris@202: /* that relevant values are stored and not held in registers, or */ Chris@202: /* are not affected by optimizers. */ Chris@202: Chris@202: /* DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. */ Chris@202: Chris@202: dlamc1_(&lbeta, <, &lrnd, &lieee1); Chris@202: Chris@202: /* Start to find EPS. */ Chris@202: Chris@202: b = (doublereal) lbeta; Chris@202: i__1 = -lt; Chris@202: a = pow_di(&b, &i__1); Chris@202: leps = a; Chris@202: Chris@202: /* Try some tricks to see whether or not this is the correct EPS. */ Chris@202: Chris@202: b = two / 3; Chris@202: half = one / 2; Chris@202: d__1 = -half; Chris@202: sixth = dlamc3_(&b, &d__1); Chris@202: third = dlamc3_(&sixth, &sixth); Chris@202: d__1 = -half; Chris@202: b = dlamc3_(&third, &d__1); Chris@202: b = dlamc3_(&b, &sixth); Chris@202: b = abs(b); Chris@202: if (b < leps) { Chris@202: b = leps; Chris@202: } Chris@202: Chris@202: leps = 1.; Chris@202: Chris@202: /* + WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP */ Chris@202: L10: Chris@202: if (leps > b && b > zero) { Chris@202: leps = b; Chris@202: d__1 = half * leps; Chris@202: /* Computing 5th power */ Chris@202: d__3 = two, d__4 = d__3, d__3 *= d__3; Chris@202: /* Computing 2nd power */ Chris@202: d__5 = leps; Chris@202: d__2 = d__4 * (d__3 * d__3) * (d__5 * d__5); Chris@202: c__ = dlamc3_(&d__1, &d__2); Chris@202: d__1 = -c__; Chris@202: c__ = dlamc3_(&half, &d__1); Chris@202: b = dlamc3_(&half, &c__); Chris@202: d__1 = -b; Chris@202: c__ = dlamc3_(&half, &d__1); Chris@202: b = dlamc3_(&half, &c__); Chris@202: goto L10; Chris@202: } Chris@202: /* + END WHILE */ Chris@202: Chris@202: if (a < leps) { Chris@202: leps = a; Chris@202: } Chris@202: Chris@202: /* Computation of EPS complete. */ Chris@202: Chris@202: /* Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)). */ Chris@202: /* Keep dividing A by BETA until (gradual) underflow occurs. This */ Chris@202: /* is detected when we cannot recover the previous A. */ Chris@202: Chris@202: rbase = one / lbeta; Chris@202: small = one; Chris@202: for (i__ = 1; i__ <= 3; ++i__) { Chris@202: d__1 = small * rbase; Chris@202: small = dlamc3_(&d__1, &zero); Chris@202: /* L20: */ Chris@202: } Chris@202: a = dlamc3_(&one, &small); Chris@202: dlamc4_(&ngpmin, &one, &lbeta); Chris@202: d__1 = -one; Chris@202: dlamc4_(&ngnmin, &d__1, &lbeta); Chris@202: dlamc4_(&gpmin, &a, &lbeta); Chris@202: d__1 = -a; Chris@202: dlamc4_(&gnmin, &d__1, &lbeta); Chris@202: ieee = FALSE_; Chris@202: Chris@202: if (ngpmin == ngnmin && gpmin == gnmin) { Chris@202: if (ngpmin == gpmin) { Chris@202: lemin = ngpmin; Chris@202: /* ( Non twos-complement machines, no gradual underflow; */ Chris@202: /* e.g., VAX ) */ Chris@202: } else if (gpmin - ngpmin == 3) { Chris@202: lemin = ngpmin - 1 + lt; Chris@202: ieee = TRUE_; Chris@202: /* ( Non twos-complement machines, with gradual underflow; */ Chris@202: /* e.g., IEEE standard followers ) */ Chris@202: } else { Chris@202: lemin = min(ngpmin,gpmin); Chris@202: /* ( A guess; no known machine ) */ Chris@202: iwarn = TRUE_; Chris@202: } Chris@202: Chris@202: } else if (ngpmin == gpmin && ngnmin == gnmin) { Chris@202: if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1) { Chris@202: lemin = max(ngpmin,ngnmin); Chris@202: /* ( Twos-complement machines, no gradual underflow; */ Chris@202: /* e.g., CYBER 205 ) */ Chris@202: } else { Chris@202: lemin = min(ngpmin,ngnmin); Chris@202: /* ( A guess; no known machine ) */ Chris@202: iwarn = TRUE_; Chris@202: } Chris@202: Chris@202: } else if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1 && gpmin == gnmin) Chris@202: { Chris@202: if (gpmin - min(ngpmin,ngnmin) == 3) { Chris@202: lemin = max(ngpmin,ngnmin) - 1 + lt; Chris@202: /* ( Twos-complement machines with gradual underflow; */ Chris@202: /* no known machine ) */ Chris@202: } else { Chris@202: lemin = min(ngpmin,ngnmin); Chris@202: /* ( A guess; no known machine ) */ Chris@202: iwarn = TRUE_; Chris@202: } Chris@202: Chris@202: } else { Chris@202: /* Computing MIN */ Chris@202: i__1 = min(ngpmin,ngnmin), i__1 = min(i__1,gpmin); Chris@202: lemin = min(i__1,gnmin); Chris@202: /* ( A guess; no known machine ) */ Chris@202: iwarn = TRUE_; Chris@202: } Chris@202: first = FALSE_; Chris@202: /* ** */ Chris@202: /* Comment out this if block if EMIN is ok */ Chris@202: if (iwarn) { Chris@202: first = TRUE_; Chris@202: s_wsfe(&io___58); Chris@202: do_fio(&c__1, (char *)&lemin, (ftnlen)sizeof(integer)); Chris@202: e_wsfe(); Chris@202: } Chris@202: /* ** */ Chris@202: Chris@202: /* Assume IEEE arithmetic if we found denormalised numbers above, */ Chris@202: /* or if arithmetic seems to round in the IEEE style, determined */ Chris@202: /* in routine DLAMC1. A true IEEE machine should have both things */ Chris@202: /* true; however, faulty machines may have one or the other. */ Chris@202: Chris@202: ieee = ieee || lieee1; Chris@202: Chris@202: /* Compute RMIN by successive division by BETA. We could compute */ Chris@202: /* RMIN as BASE**( EMIN - 1 ), but some machines underflow during */ Chris@202: /* this computation. */ Chris@202: Chris@202: lrmin = 1.; Chris@202: i__1 = 1 - lemin; Chris@202: for (i__ = 1; i__ <= i__1; ++i__) { Chris@202: d__1 = lrmin * rbase; Chris@202: lrmin = dlamc3_(&d__1, &zero); Chris@202: /* L30: */ Chris@202: } Chris@202: Chris@202: /* Finally, call DLAMC5 to compute EMAX and RMAX. */ Chris@202: Chris@202: dlamc5_(&lbeta, <, &lemin, &ieee, &lemax, &lrmax); Chris@202: } Chris@202: Chris@202: *beta = lbeta; Chris@202: *t = lt; Chris@202: *rnd = lrnd; Chris@202: *eps = leps; Chris@202: *emin = lemin; Chris@202: *rmin = lrmin; Chris@202: *emax = lemax; Chris@202: *rmax = lrmax; Chris@202: Chris@202: return 0; Chris@202: Chris@202: Chris@202: /* End of DLAMC2 */ Chris@202: Chris@202: } /* dlamc2_ */ Chris@202: Chris@202: Chris@202: /* *********************************************************************** */ Chris@202: Chris@202: doublereal dlamc3_(doublereal *a, doublereal *b) Chris@202: { Chris@202: /* System generated locals */ Chris@202: doublereal ret_val; Chris@202: Chris@202: Chris@202: /* -- LAPACK auxiliary routine (version 3.2) -- */ Chris@202: /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ Chris@202: /* November 2006 */ Chris@202: Chris@202: /* .. Scalar Arguments .. */ Chris@202: /* .. */ Chris@202: Chris@202: /* Purpose */ Chris@202: /* ======= */ Chris@202: Chris@202: /* DLAMC3 is intended to force A and B to be stored prior to doing */ Chris@202: /* the addition of A and B , for use in situations where optimizers */ Chris@202: /* might hold one of these in a register. */ Chris@202: Chris@202: /* Arguments */ Chris@202: /* ========= */ Chris@202: Chris@202: /* A (input) DOUBLE PRECISION */ Chris@202: /* B (input) DOUBLE PRECISION */ Chris@202: /* The values A and B. */ Chris@202: Chris@202: /* ===================================================================== */ Chris@202: Chris@202: /* .. Executable Statements .. */ Chris@202: Chris@202: ret_val = *a + *b; Chris@202: Chris@202: return ret_val; Chris@202: Chris@202: /* End of DLAMC3 */ Chris@202: Chris@202: } /* dlamc3_ */ Chris@202: Chris@202: Chris@202: /* *********************************************************************** */ Chris@202: Chris@202: /* Subroutine */ int dlamc4_(integer *emin, doublereal *start, integer *base) Chris@202: { Chris@202: /* System generated locals */ Chris@202: integer i__1; Chris@202: doublereal d__1; Chris@202: Chris@202: /* Local variables */ Chris@202: doublereal a; Chris@202: integer i__; Chris@202: doublereal b1, b2, c1, c2, d1, d2, one, zero, rbase; Chris@202: extern doublereal dlamc3_(doublereal *, doublereal *); Chris@202: Chris@202: Chris@202: /* -- LAPACK auxiliary routine (version 3.2) -- */ Chris@202: /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ Chris@202: /* November 2006 */ Chris@202: Chris@202: /* .. Scalar Arguments .. */ Chris@202: /* .. */ Chris@202: Chris@202: /* Purpose */ Chris@202: /* ======= */ Chris@202: Chris@202: /* DLAMC4 is a service routine for DLAMC2. */ Chris@202: Chris@202: /* Arguments */ Chris@202: /* ========= */ Chris@202: Chris@202: /* EMIN (output) INTEGER */ Chris@202: /* The minimum exponent before (gradual) underflow, computed by */ Chris@202: /* setting A = START and dividing by BASE until the previous A */ Chris@202: /* can not be recovered. */ Chris@202: Chris@202: /* START (input) DOUBLE PRECISION */ Chris@202: /* The starting point for determining EMIN. */ Chris@202: Chris@202: /* BASE (input) INTEGER */ Chris@202: /* The base of the machine. */ Chris@202: Chris@202: /* ===================================================================== */ Chris@202: Chris@202: /* .. Local Scalars .. */ Chris@202: /* .. */ Chris@202: /* .. External Functions .. */ Chris@202: /* .. */ Chris@202: /* .. Executable Statements .. */ Chris@202: Chris@202: a = *start; Chris@202: one = 1.; Chris@202: rbase = one / *base; Chris@202: zero = 0.; Chris@202: *emin = 1; Chris@202: d__1 = a * rbase; Chris@202: b1 = dlamc3_(&d__1, &zero); Chris@202: c1 = a; Chris@202: c2 = a; Chris@202: d1 = a; Chris@202: d2 = a; Chris@202: /* + WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. */ Chris@202: /* $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP */ Chris@202: L10: Chris@202: if (c1 == a && c2 == a && d1 == a && d2 == a) { Chris@202: --(*emin); Chris@202: a = b1; Chris@202: d__1 = a / *base; Chris@202: b1 = dlamc3_(&d__1, &zero); Chris@202: d__1 = b1 * *base; Chris@202: c1 = dlamc3_(&d__1, &zero); Chris@202: d1 = zero; Chris@202: i__1 = *base; Chris@202: for (i__ = 1; i__ <= i__1; ++i__) { Chris@202: d1 += b1; Chris@202: /* L20: */ Chris@202: } Chris@202: d__1 = a * rbase; Chris@202: b2 = dlamc3_(&d__1, &zero); Chris@202: d__1 = b2 / rbase; Chris@202: c2 = dlamc3_(&d__1, &zero); Chris@202: d2 = zero; Chris@202: i__1 = *base; Chris@202: for (i__ = 1; i__ <= i__1; ++i__) { Chris@202: d2 += b2; Chris@202: /* L30: */ Chris@202: } Chris@202: goto L10; Chris@202: } Chris@202: /* + END WHILE */ Chris@202: Chris@202: return 0; Chris@202: Chris@202: /* End of DLAMC4 */ Chris@202: Chris@202: } /* dlamc4_ */ Chris@202: Chris@202: Chris@202: /* *********************************************************************** */ Chris@202: Chris@202: /* Subroutine */ int dlamc5_(integer *beta, integer *p, integer *emin, Chris@202: logical *ieee, integer *emax, doublereal *rmax) Chris@202: { Chris@202: /* System generated locals */ Chris@202: integer i__1; Chris@202: doublereal d__1; Chris@202: Chris@202: /* Local variables */ Chris@202: integer i__; Chris@202: doublereal y, z__; Chris@202: integer try__, lexp; Chris@202: doublereal oldy; Chris@202: integer uexp, nbits; Chris@202: extern doublereal dlamc3_(doublereal *, doublereal *); Chris@202: doublereal recbas; Chris@202: integer exbits, expsum; Chris@202: Chris@202: Chris@202: /* -- LAPACK auxiliary routine (version 3.2) -- */ Chris@202: /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ Chris@202: /* November 2006 */ Chris@202: Chris@202: /* .. Scalar Arguments .. */ Chris@202: /* .. */ Chris@202: Chris@202: /* Purpose */ Chris@202: /* ======= */ Chris@202: Chris@202: /* DLAMC5 attempts to compute RMAX, the largest machine floating-point */ Chris@202: /* number, without overflow. It assumes that EMAX + abs(EMIN) sum */ Chris@202: /* approximately to a power of 2. It will fail on machines where this */ Chris@202: /* assumption does not hold, for example, the Cyber 205 (EMIN = -28625, */ Chris@202: /* EMAX = 28718). It will also fail if the value supplied for EMIN is */ Chris@202: /* too large (i.e. too close to zero), probably with overflow. */ Chris@202: Chris@202: /* Arguments */ Chris@202: /* ========= */ Chris@202: Chris@202: /* BETA (input) INTEGER */ Chris@202: /* The base of floating-point arithmetic. */ Chris@202: Chris@202: /* P (input) INTEGER */ Chris@202: /* The number of base BETA digits in the mantissa of a */ Chris@202: /* floating-point value. */ Chris@202: Chris@202: /* EMIN (input) INTEGER */ Chris@202: /* The minimum exponent before (gradual) underflow. */ Chris@202: Chris@202: /* IEEE (input) LOGICAL */ Chris@202: /* A logical flag specifying whether or not the arithmetic */ Chris@202: /* system is thought to comply with the IEEE standard. */ Chris@202: Chris@202: /* EMAX (output) INTEGER */ Chris@202: /* The largest exponent before overflow */ Chris@202: Chris@202: /* RMAX (output) DOUBLE PRECISION */ Chris@202: /* The largest machine floating-point number. */ Chris@202: Chris@202: /* ===================================================================== */ Chris@202: Chris@202: /* .. Parameters .. */ Chris@202: /* .. */ Chris@202: /* .. Local Scalars .. */ Chris@202: /* .. */ Chris@202: /* .. External Functions .. */ Chris@202: /* .. */ Chris@202: /* .. Intrinsic Functions .. */ Chris@202: /* .. */ Chris@202: /* .. Executable Statements .. */ Chris@202: Chris@202: /* First compute LEXP and UEXP, two powers of 2 that bound */ Chris@202: /* abs(EMIN). We then assume that EMAX + abs(EMIN) will sum */ Chris@202: /* approximately to the bound that is closest to abs(EMIN). */ Chris@202: /* (EMAX is the exponent of the required number RMAX). */ Chris@202: Chris@202: lexp = 1; Chris@202: exbits = 1; Chris@202: L10: Chris@202: try__ = lexp << 1; Chris@202: if (try__ <= -(*emin)) { Chris@202: lexp = try__; Chris@202: ++exbits; Chris@202: goto L10; Chris@202: } Chris@202: if (lexp == -(*emin)) { Chris@202: uexp = lexp; Chris@202: } else { Chris@202: uexp = try__; Chris@202: ++exbits; Chris@202: } Chris@202: Chris@202: /* Now -LEXP is less than or equal to EMIN, and -UEXP is greater */ Chris@202: /* than or equal to EMIN. EXBITS is the number of bits needed to */ Chris@202: /* store the exponent. */ Chris@202: Chris@202: if (uexp + *emin > -lexp - *emin) { Chris@202: expsum = lexp << 1; Chris@202: } else { Chris@202: expsum = uexp << 1; Chris@202: } Chris@202: Chris@202: /* EXPSUM is the exponent range, approximately equal to */ Chris@202: /* EMAX - EMIN + 1 . */ Chris@202: Chris@202: *emax = expsum + *emin - 1; Chris@202: nbits = exbits + 1 + *p; Chris@202: Chris@202: /* NBITS is the total number of bits needed to store a */ Chris@202: /* floating-point number. */ Chris@202: Chris@202: if (nbits % 2 == 1 && *beta == 2) { Chris@202: Chris@202: /* Either there are an odd number of bits used to store a */ Chris@202: /* floating-point number, which is unlikely, or some bits are */ Chris@202: /* not used in the representation of numbers, which is possible, */ Chris@202: /* (e.g. Cray machines) or the mantissa has an implicit bit, */ Chris@202: /* (e.g. IEEE machines, Dec Vax machines), which is perhaps the */ Chris@202: /* most likely. We have to assume the last alternative. */ Chris@202: /* If this is true, then we need to reduce EMAX by one because */ Chris@202: /* there must be some way of representing zero in an implicit-bit */ Chris@202: /* system. On machines like Cray, we are reducing EMAX by one */ Chris@202: /* unnecessarily. */ Chris@202: Chris@202: --(*emax); Chris@202: } Chris@202: Chris@202: if (*ieee) { Chris@202: Chris@202: /* Assume we are on an IEEE machine which reserves one exponent */ Chris@202: /* for infinity and NaN. */ Chris@202: Chris@202: --(*emax); Chris@202: } Chris@202: Chris@202: /* Now create RMAX, the largest machine number, which should */ Chris@202: /* be equal to (1.0 - BETA**(-P)) * BETA**EMAX . */ Chris@202: Chris@202: /* First compute 1.0 - BETA**(-P), being careful that the */ Chris@202: /* result is less than 1.0 . */ Chris@202: Chris@202: recbas = 1. / *beta; Chris@202: z__ = *beta - 1.; Chris@202: y = 0.; Chris@202: i__1 = *p; Chris@202: for (i__ = 1; i__ <= i__1; ++i__) { Chris@202: z__ *= recbas; Chris@202: if (y < 1.) { Chris@202: oldy = y; Chris@202: } Chris@202: y = dlamc3_(&y, &z__); Chris@202: /* L20: */ Chris@202: } Chris@202: if (y >= 1.) { Chris@202: y = oldy; Chris@202: } Chris@202: Chris@202: /* Now multiply by BETA**EMAX to get RMAX. */ Chris@202: Chris@202: i__1 = *emax; Chris@202: for (i__ = 1; i__ <= i__1; ++i__) { Chris@202: d__1 = y * *beta; Chris@202: y = dlamc3_(&d__1, &c_b32); Chris@202: /* L30: */ Chris@202: } Chris@202: Chris@202: *rmax = y; Chris@202: return 0; Chris@202: Chris@202: /* End of DLAMC5 */ Chris@202: Chris@202: } /* dlamc5_ */