diff --git a/SRC/cgels.f b/SRC/cgels.f
index 39d94946c..08b2a8d80 100644
--- a/SRC/cgels.f
+++ b/SRC/cgels.f
@@ -493,21 +493,30 @@ SUBROUTINE CGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK,
 *
       END IF
 *
-*     Undo scaling
-*
-      IF( IASCL.EQ.1 ) THEN
-         CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IASCL.EQ.2 ) THEN
-         CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      END IF
-      IF( IBSCL.EQ.1 ) THEN
-         CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-         CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+         CALL CLASCL( 'G', 0, 0, ANRM, BNRM, SCLLEN, NRHS, B, LDB,
      $                INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+            CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+            CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+            CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+            CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         END IF
       END IF
 *
    50 CONTINUE
diff --git a/SRC/cgelsd.f b/SRC/cgelsd.f
index 11ab07ac2..d12c293f0 100644
--- a/SRC/cgelsd.f
+++ b/SRC/cgelsd.f
@@ -642,26 +642,38 @@ SUBROUTINE CGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
 *
       END IF
 *
-*     Undo scaling.
-*
-      IF( IASCL.EQ.1 ) THEN
-         CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+         CALL CLASCL( 'G', 0, 0, ANRM, BNRM, N, NRHS, B, LDB,
      $                INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+            CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+            CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+            CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+            CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+      END IF
+      IF( IASCL.EQ.1 ) THEN
          CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
      $                INFO )
       ELSE IF( IASCL.EQ.2 ) THEN
-         CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
-     $                INFO )
          CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
      $                INFO )
       END IF
-      IF( IBSCL.EQ.1 ) THEN
-         CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-         CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      END IF
 *
    10 CONTINUE
       WORK( 1 ) = SROUNDUP_LWORK(MAXWRK)
diff --git a/SRC/cgelss.f b/SRC/cgelss.f
index fbadfa09e..0a13e03c2 100644
--- a/SRC/cgelss.f
+++ b/SRC/cgelss.f
@@ -764,26 +764,38 @@ SUBROUTINE CGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
          END IF
       END IF
 *
-*     Undo scaling
-*
-      IF( IASCL.EQ.1 ) THEN
-         CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+         CALL CLASCL( 'G', 0, 0, ANRM, BNRM, N, NRHS, B, LDB,
      $                INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+            CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+            CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+            CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+            CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+      END IF
+      IF( IASCL.EQ.1 ) THEN
          CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
      $                INFO )
       ELSE IF( IASCL.EQ.2 ) THEN
-         CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
-     $                INFO )
          CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
      $                INFO )
       END IF
-      IF( IBSCL.EQ.1 ) THEN
-         CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-         CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      END IF
    70 CONTINUE
       WORK( 1 ) = SROUNDUP_LWORK(MAXWRK)
       RETURN
diff --git a/SRC/cgelst.f b/SRC/cgelst.f
index 6a938ec09..25c7b8506 100644
--- a/SRC/cgelst.f
+++ b/SRC/cgelst.f
@@ -523,21 +523,30 @@ SUBROUTINE CGELST( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK,
 *
       END IF
 *
-*     Undo scaling
-*
-      IF( IASCL.EQ.1 ) THEN
-         CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IASCL.EQ.2 ) THEN
-         CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      END IF
-      IF( IBSCL.EQ.1 ) THEN
-         CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-         CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+         CALL CLASCL( 'G', 0, 0, ANRM, BNRM, SCLLEN, NRHS, B, LDB,
      $                INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+            CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+            CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+            CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+            CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         END IF
       END IF
 *
       WORK( 1 ) = SROUNDUP_LWORK( LWOPT )
diff --git a/SRC/cgelsy.f b/SRC/cgelsy.f
index e8dfe3bbc..b1f28afd9 100644
--- a/SRC/cgelsy.f
+++ b/SRC/cgelsy.f
@@ -452,26 +452,38 @@ SUBROUTINE CGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND,
 *
 *     complex workspace: N.
 *
-*     Undo scaling
-*
-      IF( IASCL.EQ.1 ) THEN
-         CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+         CALL CLASCL( 'G', 0, 0, ANRM, BNRM, N, NRHS, B, LDB,
      $                INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+            CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+            CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+            CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+            CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+      END IF
+      IF( IASCL.EQ.1 ) THEN
          CALL CLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA,
      $                INFO )
       ELSE IF( IASCL.EQ.2 ) THEN
-         CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
-     $                INFO )
          CALL CLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA,
      $                INFO )
       END IF
-      IF( IBSCL.EQ.1 ) THEN
-         CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-         CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      END IF
 *
    70 CONTINUE
       WORK( 1 ) = CMPLX( LWKOPT )
diff --git a/SRC/cgetsls.f b/SRC/cgetsls.f
index dcd9a0dd9..6ed933734 100644
--- a/SRC/cgetsls.f
+++ b/SRC/cgetsls.f
@@ -480,21 +480,30 @@ SUBROUTINE CGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
 *
       END IF
 *
-*     Undo scaling
-*
-      IF( IASCL.EQ.1 ) THEN
-        CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
-     $               INFO )
-      ELSE IF( IASCL.EQ.2 ) THEN
-        CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
-     $               INFO )
-      END IF
-      IF( IBSCL.EQ.1 ) THEN
-        CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
-     $               INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-        CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+        CALL CLASCL( 'G', 0, 0, ANRM, BNRM, SCLLEN, NRHS, B, LDB,
      $               INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+           CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
+     $                  INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+           CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
+     $                  INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+           CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                  INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+           CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                  INFO )
+         END IF
       END IF
 *
    50 CONTINUE
