Guard against a NaN pivot in the packed, banded and tridiagonal Cholesky factorizations

[rmlarsen]

Disclaimer: This PR was prepared using Claude Code.

Summary

The packed (xPPTRF), banded (xPBTF2, and through it xPBTRF) and tridiagonal (xPTTRF) Cholesky factorizations test the candidate pivot with AJJ.LE.ZERO alone. A NaN is not .LE.ZERO, so a matrix containing a NaN is factored to a NaN factor and returned with INFO = 0; the drivers xPPSV, xPBSV and xPTSV then return a NaN solution as success. The dense routines xPOTF2, xPOTRF2 and xPSTF2 test AJJ.LE.ZERO.OR.DISNAN( AJJ ) and report the column. This PR adds the same DISNAN / SISNAN term 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:

	INFO
DPOTRF, ZPOTRF, DPSTRF	the NaN column
DPPTRF, DPBTRF, DPTTRF, ZPPTRF, ZPBTRF	0, factor full of NaN
DPOSV	the NaN column
DPPSV, DPBSV, DPTSV	0, NaN solution

Unlike 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.

xPBTRF is also inconsistent with itself. For KD <= 64 ILAENV returns NB = 1 and the whole factorization runs through xPBTF2, which misses the NaN; for KD > 64 the blocked path factors the diagonal blocks with xPOTF2, which catches it. The same matrix with the same NaN gets INFO = 0 or INFO = N depending on the bandwidth.

Fix. AJJ.LE.ZERO becomes AJJ.LE.ZERO.OR.DISNAN( AJJ ) at both sites of xPPTRF and xPBTF2, and D( I ).LE.ZERO becomes D( I ).LE.ZERO.OR.DISNAN( D( I ) ) at the six sites of xPTTRF (the 4-way unrolled loop, its remainder, and D( N )). SISNAN in the single-precision routines. Twelve files.

Minimal reproducer

program minimal
  implicit none
  integer, parameter :: n = 6
  double precision :: d(n), e(n-1), z
  integer :: info
  z = 0
  d = 4
  e = 1
  d(n) = z / z
  call dpttrf(n, d, e, info)
  print '(a,i0,a)', 'DPTTRF: info = ', info, '   (expected 6, as DPOTRF reports on the same matrix)'
end program

BEFORE (master):   DPTTRF: info = 0
AFTER (this branch): DPTTRF: info = 6

The fuller reproducer in the validation section exercises DPOTRF, DPPTRF, DPBTRF and DPTTRF side by side.

Validation

Storage-format sweep against dense xPOTRF, 6864 cases per build

xPPTRF (both UPLO), xPBTRF (both UPLO, KD = 0..3) and xPTTRF against xPOTRF on 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 (for KD < 3 the 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's INFO, the dense INFO, and an FNV-1a hash of the factor with NaNs hashed as a constant.

	cases	INFO differs from dense, master	INFO differs from dense, this branch
xPPTRF, finite	96	0	0
xPPTRF, NaN	1152	1152	0
xPBTRF, finite	384	0	0
xPBTRF, NaN	4608	4080	0
xPTTRF, finite	48	0	0
xPTTRF, NaN	576	576	0

(The 528 xPBTRF NaN 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 / DPTTRF on the same matrix with a NaN at (n, n):

BEFORE (master):            AFTER (this branch):
DPOTRF : info = 6           DPOTRF : info = 6
DPPTRF : info = 0           DPPTRF : info = 6
DPBTRF : info = 0           DPBTRF : info = 6
DPTTRF : info = 0           DPTTRF : info = 6

Regression test. Each of the three test paths gets a matrix type carrying a NaN on the last diagonal entry, with IZERO = N: ?PT type 13, ?PP type 10 and ?PB type 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 either UPLO and AB( KD+1, N ) or AB( 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 INFO with IZERO, which called ALAERH. ALAERH returns at its first statement when INFO is 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 unexpected INFO = 0 themselves and leave every other mismatch to ALAERH. On the parent commit the new types fail 6 times per precision for ?PT, 12 for ?PP and 88 for ?PB:

 *** DPTTRF returned INFO = 0 instead of     3 for N =    3, type 13
 *** DPPTRF returned INFO = 0 instead of     5 for UPLO = 'U', N =    5, type 10
 *** DPBTRF returned INFO = 0 instead of    10 for UPLO='U', N=   10, KD=    1, type  9

The SQRT( -ONE ) that builds the NaN uses a runtime variable, following xERRCXX, and the four xCHKPT files join the -Onopropagate list 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 f96546fc9 are the new ?PT type, the only one of the three that runs test ratios. That includes the ?PP, ?PB and ?PT routine and driver families in all four precisions (DPP 1332 + 1910, DPB 3458 + 4750, DPT 959 + 788 tests) and the _64 extended-API variants. Built with GCC 13.3, CMAKE_BUILD_TYPE=Release, BUILD_INDEX64_EXT_API=ON.

Checklist

• The documentation has been updated. (The INFO > 0 description, "the leading principal minor of order i is not positive definite", already covers a NaN; no interface changes.)
• If the PR solves a specific issue, it is set to be closed on merge. (No tracking issue; happy to open one.)

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 69.37%. Comparing base (a6c6e74) to head (edca5c7).
⚠️ Report is 2 commits behind head on master.
✅ All tests successful. No failed tests found.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1380      +/-   ##
==========================================
+ Coverage   69.36%   69.37%   +0.01%     
==========================================
  Files        6122     6123       +1     
  Lines      486711   486813     +102     
  Branches    23268    23272       +4     
==========================================
+ Hits       337584   337728     +144     
+ Misses     148689   148647      -42     
  Partials      438      438              

Components	Coverage Δ
BLAS	97.94% <ø> (ø)
CBLAS	96.98% <ø> (ø)
LAPACK	82.38% <100.00%> (ø)
LAPACKE	2.17% <ø> (ø)
TMGLIB	55.69% <ø> (ø)
BLAS testing	88.33% <ø> (ø)
CBLAS testing	89.63% <ø> (ø)
LAPACK testing	82.25% <100.00%> (+0.04%)	⬆️
LAPACKE testing	∅ <ø> (∅)

Files with missing lines	Coverage Δ
SRC/cpbtf2.f	97.87% <100.00%> (ø)
SRC/cpptrf.f	100.00% <100.00%> (ø)
SRC/cpttrf.f	96.72% <100.00%> (ø)
SRC/dpbtf2.f	97.67% <100.00%> (ø)
SRC/dpptrf.f	100.00% <100.00%> (ø)
SRC/dpttrf.f	95.65% <100.00%> (ø)
SRC/spbtf2.f	97.67% <100.00%> (ø)
SRC/spptrf.f	100.00% <100.00%> (ø)
SRC/spttrf.f	95.65% <100.00%> (ø)
SRC/zpbtf2.f	97.87% <100.00%> (ø)
... and 21 more

---

Continue to review full report in Codecov by Harness.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update a6c6e74...edca5c7. Read the comment docs.

[rmlarsen]

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.