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

1276 lines
34 KiB
C
Raw Permalink 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 double dpow_ui(double x, integer n) {
double 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 doublereal c_b32 = 0.;
/* > \brief \b DLAMCHF77 deprecated */
/* =========== DOCUMENTATION =========== */
/* Online html documentation available at */
/* http://www.netlib.org/lapack/explore-html/ */
/* Definition: */
/* =========== */
/* DOUBLE PRECISION FUNCTION DLAMCH( CMACH ) */
/* CHARACTER CMACH */
/* > \par Purpose: */
/* ============= */
/* > */
/* > \verbatim */
/* > */
/* > DLAMCHF77 determines double precision machine parameters. */
/* > \endverbatim */
/* Arguments: */
/* ========== */
/* > \param[in] CMACH */
/* > \verbatim */
/* > 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) */
/* > \endverbatim */
/* Authors: */
/* ======== */
/* > \author Univ. of Tennessee */
/* > \author Univ. of California Berkeley */
/* > \author Univ. of Colorado Denver */
/* > \author NAG Ltd. */
/* > \date April 2012 */
/* > \ingroup auxOTHERauxiliary */
/* ===================================================================== */
doublereal dlamch_(char *cmach)
{
/* Initialized data */
static logical first = TRUE_;
/* System generated locals */
integer i__1;
doublereal ret_val;
/* Local variables */
static doublereal base;
integer beta;
static doublereal emin, prec, emax;
integer imin, imax;
logical lrnd;
static doublereal rmin, rmax, t;
doublereal rmach;
extern logical lsame_(char *, char *);
doublereal small;
static doublereal sfmin;
extern /* Subroutine */ int dlamc2_(integer *, integer *, logical *,
doublereal *, integer *, doublereal *, integer *, doublereal *);
integer it;
static doublereal 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) {
dlamc2_(&beta, &it, &lrnd, &eps, &imin, &rmin, &imax, &rmax);
base = (doublereal) beta;
t = (doublereal) it;
if (lrnd) {
rnd = 1.;
i__1 = 1 - it;
eps = pow_di(&base, &i__1) / 2;
} else {
rnd = 0.;
i__1 = 1 - it;
eps = pow_di(&base, &i__1);
}
prec = eps * base;
emin = (doublereal) imin;
emax = (doublereal) imax;
sfmin = rmin;
small = 1. / 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.);
}
}
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 DLAMCH */
} /* dlamch_ */
/* *********************************************************************** */
/* > \brief \b DLAMC1 */
/* > \details */
/* > \b Purpose: */
/* > \verbatim */
/* > DLAMC1 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 dlamc1_(integer *beta, integer *t, logical *rnd, logical
*ieee1)
{
/* Initialized data */
static logical first = TRUE_;
/* System generated locals */
doublereal d__1, d__2;
/* Local variables */
static logical lrnd;
doublereal a, b, c__, f;
static integer lbeta;
doublereal savec;
extern doublereal dlamc3_(doublereal *, doublereal *);
static logical lieee1;
doublereal t1, t2;
static integer lt;
doublereal 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.;
/* 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.;
c__ = 1.;
/* + WHILE( C.EQ.ONE )LOOP */
L10:
if (c__ == one) {
a *= 2;
c__ = dlamc3_(&a, &one);
d__1 = -a;
c__ = dlamc3_(&c__, &d__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.;
c__ = dlamc3_(&a, &b);
/* + WHILE( C.EQ.A )LOOP */
L20:
if (c__ == a) {
b *= 2;
c__ = dlamc3_(&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__;
d__1 = -a;
c__ = dlamc3_(&c__, &d__1);
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 = (doublereal) lbeta;
d__1 = b / 2;
d__2 = -b / 100;
f = dlamc3_(&d__1, &d__2);
c__ = dlamc3_(&f, &a);
if (c__ == a) {
lrnd = TRUE_;
} else {
lrnd = FALSE_;
}
d__1 = b / 2;
d__2 = b / 100;
f = dlamc3_(&d__1, &d__2);
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. */
d__1 = b / 2;
t1 = dlamc3_(&d__1, &a);
d__1 = b / 2;
t2 = dlamc3_(&d__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.;
c__ = 1.;
/* + WHILE( C.EQ.ONE )LOOP */
L30:
if (c__ == one) {
++lt;
a *= lbeta;
c__ = dlamc3_(&a, &one);
d__1 = -a;
c__ = dlamc3_(&c__, &d__1);
goto L30;
}
/* + END WHILE */
}
*beta = lbeta;
*t = lt;
*rnd = lrnd;
*ieee1 = lieee1;
first = FALSE_;
return 0;
/* End of DLAMC1 */
} /* dlamc1_ */
/* *********************************************************************** */
/* > \brief \b DLAMC2 */
/* > \details */
/* > \b Purpose: */
/* > \verbatim */
/* > DLAMC2 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 dlamc2_(integer *beta, integer *t, logical *rnd,
doublereal *eps, integer *emin, doublereal *rmin, integer *emax,
doublereal *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 DLAM"
"C2,\002,/\002 otherwise supply EMIN explicitly.\002,/)";
/* System generated locals */
integer i__1;
doublereal d__1, d__2, d__3, d__4, d__5;
/* Local variables */
logical ieee;
doublereal half;
logical lrnd;
static doublereal leps;
doublereal zero, a, b, c__;
integer i__;
static integer lbeta;
doublereal rbase;
static integer lemin, lemax;
integer gnmin;
doublereal small;
integer gpmin;
doublereal third;
static doublereal lrmin, lrmax;
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 *);
static integer lt;
integer ngnmin, ngpmin;
doublereal 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.;
one = 1.;
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. */
dlamc1_(&lbeta, &lt, &lrnd, &lieee1);
/* Start to find EPS. */
b = (doublereal) lbeta;
i__1 = -lt;
a = pow_di(&b, &i__1);
leps = a;
/* Try some tricks to see whether or not this is the correct EPS. */
b = two / 3;
half = one / 2;
d__1 = -half;
sixth = dlamc3_(&b, &d__1);
third = dlamc3_(&sixth, &sixth);
d__1 = -half;
b = dlamc3_(&third, &d__1);
b = dlamc3_(&b, &sixth);
b = abs(b);
if (b < leps) {
b = leps;
}
leps = 1.;
/* + WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP */
L10:
if (leps > b && b > zero) {
leps = b;
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);
d__1 = -c__;
c__ = dlamc3_(&half, &d__1);
b = dlamc3_(&half, &c__);
d__1 = -b;
c__ = dlamc3_(&half, &d__1);
b = dlamc3_(&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__) {
d__1 = small * rbase;
small = dlamc3_(&d__1, &zero);
/* L20: */
}
a = dlamc3_(&one, &small);
dlamc4_(&ngpmin, &one, &lbeta);
d__1 = -one;
dlamc4_(&ngnmin, &d__1, &lbeta);
dlamc4_(&gpmin, &a, &lbeta);
d__1 = -a;
dlamc4_(&gnmin, &d__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 DLAMC1. 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.;
i__1 = 1 - lemin;
for (i__ = 1; i__ <= i__1; ++i__) {
d__1 = lrmin * rbase;
lrmin = dlamc3_(&d__1, &zero);
/* L30: */
}
/* Finally, call DLAMC5 to compute EMAX and RMAX. */
dlamc5_(&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 DLAMC2 */
} /* dlamc2_ */
/* *********************************************************************** */
/* > \brief \b DLAMC3 */
/* > \details */
/* > \b Purpose: */
/* > \verbatim */
/* > DLAMC3 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 */
doublereal dlamc3_(doublereal *a, doublereal *b)
{
/* System generated locals */
doublereal 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 DLAMC3 */
} /* dlamc3_ */
/* *********************************************************************** */
/* > \brief \b DLAMC4 */
/* > \details */
/* > \b Purpose: */
/* > \verbatim */
/* > DLAMC4 is a service routine for DLAMC2. */
/* > \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 dlamc4_(integer *emin, doublereal *start, integer *base)
{
/* System generated locals */
integer i__1;
doublereal d__1;
/* Local variables */
doublereal zero, a;
integer i__;
doublereal rbase, b1, b2, c1, c2, d1, d2;
extern doublereal dlamc3_(doublereal *, doublereal *);
doublereal one;
/* -- LAPACK auxiliary routine (version 3.7.0) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2010 */
/* ===================================================================== */
a = *start;
one = 1.;
rbase = one / *base;
zero = 0.;
*emin = 1;
d__1 = a * rbase;
b1 = dlamc3_(&d__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;
d__1 = a / *base;
b1 = dlamc3_(&d__1, &zero);
d__1 = b1 * *base;
c1 = dlamc3_(&d__1, &zero);
d1 = zero;
i__1 = *base;
for (i__ = 1; i__ <= i__1; ++i__) {
d1 += b1;
/* L20: */
}
d__1 = a * rbase;
b2 = dlamc3_(&d__1, &zero);
d__1 = b2 / rbase;
c2 = dlamc3_(&d__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 DLAMC4 */
} /* dlamc4_ */
/* *********************************************************************** */
/* > \brief \b DLAMC5 */
/* > \details */
/* > \b Purpose: */
/* > \verbatim */
/* > DLAMC5 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 dlamc5_(integer *beta, integer *p, integer *emin,
logical *ieee, integer *emax, doublereal *rmax)
{
/* System generated locals */
integer i__1;
doublereal d__1;
/* Local variables */
integer lexp;
doublereal oldy;
integer uexp, i__;
doublereal y, z__;
integer nbits;
extern doublereal dlamc3_(doublereal *, doublereal *);
doublereal 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. / *beta;
z__ = *beta - 1.;
y = 0.;
i__1 = *p;
for (i__ = 1; i__ <= i__1; ++i__) {
z__ *= recbas;
if (y < 1.) {
oldy = y;
}
y = dlamc3_(&y, &z__);
/* L20: */
}
if (y >= 1.) {
y = oldy;
}
/* Now multiply by BETA**EMAX to get RMAX. */
i__1 = *emax;
for (i__ = 1; i__ <= i__1; ++i__) {
d__1 = y * *beta;
y = dlamc3_(&d__1, &c_b32);
/* L30: */
}
*rmax = y;
return 0;
/* End of DLAMC5 */
} /* dlamc5_ */
Reference in New Issue View Git Blame Copy Permalink