diff --git a/SRC/dgels.f b/SRC/dgels.f
index 21900fe59..3a3ff90eb 100644
--- a/SRC/dgels.f
+++ b/SRC/dgels.f
@@ -489,21 +489,30 @@ SUBROUTINE DGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK,
 *
       END IF
 *
-*     Undo scaling
-*
-      IF( IASCL.EQ.1 ) THEN
-         CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IASCL.EQ.2 ) THEN
-         CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      END IF
-      IF( IBSCL.EQ.1 ) THEN
-         CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-         CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+         CALL DLASCL( 'G', 0, 0, ANRM, BNRM, SCLLEN, NRHS, B, LDB,
      $                INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+            CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+            CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+            CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+            CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         END IF
       END IF
 *
    50 CONTINUE
diff --git a/SRC/dgelsd.f b/SRC/dgelsd.f
index 0043f2593..af8fd6f72 100644
--- a/SRC/dgelsd.f
+++ b/SRC/dgelsd.f
@@ -605,26 +605,38 @@ SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
 *
       END IF
 *
-*     Undo scaling.
-*
-      IF( IASCL.EQ.1 ) THEN
-         CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+         CALL DLASCL( 'G', 0, 0, ANRM, BNRM, N, NRHS, B, LDB,
      $                INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+            CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+            CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+            CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+            CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+      END IF
+      IF( IASCL.EQ.1 ) THEN
          CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
      $                INFO )
       ELSE IF( IASCL.EQ.2 ) THEN
-         CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
-     $                INFO )
          CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
      $                INFO )
       END IF
-      IF( IBSCL.EQ.1 ) THEN
-         CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-         CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      END IF
 *
    10 CONTINUE
       WORK( 1 ) = MAXWRK
diff --git a/SRC/dgelss.f b/SRC/dgelss.f
index 740820763..1ad1de72e 100644
--- a/SRC/dgelss.f
+++ b/SRC/dgelss.f
@@ -734,26 +734,38 @@ SUBROUTINE DGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
          END IF
       END IF
 *
-*     Undo scaling
-*
-      IF( IASCL.EQ.1 ) THEN
-         CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+         CALL DLASCL( 'G', 0, 0, ANRM, BNRM, N, NRHS, B, LDB,
      $                INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+            CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+            CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+            CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+            CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+      END IF
+      IF( IASCL.EQ.1 ) THEN
          CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
      $                INFO )
       ELSE IF( IASCL.EQ.2 ) THEN
-         CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
-     $                INFO )
          CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
      $                INFO )
       END IF
-      IF( IBSCL.EQ.1 ) THEN
-         CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-         CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      END IF
 *
    70 CONTINUE
       WORK( 1 ) = MAXWRK
diff --git a/SRC/dgelst.f b/SRC/dgelst.f
index 6e69a7159..0ff05ee6f 100644
--- a/SRC/dgelst.f
+++ b/SRC/dgelst.f
@@ -519,21 +519,30 @@ SUBROUTINE DGELST( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK,
 *
       END IF
 *
-*     Undo scaling
-*
-      IF( IASCL.EQ.1 ) THEN
-         CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IASCL.EQ.2 ) THEN
-         CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      END IF
-      IF( IBSCL.EQ.1 ) THEN
-         CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-         CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+         CALL DLASCL( 'G', 0, 0, ANRM, BNRM, SCLLEN, NRHS, B, LDB,
      $                INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+            CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+            CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+            CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+            CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         END IF
       END IF
 *
       WORK( 1 ) = DBLE( LWOPT )
diff --git a/SRC/dgelsy.f b/SRC/dgelsy.f
index 7f51af280..0a087550c 100644
--- a/SRC/dgelsy.f
+++ b/SRC/dgelsy.f
@@ -454,26 +454,38 @@ SUBROUTINE DGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND,
 *
 *     workspace: N.
 *
-*     Undo scaling
-*
-      IF( IASCL.EQ.1 ) THEN
-         CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+         CALL DLASCL( 'G', 0, 0, ANRM, BNRM, N, NRHS, B, LDB,
      $                INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+            CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+            CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+            CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+            CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+      END IF
+      IF( IASCL.EQ.1 ) THEN
          CALL DLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA,
      $                INFO )
       ELSE IF( IASCL.EQ.2 ) THEN
-         CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
-     $                INFO )
          CALL DLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA,
      $                INFO )
       END IF
-      IF( IBSCL.EQ.1 ) THEN
-         CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-         CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      END IF
 *
    70 CONTINUE
       WORK( 1 ) = LWKOPT
diff --git a/SRC/dgetsls.f b/SRC/dgetsls.f
index 5b61d5888..b5e29593d 100644
--- a/SRC/dgetsls.f
+++ b/SRC/dgetsls.f
@@ -476,21 +476,30 @@ SUBROUTINE DGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
 *
       END IF
 *
