phynexis-deps/OpenBLAS
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
Files
develop
OpenBLAS/lapack-netlib/INSTALL/dlamch.c
T
David Heffernan ffe8b242c7 C_LAPACK: take the float format from <float.h> in dlamch/slamch 
INSTALL/dlamch.c and INSTALL/slamch.c are f2c translations of the deprecated
dlamchf77.f and slamchf77.f, which determine the floating point format at run
time by probing in dlamc1/dlamc2 rather than reading it from the environment.
The current dlamch.f and slamch.f use the Fortran 90 inquiry intrinsics
(EPSILON, TINY, HUGE, DIGITS, MINEXPONENT, MAXEXPONENT, RADIX) instead, but f2c
cannot translate those, so the C LAPACK selected by NOFORTRAN=1 has been left
with the older probing implementation.

The probe is only correct if double intermediates are genuinely rounded to
double. That does not hold on x87. Building 32 bit for a target without SSE2,
so with -mfpmath=387, and with gcc 16, the intermediates stay in registers and
the probe measures the 80 bit register format: it reports emin/emax as
-16381/16384, and a mantissa width that follows the caller's x87 precision
control bits (64 at extended precision, 53 at double). Written back as doubles,
rmin underflows to 0 and rmax overflows to +Inf, so dlamch('S') and dlamch('U')
return 0 and dlamch('O') returns +Inf.

Everything that scales by those values is then wrong, mostly silently. The
first symptom to surface was a floating point exception rather than a wrong
answer: dsbevx computes

    safmin = dlamch('S')
    eps    = dlamch('P')
    smlnum = safmin / eps
    bignum = 1 / smlnum

so a zero safe minimum makes smlnum zero and the next line divides by zero.
Callers that unmask the divide by zero exception get a hard failure there;
callers that do not get whatever the wrong scaling produces.

Replace the probe with the <float.h> constants, mirroring the values the
current dlamch.f and slamch.f return. This fixes two lesser problems at the
same time: rmach was left uninitialised when cmach matched nothing, where
dlamch.f returns zero; and the cached static results made both routines
unsafe to call concurrently on first use.

dlamc1-dlamc5 and slamc1-slamc5 are left in place. They become unreachable
from dlamch/slamch, but dlamc3 and slamc3 have callers of their own in dlaed3,
dlaed9, dlals0, dlasd3, dlasd8 and their complex equivalents, where they serve
as optimiser barriers.
2026-08-17 10:55:05 +01:00

1256 lines
34 KiB
C
Raw Permalink Blame History

#include <float.h>
#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)
{
/* System generated locals */
doublereal ret_val;
/* Local variables */
doublereal rmach;
extern logical lsame_(char *, char *);
doublereal small, sfmin, 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 */
/* The values below are those returned by the current dlamch.f, which */
/* obtains them from the Fortran 90 inquiry intrinsics. They replace the */
/* dlamc1/dlamc2 probe of the deprecated dlamchf77.f, which measures the */
/* format at run time and is only correct if double intermediates are */
/* genuinely rounded to double. That does not hold on x87: the probe */
/* measures the 80 bit register format instead, and dlamch then returns 0 */
/* for the safe minimum and +Inf for the overflow threshold. Reading the */
/* constants also makes this routine thread safe, which the cached, */
/* probed version was not. */
eps = DBL_EPSILON * 0.5;
if (lsame_(cmach, "E")) {
rmach = eps;
} else if (lsame_(cmach, "S")) {
sfmin = DBL_MIN;
small = 1. / DBL_MAX;
if (small >= sfmin) {
/* Use SMALL plus a bit, to avoid the possibility of rounding */
/* causing overflow when computing 1/sfmin. */
sfmin = small * (eps + 1.);
}
rmach = sfmin;
} else if (lsame_(cmach, "B")) {
rmach = FLT_RADIX;
} else if (lsame_(cmach, "P")) {
rmach = eps * FLT_RADIX;
} else if (lsame_(cmach, "N")) {
rmach = DBL_MANT_DIG;
} else if (lsame_(cmach, "R")) {
rmach = 1.;
} else if (lsame_(cmach, "M")) {
rmach = DBL_MIN_EXP;
} else if (lsame_(cmach, "U")) {
rmach = DBL_MIN;
} else if (lsame_(cmach, "L")) {
rmach = DBL_MAX_EXP;
} else if (lsame_(cmach, "O")) {
rmach = DBL_MAX;
} else {
rmach = 0.;
}
ret_val = rmach;
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