From 63ed409c31f4adbc7bcfc64c6cf5a02f65f04636 Mon Sep 17 00:00:00 2001
From: Rasmus Munk Larsen <rmlarsen@gmail.com>
Date: Mon, 7 Sep 2026 21:27:23 -0700
Subject: [PATCH] Rescale every eigenvalue after a partial failure in the
 xSTEBZ-based drivers

The 32 symmetric eigensolver drivers that reduce to xSTEBZ + xSTEIN
(xSTEVX, xSYEVX, xHEEVX, xSPEVX, xHPEVX, xSBEVX, xHBEVX, xSTEVR,
xSYEVR, xHEEVR and the _2STAGE variants) scale a matrix whose largest
entry lies outside [RMIN, RMAX] and undo the scaling of W at the end
with IMAX = M if INFO = 0 and IMAX = INFO - 1 otherwise.  That block
was copied from xSTEV and xSYEV, where a nonzero INFO comes from
xSTEQR.  Here a nonzero INFO means that INFO eigenvectors failed to
converge in xSTEIN, in which case the eigenvalues are all valid and
W(INFO:M) are returned still multiplied by SIGMA, or that xSTEBZ
returned INFO - N, in which case IMAX = N + i - 1 and xSCAL writes past
the end of W.  A failure of xSTERF, xSTEQR or xSTEMR earlier in the
driver resets INFO and falls through to xSTEBZ, so at the rescale label
M is the right count in every case.

Two related defects sit in the same drivers.  A matrix containing an
Inf has ANRM = Inf and SIGMA = RMAX / ANRM = 0, so for RANGE = 'V' the
scaled interval collapses to VLL = VUU = 0 and xSTEBZ stops the process
in XERBLA; compute SIGMA first and scale only if it is positive, which
sends an infinite matrix down the same path as one containing a NaN.
And xSTEIN rejects with INFO = -6, again XERBLA, or silently skips the
negative block numbers with which xSTEBZ flags eigenvalues that did not
converge, although every driver deliberately continues into xSTEIN
when xSTEBZ returns INFO = 1; take ABS( IBLOCK( . ) ) at its four uses
and document it.

DSTEVX on 1e-150 tridiag(1, 0.5, 1) with ABSTOL = 0.3 |T| reports
INFO = 6 and returns the sixth eigenvalue 1e4 times too large on the
parent, all six correctly on this branch.  xERRST gets the infinite matrix as a regression test: xSTEVX and
xSYEVX with RANGE = 'V' on a tridiagonal whose first diagonal entry is
an infinity, which must return without reaching XERBLA.  On the parent
both report the illegal fifth argument of xSTEBZ.  The other two
changes are not reachable from the test suite: a partial xSTEIN failure
needs a matrix outside the range the drivers document, and the negative
block numbers likewise.

Over 1584 NaN/Inf cases
(8 routines, 3 special values, 5 positions, 5 orders, 3 ranges, one
process each with XERBLA overridden and canaries around every array)
the parent aborts 280 times and writes past the end of W 200 times;
this branch never does either.  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>
---
 SRC/chbevx.f         | 17 +++++++---------
 SRC/chbevx_2stage.f  | 17 +++++++---------
 SRC/cheevr.f         | 17 +++++++---------
 SRC/cheevr_2stage.f  | 17 +++++++---------
 SRC/cheevx.f         | 17 +++++++---------
 SRC/cheevx_2stage.f  | 17 +++++++---------
 SRC/chpevx.f         | 17 +++++++---------
 SRC/cstein.f         | 11 +++++++----
 SRC/dsbevx.f         | 17 +++++++---------
 SRC/dsbevx_2stage.f  | 17 +++++++---------
 SRC/dspevx.f         | 17 +++++++---------
 SRC/dstein.f         | 11 +++++++----
 SRC/dstevr.f         | 17 +++++++---------
 SRC/dstevx.f         | 17 +++++++---------
 SRC/dsyevr.f         | 17 +++++++---------
 SRC/dsyevr_2stage.f  | 17 +++++++---------
 SRC/dsyevx.f         | 17 +++++++---------
 SRC/dsyevx_2stage.f  | 17 +++++++---------
 SRC/ssbevx.f         | 17 +++++++---------
 SRC/ssbevx_2stage.f  | 17 +++++++---------
 SRC/sspevx.f         | 17 +++++++---------
 SRC/sstein.f         | 11 +++++++----
 SRC/sstevr.f         | 17 +++++++---------
 SRC/sstevx.f         | 17 +++++++---------
 SRC/ssyevr.f         | 17 +++++++---------
 SRC/ssyevr_2stage.f  | 17 +++++++---------
 SRC/ssyevx.f         | 17 +++++++---------
 SRC/ssyevx_2stage.f  | 17 +++++++---------
 SRC/zhbevx.f         | 17 +++++++---------
 SRC/zhbevx_2stage.f  | 17 +++++++---------
 SRC/zheevr.f         | 17 +++++++---------
 SRC/zheevr_2stage.f  | 17 +++++++---------
 SRC/zheevx.f         | 17 +++++++---------
 SRC/zheevx_2stage.f  | 17 +++++++---------
 SRC/zhpevx.f         | 17 +++++++---------
 SRC/zstein.f         | 11 +++++++----
 TESTING/EIG/derrst.f | 46 ++++++++++++++++++++++++++++++++++++++++++++
 TESTING/EIG/serrst.f | 46 ++++++++++++++++++++++++++++++++++++++++++++
 38 files changed, 344 insertions(+), 336 deletions(-)

