C_LAPACK: take the float format from <float.h> in dlamch/slamch

[davidheff]

lapack-netlib/INSTALL/dlamch.c and slamch.c are f2c translations of the deprecated
dlamchf77.f/slamchf77.f, which determine the floating point format at run time by
probing in dlamc1/dlamc2. The current dlamch.f uses the Fortran 90 inquiry
intrinsics instead, but f2c cannot translate those, so the C LAPACK selected by
NOFORTRAN=1 was 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, where they stay in 80 bit registers and the probe measures the
register format instead of the storage format.

This affects stock 32 bit builds, not just unusual ones

Makefile.x86 adds -msse2 but never -mfpmath=sse, and on 32 bit x86 -msse2 enables
the instructions without changing the FP math default away from 387. With gcc 16.1.0:

-msse2                 ->  fldl / fmull / faddl     (x87)
-msse2 -mfpmath=sse    ->  movsd / mulsd / addsd    (SSE)

So a stock 32 bit NOFORTRAN=1 build compiles dlamch.c with x87 double arithmetic,
which is exactly the condition that breaks the probe.

Reproduction

Linking INSTALL/dlamch.c and INSTALL/slamch.c from current develop against a test
program that compares every cmach against <float.h>, built with the stock 32 bit
flags -O2 -msse2, 16 of the 20 constants are wrong:

	before	expected
dlamch('E')	5.4210108624275222e-20	1.1102230246251565e-16
dlamch('S')	0	2.2250738585072014e-308
dlamch('N')	64	53
dlamch('M')	-16381	-1021
dlamch('L')	16385	1024
dlamch('O')	inf	1.7976931348623157e+308
slamch('S')	0	1.1754943508222875e-38
slamch('N')	64	24
slamch('O')	inf	3.4028234663852886e+38

Note that slamch reports the 80 bit format too, and that the mantissa width follows the
caller's x87 precision control bits rather than being a property of the type at all.

It is not always a silent wrong answer. With the x87 control word 0x1332 (extended
precision, overflow unmasked, which is what Delphi sets) the unpatched code terminates
with STATUS_FLOAT_OVERFLOW (0xC0000091) inside dlamch itself.

Where the values are merely wrong rather than fatal, everything that scales by them is
wrong. The first symptom we saw was dsbevx, which 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.

With this patch all 20 constants are correct, at the default control word and at 0x1332
and 0x037F.

x86_64 is not affected by this part, since -mfpmath=sse is the default there and the
probe measures the right format. I have verified that separately.

Two further problems fixed at the same time, on all architectures

These are not x87 specific and apply to every NOFORTRAN=1 build, x86_64 included:

• rmach was left uninitialised when cmach matched nothing, where dlamch.f returns
  zero.
• The cached static results made both routines unsafe to call concurrently on first
  use.

What the patch does

Replaces the probe with the <float.h> constants, mirroring the values the current
dlamch.f and slamch.f return via the F90 intrinsics.

dlamc1-dlamc5 and slamc1-slamc5 are deliberately left in place. They become
unreachable from dlamch/slamch, but dlamc3/slamc3 have callers of their own in
dlaed3, dlaed9, dlals0, dlasd3, dlasd8 and their complex equivalents, where they
serve as optimiser barriers. Happy to remove the genuinely dead dlamc1/2/4/5 in a
follow-up if you would prefer.

Regenerating these files from the modern dlamch.f is not an option, since f2c cannot
translate the F90 inquiry intrinsics — which is how the C path came to be stuck on the
deprecated version in the first place.

These two .c files are OpenBLAS's own artifacts, added in #3539 and hand maintained
since (#3605, 4041b7fb4); Reference-LAPACK ships only the .f, and its dlamch.f is
already correct. So this does not need to go to Reference-LAPACK, and will not be
clobbered by a future re-sync.

gcc -Wall reports no new warnings; the patch removes one of the five that current
develop produces for this file.

This has been running in production in our fork since 0.3.34.

🤖 Generated with Claude Code

[martin-frbg]

Ah right, thanks - the big f2c'ing was mostly a hack to keep the LAPACK parts available on Android after the NDK switched away from GCC (and also for the people who insist on compiling with plain MSVC under Windows...). I had meant to revise it at some point, if only to remove the remaining compiler warnings.

In the not too distant future, it might make sense to replace it with ilayn 's semicolon-lapack that is a complete C11 rewrite of the latest Reference-LAPACK rather than a half-baked machine translation of something resembling 3.9.0

[davidheff]

Thanks for merging, and yes replacing the f2c version with semicolon-lapack does sound like a better option.