Scale down in xLARFG and xLARFGP when |BETA| exceeds half the overflow threshold - #1381
Open
rmlarsen wants to merge 1 commit into
Open
Scale down in xLARFG and xLARFGP when |BETA| exceeds half the overflow threshold#1381rmlarsen wants to merge 1 commit into
rmlarsen wants to merge 1 commit into
Conversation
❌ 12 Tests Failed:
View the top 3 failed test(s) by shortest run time
View the full list of 77 ❄️ flaky test(s)
To view more test analytics, go to the Test Analytics Dashboard |
2 tasks
…w threshold The Householder generators rescale their input only at the small end, when |BETA| < SAFMIN. At the large end they form ALPHA-BETA (xLARFG) and ALPHA+BETA (xLARFGP), sums of two like-signed terms each bounded by |BETA|, which overflow whenever |BETA| > OVFL/2 although BETA itself and the reflector are representable. xLARFG then returns TAU = Inf and v = 0; xLARFGP gets ALPHA = Inf, TAU = 0, takes its flush branch and returns H = I with the tail of the column not annihilated. Every QR-type factorization built on them (xGEQRF, xGEQRFP, xGELQF, xGEQLF, xGERQF, xGEQRT, xGEQR, xGEQP3, xGEQP3RK, xTZRZF) returns a wrong R or Inf/NaN with INFO = 0 for a column whose leading entry exceeds about 9e307 (1.7e38 in single precision). The least-squares drivers pre-scale A and are not affected. When |BETA| > OVFL/2, scale X and ALPHA by SAFMIN, recompute XNORM and BETA, and multiply BETA back by 1/SAFMIN on exit, mirroring the existing small-end branch. The factors are powers of two, so TAU and v are the same reflector. Since |ALPHA| <= |BETA|, the sums cannot overflow for |BETA| <= OVFL/2, and the test is placed there so that no input the old code handled takes the new branch. The overflow threshold is fetched from xLAMCH( 'O' ) only when |BETA|*SAFMIN > 1, which keeps the extra cost on the common path to one multiply and one compare. The QR, LQ, QL and RQ test paths get a matrix type for it: type 9 is the random matrix of type 4 with the entry the first reflector works on raised to three quarters of the overflow threshold, A( 1, 1 ) for QR and LQ and A( M, N ) for QL and RQ. The generator cannot produce such a matrix, because it scales what it generates by the requested norm. The existing test ratios detect the failure, but only once the comparison with the threshold is guarded with xISNAN: a NaN is not .GE. THRESH, so without the guard every ratio of the new type passes on the parent commit. With it, the parent fails 5670 ratios per precision and this branch none. Over a sweep of 856704 (precision, generator, n, INCX, exponent of ALPHA, exponent of X, sign, pattern) cases the outputs are bit-identical to the parent commit whenever the true |BETA| <= OVFL/2; the parent returns Inf, NaN or a reflector with residual above tolerance in 12622 cases with a representable BETA, this branch in none, and each reflector was verified against its defining relation in quadruple precision. The full LAPACK test suite passes: 5506497 LAPACK and 315872 BLAS tests, 0 numerical errors, 0 other errors; the 64596 tests above the parent are the four new matrix types. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
rmlarsen
force-pushed
the
larfg-top-end-scaling
branch
from
September 8, 2026 03:55
d957ab3 to
9b3434c
Compare
Contributor
Author
|
Verified on an Apple M4 (macOS, Homebrew gfortran 16.2, Release build with the CI flags). With this branch merged onto current master, the full test suite passes, the new tests fail without the fix, and the reproducer behaves as described above. |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Disclaimer: This PR was prepared using Claude Code.
Summary
xLARFGandxLARFGPrescale their input only at the small end of the exponent range. At the large end they formALPHA-BETA(xLARFG) andALPHA+BETA(xLARFGP), sums of two like-signed terms each bounded by|BETA|, so the sum overflows whenever|BETA| > OVFL/2even thoughBETAand the reflector are representable.xLARFGthen returnsTAU = Infandv = 0;xLARFGPgetsALPHA = Inf, computesTAU = 0, takes its flush branch and returnsH = Iwith the tail of the column not annihilated. Every QR-type factorization built on the generators returns a wrongR, orInf/NaN, withINFO = 0for any column whose leading entry exceeds about 9e307 (1.7e38 in single precision). This PR adds the missing large-end branch to the eight generators. Nothing changes for an input the old code handled.Description
For a well-conditioned 6-by-4 matrix scaled by 2^1021 (largest entry 1.03e308), compared with the factorization of the unscaled matrix (Householder QR is exactly invariant under a power-of-two scaling):
DGEQRF,DGELQF,DTZRZF,ZGEQRFInfin the factor,INFO = 0DGEQRT,DGEQR,DGEQP3,DGEQP3RKNaN,INFO = 0DGEQRFP,DGEQLFINFO = 0At 2^1019 (largest entry 2.6e307) every routine above is bit-identical between master and this branch.
DGELS,DGELST,DGETSLS,DGELSY,DGELSDandDGELSSpre-scaleAinto[SMLNUM, BIGNUM]and are not affected;DGGLSEandDGGGLMdo not pre-scale and are the subject of a separate PR.LAWN 203 (section 2) analyzes the like-signed rearrangement and the small-end rescaling of these generators but does not consider the large end; #938 fixed the small-end
1/ALPHAoverflow ofxLARFGPonly.Fix. When
|BETA| > OVFL/2, scaleXandALPHAbySAFMIN(2^-969), recomputeXNORMandBETA, and multiplyBETAback by1/SAFMINon exit, mirroring the existing small-end branch (KNT = -1marks the case). Since|ALPHA| <= |BETA|, the sum cannot overflow when|BETA| <= OVFL/2, so the test is placed exactly there and no input that did not overflow before takes the new branch. The threshold is writtenHALF*HUGE( ZERO );HUGE( ZERO )is whatxLAMCH( 'O' )returns, and as a compile-time constant it costs the common path one compare and no arithmetic that could raise a floating-point exception. (The first revision of this PR guarded the branch with|BETA|*SAFMIN > 1to keepxLAMCH( 'O' )off the common path; that product underflows for every|BETA| < 1, so ordinary inputs such asALPHA = X(1) = 1e-100set the IEEE underflow flag and aborted under-ffpe-trap=underflow.repro/precheck_underflow.f90checks the flag after each of the eight generators.) Eight files:{s,d,c,z}larfg.f,{s,d,c,z}larfgp.f.Minimal reproducer
The same input gives
tau = InfinityfromSLARFG(with0.6*huge(1.0)) andZLARFG. ThroughDGEQRF, on the 6-by-4 matrix above:Validation
Exponent sweep of the eight generators, 856704 cases, verified in quadruple precision
repro/larfg_sweep.f90runsxLARFGandxLARFGPin all four precisions for n in {1, 2, 3, 4, 8, 33},INCXin {1, 2}, two signs ofALPHA, three tail patterns, and every pair of exponents from a grid over the whole range (every 64 binades plus the binades next to the underflow, subnormal,SAFMIN,OVFL/2andOVFLboundaries). Every output is printed in hex, and every reflector is checked against its defining relationxLARFGPreturns a negativeBETA. Cases whose trueBETAis subnormal or aboveOVFLare not judged for accuracy (the routines document the loss of accuracy for a subnormal result) but still must not produce a non-finite output whenBETAis representable.|BETA| > OVFL/2|BETA| <= OVFL/2SLARFGSLARFGPDLARFGDLARFGPCLARFGCLARFGPZLARFGZLARFGPThis branch never returns a non-finite output when the true
BETAis representable.xLARFGPfails less often on master because its flush branch happens to give an acceptableH = Iwhen the tail is negligible relative toALPHA.Among the changed cases where master was also accurate, the two libraries compute the same reflector to rounding; the maximum residual over those cases is smaller on this branch in every precision (double: 1.1e-15 on master vs 6.9e-16 here; complex double: 1.8e-15 vs 9.4e-16; single: 7.2e-7 vs 3.5e-7; complex single: 1.0e-6 vs 4.1e-7). In single precision a handful of cases (8 to 20 per generator) move by one or two ulps in the other direction, all within tolerance. The changes come from
xNRM2rounding the scaled tail differently, and they occur only for|BETA| > OVFL/2, where master's result is either these same bits or an overflow.run_sweep.shreproduces the table;classify_sweep.pyproduces it from the two output files.Regression test. The$\lvert\alpha\rvert \approx \|a\|_2/\sqrt{M}$ , so the sum that overflows needs a single entry close to the norm of its column. The one-norm of the new matrix stays finite, so the orthogonality ratio and the
?QR,?LQ,?QLand?RQpaths get a matrix type 9: the random matrix of type 4 with the entry the first reflector works on raised to three quarters of the overflow threshold,A( 1, 1 )for QR and LQ andA( M, N )for QL and RQ. The generator cannot produce such a matrix, sincexLATMSscales what it generates to the requested norm, and a whole matrix scaled to the top of the range does not trigger the defect either: for a random column,xORGQRandxORMQRratios remain meaningful.The existing ratios detect the failure only once their comparison with the threshold is guarded with
xISNAN, because a NaN is not.GE. THRESHand every ratio of the new type would otherwise pass on the parent commit. With the guard, the parent fails 5670 ratios per precision and this branch none. The guard is on the print loop of the four checkers, so it also covers the existing types.Test suite. The full LAPACK test suite passes on this branch: 215 of 215 CTest entries, 5506497 LAPACK tests and 315872 BLAS tests with 0 numerical errors and 0 other errors, built and run the same way as the parent commit
f96546fc9(GCC 13.3,CMAKE_BUILD_TYPE=Release,BUILD_INDEX64_EXT_API=ON). The 64596 tests above the parent are the four new matrix types; no existing type produced a NaN ratio that the new guard turned into a failure.Performance. Timed on a 13th Gen Intel(R) Core(TM) i7-13700HX under WSL2 with the reference BLAS, GCC 13.3,
-O2. To separate the change from code-placement effects (which move untouched routines by up to 28% between two separately linked static libraries on this machine), the parent library is a shared object shared by both sides, and each benchmark binary carries its own copy of only the changed routines, parent or branch, which interposes over the library's; everything else is byte-identical. Four rounds in alternating order, one core, one process per run, on an idle machine; medians of the per-round medians, withDPOTRFas an untouched control. The benchmark driver and raw output are available on request.DLARFG(ns/call)DGEQRF(us/call)DPOTRFcontrol (us/call)No measurable difference: the common path gains one compare per reflector, below the resolution of the measurement (the round-to-round spread of
DLARFGis 2 ns). The timing was taken on the first revision, whose common path also carried a multiply; the current one does strictly less work there.Checklist