Correct bug0096 reported by Joseph Young from Sandia.

Followed recommendation, use the existing sorting code.

Report sent to LAPACK mailing list on June 26th 2012

From Joseph:
>     There appears to be an inconsistency and possible bug in the dstemr 
> implementation.  When calculating the eigenvalues of a matrix, the 
> returned eigenvalues are supposed to be returned in ascending order.  
> Although this appears to be the case for N >= 3, it does not appear to 
> be the case for N=2.  I believe this happens because the dstemr routine 
> has special cases for N=0,1, and 2, which immediately return after their 
> computation.  Because these cases return immediately, they do not call 
> the sorting routines around line 723 (in LAPACK version 3.4.1).  As 
> such, a simple fix would be to have the N=2 case call this sorting code 
> rather than returning. 

4 files changed, 560 insertions, 557 deletions

diff --git a/SRC/cstemr.f b/SRC/cstemr.f
index c38da6ba..7137fae9 100644
--- a/SRC/cstemr.f
+++ b/SRC/cstemr.f
@@ -560,184 +560,184 @@
                END IF
             ENDIF
          ENDIF
-         RETURN
-      END IF
+      ELSE
 
-*     Continue with general N
+*        Continue with general N
 
-      INDGRS = 1
-      INDERR = 2*N + 1
-      INDGP = 3*N + 1
-      INDD = 4*N + 1
-      INDE2 = 5*N + 1
-      INDWRK = 6*N + 1
-*
-      IINSPL = 1
-      IINDBL = N + 1
-      IINDW = 2*N + 1
-      IINDWK = 3*N + 1
-*
-*     Scale matrix to allowable range, if necessary.
-*     The allowable range is related to the PIVMIN parameter; see the
-*     comments in SLARRD.  The preference for scaling small values
-*     up is heuristic; we expect users' matrices not to be close to the
-*     RMAX threshold.
-*
-      SCALE = ONE
-      TNRM = SLANST( 'M', N, D, E )
-      IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
-         SCALE = RMIN / TNRM
-      ELSE IF( TNRM.GT.RMAX ) THEN
-         SCALE = RMAX / TNRM
-      END IF
-      IF( SCALE.NE.ONE ) THEN
-         CALL SSCAL( N, SCALE, D, 1 )
-         CALL SSCAL( N-1, SCALE, E, 1 )
-         TNRM = TNRM*SCALE
-         IF( VALEIG ) THEN
-*           If eigenvalues in interval have to be found,
-*           scale (WL, WU] accordingly
-            WL = WL*SCALE
-            WU = WU*SCALE
-         ENDIF
-      END IF
+         INDGRS = 1
+         INDERR = 2*N + 1
+         INDGP = 3*N + 1
+         INDD = 4*N + 1
+         INDE2 = 5*N + 1
+         INDWRK = 6*N + 1
+*
+         IINSPL = 1
+         IINDBL = N + 1
+         IINDW = 2*N + 1
+         IINDWK = 3*N + 1
+*
+*        Scale matrix to allowable range, if necessary.
+*        The allowable range is related to the PIVMIN parameter; see the
+*        comments in SLARRD.  The preference for scaling small values
+*        up is heuristic; we expect users' matrices not to be close to the
+*        RMAX threshold.
+*
+         SCALE = ONE
+         TNRM = SLANST( 'M', N, D, E )
+         IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
+            SCALE = RMIN / TNRM
+         ELSE IF( TNRM.GT.RMAX ) THEN
+            SCALE = RMAX / TNRM
+         END IF
+         IF( SCALE.NE.ONE ) THEN
+            CALL SSCAL( N, SCALE, D, 1 )
+            CALL SSCAL( N-1, SCALE, E, 1 )
+            TNRM = TNRM*SCALE
+            IF( VALEIG ) THEN
+*              If eigenvalues in interval have to be found,
+*              scale (WL, WU] accordingly
+               WL = WL*SCALE
+               WU = WU*SCALE
+            ENDIF
+         END IF
 *
-*     Compute the desired eigenvalues of the tridiagonal after splitting
-*     into smaller subblocks if the corresponding off-diagonal elements
-*     are small
-*     THRESH is the splitting parameter for SLARRE
-*     A negative THRESH forces the old splitting criterion based on the
-*     size of the off-diagonal. A positive THRESH switches to splitting
-*     which preserves relative accuracy.
-*
-      IF( TRYRAC ) THEN
-*        Test whether the matrix warrants the more expensive relative approach.
-         CALL SLARRR( N, D, E, IINFO )
-      ELSE
-*        The user does not care about relative accurately eigenvalues
-         IINFO = -1
-      ENDIF
-*     Set the splitting criterion
-      IF (IINFO.EQ.0) THEN
-         THRESH = EPS
-      ELSE
-         THRESH = -EPS
-*        relative accuracy is desired but T does not guarantee it
-         TRYRAC = .FALSE.
-      ENDIF
+*        Compute the desired eigenvalues of the tridiagonal after splitting
+*        into smaller subblocks if the corresponding off-diagonal elements
+*        are small
+*        THRESH is the splitting parameter for SLARRE
+*        A negative THRESH forces the old splitting criterion based on the
+*        size of the off-diagonal. A positive THRESH switches to splitting
+*        which preserves relative accuracy.
+*
+         IF( TRYRAC ) THEN
+*           Test whether the matrix warrants the more expensive relative approach.
+            CALL SLARRR( N, D, E, IINFO )
+         ELSE
+*           The user does not care about relative accurately eigenvalues
+            IINFO = -1
+         ENDIF
+*        Set the splitting criterion
+         IF (IINFO.EQ.0) THEN
+            THRESH = EPS
+         ELSE
+            THRESH = -EPS
+*           relative accuracy is desired but T does not guarantee it
+            TRYRAC = .FALSE.
+         ENDIF
 *
-      IF( TRYRAC ) THEN
-*        Copy original diagonal, needed to guarantee relative accuracy
-         CALL SCOPY(N,D,1,WORK(INDD),1)
-      ENDIF
-*     Store the squares of the offdiagonal values of T
-      DO 5 J = 1, N-1
-         WORK( INDE2+J-1 ) = E(J)**2
+         IF( TRYRAC ) THEN
+*           Copy original diagonal, needed to guarantee relative accuracy
+            CALL SCOPY(N,D,1,WORK(INDD),1)
+         ENDIF
+*        Store the squares of the offdiagonal values of T
+         DO 5 J = 1, N-1
+            WORK( INDE2+J-1 ) = E(J)**2
  5    CONTINUE
 
