floraachy
/
OpenBLAS
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
OpenBLAS/lapack-netlib/INSTALL/slamch.c

1274 lines
33 KiB
C
Raw Blame History

#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <complex.h>
#ifdef complex
#undef complex
#endif
#ifdef I
#undef I
#endif
#if defined(_WIN64)
typedef long long BLASLONG;
typedef unsigned long long BLASULONG;
#else
typedef long BLASLONG;
typedef unsigned long BLASULONG;
#endif
#ifdef LAPACK_ILP64
typedef BLASLONG blasint;
#if defined(_WIN64)
#define blasabs(x) llabs(x)
#else
#define blasabs(x) labs(x)
#endif
#else
typedef int blasint;
#define blasabs(x) abs(x)
#endif
typedef blasint integer;
typedef unsigned int uinteger;
typedef char *address;
typedef short int shortint;
typedef float real;
typedef double doublereal;
typedef struct { real r, i; } complex;
typedef struct { doublereal r, i; } doublecomplex;
#ifdef _MSC_VER
static inline _Fcomplex Cf(complex *z) {_Fcomplex zz={z->r , z->i}; return zz;}
static inline _Dcomplex Cd(doublecomplex *z) {_Dcomplex zz={z->r , z->i};return zz;}
static inline _Fcomplex * _pCf(complex *z) {return (_Fcomplex*)z;}
static inline _Dcomplex * _pCd(doublecomplex *z) {return (_Dcomplex*)z;}
#else
static inline _Complex float Cf(complex *z) {return z->r + z->i*_Complex_I;}
static inline _Complex double Cd(doublecomplex *z) {return z->r + z->i*_Complex_I;}
static inline _Complex float * _pCf(complex *z) {return (_Complex float*)z;}
static inline _Complex double * _pCd(doublecomplex *z) {return (_Complex double*)z;}
#endif
#define pCf(z) (*_pCf(z))
#define pCd(z) (*_pCd(z))
typedef int logical;
typedef short int shortlogical;
typedef char logical1;
typedef char integer1;
#define TRUE_ (1)
#define FALSE_ (0)
/* Extern is for use with -E */
#ifndef Extern
#define Extern extern
#endif
/* I/O stuff */
typedef int flag;
typedef int ftnlen;
typedef int ftnint;
/*external read, write*/
typedef struct
{ flag cierr;
ftnint ciunit;
flag ciend;
char *cifmt;
ftnint cirec;
} cilist;
/*internal read, write*/
typedef struct
{ flag icierr;
char *iciunit;
flag iciend;
char *icifmt;
ftnint icirlen;
ftnint icirnum;
} icilist;
/*open*/
typedef struct
{ flag oerr;
ftnint ounit;
char *ofnm;
ftnlen ofnmlen;
char *osta;
char *oacc;
char *ofm;
ftnint orl;
char *oblnk;
} olist;
/*close*/
typedef struct
{ flag cerr;
ftnint cunit;
char *csta;
} cllist;
/*rewind, backspace, endfile*/
typedef struct
{ flag aerr;
ftnint aunit;
} alist;
/* inquire */
typedef struct
{ flag inerr;
ftnint inunit;
char *infile;
ftnlen infilen;
ftnint *inex; /*parameters in standard's order*/
ftnint *inopen;
ftnint *innum;
ftnint *innamed;
char *inname;
ftnlen innamlen;
char *inacc;
ftnlen inacclen;
char *inseq;
ftnlen inseqlen;
char *indir;
ftnlen indirlen;
char *infmt;
ftnlen infmtlen;
char *inform;
ftnint informlen;
char *inunf;
ftnlen inunflen;
ftnint *inrecl;
ftnint *innrec;
char *inblank;
ftnlen inblanklen;
} inlist;
#define VOID void
union Multitype { /* for multiple entry points */
integer1 g;
shortint h;
integer i;
/* longint j; */
real r;
doublereal d;
complex c;
doublecomplex z;
};
typedef union Multitype Multitype;
struct Vardesc { /* for Namelist */
char *name;
char *addr;
ftnlen *dims;
int type;
};
typedef struct Vardesc Vardesc;
struct Namelist {
char *name;
Vardesc **vars;
int nvars;
};
typedef struct Namelist Namelist;
#define abs(x) ((x) >= 0 ? (x) : -(x))
#define dabs(x) (fabs(x))
#define f2cmin(a,b) ((a) <= (b) ? (a) : (b))
#define f2cmax(a,b) ((a) >= (b) ? (a) : (b))
#define dmin(a,b) (f2cmin(a,b))
#define dmax(a,b) (f2cmax(a,b))
#define bit_test(a,b) ((a) >> (b) & 1)
#define bit_clear(a,b) ((a) & ~((uinteger)1 << (b)))
#define bit_set(a,b) ((a) | ((uinteger)1 << (b)))
#define abort_() { sig_die("Fortran abort routine called", 1); }
#define c_abs(z) (cabsf(Cf(z)))
#define c_cos(R,Z) { pCf(R)=ccos(Cf(Z)); }
#ifdef _MSC_VER
#define c_div(c, a, b) {Cf(c)._Val[0] = (Cf(a)._Val[0]/Cf(b)._Val[0]); Cf(c)._Val[1]=(Cf(a)._Val[1]/Cf(b)._Val[1]);}
#define z_div(c, a, b) {Cd(c)._Val[0] = (Cd(a)._Val[0]/Cd(b)._Val[0]); Cd(c)._Val[1]=(Cd(a)._Val[1]/df(b)._Val[1]);}
#else
#define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);}
#define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);}
#endif
#define c_exp(R, Z) {pCf(R) = cexpf(Cf(Z));}
#define c_log(R, Z) {pCf(R) = clogf(Cf(Z));}
#define c_sin(R, Z) {pCf(R) = csinf(Cf(Z));}
//#define c_sqrt(R, Z) {*(R) = csqrtf(Cf(Z));}
#define c_sqrt(R, Z) {pCf(R) = csqrtf(Cf(Z));}
#define d_abs(x) (fabs(*(x)))
#define d_acos(x) (acos(*(x)))
#define d_asin(x) (asin(*(x)))
#define d_atan(x) (atan(*(x)))
#define d_atn2(x, y) (atan2(*(x),*(y)))
#define d_cnjg(R, Z) { pCd(R) = conj(Cd(Z)); }
#define r_cnjg(R, Z) { pCf(R) = conjf(Cf(Z)); }
#define d_cos(x) (cos(*(x)))
#define d_cosh(x) (cosh(*(x)))
#define d_dim(__a, __b) ( *(__a) > *(__b) ? *(__a) - *(__b) : 0.0 )
#define d_exp(x) (exp(*(x)))
#define d_imag(z) (cimag(Cd(z)))
#define r_imag(z) (cimagf(Cf(z)))
#define d_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
#define r_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
#define d_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
#define r_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
#define d_log(x) (log(*(x)))
#define d_mod(x, y) (fmod(*(x), *(y)))
#define u_nint(__x) ((__x)>=0 ? floor((__x) + .5) : -floor(.5 - (__x)))
#define d_nint(x) u_nint(*(x))
#define u_sign(__a,__b) ((__b) >= 0 ? ((__a) >= 0 ? (__a) : -(__a)) : -((__a) >= 0 ? (__a) : -(__a)))
#define d_sign(a,b) u_sign(*(a),*(b))
#define r_sign(a,b) u_sign(*(a),*(b))
#define d_sin(x) (sin(*(x)))
#define d_sinh(x) (sinh(*(x)))
#define d_sqrt(x) (sqrt(*(x)))
#define d_tan(x) (tan(*(x)))
#define d_tanh(x) (tanh(*(x)))
#define i_abs(x) abs(*(x))
#define i_dnnt(x) ((integer)u_nint(*(x)))
#define i_len(s, n) (n)
#define i_nint(x) ((integer)u_nint(*(x)))
#define i_sign(a,b) ((integer)u_sign((integer)*(a),(integer)*(b)))
#define pow_dd(ap, bp) ( pow(*(ap), *(bp)))
#define pow_si(B,E) spow_ui(*(B),*(E))
#define pow_ri(B,E) spow_ui(*(B),*(E))
#define pow_di(B,E) dpow_ui(*(B),*(E))
#define pow_zi(p, a, b) {pCd(p) = zpow_ui(Cd(a), *(b));}
#define pow_ci(p, a, b) {pCf(p) = cpow_ui(Cf(a), *(b));}
#define pow_zz(R,A,B) {pCd(R) = cpow(Cd(A),*(B));}
#define s_cat(lpp, rpp, rnp, np, llp) { ftnlen i, nc, ll; char *f__rp, *lp; ll = (llp); lp = (lpp); for(i=0; i < (int)*(np); ++i) { nc = ll; if((rnp)[i] < nc) nc = (rnp)[i]; ll -= nc; f__rp = (rpp)[i]; while(--nc >= 0) *lp++ = *(f__rp)++; } while(--ll >= 0) *lp++ = ' '; }
#define s_cmp(a,b,c,d) ((integer)strncmp((a),(b),f2cmin((c),(d))))
#define s_copy(A,B,C,D) { int __i,__m; for (__i=0, __m=f2cmin((C),(D)); __i<__m && (B)[__i] != 0; ++__i) (A)[__i] = (B)[__i]; }
#define sig_die(s, kill) { exit(1); }
#define s_stop(s, n) {exit(0);}
#define z_abs(z) (cabs(Cd(z)))
#define z_exp(R, Z) {pCd(R) = cexp(Cd(Z));}
#define z_sqrt(R, Z) {pCd(R) = csqrt(Cd(Z));}
#define myexit_() break;
#define mycycle() continue;
#define myceiling(w) {ceil(w)}
#define myhuge(w) {HUGE_VAL}
//#define mymaxloc_(w,s,e,n) {if (sizeof(*(w)) == sizeof(double)) dmaxloc_((w),*(s),*(e),n); else dmaxloc_((w),*(s),*(e),n);}
#define mymaxloc(w,s,e,n) {dmaxloc_(w,*(s),*(e),n)}
/* procedure parameter types for -A and -C++ */
#define F2C_proc_par_types 1
static float spow_ui(float x, integer n) {
float pow=1.0; unsigned long int u;
if(n != 0) {
if(n < 0) n = -n, x = 1/x;
for(u = n; ; ) {
if(u & 01) pow *= x;
if(u >>= 1) x *= x;
else break;
}
}
return pow;
}
/* -- translated by f2c (version 20000121).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Table of constant values */
static integer c__1 = 1;
static real c_b32 = 0.f;
/* > \brief \b SLAMCHF77 deprecated */
/* =========== DOCUMENTATION =========== */
/* Online html documentation available at */
/* http://www.netlib.org/lapack/explore-html/ */
/* Definition: */
/* =========== */
/* REAL FUNCTION SLAMCH( CMACH ) */
/* CHARACTER CMACH */
/* > \par Purpose: */
/* ============= */
/* > */
/* > \verbatim */
/* > */
/* > SLAMCH determines single precision machine parameters. */
/* > \endverbatim */
/* Arguments: */
/* ========== */
/* > \param[in] CMACH */
/* > \verbatim */
/* > Specifies the value to be returned by SLAMCH: */
/* > = 'E' or 'e', SLAMCH := eps */
/* > = 'S' or 's , SLAMCH := sfmin */
/* > = 'B' or 'b', SLAMCH := base */
/* > = 'P' or 'p', SLAMCH := eps*base */
/* > = 'N' or 'n', SLAMCH := t */
/* > = 'R' or 'r', SLAMCH := rnd */
/* > = 'M' or 'm', SLAMCH := emin */
/* > = 'U' or 'u', SLAMCH := rmin */
/* > = 'L' or 'l', SLAMCH := emax */
/* > = 'O' or 'o', SLAMCH := 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) */
/* > \endverbatim */
/* Authors: */
/* ======== */
/* > \author Univ. of Tennessee */
/* > \author Univ. of California Berkeley */
/* > \author Univ. of Colorado Denver */
/* > \author NAG Ltd. */
/* > \date April 2012 */
/* > \ingroup auxOTHERauxiliary */
/* ===================================================================== */
real slamch_(char *cmach)
{
/* Initialized data */
static logical first = TRUE_;
/* System generated locals */
integer i__1;
real ret_val;
/* Local variables */
static real base;
integer beta;
static real emin, prec, emax;
integer imin, imax;
logical lrnd;
static real rmin, rmax, t;
real rmach;
extern logical lsame_(char *, char *);
real small;
static real sfmin;
extern /* Subroutine */ int slamc2_(integer *, integer *, logical *, real
*, integer *, real *, integer *, real *);
integer it;
static real rnd, eps;
/* -- LAPACK auxiliary routine (version 3.7.0) -- */
/* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
/* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
/* April 2012 */
if (first) {
slamc2_(&beta, &it, &lrnd, &eps, &imin, &rmin, &imax, &rmax);
base = (real) beta;
t = (real) it;
if (lrnd) {
rnd = 1.f;
i__1 = 1 - it;
eps = pow_ri(&base, &i__1) / 2;
} else {
rnd = 0.f;
i__1 = 1 - it;
eps = pow_ri(&base, &i__1);
}
prec = eps * base;
emin = (real) imin;
emax = (real) imax;
sfmin = rmin;
small = 1.f / rmax;
if (small >= sfmin) {
/* Use SMALL plus a bit, to avoid the possibility of rounding */
/* causing overflow when computing 1/sfmin. */
sfmin = small * (eps + 1.f);
}
}
if (lsame_(cmach, "E")) {
rmach = eps;
} else if (lsame_(cmach, "S")) {
rmach = sfmin;
} else if (lsame_(cmach, "B")) {
rmach = base;
} else if (lsame_(cmach, "P")) {
rmach = prec;
} else if (lsame_(cmach, "N")) {
rmach = t;
} else if (lsame_(cmach, "R")) {
rmach = rnd;
} else if (lsame_(cmach, "M")) {
rmach = emin;
} else if (lsame_(cmach, "U")) {
rmach = rmin;
} else if (lsame_(cmach, "L")) {
rmach = emax;
} else if (lsame_(cmach, "O")) {
rmach = rmax;
}
ret_val = rmach;
first = FALSE_;
return ret_val;
/* End of SLAMCH */
} /* slamch_ */
/* *********************************************************************** */
/* > \brief \b SLAMC1 */
/* > \details */
/* > \b Purpose: */
/* > \verbatim */
/* > SLAMC1 determines the machine parameters given by BETA, T, RND, and */
/* > IEEE1. */
/* > \endverbatim */
/* > */
/* > \param[out] BETA */
/* > \verbatim */
/* > The base of the machine. */
/* > \endverbatim */
/* > */
/* > \param[out] T */
/* > \verbatim */
/* > The number of ( BETA ) digits in the mantissa. */
/* > \endverbatim */
/* > */
/* > \param[out] RND */
/* > \verbatim */
/* > 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. */
/* > \endverbatim */
/* > */
/* > \param[out] IEEE1 */
/* > \verbatim */
/* > Specifies whether rounding appears to be done in the IEEE */
/* > 'round to nearest' style. */
/* > \endverbatim */
/* > \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ.
of Colorado Denver and NAG Ltd.. */
/* > \date April 2012 */
/* > \ingroup auxOTHERauxiliary */
/* > */
/* > \details \b Further \b Details */
/* > \verbatim */
/* > */
/* > 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. */
/* > \endverbatim */
/* > */
/* Subroutine */ int slamc1_(integer *beta, integer *t, logical *rnd, logical
*ieee1)
{
/* Initialized data */
static logical first = TRUE_;
/* System generated locals */
real r__1, r__2;
/* Local variables */
static logical lrnd;
real a, b, c__, f;
static integer lbeta;
real savec;
static logical lieee1;
real t1, t2;
extern real slamc3_(real *, real *);
static integer lt;
real one, qtr;
/* -- LAPACK auxiliary routine (version 3.7.0) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2010 */
/* ===================================================================== */
if (first) {
one = 1.f;
/* LBETA, LIEEE1, LT and LRND are the local values of BETA, */
/* IEEE1, T and RND. */
/* 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. */
/* Compute a = 2.0**m with the smallest positive integer m such */
/* that */
/* fl( a + 1.0 ) = a. */
a = 1.f;
c__ = 1.f;
/* + WHILE( C.EQ.ONE )LOOP */
L10:
if (c__ == one) {
a *= 2;
c__ = slamc3_(&a, &one);
r__1 = -a;
c__ = slamc3_(&c__, &r__1);
goto L10;
}
/* + END WHILE */
/* Now compute b = 2.0**m with the smallest positive integer m */
/* such that */
/* fl( a + b ) .gt. a. */
b = 1.f;
c__ = slamc3_(&a, &b);
/* + WHILE( C.EQ.A )LOOP */
L20:
if (c__ == a) {
b *= 2;
c__ = slamc3_(&a, &b);
goto L20;
}
/* + 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;
savec = c__;
r__1 = -a;
c__ = slamc3_(&c__, &r__1);
lbeta = 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 = (real) lbeta;
r__1 = b / 2;
r__2 = -b / 100;
f = slamc3_(&r__1, &r__2);
c__ = slamc3_(&f, &a);
if (c__ == a) {
lrnd = TRUE_;
} else {
lrnd = FALSE_;
}
r__1 = b / 2;
r__2 = b / 100;
f = slamc3_(&r__1, &r__2);
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. */
r__1 = b / 2;
t1 = slamc3_(&r__1, &a);
r__1 = b / 2;
t2 = slamc3_(&r__1, &savec);
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;
a = 1.f;
c__ = 1.f;
/* + WHILE( C.EQ.ONE )LOOP */
L30:
if (c__ == one) {
++lt;
a *= lbeta;
c__ = slamc3_(&a, &one);
r__1 = -a;
c__ = slamc3_(&c__, &r__1);
goto L30;
}
/* + END WHILE */
}
*beta = lbeta;
*t = lt;
*rnd = lrnd;
*ieee1 = lieee1;
first = FALSE_;
return 0;
/* End of SLAMC1 */
} /* slamc1_ */
/* *********************************************************************** */
/* > \brief \b SLAMC2 */
/* > \details */
/* > \b Purpose: */
/* > \verbatim */
/* > SLAMC2 determines the machine parameters specified in its argument */
/* > list. */
/* > \endverbatim */
/* > \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ.
of Colorado Denver and NAG Ltd.. */
/* > \date April 2012 */
/* > \ingroup auxOTHERauxiliary */
/* > */
/* > \param[out] BETA */
/* > \verbatim */
/* > The base of the machine. */
/* > \endverbatim */
/* > */
/* > \param[out] T */
/* > \verbatim */
/* > The number of ( BETA ) digits in the mantissa. */
/* > \endverbatim */
/* > */
/* > \param[out] RND */
/* > \verbatim */
/* > 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. */
/* > \endverbatim */
/* > */
/* > \param[out] EPS */
/* > \verbatim */
/* > The smallest positive number such that */
/* > fl( 1.0 - EPS ) .LT. 1.0, */
/* > where fl denotes the computed value. */
/* > \endverbatim */
/* > */
/* > \param[out] EMIN */
/* > \verbatim */
/* > The minimum exponent before (gradual) underflow occurs. */
/* > \endverbatim */
/* > */
/* > \param[out] RMIN */
/* > \verbatim */
/* > The smallest normalized number for the machine, given by */
/* > BASE**( EMIN - 1 ), where BASE is the floating point value */
/* > of BETA. */
/* > \endverbatim */
/* > */
/* > \param[out] EMAX */
/* > \verbatim */
/* > The maximum exponent before overflow occurs. */
/* > \endverbatim */
/* > */
/* > \param[out] RMAX */
/* > \verbatim */
/* > The largest positive number for the machine, given by */
/* > BASE**EMAX * ( 1 - EPS ), where BASE is the floating point */
/* > value of BETA. */
/* > \endverbatim */
/* > */
/* > \details \b Further \b Details */
/* > \verbatim */
/* > */
/* > The computation of EPS is based on a routine PARANOIA by */
/* > W. Kahan of the University of California at Berkeley. */
/* > \endverbatim */
/* 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_;
static logical iwarn = FALSE_;
/* Format strings */
static char fmt_9999[] = "(//\002 WARNING. The value EMIN may be incorre"
"ct:-\002,\002 EMIN = \002,i8,/\002 If, after inspection, the va"
"lue EMIN looks\002,\002 acceptable please comment out \002,/\002"
" the IF block as marked within the code of routine\002,\002 SLAM"
"C2,\002,/\002 otherwise supply EMIN explicitly.\002,/)";
/* System generated locals */
integer i__1;
real r__1, r__2, r__3, r__4, r__5;
/* Local variables */
logical ieee;
real half;
logical lrnd;
static real leps;
real zero, a, b, c__;
integer i__;
static integer lbeta;
real rbase;
static integer lemin, lemax;
integer gnmin;
real small;
integer gpmin;
real third;
static real lrmin, lrmax;
real sixth;
logical lieee1;
extern /* Subroutine */ int slamc1_(integer *, integer *, logical *,
logical *);
extern real slamc3_(real *, real *);
extern /* Subroutine */ int slamc4_(integer *, real *, integer *),
slamc5_(integer *, integer *, integer *, logical *, integer *,
real *);
static integer lt;
integer ngnmin, ngpmin;
real one, two;
/* Fortran I/O blocks */
static cilist io___58 = { 0, 6, 0, fmt_9999, 0 };
/* -- LAPACK auxiliary routine (version 3.7.0) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2010 */
/* ===================================================================== */
if (first) {
zero = 0.f;
one = 1.f;
two = 2.f;
/* 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. */
slamc1_(&lbeta, &lt, &lrnd, &lieee1);
/* Start to find EPS. */
b = (real) lbeta;
i__1 = -lt;
a = pow_ri(&b, &i__1);
leps = a;
/* Try some tricks to see whether or not this is the correct EPS. */
b = two / 3;
half = one / 2;
r__1 = -half;
sixth = slamc3_(&b, &r__1);
third = slamc3_(&sixth, &sixth);
r__1 = -half;
b = slamc3_(&third, &r__1);
b = slamc3_(&b, &sixth);
b = abs(b);
if (b < leps) {
b = leps;
}
leps = 1.f;
/* + WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP */
L10:
if (leps > b && b > zero) {
leps = b;
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);
r__1 = -c__;
c__ = slamc3_(&half, &r__1);
b = slamc3_(&half, &c__);
r__1 = -b;
c__ = slamc3_(&half, &r__1);
b = slamc3_(&half, &c__);
goto L10;
}
/* + 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;
small = one;
for (i__ = 1; i__ <= 3; ++i__) {
r__1 = small * rbase;
small = slamc3_(&r__1, &zero);
/* L20: */
}
a = slamc3_(&one, &small);
slamc4_(&ngpmin, &one, &lbeta);
r__1 = -one;
slamc4_(&ngnmin, &r__1, &lbeta);
slamc4_(&gpmin, &a, &lbeta);
r__1 = -a;
slamc4_(&gnmin, &r__1, &lbeta);
ieee = FALSE_;
if (ngpmin == ngnmin && gpmin == gnmin) {
if (ngpmin == gpmin) {
lemin = ngpmin;
/* ( Non twos-complement machines, no gradual underflow; */
/* e.g., VAX ) */
} else if (gpmin - ngpmin == 3) {
lemin = ngpmin - 1 + lt;
ieee = TRUE_;
/* ( Non twos-complement machines, with gradual underflow; */
/* e.g., IEEE standard followers ) */
} else {
lemin = f2cmin(ngpmin,gpmin);
/* ( A guess; no known machine ) */
iwarn = TRUE_;
}
} else if (ngpmin == gpmin && ngnmin == gnmin) {
if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1) {
lemin = f2cmax(ngpmin,ngnmin);
/* ( Twos-complement machines, no gradual underflow; */
/* e.g., CYBER 205 ) */
} else {
lemin = f2cmin(ngpmin,ngnmin);
/* ( A guess; no known machine ) */
iwarn = TRUE_;
}
} else if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1 && gpmin == gnmin)
{
if (gpmin - f2cmin(ngpmin,ngnmin) == 3) {
lemin = f2cmax(ngpmin,ngnmin) - 1 + lt;
/* ( Twos-complement machines with gradual underflow; */
/* no known machine ) */
} else {
lemin = f2cmin(ngpmin,ngnmin);
/* ( A guess; no known machine ) */
iwarn = TRUE_;
}
} else {
/* Computing MIN */
i__1 = f2cmin(ngpmin,ngnmin), i__1 = f2cmin(i__1,gpmin);
lemin = f2cmin(i__1,gnmin);
/* ( A guess; no known machine ) */
iwarn = TRUE_;
}
first = FALSE_;
/* ** */
/* Comment out this if block if EMIN is ok */
/*
if (iwarn) {
first = TRUE_;
s_wsfe(&io___58);
do_fio(&c__1, (char *)&lemin, (ftnlen)sizeof(integer));
e_wsfe();
}
*/
/* ** */
/* 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 || 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.f;
i__1 = 1 - lemin;
for (i__ = 1; i__ <= i__1; ++i__) {
r__1 = lrmin * rbase;
lrmin = slamc3_(&r__1, &zero);
/* L30: */
}
/* Finally, call SLAMC5 to compute EMAX and RMAX. */
slamc5_(&lbeta, &lt, &lemin, &ieee, &lemax, &lrmax);
}
*beta = lbeta;
*t = lt;
*rnd = lrnd;
*eps = leps;
*emin = lemin;
*rmin = lrmin;
*emax = lemax;
*rmax = lrmax;
return 0;
/* End of SLAMC2 */
} /* slamc2_ */
/* *********************************************************************** */
/* > \brief \b SLAMC3 */
/* > \details */
/* > \b Purpose: */
/* > \verbatim */
/* > 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. */
/* > \endverbatim */
/* > */
/* > \param[in] A */
/* > */
/* > \param[in] B */
/* > \verbatim */
/* > The values A and B. */
/* > \endverbatim */
real slamc3_(real *a, real *b)
{
/* System generated locals */
real ret_val;
/* -- LAPACK auxiliary routine (version 3.7.0) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2010 */
/* ===================================================================== */
ret_val = *a + *b;
return ret_val;
/* End of SLAMC3 */
} /* slamc3_ */
/* *********************************************************************** */
/* > \brief \b SLAMC4 */
/* > \details */
/* > \b Purpose: */
/* > \verbatim */
/* > SLAMC4 is a service routine for SLAMC2. */
/* > \endverbatim */
/* > */
/* > \param[out] EMIN */
/* > \verbatim */
/* > The minimum exponent before (gradual) underflow, computed by */
/* > setting A = START and dividing by BASE until the previous A */
/* > can not be recovered. */
/* > \endverbatim */
/* > */
/* > \param[in] START */
/* > \verbatim */
/* > The starting point for determining EMIN. */
/* > \endverbatim */
/* > */
/* > \param[in] BASE */
/* > \verbatim */
/* > The base of the machine. */
/* > \endverbatim */
/* > */
/* Subroutine */ int slamc4_(integer *emin, real *start, integer *base)
{
/* System generated locals */
integer i__1;
real r__1;
/* Local variables */
real zero, a;
integer i__;
real rbase, b1, b2, c1, c2, d1, d2;
extern real slamc3_(real *, real *);
real one;
/* -- LAPACK auxiliary routine (version 3.7.0) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2010 */
/* ===================================================================== */
a = *start;
one = 1.f;
rbase = one / *base;
zero = 0.f;
*emin = 1;
r__1 = a * rbase;
b1 = slamc3_(&r__1, &zero);
c1 = a;
c2 = a;
d1 = a;
d2 = a;
/* + WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. */
/* $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP */
L10:
if (c1 == a && c2 == a && d1 == a && d2 == a) {
--(*emin);
a = b1;
r__1 = a / *base;
b1 = slamc3_(&r__1, &zero);
r__1 = b1 * *base;
c1 = slamc3_(&r__1, &zero);
d1 = zero;
i__1 = *base;
for (i__ = 1; i__ <= i__1; ++i__) {
d1 += b1;
/* L20: */
}
r__1 = a * rbase;
b2 = slamc3_(&r__1, &zero);
r__1 = b2 / rbase;
c2 = slamc3_(&r__1, &zero);
d2 = zero;
i__1 = *base;
for (i__ = 1; i__ <= i__1; ++i__) {
d2 += b2;
/* L30: */
}
goto L10;
}
/* + END WHILE */
return 0;
/* End of SLAMC4 */
} /* slamc4_ */
/* *********************************************************************** */
/* > \brief \b SLAMC5 */
/* > \details */
/* > \b Purpose: */
/* > \verbatim */
/* > 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. */
/* > \endverbatim */
/* > */
/* > \param[in] BETA */
/* > \verbatim */
/* > The base of floating-point arithmetic. */
/* > \endverbatim */
/* > */
/* > \param[in] P */
/* > \verbatim */
/* > The number of base BETA digits in the mantissa of a */
/* > floating-point value. */
/* > \endverbatim */
/* > */
/* > \param[in] EMIN */
/* > \verbatim */
/* > The minimum exponent before (gradual) underflow. */
/* > \endverbatim */
/* > */
/* > \param[in] IEEE */
/* > \verbatim */
/* > A logical flag specifying whether or not the arithmetic */
/* > system is thought to comply with the IEEE standard. */
/* > \endverbatim */
/* > */
/* > \param[out] EMAX */
/* > \verbatim */
/* > The largest exponent before overflow */
/* > \endverbatim */
/* > */
/* > \param[out] RMAX */
/* > \verbatim */
/* > The largest machine floating-point number. */
/* > \endverbatim */
/* > */
/* 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 lexp;
real oldy;
integer uexp, i__;
real y, z__;
integer nbits;
extern real slamc3_(real *, real *);
real recbas;
integer exbits, expsum, try__;
/* -- LAPACK auxiliary routine (version 3.7.0) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2010 */
/* ===================================================================== */
/* 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;
exbits = 1;
L10:
try__ = lexp << 1;
if (try__ <= -(*emin)) {
lexp = try__;
++exbits;
goto L10;
}
if (lexp == -(*emin)) {
uexp = lexp;
} else {
uexp = try__;
++exbits;
}
/* 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 > -lexp - *emin) {
expsum = lexp << 1;
} else {
expsum = uexp << 1;
}
/* EXPSUM is the exponent range, approximately equal to */
/* EMAX - EMIN + 1 . */
*emax = expsum + *emin - 1;
nbits = exbits + 1 + *p;
/* NBITS is the total number of bits needed to store a */
/* floating-point number. */
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);
}
if (*ieee) {
/* Assume we are on an IEEE machine which reserves one exponent */
/* for infinity and NaN. */
--(*emax);
}
/* 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 = 1.f / *beta;
z__ = *beta - 1.f;
y = 0.f;
i__1 = *p;
for (i__ = 1; i__ <= i__1; ++i__) {
z__ *= recbas;
if (y < 1.f) {
oldy = y;
}
y = slamc3_(&y, &z__);
/* L20: */
}
if (y >= 1.f) {
y = oldy;
}
/* Now multiply by BETA**EMAX to get RMAX. */
i__1 = *emax;
for (i__ = 1; i__ <= i__1; ++i__) {
r__1 = y * *beta;
y = slamc3_(&r__1, &c_b32);
/* L30: */
}
*rmax = y;
return 0;
/* End of SLAMC5 */
} /* slamc5_ */
Reference in New Issue View Git Blame Copy Permalink