Keep a NaN inside its block in xSTEDC and report it through INFO - #1403
Open
rmlarsen wants to merge 2 commits into
Open
Keep a NaN inside its block in xSTEDC and report it through INFO#1403rmlarsen wants to merge 2 commits into
rmlarsen wants to merge 2 commits into
Conversation
xSTEDC splits the tridiagonal matrix into independent blocks where an off-diagonal entry is negligible against its diagonal neighbours: TINY = EPS*SQRT( ABS( D( FINISH ) ) )*SQRT( ABS( D( FINISH+1 ) ) ) IF( ABS( E( FINISH ) ).GT.TINY ) THEN extend the block Both comparisons are false for a NaN, so a NaN off-diagonal entry split the matrix and was dropped, and a NaN diagonal entry made TINY a NaN, ended the block before it and left it as a 1 by 1 block. The routine then returned INFO = 0 with a finite spectrum for a matrix with a NaN off the diagonal, or with the NaN reported as an eigenvalue and the coupling to its neighbours ignored. The other tridiagonal solvers, xSTEQR, xSTERF, xSTEBZ and xSTEMR, keep a NaN in its block. Extend the block unless the off-diagonal entry is known to be small, so that a NaN stays with its neighbours. A block that reaches the divide and conquer recursion is scaled by its max-norm with xLASCL, which stops in XERBLA for a NaN, so the norm is tested first and a NaN block is reported as a failure on that block, INFO = START*(N+1) + FINISH, the encoding the QR fallback already uses. A block small enough for xSTEQR gets the NaN through the QR iteration, which returns INFO > 0 or, for a 2 by 2 block, NaN eigenvalues. The same code sits in cSTEDC and zSTEDC; their COMPZ = 'I' path calls the real routine. xCHKST gets the case as a regression test for the divide and conquer path, which the sizes in sep.in (N <= 20, below SMLSIZ = 25) never reach: after the size and type loops it calls xSTEDC with COMPZ = 'I' and 'V' on an N = LDU tridiagonal matrix, once with a NaN in the middle of E and once with a NaN at the end of D, and reports INFO = 0 as a failure. On the parent commit all four calls return INFO = 0 in every precision; here they return INFO > 0. Over a NaN and Inf sweep of xSTEDC and xSTEVD (six positions, n = 1 to 64, every COMPZ and JOBZ, real and complex) the parent returned INFO = 0 for a NaN matrix with n >= 3 in 90 cases (DSTEDC 36, ZSTEDC 36, DSTEVD 18) and this branch in none of the 280 such cases; every finite case is bit-identical. The full LAPACK test suite passes: 0 numerical errors, 0 other errors, 80 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 #1403 +/- ##
==========================================
+ Coverage 69.36% 69.39% +0.03%
==========================================
Files 6122 6122
Lines 486337 486457 +120
Branches 23268 23268
==========================================
+ Hits 337330 337572 +242
+ Misses 148569 148447 -122
Partials 438 438
Continue to review full report in Codecov by Harness.
|
Use local arrays sized above SMLSIZ instead of treating LDU as the capacity of the caller's vectors and matrix columns. Preserve both NaN positions and both eigenvector options. Validation: 8 sep/se2 driver runs, 32 AddressSanitizer capacity cases, and 4 runs against the parent kernels that retain the expected failures.
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: The initial changes were prepared using Claude Code; review and subsequent test fixes used Codex.
Summary
xSTEDCsplits the tridiagonal matrix into independent blocks wherever an off-diagonal entry is negligible against its diagonal neighbours, and both halves of that test are false for a NaN: a NaN off-diagonal entry split the matrix and was dropped, and a NaN diagonal entry ended the block before it and was left as a 1 by 1 block. Whenever the block left over was large enough for the divide and conquer recursion (N > SMLSIZ = 25, which the sizes in the test inputs never reach), the routine returnedINFO = 0with a finite spectrum for a matrix with a NaN off the diagonal, or with the NaN reported as an eigenvalue and its coupling to the neighbours ignored. This PR keeps a NaN inside its block, asxSTEQR,xSTERF,xSTEBZandxSTEMRdo, and reports a NaN block throughINFOwith the encoding the QR fallback already uses. The numerical changes are in{s,d,c,z}stedc.f, with regression coverage in{s,d,c,z}chkst.f.Description
The block search is
A NaN in
E( FINISH )fails the comparison and the block ends there, without the entry; a NaN inD( FINISH+1 )makesTINYa NaN, and the block ends before it. The 1 by 1 block is skipped, so a NaN diagonal entry comes back as its own "eigenvalue" while the rest of the matrix is solved as if the neighbouring off-diagonal entries were zero. The complex routines have the same loop, and theirCOMPZ = 'I'path calls the real routine.Fix. The block is extended unless the off-diagonal entry is known to be small,
.NOT.( ABS( E( FINISH ) ).LE.TINY ), so a NaN stays with its neighbours. A block that reaches the recursion is scaled by its max-norm withxLASCL, which stops inXERBLAfor a NaN, so the norm is tested first and a NaN block is reported as a failure on that block,INFO = START*( N+1 ) + FINISH, the code that thexSTEQRfallback on a small block already returns, and that the documentation describes as the submatrix in rows and columnsINFO/(N+1)throughmod(INFO,N+1). A block small enough forxSTEQRgets the NaN through the QR iteration, which returnsINFO > 0from its iteration limit, or NaN eigenvalues for a 2 by 2 block. Finite matrices take exactly the path they took before: the negated comparison differs from the original only for a NaN.Minimal reproducer
INFO = 53 = 1*(N+1) + 26: the failure is reported on rows 1 through 26, the whole matrix. WithCOMPZ = 'N'the routine callsxSTERF, which already returnedINFO = 25for the same matrices.Regression test.
xCHKSTcallsxSTEDCwithCOMPZ = 'I'and'V', once with a NaN in the middle ofEand once with a NaN at the end ofD, and reportsINFO = 0as a failure. The case owns its diagonal, off-diagonal, and eigenvector arrays, withN = max(2, SMLSIZ+1)(26 for the defaultSMLSIZ = 25), so it reaches the divide and conquer path even when the input sizes are small.LDUis only a row stride and is not used to infer the capacity of caller arrays. On the parent commit all four calls returnINFO = 0in every precision; with the fix they returnINFO > 0.Validation
-O2 -fcheck=all; all 8 focusedsep.in/se2.indriver runs pass.SMLSIZ = 25and40, and verify that the dedicated regression leaves caller matrix and eigenvalue arrays untouched.Found while auditing the symmetric tridiagonal eigensolvers for the NaN and overflow handling of #1377-#1391.
Update: Codex review identified that LDU does not establish the capacity of the caller arrays. The regression now owns arrays sized above SMLSIZ; sanitizer checks cover small caller arrays and padded leading dimensions.