📄 dlamch.c
字号:
/* blas/dlamch.f -- translated by f2c (version 20050501).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
/* This code expects correct IEEE rounding behaviour which is not
always provided. The source should be built with -ffloat-store.
A note from the GCC man page:
-ffloat-store
Do not store floating point variables in registers. This pre-
vents undesirable excess precision on machines such as the 68000
where the floating registers (of the 68881) keep more precision
than a double is supposed to have.
For most programs, the excess precision does only good, but a
few programs rely on the precise definition of IEEE floating
point. Use `-ffloat-store' for such programs. */
/* Disable "global optimizations" to avoid optimizer bugs in this file.
For GCC the file should be compiled with the -fno-gcse option. Here
is a note from the GCC man page:
Note: When compiling a program using computed gotos, a GCC exten-
sion, you may get better runtime performance if you disable the
global common subexpression elimination pass by adding -fno-gcse to
the command line.
*/
#ifdef _MSC_VER
# pragma optimize ("g",off)
#endif
#ifdef __cplusplus
extern "C" {
#endif
#include "v3p_netlib.h"
#include <stdio.h>
/* Initialization function just calls the function once so that its
runtime-initialized constants are initialized. After the first
call it is safe to call the function from multiple threads at
once. */
void v3p_netlib_dlamch_init()
{
dlamch_(" ", 1);
}
/* Table of constant values */
static doublereal c_b32 = 0.;
/*< DOUBLE PRECISION FUNCTION DLAMCH( CMACH ) >*/
doublereal dlamch_(char *cmach, ftnlen cmach_len)
{
/* Initialized data */
static logical first = TRUE_; /* runtime-initialized constant */
/* System generated locals */
integer i__1;
doublereal ret_val;
/* Builtin functions */
double pow_di(doublereal *, integer *);
/* Local variables */
static doublereal t; /* runtime-initialized constant */
integer it;
static doublereal rnd, eps, base; /* runtime-initialized constant */
integer beta;
static doublereal emin, prec, emax; /* runtime-initialized constant */
integer imin, imax;
logical lrnd;
static doublereal rmin, rmax; /* runtime-initialized constant */
doublereal rmach=0;
extern logical lsame_(char *, char *, ftnlen, ftnlen);
doublereal small;
static doublereal sfmin; /* runtime-initialized constant */
extern /* Subroutine */ int dlamc2_(integer *, integer *, logical *,
doublereal *, integer *, doublereal *, integer *, doublereal *);
(void)cmach_len;
/* -- 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 .. */
/*< CHARACTER CMACH >*/
/* .. */
/* Purpose */
/* ======= */
/* DLAMCH determines double precision machine parameters. */
/* Arguments */
/* ========= */
/* CMACH (input) CHARACTER*1 */
/* Specifies the value to be returned by DLAMCH: */
/* = 'E' or 'e', DLAMCH := eps */
/* = 'S' or 's , DLAMCH := sfmin */
/* = 'B' or 'b', DLAMCH := base */
/* = 'P' or 'p', DLAMCH := eps*base */
/* = 'N' or 'n', DLAMCH := t */
/* = 'R' or 'r', DLAMCH := rnd */
/* = 'M' or 'm', DLAMCH := emin */
/* = 'U' or 'u', DLAMCH := rmin */
/* = 'L' or 'l', DLAMCH := emax */
/* = 'O' or 'o', DLAMCH := rmax */
/* where */
/* eps = relative machine precision */
/* sfmin = safe minimum, such that 1/sfmin does not overflow */
/* base = base of the machine */
/* prec = eps*base */
/* t = number of (base) digits in the mantissa */
/* rnd = 1.0 when rounding occurs in addition, 0.0 otherwise */
/* emin = minimum exponent before (gradual) underflow */
/* rmin = underflow threshold - base**(emin-1) */
/* emax = largest exponent before overflow */
/* rmax = overflow threshold - (base**emax)*(1-eps) */
/* ===================================================================== */
/* .. Parameters .. */
/*< DOUBLE PRECISION ONE, ZERO >*/
/*< PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) >*/
/* .. */
/* .. Local Scalars .. */
/*< LOGICAL FIRST, LRND >*/
/*< INTEGER BETA, IMAX, IMIN, IT >*/
/*< >*/
/* .. */
/* .. External Functions .. */
/*< LOGICAL LSAME >*/
/*< EXTERNAL LSAME >*/
/* .. */
/* .. External Subroutines .. */
/*< EXTERNAL DLAMC2 >*/
/* .. */
/* .. Save statement .. */
/*< >*/
/* .. */
/* .. Data statements .. */
/*< DATA FIRST / .TRUE. / >*/
/* .. */
/* .. Executable Statements .. */
/*< IF( FIRST ) THEN >*/
if (first) {
/*< FIRST = .FALSE. >*/
first = FALSE_;
/*< CALL DLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX ) >*/
dlamc2_(&beta, &it, &lrnd, &eps, &imin, &rmin, &imax, &rmax);
/*< BASE = BETA >*/
base = (doublereal) beta;
/*< T = IT >*/
t = (doublereal) it;
/*< IF( LRND ) THEN >*/
if (lrnd) {
/*< RND = ONE >*/
rnd = 1.;
/*< EPS = ( BASE**( 1-IT ) ) / 2 >*/
i__1 = 1 - it;
eps = pow_di(&base, &i__1) / 2;
/*< ELSE >*/
} else {
/*< RND = ZERO >*/
rnd = 0.;
/*< EPS = BASE**( 1-IT ) >*/
i__1 = 1 - it;
eps = pow_di(&base, &i__1);
/*< END IF >*/
}
/*< PREC = EPS*BASE >*/
prec = eps * base;
/*< EMIN = IMIN >*/
emin = (doublereal) imin;
/*< EMAX = IMAX >*/
emax = (doublereal) imax;
/*< SFMIN = RMIN >*/
sfmin = rmin;
/*< SMALL = ONE / RMAX >*/
small = 1. / rmax;
/*< IF( SMALL.GE.SFMIN ) THEN >*/
if (small >= sfmin) {
/* Use SMALL plus a bit, to avoid the possibility of rounding */
/* causing overflow when computing 1/sfmin. */
/*< SFMIN = SMALL*( ONE+EPS ) >*/
sfmin = small * (eps + 1.);
/*< END IF >*/
}
/*< END IF >*/
}
/*< IF( LSAME( CMACH, 'E' ) ) THEN >*/
if (lsame_(cmach, "E", (ftnlen)1, (ftnlen)1)) {
/*< RMACH = EPS >*/
rmach = eps;
/*< ELSE IF( LSAME( CMACH, 'S' ) ) THEN >*/
} else if (lsame_(cmach, "S", (ftnlen)1, (ftnlen)1)) {
/*< RMACH = SFMIN >*/
rmach = sfmin;
/*< ELSE IF( LSAME( CMACH, 'B' ) ) THEN >*/
} else if (lsame_(cmach, "B", (ftnlen)1, (ftnlen)1)) {
/*< RMACH = BASE >*/
rmach = base;
/*< ELSE IF( LSAME( CMACH, 'P' ) ) THEN >*/
} else if (lsame_(cmach, "P", (ftnlen)1, (ftnlen)1)) {
/*< RMACH = PREC >*/
rmach = prec;
/*< ELSE IF( LSAME( CMACH, 'N' ) ) THEN >*/
} else if (lsame_(cmach, "N", (ftnlen)1, (ftnlen)1)) {
/*< RMACH = T >*/
rmach = t;
/*< ELSE IF( LSAME( CMACH, 'R' ) ) THEN >*/
} else if (lsame_(cmach, "R", (ftnlen)1, (ftnlen)1)) {
/*< RMACH = RND >*/
rmach = rnd;
/*< ELSE IF( LSAME( CMACH, 'M' ) ) THEN >*/
} else if (lsame_(cmach, "M", (ftnlen)1, (ftnlen)1)) {
/*< RMACH = EMIN >*/
rmach = emin;
/*< ELSE IF( LSAME( CMACH, 'U' ) ) THEN >*/
} else if (lsame_(cmach, "U", (ftnlen)1, (ftnlen)1)) {
/*< RMACH = RMIN >*/
rmach = rmin;
/*< ELSE IF( LSAME( CMACH, 'L' ) ) THEN >*/
} else if (lsame_(cmach, "L", (ftnlen)1, (ftnlen)1)) {
/*< RMACH = EMAX >*/
rmach = emax;
/*< ELSE IF( LSAME( CMACH, 'O' ) ) THEN >*/
} else if (lsame_(cmach, "O", (ftnlen)1, (ftnlen)1)) {
/*< RMACH = RMAX >*/
rmach = rmax;
/*< END IF >*/
}
/*< DLAMCH = RMACH >*/
ret_val = rmach;
/*< RETURN >*/
return ret_val;
/* End of DLAMCH */
/*< END >*/
} /* dlamch_ */
/* *********************************************************************** */
/*< SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 ) >*/
/* Subroutine */ int dlamc1_(integer *beta, integer *t, logical *rnd, logical
*ieee1)
{
/* Initialized data */
static logical first = TRUE_; /* runtime-initialized constant */
/* System generated locals */
doublereal d__1, d__2;
/* Local variables */
doublereal a, b, c__, f, t1, t2;
static integer lt; /* runtime-initialized constant */
doublereal one, qtr;
static logical lrnd; /* runtime-initialized constant */
static integer lbeta; /* runtime-initialized constant */
doublereal savec;
extern doublereal dlamc3_(doublereal *, doublereal *);
static logical lieee1; /* runtime-initialized constant */
/* -- 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 IEEE1, RND >*/
/*< INTEGER BETA, T >*/
/* .. */
/* Purpose */
/* ======= */
/* DLAMC1 determines the machine parameters given by BETA, T, RND, and */
/* IEEE1. */
/* Arguments */
/* ========= */
/* BETA (output) INTEGER */
/* The base of the machine. */
/* T (output) INTEGER */
/* The number of ( BETA ) digits in the mantissa. */
/* RND (output) LOGICAL */
/* Specifies whether proper rounding ( RND = .TRUE. ) or */
/* chopping ( RND = .FALSE. ) occurs in addition. This may not */
/* be a reliable guide to the way in which the machine performs */
/* its arithmetic. */
/* IEEE1 (output) LOGICAL */
/* Specifies whether rounding appears to be done in the IEEE */
/* 'round to nearest' style. */
/* Further Details */
/* =============== */
/* The routine is based on the routine ENVRON by Malcolm and */
/* incorporates suggestions by Gentleman and Marovich. See */
/* Malcolm M. A. (1972) Algorithms to reveal properties of */
/* floating-point arithmetic. Comms. of the ACM, 15, 949-951. */
/* Gentleman W. M. and Marovich S. B. (1974) More on algorithms */
/* that reveal properties of floating point arithmetic units. */
/* Comms. of the ACM, 17, 276-277. */
/* ===================================================================== */
/* .. Local Scalars .. */
/*< LOGICAL FIRST, LIEEE1, LRND >*/
/*< INTEGER LBETA, LT >*/
/*< DOUBLE PRECISION A, B, C, F, ONE, QTR, SAVEC, T1, T2 >*/
/* .. */
/* .. External Functions .. */
/*< DOUBLE PRECISION DLAMC3 >*/
/*< EXTERNAL DLAMC3 >*/
/* .. */
/* .. Save statement .. */
/*< SAVE FIRST, LIEEE1, LBETA, LRND, LT >*/
/* .. */
/* .. Data statements .. */
/*< DATA FIRST / .TRUE. / >*/
/* .. */
/* .. Executable Statements .. */
/*< IF( FIRST ) THEN >*/
if (first) {
/*< FIRST = .FALSE. >*/
first = FALSE_;
/*< ONE = 1 >*/
one = 1.;
/* LBETA, LIEEE1, LT and LRND are the local values of BETA, */
/* IEEE1, T and RND. */
/* Throughout this routine we use the function DLAMC3 to ensure */
/* that relevant values are stored and not held in registers, or */
/* are not affected by optimizers. */
/* Compute a = 2.0**m with the smallest positive integer m such */
/* that */
/* fl( a + 1.0 ) = a. */
/*< A = 1 >*/
a = 1.;
/*< C = 1 >*/
c__ = 1.;
/* + WHILE( C.EQ.ONE )LOOP */
/*< 10 CONTINUE >*/
L10:
/*< IF( C.EQ.ONE ) THEN >*/
if (c__ == one) {
/*< A = 2*A >*/
a *= 2;
/*< C = DLAMC3( A, ONE ) >*/
c__ = dlamc3_(&a, &one);
/*< C = DLAMC3( C, -A ) >*/
d__1 = -a;
c__ = dlamc3_(&c__, &d__1);
/*< GO TO 10 >*/
goto L10;
/*< END IF >*/
}
/* + END WHILE */
/* Now compute b = 2.0**m with the smallest positive integer m */
/* such that */
/* fl( a + b ) .gt. a. */
/*< B = 1 >*/
b = 1.;
/*< C = DLAMC3( A, B ) >*/
c__ = dlamc3_(&a, &b);
/* + WHILE( C.EQ.A )LOOP */
/*< 20 CONTINUE >*/
L20:
/*< IF( C.EQ.A ) THEN >*/
if (c__ == a) {
/*< B = 2*B >*/
b *= 2;
/*< C = DLAMC3( A, B ) >*/
c__ = dlamc3_(&a, &b);
/*< GO TO 20 >*/
goto L20;
/*< END IF >*/
}
/* + END WHILE */
/* Now compute the base. a and c are neighbouring floating point */
/* numbers in the interval ( beta**t, beta**( t + 1 ) ) and so */
/* their difference is beta. Adding 0.25 to c is to ensure that it */
/* is truncated to beta and not ( beta - 1 ). */
/*< QTR = ONE / 4 >*/
qtr = one / 4;
/*< SAVEC = C >*/
savec = c__;
/*< C = DLAMC3( C, -A ) >*/
d__1 = -a;
c__ = dlamc3_(&c__, &d__1);
/*< LBETA = C + QTR >*/
lbeta = (integer) (c__ + qtr);
/* Now determine whether rounding or chopping occurs, by adding a */
/* bit less than beta/2 and a bit more than beta/2 to a. */
/*< B = LBETA >*/
b = (doublereal) lbeta;
/*< F = DLAMC3( B / 2, -B / 100 ) >*/
d__1 = b / 2;
d__2 = -b / 100;
f = dlamc3_(&d__1, &d__2);
/*< C = DLAMC3( F, A ) >*/
c__ = dlamc3_(&f, &a);
/*< IF( C.EQ.A ) THEN >*/
if (c__ == a) {
/*< LRND = .TRUE. >*/
lrnd = TRUE_;
/*< ELSE >*/
} else {
/*< LRND = .FALSE. >*/
lrnd = FALSE_;
/*< END IF >*/
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -