Divide by a subnormal pivot in the banded LU factorizations instead of scaling by its reciprocal - #1379
Open
rmlarsen wants to merge 1 commit into
Open
Divide by a subnormal pivot in the banded LU factorizations instead of scaling by its reciprocal#1379rmlarsen 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 #1379 +/- ##
==========================================
+ Coverage 69.01% 69.36% +0.34%
==========================================
Files 6122 6122
Lines 486123 486409 +286
Branches 23286 23268 -18
==========================================
+ Hits 335514 337398 +1884
+ Misses 150420 148573 -1847
- Partials 189 438 +249
... and 144 files with indirect coverage changes Continue to review full report in Codecov by Harness.
|
rmlarsen
force-pushed
the
band-lu-sfmin
branch
from
September 8, 2026 04:06
d32bd02 to
a0646b5
Compare
…f scaling by its reciprocal
xGBTF2 and xGBTRF form the multipliers of column J as
CALL xSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), 1 )
When the pivot is subnormal its reciprocal is not representable, so the
whole column of L becomes Inf, or NaN where an entry is zero, the xGER
update spreads that through the band, and the factorization completes
with INFO = 0. Partial pivoting cannot avoid it: the pivot is the
largest entry of the column, and when the matrix is small every column
is. A well-conditioned banded system scaled to 2^-1030 is solved by
xGESV to 1e-14 and by xGBSV to NaN.
The dense routines xGETF2 and xGETRF2 have guarded this since at least
LAPACK 3.2 by dividing element-wise when |pivot| < SFMIN. Apply the
same test at both banded sites, in all four precisions. SFMIN comes
from xLAMCH('S'), computed once after the quick return as in xGETF2.
The GB test path gets a matrix type for it: type 9 is the type 1
matrix scaled into the subnormal range, one eighth of the safe minimum,
which xLATMS cannot generate because it scales its output to the
requested norm. At that scale the matrix carries too few bits to
reconstruct, so the type tests what the guarded division promises, that
the factor is finite, instead of a residual; the remaining ratios,
which estimate a condition number from a subnormal norm, are skipped as
for a block size other than the first. On the parent commit the type
fails 798 times per precision.
The new branch is taken only when |pivot| < SFMIN. Over a sweep of
2124 (precision, n, KL, KU, scale) cases the banded factor and IPIV are
bit-identical to the parent commit for every pivot of normal magnitude,
and for the subnormal cases the banded solution error now equals the
dense one in every case, where before 560 of 708 came back NaN or Inf.
The full LAPACK test suite passes: 5447193 LAPACK and 315872 BLAS
tests, 0 numerical errors, 0 other errors; the 5292 tests above the
parent are the new type.
The test files declare the new xSCAL calls EXTERNAL: the extended-API
build renames only the routines a file declares, so without the
declaration the xlintst*_64 executables failed to link against the
64-bit BLAS.
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
rmlarsen
force-pushed
the
band-lu-sfmin
branch
from
September 10, 2026 07:18
a0646b5 to
8ad354a
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. |
2 tasks
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 banded LU factorizations
xGBTF2andxGBTRF(all eight:S/D/C/Z× unblocked/blocked) form the multipliers of columnJasxSCAL( KM, ONE / pivot, ... ). When the pivot is subnormal,ONE / pivotoverflows toInf, the column ofLbecomesInf/NaN, and the factorization completes withINFO = 0.xGBSVthen returns aNaNsolution as success for a system that is perfectly well conditioned. The dense LU (xGETF2,xGETRF2) has guarded this since at least LAPACK 3.2 (it is present in the 2008 import of the trunk) by dividing element-wise when|pivot| < SFMIN; this PR adds the same guard to the banded routines. Nothing changes for a pivot of normal magnitude.Description
xGETF2computes the column of multipliers asand
xGETRF2the same at the base of its recursion. The banded routines have only the first branch:For a pivot below
2^-1022(double) the reciprocal is not representable, so every multiplier in the column isInf, orNaNwhere the entry is zero; thexGERtrailing update then spreads that through the rest of the band. Partial pivoting does not help: it picks the largest entry of the column, and when the whole matrix is small every column is.xGTTRFandxGETC2are unaffected because they divide. The tridiagonal driverxGTSVand the dense driverxGESVboth solve the reproducer below to full accuracy.Fix. The
xGETF2form at both sites: keep the reciprocal scaling when|pivot| >= SFMIN, divide element-wise otherwise. Eight files.SFMIN = xLAMCH('S')is computed once after the quick return, as inxGETF2; neither banded routine referencedxLAMCHbefore.Minimal reproducer
A banded, diagonally dominant matrix with every entry near
2^-1030;b = A * ones, so the exact solution is the vector of ones.Validation
Band versus dense sweep, 2124 cases per build
xGBSVagainstxGESVon the same non-symmetric, diagonally dominant banded matrix withb = A * ones: four precisions,nin {1, 2, 3, 5, 8, 13, 40, 64, 130},KLandKUeach in {0, 1, 2, 3, 5, 40} (so both the unblockedxGBTF2path and, forKL = 40 >= NB, the blockedxGBTRFpath), and three scalings: normal (s = 1), subnormal (s = 2^-1030in double,2^-140in single), and near the overflow threshold (s = huge/64). Each line recordsINFOandmax|x - 1|for both drivers plus FNV-1a hashes of the banded factor andIPIV.< 1e-3, master< 1e-3, this branch< 1e-3NaN/Inf, masterNaN/Inf, this branchOn master, 560 of the 708 subnormal cases come back
NaNorInfwithINFO = 0; the 148 that survive haveKL = 0orn = 1, where there is no multiplier to scale. On this branch the banded error equals the dense error in every case (the largest ratioerr_band / err_denseover all finite pairs is 1.00 in all four precisions), andINFOagrees with the dense driver in 708 of 708. The 54 subnormal cases where both drivers exceed1e-3are all single precision: at2^-140aREALcarries about 9 significant bits, so that is the input, not the solver.For the normal and near-max scalings all 1416 lines, hashes included, are byte-identical between master and this branch: the new branch is reached only when
|pivot| < SFMIN.Regression test. The
?GBpath gets a matrix type 9: the type 1 matrix scaled in place into the subnormal range, one eighth of the safe minimum, which is below the reciprocal of the overflow threshold in every precision.xLATMScannot generate such a matrix, since it scales its output to the requested norm, and it returns Inf or NaN entries for a subnormalANORMon the shapes that go through its Givens chase.At that scale the stored matrix carries too few bits for a reconstruction residual to mean anything: in single precision an entry near 2^-134 keeps about fifteen, and the ratio of a correct factorization lands near 300. The type therefore checks what the guarded division promises, that the factor is finite, through the max-norm of the factor rather than
xGBT01. The condition-number and solve ratios, which would form the reciprocal of a subnormal norm, are skipped for it as they already are for a block size other than the first.On the parent commit the type fails 798 times per precision, one for each shape and block size, with the factor holding an infinity; with the fix all four precisions pass.
Test suite. The full LAPACK test suite passes on this branch: 215 of 215 CTest entries, 5447193 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. That includes the?GBroutine and driver families (DGB: 30261 routine tests, 36567 driver tests) and the_64extended-API variants. Built with GCC 13.3,CMAKE_BUILD_TYPE=Release,BUILD_INDEX64_EXT_API=ON.Checklist