Guard against a NaN pivot in the packed, banded and tridiagonal Cholesky factorizations - #1380
Open
rmlarsen wants to merge 2 commits into
Open
Guard against a NaN pivot in the packed, banded and tridiagonal Cholesky factorizations#1380rmlarsen wants to merge 2 commits into
rmlarsen wants to merge 2 commits into
Conversation
❌ 10 Tests Failed:
View the top 3 failed test(s) by shortest run time
View the full list of 74 ❄️ flaky test(s)
To view more test analytics, go to the Test Analytics Dashboard |
…sky factorizations xPOTF2, xPOTRF2 and xPSTF2 reject a pivot with IF( AJJ.LE.ZERO.OR.DISNAN( AJJ ) ) THEN Their packed, banded and tridiagonal counterparts xPPTRF, xPBTF2 and xPTTRF test AJJ.LE.ZERO (or D( I ).LE.ZERO) alone. A NaN is not .LE.ZERO, so a matrix containing one is factored into a NaN factor and returned with INFO = 0, and the drivers xPPSV, xPBSV and xPTSV return a NaN solution as success, where xPOSV reports the column. xPBTRF is also inconsistent with itself: for KD <= 64 ILAENV selects the unblocked xPBTF2 and the NaN passes; for KD > 64 the blocked path factors the diagonal blocks with xPOTF2 and catches it. No index depends on this test, so unlike the packed Bunch-Kaufman case (Reference-LAPACK#1378) the routines stay in bounds; the defect is the silent INFO = 0. Add the DISNAN / SISNAN term at both sites of xPPTRF and xPBTF2 and at the six sites of xPTTRF, in all four precisions. The three test paths get a matrix type for it: PT 13, PP 10 and PB 9, each the generated matrix of the preceding type with a NaN written on its last diagonal entry and IZERO set to N, so that the existing check of INFO against IZERO covers it. That check reached ALAERH, which returns without a message when INFO is zero, so the case it has to report was the one it dropped; the three checkers now report an unexpected INFO = 0 themselves and leave the rest to ALAERH. On the parent commit the new types give 6 failures per precision for PT, 12 for PP and 88 for PB. Behavior is unchanged on matrices that contain no NaN: over a sweep of 6864 (precision, format, UPLO, KD, n, NaN position) cases the 528 finite-input factors are bit-identical to the parent commit, and INFO now agrees with xPOTRF on the same matrix in every NaN case, where before 5808 of 6336 disagreed. The full LAPACK test suite passes: 5441925 LAPACK and 315872 BLAS tests, 0 numerical errors, 0 other errors, 24 tests more than the parent commit from the new PT type. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
rmlarsen
force-pushed
the
cholesky-nan-guard
branch
from
September 8, 2026 03:54
30b57b3 to
c09328c
Compare
xCHKPP and xCHKPB build their NaN type the same way xCHKPT does, with SQRT of a negative variable, but only xCHKPT was listed, and the comment above the _64 branch lost its indentation. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
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
The packed (
xPPTRF), banded (xPBTF2, and through itxPBTRF) and tridiagonal (xPTTRF) Cholesky factorizations test the candidate pivot withAJJ.LE.ZEROalone. A NaN is not.LE.ZERO, so a matrix containing a NaN is factored to a NaN factor and returned withINFO = 0; the driversxPPSV,xPBSVandxPTSVthen return a NaN solution as success. The dense routinesxPOTF2,xPOTRF2andxPSTF2testAJJ.LE.ZERO.OR.DISNAN( AJJ )and report the column. This PR adds the sameDISNAN/SISNANterm to the twelve packed, banded and tridiagonal routines. Nothing changes for a matrix without a NaN.Description
With a NaN on the last diagonal entry of an otherwise SPD matrix:
INFODPOTRF,ZPOTRF,DPSTRFDPPTRF,DPBTRF,DPTTRF,ZPPTRF,ZPBTRF0, factor full of NaNDPOSVDPPSV,DPBSV,DPTSV0, NaN solutionUnlike the corresponding gap in the packed Bunch-Kaufman routines (#1378), no index depends on this test, so the routines stay in bounds; the defect is the silent
INFO = 0.xPBTRFis also inconsistent with itself. ForKD <= 64ILAENVreturnsNB = 1and the whole factorization runs throughxPBTF2, which misses the NaN; forKD > 64the blocked path factors the diagonal blocks withxPOTF2, which catches it. The same matrix with the same NaN getsINFO = 0orINFO = Ndepending on the bandwidth.Fix.
AJJ.LE.ZERObecomesAJJ.LE.ZERO.OR.DISNAN( AJJ )at both sites ofxPPTRFandxPBTF2, andD( I ).LE.ZERObecomesD( I ).LE.ZERO.OR.DISNAN( D( I ) )at the six sites ofxPTTRF(the 4-way unrolled loop, its remainder, andD( N )).SISNANin the single-precision routines. Twelve files.Minimal reproducer
The fuller reproducer in the validation section exercises
DPOTRF,DPPTRF,DPBTRFandDPTTRFside by side.Validation
Storage-format sweep against dense xPOTRF, 6864 cases per build
xPPTRF(bothUPLO),xPBTRF(bothUPLO,KD= 0..3) andxPTTRFagainstxPOTRFon the same SPD or HPD matrix, four precisions,n= 1..12, with the matrix unperturbed, with a NaN on each diagonal entry in turn, and with a NaN on each sub-diagonal pair in turn (forKD < 3the banded routine and its dense reference both see the matrix truncated to the band, so a NaN outside it is dropped by both). Each line records the variant'sINFO, the denseINFO, and an FNV-1a hash of the factor with NaNs hashed as a constant.INFOdiffers from dense, masterINFOdiffers from dense, this branchxPPTRF, finitexPPTRF, NaNxPBTRF, finitexPBTRF, NaNxPTTRF, finitexPTTRF, NaN(The 528
xPBTRFNaN cases that already agreed on master are the ones where the NaN falls outside the band.) All 528 finite-input lines, factor hashes included, are byte-identical between master and this branch.The fuller reproducer,
DPOTRF/DPPTRF/DPBTRF/DPTTRFon the same matrix with a NaN at(n, n):Regression test. Each of the three test paths gets a matrix type carrying a NaN on the last diagonal entry, with
IZERO = N:?PTtype 13,?PPtype 10 and?PBtype 9. The matrix is generated exactly as for the preceding type and one entry is overwritten, which is the last entry of the packed array for eitherUPLOandAB( KD+1, N )orAB( 1, N )in band storage. No test ratio is computed for it, as for the existing zero-row types; the factorization has to report the pivot.That report goes through the existing comparison of
INFOwithIZERO, which calledALAERH.ALAERHreturns at its first statement whenINFOis zero, so the one outcome the new type has to catch was the one it dropped silently, and the same is true of the existing zero-row types. The three checkers now print an unexpectedINFO = 0themselves and leave every other mismatch toALAERH. On the parent commit the new types fail 6 times per precision for?PT, 12 for?PPand 88 for?PB:The
SQRT( -ONE )that builds the NaN uses a runtime variable, followingxERRCXX, and the fourxCHKPTfiles join the-Onopropagatelist that keeps the NAG compiler from folding it.Test suite. The full LAPACK test suite passes on this branch: 215 of 215 CTest entries, 5441925 LAPACK tests and 315872 BLAS tests with 0 numerical errors and 0 other errors. The 24 tests above the parent commit
f96546fc9are the new?PTtype, the only one of the three that runs test ratios. That includes the?PP,?PBand?PTroutine and driver families in all four precisions (DPP1332 + 1910,DPB3458 + 4750,DPT959 + 788 tests) and the_64extended-API variants. Built with GCC 13.3,CMAKE_BUILD_TYPE=Release,BUILD_INDEX64_EXT_API=ON.Checklist
INFO > 0description, "the leading principal minor of order i is not positive definite", already covers a NaN; no interface changes.)