-*     Undo scaling
-*
-      IF( IASCL.EQ.1 ) THEN
-        CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
-     $               INFO )
-      ELSE IF( IASCL.EQ.2 ) THEN
-        CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
-     $               INFO )
-      END IF
-      IF( IBSCL.EQ.1 ) THEN
-        CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
-     $               INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-        CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+        CALL DLASCL( 'G', 0, 0, ANRM, BNRM, SCLLEN, NRHS, B, LDB,
      $               INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+           CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
+     $                  INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+           CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
+     $                  INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+           CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                  INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+           CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                  INFO )
+         END IF
       END IF
 *
    50 CONTINUE
diff --git a/SRC/sgels.f b/SRC/sgels.f
index e87e1fa94..751e67d4f 100644
--- a/SRC/sgels.f
+++ b/SRC/sgels.f
@@ -490,21 +490,30 @@ SUBROUTINE SGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK,
 *
       END IF
 *
-*     Undo scaling
-*
-      IF( IASCL.EQ.1 ) THEN
-         CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IASCL.EQ.2 ) THEN
-         CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      END IF
-      IF( IBSCL.EQ.1 ) THEN
-         CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-         CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+         CALL SLASCL( 'G', 0, 0, ANRM, BNRM, SCLLEN, NRHS, B, LDB,
      $                INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+            CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+            CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+            CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+            CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         END IF
       END IF
 *
    50 CONTINUE
diff --git a/SRC/sgelsd.f b/SRC/sgelsd.f
index 423450e1c..051cfd871 100644
--- a/SRC/sgelsd.f
+++ b/SRC/sgelsd.f
@@ -610,26 +610,38 @@ SUBROUTINE SGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND,
 *
       END IF
 *
-*     Undo scaling.
-*
-      IF( IASCL.EQ.1 ) THEN
-         CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+         CALL SLASCL( 'G', 0, 0, ANRM, BNRM, N, NRHS, B, LDB,
      $                INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+            CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+            CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+            CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+            CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+      END IF
+      IF( IASCL.EQ.1 ) THEN
          CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
      $                INFO )
       ELSE IF( IASCL.EQ.2 ) THEN
-         CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
-     $                INFO )
          CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
      $                INFO )
       END IF
-      IF( IBSCL.EQ.1 ) THEN
-         CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-         CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      END IF
 *
    10 CONTINUE
       WORK( 1 ) = SROUNDUP_LWORK(MAXWRK)
diff --git a/SRC/sgelss.f b/SRC/sgelss.f
index 99371e324..b0b3c92e6 100644
--- a/SRC/sgelss.f
+++ b/SRC/sgelss.f
@@ -731,26 +731,38 @@ SUBROUTINE SGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
          END IF
       END IF
 *
-*     Undo scaling
-*
-      IF( IASCL.EQ.1 ) THEN
-         CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+         CALL SLASCL( 'G', 0, 0, ANRM, BNRM, N, NRHS, B, LDB,
      $                INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+            CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+            CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+            CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+            CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+      END IF
+      IF( IASCL.EQ.1 ) THEN
          CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
      $                INFO )
       ELSE IF( IASCL.EQ.2 ) THEN
-         CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
-     $                INFO )
          CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
      $                INFO )
       END IF
-      IF( IBSCL.EQ.1 ) THEN
-         CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-         CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      END IF
 *
    70 CONTINUE
       WORK( 1 ) = SROUNDUP_LWORK(MAXWRK)
diff --git a/SRC/sgelst.f b/SRC/sgelst.f
index 92c37a370..f7b421807 100644
--- a/SRC/sgelst.f
+++ b/SRC/sgelst.f
@@ -519,21 +519,30 @@ SUBROUTINE SGELST( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK,
 *
       END IF
 *
-*     Undo scaling
-*
-      IF( IASCL.EQ.1 ) THEN
-         CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IASCL.EQ.2 ) THEN
-         CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      END IF
-      IF( IBSCL.EQ.1 ) THEN
-         CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-         CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+         CALL SLASCL( 'G', 0, 0, ANRM, BNRM, SCLLEN, NRHS, B, LDB,
      $                INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+            CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+            CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+            CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+            CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         END IF
       END IF
 *
       WORK( 1 ) = SROUNDUP_LWORK( LWOPT )
diff --git a/SRC/sgelsy.f b/SRC/sgelsy.f
index 22f65eac1..5582f02a3 100644
--- a/SRC/sgelsy.f
+++ b/SRC/sgelsy.f
@@ -455,26 +455,38 @@ SUBROUTINE SGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND,
 *
 *     workspace: N.
 *
-*     Undo scaling
-*
-      IF( IASCL.EQ.1 ) THEN
-         CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+         CALL SLASCL( 'G', 0, 0, ANRM, BNRM, N, NRHS, B, LDB,
      $                INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+            CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+            CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+            CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+            CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+      END IF
+      IF( IASCL.EQ.1 ) THEN
          CALL SLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA,
      $                INFO )
       ELSE IF( IASCL.EQ.2 ) THEN
-         CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
-     $                INFO )
          CALL SLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA,
      $                INFO )
       END IF
-      IF( IBSCL.EQ.1 ) THEN
-         CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-         CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      END IF
 *
    70 CONTINUE
       WORK( 1 ) = SROUNDUP_LWORK(LWKOPT)
diff --git a/SRC/sgetsls.f b/SRC/sgetsls.f
index 246d544f6..fb6841647 100644
--- a/SRC/sgetsls.f
+++ b/SRC/sgetsls.f
@@ -477,21 +477,30 @@ SUBROUTINE SGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
 *
       END IF
 *