-*     Set the tolerance parameters for bisection
-      IF( .NOT.WANTZ ) THEN
-*        SLARRE computes the eigenvalues to full precision.
-         RTOL1 = FOUR * EPS
-         RTOL2 = FOUR * EPS
-      ELSE
-*        SLARRE computes the eigenvalues to less than full precision.
-*        CLARRV will refine the eigenvalue approximations, and we only
-*        need less accurate initial bisection in SLARRE.
-*        Note: these settings do only affect the subset case and SLARRE
-         RTOL1 = MAX( SQRT(EPS)*5.0E-2, FOUR * EPS )
-         RTOL2 = MAX( SQRT(EPS)*5.0E-3, FOUR * EPS )
-      ENDIF
-      CALL SLARRE( RANGE, N, WL, WU, IIL, IIU, D, E,
+*        Set the tolerance parameters for bisection
+         IF( .NOT.WANTZ ) THEN
+*           SLARRE computes the eigenvalues to full precision.
+            RTOL1 = FOUR * EPS
+            RTOL2 = FOUR * EPS
+         ELSE
+*           SLARRE computes the eigenvalues to less than full precision.
+*           CLARRV will refine the eigenvalue approximations, and we only
+*           need less accurate initial bisection in SLARRE.
+*           Note: these settings do only affect the subset case and SLARRE
+            RTOL1 = MAX( SQRT(EPS)*5.0E-2, FOUR * EPS )
+            RTOL2 = MAX( SQRT(EPS)*5.0E-3, FOUR * EPS )
+         ENDIF
+         CALL SLARRE( RANGE, N, WL, WU, IIL, IIU, D, E,
      $             WORK(INDE2), RTOL1, RTOL2, THRESH, NSPLIT,
      $             IWORK( IINSPL ), M, W, WORK( INDERR ),
      $             WORK( INDGP ), IWORK( IINDBL ),
      $             IWORK( IINDW ), WORK( INDGRS ), PIVMIN,
      $             WORK( INDWRK ), IWORK( IINDWK ), IINFO )
-      IF( IINFO.NE.0 ) THEN
-         INFO = 10 + ABS( IINFO )
-         RETURN
-      END IF
-*     Note that if RANGE .NE. 'V', SLARRE computes bounds on the desired
-*     part of the spectrum. All desired eigenvalues are contained in
-*     (WL,WU]
+         IF( IINFO.NE.0 ) THEN
+            INFO = 10 + ABS( IINFO )
+            RETURN
+         END IF
+*        Note that if RANGE .NE. 'V', SLARRE computes bounds on the desired
+*        part of the spectrum. All desired eigenvalues are contained in
+*        (WL,WU]
 
 
-      IF( WANTZ ) THEN
+         IF( WANTZ ) THEN
 *
-*        Compute the desired eigenvectors corresponding to the computed
-*        eigenvalues
+*           Compute the desired eigenvectors corresponding to the computed
+*           eigenvalues
 *
-         CALL CLARRV( N, WL, WU, D, E,
+            CALL CLARRV( N, WL, WU, D, E,
      $                PIVMIN, IWORK( IINSPL ), M,
      $                1, M, MINRGP, RTOL1, RTOL2,
      $                W, WORK( INDERR ), WORK( INDGP ), IWORK( IINDBL ),
      $                IWORK( IINDW ), WORK( INDGRS ), Z, LDZ,
      $                ISUPPZ, WORK( INDWRK ), IWORK( IINDWK ), IINFO )
-         IF( IINFO.NE.0 ) THEN
-            INFO = 20 + ABS( IINFO )
-            RETURN
-         END IF
-      ELSE
-*        SLARRE computes eigenvalues of the (shifted) root representation
-*        CLARRV returns the eigenvalues of the unshifted matrix.
-*        However, if the eigenvectors are not desired by the user, we need
-*        to apply the corresponding shifts from SLARRE to obtain the
-*        eigenvalues of the original matrix.
-         DO 20 J = 1, M
-            ITMP = IWORK( IINDBL+J-1 )
-            W( J ) = W( J ) + E( IWORK( IINSPL+ITMP-1 ) )
+            IF( IINFO.NE.0 ) THEN
+               INFO = 20 + ABS( IINFO )
+               RETURN
+            END IF
+         ELSE
+*           SLARRE computes eigenvalues of the (shifted) root representation
+*           CLARRV returns the eigenvalues of the unshifted matrix.
+*           However, if the eigenvectors are not desired by the user, we need
+*           to apply the corresponding shifts from SLARRE to obtain the
+*           eigenvalues of the original matrix.
+            DO 20 J = 1, M
+               ITMP = IWORK( IINDBL+J-1 )
+               W( J ) = W( J ) + E( IWORK( IINSPL+ITMP-1 ) )
  20      CONTINUE
-      END IF
+         END IF
 *
 
-      IF ( TRYRAC ) THEN
-*        Refine computed eigenvalues so that they are relatively accurate
-*        with respect to the original matrix T.
-         IBEGIN = 1
-         WBEGIN = 1
-         DO 39  JBLK = 1, IWORK( IINDBL+M-1 )
-            IEND = IWORK( IINSPL+JBLK-1 )
-            IN = IEND - IBEGIN + 1
-            WEND = WBEGIN - 1
-*           check if any eigenvalues have to be refined in this block
+         IF ( TRYRAC ) THEN
+*           Refine computed eigenvalues so that they are relatively accurate
+*           with respect to the original matrix T.
+            IBEGIN = 1
+            WBEGIN = 1
+            DO 39  JBLK = 1, IWORK( IINDBL+M-1 )
+               IEND = IWORK( IINSPL+JBLK-1 )
+               IN = IEND - IBEGIN + 1
+               WEND = WBEGIN - 1
+*              check if any eigenvalues have to be refined in this block
  36         CONTINUE
-            IF( WEND.LT.M ) THEN
-               IF( IWORK( IINDBL+WEND ).EQ.JBLK ) THEN
-                  WEND = WEND + 1
-                  GO TO 36
+               IF( WEND.LT.M ) THEN
+                  IF( IWORK( IINDBL+WEND ).EQ.JBLK ) THEN
+                     WEND = WEND + 1
+                     GO TO 36
+                  END IF
+               END IF
+               IF( WEND.LT.WBEGIN ) THEN
+                  IBEGIN = IEND + 1
+                  GO TO 39
                END IF
-            END IF
-            IF( WEND.LT.WBEGIN ) THEN
-               IBEGIN = IEND + 1
-               GO TO 39
-            END IF
 
-            OFFSET = IWORK(IINDW+WBEGIN-1)-1
-            IFIRST = IWORK(IINDW+WBEGIN-1)
-            ILAST = IWORK(IINDW+WEND-1)
-            RTOL2 = FOUR * EPS
-            CALL SLARRJ( IN,
+               OFFSET = IWORK(IINDW+WBEGIN-1)-1
+               IFIRST = IWORK(IINDW+WBEGIN-1)
+               ILAST = IWORK(IINDW+WEND-1)
+               RTOL2 = FOUR * EPS
+               CALL SLARRJ( IN,
      $                   WORK(INDD+IBEGIN-1), WORK(INDE2+IBEGIN-1),
      $                   IFIRST, ILAST, RTOL2, OFFSET, W(WBEGIN),
      $                   WORK( INDERR+WBEGIN-1 ),
      $                   WORK( INDWRK ), IWORK( IINDWK ), PIVMIN,
      $                   TNRM, IINFO )
-            IBEGIN = IEND + 1
-            WBEGIN = WEND + 1
+               IBEGIN = IEND + 1
+               WBEGIN = WEND + 1
  39      CONTINUE
-      ENDIF
+         ENDIF
 *
-*     If matrix was scaled, then rescale eigenvalues appropriately.
+*        If matrix was scaled, then rescale eigenvalues appropriately.
 *
-      IF( SCALE.NE.ONE ) THEN
-         CALL SSCAL( M, ONE / SCALE, W, 1 )
+         IF( SCALE.NE.ONE ) THEN
+            CALL SSCAL( M, ONE / SCALE, W, 1 )
+         END IF
       END IF
 *
 *     If eigenvalues are not in increasing order, then sort them,
 *     possibly along with eigenvectors.
 *
-      IF( NSPLIT.GT.1 ) THEN
+      IF( NSPLIT.GT.1 .OR. N.EQ.2 ) THEN
          IF( .NOT. WANTZ ) THEN
             CALL SLASRT( 'I', M, W, IINFO )
             IF( IINFO.NE.0 ) THEN
diff --git a/SRC/dstemr.f b/SRC/dstemr.f
index 24e0a4b4..58d4299b 100644
--- a/SRC/dstemr.f
+++ b/SRC/dstemr.f
@@ -543,184 +543,187 @@
                END IF
             ENDIF
          ENDIF
-         RETURN
-      END IF
+
+      ELSE
 
 *     Continue with general N
 
-      INDGRS = 1
-      INDERR = 2*N + 1
-      INDGP = 3*N + 1
-      INDD = 4*N + 1
-      INDE2 = 5*N + 1
-      INDWRK = 6*N + 1
-*
-      IINSPL = 1
-      IINDBL = N + 1
-      IINDW = 2*N + 1
-      IINDWK = 3*N + 1
-*
-*     Scale matrix to allowable range, if necessary.
-*     The allowable range is related to the PIVMIN parameter; see the
-*     comments in DLARRD.  The preference for scaling small values
-*     up is heuristic; we expect users' matrices not to be close to the
-*     RMAX threshold.
-*
-      SCALE = ONE
-      TNRM = DLANST( 'M', N, D, E )
-      IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
-         SCALE = RMIN / TNRM
-      ELSE IF( TNRM.GT.RMAX ) THEN
-         SCALE = RMAX / TNRM
-      END IF
-      IF( SCALE.NE.ONE ) THEN
-         CALL DSCAL( N, SCALE, D, 1 )
-         CALL DSCAL( N-1, SCALE, E, 1 )
-         TNRM = TNRM*SCALE
-         IF( VALEIG ) THEN
-*           If eigenvalues in interval have to be found,
-*           scale (WL, WU] accordingly
-            WL = WL*SCALE
-            WU = WU*SCALE
-         ENDIF
-      END IF
+         INDGRS = 1
+         INDERR = 2*N + 1
+         INDGP = 3*N + 1
+         INDD = 4*N + 1
+         INDE2 = 5*N + 1
+         INDWRK = 6*N + 1
+*
+         IINSPL = 1
+         IINDBL = N + 1
+         IINDW = 2*N + 1
+         IINDWK = 3*N + 1
+*
+*        Scale matrix to allowable range, if necessary.
+*        The allowable range is related to the PIVMIN parameter; see the
+*        comments in DLARRD.  The preference for scaling small values
+*        up is heuristic; we expect users' matrices not to be close to the
+*        RMAX threshold.
+*
+         SCALE = ONE
+         TNRM = DLANST( 'M', N, D, E )
+         IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
+            SCALE = RMIN / TNRM
+         ELSE IF( TNRM.GT.RMAX ) THEN
+            SCALE = RMAX / TNRM
+         END IF
+         IF( SCALE.NE.ONE ) THEN
+            CALL DSCAL( N, SCALE, D, 1 )
+            CALL DSCAL( N-1, SCALE, E, 1 )
+            TNRM = TNRM*SCALE
+            IF( VALEIG ) THEN
+*              If eigenvalues in interval have to be found,
+*              scale (WL, WU] accordingly
+               WL = WL*SCALE
+               WU = WU*SCALE
+            ENDIF
+         END IF
 *
-*     Compute the desired eigenvalues of the tridiagonal after splitting
-*     into smaller subblocks if the corresponding off-diagonal elements
-*     are small
-*     THRESH is the splitting parameter for DLARRE
-*     A negative THRESH forces the old splitting criterion based on the
-*     size of the off-diagonal. A positive THRESH switches to splitting
-*     which preserves relative accuracy.
-*
-      IF( TRYRAC ) THEN
-*        Test whether the matrix warrants the more expensive relative approach.
-         CALL DLARRR( N, D, E, IINFO )
-      ELSE
-*        The user does not care about relative accurately eigenvalues
-         IINFO = -1
-      ENDIF
-*     Set the splitting criterion
-      IF (IINFO.EQ.0) THEN
-         THRESH = EPS
-      ELSE
-         THRESH = -EPS
-*        relative accuracy is desired but T does not guarantee it
-         TRYRAC = .FALSE.
-      ENDIF
+*        Compute the desired eigenvalues of the tridiagonal after splitting
+*        into smaller subblocks if the corresponding off-diagonal elements
+*        are small
+*        THRESH is the splitting parameter for DLARRE
+*        A negative THRESH forces the old splitting criterion based on the
+*        size of the off-diagonal. A positive THRESH switches to splitting
+*        which preserves relative accuracy.
+*
+         IF( TRYRAC ) THEN
+*           Test whether the matrix warrants the more expensive relative approach.
+            CALL DLARRR( N, D, E, IINFO )
+         ELSE
+*           The user does not care about relative accurately eigenvalues
+            IINFO = -1
+         ENDIF
+*        Set the splitting criterion
+         IF (IINFO.EQ.0) THEN
+            THRESH = EPS
+         ELSE
+            THRESH = -EPS
+*           relative accuracy is desired but T does not guarantee it
+            TRYRAC = .FALSE.
+         ENDIF
 *
-      IF( TRYRAC ) THEN
-*        Copy original diagonal, needed to guarantee relative accuracy
-         CALL DCOPY(N,D,1,WORK(INDD),1)
-      ENDIF
-*     Store the squares of the offdiagonal values of T
-      DO 5 J = 1, N-1
-         WORK( INDE2+J-1 ) = E(J)**2
- 5    CONTINUE
+         IF( TRYRAC ) THEN
+*           Copy original diagonal, needed to guarantee relative accuracy
+            CALL DCOPY(N,D,1,WORK(INDD),1)
+         ENDIF
+*        Store the squares of the offdiagonal values of T
+         DO 5 J = 1, N-1
+            WORK( INDE2+J-1 ) = E(J)**2
+ 5       CONTINUE
 
-*     Set the tolerance parameters for bisection
-      IF( .NOT.WANTZ ) THEN
-*        DLARRE computes the eigenvalues to full precision.
-         RTOL1 = FOUR * EPS
-         RTOL2 = FOUR * EPS
-      ELSE
-*        DLARRE computes the eigenvalues to less than full precision.
-*        DLARRV will refine the eigenvalue approximations, and we can
-*        need less accurate initial bisection in DLARRE.
-*        Note: these settings do only affect the subset case and DLARRE
-         RTOL1 = SQRT(EPS)
-         RTOL2 = MAX( SQRT(EPS)*5.0D-3, FOUR * EPS )
-      ENDIF
-      CALL DLARRE( RANGE, N, WL, WU, IIL, IIU, D, E,
+*        Set the tolerance parameters for bisection
+         IF( .NOT.WANTZ ) THEN
+*           DLARRE computes the eigenvalues to full precision.
+            RTOL1 = FOUR * EPS
+            RTOL2 = FOUR * EPS
+         ELSE
+*           DLARRE computes the eigenvalues to less than full precision.
+*           DLARRV will refine the eigenvalue approximations, and we can
+*           need less accurate initial bisection in DLARRE.
+*           Note: these settings do only affect the subset case and DLARRE
+            RTOL1 = SQRT(EPS)
+            RTOL2 = MAX( SQRT(EPS)*5.0D-3, FOUR * EPS )
+         ENDIF
+         CALL DLARRE( RANGE, N, WL, WU, IIL, IIU, D, E,
      $             WORK(INDE2), RTOL1, RTOL2, THRESH, NSPLIT,
      $             IWORK( IINSPL ), M, W, WORK( INDERR ),
      $             WORK( INDGP ), IWORK( IINDBL ),
      $             IWORK( IINDW ), WORK( INDGRS ), PIVMIN,
      $             WORK( INDWRK ), IWORK( IINDWK ), IINFO )
-      IF( IINFO.NE.0 ) THEN
-         INFO = 10 + ABS( IINFO )
-         RETURN
-      END IF
-*     Note that if RANGE .NE. 'V', DLARRE computes bounds on the desired
-*     part of the spectrum. All desired eigenvalues are contained in
-*     (WL,WU]
+         IF( IINFO.NE.0 ) THEN
+            INFO = 10 + ABS( IINFO )
+            RETURN
+         END IF
+*        Note that if RANGE .NE. 'V', DLARRE computes bounds on the desired
+*        part of the spectrum. All desired eigenvalues are contained in
+*        (WL,WU]
 
 
-      IF( WANTZ ) THEN
+         IF( WANTZ ) THEN
 *
-*        Compute the desired eigenvectors corresponding to the computed
-*        eigenvalues
+*           Compute the desired eigenvectors corresponding to the computed
+*           eigenvalues
 *
-         CALL DLARRV( N, WL, WU, D, E,
+            CALL DLARRV( N, WL, WU, D, E,
      $                PIVMIN, IWORK( IINSPL ), M,
      $                1, M, MINRGP, RTOL1, RTOL2,
      $                W, WORK( INDERR ), WORK( INDGP ), IWORK( IINDBL ),
      $                IWORK( IINDW ), WORK( INDGRS ), Z, LDZ,
      $                ISUPPZ, WORK( INDWRK ), IWORK( IINDWK ), IINFO )
-         IF( IINFO.NE.0 ) THEN
-            INFO = 20 + ABS( IINFO )
-            RETURN
+            IF( IINFO.NE.0 ) THEN
+               INFO = 20 + ABS( IINFO )
+               RETURN
+            END IF
+         ELSE
+*           DLARRE computes eigenvalues of the (shifted) root representation
+*           DLARRV returns the eigenvalues of the unshifted matrix.
+*           However, if the eigenvectors are not desired by the user, we need
+*           to apply the corresponding shifts from DLARRE to obtain the
+*           eigenvalues of the original matrix.
+            DO 20 J = 1, M
+               ITMP = IWORK( IINDBL+J-1 )
+               W( J ) = W( J ) + E( IWORK( IINSPL+ITMP-1 ) )
+ 20         CONTINUE
          END IF
-      ELSE
-*        DLARRE computes eigenvalues of the (shifted) root representation
-*        DLARRV returns the eigenvalues of the unshifted matrix.
-*        However, if the eigenvectors are not desired by the user, we need
-*        to apply the corresponding shifts from DLARRE to obtain the
-*        eigenvalues of the original matrix.
-         DO 20 J = 1, M
-            ITMP = IWORK( IINDBL+J-1 )
-            W( J ) = W( J ) + E( IWORK( IINSPL+ITMP-1 ) )
- 20      CONTINUE
-      END IF
 *
 
-      IF ( TRYRAC ) THEN
-*        Refine computed eigenvalues so that they are relatively accurate
-*        with respect to the original matrix T.
-         IBEGIN = 1
-         WBEGIN = 1
-         DO 39  JBLK = 1, IWORK( IINDBL+M-1 )
-            IEND = IWORK( IINSPL+JBLK-1 )
-            IN = IEND - IBEGIN + 1
-            WEND = WBEGIN - 1
-*           check if any eigenvalues have to be refined in this block
- 36         CONTINUE
-            IF( WEND.LT.M ) THEN
-               IF( IWORK( IINDBL+WEND ).EQ.JBLK ) THEN
-                  WEND = WEND + 1
-                  GO TO 36
+         IF ( TRYRAC ) THEN
+*           Refine computed eigenvalues so that they are relatively accurate
+*           with respect to the original matrix T.
+            IBEGIN = 1
+            WBEGIN = 1
+            DO 39  JBLK = 1, IWORK( IINDBL+M-1 )
+               IEND = IWORK( IINSPL+JBLK-1 )
+               IN = IEND - IBEGIN + 1
+               WEND = WBEGIN - 1
+*              check if any eigenvalues have to be refined in this block
+ 36            CONTINUE
+               IF( WEND.LT.M ) THEN
+                  IF( IWORK( IINDBL+WEND ).EQ.JBLK ) THEN
+                     WEND = WEND + 1
+                     GO TO 36
+                  END IF
+               END IF
+               IF( WEND.LT.WBEGIN ) THEN
+                  IBEGIN = IEND + 1
+                  GO TO 39
                END IF
-            END IF
-            IF( WEND.LT.WBEGIN ) THEN
-               IBEGIN = IEND + 1
-               GO TO 39
-            END IF
 
-            OFFSET = IWORK(IINDW+WBEGIN-1)-1
-            IFIRST = IWORK(IINDW+WBEGIN-1)
-            ILAST = IWORK(IINDW+WEND-1)
-            RTOL2 = FOUR * EPS
-            CALL DLARRJ( IN,
+               OFFSET = IWORK(IINDW+WBEGIN-1)-1
+               IFIRST = IWORK(IINDW+WBEGIN-1)
+               ILAST = IWORK(IINDW+WEND-1)
+               RTOL2 = FOUR * EPS
+               CALL DLARRJ( IN,
      $                   WORK(INDD+IBEGIN-1), WORK(INDE2+IBEGIN-1),
      $                   IFIRST, ILAST, RTOL2, OFFSET, W(WBEGIN),
      $                   WORK( INDERR+WBEGIN-1 ),
      $                   WORK( INDWRK ), IWORK( IINDWK ), PIVMIN,
      $                   TNRM, IINFO )
-            IBEGIN = IEND + 1
-            WBEGIN = WEND + 1
- 39      CONTINUE
-      ENDIF
+               IBEGIN = IEND + 1
+               WBEGIN = WEND + 1
+ 39         CONTINUE
+         ENDIF
 *
-*     If matrix was scaled, then rescale eigenvalues appropriately.
+*        If matrix was scaled, then rescale eigenvalues appropriately.
 *
-      IF( SCALE.NE.ONE ) THEN
-         CALL DSCAL( M, ONE / SCALE, W, 1 )
+         IF( SCALE.NE.ONE ) THEN
+            CALL DSCAL( M, ONE / SCALE, W, 1 )
+         END IF
+  
       END IF
+    
 *
 *     If eigenvalues are not in increasing order, then sort them,
 *     possibly along with eigenvectors.
 *
-      IF( NSPLIT.GT.1 ) THEN
+      IF( NSPLIT.GT.1 .OR. N.EQ.2 ) THEN
          IF( .NOT. WANTZ ) THEN
             CALL DLASRT( 'I', M, W, IINFO )
             IF( IINFO.NE.0 ) THEN
diff --git a/SRC/sstemr.f b/SRC/sstemr.f
index 7b1ab2eb..a06f6b44 100644
--- a/SRC/sstemr.f
+++ b/SRC/sstemr.f
@@ -541,184 +541,184 @@
                END IF
             ENDIF
          ENDIF
-         RETURN
-      END IF
+      ELSE
 
 *     Continue with general N
 
-      INDGRS = 1
-      INDERR = 2*N + 1
-      INDGP = 3*N + 1
-      INDD = 4*N + 1
-      INDE2 = 5*N + 1
-      INDWRK = 6*N + 1
-*
-      IINSPL = 1
-      IINDBL = N + 1
-      IINDW = 2*N + 1
-      IINDWK = 3*N + 1
-*
-*     Scale matrix to allowable range, if necessary.
-*     The allowable range is related to the PIVMIN parameter; see the
-*     comments in SLARRD.  The preference for scaling small values
-*     up is heuristic; we expect users' matrices not to be close to the
-*     RMAX threshold.
-*
-      SCALE = ONE
-      TNRM = SLANST( 'M', N, D, E )
-      IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
-         SCALE = RMIN / TNRM
-      ELSE IF( TNRM.GT.RMAX ) THEN
-         SCALE = RMAX / TNRM
-      END IF
-      IF( SCALE.NE.ONE ) THEN
-         CALL SSCAL( N, SCALE, D, 1 )
-         CALL SSCAL( N-1, SCALE, E, 1 )
-         TNRM = TNRM*SCALE
-         IF( VALEIG ) THEN
-*           If eigenvalues in interval have to be found,
-*           scale (WL, WU] accordingly
-            WL = WL*SCALE
-            WU = WU*SCALE
-         ENDIF
-      END IF
+         INDGRS = 1
+         INDERR = 2*N + 1
+         INDGP = 3*N + 1
+         INDD = 4*N + 1
+         INDE2 = 5*N + 1
+         INDWRK = 6*N + 1
+*
+         IINSPL = 1
+         IINDBL = N + 1
+         IINDW = 2*N + 1
+         IINDWK = 3*N + 1
+*
+*        Scale matrix to allowable range, if necessary.
+*        The allowable range is related to the PIVMIN parameter; see the
+*        comments in SLARRD.  The preference for scaling small values
+*        up is heuristic; we expect users' matrices not to be close to the
+*        RMAX threshold.
+*
+         SCALE = ONE
+         TNRM = SLANST( 'M', N, D, E )
+         IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
+            SCALE = RMIN / TNRM
+         ELSE IF( TNRM.GT.RMAX ) THEN
+            SCALE = RMAX / TNRM
+         END IF
+         IF( SCALE.NE.ONE ) THEN
+            CALL SSCAL( N, SCALE, D, 1 )
+            CALL SSCAL( N-1, SCALE, E, 1 )
+            TNRM = TNRM*SCALE
+            IF( VALEIG ) THEN
+*              If eigenvalues in interval have to be found,
+*              scale (WL, WU] accordingly
+               WL = WL*SCALE
+               WU = WU*SCALE
+            ENDIF
+         END IF
 *
-*     Compute the desired eigenvalues of the tridiagonal after splitting
-*     into smaller subblocks if the corresponding off-diagonal elements
-*     are small
-*     THRESH is the splitting parameter for SLARRE
-*     A negative THRESH forces the old splitting criterion based on the
-*     size of the off-diagonal. A positive THRESH switches to splitting
-*     which preserves relative accuracy.
-*
-      IF( TRYRAC ) THEN
-*        Test whether the matrix warrants the more expensive relative approach.
-         CALL SLARRR( N, D, E, IINFO )
-      ELSE
-*        The user does not care about relative accurately eigenvalues
-         IINFO = -1
-      ENDIF
-*     Set the splitting criterion
-      IF (IINFO.EQ.0) THEN
-         THRESH = EPS
-      ELSE
-         THRESH = -EPS
-*        relative accuracy is desired but T does not guarantee it
-         TRYRAC = .FALSE.
-      ENDIF
+*        Compute the desired eigenvalues of the tridiagonal after splitting
+*        into smaller subblocks if the corresponding off-diagonal elements
+*        are small
+*        THRESH is the splitting parameter for SLARRE
+*        A negative THRESH forces the old splitting criterion based on the
+*        size of the off-diagonal. A positive THRESH switches to splitting
+*        which preserves relative accuracy.
+*
+         IF( TRYRAC ) THEN
+*           Test whether the matrix warrants the more expensive relative approach.
+            CALL SLARRR( N, D, E, IINFO )
+         ELSE
+*           The user does not care about relative accurately eigenvalues
+            IINFO = -1
+         ENDIF
+*        Set the splitting criterion
+         IF (IINFO.EQ.0) THEN
+            THRESH = EPS
+         ELSE
+            THRESH = -EPS
+*           relative accuracy is desired but T does not guarantee it
+            TRYRAC = .FALSE.
+         ENDIF
 *
-      IF( TRYRAC ) THEN
-*        Copy original diagonal, needed to guarantee relative accuracy
-         CALL SCOPY(N,D,1,WORK(INDD),1)
-      ENDIF
-*     Store the squares of the offdiagonal values of T
-      DO 5 J = 1, N-1
-         WORK( INDE2+J-1 ) = E(J)**2
+         IF( TRYRAC ) THEN
+*           Copy original diagonal, needed to guarantee relative accuracy
+            CALL SCOPY(N,D,1,WORK(INDD),1)
+         ENDIF
+*        Store the squares of the offdiagonal values of T
+         DO 5 J = 1, N-1
+            WORK( INDE2+J-1 ) = E(J)**2
  5    CONTINUE
 
-*     Set the tolerance parameters for bisection
-      IF( .NOT.WANTZ ) THEN
-*        SLARRE computes the eigenvalues to full precision.
-         RTOL1 = FOUR * EPS
-         RTOL2 = FOUR * EPS
-      ELSE
-*        SLARRE computes the eigenvalues to less than full precision.
-*        SLARRV will refine the eigenvalue approximations, and we can
-*        need less accurate initial bisection in SLARRE.
-*        Note: these settings do only affect the subset case and SLARRE
-         RTOL1 = MAX( SQRT(EPS)*5.0E-2, FOUR * EPS )
-         RTOL2 = MAX( SQRT(EPS)*5.0E-3, FOUR * EPS )
-      ENDIF
-      CALL SLARRE( RANGE, N, WL, WU, IIL, IIU, D, E,
+*        Set the tolerance parameters for bisection
+         IF( .NOT.WANTZ ) THEN
+*           SLARRE computes the eigenvalues to full precision.
+            RTOL1 = FOUR * EPS
+            RTOL2 = FOUR * EPS
+         ELSE
+*           SLARRE computes the eigenvalues to less than full precision.
+*           SLARRV will refine the eigenvalue approximations, and we can
+*           need less accurate initial bisection in SLARRE.
+*           Note: these settings do only affect the subset case and SLARRE
+            RTOL1 = MAX( SQRT(EPS)*5.0E-2, FOUR * EPS )
+            RTOL2 = MAX( SQRT(EPS)*5.0E-3, FOUR * EPS )
+         ENDIF
+         CALL SLARRE( RANGE, N, WL, WU, IIL, IIU, D, E,
      $             WORK(INDE2), RTOL1, RTOL2, THRESH, NSPLIT,
      $             IWORK( IINSPL ), M, W, WORK( INDERR ),
      $             WORK( INDGP ), IWORK( IINDBL ),
      $             IWORK( IINDW ), WORK( INDGRS ), PIVMIN,
      $             WORK( INDWRK ), IWORK( IINDWK ), IINFO )
-      IF( IINFO.NE.0 ) THEN
-         INFO = 10 + ABS( IINFO )
-         RETURN
-      END IF
-*     Note that if RANGE .NE. 'V', SLARRE computes bounds on the desired
-*     part of the spectrum. All desired eigenvalues are contained in
-*     (WL,WU]
+         IF( IINFO.NE.0 ) THEN
+            INFO = 10 + ABS( IINFO )
+            RETURN
+         END IF
+*        Note that if RANGE .NE. 'V', SLARRE computes bounds on the desired
+*        part of the spectrum. All desired eigenvalues are contained in
+*        (WL,WU]
 
 
-      IF( WANTZ ) THEN
+         IF( WANTZ ) THEN
 *
-*        Compute the desired eigenvectors corresponding to the computed
-*        eigenvalues
+*           Compute the desired eigenvectors corresponding to the computed
+*           eigenvalues
 *
-         CALL SLARRV( N, WL, WU, D, E,
+            CALL SLARRV( N, WL, WU, D, E,
      $                PIVMIN, IWORK( IINSPL ), M,
      $                1, M, MINRGP, RTOL1, RTOL2,
      $                W, WORK( INDERR ), WORK( INDGP ), IWORK( IINDBL ),
      $                IWORK( IINDW ), WORK( INDGRS ), Z, LDZ,
      $                ISUPPZ, WORK( INDWRK ), IWORK( IINDWK ), IINFO )
-         IF( IINFO.NE.0 ) THEN
-            INFO = 20 + ABS( IINFO )
-            RETURN
-         END IF
-      ELSE
-*        SLARRE computes eigenvalues of the (shifted) root representation
-*        SLARRV returns the eigenvalues of the unshifted matrix.
-*        However, if the eigenvectors are not desired by the user, we need
-*        to apply the corresponding shifts from SLARRE to obtain the
-*        eigenvalues of the original matrix.
-         DO 20 J = 1, M
-            ITMP = IWORK( IINDBL+J-1 )
-            W( J ) = W( J ) + E( IWORK( IINSPL+ITMP-1 ) )
+            IF( IINFO.NE.0 ) THEN
+               INFO = 20 + ABS( IINFO )
+               RETURN
+            END IF
+         ELSE
+*           SLARRE computes eigenvalues of the (shifted) root representation
+*           SLARRV returns the eigenvalues of the unshifted matrix.
+*           However, if the eigenvectors are not desired by the user, we need
+*           to apply the corresponding shifts from SLARRE to obtain the
+*           eigenvalues of the original matrix.
+            DO 20 J = 1, M
+               ITMP = IWORK( IINDBL+J-1 )
+               W( J ) = W( J ) + E( IWORK( IINSPL+ITMP-1 ) )
  20      CONTINUE
-      END IF
+         END IF
 *
 
-      IF ( TRYRAC ) THEN
-*        Refine computed eigenvalues so that they are relatively accurate
-*        with respect to the original matrix T.
-         IBEGIN = 1
-         WBEGIN = 1
-         DO 39  JBLK = 1, IWORK( IINDBL+M-1 )
-            IEND = IWORK( IINSPL+JBLK-1 )
-            IN = IEND - IBEGIN + 1
-            WEND = WBEGIN - 1
-*           check if any eigenvalues have to be refined in this block
+         IF ( TRYRAC ) THEN
+*           Refine computed eigenvalues so that they are relatively accurate
+*           with respect to the original matrix T.
+            IBEGIN = 1
+            WBEGIN = 1
+            DO 39  JBLK = 1, IWORK( IINDBL+M-1 )
+               IEND = IWORK( IINSPL+JBLK-1 )
+               IN = IEND - IBEGIN + 1
+               WEND = WBEGIN - 1
+*              check if any eigenvalues have to be refined in this block
  36         CONTINUE
-            IF( WEND.LT.M ) THEN
-               IF( IWORK( IINDBL+WEND ).EQ.JBLK ) THEN
-                  WEND = WEND + 1
-                  GO TO 36
+               IF( WEND.LT.M ) THEN
+                  IF( IWORK( IINDBL+WEND ).EQ.JBLK ) THEN
+                     WEND = WEND + 1
+                     GO TO 36
+                  END IF
+               END IF
+               IF( WEND.LT.WBEGIN ) THEN
+                  IBEGIN = IEND + 1
+                  GO TO 39
                END IF
-            END IF
-            IF( WEND.LT.WBEGIN ) THEN
-               IBEGIN = IEND + 1
-               GO TO 39
-            END IF
 
-            OFFSET = IWORK(IINDW+WBEGIN-1)-1
-            IFIRST = IWORK(IINDW+WBEGIN-1)
-            ILAST = IWORK(IINDW+WEND-1)
-            RTOL2 = FOUR * EPS
-            CALL SLARRJ( IN,
+               OFFSET = IWORK(IINDW+WBEGIN-1)-1
+               IFIRST = IWORK(IINDW+WBEGIN-1)
+               ILAST = IWORK(IINDW+WEND-1)
+               RTOL2 = FOUR * EPS
+               CALL SLARRJ( IN,
      $                   WORK(INDD+IBEGIN-1), WORK(INDE2+IBEGIN-1),
      $                   IFIRST, ILAST, RTOL2, OFFSET, W(WBEGIN),
      $                   WORK( INDERR+WBEGIN-1 ),
      $                   WORK( INDWRK ), IWORK( IINDWK ), PIVMIN,
      $                   TNRM, IINFO )
-            IBEGIN = IEND + 1
-            WBEGIN = WEND + 1
+               IBEGIN = IEND + 1
+               WBEGIN = WEND + 1
  39      CONTINUE
-      ENDIF
+         ENDIF
 *
-*     If matrix was scaled, then rescale eigenvalues appropriately.
+*        If matrix was scaled, then rescale eigenvalues appropriately.
 *
-      IF( SCALE.NE.ONE ) THEN
-         CALL SSCAL( M, ONE / SCALE, W, 1 )
+         IF( SCALE.NE.ONE ) THEN
+            CALL SSCAL( M, ONE / SCALE, W, 1 )
+         END IF
       END IF
 *
 *     If eigenvalues are not in increasing order, then sort them,
 *     possibly along with eigenvectors.
 *
-      IF( NSPLIT.GT.1 ) THEN
+      IF( NSPLIT.GT.1 .OR. N.EQ.2 ) THEN
          IF( .NOT. WANTZ ) THEN
             CALL SLASRT( 'I', M, W, IINFO )
             IF( IINFO.NE.0 ) THEN
diff --git a/SRC/zstemr.f b/SRC/zstemr.f
index 6d0833fa..8cee9316 100644
--- a/SRC/zstemr.f
+++ b/SRC/zstemr.f
@@ -560,184 +560,184 @@
                END IF
             ENDIF
          ENDIF
-         RETURN
-      END IF
+      ELSE
 
-*     Continue with general N
+*        Continue with general N
 
-      INDGRS = 1
-      INDERR = 2*N + 1
-      INDGP = 3*N + 1
-      INDD = 4*N + 1
-      INDE2 = 5*N + 1
-      INDWRK = 6*N + 1
-*
-      IINSPL = 1
-      IINDBL = N + 1
-      IINDW = 2*N + 1
-      IINDWK = 3*N + 1
-*
-*     Scale matrix to allowable range, if necessary.
-*     The allowable range is related to the PIVMIN parameter; see the
-*     comments in DLARRD.  The preference for scaling small values
-*     up is heuristic; we expect users' matrices not to be close to the
-*     RMAX threshold.
-*
-      SCALE = ONE
-      TNRM = DLANST( 'M', N, D, E )
-      IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
-         SCALE = RMIN / TNRM
-      ELSE IF( TNRM.GT.RMAX ) THEN
-         SCALE = RMAX / TNRM
-      END IF
-      IF( SCALE.NE.ONE ) THEN
-         CALL DSCAL( N, SCALE, D, 1 )
-         CALL DSCAL( N-1, SCALE, E, 1 )
-         TNRM = TNRM*SCALE
-         IF( VALEIG ) THEN
-*           If eigenvalues in interval have to be found,
-*           scale (WL, WU] accordingly
-            WL = WL*SCALE
-            WU = WU*SCALE
-         ENDIF
-      END IF
+         INDGRS = 1
+         INDERR = 2*N + 1
+         INDGP = 3*N + 1
+         INDD = 4*N + 1
+         INDE2 = 5*N + 1
+         INDWRK = 6*N + 1
+*
+         IINSPL = 1
+         IINDBL = N + 1
+         IINDW = 2*N + 1
+         IINDWK = 3*N + 1
+*
+*        Scale matrix to allowable range, if necessary.
+*        The allowable range is related to the PIVMIN parameter; see the
+*        comments in DLARRD.  The preference for scaling small values
+*        up is heuristic; we expect users' matrices not to be close to the
+*        RMAX threshold.
+*
+         SCALE = ONE
+         TNRM = DLANST( 'M', N, D, E )
+         IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
+            SCALE = RMIN / TNRM
+         ELSE IF( TNRM.GT.RMAX ) THEN
+            SCALE = RMAX / TNRM
+         END IF
+         IF( SCALE.NE.ONE ) THEN
+            CALL DSCAL( N, SCALE, D, 1 )
+            CALL DSCAL( N-1, SCALE, E, 1 )
+            TNRM = TNRM*SCALE
+            IF( VALEIG ) THEN
+*              If eigenvalues in interval have to be found,
+*              scale (WL, WU] accordingly
+               WL = WL*SCALE
+               WU = WU*SCALE
+            ENDIF
+         END IF
 *
-*     Compute the desired eigenvalues of the tridiagonal after splitting
-*     into smaller subblocks if the corresponding off-diagonal elements
-*     are small
-*     THRESH is the splitting parameter for DLARRE
-*     A negative THRESH forces the old splitting criterion based on the
-*     size of the off-diagonal. A positive THRESH switches to splitting
-*     which preserves relative accuracy.
-*
-      IF( TRYRAC ) THEN
-*        Test whether the matrix warrants the more expensive relative approach.
-         CALL DLARRR( N, D, E, IINFO )
-      ELSE
-*        The user does not care about relative accurately eigenvalues
-         IINFO = -1
-      ENDIF
-*     Set the splitting criterion
-      IF (IINFO.EQ.0) THEN
-         THRESH = EPS
-      ELSE
-         THRESH = -EPS
-*        relative accuracy is desired but T does not guarantee it
-         TRYRAC = .FALSE.
-      ENDIF
+*        Compute the desired eigenvalues of the tridiagonal after splitting
+*        into smaller subblocks if the corresponding off-diagonal elements
+*        are small
+*        THRESH is the splitting parameter for DLARRE
+*        A negative THRESH forces the old splitting criterion based on the
+*        size of the off-diagonal. A positive THRESH switches to splitting
+*        which preserves relative accuracy.
+*
+         IF( TRYRAC ) THEN
+*           Test whether the matrix warrants the more expensive relative approach.
+            CALL DLARRR( N, D, E, IINFO )
+         ELSE
+*           The user does not care about relative accurately eigenvalues
+            IINFO = -1
+         ENDIF
+*        Set the splitting criterion
+         IF (IINFO.EQ.0) THEN
+            THRESH = EPS
+         ELSE
+            THRESH = -EPS
+*           relative accuracy is desired but T does not guarantee it
+            TRYRAC = .FALSE.
+         ENDIF
 *
-      IF( TRYRAC ) THEN
-*        Copy original diagonal, needed to guarantee relative accuracy
-         CALL DCOPY(N,D,1,WORK(INDD),1)
-      ENDIF
-*     Store the squares of the offdiagonal values of T
-      DO 5 J = 1, N-1
-         WORK( INDE2+J-1 ) = E(J)**2
+         IF( TRYRAC ) THEN
+*           Copy original diagonal, needed to guarantee relative accuracy
+            CALL DCOPY(N,D,1,WORK(INDD),1)
+         ENDIF
+*        Store the squares of the offdiagonal values of T
+         DO 5 J = 1, N-1
+            WORK( INDE2+J-1 ) = E(J)**2
  5    CONTINUE
 
-*     Set the tolerance parameters for bisection
-      IF( .NOT.WANTZ ) THEN
-*        DLARRE computes the eigenvalues to full precision.
-         RTOL1 = FOUR * EPS
-         RTOL2 = FOUR * EPS
-      ELSE
-*        DLARRE computes the eigenvalues to less than full precision.
-*        ZLARRV will refine the eigenvalue approximations, and we only
-*        need less accurate initial bisection in DLARRE.
-*        Note: these settings do only affect the subset case and DLARRE
-         RTOL1 = SQRT(EPS)
-         RTOL2 = MAX( SQRT(EPS)*5.0D-3, FOUR * EPS )
-      ENDIF
-      CALL DLARRE( RANGE, N, WL, WU, IIL, IIU, D, E,
+*        Set the tolerance parameters for bisection
+         IF( .NOT.WANTZ ) THEN
+*           DLARRE computes the eigenvalues to full precision.
+            RTOL1 = FOUR * EPS
+            RTOL2 = FOUR * EPS
+         ELSE
+*           DLARRE computes the eigenvalues to less than full precision.
+*           ZLARRV will refine the eigenvalue approximations, and we only
+*           need less accurate initial bisection in DLARRE.
+*           Note: these settings do only affect the subset case and DLARRE
+            RTOL1 = SQRT(EPS)
+            RTOL2 = MAX( SQRT(EPS)*5.0D-3, FOUR * EPS )
+         ENDIF
+         CALL DLARRE( RANGE, N, WL, WU, IIL, IIU, D, E,
      $             WORK(INDE2), RTOL1, RTOL2, THRESH, NSPLIT,
      $             IWORK( IINSPL ), M, W, WORK( INDERR ),
      $             WORK( INDGP ), IWORK( IINDBL ),
      $             IWORK( IINDW ), WORK( INDGRS ), PIVMIN,
      $             WORK( INDWRK ), IWORK( IINDWK ), IINFO )
-      IF( IINFO.NE.0 ) THEN
-         INFO = 10 + ABS( IINFO )
-         RETURN
-      END IF
-*     Note that if RANGE .NE. 'V', DLARRE computes bounds on the desired
-*     part of the spectrum. All desired eigenvalues are contained in
-*     (WL,WU]
+         IF( IINFO.NE.0 ) THEN
+            INFO = 10 + ABS( IINFO )
+            RETURN
+         END IF
+*        Note that if RANGE .NE. 'V', DLARRE computes bounds on the desired
+*        part of the spectrum. All desired eigenvalues are contained in
+*        (WL,WU]
 
 
-      IF( WANTZ ) THEN
+         IF( WANTZ ) THEN
 *
-*        Compute the desired eigenvectors corresponding to the computed
-*        eigenvalues
+*           Compute the desired eigenvectors corresponding to the computed
+*           eigenvalues
 *
-         CALL ZLARRV( N, WL, WU, D, E,
+            CALL ZLARRV( N, WL, WU, D, E,
      $                PIVMIN, IWORK( IINSPL ), M,
      $                1, M, MINRGP, RTOL1, RTOL2,
      $                W, WORK( INDERR ), WORK( INDGP ), IWORK( IINDBL ),
      $                IWORK( IINDW ), WORK( INDGRS ), Z, LDZ,
      $                ISUPPZ, WORK( INDWRK ), IWORK( IINDWK ), IINFO )
-         IF( IINFO.NE.0 ) THEN
-            INFO = 20 + ABS( IINFO )
-            RETURN
-         END IF
-      ELSE
-*        DLARRE computes eigenvalues of the (shifted) root representation
-*        ZLARRV returns the eigenvalues of the unshifted matrix.
-*        However, if the eigenvectors are not desired by the user, we need
-*        to apply the corresponding shifts from DLARRE to obtain the
-*        eigenvalues of the original matrix.
-         DO 20 J = 1, M
-            ITMP = IWORK( IINDBL+J-1 )
-            W( J ) = W( J ) + E( IWORK( IINSPL+ITMP-1 ) )
+            IF( IINFO.NE.0 ) THEN
+               INFO = 20 + ABS( IINFO )
+               RETURN
+            END IF
+         ELSE
+*           DLARRE computes eigenvalues of the (shifted) root representation
+*           ZLARRV returns the eigenvalues of the unshifted matrix.
+*           However, if the eigenvectors are not desired by the user, we need
+*           to apply the corresponding shifts from DLARRE to obtain the
+*           eigenvalues of the original matrix.
+            DO 20 J = 1, M
+               ITMP = IWORK( IINDBL+J-1 )
+               W( J ) = W( J ) + E( IWORK( IINSPL+ITMP-1 ) )
  20      CONTINUE
-      END IF
+         END IF
 *
 
-      IF ( TRYRAC ) THEN
-*        Refine computed eigenvalues so that they are relatively accurate
-*        with respect to the original matrix T.
-         IBEGIN = 1
-         WBEGIN = 1
-         DO 39  JBLK = 1, IWORK( IINDBL+M-1 )
-            IEND = IWORK( IINSPL+JBLK-1 )
-            IN = IEND - IBEGIN + 1
-            WEND = WBEGIN - 1
-*           check if any eigenvalues have to be refined in this block
+         IF ( TRYRAC ) THEN
+*           Refine computed eigenvalues so that they are relatively accurate
+*           with respect to the original matrix T.
+            IBEGIN = 1
+            WBEGIN = 1
+            DO 39  JBLK = 1, IWORK( IINDBL+M-1 )
+               IEND = IWORK( IINSPL+JBLK-1 )
+               IN = IEND - IBEGIN + 1
+               WEND = WBEGIN - 1
+*              check if any eigenvalues have to be refined in this block
  36         CONTINUE
-            IF( WEND.LT.M ) THEN
-               IF( IWORK( IINDBL+WEND ).EQ.JBLK ) THEN
-                  WEND = WEND + 1
-                  GO TO 36
+               IF( WEND.LT.M ) THEN
+                  IF( IWORK( IINDBL+WEND ).EQ.JBLK ) THEN
+                     WEND = WEND + 1
+                     GO TO 36
+                  END IF
+               END IF
+               IF( WEND.LT.WBEGIN ) THEN
+                  IBEGIN = IEND + 1
+                  GO TO 39
                END IF
-            END IF
-            IF( WEND.LT.WBEGIN ) THEN
-               IBEGIN = IEND + 1
-               GO TO 39
-            END IF
 
-            OFFSET = IWORK(IINDW+WBEGIN-1)-1
-            IFIRST = IWORK(IINDW+WBEGIN-1)
-            ILAST = IWORK(IINDW+WEND-1)
-            RTOL2 = FOUR * EPS
-            CALL DLARRJ( IN,
+               OFFSET = IWORK(IINDW+WBEGIN-1)-1
+               IFIRST = IWORK(IINDW+WBEGIN-1)
+               ILAST = IWORK(IINDW+WEND-1)
+               RTOL2 = FOUR * EPS
+               CALL DLARRJ( IN,
      $                   WORK(INDD+IBEGIN-1), WORK(INDE2+IBEGIN-1),
      $                   IFIRST, ILAST, RTOL2, OFFSET, W(WBEGIN),
      $                   WORK( INDERR+WBEGIN-1 ),
      $                   WORK( INDWRK ), IWORK( IINDWK ), PIVMIN,
      $                   TNRM, IINFO )
-            IBEGIN = IEND + 1
-            WBEGIN = WEND + 1
+               IBEGIN = IEND + 1
+               WBEGIN = WEND + 1
  39      CONTINUE
-      ENDIF
+         ENDIF
 *
-*     If matrix was scaled, then rescale eigenvalues appropriately.
+*        If matrix was scaled, then rescale eigenvalues appropriately.
 *
-      IF( SCALE.NE.ONE ) THEN
-         CALL DSCAL( M, ONE / SCALE, W, 1 )
+         IF( SCALE.NE.ONE ) THEN
+            CALL DSCAL( M, ONE / SCALE, W, 1 )
+         END IF
       END IF
 *
 *     If eigenvalues are not in increasing order, then sort them,
 *     possibly along with eigenvectors.
 *
-      IF( NSPLIT.GT.1 ) THEN
+      IF( NSPLIT.GT.1 .OR. N.EQ.2 ) THEN
          IF( .NOT. WANTZ ) THEN
             CALL DLASRT( 'I', M, W, IINFO )
             IF( IINFO.NE.0 ) THEN