Bound the bracket widening in xLARRB and xLARRJ so that a NaN matrix 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>

Author: rmlarsen

SRC/clarrv.f

@@ -702,6 +702,10 @@ SUBROUTINE CLARRV( N, VL, VU, D, L, PIVMIN,
      $                       WERR(WBEGIN),WORK( INDWRK ),
      $                       IWORK( IINDWK ), PIVMIN, SPDIAM,
      $                       IN, IINFO )
+                        IF( IINFO.NE.0 ) THEN
+                           INFO = -1
+                           RETURN
+                        END IF
  55                  CONTINUE
 *
                      IF((WBEGIN+NEWLST-1.LT.DOL).OR.

SRC/cstemr.f

@@ -304,10 +304,11 @@
 *>          = 0:  successful exit
 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 *>          > 0:  if INFO = 1X, internal error in SLARRE,
-*>                if INFO = 2X, internal error in CLARRV.
+*>                if INFO = 2X, internal error in CLARRV,
+*>                if INFO = 3X, internal error in SLARRJ.
 *>                Here, the digit X = ABS( IINFO ) < 10, where IINFO is
-*>                the nonzero error code returned by SLARRE or
-*>                CLARRV, respectively.
+*>                the nonzero error code returned by SLARRE, CLARRV
+*>                or SLARRJ, respectively.
 *> \endverbatim
 *
 *  Authors:
@@ -752,6 +753,10 @@ SUBROUTINE CSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
      $                   WORK( INDERR+WBEGIN-1 ),
      $                   WORK( INDWRK ), IWORK( IINDWK ), PIVMIN,
      $                   TNRM, IINFO )
+               IF( IINFO.NE.0 ) THEN
+                  INFO = 30 + ABS( IINFO )
+                  RETURN
+               END IF
                IBEGIN = IEND + 1
                WBEGIN = WEND + 1
  39      CONTINUE

SRC/dlarrb.f

@@ -165,7 +165,11 @@
 *> \param[out] INFO
 *> \verbatim
 *>          INFO is INTEGER
-*>          Error flag.
+*>          = 0:  successful exit
+*>          = 1:  the interval around one of the eigenvalues could not
+*>                be widened to contain it: the Sturm counts of the
+*>                representation are inconsistent, as they are when it
+*>                contains a NaN or an Inf.
 *> \endverbatim
 *
 *  Authors:
