Scale the inputs of xGGLSE and xGGGLM into the safe range, as the other least-squares drivers do - #1383
Open
rmlarsen wants to merge 1 commit into
Open
Scale the inputs of xGGLSE and xGGGLM into the safe range, as the other least-squares drivers do#1383rmlarsen wants to merge 1 commit into
rmlarsen wants to merge 1 commit into
Conversation
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #1383 +/- ##
==========================================
+ Coverage 69.01% 69.38% +0.36%
==========================================
Files 6122 6122
Lines 486123 486809 +686
Branches 23286 23268 -18
==========================================
+ Hits 335514 337766 +2252
+ Misses 150420 148605 -1815
- Partials 189 438 +249
... and 148 files with indirect coverage changes Continue to review full report in Codecov by Harness.
|
2 tasks
rmlarsen
force-pushed
the
gglse-ggglm-scaling
branch
from
September 8, 2026 04:38
996fd4a to
89c4bea
Compare
…y, by exact powers of two Review of the previous revision found three regressions against the unscaled drivers. A and B shared one factor from their combined max-norm, so A = 2^1023, B = 2^-1023, d = 2^-1023 (x = 1) had B rounded to zero and returned INFO = 1. The two undo scalings of x cancelled through an intermediate that could flush, so x = (1, 2^-1050) came back as (1, 0). All of c was rescaled although only c(N-P+1:M) holds the residual, and the intermediate left in c(1:N-P) could overflow. Scale each operand by a power of two chosen from the EXPONENT of its max-norm instead. A and B get independent factors 2^KA and 2^KB. In xGGLSE, c carries 2^(KA+KT) and d carries 2^(KB+KT), which keeps the constraint and the residual consistent and scales x by 2^KT, with KT chosen to bring the larger right-hand side into range; on exit x is rescaled by 2^-KT and the residual in c(N-P+1:M) by 2^-(KA+KT). In xGGGLM, A, B and d are scaled independently and x and y are rescaled by 2^(KA-KD) and 2^(KB-KD). The factors are exact, so a problem the unscaled computation handles returns the same solution to the last bit, and each undo is one exact operation. xLASCL is called with both endpoints at or above one so that forming the factor raises no underflow; vectors are scaled with xSCAL. A zero, infinite or NaN norm takes no part. The three review cases and two neighbouring corners agree with the parent bit for bit and raise the same IEEE flags. Over the exponent sweep of 2860 cases the branch fails in none (the parent in 160) and, in the 1480 scaled cases with normal-range inputs, returns the solution of the unscaled twin bit for bit in 1388; the rest have entries beyond the thresholds at which xNRM2 switches accumulators and agree with the twin to rounding. xGLMTS and xLSETS solve their problem a second time with every operand scaled by one power of two, so that the largest entry sits just below the overflow threshold. The problem is exactly invariant under that scaling, so the residual of the second solution against the original data has to match the first; a NaN counts as a failure, which the comparison with the threshold would otherwise pass over. On the parent commit the GLM twin fails one of the 48 ratios in every precision. The LSE twin passes on both, because the right-hand side is the largest operand there and scaling it to the top of the range leaves the matrix below the exponent at which the reflector generator overflows; it covers the new scaling path rather than the defect. The full LAPACK test suite passes with the same totals as the parent. xGLMTS and xLSETS declare the new xSCAL calls EXTERNAL: the extended-API build renames only the routines a file declares, so without the declaration the xeigtst*_64 executables failed to link against the 64-bit BLAS. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
rmlarsen
force-pushed
the
gglse-ggglm-scaling
branch
from
September 10, 2026 07:18
89c4bea to
c7a7e38
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. Two differences from x86:
|
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
xGELS,xGELST,xGETSLS,xGELSY,xGELSDandxGELSSscaleAand the right-hand side into[SMLNUM, BIGNUM]before factoring and undo the scaling on the solution.xGGLSEandxGGGLMare the only least-squares drivers that do not, so a problem whose entries approach the overflow threshold is factored and solved as given, and a well-conditioned problem scaled by 2^1021 comes back asNaN, or as a finite wrongx, withINFO = 0. This PR adds the scaling block ofxGELSto both drivers. Nothing changes for a problem whose entries already lie in[SMLNUM, BIGNUM].Description
For
min ||c - A x||subject toB x = dwith a well-conditioned 6-by-4A, 2-by-4Band consistentc,d, all scaled by 2^1021 (largest entry 1e308); the solution is invariant under that scaling:INFOxNaN NaN NaN NaN0.1 0.2 0.3 0.4(exact solution)DGGLSEcomputes the GRQ factorization of(B, A)in place, and its Householder generators overflow first (a separate PR); with that fixed, the reflector applications inxORMQR, the triangular solves and the matrix-vector products still overflow, so the drivers need the same protection the other least-squares drivers have. The documentation of both routines describes their rank checks as rudimentary but says nothing about range.Fix. After the quick return, each driver scales
A,Band the right-hand sides by exact powers of two so that their largest entries lie in[SMLNUM, BIGNUM): the shift of each operand is computed from theEXPONENTof its max-norm, applied withxLASCL(matrices) orxSCAL(vectors), and undone on exit by the opposite shift.AandBget independent factors2^KAand2^KB. InxGGLSE,cis scaled by2^(KA+KT)anddby2^(KB+KT), which keeps the constraintB x = dand the residualc - A xconsistent and scalesxby2^KT;KTbrings the larger of the two right-hand sides, each measured against its matrix, into range. On exitxis rescaled by2^-KTand the residual left in elementsN-P+1toMofcby2^-(KA+KT); the intermediate left inc(1:N-P)is not touched. InxGGGLMthe three operands are independent:dis scaled by2^KD, andxandyare rescaled by2^(KA-KD)and2^(KB-KD). Scaling by a power of two is exact, so a problem the unscaled computation can handle returns the same solution to the last bit (sweep below), and each undo is one exact operation. A norm that is zero, infinite or NaN takes no part.xLASCLis called with both endpoints at or above one, so forming the factor raises no spurious underflow, and a problem whose entries lie in[SMLNUM, BIGNUM)calls neitherxLASCLnorxSCAL. As inxGELS, the factorizations returned inAandBare those of the scaled matrices when scaling took place.xLANGE( 'M' )does not touch its work argument, so no workspace is needed and the workspace query is unchanged. Eight files:{s,d,c,z}gglse.f,{s,d,c,z}ggglm.f.Known limit, the same as in
xGELS: when a right-hand side has to be scaled down (largest entry aboveBIGNUM), its entries more than about 2^1990 below that entry lose precision and entries more than 2^2044 below it are lost; inxGGLSEthe same holds between2^KA cand2^KB d. InxGGGLM, whered = A x + B ycarries the scales of bothAandB, adwhoseB ypart is that far below itsA xpart loses it:A = 2^1023,B = 2^-1023,d = (2^1023, 2^-1023)returnsy = 0instead of1. Such adis representable only when the two parts occupy separate entries, as here; master returnsy = 1because the factorization of thatAis trivial.Revision 2. A review of the first revision found three regressions against master, all reproduced by
repro/review_cases.f90:AandBshared one factor from their combined max-norm, so forA = 2^1023,B = 2^-1023,c = 1,d = 2^-1023(exact solutionx = 1) the scaling roundedBto zero and the driver returnedINFO = 1.xcancelled through an intermediate that could flush:x = (1, 2^-1050)came back as(1, 0).cwas rescaled although onlyc(N-P+1:M)holds the residual; the intermediate inc(1:N-P)could overflow, and trap, on a problem master solves.The revision fixes the three by construction (independent factors for
AandB; exact factors undone in one operation; only the residual part ofcrescaled). The fivexGGLSEcases ofreview_cases.f90, the three findings and two neighbouring corners, agree with master bit for bit and raise the same IEEE flags. The same two-step undo exists in the sixxGELS-family drivers (xGELS,xGELST,xGETSLS,xGELSY,xGELSD,xGELSS), whose first step flushes for the same input; #1391 fixes it there. The two PRs are independent and can land in either order.Minimal reproducer
Validation
Exponent sweep of both drivers, all precisions, 2860 cases
repro/gg_sweep.f90runsxGGLSEon the shapes(M,N,P)= (6,4,2), (4,4,4), (5,4,0), (2,4,4), (3,3,1) andxGGGLMon(N,M,P)= (6,4,6), (6,4,2), (5,5,0), (6,0,6), (3,1,3) in all four precisions, with the matrices scaled by 2^ka and the right-hand sides by 2^kb over a grid of exponents from the subnormal range to the overflow threshold (double: -1070 to 1021; single: -148 to 125), skipping pairs whose exact solution is not representable. The solution of each case is compared with the unscaled twin after rescaling, and printed in hex for a bit-for-bit comparison with master. Cases whose inputs carry fewer than about 12 significant bits (deep subnormal) are excluded from the accuracy judgment.SGGLSE/CGGLSEDGGLSE/ZGGLSESGGGLM/CGGGLMDGGGLM/ZGGGLMA case fails when
INFO /= 0, the solution containsNaNorInf, or it deviates from the unscaled twin by more than 5e-13 (2e-5 single; 2e-11 and 2e-3 when the inputs are subnormal and carry fewer bits). Of the 1480 cases with inputs in the normal range in which scaling is applied, this branch returns the solution of the unscaled twin bit for bit in 1388; the other 92 (single precision, andZGGGLMwith an exponent beyond 486) have entries beyond the thresholds at whichxNRM2switches accumulators, and agree with the twin to 1.6e-6 (single) and 7e-22 (double) relative to the largest solution entry. Master returns the twin bit for bit in 1124 of the 1480 and fails in 160. All of master's failures are at the largest matrix exponent (2^1021, or 2^125 in single precision); thexGGGLMfailures there include finite, wrong solutions returned withINFO = 0.run_sweep.shreproduces the table.Regression test.
xGLMTSandxLSETSsolve their problem a second time with every operand scaled by one power of two, so that the largest entry sits just below the overflow threshold. The problem is exactly invariant under that scaling, so the residual of the second solution against the original data has to match the first, and a NaN result counts as a failure, which the plain comparison with the threshold would pass over.On the parent commit the GLM twin fails one of the 48 ratios in every precision, and with the fix all of them pass. The LSE twin passes on both: there the right-hand side is the largest operand, so scaling it to the top of the range leaves the matrix below the exponent at which the reflector generator overflows. It covers the new scaling path rather than the defect, which is what the existing types could not reach at all, since
xLATB9fixes the norms ofAandBat 10 and 1000 for these two paths.Test suite. The full LAPACK test suite passes on this branch, including the
xLSEandxGLMtests of the eigenvalue suite: 215 of 215 CTest entries, 5441901 LAPACK tests and 315872 BLAS tests with 0 numerical errors and 0 other errors, the same totals as the parent commitf96546fc9built and run the same way (GCC 13.3,CMAKE_BUILD_TYPE=Release,BUILD_INDEX64_EXT_API=ON).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.DGGLSE(us/call)DPOTRFcontrol (us/call)The cost is the four max-norm scans over
A,B,candd(O(mn + pn)against theO(mn^2)factorization), about 1.4 us on the 32-by-16 problem and within the measurement noise from n = 256 on; it is the same workxGELSdoes. Nothing else runs on a problem whose entries lie in[SMLNUM, BIGNUM). (Revision 2 numbers; revision 1 measured 1.17 / 1.06 / 1.01 / 1.01.)Checklist
xGELS.)