-*     Undo scaling
-*
-      IF( IASCL.EQ.1 ) THEN
-        CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
-     $               INFO )
-      ELSE IF( IASCL.EQ.2 ) THEN
-        CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
-     $               INFO )
-      END IF
-      IF( IBSCL.EQ.1 ) THEN
-        CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
-     $               INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-        CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+        CALL SLASCL( 'G', 0, 0, ANRM, BNRM, SCLLEN, NRHS, B, LDB,
      $               INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+           CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
+     $                  INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+           CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
+     $                  INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+           CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                  INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+           CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                  INFO )
+         END IF
       END IF
 *
    50 CONTINUE
diff --git a/SRC/zgels.f b/SRC/zgels.f
index fb5351c18..7c0d9a536 100644
--- a/SRC/zgels.f
+++ b/SRC/zgels.f
@@ -493,21 +493,30 @@ SUBROUTINE ZGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK,
 *
       END IF
 *
-*     Undo scaling
-*
-      IF( IASCL.EQ.1 ) THEN
-         CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IASCL.EQ.2 ) THEN
-         CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      END IF
-      IF( IBSCL.EQ.1 ) THEN
-         CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-         CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+         CALL ZLASCL( 'G', 0, 0, ANRM, BNRM, SCLLEN, NRHS, B, LDB,
      $                INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+            CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+            CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+            CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+            CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         END IF
       END IF
 *
    50 CONTINUE
diff --git a/SRC/zgelsd.f b/SRC/zgelsd.f
index c5adf6c9e..7a16c2f89 100644
--- a/SRC/zgelsd.f
+++ b/SRC/zgelsd.f
@@ -641,26 +641,38 @@ SUBROUTINE ZGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
 *
       END IF
 *
-*     Undo scaling.
-*
-      IF( IASCL.EQ.1 ) THEN
-         CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+         CALL ZLASCL( 'G', 0, 0, ANRM, BNRM, N, NRHS, B, LDB,
      $                INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+            CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+            CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+            CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+            CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+      END IF
+      IF( IASCL.EQ.1 ) THEN
          CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
      $                INFO )
       ELSE IF( IASCL.EQ.2 ) THEN
-         CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
-     $                INFO )
          CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
      $                INFO )
       END IF
-      IF( IBSCL.EQ.1 ) THEN
-         CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-         CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      END IF
 *
    10 CONTINUE
       WORK( 1 ) = MAXWRK
diff --git a/SRC/zgelss.f b/SRC/zgelss.f
index 316133723..a1162a608 100644
--- a/SRC/zgelss.f
+++ b/SRC/zgelss.f
@@ -764,26 +764,38 @@ SUBROUTINE ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
          END IF
       END IF
 *
-*     Undo scaling
-*
-      IF( IASCL.EQ.1 ) THEN
-         CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+         CALL ZLASCL( 'G', 0, 0, ANRM, BNRM, N, NRHS, B, LDB,
      $                INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+            CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+            CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+            CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+            CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+      END IF
+      IF( IASCL.EQ.1 ) THEN
          CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
      $                INFO )
       ELSE IF( IASCL.EQ.2 ) THEN
-         CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
-     $                INFO )
          CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
      $                INFO )
       END IF
-      IF( IBSCL.EQ.1 ) THEN
-         CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-         CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      END IF
    70 CONTINUE
       WORK( 1 ) = MAXWRK
       RETURN
diff --git a/SRC/zgelst.f b/SRC/zgelst.f
index 1a5de8d68..894e40b69 100644
--- a/SRC/zgelst.f
+++ b/SRC/zgelst.f
@@ -523,21 +523,30 @@ SUBROUTINE ZGELST( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK,
 *
       END IF
 *
-*     Undo scaling
-*
-      IF( IASCL.EQ.1 ) THEN
-         CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IASCL.EQ.2 ) THEN
-         CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      END IF
-      IF( IBSCL.EQ.1 ) THEN
-         CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-         CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+         CALL ZLASCL( 'G', 0, 0, ANRM, BNRM, SCLLEN, NRHS, B, LDB,
      $                INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+            CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+            CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+            CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+            CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                   INFO )
+         END IF
       END IF
 *
       WORK( 1 ) = DBLE( LWOPT )
diff --git a/SRC/zgelsy.f b/SRC/zgelsy.f
index 4a0cdb935..091c31325 100644
--- a/SRC/zgelsy.f
+++ b/SRC/zgelsy.f
@@ -452,26 +452,38 @@ SUBROUTINE ZGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND,
 *
 *     complex workspace: N.
 *
-*     Undo scaling
-*
-      IF( IASCL.EQ.1 ) THEN
-         CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+         CALL ZLASCL( 'G', 0, 0, ANRM, BNRM, N, NRHS, B, LDB,
      $                INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+            CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+            CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+            CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+            CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
+     $                   INFO )
+         END IF
+      END IF
+      IF( IASCL.EQ.1 ) THEN
          CALL ZLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA,
      $                INFO )
       ELSE IF( IASCL.EQ.2 ) THEN
-         CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB,
-     $                INFO )
          CALL ZLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA,
      $                INFO )
       END IF
