Bound the bracket widening in xLARRB and xLARRJ so that a NaN matrix cannot hang xSTEMR by rmlarsen · Pull Request #1389 · Reference-LAPACK/lapack

rmlarsen · 2026-09-07T22:45:10Z

Disclaimer: This PR was prepared using Claude Code.

Summary

xLARRB refines eigenvalues of a relatively robust representation $LDL^T$ by bisection, and xLARRJ does the same on the tridiagonal matrix itself when only eigenvalues are wanted. Before bisecting each makes sure the initial interval [W - WERR, W + WERR] contains the wanted eigenvalue by widening it: it moves the endpoint by BACK, doubles BACK, and repeats until the Sturm count agrees with the index. In exact arithmetic that ends quickly, because the count is 0 below the spectrum and N above it. A NaN pivot is never counted, so with a NaN in the representation the count stays short of the index and the loop runs forever: DSTEMR, and through it xSTEVR, xSYEVR and xHEEVR with RANGE = 'A', never return for a symmetric tridiagonal matrix with a NaN in one of its first two rows, with JOBZ = 'V' through xLARRB and with JOBZ = 'N' through xLARRJ. This PR stops the widening in both routines once BACK has overflowed, at which point no further step can move the endpoint, and returns INFO = 1, which the xLARRB callers already propagate (xLARRV returns -1, xLARRE -4, xSTEMR 21 or 11, and the drivers fall back to xSTEBZ + xSTEIN). It also checks xLARRB's return code at the one call site in xLARRV that ignored it, and xLARRJ's return code in xSTEMR, which never looked at it and now returns INFO = 3X for a failure there. Twelve files: {s,d}larrb.f, {s,d}larrj.f, {s,d,c,z}larrv.f, {s,d,c,z}stemr.f.

Description

The loops are

         BACK = WERR( II )
 20      CONTINUE
         NEGCNT = DLANEG( N, D, LLD, LEFT, PIVMIN, R )
         IF( NEGCNT.GT.I-1 ) THEN
            LEFT = LEFT - BACK
            BACK = TWO*BACK
            GO TO 20
         END IF

and its mirror image for RIGHT. DLANEG deliberately does not count a NaN pivot (DPLUS.LT.ZERO is false) and replaces a NaN quotient by 1 to continue the recurrence, so for a representation with a NaN in D(1) the count never exceeds N - 1 for any shift, and the RIGHT loop for the last eigenvalue widens through +Inf without end. The widening also stalls when WERR is exactly zero, which xLARRE produces for an eigenvalue that the dqds algorithm returns as exactly zero (WERR(I) = RTOL*ABS(W(I))): then BACK = 0, the endpoint never moves, and the loop terminates only if the Sturm count at W already agrees with the index.

xLARRJ has the same two loops, written with a factor FAC that doubles and a step WERR( II )*FAC, and the same Sturm count (DPLUS.LT.ZERO), so DSTEMR( 'N', 'A', ... ) on the same matrices never returns either: it skips xLARRV and refines the xLARRE eigenvalues with xLARRJ instead. xSTEMR ignores the INFO of xLARRJ.

Fix. In both routines BACK starts at MAX( WERR( II ), MNWDTH ), where MNWDTH = 2*PIVMIN is the minimum interval width the bisection already uses, so a zero WERR still makes progress; and each loop gives up with INFO = 1 when BACK.LT.TWO*BACK is false, i.e. when BACK is infinite (or NaN), because no further step can move the endpoint. The test is false only after BACK has already overflowed, which a finite representation never reaches: the count is 0 once LEFT is below the Gershgorin bound and N once RIGHT is above it, and BACK doubles from at least 2*PIVMIN to beyond the spectral diameter in at most MAXITR steps. INFO is documented in both (= 1: the interval could not be widened to contain the eigenvalue; the counts are inconsistent, as they are for a NaN or an Inf). xSTEMR now checks xLARRJ's INFO and returns INFO = 30 + ABS( IINFO ), documented next to the existing 1X (xLARRE) and 2X (xLARRV) codes. The second xLARRB call in xLARRV (refining the extremal eigenvalues of a child cluster) now returns INFO = -1 on failure like the first one instead of continuing with an unrefined bracket, in all four precisions: CLARRV and ZLARRV call the real xLARRB at the same three places as SLARRV and DLARRV, and without the check CSTEMR and ZSTEMR on the NaN matrices below go on into the representation construction and report the xLARRF failure (INFO = 22) where SSTEMR and DSTEMR report the refinement failure (INFO = 21).

Minimal reproducer

