From ffe8b242c765c1b7180d949ca9dd5e1cc164392c Mon Sep 17 00:00:00 2001 From: David Heffernan Date: Tue, 11 Aug 2026 11:26:11 +0100 Subject: [PATCH] C_LAPACK: take the float format from 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 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. --- lapack-netlib/INSTALL/dlamch.c | 78 +++++++++++++--------------------- lapack-netlib/INSTALL/slamch.c | 78 +++++++++++++--------------------- 2 files changed, 58 insertions(+), 98 deletions(-) diff --git a/lapack-netlib/INSTALL/dlamch.c b/lapack-netlib/INSTALL/dlamch.c index ce6b76a32a..5eaa1e595d 100644 --- a/lapack-netlib/INSTALL/dlamch.c +++ b/lapack-netlib/INSTALL/dlamch.c @@ -1,3 +1,4 @@ +#include #include #include #include @@ -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,24 +370,23 @@ 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 */ @@ -410,32 +394,28 @@ doublereal dlamch_(char *cmach) 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 */ diff --git a/lapack-netlib/INSTALL/slamch.c b/lapack-netlib/INSTALL/slamch.c index 2def9337a7..3ffc715d17 100644 --- a/lapack-netlib/INSTALL/slamch.c +++ b/lapack-netlib/INSTALL/slamch.c @@ -1,3 +1,4 @@ +#include #include #include #include @@ -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,24 +369,23 @@ 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 */ @@ -409,32 +393,28 @@ real slamch_(char *cmach) 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 */