Fix rescaling after a partial failure in the xSTEBZ-based drivers; accept Inf input and xSTEBZ's negative block numbers - #1384
Open
rmlarsen wants to merge 1 commit into
Conversation
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #1384 +/- ##
=======================================
Coverage 69.36% 69.36%
=======================================
Files 6122 6122
Lines 486711 486715 +4
Branches 23268 23268
=======================================
+ Hits 337584 337614 +30
+ Misses 148689 148663 -26
Partials 438 438
... and 1 file with indirect coverage changes Continue to review full report in Codecov by Harness.
|
rmlarsen
force-pushed
the
symeig-driver-scaling
branch
from
September 8, 2026 04:27
c9bdaf6 to
63ed409
Compare
rmlarsen
changed the title
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. |
…drivers The 32 symmetric eigensolver drivers that reduce to xSTEBZ + xSTEIN (xSTEVX, xSYEVX, xHEEVX, xSPEVX, xHPEVX, xSBEVX, xHBEVX, xSTEVR, xSYEVR, xHEEVR and the _2STAGE variants) scale a matrix whose largest entry lies outside [RMIN, RMAX] and undo the scaling of W at the end with IMAX = M if INFO = 0 and IMAX = INFO - 1 otherwise. That block was copied from xSTEV and xSYEV, where a nonzero INFO comes from xSTEQR. Here a nonzero INFO means that INFO eigenvectors failed to converge in xSTEIN, in which case the eigenvalues are all valid and W(INFO:M) are returned still multiplied by SIGMA, or that xSTEBZ returned INFO - N, in which case IMAX = N + i - 1 and xSCAL writes past the end of W. A failure of xSTERF, xSTEQR or xSTEMR earlier in the driver resets INFO and falls through to xSTEBZ, so at the rescale label M is the right count in every case. Two related defects sit in the same drivers. A matrix containing an Inf has ANRM = Inf and SIGMA = RMAX / ANRM = 0, so for RANGE = 'V' the scaled interval collapses to VLL = VUU = 0 and xSTEBZ stops the process in XERBLA; compute SIGMA first and scale only if it is positive, which sends an infinite matrix down the same path as one containing a NaN. And xSTEIN rejects with INFO = -6, again XERBLA, or silently skips the negative block numbers with which xSTEBZ flags eigenvalues that did not converge, although every driver deliberately continues into xSTEIN when xSTEBZ returns INFO = 1; take ABS( IBLOCK( . ) ) at its four uses and document it. DSTEVX on 1e-150 tridiag(1, 0.5, 1) with ABSTOL = 0.3 |T| reports INFO = 6 and returns the sixth eigenvalue 1e4 times too large on the parent, all six correctly on this branch. xERRST gets the infinite matrix as a regression test: xSTEVX and xSYEVX with RANGE = 'V' on a tridiagonal whose first diagonal entry is an infinity, which must return without reaching XERBLA. On the parent both report the illegal fifth argument of xSTEBZ. The other two changes are not reachable from the test suite: a partial xSTEIN failure needs a matrix outside the range the drivers document, and the negative block numbers likewise. Over 1584 NaN/Inf cases (8 routines, 3 special values, 5 positions, 5 orders, 3 ranges, one process each with XERBLA overridden and canaries around every array) the parent aborts 280 times and writes past the end of W 200 times; this branch never does either. The full LAPACK test suite passes: 5441901 LAPACK tests, 0 numerical errors, 0 other errors, the same totals as the parent. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
rmlarsen
force-pushed
the
symeig-driver-scaling
branch
from
September 15, 2026 20:28
63ed409 to
4f91c6d
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 32 symmetric eigensolver drivers that reduce to
xSTEBZ+xSTEIN(xSTEVX,xSYEVX,xHEEVX,xSPEVX,xHPEVX,xSBEVX,xHBEVX,xSTEVR,xSYEVR,xHEEVRand the_2STAGEvariants) scale a matrix whose largest entry lies outside[RMIN, RMAX]and undo the scaling at the end withThat block was copied from
xSTEV/xSYEV, where a nonzeroINFOcomes fromxSTEQR. In these drivers a nonzeroINFOat that point means either thatINFOeigenvectors failed to converge inxSTEIN(INFO <= N; the documentation says the eigenvalues are valid) or thatxSTEBZreturnedINFO - N. In the first caseW(INFO:M)are returned still multiplied bySIGMA, that is wrong by the factorRMIN/|T|orRMAX/|T|; in the secondIMAX = N + i - 1 > NandxSCALwrites past the end ofW. Two related robustness defects sit in the same drivers: when the matrix contains anInf,SIGMA = RMAX/ANRM = 0, so forRANGE = 'V'the scaled interval isVLL = VUU = 0andxSTEBZstops the process inXERBLA(parameter number 5 had an illegal value); andxSTEINrejects withINFO = -6(againXERBLA) or silently skips the negative block numbers with whichxSTEBZflags eigenvalues that did not converge, although every driver deliberately continues intoxSTEINwhenxSTEBZreturnsINFO = 1. This PR rescales allMeigenvalues, leaves a matrix with an infinite norm unscaled, and makesxSTEINacceptxSTEBZ's negative block numbers. 32 driver files plus{s,d,c,z}stein.f.Description
Rescaling. At the rescale label
INFOcan be nonzero only because ofxSTEBZorxSTEIN: a failure ofxSTERF,xSTEQRorxSTEMRearlier in the driver resetsINFOand falls through toxSTEBZ. Both leaveMvalid eigenvalues inW(xSTEBZwithINFO = 1returns its best estimates for the unconverged ones, withINFO = 2, 3theMit found, withINFO = 4none), so the right count isMin every case andIMAXdisappears. With finite input the wrong branch is reached whenever inverse iteration fails for at least one vector, which the drivers report asINFO = i <= N, "i eigenvectors failed to converge, their indices are stored in IFAIL", with no hint that eigenvaluesi .. Mare off by up to 1e146:Infinite norm.
SIGMAis now computed first and the scaling is applied only if it is positive, which fails exactly forANRM = Inf. An infinite matrix then takes the same path as one containing aNaN(which hasANRM = NaNand is never scaled): no abort,INFO = 0orN + 4depending onRANGE,NaNin the output, nothing written beyondW(M).Negative block numbers.
xSTEBZdocumentsIBLOCK(i) = -jfor an eigenvalue of blockjthat did not converge, and the drivers rely onxSTEINto accept that array (IF( IINFO.NE.1 ) GO TO 20skipsxSTEINfor every otherxSTEBZfailure).xSTEINnow usesABS( IBLOCK( . ) )at its four uses and documents it; the drivers then returnINFO = N + 1together with eigenvectors, as intended, instead of stopping inXERBLAwhen a negative entry follows a positive one or returning an untouchedZwhen the last entry is negative.Regression test.
xERRSTgets the second case:xSTEVXandxSYEVXwithRANGE = 'V'on a tridiagonal matrix whose first diagonal entry is an infinity, which must return without reachingXERBLA. On the parent commit both report the illegal fifth argument ofxSTEBZ, since the scale factorRMAX / ANRMis zero and both interval ends are scaled to zero; with the fix the matrix goes through unscaled and the drivers return.The other two changes cannot be reached from the test suite. A partial
xSTEINfailure needs a matrix whose norm is below the range the drivers document, and the negative block numbers needxSTEBZto returnINFO = 1, which finite in-range input cannot produce. The sweep below covers both.Validation
ctest: 100% of 215 tests passed).DSTEBZ+DSTEINdirect,DSTEVX,DSTEVR,DSYEVX,DSYEVR,ZHEEVX,DSPEVX,DSBEVX;NaN,+Inf,-Infin each of five positions ofDorE;n = 1, 2, 3, 5, 8;RANGE = 'A', 'V', 'I'), one process per case withXERBLAoverridden to exit, canaries around every output and work array and a 10 s timeout. Master: 280 aborts (every driver family,Inf,RANGE = 'V'), 200 cases writing three entries past the end ofW(RANGE = 'I',Inf,INFO = N + 4), 18 hangs. This branch: 0 aborts, 0 out-of-bounds writes, the same 18 hangs. The hangs areDSTEVR/DSYEVRwithRANGE = 'A'and aNaNin the first two rows, insideDLARRB(its bracket-widening loops have no iteration bound); they are an MRRR issue independent of this PR; Bound the bracket widening in xLARRB and xLARRJ so that a NaN matrix cannot hang xSTEMR #1389 fixes them.xSTEBZ's documented range (T * 2**512, wheree(j)**2overflows) aborted inDSTEINon master and returnsINFO = 1fromDSTEBZon this branch.[RMIN, RMAX]the drivers do not scale and this PR changes nothing; for a scaled matrix withINFO = 0the rescaling is the samexSCAL( M, ... )as before.