Return INFO = 1 from xLALSD and xBDSDC instead of stopping in XERBLA when the bidiagonal contains a NaN - #1382
Open
rmlarsen wants to merge 1 commit into
Open
Conversation
❌ 4 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 |
This was referenced Sep 7, 2026
…when the bidiagonal contains a NaN xLALSD and xBDSDC scale the bidiagonal by its largest entry before the divide and conquer phase, guarded only by ORGNRM.EQ.ZERO. xLANST propagates a NaN, and xLASCL rejects CFROM = NaN as an illegal argument, so the call ends in XERBLA and the reference XERBLA stops the program. A NaN anywhere in A reaches xGELSD's bidiagonal, and so does an infinity, because the driver's scaling multiplies by zero and leaves NaN where the infinity was. xGESDD guards its own input against a NaN (Reference-LAPACK#469) but an infinity still reaches xBDSDC the same way. Every other least-squares driver returns from such input (xGELSS with INFO > 0 from xBDSQR, the rest with INFO = 0); xGELSD is the only one that kills the process. Test the norm with xISNAN at the scaling sites of xLALSD (including the N = 1 branch) and xBDSDC and return INFO = 1, the value xLASDA already uses when a singular value fails to converge, in all precisions. Document the new return in xLALSD, xBDSDC and xGELSD. xERRLS gets the case as a regression test: a 2 by 2 matrix with a NaN, passed to xGELSD, which must return INFO = 1 without calling XERBLA. The error-exit tests are where it belongs, since the driver returns no solution for such a matrix and the point of the test is that it returns at all; on the parent it reports the illegal xLASCL argument and returns INFO = -4 with a NaN solution. The added test changes nothing for finite input: over the finite cases of the sweep the solutions are bit-identical to the parent commit. Over 2296 (precision, m, n, NaN or infinity, position) xGELSD cases and 810 xBDSDC cases with a NaN, this branch returns from every call; xGELSD returns INFO = 1 in every case, xBDSDC in every divide and conquer case, and no call reaches XERBLA. The full LAPACK test suite passes: 5441901 LAPACK and 315872 BLAS tests, 0 numerical errors, 0 other errors, the same totals as the parent. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
rmlarsen
force-pushed
the
dc-nan-info-not-xerbla
branch
from
September 8, 2026 03:20
9e8fed1 to
24e080d
Compare
rmlarsen
added a commit
to rmlarsen/lapack
that referenced
this pull request
Sep 8, 2026
…is a NaN xLASQ1 scales the bidiagonal matrix by its largest entry SIGMX before running dqds. When D or E contains a NaN that survives the MAX reductions computing SIGMX (with gfortran, a NaN in D(N)), SIGMX is a NaN and xLASCL rejects it by stopping the process in XERBLA. xBDSQR without singular vectors and xGESVD with JOBU = JOBVT = 'N' take this path; a NaN elsewhere in D returns INFO = 0 with NaN output. Test SIGMX with xISNAN and return INFO = 1, the mechanism Reference-LAPACK#1382 uses for xLALSD and xBDSDC. xBDSQR treats every nonzero INFO from xLASQ1 as a request to finish with the QR algorithm, so it then reports the non-convergence on the NaN data through its own INFO like it does for other non-finite input. Finite input never takes the new branch. xERRBD gets the case as a regression test: a 4 by 4 bidiagonal with a NaN in D(N), passed to xBDSQR without singular vectors, which must return without calling XERBLA. The error-exit tests are where it belongs, since the routine returns no meaningful output for such a matrix and the point of the test is that it returns at all. All four xERRBD files carry it, because the complex xBDSQR takes the same path through the real xLASQ1. Over 120 NaN/Inf cases of DBDSQR without vectors and DGESVD the parent stops in XERBLA 3 times and this branch never; every other case returns the same INFO on both. The full LAPACK test suite passes: 5441901 LAPACK tests, 0 numerical errors, 0 other errors, the same totals as the parent. 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
xLALSD(the solver behindxGELSD) andxBDSDCscale the bidiagonal by its largest entry before the divide and conquer phase, guarded only byORGNRM.EQ.ZERO.xLANSTpropagates aNaN, andxLASCLrejectsCFROM = NaNas an illegal argument, so the call ends inXERBLA, and the referenceXERBLAstops the program. ANaNanywhere inAreaches the bidiagonal ofxGELSD, and so does an infinity, because the driver's own scaling multiplies by zero and leaves aNaNwhere the infinity was. This PR tests the norm withxISNANat those sites and returnsINFO = 1, the valuexLASDAalready uses when a singular value fails to converge. Nothing changes for finite input.Description
With a
NaNatA(2,2)of a 6-by-4 least-squares problem, one process per call:DGELS,DGELST,DGETSLS,DGELSYINFO = 0,NaNsolutionDGELSSINFO = 3(DBDSQRdid not converge)DGELSD,ZGELSD** On entry to DLASCL parameter number 4 had an illegal valueThe backtrace is
DGELSD->DLALSD->DLASCL->XERBLA, from the unconditionalDLASCL( 'G', 0, 0, ORGNRM, ONE, ... )afterORGNRM = DLANST( 'M', N, D, E ). An infinity inAends the same way:ANRM > BIGNUMmakesDGELSDcallDLASCL( 'G', 0, 0, ANRM, BIGNUM, ... ), which multiplies by zero and turns the infinity into aNaN.xGESDDguards its own input against aNaN(returnsINFO = -4, #469), but an infinity passes that guard, becomes aNaNin the same way, and stops the process insidexBDSDC. Callers that replaceXERBLAwith an exception see an illegal-argument error fromDLASCLfor what is a data problem.Whether a
NaNis detected in a scan that usesMAXand comparisons is compiler dependent; the sites changed here usexISNAN, which is not.Fix. At the scaling sites of
{s,d,c,z}lalsd.f(theN = 1branch and the general one) and{s,d}bdsdc.f, addELSE IF( xISNAN( ... ) ) THEN INFO = 1; RETURN. TheINFOdescriptions ofxLALSD,xBDSDCandxGELSDdocument the new return. Ten files, six with code changes.xBDSDCreaches the changed site only on its divide and conquer path (COMPQ = 'I'or'P'andN > SMLSIZ). Its small-problem path and its singular-values-only path go throughxLASDQas before and are not changed by this PR: forN <= 2they returnINFO = 0with aNaNresult, for largerNa non-convergence count.Minimal reproducer
Validation
NaN and infinity at every position, all precisions, both sides of the divide and conquer threshold
repro/nan_sweep.f90runs{s,d,c,z}GELSDon shapes (1,1), (2,2), (5,3), (3,5), (25,25), (26,26), (40,25), (25,40) and (80,60) with aNaNor an infinity at every position (a lattice of 40 positions for the larger shapes), and{s,d}BDSDCwith aNaNon every diagonal and off-diagonal entry for the same sizes.XERBLAis replaced by a routine that stops the process with status 99, so a completed run proves that no case reached it. The same program with finite input prints the solutions in hex for a bit-for-bit comparison with master.INFO = 1SGELSD,DGELSD,CGELSD,ZGELSDSBDSDC,DBDSDC,N > 25SBDSDC,DBDSDC,N <= 25xLASDQbehavior (see above)Finite input: 27 (routine, shape) cases, bit-identical between master and this branch. Master stops on the first
NaNcase.Regression test.
xERRLSgets the case: a 2 by 2 matrix with a NaN in one entry, passed toxGELSD, which must returnINFO = 1without callingXERBLA. It belongs with the error-exit tests rather than inxDRVLS, because the driver returns no solution for such a matrix and the point of the test is that it returns at all; a NaN matrix type inxDRVLSwould report the fixed behaviour as a failure (INFO .NE. 0goes toALAERH) and let the master behaviour through, since its residual ratios are NaN andRESULT( J ) .GE. THRESHis false for a NaN. On the parent commit the fourxerrlssections report*** xGELSD on a matrix with a NaN returned INFO = -4 instead of 1 ***preceded by thexLASCLillegal-argument line, and*** xLS routines failed the tests of the error exits ***; with the fix they pass.xGELSDneeds a real workspace of about 800 for a 2 by 2 problem, so the arrays inxERRLSgrow accordingly, and the four files join the-Onopropagatelist that already keeps the NAG compiler from folding theSQRT( -ONE )that builds the NaN.Test suite. The full LAPACK test suite passes on this branch: 215 of 215 CTest entries, 5441901 LAPACK tests and 315872 BLAS tests with 0 numerical errors and 0 other errors, the same totals as the parent commit
f96546fc9built and run the same way (GCC 13.3,CMAKE_BUILD_TYPE=Release,BUILD_INDEX64_EXT_API=ON).Performance. Timed on a 13th Gen Intel(R) Core(TM) i7-13700HX under WSL2 with the reference BLAS, GCC 13.3,
-O2. To separate the change from code-placement effects (which move untouched routines by up to 28% between two separately linked static libraries on this machine), the parent library is a shared object shared by both sides, and each benchmark binary carries its own copy of only the changed routines, parent or branch, which interposes over the library's; everything else is byte-identical. Four rounds in alternating order, one core, one process per run, on an idle machine; medians of the per-round medians, withDPOTRFas an untouched control. The benchmark driver and raw output are available on request.DGELSD(us/call)DBDSDC(us/call)DPOTRFcontrol (us/call)No measurable difference; the change is one
xISNANtest per call.Not changed here.
xBDSQRwith singular vectors loops indefinitely on an infinity inD(reachable throughxBDSDCwithN <= SMLSIZ);xGELSDcannot pass an infinity down because its scaling turns it into aNaNfirst. That is a separate defect.Checklist
INFOdescriptions ofxLALSD,xBDSDCandxGELSD.)