Scale ABSTOL with the matrix in xSTEVX and xSTEVR - #1386
Open
rmlarsen wants to merge 1 commit into
Open
Conversation
❌ 5 Tests Failed:
View the top 3 failed test(s) by shortest run time
View the full list of 77 ❄️ flaky test(s)
To view more test analytics, go to the Test Analytics Dashboard |
xSTEVX and xSTEVR scale the tridiagonal matrix into [RMIN, RMAX] before calling xSTEBZ but, unlike the other 28 drivers with the same scaling block, pass the caller's ABSTOL unscaled. A positive ABSTOL is then applied to a matrix multiplied by SIGMA = RMAX / |T| whenever |T| > RMAX = SAFMIN**(-1/4), which is 1.6e77 in double precision and 3e9 in single, so bisection stops as soon as an interval is narrower than ABSTOL / SIGMA in the caller's units, a width that can exceed the norm of the scaled matrix. SSTEVX on tridiag(1e12, 5e11, 1e12) with ABSTOL = 1e-6 |T| returns eigenvalues with relative errors of 2e-4, 3e-2 at 1e14, and DSTEVX at |T| = 1e100 with ABSTOL = 1e-10 |T| returns 0.4, all with INFO = 0, while SSYEVX and DSYEVX honour the tolerance. Below RMIN the tolerance is only applied more tightly than requested. Scale a positive ABSTOL with the matrix and pass ABSTLL = MAX( MIN( ABSTOL, |T| )*SIGMA, SAFMIN ) to xSTEBZ. The plain product of xSYEVX is not enough: for the documented "most accurate" ABSTOL = 2*SAFMIN it underflows to zero once |T| > RMAX/EPS, and xSTEBZ reads a zero tolerance as a request for its default EPS*|T|, so DSTEVX on D = (1e100, 1), E = 5e49 would return 1.1e83 for the eigenvalue 0.75. Every positive tolerance below PIVMIN >= SAFMIN is equivalent in xLAEBZ, so the floor changes nothing that did not underflow; a tolerance above |T| is met by the Gershgorin interval, so the cap changes nothing that could not overflow when a tiny matrix is scaled up. Nothing changes for ABSTOL <= 0 or for a matrix the driver does not scale. xDRVST gets a matrix type for it: type 19 is a tridiagonal matrix with equal diagonal entries and half that off the diagonal, scaled to the "large" magnitude of the existing types, which is above the norm the drivers scale down from, and it is the only type that asks for a positive ABSTOL, ten ulps of the matrix norm. Every other type passes 2*SAFMIN, which is inert in the bisection and so cannot show the defect, and the tridiagonal drivers were tested on diagonal matrices alone, where each block is one by one and the tolerance never matters. On the parent commit the type fails 18 ratios in single precision and 40 in double. With the change xSTEVX and xSYEVX agree to the last digit on the matrices above at every scale, and ABSTOL = 2*SAFMIN gives the same eigenvalues as on master. The full LAPACK test suite passes: 5449821 LAPACK tests, 0 numerical errors, 0 other errors; the 7920 tests above the parent are the new type. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
rmlarsen
force-pushed
the
stevx-abstol-scaling
branch
from
September 8, 2026 04:19
7885655 to
b364dad
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
xSTEVXandxSTEVRscale the tridiagonal matrix into[RMIN, RMAX]before callingxSTEBZbut, unlike the other 28 drivers with the same scaling block, pass the caller'sABSTOLunscaled. A positiveABSTOLis then applied to a matrix that has been multiplied bySIGMA = RMAX/|T|when|T| > RMAX = SAFMIN**(-1/4)(1.6e77 in double precision, 3e9 in single), or byRMIN/|T|when|T| < RMIN. AboveRMAXbisection therefore stops as soon as an interval is narrower thanABSTOL, which in the scaled units isABSTOL/SIGMAand can exceed the norm of the scaled matrix:SSTEVXontridiag(1e12, 5e11, 1e12)withABSTOL = 1e-6 |T|returns eigenvalues with relative errors of 2e-4 (3e-2 at 1e14),DSTEVXat|T| = 1e100withABSTOL = 1e-10 |T|returns 0.4, all withINFO = 0;SSYEVX/DSYEVXon the same matrices honour the tolerance. BelowRMINthe tolerance is only applied more tightly than requested. This PR scalesABSTOLwith the matrix in the four tridiagonal drivers, asABSTLL = MAX( MIN( ABSTOL, |T| )*SIGMA, SAFMIN ), so that the scaled tolerance can neither underflow to zero nor overflow.Description
The scaling in these drivers exists so that the caller does not have to (the drivers'$s\,(1 + \cos(k\pi/7))$ , with
RMAXis exactly the largest entryxSTEBZdocuments as safe), so a matrix with entries above 3e9 in single precision or 1.6e77 in double is an ordinary input, and the routines' contract is that each eigenvalue is found to withinABSTOL(plusEPS |T|). Maximum relative eigenvalue error forT = s * tridiag(1, 0.5, 1),n = 6, whose eigenvalues areABSTOL = rel * |T|:sABSTOL / |T|xSTEVXmasterxSTEVXthis branchxSYEVXxSTEVRreachesxSTEBZforRANGE = 'V'and wheneverxSTEMRfails, with the same unscaledABSTOL; it is fixed the same way. ItsTRYRACdecision (ABSTOL .LE. 2 N EPS) is left as inxSYEVR.The plain
ABSTLL = ABSTOL*SIGMAofxSYEVXis not enough. For the documented "most accurate"ABSTOL = 2*SAFMINthe product underflows to zero once|T|exceeds roughlyRMAX/EPS(1e93 in double precision, 5e16 in single), andxSTEBZreads a zero tolerance as a request for its defaultEPS*|T|, which for a graded matrix is far worse than the unscaled2*SAFMINmaster passes:DSTEVXonD = (1e100, 1),E = 5e49(eigenvalues1e100and0.75) withABSTOL = 2*DLAMCH('S')returns0.75on master and1.1e83with the plain product,INFO = 0(SSTEVXonD = (1e20, 1),E = 5e9:1.8e12). InxLAEBZan interval has converged once its width is belowMAX( ABSTOL, PIVMIN, RELTOL*|lambda| ), andxSTEBZsetsPIVMIN >= SAFMIN, so every positive tolerance belowSAFMINacts exactly likeSAFMIN: the floor changes nothing that did not underflow. When a tiny matrix is scaled up (|T| < RMIN,SIGMAup to about 1e154) the product can instead overflow, which traps under-ffpe-trap=overflow; a tolerance above|T|is already met by the Gershgorin interval, so cappingABSTOLat|T|before scaling changes nothing that could not overflow.xSYEVXand the other 27 drivers with theABSTLL = ABSTOL*SIGMAblock have both defects on master (DSYEVXreturns1.1e83on the matrix above); they are left for a separate PR.Minimal reproducer
Regression test.
xDRVSTgets a matrix type 19: a tridiagonal matrix with equal diagonal entries and half that off the diagonal, at the "large" magnitude the existing types already use, which is above the norm from which the drivers scale down. It is the only type that asks for a positiveABSTOL, ten ulps of the matrix norm.Nothing in the suite could have caught this. Every type passes
ABSTOL = 2*SAFMIN, which is at or below the smallest pivot and therefore inert insidexLAEBZ, and the tridiagonal drivers were exercised only on the diagonal matrices of types 1 to 7, where every block of the bisection is one by one and the tolerance never enters. The new type lifts that restriction for itself alone, soxSTEV,xSTEVX,xSTEVRandxSTEVDsee a genuine tridiagonal for the first time.On the parent commit the type fails 18 ratios in single precision and 40 in double, comparing
xSTEVXandxSTEVRagainstxSTEV; with the fix all of them pass.Validation
ctest: 100% of 215 tests passed). The 7920 tests above the parent commit are the new type.xSTEVXandxSYEVXagree to the last digit at every scale.ABSTOL = 2*SAFMINon the graded matrices above:DSTEVX(RANGE = 'A'and'V'),DSTEVR('V'),SSTEVXandSSTEVRreturn0.75as on master. Under-ffpe-trap=overflow,D = (1e-300, 1e-301),E = 1e-301withABSTOL = 1e200runs through the four drivers without a trap.ABSTOL <= 0, and for any matrix the driver does not scale, the call toxSTEBZis unchanged.