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.

Author: davidheff

lapack-netlib/INSTALL/dlamch.c

@@ -1,3 +1,4 @@
+#include <float.h>
 #include <math.h>
 #include <stdlib.h>
 #include <string.h>
@@ -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) -- */
@@ -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 */

lapack-netlib/INSTALL/slamch.c

@@ -1,3 +1,4 @@
+#include <float.h>
 #include <math.h>
 #include <stdlib.h>
 #include <string.h>
@@ -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) -- */
@@ -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 */