Let an exact zero always split the matrix in xBDSQR so that an Inf cannot hang it - #1390
Open
rmlarsen wants to merge 1 commit into
Open
Let an exact zero always split the matrix in xBDSQR so that an Inf cannot hang it#1390rmlarsen wants to merge 1 commit into
rmlarsen wants to merge 1 commit 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 |
…nnot hang it xBDSQR with singular vectors never returns for a bidiagonal matrix with an Inf in one of its first two diagonal entries, and neither do xBDSDC with COMPQ = 'I' and the drivers xGESVD and xGESDD built on them. The relative-accuracy threshold THRESH comes from a recurrence that evaluates Inf/Inf for such input and becomes NaN; the block search IF( ABSE.LE.THRESH ) then never recognises a split, while the convergence test zeroes E( M-1 ) and jumps back to the top of the loop without advancing the iteration counter. The two chase each other forever. Without vectors xBDSQR calls xLASQ1 instead and returns. Make an exact zero in E a split regardless of THRESH. For a finite THRESH, which is at least MAXITR*N*N*UNFL > 0, the extra clause is implied by the existing test, so nothing changes for finite input; for a NaN threshold every path back to the loop head now shrinks the block or advances ITER, and the routine returns within its iteration bound like it does for other non-finite input. xERRBD gets the case as a regression test: a 4 by 4 bidiagonal with an infinite D(1), passed to xBDSQR with both sets of 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; on the parent, all four xeigtst binaries stop responding at that call. Over 480 Inf/NaN cases (the four xBDSQR, DBDSDC, DGESVD, DGESDD, ZGESVD; five positions; n = 2, 3, 5, 8) the parent hangs in 66, this branch returns in all of them. 160 finite bidiagonal matrices give bit-identical singular values and vectors 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>
rmlarsen
force-pushed
the
bdsqr-zero-split
branch
from
September 8, 2026 03:23
e1c57a5 to
941f78c
Compare
This was referenced Sep 8, 2026
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. Difference from x86: the hang does not occur on this machine. gfortran compiles MIN to |
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
xBDSQRwith singular vectors never returns when the bidiagonal matrix contains anInfin one of its first two diagonal entries, and neither doxBDSDCwithCOMPQ = 'I'and the SVD driversxGESVDandxGESDDbuilt on them. The relative-accuracy thresholdTHRESHis computed from a recurrence that producesInf/Inf = NaNfor such input, and withTHRESH = NaNthe block searchIF( ABSE.LE.THRESH ) GO TO 80never recognises a split. The convergence test further down zeroesE( M-1 )and jumps back to the top of the loop without touching the iteration counter, the block search does not see the zero, the test zeroes it again, and so on. This PR makes an exact zero inEalways a split, which is what it already is for every finiteTHRESH(THRESH >= MAXITR*N*N*UNFL > 0), so the change is a no-op for finite input and closes the cycle for a NaN threshold. Four files:{s,d,c,z}bdsqr.f.Description
With
TOL >= 0(relative accuracy, the default), the threshold isFor
D(1) = Infthe first step evaluatesInf/Inf,SMINOAbecomes NaN and so doesTHRESH. In the main loop every test againstTHRESHis then false. The block search (DO 70) therefore never finds a split, soLL = 1and the whole matrix is the working block; the standard testABS( E( M-1 ) ).LE.ABS( TOL )*ABS( D( M ) )(or its forward/backward relative variants) zeroes an off-diagonal entry and executesGO TO 60, and the loop returns to the block search, which still does not see the zero.ITERis only advanced by a QR sweep, so the iteration limit is never reached. Without singular vectorsxBDSQRcallsxLASQ1instead and returns.The fix changes the split test to
IF( ABSE.LE.THRESH .OR. ABSE.EQ.ZERO ). For a finiteTHRESHthe added clause is implied by the first (THRESHis at leastMAXITR*N*N*UNFL, which is positive), so the singular values and vectors of every finite matrix are unchanged, bit for bit. For a NaNTHRESHevery path back to the loop head now makes progress: a split shrinks the working block, a 2-by-2 block is finished byxLASV2, and a QR sweep advancesITER, so the loop ends within the existingMAXITR*N*Nbound and the routine returns like it does for any other non-finite input.Minimal reproducer
Regression test.
xERRBDgets the case: a 4 by 4 bidiagonal with an infiniteD(1), passed toxBDSQRwith both sets of singular vectors, which must return without callingXERBLA. It belongs with the error-exit tests rather than inxCHKBD, because the routine returns no meaningful output for such a matrix and the point of the test is that it returns at all; a matrix type would additionally feed the infinity toxBDSDCandxBDSVDX, and its residual ratios would be NaN, which theRESULT( J ) .GE. THRESHcomparison passes silently. On the parent commit all fourxeigtstbinaries stop responding at that call, soLAPACK-xeigtst{s,d,c,z}_svd_infail by timeout; with the fix the four error-exit sections pass (SBD/DBD56 tests,CBD/ZBD39).Validation
ctest: 100% of 215 tests passed).DBDSQRwith and without vectors,SBDSQR,ZBDSQR,DBDSDC('I'),DGESVD,DGESDD,ZGESVD;+Inf,-Inf, NaN inD(1),D(2),D(N),E(1),E(N-1);N = 2, 3, 5, 8), one process per case with a 3 s timeout: master hangs in 66 (the threexBDSQRandDBDSDCfor+-InfinD(1)orD(2)withN >= 3, the three drivers for+-InfinD(1)); this branch returns in all 480 with the sameINFOvalues master produces for the non-finite cases it does return from (0, or the number of unconverged off-diagonals).N = 1 ... 40; random, scaled to 1e-150 and to 1e100, nearly diagonal with an exact zero off-diagonal) give bit-identical singular values and vectors fromDBDSQRon master and on this branch.xBDSQRwithout vectors takes thexLASQ1path; that path stops the process inXERBLA(DLASCL, parameter 4) for a NaN inD(N)and is not changed here; #1387 fixes it.