📄 slamch.c
字号:
/*< END IF >*/
}
/*< F = SLAMC3( B / 2, B / 100 ) >*/
r__1 = b / 2;
r__2 = b / 100;
f = slamc3_(&r__1, &r__2);
/*< C = SLAMC3( F, A ) >*/
c__ = slamc3_(&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 = SLAMC3( B / 2, A ) >*/
r__1 = b / 2;
t1 = slamc3_(&r__1, &a);
/*< T2 = SLAMC3( B / 2, SAVEC ) >*/
r__1 = b / 2;
t2 = slamc3_(&r__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 = (float)1.;
/*< C = 1 >*/
c__ = (float)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 = SLAMC3( A, ONE ) >*/
c__ = slamc3_(&a, &one);
/*< C = SLAMC3( C, -A ) >*/
r__1 = -a;
c__ = slamc3_(&c__, &r__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 SLAMC1 */
/*< END >*/
} /* slamc1_ */
/* *********************************************************************** */
/*< SUBROUTINE SLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX ) >*/
/* Subroutine */ int slamc2_(integer *beta, integer *t, logical *rnd, real *
eps, integer *emin, real *rmin, integer *emax, real *rmax)
{
/* Initialized data */
static logical first = TRUE_; /* runtime-initialized constant */
static logical iwarn = FALSE_; /* runtime-initialized constant */
/* System generated locals */
integer i__1;
real r__1, r__2, r__3, r__4, r__5;
/* Builtin functions */
double pow_ri(real *, integer *);
/* Local variables */
real a, b, c__;
integer i__;
static integer lt; /* runtime-initialized constant */
real one, two;
logical ieee;
real half;
logical lrnd;
static real leps; /* runtime-initialized constant */
real zero;
static integer lbeta; /* runtime-initialized constant */
real rbase;
static integer lemin, lemax; /* runtime-initialized constant */
integer gnmin;
real small;
integer gpmin;
real third;
static real lrmin, lrmax; /* runtime-initialized constant */
real sixth;
logical lieee1;
extern /* Subroutine */ int slamc1_(integer *, integer *, logical *,
logical *);
extern doublereal slamc3_(real *, real *);
extern /* Subroutine */ int slamc4_(integer *, real *, integer *),
slamc5_(integer *, integer *, integer *, logical *, integer *,
real *);
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 >*/
/*< REAL EPS, RMAX, RMIN >*/
/* .. */
/* Purpose */
/* ======= */
/* SLAMC2 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) REAL */
/* 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) REAL */
/* 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) REAL */
/* 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 .. */
/*< REAL SLAMC3 >*/
/*< EXTERNAL SLAMC3 >*/
/* .. */
/* .. External Subroutines .. */
/*< EXTERNAL SLAMC1, SLAMC4, SLAMC5 >*/
/* .. */
/* .. 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 = (float)0.;
/*< ONE = 1 >*/
one = (float)1.;
/*< TWO = 2 >*/
two = (float)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 SLAMC3 to ensure */
/* 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. */
/*< CALL SLAMC1( LBETA, LT, LRND, LIEEE1 ) >*/
slamc1_(&lbeta, <, &lrnd, &lieee1);
/* Start to find EPS. */
/*< B = LBETA >*/
b = (real) lbeta;
/*< A = B**( -LT ) >*/
i__1 = -lt;
a = pow_ri(&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 = SLAMC3( B, -HALF ) >*/
r__1 = -half;
sixth = slamc3_(&b, &r__1);
/*< THIRD = SLAMC3( SIXTH, SIXTH ) >*/
third = slamc3_(&sixth, &sixth);
/*< B = SLAMC3( THIRD, -HALF ) >*/
r__1 = -half;
b = slamc3_(&third, &r__1);
/*< B = SLAMC3( B, SIXTH ) >*/
b = slamc3_(&b, &sixth);
/*< B = ABS( B ) >*/
b = dabs(b);
/*< >*/
if (b < leps) {
b = leps;
}
/*< LEPS = 1 >*/
leps = (float)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 = SLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) ) >*/
r__1 = half * leps;
/* Computing 5th power */
r__3 = two, r__4 = r__3, r__3 *= r__3;
/* Computing 2nd power */
r__5 = leps;
r__2 = r__4 * (r__3 * r__3) * (r__5 * r__5);
c__ = slamc3_(&r__1, &r__2);
/*< C = SLAMC3( HALF, -C ) >*/
r__1 = -c__;
c__ = slamc3_(&half, &r__1);
/*< B = SLAMC3( HALF, C ) >*/
b = slamc3_(&half, &c__);
/*< C = SLAMC3( HALF, -B ) >*/
r__1 = -b;
c__ = slamc3_(&half, &r__1);
/*< B = SLAMC3( HALF, C ) >*/
b = slamc3_(&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 = SLAMC3( SMALL*RBASE, ZERO ) >*/
r__1 = small * rbase;
small = slamc3_(&r__1, &zero);
/*< 20 CONTINUE >*/
/* L20: */
}
/*< A = SLAMC3( ONE, SMALL ) >*/
a = slamc3_(&one, &small);
/*< CALL SLAMC4( NGPMIN, ONE, LBETA ) >*/
slamc4_(&ngpmin, &one, &lbeta);
/*< CALL SLAMC4( NGNMIN, -ONE, LBETA ) >*/
r__1 = -one;
slamc4_(&ngnmin, &r__1, &lbeta);
/*< CALL SLAMC4( GPMIN, A, LBETA ) >*/
slamc4_(&gpmin, &a, &lbeta);
/*< CALL SLAMC4( GNMIN, -A, LBETA ) >*/
r__1 = -a;
slamc4_(&gnmin, &r__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 >*/
if (iwarn) {
/*< FIRST = .TRUE. >*/
first = TRUE_;
/*< WRITE( 6, FMT = 9999 )LEMIN >*/
printf("\n\n WARNING. The value EMIN may be incorrect: - ");
printf("EMIN = %8li\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");
/*< END IF >*/
}
/* ** */
/* 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 .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 = (float)1.;
/*< DO 30 I = 1, 1 - LEMIN >*/
i__1 = 1 - lemin;
for (i__ = 1; i__ <= i__1; ++i__) {
/*< LRMIN = SLAMC3( LRMIN*RBASE, ZERO ) >*/
r__1 = lrmin * rbase;
lrmin = slamc3_(&r__1, &zero);
/*< 30 CONTINUE >*/
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -