Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
78 changes: 29 additions & 49 deletions lapack-netlib/INSTALL/dlamch.c
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
#include <float.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
Expand Down Expand Up @@ -354,29 +355,13 @@ static doublereal c_b32 = 0.;
/* ===================================================================== */
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;
doublereal small, sfmin, eps;


/* -- LAPACK auxiliary routine (version 3.7.0) -- */
Expand All @@ -385,57 +370,52 @@ doublereal dlamch_(char *cmach)
/* 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;
/* 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.);
}
}

if (lsame_(cmach, "E")) {
rmach = eps;
} else if (lsame_(cmach, "S")) {
rmach = sfmin;
} else if (lsame_(cmach, "B")) {
rmach = base;
rmach = FLT_RADIX;
} else if (lsame_(cmach, "P")) {
rmach = prec;
rmach = eps * FLT_RADIX;
} else if (lsame_(cmach, "N")) {
rmach = t;
rmach = DBL_MANT_DIG;
} else if (lsame_(cmach, "R")) {
rmach = rnd;
rmach = 1.;
} else if (lsame_(cmach, "M")) {
rmach = emin;
rmach = DBL_MIN_EXP;
} else if (lsame_(cmach, "U")) {
rmach = rmin;
rmach = DBL_MIN;
} else if (lsame_(cmach, "L")) {
rmach = emax;
rmach = DBL_MAX_EXP;
} else if (lsame_(cmach, "O")) {
rmach = rmax;
rmach = DBL_MAX;
} else {
rmach = 0.;
}

ret_val = rmach;
first = FALSE_;
return ret_val;

/* End of DLAMCH */
Expand Down
78 changes: 29 additions & 49 deletions lapack-netlib/INSTALL/slamch.c
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
#include <float.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
Expand Down Expand Up @@ -353,29 +354,13 @@ static real c_b32 = 0.f;
/* ===================================================================== */
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;
real small, sfmin, eps;


/* -- LAPACK auxiliary routine (version 3.7.0) -- */
Expand All @@ -384,57 +369,52 @@ real slamch_(char *cmach)
/* 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;
/* The values below are those returned by the current slamch.f, which */
/* obtains them from the Fortran 90 inquiry intrinsics. They replace the */
/* slamc1/slamc2 probe of the deprecated slamchf77.f, which measures the */
/* format at run time and is only correct if intermediates are genuinely */
/* rounded to single. That does not hold on x87: the probe measures the */
/* 80 bit register format instead, and slamch 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 = FLT_EPSILON * 0.5f;

if (lsame_(cmach, "E")) {
rmach = eps;
} else if (lsame_(cmach, "S")) {
sfmin = FLT_MIN;
small = 1.f / FLT_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.f);
}
}

if (lsame_(cmach, "E")) {
rmach = eps;
} else if (lsame_(cmach, "S")) {
rmach = sfmin;
} else if (lsame_(cmach, "B")) {
rmach = base;
rmach = FLT_RADIX;
} else if (lsame_(cmach, "P")) {
rmach = prec;
rmach = eps * FLT_RADIX;
} else if (lsame_(cmach, "N")) {
rmach = t;
rmach = FLT_MANT_DIG;
} else if (lsame_(cmach, "R")) {
rmach = rnd;
rmach = 1.f;
} else if (lsame_(cmach, "M")) {
rmach = emin;
rmach = FLT_MIN_EXP;
} else if (lsame_(cmach, "U")) {
rmach = rmin;
rmach = FLT_MIN;
} else if (lsame_(cmach, "L")) {
rmach = emax;
rmach = FLT_MAX_EXP;
} else if (lsame_(cmach, "O")) {
rmach = rmax;
rmach = FLT_MAX;
} else {
rmach = 0.f;
}

ret_val = rmach;
first = FALSE_;
return ret_val;

/* End of SLAMCH */
Expand Down
Loading