Drop the absolute floor from the splitting test in xSTEBZ - #1385
Open
rmlarsen wants to merge 1 commit into
Open
Conversation
❌ 10 Tests Failed:
View the top 3 failed test(s) by shortest run time
View the full list of 75 ❄️ flaky test(s)
To view more test analytics, go to the Test Analytics Dashboard |
xSTEBZ declares an off-diagonal entry negligible when e(j)**2 < ulp**2 |d(j) d(j+1)| + SAFEMN. The SAFEMN term makes the test absolute: any |e(j)| below sqrt(SAFEMN) (1.5e-154 in double precision) splits the matrix whatever its neighbours are. The drivers built on xSTEBZ (xSTEVX, xSYEVX, xHEEVX, xSPEVX, xHPEVX, xSBEVX, xHBEVX, and xSTEVR/xSYEVR/xHEEVR for RANGE = 'V') scale a small matrix up only to RMIN = sqrt(SAFMIN/ULP), about 1e-146, where that floor corresponds to |e(j)|/|T| < sqrt(ULP) = 1.5e-8. So for a matrix whose largest entry is below sqrt(SAFMIN)/ULP (6.7e-139) the drivers drop off-diagonals of relative size between ULP and 1.5e-8: the eigenvalues are wrong by up to 1.5e-8 |T| and the eigenvectors have residuals of the same size, with INFO = 0, while xSTEV and xSYEV return them to working precision. The test suite never reaches this regime; its smallest matrix scale is RTUNFL*N/ULP, exactly where the floor stops mattering. Use Kahan's relative criterion e(j)**2 < ulp**2 |d(j) d(j+1)|, or e(j) = 0, as xLARRA already does, with the threshold formed as (|d(j)| ulp)(|d(j+1)| ulp) so that it cannot overflow before e(j)**2 does. For a matrix in the normal range this changes the treatment only of an off-diagonal below 1.5e-154 sitting between two diagonal entries below 2e-138, which now stays in its block. Also correct "underflow" to "overflow" in the scaling advice of the xLAEBZ documentation. For T = s tridiag(1, r, 1) with the eigenvalues known in closed form, every DSTEBZ-based driver returns eigenvalue errors of r |T| for r between 1e-14 and 1e-8 at s <= 1e-140 on the parent commit and errors below n ulp |T| on this branch, at every scale. xCHKST gets a matrix type for it: type 22 is a tridiagonal matrix with equal diagonal entries and an off-diagonal entry one tenth of the square root of an ulp of them, at the smallest norm the drivers scale a matrix up to. Its off-diagonal square underflows, so the absolute term of the old test hid it, while the relative test does not. The existing comparison of xSTEBZ against xSTERF, and the residual of the eigenvectors xSTEIN builds from its blocks, both catch the difference: on the parent commit the type fails 7 ratios in single precision and 8 in double at ratios of 1e6 and above, against a threshold of 60. The 'I' against 'V' comparison is skipped for the type, because the interval the test builds is floored at twice the square root of the underflow threshold, which is wider than the whole spectrum here. The complex checkers carry the same type, since they call the real xSTEBZ. A power-of-two twin sweep of DSTEBZ + DSTEIN over 42 exponents from -1070 to 1020 is identical to the parent wherever the parent was correct. The full LAPACK test suite passes: 5446341 LAPACK tests, 0 numerical errors, 0 other errors; the 4440 tests above the parent are the new type. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
rmlarsen
force-pushed
the
stebz-relative-split
branch
from
September 8, 2026 04:20
12b5af1 to
66f0863
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 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
xSTEBZtreats an off-diagonal entry as negligible whene(j)**2 < ulp**2 * |d(j) d(j+1)| + SAFEMN. TheSAFEMNterm makes the test absolute: any|e(j)|belowsqrt(SAFEMN)(1.5e-154 in double precision) splits the matrix whatever its neighbours are. The drivers that callxSTEBZ(xSTEVX,xSYEVX,xHEEVX,xSPEVX,xHPEVX,xSBEVX,xHBEVX, andxSTEVR/xSYEVR/xHEEVRforRANGE = 'V'or after anxSTEMRfailure) scale a small matrix up only toRMIN = sqrt(SAFMIN/ULP)(1e-146), where that floor corresponds to|e(j)|/|T| < sqrt(ULP)= 1.5e-8. So for every matrix whose largest entry is belowsqrt(SAFMIN)/ULP(6.7e-139) the drivers silently drop off-diagonal entries of relative size betweenULPand 1.5e-8: the eigenvalues are wrong by up to1.5e-8 * |T|and the eigenvectors have residuals of the same size, all withINFO = 0, whilexSTEV/xSYEV(QR) return them to working precision. This PR replaces the test by Kahan's relative criterione(j)**2 < ulp**2 * |d(j) d(j+1)|(ore(j) = 0), which is whatxLARRAalready uses for the MRRR drivers, and fixes an "underflow"/"overflow" typo in thexLAEBZdocumentation. Two code files, two documentation files.Description
For$s\,(1 + 2r\cos(k\pi/7))$ exactly. Maximum eigenvalue error divided by
T = s * tridiag(1, r, 1)of order 6 the eigenvalues aren * ulp * |T|on master;DSTEVX(A)uses the documented "most accurate"ABSTOL = 2*DLAMCH('S'):srDSTEVDSTEVX(A)DSTEVX(I)DSTEVX(V)DSTEVR(V)DSYEVX(I)DSYEVR(V)ZHEEVX(I)DSPEVX(I)DSBEVX(I)The eigenvector residuals$\|Tz - wz\| / (\|T\|\,\|z\|\,n\,\mathrm{ulp})$ are the same numbers to within 25%. With this branch every entry of the table is below 0.25. At
s = 1e-140the drivers do not scale and the floor still discards|e| = 1e-154; atr = 1e-6, and at everyrfors >= 1e-138, master is already accurate, which locates the affected region exactly at|e(j)| < sqrt(SAFEMN)next to diagonal entries belowsqrt(SAFEMN)/ULP.The test suite cannot see this: the smallest matrix scale in
xCHKST/xDRVSTisRTUNFL*N*ULPINV = 6.7e-139 * N, precisely the scale at which the floor stops mattering.Fix. The threshold is formed as
(|d(j)| ulp) (|d(j+1)| ulp), which is bit-identical to|d(j) d(j+1)| ulp**2whenever the latter is finite and does not overflow beforee(j)**2does, so an off-diagonal next to two entries above 1.3e154 is no longer declared negligible by an infinite threshold (that case is outside the rangexSTEBZdocuments, but the old test turned it into a wrong answer withINFO = 0).e(j) = 0splits as before. For a matrix in the normal range the only entries whose treatment changes are off-diagonals below 1.5e-154 sitting between two diagonal entries below 2e-138: they now stay in their block instead of being zeroed, their squares take part in the Sturm counts as subnormal numbers, and the twin sweep below shows the counts unaffected. A nonzeroe(j)whose square underflows to zero still splits, which the Sturm sequence could not distinguish from zero anyway.Minimal reproducer
The remaining differences are one
ulpof|T|.Regression test.
xCHKSTgets a matrix type 22: a tridiagonal matrix with equal diagonal entries and an off-diagonal entry one tenth of the square root of an ulp of them, atsqrt( SAFMIN / EPS ), the smallest norm the eigenvalue drivers scale a matrix up to. The square of that off-diagonal entry underflows, which is exactly what the absolute term of the old test hid and what the relative test sees.Two existing ratios catch it, the comparison of
xSTEBZagainstxSTERFand the residual of the eigenvectorsxSTEINbuilds from its blocks: on the parent commit the type fails 7 ratios per run in single precision and 8 in double, at 1e6 and above against a threshold of 60, and with the fix all of them pass. The'I'against'V'comparison is skipped for this type, because the interval that test builds is floored at twice the square root of the underflow threshold, which is wider than the whole spectrum of such a matrix, so the two ranges cannot agree either way. All fourxCHKSTfiles carry the type, since the complex ones call the realxSTEBZas well.Adding the type raises the number of types in
sep.inpast the count that madeALAREQread its list, soMAXTYPfor the symmetric problem goes to 22 in all fourxCHKEEfiles; the two-stage checker clamps to its own 21 and ignores the new type.Validation
ctest: 100% of 215 tests passed). The 4440 tests above the parent commit are the new type.DSTEBZcalls at|T| = 1e-290, outside the rangexSTEBZdocuments (e(j)**2underflows to zero there), which the drivers handle by scaling.DSTEBZ+DSTEINonT * 2**kfor 42 values ofkfrom -1070 to 1020, compared line by line with master: eigenvalues, block structure,INFOand eigenvectors are identical for everykat which master was correct (-508 <= k <= 507); atk = -510 ... -520master splits the matrix completely and returns the diagonal, this branch returns the eigenvalues to working precision; atk = 512, 513master aborts inDSTEIN(XERBLA, illegalIBLOCK), this branch returnsINFO = 1fromDSTEBZlike master does for largerk.Found while auditing the tridiagonal bisection and inverse-iteration routines; the routines themselves are correct within their documented range, this is the one place where the drivers' scaling and
xSTEBZ's splitting disagree.