-      IF( IBSCL.EQ.1 ) THEN
-         CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-         CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB,
-     $                INFO )
-      END IF
 *
    70 CONTINUE
       WORK( 1 ) = DCMPLX( LWKOPT )
diff --git a/SRC/zgetsls.f b/SRC/zgetsls.f
index 86784da4c..c94f94edd 100644
--- a/SRC/zgetsls.f
+++ b/SRC/zgetsls.f
@@ -479,21 +479,30 @@ SUBROUTINE ZGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
 *
       END IF
 *
-*     Undo scaling
-*
-      IF( IASCL.EQ.1 ) THEN
-        CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
-     $               INFO )
-      ELSE IF( IASCL.EQ.2 ) THEN
-        CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
-     $                INFO )
-      END IF
-      IF( IBSCL.EQ.1 ) THEN
-        CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
-     $               INFO )
-      ELSE IF( IBSCL.EQ.2 ) THEN
-        CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+*     Undo scaling.  The solution carries the factor of B divided by
+*     the factor of A.  When both were scaled to the same end of the
+*     range the constants cancel and BNRM/ANRM is applied in one step:
+*     applied in two, the first step could flush or overflow an entry
+*     that the second would have brought back into range.
+*
+      IF( IASCL.EQ.IBSCL .AND. IASCL.NE.0 ) THEN
+        CALL ZLASCL( 'G', 0, 0, ANRM, BNRM, SCLLEN, NRHS, B, LDB,
      $               INFO )
+      ELSE
+         IF( IASCL.EQ.1 ) THEN
+           CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
+     $                  INFO )
+         ELSE IF( IASCL.EQ.2 ) THEN
+           CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
+     $                  INFO )
+         END IF
+         IF( IBSCL.EQ.1 ) THEN
+           CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                  INFO )
+         ELSE IF( IBSCL.EQ.2 ) THEN
+           CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+     $                  INFO )
+         END IF
       END IF
 *
    50 CONTINUE
diff --git a/TESTING/LIN/cdrvls.f b/TESTING/LIN/cdrvls.f
index ac27456a4..8f96c8943 100644
--- a/TESTING/LIN/cdrvls.f
+++ b/TESTING/LIN/cdrvls.f
@@ -200,6 +200,7 @@ SUBROUTINE CDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
       LOGICAL            TSTERR
       INTEGER            NM, NN, NNB, NNS, NOUT
       REAL               THRESH
+      REAL               SMLX
 *     ..
 *     .. Array Arguments ..
       LOGICAL            DOTYPE( * )
@@ -283,6 +284,7 @@ SUBROUTINE CDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
          ISEED( I ) = ISEEDY( I )
    10 CONTINUE
       EPS = SLAMCH( 'Epsilon' )
+      SMLX = 256*SLAMCH( 'Safe minimum' )
 *
 *     Threshold for rank estimation
 *
@@ -471,6 +473,14 @@ SUBROUTINE CDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
                                  CALL CSSCAL( NCOLS*NRHS,
      $                                        ONE / REAL( NCOLS ), WORK,
      $                                        1 )