! DSYEVR('V', 'A') on a 3x3 symmetric tridiagonal matrix with A(1,1) = NaN.
program minimal
  implicit none
  integer, parameter :: n = 3
  double precision :: a(n,n), w(n), z(n,n), work(26*n), zero, d(n), e(n), dd(n), ee(n)
  integer :: iwork(10*n), isuppz(2*n), info, m, i
  logical :: tryrac
  zero = 0
  a = 0
  do i = 1, n
    a(i,i) = i
  end do
  do i = 1, n-1
    a(i,i+1) = 0.3d0; a(i+1,i) = 0.3d0
  end do
  a(1,1) = zero/zero
  d = [a(1,1), a(2,2), a(3,3)]
  e = [0.3d0, 0.3d0, 0d0]
  dd = d; ee = e; tryrac = .true.
  call dstemr('V', 'A', n, dd, ee, 0d0, 0d0, 1, n, m, w, z, n, n, isuppz, tryrac, work, 26*n, iwork, 10*n, info)
  print '(a,i3,a,i2)', 'DSTEMR(V): info = ', info, '  m = ', m
  dd = d; ee = e; tryrac = .true.
  call dstemr('N', 'A', n, dd, ee, 0d0, 0d0, 1, n, m, w, z, n, n, isuppz, tryrac, work, 26*n, iwork, 10*n, info)
  print '(a,i3,a,i2)', 'DSTEMR(N): info = ', info, '  m = ', m
  call dsyevr('V', 'A', 'U', n, a, n, 0d0, 0d0, 1, n, 0d0, m, w, z, n, isuppz, work, 26*n, iwork, 10*n, info)
  print '(a,i3,a,i2)', 'DSYEVR: info = ', info, '  m = ', m
end program

BEFORE (master):     (no output; the first DSTEMR call never returns)
AFTER (this branch): DSTEMR(V): info =  21  m =  3
                     DSTEMR(N): info =  31  m =  3
                     DSYEVR: info =   0  m =  0

DSYEVR returns through its DSTEBZ + DSTEIN fallback, which reports no eigenvalue for the NaN matrix, the same result it gives on master for a NaN that is not in the first two rows.

Regression test. xCHKST gets the case: after the size and type loops it calls xSTEMR with RANGE = 'A' on a tridiagonal matrix whose first diagonal entry is a NaN, once with JOBZ = 'V' and once with JOBZ = 'N', which covers xLARRB and xLARRJ respectively. Only a negative INFO counts as a failure, 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 commit all four xeigtst binaries stop responding at that call, so LAPACK-xeigtst{s,d,c,z}_sep_in fail by timeout; with the fix sep.in passes unchanged (13464 tests for the real precisions, 11016 for the complex).

Validation

