📄 dlamch.c
字号:
}
/*< 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 (lrnd && c__ == a) {
lrnd = FALSE_;
}
/* Try and decide whether rounding is done in the IEEE 'round to */
/* nearest' style. B/2 is half a unit in the last place of the two */
/* numbers A and SAVEC. Furthermore, A is even, i.e. has last bit */
/* zero, and SAVEC is odd. Thus adding B/2 to A should not change */
/* A, but adding B/2 to SAVEC should change SAVEC. */
/*< T1 = DLAMC3( B / 2, A ) >*/
d__1 = b / 2;
t1 = dlamc3_(&d__1, &a);
/*< T2 = DLAMC3( B / 2, SAVEC ) >*/
d__1 = b / 2;
t2 = dlamc3_(&d__1, &savec);
/*< LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND >*/
lieee1 = t1 == a && t2 > savec && lrnd;
/* Now find the mantissa, t. It should be the integer part of */
/* log to the base beta of a, however it is safer to determine t */
/* by powering. So we find t as the smallest positive integer for */
/* which */
/* fl( beta**t + 1.0 ) = 1.0. */
/*< LT = 0 >*/
lt = 0;
/*< A = 1 >*/
a = 1.;
/*< C = 1 >*/
c__ = 1.;
/* + WHILE( C.EQ.ONE )LOOP */
/*< 30 CONTINUE >*/
L30:
/*< IF( C.EQ.ONE ) THEN >*/
if (c__ == one) {
/*< LT = LT + 1 >*/
++lt;
/*< A = A*LBETA >*/
a *= lbeta;
/*< C = DLAMC3( A, ONE ) >*/
c__ = dlamc3_(&a, &one);
/*< C = DLAMC3( C, -A ) >*/
d__1 = -a;
c__ = dlamc3_(&c__, &d__1);
/*< GO TO 30 >*/
goto L30;
/*< END IF >*/
}
/* + END WHILE */
/*< END IF >*/
}
/*< BETA = LBETA >*/
*beta = lbeta;
/*< T = LT >*/
*t = lt;
/*< RND = LRND >*/
*rnd = lrnd;
/*< IEEE1 = LIEEE1 >*/
*ieee1 = lieee1;
/*< RETURN >*/
return 0;
/* End of DLAMC1 */
/*< END >*/
} /* dlamc1_ */
/* *********************************************************************** */
/*< SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX ) >*/
/* Subroutine */ int dlamc2_(integer *beta, integer *t, logical *rnd,
doublereal *eps, integer *emin, doublereal *rmin, integer *emax,
doublereal *rmax)
{
/* Initialized data */
static logical first = TRUE_; /* runtime-initialized constant */
static logical iwarn = FALSE_; /* runtime-initialized constant */
/* System generated locals */
integer i__1;
doublereal d__1, d__2, d__3, d__4, d__5;
/* Builtin functions */
double pow_di(doublereal *, integer *);
integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe();
/* Local variables */
doublereal a, b, c__;
integer i__;
static integer lt; /* runtime-initialized constant */
doublereal one, two;
logical ieee;
doublereal half;
logical lrnd;
static doublereal leps; /* runtime-initialized constant */
doublereal zero;
static integer lbeta; /* runtime-initialized constant */
doublereal rbase;
static integer lemin, lemax; /* runtime-initialized constant */
integer gnmin;
doublereal small;
integer gpmin;
doublereal third;
static doublereal lrmin, lrmax; /* runtime-initialized constant */
doublereal sixth;
extern /* Subroutine */ int dlamc1_(integer *, integer *, logical *,
logical *);
extern doublereal dlamc3_(doublereal *, doublereal *);
logical lieee1;
extern /* Subroutine */ int dlamc4_(integer *, doublereal *, integer *),
dlamc5_(integer *, integer *, integer *, logical *, integer *,
doublereal *);
integer ngnmin, ngpmin;
/* -- 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 RND >*/
/*< INTEGER BETA, EMAX, EMIN, T >*/
/*< DOUBLE PRECISION EPS, RMAX, RMIN >*/
/* .. */
/* Purpose */
/* ======= */
/* DLAMC2 determines the machine parameters specified in its argument */
/* list. */
/* 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. */
/* EPS (output) DOUBLE PRECISION */
/* The smallest positive number such that */
/* fl( 1.0 - EPS ) .LT. 1.0, */
/* where fl denotes the computed value. */
/* EMIN (output) INTEGER */
/* The minimum exponent before (gradual) underflow occurs. */
/* RMIN (output) DOUBLE PRECISION */
/* The smallest normalized number for the machine, given by */
/* BASE**( EMIN - 1 ), where BASE is the floating point value */
/* of BETA. */
/* EMAX (output) INTEGER */
/* The maximum exponent before overflow occurs. */
/* RMAX (output) DOUBLE PRECISION */
/* The largest positive number for the machine, given by */
/* BASE**EMAX * ( 1 - EPS ), where BASE is the floating point */
/* value of BETA. */
/* Further Details */
/* =============== */
/* The computation of EPS is based on a routine PARANOIA by */
/* W. Kahan of the University of California at Berkeley. */
/* ===================================================================== */
/* .. Local Scalars .. */
/*< LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND >*/
/*< >*/
/*< >*/
/* .. */
/* .. External Functions .. */
/*< DOUBLE PRECISION DLAMC3 >*/
/*< EXTERNAL DLAMC3 >*/
/* .. */
/* .. External Subroutines .. */
/*< EXTERNAL DLAMC1, DLAMC4, DLAMC5 >*/
/* .. */
/* .. Intrinsic Functions .. */
/*< INTRINSIC ABS, MAX, MIN >*/
/* .. */
/* .. Save statement .. */
/*< >*/
/* .. */
/* .. Data statements .. */
/*< DATA FIRST / .TRUE. / , IWARN / .FALSE. / >*/
/* .. */
/* .. Executable Statements .. */
/*< IF( FIRST ) THEN >*/
if (first) {
/*< FIRST = .FALSE. >*/
first = FALSE_;
/*< ZERO = 0 >*/
zero = 0.;
/*< ONE = 1 >*/
one = 1.;
/*< TWO = 2 >*/
two = 2.;
/* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of */
/* BETA, T, RND, EPS, EMIN and RMIN. */
/* 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. */
/* DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. */
/*< CALL DLAMC1( LBETA, LT, LRND, LIEEE1 ) >*/
dlamc1_(&lbeta, <, &lrnd, &lieee1);
/* Start to find EPS. */
/*< B = LBETA >*/
b = (doublereal) lbeta;
/*< A = B**( -LT ) >*/
i__1 = -lt;
a = pow_di(&b, &i__1);
/*< LEPS = A >*/
leps = a;
/* Try some tricks to see whether or not this is the correct EPS. */
/*< B = TWO / 3 >*/
b = two / 3;
/*< HALF = ONE / 2 >*/
half = one / 2;
/*< SIXTH = DLAMC3( B, -HALF ) >*/
d__1 = -half;
sixth = dlamc3_(&b, &d__1);
/*< THIRD = DLAMC3( SIXTH, SIXTH ) >*/
third = dlamc3_(&sixth, &sixth);
/*< B = DLAMC3( THIRD, -HALF ) >*/
d__1 = -half;
b = dlamc3_(&third, &d__1);
/*< B = DLAMC3( B, SIXTH ) >*/
b = dlamc3_(&b, &sixth);
/*< B = ABS( B ) >*/
b = abs(b);
/*< >*/
if (b < leps) {
b = leps;
}
/*< LEPS = 1 >*/
leps = 1.;
/* + WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP */
/*< 10 CONTINUE >*/
L10:
/*< IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN >*/
if (leps > b && b > zero) {
/*< LEPS = B >*/
leps = b;
/*< C = DLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) ) >*/
d__1 = half * leps;
/* Computing 5th power */
d__3 = two, d__4 = d__3, d__3 *= d__3;
/* Computing 2nd power */
d__5 = leps;
d__2 = d__4 * (d__3 * d__3) * (d__5 * d__5);
c__ = dlamc3_(&d__1, &d__2);
/*< C = DLAMC3( HALF, -C ) >*/
d__1 = -c__;
c__ = dlamc3_(&half, &d__1);
/*< B = DLAMC3( HALF, C ) >*/
b = dlamc3_(&half, &c__);
/*< C = DLAMC3( HALF, -B ) >*/
d__1 = -b;
c__ = dlamc3_(&half, &d__1);
/*< B = DLAMC3( HALF, C ) >*/
b = dlamc3_(&half, &c__);
/*< GO TO 10 >*/
goto L10;
/*< END IF >*/
}
/* + END WHILE */
/*< >*/
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 >*/
rbase = one / lbeta;
/*< SMALL = ONE >*/
small = one;
/*< DO 20 I = 1, 3 >*/
for (i__ = 1; i__ <= 3; ++i__) {
/*< SMALL = DLAMC3( SMALL*RBASE, ZERO ) >*/
d__1 = small * rbase;
small = dlamc3_(&d__1, &zero);
/*< 20 CONTINUE >*/
/* L20: */
}
/*< A = DLAMC3( ONE, SMALL ) >*/
a = dlamc3_(&one, &small);
/*< CALL DLAMC4( NGPMIN, ONE, LBETA ) >*/
dlamc4_(&ngpmin, &one, &lbeta);
/*< CALL DLAMC4( NGNMIN, -ONE, LBETA ) >*/
d__1 = -one;
dlamc4_(&ngnmin, &d__1, &lbeta);
/*< CALL DLAMC4( GPMIN, A, LBETA ) >*/
dlamc4_(&gpmin, &a, &lbeta);
/*< CALL DLAMC4( GNMIN, -A, LBETA ) >*/
d__1 = -a;
dlamc4_(&gnmin, &d__1, &lbeta);
/*< IEEE = .FALSE. >*/
ieee = FALSE_;
/*< IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN >*/
if (ngpmin == ngnmin && gpmin == gnmin) {
/*< IF( NGPMIN.EQ.GPMIN ) THEN >*/
if (ngpmin == gpmin) {
/*< LEMIN = NGPMIN >*/
lemin = ngpmin;
/* ( Non twos-complement machines, no gradual underflow; */
/* e.g., VAX ) */
/*< ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN >*/
} else if (gpmin - ngpmin == 3) {
/*< LEMIN = NGPMIN - 1 + LT >*/
lemin = ngpmin - 1 + lt;
/*< IEEE = .TRUE. >*/
ieee = TRUE_;
/* ( Non twos-complement machines, with gradual underflow; */
/* e.g., IEEE standard followers ) */
/*< ELSE >*/
} else {
/*< LEMIN = MIN( NGPMIN, GPMIN ) >*/
lemin = min(ngpmin,gpmin);
/* ( A guess; no known machine ) */
/*< IWARN = .TRUE. >*/
iwarn = TRUE_;
/*< END IF >*/
}
/*< ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN >*/
} else if (ngpmin == gpmin && ngnmin == gnmin) {
/*< IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN >*/
if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1) {
/*< LEMIN = MAX( NGPMIN, NGNMIN ) >*/
lemin = max(ngpmin,ngnmin);
/* ( Twos-complement machines, no gradual underflow; */
/* e.g., CYBER 205 ) */
/*< ELSE >*/
} else {
/*< LEMIN = MIN( NGPMIN, NGNMIN ) >*/
lemin = min(ngpmin,ngnmin);
/* ( A guess; no known machine ) */
/*< IWARN = .TRUE. >*/
iwarn = TRUE_;
/*< END IF >*/
}
/*< >*/
} else if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1 && gpmin == gnmin)
{
/*< IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN >*/
if (gpmin - min(ngpmin,ngnmin) == 3) {
/*< LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT >*/
lemin = max(ngpmin,ngnmin) - 1 + lt;
/* ( Twos-complement machines with gradual underflow; */
/* no known machine ) */
/*< ELSE >*/
} else {
/*< LEMIN = MIN( NGPMIN, NGNMIN ) >*/
lemin = min(ngpmin,ngnmin);
/* ( A guess; no known machine ) */
/*< IWARN = .TRUE. >*/
iwarn = TRUE_;
/*< END IF >*/
}
/*< ELSE >*/
} else {
/*< LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN ) >*/
/* Computing MIN */
i__1 = min(ngpmin,ngnmin), i__1 = min(i__1,gpmin);
lemin = min(i__1,gnmin);
/* ( A guess; no known machine ) */
/*< IWARN = .TRUE. >*/
iwarn = TRUE_;
/*< END IF >*/
}
/* ** */
/* Comment out this if block if EMIN is ok */
/*< IF( IWARN ) THEN >*/
/*< FIRST = .TRUE. >*/
/*< WRITE( 6, FMT = 9999 )LEMIN >*/
/*< END IF >*/
if (iwarn) {
first = TRUE_;
printf("\n\n WARNING. The value EMIN may be incorrect: - ");
printf("EMIN = %8ld\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 .OR. LIEEE1 >*/
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 = 1 >*/
lrmin = 1.;
/*< DO 30 I = 1, 1 - LEMIN >*/
i__1 = 1 - lemin;
for (i__ = 1; i__ <= i__1; ++i__) {
/*< LRMIN = DLAMC3( LRMIN*RBASE, ZERO ) >*/
d__1 = lrmin * rbase;
lrmin = dlamc3_(&d__1, &zero);
/*< 30 CONTINUE >*/
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -