From 4bc258e12aab904fe15095345e3450f605995cbd Mon Sep 17 00:00:00 2001
From: Rasmus Munk Larsen <rmlarsen@gmail.com>
Date: Mon, 7 Sep 2026 21:11:25 -0700
Subject: [PATCH] Undo the scaling of the least-squares solution in one step
 when A and B were scaled to the same end of the range

xGELS, xGELST, xGETSLS, xGELSY, xGELSD and xGELSS scale A and B into
[SMLNUM, BIGNUM] before solving and undo the two scalings on the
solution one after the other.  When both were scaled to the same end
the factors share a constant that cancels, but the first step runs on
its own: with both above BIGNUM it multiplies by BIGNUM/ANRM, as small
as 2^-54, and flushes any solution entry below 2^-1020 that the second
step, BNRM/BIGNUM, would have restored.  For A = 2^1023 I and
b = (2^1023, 2^-27), whose solution is (1, 2^-1050), all 24 drivers
return (1, 0) with INFO = 0.  The mirror case below SMLNUM applies the
large factor first and could overflow, though a full-rank problem
cannot reach it.

When IASCL = IBSCL /= 0, apply the quotient BNRM/ANRM in one xLASCL
call, which never overshoots its target; the other combinations apply
two factors in the same direction and are unchanged.  In xGELSY,
xGELSD and xGELSS the rescaling of R or S by the factor of A alone
moves into its own IF and is unchanged.  Where both factors apply the
solution is rounded once instead of twice and can differ from the
previous result in the last bit.

Found while reviewing the first revision of the xGGLSE/xGGGLM scaling
change, which copied this block.

The test suite never ran the block: xQRT13 scales its "scaled up"
matrix to 1/(SFMIN/EPS) = 2^969, one binade inside the drivers'
threshold of 2^970, so IASCL and IBSCL were zero for every matrix it
generates.  xQRT13 now scales up to a fixed 2^1016 instead, and xDRVLS
makes the last column of the exact solution 256*SFMIN in the scaled-up
types, which is small enough for the first of the two undo steps to
flush it.  xQRT16 normalizes its residual per right-hand side, so the
lost column is visible: the parent fails between 2395 and 5041 ratios
per precision, this branch none.

Over a sweep of 2028 (precision,
driver, exponent of A, exponent of b) cases every case whose code path
is unchanged is bit-identical to the parent, and the merged cases agree
with the unscaled twin as closely as before.  The full LAPACK test
suite passes with the same totals as the parent.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
---
 SRC/cgels.f          | 37 +++++++++++++++++++++++--------------
 SRC/cgelsd.f         | 38 +++++++++++++++++++++++++-------------
 SRC/cgelss.f         | 38 +++++++++++++++++++++++++-------------
 SRC/cgelst.f         | 37 +++++++++++++++++++++++--------------
 SRC/cgelsy.f         | 38 +++++++++++++++++++++++++-------------
 SRC/cgetsls.f        | 37 +++++++++++++++++++++++--------------
 SRC/dgels.f          | 37 +++++++++++++++++++++++--------------
 SRC/dgelsd.f         | 38 +++++++++++++++++++++++++-------------
 SRC/dgelss.f         | 38 +++++++++++++++++++++++++-------------
 SRC/dgelst.f         | 37 +++++++++++++++++++++++--------------
 SRC/dgelsy.f         | 38 +++++++++++++++++++++++++-------------
 SRC/dgetsls.f        | 37 +++++++++++++++++++++++--------------
 SRC/sgels.f          | 37 +++++++++++++++++++++++--------------
 SRC/sgelsd.f         | 38 +++++++++++++++++++++++++-------------
 SRC/sgelss.f         | 38 +++++++++++++++++++++++++-------------
 SRC/sgelst.f         | 37 +++++++++++++++++++++++--------------
 SRC/sgelsy.f         | 38 +++++++++++++++++++++++++-------------
 SRC/sgetsls.f        | 37 +++++++++++++++++++++++--------------
 SRC/zgels.f          | 37 +++++++++++++++++++++++--------------
 SRC/zgelsd.f         | 38 +++++++++++++++++++++++++-------------
 SRC/zgelss.f         | 38 +++++++++++++++++++++++++-------------
 SRC/zgelst.f         | 37 +++++++++++++++++++++++--------------
 SRC/zgelsy.f         | 38 +++++++++++++++++++++++++-------------
 SRC/zgetsls.f        | 37 +++++++++++++++++++++++--------------
 TESTING/LIN/cdrvls.f | 26 ++++++++++++++++++++++++++
 TESTING/LIN/cqrt13.f | 10 ++++++++--
 TESTING/LIN/ddrvls.f | 26 ++++++++++++++++++++++++++
 TESTING/LIN/dqrt13.f | 10 ++++++++--
 TESTING/LIN/sdrvls.f | 26 ++++++++++++++++++++++++++
 TESTING/LIN/sqrt13.f | 10 ++++++++--
 TESTING/LIN/zdrvls.f | 26 ++++++++++++++++++++++++++
 TESTING/LIN/zqrt13.f | 10 ++++++++--
 32 files changed, 712 insertions(+), 332 deletions(-)

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
 *