Undo the scaling of the least-squares solution in one step when A and B were scaled to the same end of the range - #1391
Open
rmlarsen wants to merge 1 commit into
Conversation
2 tasks
Contributor
Author
|
For the reviewers, on ordering: this PR is independent of #1383 (it was found while reviewing that PR's first revision, which had copied the two-step undo) and of #1381 and #1382, and merges cleanly with each of them in either order. It shares |
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #1391 +/- ##
========================================
Coverage 69.36% 69.36%
========================================
Files 6122 6122
Lines 486711 486815 +104
Branches 23268 23268
========================================
+ Hits 337584 337702 +118
+ Misses 148689 148675 -14
Partials 438 438
Continue to review full report in Codecov by Harness.
|
rmlarsen
force-pushed
the
ls-drivers-undo-scaling
branch
from
September 8, 2026 04:11
5f4988f to
4bc258e
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 was referenced Sep 10, 2026
… B were scaled to the same end of the range xGELS, xGELST, xGETSLS, xGELSY, xGELSD and xGELSS scale A and B into [SMLNUM, BIGNUM] before solving and undo the two scalings on the solution one after the other. When both were scaled to the same end the factors share a constant that cancels, but the first step runs on its own: with both above BIGNUM it multiplies by BIGNUM/ANRM, as small as 2^-54, and flushes any solution entry below 2^-1020 that the second step, BNRM/BIGNUM, would have restored. For A = 2^1023 I and b = (2^1023, 2^-27), whose solution is (1, 2^-1050), all 24 drivers return (1, 0) with INFO = 0. The mirror case below SMLNUM applies the large factor first and could overflow, though a full-rank problem cannot reach it. When IASCL = IBSCL /= 0, apply the quotient BNRM/ANRM in one xLASCL call, which never overshoots its target; the other combinations apply two factors in the same direction and are unchanged. In xGELSY, xGELSD and xGELSS the rescaling of R or S by the factor of A alone moves into its own IF and is unchanged. Where both factors apply the solution is rounded once instead of twice and can differ from the previous result in the last bit. Found while reviewing the first revision of the xGGLSE/xGGGLM scaling change, which copied this block. The test suite never ran the block: xQRT13 scales its "scaled up" matrix to 1/(SFMIN/EPS) = 2^969, one binade inside the drivers' threshold of 2^970, so IASCL and IBSCL were zero for every matrix it generates. xQRT13 now scales up to a fixed 2^1016 instead, and xDRVLS makes the last column of the exact solution 256*SFMIN in the scaled-up types, which is small enough for the first of the two undo steps to flush it. xQRT16 normalizes its residual per right-hand side, so the lost column is visible: the parent fails between 2395 and 5041 ratios per precision, this branch none. Over a sweep of 2028 (precision, driver, exponent of A, exponent of b) cases every case whose code path is unchanged is bit-identical to the parent, and the merged cases agree with the unscaled twin as closely as before. The full LAPACK test suite passes with the same totals as the parent. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
rmlarsen
force-pushed
the
ls-drivers-undo-scaling
branch
from
September 15, 2026 20:28
4bc258e to
7668d18
Compare
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
The six least-squares drivers
xGELS,xGELST,xGETSLS,xGELSY,xGELSDandxGELSSscaleAandBinto[SMLNUM, BIGNUM]before solving and undo the two scalings on the solution one after the other. WhenAandBwere both scaled to the same end of the range the two factors share a constant that cancels, but the first step is applied on its own and can flush a solution entry to zero (both aboveBIGNUM) or overflow it (both belowSMLNUM) before the second step would have brought it back. ForA = 2^1023 I,b = (2^1023, 2^-27), whose solution isx = (1, 2^-1050), all 24 drivers returnx = (1, 0)withINFO = 0. This PR applies the quotient of the two factors in one step in that case; nothing changes when at most one ofA,Bis scaled or they are scaled to opposite ends.Description
With
IASCL = IBSCL = 2the drivers computex = x' * (BIGNUM/ANRM) * (BNRM/BIGNUM). The first factor is as small as 2^-54, so any entry ofx'below 2^-1020 (2^-125 in single precision) is rounded to zero, and the second factor cannot recover it. The mirror caseIASCL = IBSCL = 1multiplies bySMLNUM/ANRM, up to 2^104, first; an intermediate overflow there needsx'above 2^920 withbbelowSMLNUM, which a full-rank problem cannot produce, so it is merged for symmetry. The two mixed cases apply two factors in the same direction and are left as they are. The change was found while reviewing #1383, whose first revision copied this block.DGELSonA = 2^1023 I,b = (2^1023, 2^-27)x = (1, 0),INFO = 0x = (1, 2^-1050),INFO = 0(exact)Fix. In each driver, when
IASCL .EQ. IBSCL .AND. IASCL .NE. 0, callxLASCL( 'G', 0, 0, ANRM, BNRM, ... )once on the solution instead of the two calls; otherwise the existing per-flag calls run unchanged.xLASCLnever overshoots its target, so the single call cannot flush or overflow an entry whose final value is representable. InxGELSY,xGELSDandxGELSSthe rescaling ofRorS, which involves the factor ofAalone, is moved into its ownIFafter the solution block; it is unchanged. Where both factors apply, the solution is now rounded once instead of twice, so it can differ from master in the last bit. 24 files:{s,d,c,z}{gels,gelst,getsls,gelsy,gelsd,gelss}.f.Minimal reproducer (
repro/undo_cases.f90runs the case, and its mirror belowSMLNUM, through all 24 drivers)Validation
Exponent sweep, six drivers, four precisions, 2028 cases
repro/ls_sweep.f90runs each driver on a well-conditioned 6-by-4 problem withAscaled by 2^ka andbby 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, and compares the solution with the unscaled twin after rescaling;x(1)is printed in hex for a bit-for-bit diff against master.run_sweep.shreproduces the table.S/Cdrivers (12)D/Zdrivers (12)Every case whose code path is unchanged is bit-identical to master. In the merged cases the solution is rounded once instead of twice; the cases that still agree with master to the bit are those whose two norms have a power-of-two quotient. The twin problem has no solution entry near the underflow threshold, so the flush itself is covered by
undo_cases.f90, not by the sweep.Test suite. The full LAPACK test suite passes on this branch: 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 commit
f96546fc9built and run the same way (GCC 13.3,CMAKE_BUILD_TYPE=Release,BUILD_INDEX64_EXT_API=ON).Regression test.
xQRT13scales its "scaled up" matrix to1/(SFMIN/EPS)= 2^969, andxDRVLSnow scales it to a fixed 2^1016 instead, above the drivers' threshold, so that the scaling block runs at all. In the scaled-up typesxDRVLSalso multiplies the last column of the exact solution by256*SFMINbefore forming the right-hand side, which puts it below what the first of the two undo steps can represent.xQRT16normalizes its residual per right-hand side, so the flushed column is not hidden by the other columns: on the parent commit thexGELS,xGELSTandxGETSLSblocks fail between 2395 and 5041 ratios per precision, and with the fix all four precisions pass. Raising the scale also givesIASCLandIBSCLtheir first coverage in the six drivers.Why the test suite never saw it. The least-squares tests in
xDRVLStake their "scaled up" and "scaled down" matrices fromxQRT13, which scales the largest entry to1/(SFMIN/DLAMCH('Epsilon'))= 2^969 (down: 2^-969), while the drivers scale when the largest entry is outside[DLAMCH('S')/DLAMCH('P'), 1/that]= [2^-970, 2^970]. The test matrices sit one binade inside the window on both sides, soIASCLandIBSCLare 0 for every matrix the suite generates and the undo block at the end of the drivers has never run under the suite. Moving the test scale one binade outward and giving one right-hand side a solution below 2^-1020 makesxQRT16's per-column residual catch the flush (ratios of 10^3 to 10^5 on master against a threshold of 30); that test-suite change is prepared separately so that this fix stays a 24-file source change.Performance. No change on the common path: the merged call replaces two
xLASCLcalls by one and runs only when both matrices were scaled to the same end of the range.Checklist