@@ -274,10 +278,23 @@ SUBROUTINE DLARRB( N, D, LLD, IFIRST, ILAST, RTOL1,
 *
 *        Do while( NEGCNT(LEFT).GT.I-1 )
 *
-         BACK = WERR( II )
+*        The interval is widened by BACK, doubled at every step, until
+*        the Sturm count agrees with the index.  In exact arithmetic
+*        this terminates because the count is 0 below the spectrum and
+*        N above it.  A NaN or an Inf in the representation is never
+*        counted, so the count can stay short of the index; stop once
+*        BACK has overflowed, when no further step can move the
+*        endpoint.  BACK starts at no less than the minimum interval
+*        width so that a zero WERR still makes progress.
+*
+         BACK = MAX( WERR( II ), MNWDTH )
  20      CONTINUE
          NEGCNT = DLANEG( N, D, LLD, LEFT, PIVMIN, R )
          IF( NEGCNT.GT.I-1 ) THEN
+            IF( .NOT.( BACK.LT.TWO*BACK ) ) THEN
+               INFO = 1
+               RETURN
+            END IF
             LEFT = LEFT - BACK
             BACK = TWO*BACK
             GO TO 20
@@ -286,15 +303,18 @@ SUBROUTINE DLARRB( N, D, LLD, IFIRST, ILAST, RTOL1,
 *        Do while( NEGCNT(RIGHT).LT.I )
 *        Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT
 *
-         BACK = WERR( II )
+         BACK = MAX( WERR( II ), MNWDTH )
  50      CONTINUE
-
          NEGCNT = DLANEG( N, D, LLD, RIGHT, PIVMIN, R )
-          IF( NEGCNT.LT.I ) THEN
-             RIGHT = RIGHT + BACK
-             BACK = TWO*BACK
-             GO TO 50
-          END IF
+         IF( NEGCNT.LT.I ) THEN
+            IF( .NOT.( BACK.LT.TWO*BACK ) ) THEN
+               INFO = 1
+               RETURN
+            END IF
+            RIGHT = RIGHT + BACK
+            BACK = TWO*BACK
+            GO TO 50
+         END IF
          WIDTH = HALF*ABS( LEFT - RIGHT )
          TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
          CVRGD = MAX(RTOL1*GAP,RTOL2*TMP)

SRC/dlarrj.f

@@ -137,7 +137,11 @@
 *> \param[out] INFO
 *> \verbatim
 *>          INFO is INTEGER
-*>          Error flag.
+*>          = 0:  successful exit
+*>          = 1:  the interval around one of the eigenvalues could not
+*>                be widened to contain it: the Sturm counts of the
+*>                matrix are inconsistent, as they are when it contains
+*>                a NaN or an Inf.
 *> \endverbatim
 *
 *  Authors:
@@ -190,7 +194,7 @@ SUBROUTINE DLARRJ( N, D, E2, IFIRST, ILAST,
 *     .. Local Scalars ..
       INTEGER            CNT, I, I1, I2, II, ITER, J, K, NEXT, NINT,
      $                   OLNINT, P, PREV, SAVI1
-      DOUBLE PRECISION   DPLUS, FAC, LEFT, MID, RIGHT, S, TMP, WIDTH
+      DOUBLE PRECISION   BACK, DPLUS, LEFT, MID, RIGHT, S, TMP, WIDTH
 *
 *     ..
 *     .. Intrinsic Functions ..
@@ -249,7 +253,16 @@ SUBROUTINE DLARRJ( N, D, E2, IFIRST, ILAST,
 *
 *           Do while( CNT(LEFT).GT.I-1 )
 *
-            FAC = ONE
+*           The interval is widened by BACK, doubled at every step,
+*           until the Sturm count agrees with the index.  In exact
+*           arithmetic this terminates because the count is 0 below the
+*           spectrum and N above it.  A NaN or an Inf in the matrix is
+*           never counted, so the count can stay short of the index;
+*           stop once BACK has overflowed, when no further step can move
+*           the endpoint.  BACK starts at no less than the minimum
+*           interval width so that a zero WERR still makes progress.
+*
+            BACK = MAX( WERR( II ), TWO*PIVMIN )
  20         CONTINUE
             CNT = 0
             S = LEFT
@@ -260,14 +273,18 @@ SUBROUTINE DLARRJ( N, D, E2, IFIRST, ILAST,
                IF( DPLUS.LT.ZERO ) CNT = CNT + 1
  30         CONTINUE
             IF( CNT.GT.I-1 ) THEN
-               LEFT = LEFT - WERR( II )*FAC
-               FAC = TWO*FAC
+               IF( .NOT.( BACK.LT.TWO*BACK ) ) THEN
+                  INFO = 1
+                  RETURN
+               END IF
+               LEFT = LEFT - BACK
+               BACK = TWO*BACK
                GO TO 20
             END IF
 *
 *           Do while( CNT(RIGHT).LT.I )
 *
-            FAC = ONE
+            BACK = MAX( WERR( II ), TWO*PIVMIN )
  50         CONTINUE
             CNT = 0
             S = RIGHT
@@ -278,8 +295,12 @@ SUBROUTINE DLARRJ( N, D, E2, IFIRST, ILAST,
                IF( DPLUS.LT.ZERO ) CNT = CNT + 1
  60         CONTINUE
             IF( CNT.LT.I ) THEN
-               RIGHT = RIGHT + WERR( II )*FAC
-               FAC = TWO*FAC
+               IF( .NOT.( BACK.LT.TWO*BACK ) ) THEN
+                  INFO = 1
+                  RETURN
+               END IF
+               RIGHT = RIGHT + BACK
+               BACK = TWO*BACK
                GO TO 50
             END IF
             NINT = NINT + 1

SRC/dlarrv.f

@@ -698,6 +698,10 @@ SUBROUTINE DLARRV( N, VL, VU, D, L, PIVMIN,
      $                       WERR(WBEGIN),WORK( INDWRK ),
      $                       IWORK( IINDWK ), PIVMIN, SPDIAM,
      $                       IN, IINFO )
+                        IF( IINFO.NE.0 ) THEN
+                           INFO = -1
+                           RETURN
+                        END IF
  55                  CONTINUE
 *
                      IF((WBEGIN+NEWLST-1.LT.DOL).OR.

SRC/dstemr.f

@@ -287,10 +287,11 @@
 *>          = 0:  successful exit
 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 *>          > 0:  if INFO = 1X, internal error in DLARRE,
-*>                if INFO = 2X, internal error in DLARRV.
+*>                if INFO = 2X, internal error in DLARRV,
+*>                if INFO = 3X, internal error in DLARRJ.
 *>                Here, the digit X = ABS( IINFO ) < 10, where IINFO is
-*>                the nonzero error code returned by DLARRE or
-*>                DLARRV, respectively.
+*>                the nonzero error code returned by DLARRE, DLARRV
+*>                or DLARRJ, respectively.
 *> \endverbatim
 *
 *  Authors:
@@ -735,6 +736,10 @@ SUBROUTINE DSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
      $                   WORK( INDERR+WBEGIN-1 ),
      $                   WORK( INDWRK ), IWORK( IINDWK ), PIVMIN,
      $                   TNRM, IINFO )
+               IF( IINFO.NE.0 ) THEN
+                  INFO = 30 + ABS( IINFO )
+                  RETURN
+               END IF
                IBEGIN = IEND + 1
                WBEGIN = WEND + 1
  39         CONTINUE

SRC/slarrb.f

@@ -165,7 +165,11 @@
 *> \param[out] INFO
 *> \verbatim
 *>          INFO is INTEGER
-*>          Error flag.
+*>          = 0:  successful exit
+*>          = 1:  the interval around one of the eigenvalues could not
+*>                be widened to contain it: the Sturm counts of the
+*>                representation are inconsistent, as they are when it
+*>                contains a NaN or an Inf.
 *> \endverbatim
 *
 *  Authors:
@@ -274,10 +278,23 @@ SUBROUTINE SLARRB( N, D, LLD, IFIRST, ILAST, RTOL1,
 *
 *        Do while( NEGCNT(LEFT).GT.I-1 )
 *
-         BACK = WERR( II )
+*        The interval is widened by BACK, doubled at every step, until
+*        the Sturm count agrees with the index.  In exact arithmetic
+*        this terminates because the count is 0 below the spectrum and
+*        N above it.  A NaN or an Inf in the representation is never
+*        counted, so the count can stay short of the index; stop once
+*        BACK has overflowed, when no further step can move the
+*        endpoint.  BACK starts at no less than the minimum interval
+*        width so that a zero WERR still makes progress.
+*
+         BACK = MAX( WERR( II ), MNWDTH )
  20      CONTINUE
          NEGCNT = SLANEG( N, D, LLD, LEFT, PIVMIN, R )
          IF( NEGCNT.GT.I-1 ) THEN
+            IF( .NOT.( BACK.LT.TWO*BACK ) ) THEN
+               INFO = 1
+               RETURN
+            END IF
             LEFT = LEFT - BACK
             BACK = TWO*BACK
             GO TO 20
@@ -286,15 +303,18 @@ SUBROUTINE SLARRB( N, D, LLD, IFIRST, ILAST, RTOL1,
 *        Do while( NEGCNT(RIGHT).LT.I )
 *        Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT
 *
-         BACK = WERR( II )
+         BACK = MAX( WERR( II ), MNWDTH )
  50      CONTINUE
-
          NEGCNT = SLANEG( N, D, LLD, RIGHT, PIVMIN, R )
-          IF( NEGCNT.LT.I ) THEN
-             RIGHT = RIGHT + BACK
-             BACK = TWO*BACK
-             GO TO 50
-          END IF
+         IF( NEGCNT.LT.I ) THEN
+            IF( .NOT.( BACK.LT.TWO*BACK ) ) THEN
+               INFO = 1
+               RETURN
+            END IF
+            RIGHT = RIGHT + BACK
+            BACK = TWO*BACK
+            GO TO 50
+         END IF
          WIDTH = HALF*ABS( LEFT - RIGHT )
          TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
          CVRGD = MAX(RTOL1*GAP,RTOL2*TMP)

SRC/slarrj.f

@@ -137,7 +137,11 @@
 *> \param[out] INFO
 *> \verbatim
 *>          INFO is INTEGER
-*>          Error flag.
+*>          = 0:  successful exit
+*>          = 1:  the interval around one of the eigenvalues could not
+*>                be widened to contain it: the Sturm counts of the
+*>                matrix are inconsistent, as they are when it contains
+*>                a NaN or an Inf.
 *> \endverbatim
 *
 *  Authors:
@@ -190,7 +194,7 @@ SUBROUTINE SLARRJ( N, D, E2, IFIRST, ILAST,
 *     .. Local Scalars ..
       INTEGER            CNT, I, I1, I2, II, ITER, J, K, NEXT, NINT,
      $                   OLNINT, P, PREV, SAVI1
-      REAL               DPLUS, FAC, LEFT, MID, RIGHT, S, TMP, WIDTH
+      REAL               BACK, DPLUS, LEFT, MID, RIGHT, S, TMP, WIDTH
 *
 *     ..
 *     .. Intrinsic Functions ..
@@ -249,7 +253,16 @@ SUBROUTINE SLARRJ( N, D, E2, IFIRST, ILAST,
 *
 *           Do while( CNT(LEFT).GT.I-1 )
 *
-            FAC = ONE
+*           The interval is widened by BACK, doubled at every step,
+*           until the Sturm count agrees with the index.  In exact
+*           arithmetic this terminates because the count is 0 below the
+*           spectrum and N above it.  A NaN or an Inf in the matrix is
+*           never counted, so the count can stay short of the index;
+*           stop once BACK has overflowed, when no further step can move
+*           the endpoint.  BACK starts at no less than the minimum
+*           interval width so that a zero WERR still makes progress.
+*
+            BACK = MAX( WERR( II ), TWO*PIVMIN )
  20         CONTINUE
             CNT = 0
             S = LEFT
@@ -260,14 +273,18 @@ SUBROUTINE SLARRJ( N, D, E2, IFIRST, ILAST,
                IF( DPLUS.LT.ZERO ) CNT = CNT + 1
  30         CONTINUE
             IF( CNT.GT.I-1 ) THEN
-               LEFT = LEFT - WERR( II )*FAC
-               FAC = TWO*FAC
+               IF( .NOT.( BACK.LT.TWO*BACK ) ) THEN
+                  INFO = 1
+                  RETURN
+               END IF
+               LEFT = LEFT - BACK
+               BACK = TWO*BACK
                GO TO 20
             END IF
 *
 *           Do while( CNT(RIGHT).LT.I )
 *
-            FAC = ONE
+            BACK = MAX( WERR( II ), TWO*PIVMIN )
  50         CONTINUE
             CNT = 0
             S = RIGHT
@@ -278,8 +295,12 @@ SUBROUTINE SLARRJ( N, D, E2, IFIRST, ILAST,
                IF( DPLUS.LT.ZERO ) CNT = CNT + 1
  60         CONTINUE
             IF( CNT.LT.I ) THEN
-               RIGHT = RIGHT + WERR( II )*FAC
-               FAC = TWO*FAC
+               IF( .NOT.( BACK.LT.TWO*BACK ) ) THEN
+                  INFO = 1
+                  RETURN
+               END IF
+               RIGHT = RIGHT + BACK
+               BACK = TWO*BACK
                GO TO 50
             END IF
             NINT = NINT + 1

SRC/slarrv.f

@@ -698,6 +698,10 @@ SUBROUTINE SLARRV( N, VL, VU, D, L, PIVMIN,
      $                       WERR(WBEGIN),WORK( INDWRK ),
      $                       IWORK( IINDWK ), PIVMIN, SPDIAM,
      $                       IN, IINFO )
+                        IF( IINFO.NE.0 ) THEN
+                           INFO = -1
+                           RETURN
+                        END IF
  55                  CONTINUE
 *
                      IF((WBEGIN+NEWLST-1.LT.DOL).OR.

SRC/sstemr.f

@@ -287,10 +287,11 @@
 *>          = 0:  successful exit
 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 *>          > 0:  if INFO = 1X, internal error in SLARRE,
-*>                if INFO = 2X, internal error in SLARRV.
+*>                if INFO = 2X, internal error in SLARRV,
+*>                if INFO = 3X, internal error in SLARRJ.
 *>                Here, the digit X = ABS( IINFO ) < 10, where IINFO is
-*>                the nonzero error code returned by SLARRE or
-*>                SLARRV, respectively.
+*>                the nonzero error code returned by SLARRE, SLARRV
+*>                or SLARRJ, respectively.
 *> \endverbatim
 *
 *  Authors:
@@ -732,6 +733,10 @@ SUBROUTINE SSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
      $                   WORK( INDERR+WBEGIN-1 ),
      $                   WORK( INDWRK ), IWORK( IINDWK ), PIVMIN,
      $                   TNRM, IINFO )
+               IF( IINFO.NE.0 ) THEN
+                  INFO = 30 + ABS( IINFO )
+                  RETURN
+               END IF
                IBEGIN = IEND + 1
                WBEGIN = WEND + 1
  39      CONTINUE