The full LAPACK test suite passes on this branch: 5441941 LAPACK tests, 0 numerical errors, 0 other errors (ctest: 100% of 215 tests passed). The 40 extra tests over the parent commit are the new NaN calls, 5 sizes times 2 job options times 4 precisions.
100 NaN cases (DSTEMR, SSTEMR, DSYEVR, DSTEVR, ZHEEVR; NaN in D(1), D(2), E(1), D(N), E(N-1); N = 3, 5, 8, 20), one process per case with a 10 s timeout: master hangs in 60 of them (every case with the NaN in the first two rows); this branch returns in all 100, xSTEMR with INFO = 21 (48 cases) or INFO = 11 (32 cases, when xLARRE's own refinement hits the loop first), the drivers with INFO = 0 and no eigenvalues.
xSTEMR called directly in all four precisions on the 12 NaN positions in rows 1-2 (D(1), D(2), E(1); N = 3, 5, 8, 20) returns INFO = 21 in every case; CSTEMR and ZSTEMR returned INFO = 22 before the {c,z}larrv.f check.
xSTEMR( 'N', 'A', ... ) in all four precisions on five NaN positions (D(1), D(2), E(1), D(N), E(N-1); N = 3, 5, 8, 20), one process per case with a 10 s timeout: master hangs in 48 of the 80 cases; this branch returns in all 80, with INFO = 31 in every case that hung (repro/probe8n.f90, run8n.sh).
200 finite matrices (N = 1 ... 40; random, glued-Wilkinson clusters, scaled to 1e-150 and to 1e100, nearly diagonal with off-diagonals 1e-9) give bit-identical eigenvalues from DSTEMR on master and on this branch, with JOBZ = 'V' and with JOBZ = 'N', so the MNWDTH floor never changed a computed bracket in that sample.

Found while auditing the tridiagonal bisection routines (the same NaN sweep that produced #1384); xLARRB and xLARRJ are the only routines in the symmetric eigenvalue path that hang rather than returning garbage for a NaN. The xLARRJ half was found afterwards, when preparing a NaN test for the suite that calls DSTEMR with JOBZ = 'N' as dchkst.f does.

codecov · 2026-09-07T23:37:57Z

Codecov Report

❌ Patch coverage is 87.20930% with 22 lines in your changes missing coverage. Please review.
✅ Project coverage is 69.38%. Comparing base (a6c6e74) to head (fcc5958).
✅ All tests successful. No failed tests found.

Files with missing lines	Patch %	Lines
SRC/dlarrj.f	66.66%	4 Missing ⚠️
SRC/slarrj.f	66.66%	4 Missing ⚠️
SRC/dlarrb.f	75.00%	3 Missing ⚠️
SRC/slarrb.f	75.00%	3 Missing ⚠️
TESTING/EIG/cchkst.f	92.00%	2 Missing ⚠️
TESTING/EIG/dchkst.f	92.00%	2 Missing ⚠️
TESTING/EIG/schkst.f	92.00%	2 Missing ⚠️
TESTING/EIG/zchkst.f	92.00%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1389      +/-   ##
==========================================
+ Coverage   69.36%   69.38%   +0.02%     
==========================================
  Files        6122     6122              
  Lines      486711   486847     +136     
  Branches    23268    23268              
==========================================
+ Hits       337584   337815     +231     
+ Misses     148689   148594      -95     
  Partials      438      438

Components	Coverage Δ
BLAS	`97.94% <ø> (ø)`
CBLAS	`96.98% <ø> (ø)`
LAPACK	`82.44% <80.55%> (+0.05%)`	⬆️
LAPACKE	`2.17% <ø> (ø)`
TMGLIB	`55.69% <ø> (ø)`
BLAS testing	`88.33% <ø> (ø)`
CBLAS testing	`89.63% <ø> (ø)`
LAPACK testing	`82.21% <92.00%> (+<0.01%)`	⬆️
LAPACKE testing	`∅ <ø> (∅)`

Files with missing lines	Coverage Δ
SRC/clarrv.f	`68.77% <100.00%> (+3.21%)`	⬆️
SRC/cstemr.f	`86.49% <100.00%> (+1.02%)`	⬆️
SRC/dlarrv.f	`67.69% <100.00%> (+3.33%)`	⬆️
SRC/dstemr.f	`86.49% <100.00%> (+1.02%)`	⬆️
SRC/slarrv.f	`70.44% <100.00%> (+3.31%)`	⬆️
SRC/sstemr.f	`86.49% <100.00%> (+1.02%)`	⬆️
SRC/zlarrv.f	`66.11% <100.00%> (+3.23%)`	⬆️
SRC/zstemr.f	`86.49% <100.00%> (+1.02%)`	⬆️
TESTING/EIG/cchkst.f	`64.38% <92.00%> (+1.23%)`	⬆️
TESTING/EIG/dchkst.f	`64.40% <92.00%> (+1.25%)`	⬆️
... and 6 more

... and 2 files with indirect coverage changes

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...fcc5958. Read the comment docs.

rmlarsen · 2026-09-10T18:54:50Z

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.

…cannot hang xSTEMR Before bisecting, xLARRB widens the initial interval [W-WERR, W+WERR] until the Sturm count of the representation agrees with the eigenvalue index, moving the endpoint by BACK and doubling BACK at every step. In exact arithmetic that ends quickly, because the count is 0 below the spectrum and N above it. A NaN pivot is never counted, so with a NaN in the representation the count can stay short of the index and the loop runs forever: xSTEMR, and through it xSTEVR, xSYEVR and xHEEVR with RANGE = 'A', never return for a symmetric tridiagonal matrix with a NaN in one of its first two rows. The loop also stalls when WERR is exactly zero, which xLARRE produces for an eigenvalue that dqds returns as exactly zero. xLARRJ, which refines the eigenvalues by bisection on the matrix itself when only eigenvalues are wanted, has the same two widening loops with the same step and the same hang: xSTEMR with JOBZ = 'N' never returns for the same matrices. Start BACK at no less than the minimum interval width 2*PIVMIN, and stop widening with INFO = 1 once BACK has overflowed, when no further step can move the endpoint; a finite matrix never gets there, since BACK passes the spectral diameter within MAXITR doublings. The xLARRB callers already propagate a nonzero INFO (xLARRV returns -1, xLARRE -4, xSTEMR 21 or 11) and the drivers then fall back to xSTEBZ and xSTEIN; the one xLARRB call in xLARRV that ignored INFO now checks it like the others, in all four precisions. xSTEMR never looked at the INFO of xLARRJ; it now returns INFO = 3X for a failure there, next to the documented 1X and 2X codes of xLARRE and xLARRV. xCHKST gets the case as a regression test: after the size and type loops it calls xSTEMR with RANGE = 'A' on a tridiagonal matrix whose first diagonal entry is a NaN, once with JOBZ = 'V' and once with JOBZ = 'N', so that both refinement paths are covered. Only a negative INFO is reported as a failure; the test is that the call returns at all. On the parent all four xeigtst binaries stop responding there. Over 100 NaN cases (DSTEMR, SSTEMR, DSYEVR, DSTEVR, ZHEEVR; five NaN positions; n = 3, 5, 8, 20) the parent hangs in 60, this branch returns in all of them; with JOBZ = 'N' (xSTEMR in all four precisions, five NaN positions, n = 3, 5, 8, 20) the parent hangs in 48 of 80, this branch returns INFO = 31 in every case that hung. 200 finite matrices give bit-identical eigenvalues from DSTEMR on both, with JOBZ = 'V' and with JOBZ = 'N'. The full LAPACK test suite passes: 5441941 LAPACK tests, 0 numerical errors, 0 other errors, 40 tests more than the parent from the new calls. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

xSTEMR scales the matrix into the range xLARRE can represent: ELSE IF( TNRM.GT.RMAX ) THEN SCALE = RMAX / TNRM For a matrix with an infinite entry TNRM is infinite, SCALE underflows to zero, the scaling turns the infinite entry into a NaN and every other entry into zero, and the solver ran on that matrix: RANGE = 'A' and 'I' returned INFO = 0 with every eigenvalue a NaN, RANGE = 'V' scaled the interval to (0, 0] and returned no eigenvalues, and with the infinite entry off the diagonal the representation tree search never returned (the NaN hang of Reference-LAPACK#1389). Leave a matrix whose scale factor would be zero unscaled, and let xLARRE, which computes the base representations, reject a matrix with an infinite entry as one for which no representation exists, INFO = 2, the code it already uses when the search for a representation gives up. xSTEMR returns INFO = 12 and the drivers that call it, xSTEVR, xSYEVR and xHEEVR, fall back to xSTEBZ and xSTEIN as for any other xSTEMR failure. A NaN matrix is not touched by this change: its norm compares false against every threshold and it takes the same path as before. xCHKST gets the case as a regression test: after the size and type loops it calls xSTEMR with RANGE = 'A' on an N = NMAX tridiagonal matrix whose first diagonal entry is +Inf, once with JOBZ = 'V' and once with JOBZ = 'N', and reports INFO = 0 as a failure. On the parent commit all four calls return INFO = 0 with NaN eigenvalues; here they return INFO = 12. Over 420 placements of +Inf and -Inf (n = 3 to 64, RANGE = 'A', 'V', 'I', real and complex) the parent returned INFO = 0 with NaN eigenvalues in 52 cases, hung in 24 and returned INFO = 11 in the rest; this branch returns INFO = 12 in all of them. The 1 by 1 and 2 by 2 matrices, which xSTEMR solves in closed form before scaling, are unchanged, and every finite case is bit-identical. The full LAPACK test suite passes: 0 numerical errors, 0 other errors, 40 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>

rmlarsen mentioned this pull request Sep 7, 2026

Fix rescaling after a partial failure in the xSTEBZ-based drivers; accept Inf input and xSTEBZ's negative block numbers #1384

Open

rmlarsen force-pushed the larrb-bounded-bracket branch from a0d34af to dded2ab Compare September 8, 2026 03:11

rmlarsen changed the title ~~Bound the bracket widening in xLARRB so that a NaN representation cannot hang it~~ Bound the bracket widening in xLARRB and xLARRJ so that a NaN matrix cannot hang xSTEMR Sep 8, 2026

rmlarsen force-pushed the larrb-bounded-bracket branch from dded2ab to 111c8c5 Compare September 8, 2026 03:47

rmlarsen mentioned this pull request Sep 10, 2026

Do not scale an infinite matrix to zero in xSTEMR; report it from xLARRE #1404

Open

rmlarsen force-pushed the larrb-bounded-bracket branch from 111c8c5 to fcc5958 Compare September 15, 2026 20:28

Navigation Menu

Uh oh!

Bound the bracket widening in xLARRB and xLARRJ so that a NaN matrix cannot hang xSTEMR - #1389

Bound the bracket widening in xLARRB and xLARRJ so that a NaN matrix cannot hang xSTEMR#1389
rmlarsen wants to merge 1 commit into
Reference-LAPACK:masterfrom
rmlarsen:larrb-bounded-bracket

rmlarsen commented Sep 7, 2026 •

edited

Loading

Uh oh!

codecov Bot commented Sep 7, 2026 •

edited

Loading

Uh oh!

rmlarsen commented Sep 10, 2026 •

edited

Loading

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

Uh oh!

Conversation

rmlarsen commented Sep 7, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

codecov Bot commented Sep 7, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

rmlarsen commented Sep 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

rmlarsen commented Sep 7, 2026 •

edited

Loading

codecov Bot commented Sep 7, 2026 •

edited

Loading

rmlarsen commented Sep 10, 2026 •

edited

Loading