+*
+*                             Make the last solution column small
+*                             enough that undoing the scaling of
+*                             the solution in two steps flushes it.
+*
+                                 IF( ISCALE.EQ.2 .AND. NRHS.GT.1 )
+     $               CALL CSSCAL( NCOLS, SMLX,
+     $                          WORK( ( NRHS-1 )*LDWORK+1 ), 1 )
                               END IF
                               CALL CGEMM( TRANS, 'No transpose', NROWS,
      $                                    NRHS, NCOLS, CONE, COPYA, LDA,
@@ -588,6 +598,14 @@ SUBROUTINE CDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
                                  CALL CSSCAL( NCOLS*NRHS,
      $                                        ONE / REAL( NCOLS ), WORK,
      $                                        1 )
+*
+*                             Make the last solution column small
+*                             enough that undoing the scaling of
+*                             the solution in two steps flushes it.
+*
+                                 IF( ISCALE.EQ.2 .AND. NRHS.GT.1 )
+     $               CALL CSSCAL( NCOLS, SMLX,
+     $                          WORK( ( NRHS-1 )*LDWORK+1 ), 1 )
                               END IF
                               CALL CGEMM( TRANS, 'No transpose', NROWS,
      $                                    NRHS, NCOLS, CONE, COPYA, LDA,
@@ -711,6 +729,14 @@ SUBROUTINE CDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
                                     CALL CSCAL( NCOLS*NRHS,
      $                                          CONE / REAL( NCOLS ),
      $                                          WORK, 1 )
+*
+*                             Make the last solution column small
+*                             enough that undoing the scaling of
+*                             the solution in two steps flushes it.
+*
+                                    IF( ISCALE.EQ.2 .AND. NRHS.GT.1 )
+     $               CALL CSSCAL( NCOLS, SMLX,
+     $                          WORK( ( NRHS-1 )*LDWORK+1 ), 1 )
                                  END IF
                                  CALL CGEMM( TRANS, 'No transpose',
      $                                       NROWS, NRHS, NCOLS, CONE,
diff --git a/TESTING/LIN/cqrt13.f b/TESTING/LIN/cqrt13.f
index e87b7eafb..b94dc9991 100644
--- a/TESTING/LIN/cqrt13.f
+++ b/TESTING/LIN/cqrt13.f
@@ -111,7 +111,7 @@ SUBROUTINE CQRT13( SCALE, M, N, A, LDA, NORMA, ISEED )
 *     ..
 *     .. Local Scalars ..
       INTEGER            INFO, J
-      REAL               BIGNUM, SMLNUM
+      REAL               BIGNUM, BIGUP, SMLNUM
 *     ..
 *     .. External Functions ..
       REAL               CLANGE, SCASUM, SLAMCH
@@ -149,12 +149,18 @@ SUBROUTINE CQRT13( SCALE, M, N, A, LDA, NORMA, ISEED )
          BIGNUM = ONE / SMLNUM
          SMLNUM = SMLNUM / SLAMCH( 'Epsilon' )
          BIGNUM = ONE / SMLNUM
+*
+*        The least squares drivers scale a matrix whose largest
+*        entry lies outside [SMLNUM, BIGNUM].  Scale up past
+*        that bound, so that their scaling is exercised.
+*
+         BIGUP = SLAMCH( 'Overflow' ) / 256
 *
          IF( SCALE.EQ.2 ) THEN
 *
 *           matrix scaled up
 *
-            CALL CLASCL( 'General', 0, 0, NORMA, BIGNUM, M, N, A, LDA,
+            CALL CLASCL( 'General', 0, 0, NORMA, BIGUP, M, N, A, LDA,
      $                   INFO )
          ELSE IF( SCALE.EQ.3 ) THEN
 *
diff --git a/TESTING/LIN/ddrvls.f b/TESTING/LIN/ddrvls.f
index 9fb85e0a4..5e2f614a6 100644
--- a/TESTING/LIN/ddrvls.f
+++ b/TESTING/LIN/ddrvls.f
@@ -199,6 +199,7 @@ SUBROUTINE DDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
       LOGICAL            TSTERR
       INTEGER            NM, NN, NNB, NNS, NOUT
       DOUBLE PRECISION   THRESH
+      DOUBLE PRECISION   SMLX
 *     ..
 *     .. Array Arguments ..
       LOGICAL            DOTYPE( * )
@@ -276,6 +277,7 @@ SUBROUTINE DDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
          ISEED( I ) = ISEEDY( I )
    10 CONTINUE
       EPS = DLAMCH( 'Epsilon' )
+      SMLX = 256*DLAMCH( 'Safe minimum' )
 *
 *     Threshold for rank estimation
 *
@@ -456,6 +458,14 @@ SUBROUTINE DDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
                                  CALL DSCAL( NCOLS*NRHS,
      $                                       ONE / DBLE( NCOLS ), WORK,
      $                                       1 )
+*
+*                             Make the last solution column small
+*                             enough that undoing the scaling of
+*                             the solution in two steps flushes it.
+*
+                                 IF( ISCALE.EQ.2 .AND. NRHS.GT.1 )
+     $               CALL DSCAL( NCOLS, SMLX,
+     $                          WORK( ( NRHS-1 )*LDWORK+1 ), 1 )
                               END IF
                               CALL DGEMM( TRANS, 'No transpose', NROWS,
      $                                    NRHS, NCOLS, ONE, COPYA, LDA,
@@ -573,6 +583,14 @@ SUBROUTINE DDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
                                  CALL DSCAL( NCOLS*NRHS,
      $                                       ONE / DBLE( NCOLS ), WORK,
      $                                       1 )
+*
+*                             Make the last solution column small
+*                             enough that undoing the scaling of
+*                             the solution in two steps flushes it.
+*
+                                 IF( ISCALE.EQ.2 .AND. NRHS.GT.1 )
+     $               CALL DSCAL( NCOLS, SMLX,
+     $                          WORK( ( NRHS-1 )*LDWORK+1 ), 1 )
                               END IF
                               CALL DGEMM( TRANS, 'No transpose', NROWS,
      $                                    NRHS, NCOLS, ONE, COPYA, LDA,
@@ -697,6 +715,14 @@ SUBROUTINE DDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
                                     CALL DSCAL( NCOLS*NRHS,
      $                                          ONE / DBLE( NCOLS ),
      $                                          WORK, 1 )
+*
+*                             Make the last solution column small
+*                             enough that undoing the scaling of
+*                             the solution in two steps flushes it.
+*
+                                    IF( ISCALE.EQ.2 .AND. NRHS.GT.1 )
+     $               CALL DSCAL( NCOLS, SMLX,
+     $                          WORK( ( NRHS-1 )*LDWORK+1 ), 1 )
                                  END IF
                                  CALL DGEMM( TRANS, 'No transpose',
      $                                       NROWS, NRHS, NCOLS, ONE,
diff --git a/TESTING/LIN/dqrt13.f b/TESTING/LIN/dqrt13.f
index eaeefb307..8e4cb9a5b 100644
--- a/TESTING/LIN/dqrt13.f
+++ b/TESTING/LIN/dqrt13.f
@@ -111,7 +111,7 @@ SUBROUTINE DQRT13( SCALE, M, N, A, LDA, NORMA, ISEED )
 *     ..
 *     .. Local Scalars ..
       INTEGER            INFO, J
-      DOUBLE PRECISION   BIGNUM, SMLNUM
+      DOUBLE PRECISION   BIGNUM, BIGUP, SMLNUM
 *     ..
 *     .. External Functions ..
       DOUBLE PRECISION   DASUM, DLAMCH, DLANGE
@@ -149,12 +149,18 @@ SUBROUTINE DQRT13( SCALE, M, N, A, LDA, NORMA, ISEED )
          BIGNUM = ONE / SMLNUM
          SMLNUM = SMLNUM / DLAMCH( 'Epsilon' )
          BIGNUM = ONE / SMLNUM
+*
+*        The least squares drivers scale a matrix whose largest
+*        entry lies outside [SMLNUM, BIGNUM].  Scale up past
+*        that bound, so that their scaling is exercised.
+*
+         BIGUP = DLAMCH( 'Overflow' ) / 256
 *
          IF( SCALE.EQ.2 ) THEN
 *
 *           matrix scaled up
 *
-            CALL DLASCL( 'General', 0, 0, NORMA, BIGNUM, M, N, A, LDA,
+            CALL DLASCL( 'General', 0, 0, NORMA, BIGUP, M, N, A, LDA,
      $                   INFO )
          ELSE IF( SCALE.EQ.3 ) THEN
 *
diff --git a/TESTING/LIN/sdrvls.f b/TESTING/LIN/sdrvls.f
index bb3ea9a7c..cee0c68f9 100644
--- a/TESTING/LIN/sdrvls.f
+++ b/TESTING/LIN/sdrvls.f
@@ -199,6 +199,7 @@ SUBROUTINE SDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
       LOGICAL            TSTERR
       INTEGER            NM, NN, NNB, NNS, NOUT
       REAL               THRESH
+      REAL               SMLX
 *     ..
 *     .. Array Arguments ..
       LOGICAL            DOTYPE( * )
@@ -276,6 +277,7 @@ SUBROUTINE SDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
          ISEED( I ) = ISEEDY( I )
    10 CONTINUE
       EPS = SLAMCH( 'Epsilon' )
+      SMLX = 256*SLAMCH( 'Safe minimum' )
 *
 *     Threshold for rank estimation
 *
@@ -456,6 +458,14 @@ SUBROUTINE SDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
                                  CALL SSCAL( NCOLS*NRHS,
      $                                       ONE / REAL( NCOLS ), WORK,
      $                                       1 )
+*
+*                             Make the last solution column small
+*                             enough that undoing the scaling of
+*                             the solution in two steps flushes it.
+*
+                                 IF( ISCALE.EQ.2 .AND. NRHS.GT.1 )
+     $               CALL SSCAL( NCOLS, SMLX,
+     $                          WORK( ( NRHS-1 )*LDWORK+1 ), 1 )
                               END IF
                               CALL SGEMM( TRANS, 'No transpose', NROWS,
      $                                    NRHS, NCOLS, ONE, COPYA, LDA,
@@ -573,6 +583,14 @@ SUBROUTINE SDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
                                  CALL SSCAL( NCOLS*NRHS,
      $                                       ONE / REAL( NCOLS ), WORK,
      $                                       1 )
+*
+*                             Make the last solution column small
+*                             enough that undoing the scaling of
+*                             the solution in two steps flushes it.
+*
+                                 IF( ISCALE.EQ.2 .AND. NRHS.GT.1 )
+     $               CALL SSCAL( NCOLS, SMLX,
+     $                          WORK( ( NRHS-1 )*LDWORK+1 ), 1 )
                               END IF
                               CALL SGEMM( TRANS, 'No transpose', NROWS,
      $                                    NRHS, NCOLS, ONE, COPYA, LDA,
@@ -697,6 +715,14 @@ SUBROUTINE SDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
                                     CALL SSCAL( NCOLS*NRHS,
      $                                          ONE / REAL( NCOLS ),
      $                                          WORK, 1 )
+*
+*                             Make the last solution column small
+*                             enough that undoing the scaling of
+*                             the solution in two steps flushes it.
+*
+                                    IF( ISCALE.EQ.2 .AND. NRHS.GT.1 )
+     $               CALL SSCAL( NCOLS, SMLX,
+     $                          WORK( ( NRHS-1 )*LDWORK+1 ), 1 )
                                  END IF
                                  CALL SGEMM( TRANS, 'No transpose',
      $                                       NROWS, NRHS, NCOLS, ONE,
diff --git a/TESTING/LIN/sqrt13.f b/TESTING/LIN/sqrt13.f
index bb153a67b..d49cd9f5f 100644
--- a/TESTING/LIN/sqrt13.f
+++ b/TESTING/LIN/sqrt13.f
@@ -111,7 +111,7 @@ SUBROUTINE SQRT13( SCALE, M, N, A, LDA, NORMA, ISEED )
 *     ..
 *     .. Local Scalars ..
       INTEGER            INFO, J
-      REAL               BIGNUM, SMLNUM
+      REAL               BIGNUM, BIGUP, SMLNUM
 *     ..
 *     .. External Functions ..
       REAL               SASUM, SLAMCH, SLANGE
@@ -149,12 +149,18 @@ SUBROUTINE SQRT13( SCALE, M, N, A, LDA, NORMA, ISEED )
          BIGNUM = ONE / SMLNUM
          SMLNUM = SMLNUM / SLAMCH( 'Epsilon' )
          BIGNUM = ONE / SMLNUM
+*
+*        The least squares drivers scale a matrix whose largest
+*        entry lies outside [SMLNUM, BIGNUM].  Scale up past
+*        that bound, so that their scaling is exercised.
+*
+         BIGUP = SLAMCH( 'Overflow' ) / 256
 *
          IF( SCALE.EQ.2 ) THEN
 *
 *           matrix scaled up
 *
-            CALL SLASCL( 'General', 0, 0, NORMA, BIGNUM, M, N, A, LDA,
+            CALL SLASCL( 'General', 0, 0, NORMA, BIGUP, M, N, A, LDA,
      $                   INFO )
          ELSE IF( SCALE.EQ.3 ) THEN
 *
diff --git a/TESTING/LIN/zdrvls.f b/TESTING/LIN/zdrvls.f
index 2fab34f59..e158d6bbd 100644
--- a/TESTING/LIN/zdrvls.f
+++ b/TESTING/LIN/zdrvls.f
@@ -199,6 +199,7 @@ SUBROUTINE ZDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
       LOGICAL            TSTERR
       INTEGER            NM, NN, NNB, NNS, NOUT
       DOUBLE PRECISION   THRESH
+      DOUBLE PRECISION   SMLX
 *     ..
 *     .. Array Arguments ..
       LOGICAL            DOTYPE( * )
@@ -282,6 +283,7 @@ SUBROUTINE ZDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
          ISEED( I ) = ISEEDY( I )
    10 CONTINUE
       EPS = DLAMCH( 'Epsilon' )
+      SMLX = 256*DLAMCH( 'Safe minimum' )
 *
 *     Threshold for rank estimation
 *
@@ -470,6 +472,14 @@ SUBROUTINE ZDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
                                  CALL ZDSCAL( NCOLS*NRHS,
      $                                        ONE / DBLE( NCOLS ), WORK,
      $                                        1 )
+*
+*                             Make the last solution column small
+*                             enough that undoing the scaling of
+*                             the solution in two steps flushes it.
+*
+                                 IF( ISCALE.EQ.2 .AND. NRHS.GT.1 )
+     $               CALL ZDSCAL( NCOLS, SMLX,
+     $                          WORK( ( NRHS-1 )*LDWORK+1 ), 1 )
                               END IF
                               CALL ZGEMM( TRANS, 'No transpose', NROWS,
      $                                    NRHS, NCOLS, CONE, COPYA, LDA,
@@ -587,6 +597,14 @@ SUBROUTINE ZDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
                                  CALL ZDSCAL( NCOLS*NRHS,
      $                                        ONE / DBLE( NCOLS ), WORK,
      $                                        1 )
+*
+*                             Make the last solution column small
+*                             enough that undoing the scaling of
+*                             the solution in two steps flushes it.
+*
+                                 IF( ISCALE.EQ.2 .AND. NRHS.GT.1 )
+     $               CALL ZDSCAL( NCOLS, SMLX,
+     $                          WORK( ( NRHS-1 )*LDWORK+1 ), 1 )
                               END IF
                               CALL ZGEMM( TRANS, 'No transpose', NROWS,
      $                                    NRHS, NCOLS, CONE, COPYA, LDA,
@@ -710,6 +728,14 @@ SUBROUTINE ZDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
                                     CALL ZSCAL( NCOLS*NRHS,
      $                                          CONE / DBLE( NCOLS ),
      $                                          WORK, 1 )
+*
+*                             Make the last solution column small
+*                             enough that undoing the scaling of
+*                             the solution in two steps flushes it.
+*
+                                    IF( ISCALE.EQ.2 .AND. NRHS.GT.1 )
+     $               CALL ZDSCAL( NCOLS, SMLX,
+     $                          WORK( ( NRHS-1 )*LDWORK+1 ), 1 )
                                  END IF
                                  CALL ZGEMM( TRANS, 'No transpose',
      $                                       NROWS, NRHS, NCOLS, CONE,
diff --git a/TESTING/LIN/zqrt13.f b/TESTING/LIN/zqrt13.f
index 67c4884f7..ff21a1875 100644
--- a/TESTING/LIN/zqrt13.f
+++ b/TESTING/LIN/zqrt13.f
@@ -111,7 +111,7 @@ SUBROUTINE ZQRT13( SCALE, M, N, A, LDA, NORMA, ISEED )
 *     ..
 *     .. Local Scalars ..
       INTEGER            INFO, J
-      DOUBLE PRECISION   BIGNUM, SMLNUM
+      DOUBLE PRECISION   BIGNUM, BIGUP, SMLNUM
 *     ..
 *     .. External Functions ..
       DOUBLE PRECISION   DLAMCH, DZASUM, ZLANGE
@@ -149,12 +149,18 @@ SUBROUTINE ZQRT13( SCALE, M, N, A, LDA, NORMA, ISEED )
          BIGNUM = ONE / SMLNUM
          SMLNUM = SMLNUM / DLAMCH( 'Epsilon' )
          BIGNUM = ONE / SMLNUM
+*
+*        The least squares drivers scale a matrix whose largest
+*        entry lies outside [SMLNUM, BIGNUM].  Scale up past
+*        that bound, so that their scaling is exercised.
+*
+         BIGUP = DLAMCH( 'Overflow' ) / 256
 *
          IF( SCALE.EQ.2 ) THEN
 *
 *           matrix scaled up
 *
-            CALL ZLASCL( 'General', 0, 0, NORMA, BIGNUM, M, N, A, LDA,
+            CALL ZLASCL( 'General', 0, 0, NORMA, BIGUP, M, N, A, LDA,
      $                   INFO )
          ELSE IF( SCALE.EQ.3 ) THEN
 *