diff --git a/SRC/chbevx.f b/SRC/chbevx.f
index 439d85a70b..8ff2aa2fe9 100644
--- a/SRC/chbevx.f
+++ b/SRC/chbevx.f
@@ -296,7 +296,7 @@ SUBROUTINE CHBEVX( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ,
 *     .. Local Scalars ..
       LOGICAL            ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ
       CHARACTER          ORDER
-      INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
+      INTEGER            I, IINFO, INDD, INDE, INDEE, INDIBL,
      $                   INDISP, INDIWK, INDRWK, INDWRK, ISCALE, ITMP1,
      $                   J, JJ, NSPLIT
       REAL               ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
@@ -415,8 +415,11 @@ SUBROUTINE CHBEVX( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -517,14 +520,8 @@ SUBROUTINE CHBEVX( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    30 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL SSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/chbevx_2stage.f b/SRC/chbevx_2stage.f
index c05adc4553..9398baf1fe 100644
--- a/SRC/chbevx_2stage.f
+++ b/SRC/chbevx_2stage.f
@@ -357,7 +357,7 @@ SUBROUTINE CHBEVX_2STAGE( JOBZ, RANGE, UPLO, N, KD, AB, LDAB,
       LOGICAL            ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ,
      $                   LQUERY
       CHARACTER          ORDER
-      INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
+      INTEGER            I, IINFO, INDD, INDE, INDEE, INDIBL,
      $                   INDISP, INDIWK, INDRWK, INDWRK, ISCALE, ITMP1,
      $                   LLWORK, LWMIN, LHTRD, LWTRD, IB, INDHOUS,
      $                   J, JJ, NSPLIT
@@ -501,8 +501,11 @@ SUBROUTINE CHBEVX_2STAGE( JOBZ, RANGE, UPLO, N, KD, AB, LDAB,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -608,14 +611,8 @@ SUBROUTINE CHBEVX_2STAGE( JOBZ, RANGE, UPLO, N, KD, AB, LDAB,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    30 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL SSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/cheevr.f b/SRC/cheevr.f
index bab9353f8a..c8d8179d41 100644
--- a/SRC/cheevr.f
+++ b/SRC/cheevr.f
@@ -398,7 +398,7 @@ SUBROUTINE CHEEVR( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
       LOGICAL            ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
      $                   WANTZ, TRYRAC
       CHARACTER          ORDER
-      INTEGER            I, IEEEOK, IINFO, IMAX, INDIBL, INDIFL, INDISP,
+      INTEGER            I, IEEEOK, IINFO, INDIBL, INDIFL, INDISP,
      $                   INDIWO, INDRD, INDRDD, INDRE, INDREE, INDRWK,
      $                   INDTAU, INDWK, INDWKN, ISCALE, ITMP1, J, JJ,
      $                   LIWMIN, LLWORK, LLRWORK, LLWRKN, LRWMIN,
@@ -549,8 +549,11 @@ SUBROUTINE CHEEVR( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -700,14 +703,8 @@ SUBROUTINE CHEEVR( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    30 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL SSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/cheevr_2stage.f b/SRC/cheevr_2stage.f
index 5bd16f449c..f2ac021c73 100644
--- a/SRC/cheevr_2stage.f
+++ b/SRC/cheevr_2stage.f
@@ -433,7 +433,7 @@ SUBROUTINE CHEEVR_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
       LOGICAL            ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
      $                   WANTZ, TRYRAC
       CHARACTER          ORDER
-      INTEGER            I, IEEEOK, IINFO, IMAX, INDIBL, INDIFL, INDISP,
+      INTEGER            I, IEEEOK, IINFO, INDIBL, INDIFL, INDISP,
      $                   INDIWO, INDRD, INDRDD, INDRE, INDREE, INDRWK,
      $                   INDTAU, INDWK, INDWKN, ISCALE, ITMP1, J, JJ,
      $                   LIWMIN, LLWORK, LLRWORK, LLWRKN, LRWMIN,
@@ -590,8 +590,11 @@ SUBROUTINE CHEEVR_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -744,14 +747,8 @@ SUBROUTINE CHEEVR_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    30 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL SSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/cheevx.f b/SRC/cheevx.f
index 54324830e0..15d24ce646 100644
--- a/SRC/cheevx.f
+++ b/SRC/cheevx.f
@@ -287,7 +287,7 @@ SUBROUTINE CHEEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
       LOGICAL            ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
      $                   WANTZ
       CHARACTER          ORDER
-      INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
+      INTEGER            I, IINFO, INDD, INDE, INDEE, INDIBL,
      $                   INDISP, INDIWK, INDRWK, INDTAU, INDWRK, ISCALE,
      $                   ITMP1, J, JJ, LLWORK, LWKMIN, LWKOPT, NB,
      $                   NSPLIT
@@ -419,8 +419,11 @@ SUBROUTINE CHEEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -525,14 +528,8 @@ SUBROUTINE CHEEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    40 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL SSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/cheevx_2stage.f b/SRC/cheevx_2stage.f
index 5b93c5aab8..cde12aa1d5 100644
--- a/SRC/cheevx_2stage.f
+++ b/SRC/cheevx_2stage.f
@@ -334,7 +334,7 @@ SUBROUTINE CHEEVX_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
       LOGICAL            ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
      $                   WANTZ
       CHARACTER          ORDER
-      INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
+      INTEGER            I, IINFO, INDD, INDE, INDEE, INDIBL,
      $                   INDISP, INDIWK, INDRWK, INDTAU, INDWRK, ISCALE,
      $                   ITMP1, J, JJ, LLWORK,
      $                   NSPLIT, LWMIN, LHTRD, LWTRD, KD, IB, INDHOUS
@@ -470,8 +470,11 @@ SUBROUTINE CHEEVX_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -580,14 +583,8 @@ SUBROUTINE CHEEVX_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    40 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL SSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/chpevx.f b/SRC/chpevx.f
index 7bf5510b38..428084bf34 100644
--- a/SRC/chpevx.f
+++ b/SRC/chpevx.f
@@ -266,7 +266,7 @@ SUBROUTINE CHPEVX( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU,
 *     .. Local Scalars ..
       LOGICAL            ALLEIG, INDEIG, TEST, VALEIG, WANTZ
       CHARACTER          ORDER
-      INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE,
+      INTEGER            I, IINFO, INDD, INDE, INDEE,
      $                   INDISP, INDIWK, INDRWK, INDTAU, INDWRK, ISCALE,
      $                   ITMP1, J, JJ, NSPLIT
       REAL               ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
@@ -373,8 +373,11 @@ SUBROUTINE CHPEVX( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          CALL CSSCAL( ( N*( N+1 ) ) / 2, SIGMA, AP, 1 )
@@ -469,14 +472,8 @@ SUBROUTINE CHPEVX( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    20 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL SSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/cstein.f b/SRC/cstein.f
index 1040fe888c..974b790404 100644
--- a/SRC/cstein.f
+++ b/SRC/cstein.f
@@ -95,6 +95,8 @@
 *>          the first submatrix from the top, =2 if W(i) belongs to
 *>          the second submatrix, etc.  ( The output array IBLOCK
 *>          from SSTEBZ is expected here. )
+*>          A negative entry, with which SSTEBZ flags an eigenvalue that
+*>          did not converge, is treated as its absolute value.
 *> \endverbatim
 *>
 *> \param[in] ISPLIT
@@ -244,11 +246,12 @@ SUBROUTINE CSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
          INFO = -9
       ELSE
          DO 20 J = 2, M
-            IF( IBLOCK( J ).LT.IBLOCK( J-1 ) ) THEN
+            IF( ABS( IBLOCK( J ) ).LT.ABS( IBLOCK( J-1 ) ) ) THEN
                INFO = -6
                GO TO 30
             END IF
-            IF( IBLOCK( J ).EQ.IBLOCK( J-1 ) .AND. W( J ).LT.W( J-1 ) )
+            IF( ABS( IBLOCK( J ) ).EQ.ABS( IBLOCK( J-1 ) ) .AND.
+     $          W( J ).LT.W( J-1 ) )
      $           THEN
                INFO = -5
                GO TO 30
@@ -292,7 +295,7 @@ SUBROUTINE CSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
 *     Compute eigenvectors of matrix blocks.
 *
       J1 = 1
-      DO 180 NBLK = 1, IBLOCK( M )
+      DO 180 NBLK = 1, ABS( IBLOCK( M ) )
 *
 *        Find starting and ending indices of block nblk.
 *
@@ -324,7 +327,7 @@ SUBROUTINE CSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
    60    CONTINUE
          JBLK = 0
          DO 170 J = J1, M
-            IF( IBLOCK( J ).NE.NBLK ) THEN
+            IF( ABS( IBLOCK( J ) ).NE.NBLK ) THEN
                J1 = J
                GO TO 180
             END IF
diff --git a/SRC/dsbevx.f b/SRC/dsbevx.f
index 6b29432e01..20ea91e786 100644
--- a/SRC/dsbevx.f
+++ b/SRC/dsbevx.f
@@ -290,7 +290,7 @@ SUBROUTINE DSBEVX( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ,
 *     .. Local Scalars ..
       LOGICAL            ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ
       CHARACTER          ORDER
-      INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
+      INTEGER            I, IINFO, INDD, INDE, INDEE, INDIBL,
      $                   INDISP, INDIWO, INDWRK, ISCALE, ITMP1, J, JJ,
      $                   NSPLIT
       DOUBLE PRECISION   ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
@@ -406,8 +406,11 @@ SUBROUTINE DSBEVX( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -507,14 +510,8 @@ SUBROUTINE DSBEVX( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    30 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL DSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/dsbevx_2stage.f b/SRC/dsbevx_2stage.f
index 233a7121e5..b54675f35f 100644
--- a/SRC/dsbevx_2stage.f
+++ b/SRC/dsbevx_2stage.f
@@ -349,7 +349,7 @@ SUBROUTINE DSBEVX_2STAGE( JOBZ, RANGE, UPLO, N, KD, AB, LDAB,
       LOGICAL            ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ,
      $                   LQUERY
       CHARACTER          ORDER
-      INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
+      INTEGER            I, IINFO, INDD, INDE, INDEE, INDIBL,
      $                   INDISP, INDIWO, INDWRK, ISCALE, ITMP1, J, JJ,
      $                   LLWORK, LWMIN, LHTRD, LWTRD, IB, INDHOUS,
      $                   NSPLIT
@@ -490,8 +490,11 @@ SUBROUTINE DSBEVX_2STAGE( JOBZ, RANGE, UPLO, N, KD, AB, LDAB,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -596,14 +599,8 @@ SUBROUTINE DSBEVX_2STAGE( JOBZ, RANGE, UPLO, N, KD, AB, LDAB,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    30 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL DSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/dspevx.f b/SRC/dspevx.f
index 6eaff9b179..11c71990f6 100644
--- a/SRC/dspevx.f
+++ b/SRC/dspevx.f
@@ -257,7 +257,7 @@ SUBROUTINE DSPEVX( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU,
 *     .. Local Scalars ..
       LOGICAL            ALLEIG, INDEIG, TEST, VALEIG, WANTZ
       CHARACTER          ORDER
-      INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE,
+      INTEGER            I, IINFO, INDD, INDE, INDEE,
      $                   INDISP, INDIWO, INDTAU, INDWRK, ISCALE, ITMP1,
      $                   J, JJ, NSPLIT
       DOUBLE PRECISION   ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
@@ -364,8 +364,11 @@ SUBROUTINE DSPEVX( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          CALL DSCAL( ( N*( N+1 ) ) / 2, SIGMA, AP, 1 )
@@ -458,14 +461,8 @@ SUBROUTINE DSPEVX( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    20 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL DSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/dstein.f b/SRC/dstein.f
index 02de3a0b37..bb7eb0eb39 100644
--- a/SRC/dstein.f
+++ b/SRC/dstein.f
@@ -88,6 +88,8 @@
 *>          the first submatrix from the top, =2 if W(i) belongs to
 *>          the second submatrix, etc.  ( The output array IBLOCK
 *>          from DSTEBZ is expected here. )
+*>          A negative entry, with which DSTEBZ flags an eigenvalue that
+*>          did not converge, is treated as its absolute value.
 *> \endverbatim
 *>
 *> \param[in] ISPLIT
@@ -233,11 +235,12 @@ SUBROUTINE DSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
          INFO = -9
       ELSE
          DO 20 J = 2, M
-            IF( IBLOCK( J ).LT.IBLOCK( J-1 ) ) THEN
+            IF( ABS( IBLOCK( J ) ).LT.ABS( IBLOCK( J-1 ) ) ) THEN
                INFO = -6
                GO TO 30
             END IF
-            IF( IBLOCK( J ).EQ.IBLOCK( J-1 ) .AND. W( J ).LT.W( J-1 ) )
+            IF( ABS( IBLOCK( J ) ).EQ.ABS( IBLOCK( J-1 ) ) .AND.
+     $          W( J ).LT.W( J-1 ) )
      $           THEN
                INFO = -5
                GO TO 30
@@ -281,7 +284,7 @@ SUBROUTINE DSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
 *     Compute eigenvectors of matrix blocks.
 *
       J1 = 1
-      DO 160 NBLK = 1, IBLOCK( M )
+      DO 160 NBLK = 1, ABS( IBLOCK( M ) )
 *
 *        Find starting and ending indices of block nblk.
 *
@@ -313,7 +316,7 @@ SUBROUTINE DSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
    60    CONTINUE
          JBLK = 0
          DO 150 J = J1, M
-            IF( IBLOCK( J ).NE.NBLK ) THEN
+            IF( ABS( IBLOCK( J ) ).NE.NBLK ) THEN
                J1 = J
                GO TO 160
             END IF
diff --git a/SRC/dstevr.f b/SRC/dstevr.f
index d75afcb375..87548af613 100644
--- a/SRC/dstevr.f
+++ b/SRC/dstevr.f
@@ -326,7 +326,7 @@ SUBROUTINE DSTEVR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
       LOGICAL            ALLEIG, INDEIG, TEST, LQUERY, VALEIG, WANTZ,
      $                   TRYRAC
       CHARACTER          ORDER
-      INTEGER            I, IEEEOK, IMAX, INDIBL, INDIFL, INDISP,
+      INTEGER            I, IEEEOK, INDIBL, INDIFL, INDISP,
      $                   INDIWO, ISCALE, ITMP1, J, JJ, LIWMIN, LWMIN,
      $                   NSPLIT
       DOUBLE PRECISION   BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
@@ -450,8 +450,11 @@ SUBROUTINE DSTEVR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
          ISCALE = 1
          SIGMA = RMIN / TNRM
       ELSE IF( TNRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / TNRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          CALL DSCAL( N, SIGMA, D, 1 )
@@ -537,14 +540,8 @@ SUBROUTINE DSTEVR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    10 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL DSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/dstevx.f b/SRC/dstevx.f
index e18c032651..e9ddbd8066 100644
--- a/SRC/dstevx.f
+++ b/SRC/dstevx.f
@@ -251,7 +251,7 @@ SUBROUTINE DSTEVX( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
 *     .. Local Scalars ..
       LOGICAL            ALLEIG, INDEIG, TEST, VALEIG, WANTZ
       CHARACTER          ORDER
-      INTEGER            I, IINFO, IMAX, INDISP, INDIWO, INDWRK,
+      INTEGER            I, IINFO, INDISP, INDIWO, INDWRK,
      $                   ISCALE, ITMP1, J, JJ, NSPLIT
       DOUBLE PRECISION   BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
      $                   TMP1, TNRM, VLL, VUU
@@ -352,8 +352,11 @@ SUBROUTINE DSTEVX( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
          ISCALE = 1
          SIGMA = RMIN / TNRM
       ELSE IF( TNRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / TNRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          CALL DSCAL( N, SIGMA, D, 1 )
@@ -427,14 +430,8 @@ SUBROUTINE DSTEVX( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    20 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL DSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/dsyevr.f b/SRC/dsyevr.f
index 92166fa1bd..caa54cda31 100644
--- a/SRC/dsyevr.f
+++ b/SRC/dsyevr.f
@@ -371,7 +371,7 @@ SUBROUTINE DSYEVR( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
       LOGICAL            ALLEIG, INDEIG, LOWER, LQUERY, VALEIG, WANTZ,
      $                   TRYRAC
       CHARACTER          ORDER
-      INTEGER            I, IEEEOK, IINFO, IMAX, INDD, INDDD, INDE,
+      INTEGER            I, IEEEOK, IINFO, INDD, INDDD, INDE,
      $                   INDEE, INDIBL, INDIFL, INDISP, INDIWO, INDTAU,
      $                   INDWK, INDWKN, ISCALE, J, JJ, LIWMIN,
      $                   LLWORK, LLWRKN, LWKOPT, LWMIN, NB, NSPLIT
@@ -511,8 +511,11 @@ SUBROUTINE DSYEVR( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -655,14 +658,8 @@ SUBROUTINE DSYEVR( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
 *
 *  Jump here if DSTEMR/DSTEIN succeeded.
    30 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL DSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.  Note: We do not sort the IFAIL portion of IWORK.
diff --git a/SRC/dsyevr_2stage.f b/SRC/dsyevr_2stage.f
index 180800e265..f09a137b3e 100644
--- a/SRC/dsyevr_2stage.f
+++ b/SRC/dsyevr_2stage.f
@@ -405,7 +405,7 @@ SUBROUTINE DSYEVR_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
       LOGICAL            ALLEIG, INDEIG, LOWER, LQUERY, VALEIG, WANTZ,
      $                   TRYRAC
       CHARACTER          ORDER
-      INTEGER            I, IEEEOK, IINFO, IMAX, INDD, INDDD, INDE,
+      INTEGER            I, IEEEOK, IINFO, INDD, INDDD, INDE,
      $                   INDEE, INDIBL, INDIFL, INDISP, INDIWO, INDTAU,
      $                   INDWK, INDWKN, ISCALE, J, JJ, LIWMIN,
      $                   LLWORK, LLWRKN, LWMIN, NSPLIT,
@@ -556,8 +556,11 @@ SUBROUTINE DSYEVR_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -705,14 +708,8 @@ SUBROUTINE DSYEVR_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
 *
 *  Jump here if DSTEMR/DSTEIN succeeded.
    30 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL DSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.  Note: We do not sort the IFAIL portion of IWORK.
diff --git a/SRC/dsyevx.f b/SRC/dsyevx.f
index 84dfe289d0..ac3bf7d28f 100644
--- a/SRC/dsyevx.f
+++ b/SRC/dsyevx.f
@@ -278,7 +278,7 @@ SUBROUTINE DSYEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
       LOGICAL            ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
      $                   WANTZ
       CHARACTER          ORDER
-      INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
+      INTEGER            I, IINFO, INDD, INDE, INDEE, INDIBL,
      $                   INDISP, INDIWO, INDTAU, INDWKN, INDWRK, ISCALE,
      $                   ITMP1, J, JJ, LLWORK, LLWRKN, LWKMIN,
      $                   LWKOPT, NB, NSPLIT
@@ -407,8 +407,11 @@ SUBROUTINE DSYEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -514,14 +517,8 @@ SUBROUTINE DSYEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    40 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL DSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/dsyevx_2stage.f b/SRC/dsyevx_2stage.f
index 4b5519be91..cdbb7e2c1e 100644
--- a/SRC/dsyevx_2stage.f
+++ b/SRC/dsyevx_2stage.f
@@ -325,7 +325,7 @@ SUBROUTINE DSYEVX_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
       LOGICAL            ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
      $                   WANTZ
       CHARACTER          ORDER
-      INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
+      INTEGER            I, IINFO, INDD, INDE, INDEE, INDIBL,
      $                   INDISP, INDIWO, INDTAU, INDWKN, INDWRK, ISCALE,
      $                   ITMP1, J, JJ, LLWORK, LLWRKN,
      $                   NSPLIT, LWMIN, LHTRD, LWTRD, KD, IB, INDHOUS
@@ -459,8 +459,11 @@ SUBROUTINE DSYEVX_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -569,14 +572,8 @@ SUBROUTINE DSYEVX_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    40 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL DSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/ssbevx.f b/SRC/ssbevx.f
index 313616e7d5..c7080f906a 100644
--- a/SRC/ssbevx.f
+++ b/SRC/ssbevx.f
@@ -290,7 +290,7 @@ SUBROUTINE SSBEVX( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ,
 *     .. Local Scalars ..
       LOGICAL            ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ
       CHARACTER          ORDER
-      INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
+      INTEGER            I, IINFO, INDD, INDE, INDEE, INDIBL,
      $                   INDISP, INDIWO, INDWRK, ISCALE, ITMP1, J, JJ,
      $                   NSPLIT
       REAL               ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
@@ -406,8 +406,11 @@ SUBROUTINE SSBEVX( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -507,14 +510,8 @@ SUBROUTINE SSBEVX( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    30 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL SSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/ssbevx_2stage.f b/SRC/ssbevx_2stage.f
index 228839cb52..d70644ca22 100644
--- a/SRC/ssbevx_2stage.f
+++ b/SRC/ssbevx_2stage.f
@@ -349,7 +349,7 @@ SUBROUTINE SSBEVX_2STAGE( JOBZ, RANGE, UPLO, N, KD, AB, LDAB,
       LOGICAL            ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ,
      $                   LQUERY
       CHARACTER          ORDER
-      INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
+      INTEGER            I, IINFO, INDD, INDE, INDEE, INDIBL,
      $                   INDISP, INDIWO, INDWRK, ISCALE, ITMP1, J, JJ,
      $                   LLWORK, LWMIN, LHTRD, LWTRD, IB, INDHOUS,
      $                   NSPLIT
@@ -491,8 +491,11 @@ SUBROUTINE SSBEVX_2STAGE( JOBZ, RANGE, UPLO, N, KD, AB, LDAB,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -597,14 +600,8 @@ SUBROUTINE SSBEVX_2STAGE( JOBZ, RANGE, UPLO, N, KD, AB, LDAB,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    30 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL SSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/sspevx.f b/SRC/sspevx.f
index 375a6a223d..ca602d8c7d 100644
--- a/SRC/sspevx.f
+++ b/SRC/sspevx.f
@@ -257,7 +257,7 @@ SUBROUTINE SSPEVX( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU,
 *     .. Local Scalars ..
       LOGICAL            ALLEIG, INDEIG, TEST, VALEIG, WANTZ
       CHARACTER          ORDER
-      INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE,
+      INTEGER            I, IINFO, INDD, INDE, INDEE,
      $                   INDISP, INDIWO, INDTAU, INDWRK, ISCALE, ITMP1,
      $                   J, JJ, NSPLIT
       REAL               ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
@@ -364,8 +364,11 @@ SUBROUTINE SSPEVX( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          CALL SSCAL( ( N*( N+1 ) ) / 2, SIGMA, AP, 1 )
@@ -458,14 +461,8 @@ SUBROUTINE SSPEVX( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    20 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL SSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/sstein.f b/SRC/sstein.f
index 471d31220d..72d5507309 100644
--- a/SRC/sstein.f
+++ b/SRC/sstein.f
@@ -88,6 +88,8 @@
 *>          the first submatrix from the top, =2 if W(i) belongs to
 *>          the second submatrix, etc.  ( The output array IBLOCK
 *>          from SSTEBZ is expected here. )
+*>          A negative entry, with which SSTEBZ flags an eigenvalue that
+*>          did not converge, is treated as its absolute value.
 *> \endverbatim
 *>
 *> \param[in] ISPLIT
@@ -233,11 +235,12 @@ SUBROUTINE SSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
          INFO = -9
       ELSE
          DO 20 J = 2, M
-            IF( IBLOCK( J ).LT.IBLOCK( J-1 ) ) THEN
+            IF( ABS( IBLOCK( J ) ).LT.ABS( IBLOCK( J-1 ) ) ) THEN
                INFO = -6
                GO TO 30
             END IF
-            IF( IBLOCK( J ).EQ.IBLOCK( J-1 ) .AND. W( J ).LT.W( J-1 ) )
+            IF( ABS( IBLOCK( J ) ).EQ.ABS( IBLOCK( J-1 ) ) .AND.
+     $          W( J ).LT.W( J-1 ) )
      $           THEN
                INFO = -5
                GO TO 30
@@ -281,7 +284,7 @@ SUBROUTINE SSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
 *     Compute eigenvectors of matrix blocks.
 *
       J1 = 1
-      DO 160 NBLK = 1, IBLOCK( M )
+      DO 160 NBLK = 1, ABS( IBLOCK( M ) )
 *
 *        Find starting and ending indices of block nblk.
 *
@@ -313,7 +316,7 @@ SUBROUTINE SSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
    60    CONTINUE
          JBLK = 0
          DO 150 J = J1, M
-            IF( IBLOCK( J ).NE.NBLK ) THEN
+            IF( ABS( IBLOCK( J ) ).NE.NBLK ) THEN
                J1 = J
                GO TO 160
             END IF
diff --git a/SRC/sstevr.f b/SRC/sstevr.f
index 0b72c2f0e5..51dac1c905 100644
--- a/SRC/sstevr.f
+++ b/SRC/sstevr.f
@@ -328,7 +328,7 @@ SUBROUTINE SSTEVR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
       LOGICAL            ALLEIG, INDEIG, TEST, LQUERY, VALEIG, WANTZ,
      $                   TRYRAC
       CHARACTER          ORDER
-      INTEGER            I, IEEEOK, IMAX, INDIBL, INDIFL, INDISP,
+      INTEGER            I, IEEEOK, INDIBL, INDIFL, INDISP,
      $                   INDIWO, ISCALE, J, JJ, LIWMIN, LWMIN, NSPLIT
       REAL               BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
      $                   TMP1, TNRM, VLL, VUU
@@ -452,8 +452,11 @@ SUBROUTINE SSTEVR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
          ISCALE = 1
          SIGMA = RMIN / TNRM
       ELSE IF( TNRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / TNRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          CALL SSCAL( N, SIGMA, D, 1 )
@@ -539,14 +542,8 @@ SUBROUTINE SSTEVR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    10 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL SSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/sstevx.f b/SRC/sstevx.f
index f84e740ef2..51c563109b 100644
--- a/SRC/sstevx.f
+++ b/SRC/sstevx.f
@@ -251,7 +251,7 @@ SUBROUTINE SSTEVX( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
 *     .. Local Scalars ..
       LOGICAL            ALLEIG, INDEIG, TEST, VALEIG, WANTZ
       CHARACTER          ORDER
-      INTEGER            I, IINFO, IMAX, INDISP, INDIWO, INDWRK,
+      INTEGER            I, IINFO, INDISP, INDIWO, INDWRK,
      $                   ISCALE, ITMP1, J, JJ, NSPLIT
       REAL               BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
      $                   TMP1, TNRM, VLL, VUU
@@ -352,8 +352,11 @@ SUBROUTINE SSTEVX( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
          ISCALE = 1
          SIGMA = RMIN / TNRM
       ELSE IF( TNRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / TNRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          CALL SSCAL( N, SIGMA, D, 1 )
@@ -427,14 +430,8 @@ SUBROUTINE SSTEVX( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    20 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL SSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/ssyevr.f b/SRC/ssyevr.f
index 6a848bb1cb..cbe3e71ae3 100644
--- a/SRC/ssyevr.f
+++ b/SRC/ssyevr.f
@@ -373,7 +373,7 @@ SUBROUTINE SSYEVR( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
       LOGICAL            ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
      $                   WANTZ, TRYRAC
       CHARACTER          ORDER
-      INTEGER            I, IEEEOK, IINFO, IMAX, INDD, INDDD, INDE,
+      INTEGER            I, IEEEOK, IINFO, INDD, INDDD, INDE,
      $                   INDEE, INDIBL, INDIFL, INDISP, INDIWO, INDTAU,
      $                   INDWK, INDWKN, ISCALE, J, JJ, LIWMIN,
      $                   LLWORK, LLWRKN, LWKOPT, LWMIN, NB, NSPLIT
@@ -516,8 +516,11 @@ SUBROUTINE SSYEVR( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -665,14 +668,8 @@ SUBROUTINE SSYEVR( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
 *
 *  Jump here if SSTEMR/SSTEIN succeeded.
    30 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL SSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.  Note: We do not sort the IFAIL portion of IWORK.
diff --git a/SRC/ssyevr_2stage.f b/SRC/ssyevr_2stage.f
index 63f0023886..47d82f3084 100644
--- a/SRC/ssyevr_2stage.f
+++ b/SRC/ssyevr_2stage.f
@@ -405,7 +405,7 @@ SUBROUTINE SSYEVR_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
       LOGICAL            ALLEIG, INDEIG, LOWER, LQUERY, VALEIG, WANTZ,
      $                   TRYRAC, TEST
       CHARACTER          ORDER
-      INTEGER            I, IEEEOK, IINFO, IMAX, INDD, INDDD, INDE,
+      INTEGER            I, IEEEOK, IINFO, INDD, INDDD, INDE,
      $                   INDEE, INDIBL, INDIFL, INDISP, INDIWO, INDTAU,
      $                   INDWK, INDWKN, ISCALE, J, JJ, LIWMIN,
      $                   LLWORK, LLWRKN, LWMIN, NSPLIT,
@@ -557,8 +557,11 @@ SUBROUTINE SSYEVR_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -711,14 +714,8 @@ SUBROUTINE SSYEVR_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
 *
 *  Jump here if SSTEMR/SSTEIN succeeded.
    30 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL SSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.  Note: We do not sort the IFAIL portion of IWORK.
diff --git a/SRC/ssyevx.f b/SRC/ssyevx.f
index 9528c66fc6..857e61da01 100644
--- a/SRC/ssyevx.f
+++ b/SRC/ssyevx.f
@@ -278,7 +278,7 @@ SUBROUTINE SSYEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
       LOGICAL            ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
      $                   WANTZ
       CHARACTER          ORDER
-      INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
+      INTEGER            I, IINFO, INDD, INDE, INDEE, INDIBL,
      $                   INDISP, INDIWO, INDTAU, INDWKN, INDWRK, ISCALE,
      $                   ITMP1, J, JJ, LLWORK, LLWRKN, LWKMIN,
      $                   LWKOPT, NB, NSPLIT
@@ -408,8 +408,11 @@ SUBROUTINE SSYEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -515,14 +518,8 @@ SUBROUTINE SSYEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    40 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL SSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/ssyevx_2stage.f b/SRC/ssyevx_2stage.f
index 8a1fc2798c..a0338bac85 100644
--- a/SRC/ssyevx_2stage.f
+++ b/SRC/ssyevx_2stage.f
@@ -325,7 +325,7 @@ SUBROUTINE SSYEVX_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
       LOGICAL            ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
      $                   WANTZ
       CHARACTER          ORDER
-      INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
+      INTEGER            I, IINFO, INDD, INDE, INDEE, INDIBL,
      $                   INDISP, INDIWO, INDTAU, INDWKN, INDWRK, ISCALE,
      $                   ITMP1, J, JJ, LLWORK, LLWRKN,
      $                   NSPLIT, LWMIN, LHTRD, LWTRD, KD, IB, INDHOUS
@@ -460,8 +460,11 @@ SUBROUTINE SSYEVX_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -570,14 +573,8 @@ SUBROUTINE SSYEVX_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    40 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL SSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/zhbevx.f b/SRC/zhbevx.f
index 22487eb27e..9387893844 100644
--- a/SRC/zhbevx.f
+++ b/SRC/zhbevx.f
@@ -296,7 +296,7 @@ SUBROUTINE ZHBEVX( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ,
 *     .. Local Scalars ..
       LOGICAL            ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ
       CHARACTER          ORDER
-      INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
+      INTEGER            I, IINFO, INDD, INDE, INDEE, INDIBL,
      $                   INDISP, INDIWK, INDRWK, INDWRK, ISCALE, ITMP1,
      $                   J, JJ, NSPLIT
       DOUBLE PRECISION   ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
@@ -415,8 +415,11 @@ SUBROUTINE ZHBEVX( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -517,14 +520,8 @@ SUBROUTINE ZHBEVX( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    30 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL DSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/zhbevx_2stage.f b/SRC/zhbevx_2stage.f
index a22b589f51..bfc6e542a7 100644
--- a/SRC/zhbevx_2stage.f
+++ b/SRC/zhbevx_2stage.f
@@ -357,7 +357,7 @@ SUBROUTINE ZHBEVX_2STAGE( JOBZ, RANGE, UPLO, N, KD, AB, LDAB,
       LOGICAL            ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ,
      $                   LQUERY
       CHARACTER          ORDER
-      INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
+      INTEGER            I, IINFO, INDD, INDE, INDEE, INDIBL,
      $                   INDISP, INDIWK, INDRWK, INDWRK, ISCALE, ITMP1,
      $                   LLWORK, LWMIN, LHTRD, LWTRD, IB, INDHOUS,
      $                   J, JJ, NSPLIT
@@ -500,8 +500,11 @@ SUBROUTINE ZHBEVX_2STAGE( JOBZ, RANGE, UPLO, N, KD, AB, LDAB,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -607,14 +610,8 @@ SUBROUTINE ZHBEVX_2STAGE( JOBZ, RANGE, UPLO, N, KD, AB, LDAB,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    30 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL DSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/zheevr.f b/SRC/zheevr.f
index 038738ec8b..cf77d1c2b9 100644
--- a/SRC/zheevr.f
+++ b/SRC/zheevr.f
@@ -398,7 +398,7 @@ SUBROUTINE ZHEEVR( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
       LOGICAL            ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
      $                   WANTZ, TRYRAC
       CHARACTER          ORDER
-      INTEGER            I, IEEEOK, IINFO, IMAX, INDIBL, INDIFL, INDISP,
+      INTEGER            I, IEEEOK, IINFO, INDIBL, INDIFL, INDISP,
      $                   INDIWO, INDRD, INDRDD, INDRE, INDREE, INDRWK,
      $                   INDTAU, INDWK, INDWKN, ISCALE, ITMP1, J, JJ,
      $                   LIWMIN, LLWORK, LLRWORK, LLWRKN, LRWMIN,
@@ -548,8 +548,11 @@ SUBROUTINE ZHEEVR( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -699,14 +702,8 @@ SUBROUTINE ZHEEVR( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    30 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL DSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/zheevr_2stage.f b/SRC/zheevr_2stage.f
index 0ba5a29533..4553203a7e 100644
--- a/SRC/zheevr_2stage.f
+++ b/SRC/zheevr_2stage.f
@@ -433,7 +433,7 @@ SUBROUTINE ZHEEVR_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
       LOGICAL            ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
      $                   WANTZ, TRYRAC
       CHARACTER          ORDER
-      INTEGER            I, IEEEOK, IINFO, IMAX, INDIBL, INDIFL, INDISP,
+      INTEGER            I, IEEEOK, IINFO, INDIBL, INDIFL, INDISP,
      $                   INDIWO, INDRD, INDRDD, INDRE, INDREE, INDRWK,
      $                   INDTAU, INDWK, INDWKN, ISCALE, ITMP1, J, JJ,
      $                   LIWMIN, LLWORK, LLRWORK, LLWRKN, LRWMIN,
@@ -590,8 +590,11 @@ SUBROUTINE ZHEEVR_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -744,14 +747,8 @@ SUBROUTINE ZHEEVR_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    30 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL DSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/zheevx.f b/SRC/zheevx.f
index f7f58932f6..29699b8fb3 100644
--- a/SRC/zheevx.f
+++ b/SRC/zheevx.f
@@ -287,7 +287,7 @@ SUBROUTINE ZHEEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
       LOGICAL            ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
      $                   WANTZ
       CHARACTER          ORDER
-      INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
+      INTEGER            I, IINFO, INDD, INDE, INDEE, INDIBL,
      $                   INDISP, INDIWK, INDRWK, INDTAU, INDWRK, ISCALE,
      $                   ITMP1, J, JJ, LLWORK, LWKMIN, LWKOPT, NB,
      $                   NSPLIT
@@ -418,8 +418,11 @@ SUBROUTINE ZHEEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -524,14 +527,8 @@ SUBROUTINE ZHEEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    40 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL DSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/zheevx_2stage.f b/SRC/zheevx_2stage.f
index 1b0a22d27f..a9a7e59299 100644
--- a/SRC/zheevx_2stage.f
+++ b/SRC/zheevx_2stage.f
@@ -334,7 +334,7 @@ SUBROUTINE ZHEEVX_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
       LOGICAL            ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
      $                   WANTZ
       CHARACTER          ORDER
-      INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
+      INTEGER            I, IINFO, INDD, INDE, INDEE, INDIBL,
      $                   INDISP, INDIWK, INDRWK, INDTAU, INDWRK, ISCALE,
      $                   ITMP1, J, JJ, LLWORK,
      $                   NSPLIT, LWMIN, LHTRD, LWTRD, KD, IB, INDHOUS
@@ -469,8 +469,11 @@ SUBROUTINE ZHEEVX_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          IF( LOWER ) THEN
@@ -579,14 +582,8 @@ SUBROUTINE ZHEEVX_2STAGE( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    40 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL DSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/zhpevx.f b/SRC/zhpevx.f
index 3a60e8fd64..a3a5bcfe9a 100644
--- a/SRC/zhpevx.f
+++ b/SRC/zhpevx.f
@@ -266,7 +266,7 @@ SUBROUTINE ZHPEVX( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU,
 *     .. Local Scalars ..
       LOGICAL            ALLEIG, INDEIG, TEST, VALEIG, WANTZ
       CHARACTER          ORDER
-      INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE,
+      INTEGER            I, IINFO, INDD, INDE, INDEE,
      $                   INDISP, INDIWK, INDRWK, INDTAU, INDWRK, ISCALE,
      $                   ITMP1, J, JJ, NSPLIT
       DOUBLE PRECISION   ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
@@ -373,8 +373,11 @@ SUBROUTINE ZHPEVX( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU,
          ISCALE = 1
          SIGMA = RMIN / ANRM
       ELSE IF( ANRM.GT.RMAX ) THEN
-         ISCALE = 1
+*        An infinite norm gives SIGMA = 0; leave such a matrix
+*        unscaled so that VL and VU remain a valid interval.
          SIGMA = RMAX / ANRM
+         IF( SIGMA.GT.ZERO )
+     $      ISCALE = 1
       END IF
       IF( ISCALE.EQ.1 ) THEN
          CALL ZDSCAL( ( N*( N+1 ) ) / 2, SIGMA, AP, 1 )
@@ -469,14 +472,8 @@ SUBROUTINE ZHPEVX( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU,
 *     If matrix was scaled, then rescale eigenvalues appropriately.
 *
    20 CONTINUE
-      IF( ISCALE.EQ.1 ) THEN
-         IF( INFO.EQ.0 ) THEN
-            IMAX = M
-         ELSE
-            IMAX = INFO - 1
-         END IF
-         CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
-      END IF
+      IF( ISCALE.EQ.1 )
+     $   CALL DSCAL( M, ONE / SIGMA, W, 1 )
 *
 *     If eigenvalues are not in order, then sort them, along with
 *     eigenvectors.
diff --git a/SRC/zstein.f b/SRC/zstein.f
index ca408fe70b..a00c088c1b 100644
--- a/SRC/zstein.f
+++ b/SRC/zstein.f
@@ -95,6 +95,8 @@
 *>          the first submatrix from the top, =2 if W(i) belongs to
 *>          the second submatrix, etc.  ( The output array IBLOCK
 *>          from DSTEBZ is expected here. )
+*>          A negative entry, with which DSTEBZ flags an eigenvalue that
+*>          did not converge, is treated as its absolute value.
 *> \endverbatim
 *>
 *> \param[in] ISPLIT
@@ -244,11 +246,12 @@ SUBROUTINE ZSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
          INFO = -9
       ELSE
          DO 20 J = 2, M
-            IF( IBLOCK( J ).LT.IBLOCK( J-1 ) ) THEN
+            IF( ABS( IBLOCK( J ) ).LT.ABS( IBLOCK( J-1 ) ) ) THEN
                INFO = -6
                GO TO 30
             END IF
-            IF( IBLOCK( J ).EQ.IBLOCK( J-1 ) .AND. W( J ).LT.W( J-1 ) )
+            IF( ABS( IBLOCK( J ) ).EQ.ABS( IBLOCK( J-1 ) ) .AND.
+     $          W( J ).LT.W( J-1 ) )
      $           THEN
                INFO = -5
                GO TO 30
@@ -292,7 +295,7 @@ SUBROUTINE ZSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
 *     Compute eigenvectors of matrix blocks.
 *
       J1 = 1
-      DO 180 NBLK = 1, IBLOCK( M )
+      DO 180 NBLK = 1, ABS( IBLOCK( M ) )
 *
 *        Find starting and ending indices of block nblk.
 *
@@ -324,7 +327,7 @@ SUBROUTINE ZSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
    60    CONTINUE
          JBLK = 0
          DO 170 J = J1, M
-            IF( IBLOCK( J ).NE.NBLK ) THEN
+            IF( ABS( IBLOCK( J ) ).NE.NBLK ) THEN
                J1 = J
                GO TO 180
             END IF
diff --git a/TESTING/EIG/derrst.f b/TESTING/EIG/derrst.f
index 94018743f2..3662873a94 100644
--- a/TESTING/EIG/derrst.f
+++ b/TESTING/EIG/derrst.f
@@ -86,6 +86,7 @@ SUBROUTINE DERRST( PATH, NUNIT )
      $                   E( NMAX ), Q( NMAX, NMAX ), R( NMAX ),
      $                   TAU( NMAX ), W( LW ), X( NMAX ),
      $                   Z( NMAX, NMAX )
+      DOUBLE PRECISION   RINF, RZERO
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAMEN
@@ -439,6 +440,49 @@ SUBROUTINE DERRST( PATH, NUNIT )
          CALL CHKXER( 'DSTEIN', INFOT, NOUT, LERR, OK )
          NT = NT + 4
 *
+*        A matrix that holds an infinity must come back from the
+*        drivers that scale it, instead of reaching XERBLA with a
+*        scale factor of zero.
+*
+         RZERO = 0.0D0
+         RINF = 1.0D0 / RZERO
+         DO 50 J = 1, NMAX
+            D( J ) = DBLE( J )
+            E( J ) = 1.0D0 / DBLE( J+1 )
+            DO 40 I = 1, NMAX
+               A( I, J ) = RZERO
+   40       CONTINUE
+            A( J, J ) = D( J )
+   50    CONTINUE
+         D( 1 ) = RINF
+         A( 1, 1 ) = RINF
+         E( NMAX ) = RZERO
+         A( 1, 2 ) = E( 1 )
+         A( 2, 3 ) = E( 2 )
+         SRNAMT = 'DSTEVX'
+         INFOT = 0
+         LERR = .FALSE.
+         CALL DSTEVX( 'V', 'V', NMAX, D, E, RZERO, 1.0D0, 1, NMAX,
+     $                RZERO, M, X, Z, NMAX, W, IW, I3, INFO )
+         IF( LERR ) THEN
+            WRITE( NOUT, FMT = 9997 )'DSTEVX'
+            OK = .FALSE.
+         END IF
+         D( 1 ) = RINF
+         E( 1 ) = 1.0D0 / DBLE( 2 )
+         E( 2 ) = 1.0D0 / DBLE( 3 )
+         E( NMAX ) = RZERO
+         SRNAMT = 'DSYEVX'
+         INFOT = 0
+         LERR = .FALSE.
+         CALL DSYEVX( 'V', 'V', 'U', NMAX, A, NMAX, RZERO, 1.0D0, 1,
+     $                NMAX, RZERO, M, X, Z, NMAX, W, LW, IW, I3, INFO )
+         IF( LERR ) THEN
+            WRITE( NOUT, FMT = 9997 )'DSYEVX'
+            OK = .FALSE.
+         END IF
+         NT = NT + 2
+*
 *        DSTEQR
 *
          SRNAMT = 'DSTEQR'
@@ -1427,6 +1471,8 @@ SUBROUTINE DERRST( PATH, NUNIT )
      $      ' (', I3, ' tests done)' )
  9998 FORMAT( ' *** ', A3, ' routines failed the tests of the error ',
      $      'exits ***' )
+ 9997 FORMAT( ' *** ', A6, ' called XERBLA for a matrix with an ',
+     $      'infinity ***' )
 *
       RETURN
 *
diff --git a/TESTING/EIG/serrst.f b/TESTING/EIG/serrst.f
index 471fa35993..642a0648cb 100644
--- a/TESTING/EIG/serrst.f
+++ b/TESTING/EIG/serrst.f
@@ -86,6 +86,7 @@ SUBROUTINE SERRST( PATH, NUNIT )
      $                   E( NMAX ), Q( NMAX, NMAX ), R( NMAX ),
      $                   TAU( NMAX ), W( LW ), X( NMAX ),
      $                   Z( NMAX, NMAX )
+      REAL               RINF, RZERO
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAMEN
@@ -439,6 +440,49 @@ SUBROUTINE SERRST( PATH, NUNIT )
          CALL CHKXER( 'SSTEIN', INFOT, NOUT, LERR, OK )
          NT = NT + 4
 *
+*        A matrix that holds an infinity must come back from the
+*        drivers that scale it, instead of reaching XERBLA with a
+*        scale factor of zero.
+*
+         RZERO = 0.0E0
+         RINF = 1.0E0 / RZERO
+         DO 50 J = 1, NMAX
+            D( J ) = REAL( J )
+            E( J ) = 1.0E0 / REAL( J+1 )
+            DO 40 I = 1, NMAX
+               A( I, J ) = RZERO
+   40       CONTINUE
+            A( J, J ) = D( J )
+   50    CONTINUE
+         D( 1 ) = RINF
+         A( 1, 1 ) = RINF
+         E( NMAX ) = RZERO
+         A( 1, 2 ) = E( 1 )
+         A( 2, 3 ) = E( 2 )
+         SRNAMT = 'SSTEVX'
+         INFOT = 0
+         LERR = .FALSE.
+         CALL SSTEVX( 'V', 'V', NMAX, D, E, RZERO, 1.0E0, 1, NMAX,
+     $                RZERO, M, X, Z, NMAX, W, IW, I3, INFO )
+         IF( LERR ) THEN
+            WRITE( NOUT, FMT = 9997 )'SSTEVX'
+            OK = .FALSE.
+         END IF
+         D( 1 ) = RINF
+         E( 1 ) = 1.0E0 / REAL( 2 )
+         E( 2 ) = 1.0E0 / REAL( 3 )
+         E( NMAX ) = RZERO
+         SRNAMT = 'SSYEVX'
+         INFOT = 0
+         LERR = .FALSE.
+         CALL SSYEVX( 'V', 'V', 'U', NMAX, A, NMAX, RZERO, 1.0E0, 1,
+     $                NMAX, RZERO, M, X, Z, NMAX, W, LW, IW, I3, INFO )
+         IF( LERR ) THEN
+            WRITE( NOUT, FMT = 9997 )'SSYEVX'
+            OK = .FALSE.
+         END IF
+         NT = NT + 2
+*
 *        SSTEQR
 *
          SRNAMT = 'SSTEQR'
@@ -1425,6 +1469,8 @@ SUBROUTINE SERRST( PATH, NUNIT )
      $      ' (', I3, ' tests done)' )
  9998 FORMAT( ' *** ', A3, ' routines failed the tests of the error ',
      $      'exits ***' )
+ 9997 FORMAT( ' *** ', A6, ' called XERBLA for a matrix with an ',
+     $      'infinity ***' )
 *
       RETURN
 *