Do not scale an infinite matrix to zero in xSTEMR; report it from xLARRE - #1404
Open
rmlarsen wants to merge 1 commit into
Open
Do not scale an infinite matrix to zero in xSTEMR; report it from xLARRE#1404rmlarsen wants to merge 1 commit into
rmlarsen wants to merge 1 commit into
Conversation
xSTEMR scales the matrix into the range xLARRE can represent:
ELSE IF( TNRM.GT.RMAX ) THEN
SCALE = RMAX / TNRM
For a matrix with an infinite entry TNRM is infinite, SCALE underflows
to zero, the scaling turns the infinite entry into a NaN and every
other entry into zero, and the solver ran on that matrix: RANGE = 'A'
and 'I' returned INFO = 0 with every eigenvalue a NaN, RANGE = 'V'
scaled the interval to (0, 0] and returned no eigenvalues, and with
the infinite entry off the diagonal the representation tree search
never returned (the NaN hang of Reference-LAPACK#1389).
Leave a matrix whose scale factor would be zero unscaled, and let
xLARRE, which computes the base representations, reject a matrix with
an infinite entry as one for which no representation exists, INFO = 2,
the code it already uses when the search for a representation gives up.
xSTEMR returns INFO = 12 and the drivers that call it, xSTEVR, xSYEVR
and xHEEVR, fall back to xSTEBZ and xSTEIN as for any other xSTEMR
failure. A NaN matrix is not touched by this change: its norm
compares false against every threshold and it takes the same path as
before.
xCHKST gets the case as a regression test: after the size and type
loops it calls xSTEMR with RANGE = 'A' on an N = NMAX tridiagonal
matrix whose first diagonal entry is +Inf, once with JOBZ = 'V' and
once with JOBZ = 'N', and reports INFO = 0 as a failure. On the
parent commit all four calls return INFO = 0 with NaN eigenvalues;
here they return INFO = 12. Over 420 placements of +Inf and -Inf
(n = 3 to 64, RANGE = 'A', 'V', 'I', real and complex) the parent
returned INFO = 0 with NaN eigenvalues in 52 cases, hung in 24 and
returned INFO = 11 in the rest; this branch returns INFO = 12 in all
of them. The 1 by 1 and 2 by 2 matrices, which xSTEMR solves in closed
form before scaling, are unchanged, and every finite case is
bit-identical.
The full LAPACK test suite passes: 0 numerical errors, 0 other errors,
40 tests more than the parent from the new calls.
The regression test fails on the parent and passes with the fix, and the
reproducer prints the same before and after output, with gfortran 13
(x86-64 Release and Debug with -fcheck=all, and under QEMU on aarch64,
ppc64le, s390x and riscv64), flang-19 and Intel ifx 2025.3.
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #1404 +/- ##
========================================
Coverage 69.36% 69.36%
========================================
Files 6122 6122
Lines 486337 486447 +110
Branches 23268 23268
========================================
+ Hits 337330 337440 +110
Misses 148569 148569
Partials 438 438
Continue to review full report in Codecov by Harness.
|
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
xSTEMRscales the matrix into the range its representations can handle,SCALE = RMAX / TNRMwhen the max-normTNRMexceedsRMAX. For a matrix with an infinite entryTNRMis infinite,SCALEunderflows to zero, the scaling turns the infinite entry into a NaN and every other entry into zero, and the solver ran on that matrix:RANGE = 'A'and'I'returnedINFO = 0with every eigenvalue a NaN,RANGE = 'V'scaled the interval to(0, 0]and returned no eigenvalues, and with the infinite entry off the diagonal the representation search never returned (the NaN hang of #1389). This PR leaves such a matrix unscaled and letsxLARREreject a matrix with an infinite entry as one for which no base representation exists, soxSTEMRreturnsINFO = 12and its drivers fall back toxSTEBZandxSTEIN. Six files:{s,d,c,z}stemr.f,{s,d}larre.f.Description
RMAXis about1e77in double precision, soRMAX / TNRMis at least5e-232for every finiteTNRMand is zero exactly whenTNRMis infinite. The scaled matrix hasInf * 0 = NaNat the infinite entry and zeros elsewhere;xLARREthen computes Gershgorin bounds and representations from that, andxLARRVorxLARRJrefine them. At the end the eigenvalues are divided bySCALE,0 * Inf = NaN.Fix.
xSTEMRkeepsSCALE = ONEwhen the quotient is zero, which happens for no finite matrix.xLARREtests the max-norm of the matrix it is given against the overflow threshold right after its Gershgorin bounds and returnsINFO = 2, the code it already uses when the search for a base representation gives up (documented as "No base representation could be found in MAXTRY iterations", now "..., or the matrix has an infinite entry").xSTEMRreports that asINFO = 12under its existing1Xconvention.xSTEVR,xSYEVRandxHEEVRtreat any nonzeroINFOfromxSTEMRas a reason to fall back toxSTEBZandxSTEIN, which #1384 makes pass an infinite matrix through unscaled. A NaN matrix is not touched by this change: its norm compares false against every threshold, on both sides, and it takes the same path as before, which #1389 bounds.Minimal reproducer
With
E(1) = +Infinstead, master never returns from the same call.Regression test.
xCHKSTgets the case: after the size and type loops it callsxSTEMRwithRANGE = 'A'on anN = NMAXtridiagonal matrix whose first diagonal entry is+Inf, once withJOBZ = 'V'and once withJOBZ = 'N', and reportsINFO = 0as a failure. On the parent commit all four calls returnINFO = 0with NaN eigenvalues in every precision; with the fix they returnINFO = 12. The position was chosen so that the test fails rather than hangs on the parent.Validation
ctest: 100% of 215 tests passed). The 40 extra tests over the parent commit are the new calls, 2 job options times 5 parameter sets times 4 precisions.xSTEMRwithJOBZ = 'V'andRANGE = 'A','V','I'; +Inf and -Inf inD(1),D(2),D(N),E(1),E(N/2),E(N-1);N = 1, 2, 3, 5, 8, 26, 40, 64;DSTEMRandZSTEMR; one process per case with a 3 s timeout): forN >= 3the parent returnedINFO = 0with NaN eigenvalues in 52 cases, hung in 24 (RANGE = 'I'with the Inf off the diagonal), and returnedINFO = 11, anxLARREfailure on the zero-scaled matrix, in the remaining 344; this branch returnsINFO = 12in all 420. The 1 by 1 and 2 by 2 cases, whichxSTEMRsolves in closed form before scaling, are unchanged.-fcheck=all, and under QEMU on aarch64, ppc64le, s390x and riscv64), flang-19 and Intel ifx 2025.3 (-fp-model=strict).Found while auditing the symmetric tridiagonal eigensolvers for the NaN and overflow handling of #1377-#1391; this is the
xSTEMRinstance of the scale-to-zero defect fixed for the bisection drivers in #1384.