diff --git a/SRC/cggglm.f b/SRC/cggglm.f
index 46c4a2c73..0cf909f38 100644
--- a/SRC/cggglm.f
+++ b/SRC/cggglm.f
@@ -208,6 +208,8 @@ SUBROUTINE CGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK,
 *  ===================================================================
 *
 *     .. Parameters ..
+      REAL               ZERO, ONE
+      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
       COMPLEX            CZERO, CONE
       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
@@ -216,19 +218,24 @@ SUBROUTINE CGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK,
       LOGICAL            LQUERY
       INTEGER            I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3,
      $                   NB4, NP
+      INTEGER            IA, IB, IBIG, ID, ISML, KA, KB, KD
+      REAL               ANRM, BIGNUM, BNRM, DNRM, SMLNUM
+*     ..
+*     .. Local Arrays ..
+      REAL               RWORK( 1 )
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           CCOPY, CGEMV, CGGQRF, CTRTRS, CUNMQR,
-     $                   CUNMRQ,
-     $                   XERBLA
+      EXTERNAL           CCOPY, CGEMV, CGGQRF, CLASCL, CSSCAL, CTRTRS,
+     $                   CUNMQR, CUNMRQ, XERBLA
 *     ..
 *     .. External Functions ..
       INTEGER            ILAENV
       REAL               SROUNDUP_LWORK
-      EXTERNAL           ILAENV, SROUNDUP_LWORK
+      REAL               SLAMCH, CLANGE
+      EXTERNAL           CLANGE, ILAENV, SLAMCH, SROUNDUP_LWORK
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          INT, MAX, MIN
+      INTRINSIC          EXPONENT, HUGE, INT, MAX, MIN, SCALE
 *     ..
 *     .. Executable Statements ..
 *
@@ -290,6 +297,62 @@ SUBROUTINE CGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK,
          RETURN
       END IF
 *
+*     Get machine parameters
+*
+      SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' )
+      BIGNUM = ONE / SMLNUM
+      ISML = EXPONENT( SMLNUM )
+      IBIG = EXPONENT( BIGNUM ) - 1
+*
+*     Scale A, B and d independently by powers of two 2**KA, 2**KB and
+*     2**KD so that their largest entries lie in [SMLNUM,BIGNUM), the
+*     values whose EXPONENT lies in [ISML,IBIG].  x is then scaled by
+*     2**(KD-KA) and y by 2**(KD-KB).  Scaling by a power of two is
+*     exact.  A norm that is zero, infinite or NaN takes no part.
+*     CLASCL is called with both endpoints at or above one so that
+*     forming the factor raises no underflow.
+*
+      ANRM = CLANGE( 'M', N, M, A, LDA, RWORK )
+      KA = 0
+      IF( ANRM.GT.ZERO .AND. ANRM.LE.HUGE( ZERO ) ) THEN
+         IA = EXPONENT( ANRM )
+         IF( IA.LT.ISML ) THEN
+            KA = ISML - IA
+         ELSE IF( IA.GT.IBIG ) THEN
+            KA = IBIG - IA
+         END IF
+      END IF
+      IF( KA.NE.0 )
+     $   CALL CLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KA, 0 ) ),
+     $                SCALE( ONE, MAX( KA, 0 ) ), N, M, A, LDA, INFO )
+*
+      BNRM = CLANGE( 'M', N, P, B, LDB, RWORK )
+      KB = 0
+      IF( BNRM.GT.ZERO .AND. BNRM.LE.HUGE( ZERO ) ) THEN
+         IB = EXPONENT( BNRM )
+         IF( IB.LT.ISML ) THEN
+            KB = ISML - IB
+         ELSE IF( IB.GT.IBIG ) THEN
+            KB = IBIG - IB
+         END IF
+      END IF
+      IF( KB.NE.0 )
+     $   CALL CLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KB, 0 ) ),
+     $                SCALE( ONE, MAX( KB, 0 ) ), N, P, B, LDB, INFO )
+*
+      DNRM = CLANGE( 'M', N, 1, D, N, RWORK )
+      KD = 0
+      IF( DNRM.GT.ZERO .AND. DNRM.LE.HUGE( ZERO ) ) THEN
+         ID = EXPONENT( DNRM )
+         IF( ID.LT.ISML ) THEN
+            KD = ISML - ID
+         ELSE IF( ID.GT.IBIG ) THEN
+            KD = IBIG - ID
+         END IF
+      END IF
+      IF( KD.NE.0 )
+     $   CALL CSSCAL( N, SCALE( ONE, KD ), D, 1 )
+*
 *     Compute the GQR factorization of matrices A and B:
 *
 *          Q**H*A = ( R11 ) M,    Q**H*B*Z**H = ( T11   T12 ) M
@@ -359,6 +422,14 @@ SUBROUTINE CGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK,
       CALL CUNMRQ( 'Left', 'Conjugate transpose', P, 1, NP,
      $             B( MAX( 1, N-P+1 ), 1 ), LDB, WORK( M+1 ), Y,
      $             MAX( 1, P ), WORK( M+NP+1 ), LWORK-M-NP, INFO )
+*
+*     Undo scaling: x carries 2**(KD-KA) and y carries 2**(KD-KB)
+*
+      IF( KA.NE.KD )
+     $   CALL CSSCAL( M, SCALE( ONE, KA-KD ), X, 1 )
+      IF( KB.NE.KD )
+     $   CALL CSSCAL( P, SCALE( ONE, KB-KD ), Y, 1 )
+*
       WORK( 1 ) = CMPLX( M + NP + MAX( LOPT, INT( WORK( M+NP+1 ) ) ) )
 *
       RETURN
diff --git a/SRC/cgglse.f b/SRC/cgglse.f
index 0080717c4..1a24122bb 100644
--- a/SRC/cgglse.f
+++ b/SRC/cgglse.f
@@ -203,6 +203,8 @@ SUBROUTINE CGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK,
 *  =====================================================================
 *
 *     .. Parameters ..
+      REAL               ZERO, ONE
+      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
       COMPLEX            CONE
       PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
 *     ..
@@ -210,19 +212,24 @@ SUBROUTINE CGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK,
       LOGICAL            LQUERY
       INTEGER            LOPT, LWKMIN, LWKOPT, MN, NB, NB1, NB2, NB3,
      $                   NB4, NR
+      INTEGER            IA, IB, IBIG, IC, ID, ISML, IU, KA, KB, KT
+      REAL               ANRM, BIGNUM, BNRM, CNRM, DNRM, SMLNUM
+*     ..
+*     .. Local Arrays ..
+      REAL               RWORK( 1 )
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           CAXPY, CCOPY, CGEMV, CGGRQF, CTRMV,
-     $                   CTRTRS,
-     $                   CUNMQR, CUNMRQ, XERBLA
+      EXTERNAL           CAXPY, CCOPY, CGEMV, CGGRQF, CLASCL, CSSCAL,
+     $                   CTRMV, CTRTRS, CUNMQR, CUNMRQ, XERBLA
 *     ..
 *     .. External Functions ..
       INTEGER            ILAENV
       REAL               SROUNDUP_LWORK
-      EXTERNAL           ILAENV, SROUNDUP_LWORK
+      REAL               SLAMCH, CLANGE
+      EXTERNAL           CLANGE, ILAENV, SLAMCH, SROUNDUP_LWORK
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          INT, MAX, MIN
+      INTRINSIC          EXPONENT, HUGE, INT, MAX, MIN, SCALE
 *     ..
 *     .. Executable Statements ..
 *
@@ -277,6 +284,72 @@ SUBROUTINE CGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK,
       IF( N.EQ.0 )
      $   RETURN
 *
+*     Get machine parameters
+*
+      SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' )
+      BIGNUM = ONE / SMLNUM
+      ISML = EXPONENT( SMLNUM )
+      IBIG = EXPONENT( BIGNUM ) - 1
+*
+*     Scale A, B, c and d by powers of two so that their largest
+*     entries lie in [SMLNUM,BIGNUM), the values whose EXPONENT lies
+*     in [ISML,IBIG].  A and B are scaled independently, by 2**KA and
+*     2**KB.  c is scaled by 2**(KA+KT) and d by 2**(KB+KT), which
+*     keeps the constraint and the residual consistent and scales x by
+*     2**KT; KT brings the larger of the two right-hand sides into
+*     range.  Scaling by a power of two is exact.  A norm that is zero,
+*     infinite or NaN takes no part.  CLASCL is called with both
+*     endpoints at or above one so that forming the factor raises no
+*     underflow.
+*
+      ANRM = CLANGE( 'M', M, N, A, LDA, RWORK )
+      KA = 0
+      IF( ANRM.GT.ZERO .AND. ANRM.LE.HUGE( ZERO ) ) THEN
+         IA = EXPONENT( ANRM )
+         IF( IA.LT.ISML ) THEN
+            KA = ISML - IA
+         ELSE IF( IA.GT.IBIG ) THEN
+            KA = IBIG - IA
+         END IF
+      END IF
+      IF( KA.NE.0 )
+     $   CALL CLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KA, 0 ) ),
+     $                SCALE( ONE, MAX( KA, 0 ) ), M, N, A, LDA, INFO )
+*
+      BNRM = CLANGE( 'M', P, N, B, LDB, RWORK )
+      KB = 0
+      IF( BNRM.GT.ZERO .AND. BNRM.LE.HUGE( ZERO ) ) THEN
+         IB = EXPONENT( BNRM )
+         IF( IB.LT.ISML ) THEN
+            KB = ISML - IB
+         ELSE IF( IB.GT.IBIG ) THEN
+            KB = IBIG - IB
+         END IF
+      END IF
+      IF( KB.NE.0 )
+     $   CALL CLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KB, 0 ) ),
+     $                SCALE( ONE, MAX( KB, 0 ) ), P, N, B, LDB, INFO )
+*
+      CNRM = CLANGE( 'M', M, 1, C, MAX( 1, M ), RWORK )
+      DNRM = CLANGE( 'M', P, 1, D, MAX( 1, P ), RWORK )
+      IC = -HUGE( 0 )
+      IF( CNRM.GT.ZERO .AND. CNRM.LE.HUGE( ZERO ) )
+     $   IC = EXPONENT( CNRM ) + KA
+      ID = -HUGE( 0 )
+      IF( DNRM.GT.ZERO .AND. DNRM.LE.HUGE( ZERO ) )
+     $   ID = EXPONENT( DNRM ) + KB
+      IU = MAX( IC, ID )
+      KT = 0
+      IF( IU.GT.IBIG ) THEN
+         KT = IBIG - IU
+      ELSE IF( IU.LT.ISML .AND. IU.NE.-HUGE( 0 ) ) THEN
+         KT = ISML - IU
+      END IF
+      IF( KA+KT.NE.0 )
+     $   CALL CSSCAL( M, SCALE( ONE, KA+KT ), C, 1 )
+      IF( KB+KT.NE.0 )
+     $   CALL CSSCAL( P, SCALE( ONE, KB+KT ), D, 1 )
+*
 *     Compute the GRQ factorization of matrices B and A:
 *
 *            B*Q**H = (  0  T12 ) P   Z**H*A*Q**H = ( R11 R12 ) N-P
@@ -357,6 +430,15 @@ SUBROUTINE CGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK,
 *
       CALL CUNMRQ( 'Left', 'Conjugate Transpose', N, 1, P, B, LDB,
      $             WORK( 1 ), X, N, WORK( P+MN+1 ), LWORK-P-MN, INFO )
+*
+*     Undo scaling: x carries 2**KT and the residual in c(N-P+1:M)
+*     carries 2**(KA+KT)
+*
+      IF( KT.NE.0 )
+     $   CALL CSSCAL( N, SCALE( ONE, -KT ), X, 1 )
+      IF( KA+KT.NE.0 .AND. M+P.GT.N )
+     $   CALL CSSCAL( M+P-N, SCALE( ONE, -KA-KT ), C( N-P+1 ), 1 )
+*
       WORK( 1 ) = CMPLX( P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) ) )
 *
       RETURN
diff --git a/SRC/dggglm.f b/SRC/dggglm.f
index 32a3f8a9b..0a46b6f61 100644
--- a/SRC/dggglm.f
+++ b/SRC/dggglm.f
@@ -215,18 +215,23 @@ SUBROUTINE DGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK,
       LOGICAL            LQUERY
       INTEGER            I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3,
      $                   NB4, NP
+      INTEGER            IA, IB, IBIG, ID, ISML, KA, KB, KD
+      DOUBLE PRECISION   ANRM, BIGNUM, BNRM, DNRM, SMLNUM
+*     ..
+*     .. Local Arrays ..
+      DOUBLE PRECISION   RWORK( 1 )
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           DCOPY, DGEMV, DGGQRF, DORMQR, DORMRQ,
-     $                   DTRTRS,
-     $                   XERBLA
+      EXTERNAL           DCOPY, DGEMV, DGGQRF, DLASCL, DORMQR, DORMRQ,
+     $                   DSCAL, DTRTRS, XERBLA
 *     ..
 *     .. External Functions ..
       INTEGER            ILAENV
-      EXTERNAL           ILAENV
+      DOUBLE PRECISION   DLAMCH, DLANGE
+      EXTERNAL           DLAMCH, DLANGE, ILAENV
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          INT, MAX, MIN
+      INTRINSIC          EXPONENT, HUGE, INT, MAX, MIN, SCALE
 *     ..
 *     .. Executable Statements ..
 *
@@ -288,6 +293,62 @@ SUBROUTINE DGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK,
          RETURN
       END IF
 *
+*     Get machine parameters
+*
+      SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
+      BIGNUM = ONE / SMLNUM
+      ISML = EXPONENT( SMLNUM )
+      IBIG = EXPONENT( BIGNUM ) - 1
+*
+*     Scale A, B and d independently by powers of two 2**KA, 2**KB and
+*     2**KD so that their largest entries lie in [SMLNUM,BIGNUM), the
+*     values whose EXPONENT lies in [ISML,IBIG].  x is then scaled by
+*     2**(KD-KA) and y by 2**(KD-KB).  Scaling by a power of two is
+*     exact.  A norm that is zero, infinite or NaN takes no part.
+*     DLASCL is called with both endpoints at or above one so that
+*     forming the factor raises no underflow.
+*
+      ANRM = DLANGE( 'M', N, M, A, LDA, RWORK )
+      KA = 0
+      IF( ANRM.GT.ZERO .AND. ANRM.LE.HUGE( ZERO ) ) THEN
+         IA = EXPONENT( ANRM )
+         IF( IA.LT.ISML ) THEN
+            KA = ISML - IA
+         ELSE IF( IA.GT.IBIG ) THEN
+            KA = IBIG - IA
+         END IF
+      END IF
+      IF( KA.NE.0 )
+     $   CALL DLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KA, 0 ) ),
+     $                SCALE( ONE, MAX( KA, 0 ) ), N, M, A, LDA, INFO )
+*
+      BNRM = DLANGE( 'M', N, P, B, LDB, RWORK )
+      KB = 0
+      IF( BNRM.GT.ZERO .AND. BNRM.LE.HUGE( ZERO ) ) THEN
+         IB = EXPONENT( BNRM )
+         IF( IB.LT.ISML ) THEN
+            KB = ISML - IB
+         ELSE IF( IB.GT.IBIG ) THEN
+            KB = IBIG - IB
+         END IF
+      END IF
+      IF( KB.NE.0 )
+     $   CALL DLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KB, 0 ) ),
+     $                SCALE( ONE, MAX( KB, 0 ) ), N, P, B, LDB, INFO )
+*
+      DNRM = DLANGE( 'M', N, 1, D, N, RWORK )
+      KD = 0
+      IF( DNRM.GT.ZERO .AND. DNRM.LE.HUGE( ZERO ) ) THEN
+         ID = EXPONENT( DNRM )
+         IF( ID.LT.ISML ) THEN
+            KD = ISML - ID
+         ELSE IF( ID.GT.IBIG ) THEN
+            KD = IBIG - ID
+         END IF
+      END IF
+      IF( KD.NE.0 )
+     $   CALL DSCAL( N, SCALE( ONE, KD ), D, 1 )
+*
 *     Compute the GQR factorization of matrices A and B:
 *
 *          Q**T*A = ( R11 ) M,    Q**T*B*Z**T = ( T11   T12 ) M
@@ -355,6 +416,14 @@ SUBROUTINE DGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK,
       CALL DORMRQ( 'Left', 'Transpose', P, 1, NP,
      $             B( MAX( 1, N-P+1 ), 1 ), LDB, WORK( M+1 ), Y,
      $             MAX( 1, P ), WORK( M+NP+1 ), LWORK-M-NP, INFO )
+*
+*     Undo scaling: x carries 2**(KD-KA) and y carries 2**(KD-KB)
+*
+      IF( KA.NE.KD )
+     $   CALL DSCAL( M, SCALE( ONE, KA-KD ), X, 1 )
+      IF( KB.NE.KD )
+     $   CALL DSCAL( P, SCALE( ONE, KB-KD ), Y, 1 )
+*
       WORK( 1 ) = M + NP + MAX( LOPT, INT( WORK( M+NP+1 ) ) )
 *
       RETURN
diff --git a/SRC/dgglse.f b/SRC/dgglse.f
index 6c1198af1..23db6b0c5 100644
--- a/SRC/dgglse.f
+++ b/SRC/dgglse.f
@@ -203,25 +203,30 @@ SUBROUTINE DGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK,
 *  =====================================================================
 *
 *     .. Parameters ..
-      DOUBLE PRECISION   ONE
-      PARAMETER          ( ONE = 1.0D+0 )
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 *     ..
 *     .. Local Scalars ..
       LOGICAL            LQUERY
       INTEGER            LOPT, LWKMIN, LWKOPT, MN, NB, NB1, NB2, NB3,
      $                   NB4, NR
+      INTEGER            IA, IB, IBIG, IC, ID, ISML, IU, KA, KB, KT
+      DOUBLE PRECISION   ANRM, BIGNUM, BNRM, CNRM, DNRM, SMLNUM
+*     ..
+*     .. Local Arrays ..
+      DOUBLE PRECISION   RWORK( 1 )
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           DAXPY, DCOPY, DGEMV, DGGRQF, DORMQR,
-     $                   DORMRQ,
-     $                   DTRMV, DTRTRS, XERBLA
+      EXTERNAL           DAXPY, DCOPY, DGEMV, DGGRQF, DLASCL, DORMQR,
+     $                   DORMRQ, DSCAL, DTRMV, DTRTRS, XERBLA
 *     ..
 *     .. External Functions ..
       INTEGER            ILAENV
-      EXTERNAL           ILAENV
+      DOUBLE PRECISION   DLAMCH, DLANGE
+      EXTERNAL           DLAMCH, DLANGE, ILAENV
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          INT, MAX, MIN
+      INTRINSIC          EXPONENT, HUGE, INT, MAX, MIN, SCALE
 *     ..
 *     .. Executable Statements ..
 *
@@ -276,6 +281,72 @@ SUBROUTINE DGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK,
       IF( N.EQ.0 )
      $   RETURN
 *
+*     Get machine parameters
+*
+      SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
+      BIGNUM = ONE / SMLNUM
+      ISML = EXPONENT( SMLNUM )
+      IBIG = EXPONENT( BIGNUM ) - 1
+*
+*     Scale A, B, c and d by powers of two so that their largest
+*     entries lie in [SMLNUM,BIGNUM), the values whose EXPONENT lies
+*     in [ISML,IBIG].  A and B are scaled independently, by 2**KA and
+*     2**KB.  c is scaled by 2**(KA+KT) and d by 2**(KB+KT), which
+*     keeps the constraint and the residual consistent and scales x by
+*     2**KT; KT brings the larger of the two right-hand sides into
+*     range.  Scaling by a power of two is exact.  A norm that is zero,
+*     infinite or NaN takes no part.  DLASCL is called with both
+*     endpoints at or above one so that forming the factor raises no
+*     underflow.
+*
+      ANRM = DLANGE( 'M', M, N, A, LDA, RWORK )
+      KA = 0
+      IF( ANRM.GT.ZERO .AND. ANRM.LE.HUGE( ZERO ) ) THEN
+         IA = EXPONENT( ANRM )
+         IF( IA.LT.ISML ) THEN
+            KA = ISML - IA
+         ELSE IF( IA.GT.IBIG ) THEN
+            KA = IBIG - IA
+         END IF
+      END IF
+      IF( KA.NE.0 )
+     $   CALL DLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KA, 0 ) ),
+     $                SCALE( ONE, MAX( KA, 0 ) ), M, N, A, LDA, INFO )
+*
+      BNRM = DLANGE( 'M', P, N, B, LDB, RWORK )
+      KB = 0
+      IF( BNRM.GT.ZERO .AND. BNRM.LE.HUGE( ZERO ) ) THEN
+         IB = EXPONENT( BNRM )
+         IF( IB.LT.ISML ) THEN
+            KB = ISML - IB
+         ELSE IF( IB.GT.IBIG ) THEN
+            KB = IBIG - IB
+         END IF
+      END IF
+      IF( KB.NE.0 )
+     $   CALL DLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KB, 0 ) ),
+     $                SCALE( ONE, MAX( KB, 0 ) ), P, N, B, LDB, INFO )
+*
+      CNRM = DLANGE( 'M', M, 1, C, MAX( 1, M ), RWORK )
+      DNRM = DLANGE( 'M', P, 1, D, MAX( 1, P ), RWORK )
+      IC = -HUGE( 0 )
+      IF( CNRM.GT.ZERO .AND. CNRM.LE.HUGE( ZERO ) )
+     $   IC = EXPONENT( CNRM ) + KA
+      ID = -HUGE( 0 )
+      IF( DNRM.GT.ZERO .AND. DNRM.LE.HUGE( ZERO ) )
+     $   ID = EXPONENT( DNRM ) + KB
+      IU = MAX( IC, ID )
+      KT = 0
+      IF( IU.GT.IBIG ) THEN
+         KT = IBIG - IU
+      ELSE IF( IU.LT.ISML .AND. IU.NE.-HUGE( 0 ) ) THEN
+         KT = ISML - IU
+      END IF
+      IF( KA+KT.NE.0 )
+     $   CALL DSCAL( M, SCALE( ONE, KA+KT ), C, 1 )
+      IF( KB+KT.NE.0 )
+     $   CALL DSCAL( P, SCALE( ONE, KB+KT ), D, 1 )
+*
 *     Compute the GRQ factorization of matrices B and A:
 *
 *            B*Q**T = (  0  T12 ) P   Z**T*A*Q**T = ( R11 R12 ) N-P
@@ -357,6 +428,15 @@ SUBROUTINE DGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK,
       CALL DORMRQ( 'Left', 'Transpose', N, 1, P, B, LDB, WORK( 1 ),
      $             X,
      $             N, WORK( P+MN+1 ), LWORK-P-MN, INFO )
+*
+*     Undo scaling: x carries 2**KT and the residual in c(N-P+1:M)
+*     carries 2**(KA+KT)
+*
+      IF( KT.NE.0 )
+     $   CALL DSCAL( N, SCALE( ONE, -KT ), X, 1 )
+      IF( KA+KT.NE.0 .AND. M+P.GT.N )
+     $   CALL DSCAL( M+P-N, SCALE( ONE, -KA-KT ), C( N-P+1 ), 1 )
+*
       WORK( 1 ) = P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) )
 *
       RETURN
diff --git a/SRC/sggglm.f b/SRC/sggglm.f
index 10aea9a4f..caed7cf5d 100644
--- a/SRC/sggglm.f
+++ b/SRC/sggglm.f
@@ -215,19 +215,24 @@ SUBROUTINE SGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK,
       LOGICAL            LQUERY
       INTEGER            I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3,
      $                   NB4, NP
+      INTEGER            IA, IB, IBIG, ID, ISML, KA, KB, KD
+      REAL               ANRM, BIGNUM, BNRM, DNRM, SMLNUM
+*     ..
+*     .. Local Arrays ..
+      REAL               RWORK( 1 )
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           SCOPY, SGEMV, SGGQRF, SORMQR, SORMRQ,
-     $                   STRTRS,
-     $                   XERBLA
+      EXTERNAL           SCOPY, SGEMV, SGGQRF, SLASCL, SORMQR, SORMRQ,
+     $                   SSCAL, STRTRS, XERBLA
 *     ..
 *     .. External Functions ..
       INTEGER            ILAENV
       REAL               SROUNDUP_LWORK
-      EXTERNAL           ILAENV, SROUNDUP_LWORK
+      REAL               SLAMCH, SLANGE
+      EXTERNAL           ILAENV, SLAMCH, SLANGE, SROUNDUP_LWORK
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          INT, MAX, MIN
+      INTRINSIC          EXPONENT, HUGE, INT, MAX, MIN, SCALE
 *     ..
 *     .. Executable Statements ..
 *
@@ -289,6 +294,62 @@ SUBROUTINE SGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK,
          RETURN
       END IF
 *
+*     Get machine parameters
+*
+      SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' )
+      BIGNUM = ONE / SMLNUM
+      ISML = EXPONENT( SMLNUM )
+      IBIG = EXPONENT( BIGNUM ) - 1
+*
+*     Scale A, B and d independently by powers of two 2**KA, 2**KB and
+*     2**KD so that their largest entries lie in [SMLNUM,BIGNUM), the
+*     values whose EXPONENT lies in [ISML,IBIG].  x is then scaled by
+*     2**(KD-KA) and y by 2**(KD-KB).  Scaling by a power of two is
+*     exact.  A norm that is zero, infinite or NaN takes no part.
+*     SLASCL is called with both endpoints at or above one so that
+*     forming the factor raises no underflow.
+*
+      ANRM = SLANGE( 'M', N, M, A, LDA, RWORK )
+      KA = 0
+      IF( ANRM.GT.ZERO .AND. ANRM.LE.HUGE( ZERO ) ) THEN
+         IA = EXPONENT( ANRM )
+         IF( IA.LT.ISML ) THEN
+            KA = ISML - IA
+         ELSE IF( IA.GT.IBIG ) THEN
+            KA = IBIG - IA
+         END IF
+      END IF
+      IF( KA.NE.0 )
+     $   CALL SLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KA, 0 ) ),
+     $                SCALE( ONE, MAX( KA, 0 ) ), N, M, A, LDA, INFO )
+*
+      BNRM = SLANGE( 'M', N, P, B, LDB, RWORK )
+      KB = 0
+      IF( BNRM.GT.ZERO .AND. BNRM.LE.HUGE( ZERO ) ) THEN
+         IB = EXPONENT( BNRM )
+         IF( IB.LT.ISML ) THEN
+            KB = ISML - IB
+         ELSE IF( IB.GT.IBIG ) THEN
+            KB = IBIG - IB
+         END IF
+      END IF
+      IF( KB.NE.0 )
+     $   CALL SLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KB, 0 ) ),
+     $                SCALE( ONE, MAX( KB, 0 ) ), N, P, B, LDB, INFO )
+*
+      DNRM = SLANGE( 'M', N, 1, D, N, RWORK )
+      KD = 0
+      IF( DNRM.GT.ZERO .AND. DNRM.LE.HUGE( ZERO ) ) THEN
+         ID = EXPONENT( DNRM )
+         IF( ID.LT.ISML ) THEN
+            KD = ISML - ID
+         ELSE IF( ID.GT.IBIG ) THEN
+            KD = IBIG - ID
+         END IF
+      END IF
+      IF( KD.NE.0 )
+     $   CALL SSCAL( N, SCALE( ONE, KD ), D, 1 )
+*
 *     Compute the GQR factorization of matrices A and B:
 *
 *          Q**T*A = ( R11 ) M,    Q**T*B*Z**T = ( T11   T12 ) M
@@ -356,6 +417,14 @@ SUBROUTINE SGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK,
       CALL SORMRQ( 'Left', 'Transpose', P, 1, NP,
      $             B( MAX( 1, N-P+1 ), 1 ), LDB, WORK( M+1 ), Y,
      $             MAX( 1, P ), WORK( M+NP+1 ), LWORK-M-NP, INFO )
+*
+*     Undo scaling: x carries 2**(KD-KA) and y carries 2**(KD-KB)
+*
+      IF( KA.NE.KD )
+     $   CALL SSCAL( M, SCALE( ONE, KA-KD ), X, 1 )
+      IF( KB.NE.KD )
+     $   CALL SSCAL( P, SCALE( ONE, KB-KD ), Y, 1 )
+*
       WORK( 1 ) = REAL( M + NP + MAX( LOPT, INT( WORK( M+NP+1 ) ) ) )
 *
       RETURN
diff --git a/SRC/sgglse.f b/SRC/sgglse.f
index c14fddb3e..1f2338c57 100644
--- a/SRC/sgglse.f
+++ b/SRC/sgglse.f
@@ -203,26 +203,31 @@ SUBROUTINE SGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK,
 *  =====================================================================
 *
 *     .. Parameters ..
-      REAL               ONE
-      PARAMETER          ( ONE = 1.0E+0 )
+      REAL               ZERO, ONE
+      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 *     ..
 *     .. Local Scalars ..
       LOGICAL            LQUERY
       INTEGER            LOPT, LWKMIN, LWKOPT, MN, NB, NB1, NB2, NB3,
      $                   NB4, NR
+      INTEGER            IA, IB, IBIG, IC, ID, ISML, IU, KA, KB, KT
+      REAL               ANRM, BIGNUM, BNRM, CNRM, DNRM, SMLNUM
+*     ..
+*     .. Local Arrays ..
+      REAL               RWORK( 1 )
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           SAXPY, SCOPY, SGEMV, SGGRQF, SORMQR,
-     $                   SORMRQ,
-     $                   STRMV, STRTRS, XERBLA
+      EXTERNAL           SAXPY, SCOPY, SGEMV, SGGRQF, SLASCL, SORMQR,
+     $                   SORMRQ, SSCAL, STRMV, STRTRS, XERBLA
 *     ..
 *     .. External Functions ..
       INTEGER            ILAENV
       REAL               SROUNDUP_LWORK
-      EXTERNAL           ILAENV, SROUNDUP_LWORK
+      REAL               SLAMCH, SLANGE
+      EXTERNAL           ILAENV, SLAMCH, SLANGE, SROUNDUP_LWORK
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          INT, MAX, MIN
+      INTRINSIC          EXPONENT, HUGE, INT, MAX, MIN, SCALE
 *     ..
 *     .. Executable Statements ..
 *
@@ -277,6 +282,72 @@ SUBROUTINE SGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK,
       IF( N.EQ.0 )
      $   RETURN
 *
+*     Get machine parameters
+*
+      SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' )
+      BIGNUM = ONE / SMLNUM
+      ISML = EXPONENT( SMLNUM )
+      IBIG = EXPONENT( BIGNUM ) - 1
+*
+*     Scale A, B, c and d by powers of two so that their largest
+*     entries lie in [SMLNUM,BIGNUM), the values whose EXPONENT lies
+*     in [ISML,IBIG].  A and B are scaled independently, by 2**KA and
+*     2**KB.  c is scaled by 2**(KA+KT) and d by 2**(KB+KT), which
+*     keeps the constraint and the residual consistent and scales x by
+*     2**KT; KT brings the larger of the two right-hand sides into
+*     range.  Scaling by a power of two is exact.  A norm that is zero,
+*     infinite or NaN takes no part.  SLASCL is called with both
+*     endpoints at or above one so that forming the factor raises no
+*     underflow.
+*
+      ANRM = SLANGE( 'M', M, N, A, LDA, RWORK )
+      KA = 0
+      IF( ANRM.GT.ZERO .AND. ANRM.LE.HUGE( ZERO ) ) THEN
+         IA = EXPONENT( ANRM )
+         IF( IA.LT.ISML ) THEN
+            KA = ISML - IA
+         ELSE IF( IA.GT.IBIG ) THEN
+            KA = IBIG - IA
+         END IF
+      END IF
+      IF( KA.NE.0 )
+     $   CALL SLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KA, 0 ) ),
+     $                SCALE( ONE, MAX( KA, 0 ) ), M, N, A, LDA, INFO )
+*
+      BNRM = SLANGE( 'M', P, N, B, LDB, RWORK )
+      KB = 0
+      IF( BNRM.GT.ZERO .AND. BNRM.LE.HUGE( ZERO ) ) THEN
+         IB = EXPONENT( BNRM )
+         IF( IB.LT.ISML ) THEN
+            KB = ISML - IB
+         ELSE IF( IB.GT.IBIG ) THEN
+            KB = IBIG - IB
+         END IF
+      END IF
+      IF( KB.NE.0 )
+     $   CALL SLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KB, 0 ) ),
+     $                SCALE( ONE, MAX( KB, 0 ) ), P, N, B, LDB, INFO )
+*
+      CNRM = SLANGE( 'M', M, 1, C, MAX( 1, M ), RWORK )
+      DNRM = SLANGE( 'M', P, 1, D, MAX( 1, P ), RWORK )
+      IC = -HUGE( 0 )
+      IF( CNRM.GT.ZERO .AND. CNRM.LE.HUGE( ZERO ) )
+     $   IC = EXPONENT( CNRM ) + KA
+      ID = -HUGE( 0 )
+      IF( DNRM.GT.ZERO .AND. DNRM.LE.HUGE( ZERO ) )
+     $   ID = EXPONENT( DNRM ) + KB
+      IU = MAX( IC, ID )
+      KT = 0
+      IF( IU.GT.IBIG ) THEN
+         KT = IBIG - IU
+      ELSE IF( IU.LT.ISML .AND. IU.NE.-HUGE( 0 ) ) THEN
+         KT = ISML - IU
+      END IF
+      IF( KA+KT.NE.0 )
+     $   CALL SSCAL( M, SCALE( ONE, KA+KT ), C, 1 )
+      IF( KB+KT.NE.0 )
+     $   CALL SSCAL( P, SCALE( ONE, KB+KT ), D, 1 )
+*
 *     Compute the GRQ factorization of matrices B and A:
 *
 *            B*Q**T = (  0  T12 ) P   Z**T*A*Q**T = ( R11 R12 ) N-P
@@ -358,6 +429,15 @@ SUBROUTINE SGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK,
       CALL SORMRQ( 'Left', 'Transpose', N, 1, P, B, LDB, WORK( 1 ),
      $             X,
      $             N, WORK( P+MN+1 ), LWORK-P-MN, INFO )
+*
+*     Undo scaling: x carries 2**KT and the residual in c(N-P+1:M)
+*     carries 2**(KA+KT)
+*
+      IF( KT.NE.0 )
+     $   CALL SSCAL( N, SCALE( ONE, -KT ), X, 1 )
+      IF( KA+KT.NE.0 .AND. M+P.GT.N )
+     $   CALL SSCAL( M+P-N, SCALE( ONE, -KA-KT ), C( N-P+1 ), 1 )
+*
       WORK( 1 ) = REAL( P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) ) )
 *
       RETURN
diff --git a/SRC/zggglm.f b/SRC/zggglm.f
index 3a0d550ab..6bced0748 100644
--- a/SRC/zggglm.f
+++ b/SRC/zggglm.f
@@ -208,6 +208,8 @@ SUBROUTINE ZGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK,
 *  ===================================================================
 *
 *     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
       COMPLEX*16         CZERO, CONE
       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
@@ -216,18 +218,23 @@ SUBROUTINE ZGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK,
       LOGICAL            LQUERY
       INTEGER            I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3,
      $                   NB4, NP
+      INTEGER            IA, IB, IBIG, ID, ISML, KA, KB, KD
+      DOUBLE PRECISION   ANRM, BIGNUM, BNRM, DNRM, SMLNUM
+*     ..
+*     .. Local Arrays ..
+      DOUBLE PRECISION   RWORK( 1 )
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           XERBLA, ZCOPY, ZGEMV, ZGGQRF, ZTRTRS,
-     $                   ZUNMQR,
-     $                   ZUNMRQ
+      EXTERNAL           XERBLA, ZCOPY, ZDSCAL, ZGEMV, ZGGQRF, ZLASCL,
+     $                   ZTRTRS, ZUNMQR, ZUNMRQ
 *     ..
 *     .. External Functions ..
       INTEGER            ILAENV
-      EXTERNAL           ILAENV
+      DOUBLE PRECISION   DLAMCH, ZLANGE
+      EXTERNAL           DLAMCH, ILAENV, ZLANGE
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          INT, MAX, MIN
+      INTRINSIC          EXPONENT, HUGE, INT, MAX, MIN, SCALE
 *     ..
 *     .. Executable Statements ..
 *
@@ -289,6 +296,62 @@ SUBROUTINE ZGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK,
          RETURN
       END IF
 *
+*     Get machine parameters
+*
+      SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
+      BIGNUM = ONE / SMLNUM
+      ISML = EXPONENT( SMLNUM )
+      IBIG = EXPONENT( BIGNUM ) - 1
+*
+*     Scale A, B and d independently by powers of two 2**KA, 2**KB and
+*     2**KD so that their largest entries lie in [SMLNUM,BIGNUM), the
+*     values whose EXPONENT lies in [ISML,IBIG].  x is then scaled by
+*     2**(KD-KA) and y by 2**(KD-KB).  Scaling by a power of two is
+*     exact.  A norm that is zero, infinite or NaN takes no part.
+*     ZLASCL is called with both endpoints at or above one so that
+*     forming the factor raises no underflow.
+*
+      ANRM = ZLANGE( 'M', N, M, A, LDA, RWORK )
+      KA = 0
+      IF( ANRM.GT.ZERO .AND. ANRM.LE.HUGE( ZERO ) ) THEN
+         IA = EXPONENT( ANRM )
+         IF( IA.LT.ISML ) THEN
+            KA = ISML - IA
+         ELSE IF( IA.GT.IBIG ) THEN
+            KA = IBIG - IA
+         END IF
+      END IF
+      IF( KA.NE.0 )
+     $   CALL ZLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KA, 0 ) ),
+     $                SCALE( ONE, MAX( KA, 0 ) ), N, M, A, LDA, INFO )
+*
+      BNRM = ZLANGE( 'M', N, P, B, LDB, RWORK )
+      KB = 0
+      IF( BNRM.GT.ZERO .AND. BNRM.LE.HUGE( ZERO ) ) THEN
+         IB = EXPONENT( BNRM )
+         IF( IB.LT.ISML ) THEN
+            KB = ISML - IB
+         ELSE IF( IB.GT.IBIG ) THEN
+            KB = IBIG - IB
+         END IF
+      END IF
+      IF( KB.NE.0 )
+     $   CALL ZLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KB, 0 ) ),
+     $                SCALE( ONE, MAX( KB, 0 ) ), N, P, B, LDB, INFO )
+*
+      DNRM = ZLANGE( 'M', N, 1, D, N, RWORK )
+      KD = 0
+      IF( DNRM.GT.ZERO .AND. DNRM.LE.HUGE( ZERO ) ) THEN
+         ID = EXPONENT( DNRM )
+         IF( ID.LT.ISML ) THEN
+            KD = ISML - ID
+         ELSE IF( ID.GT.IBIG ) THEN
+            KD = IBIG - ID
+         END IF
+      END IF
+      IF( KD.NE.0 )
+     $   CALL ZDSCAL( N, SCALE( ONE, KD ), D, 1 )
+*
 *     Compute the GQR factorization of matrices A and B:
 *
 *          Q**H*A = ( R11 ) M,    Q**H*B*Z**H = ( T11   T12 ) M
@@ -358,6 +421,14 @@ SUBROUTINE ZGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK,
       CALL ZUNMRQ( 'Left', 'Conjugate transpose', P, 1, NP,
      $             B( MAX( 1, N-P+1 ), 1 ), LDB, WORK( M+1 ), Y,
      $             MAX( 1, P ), WORK( M+NP+1 ), LWORK-M-NP, INFO )
+*
+*     Undo scaling: x carries 2**(KD-KA) and y carries 2**(KD-KB)
+*
+      IF( KA.NE.KD )
+     $   CALL ZDSCAL( M, SCALE( ONE, KA-KD ), X, 1 )
+      IF( KB.NE.KD )
+     $   CALL ZDSCAL( P, SCALE( ONE, KB-KD ), Y, 1 )
+*
       WORK( 1 ) = M + NP + MAX( LOPT, INT( WORK( M+NP+1 ) ) )
 *
       RETURN
diff --git a/SRC/zgglse.f b/SRC/zgglse.f
index 00ab2cfab..1966061ca 100644
--- a/SRC/zgglse.f
+++ b/SRC/zgglse.f
@@ -203,6 +203,8 @@ SUBROUTINE ZGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK,
 *  =====================================================================
 *
 *     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
       COMPLEX*16         CONE
       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
 *     ..
@@ -210,18 +212,23 @@ SUBROUTINE ZGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK,
       LOGICAL            LQUERY
       INTEGER            LOPT, LWKMIN, LWKOPT, MN, NB, NB1, NB2, NB3,
      $                   NB4, NR
+      INTEGER            IA, IB, IBIG, IC, ID, ISML, IU, KA, KB, KT
+      DOUBLE PRECISION   ANRM, BIGNUM, BNRM, CNRM, DNRM, SMLNUM
+*     ..
+*     .. Local Arrays ..
+      DOUBLE PRECISION   RWORK( 1 )
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           XERBLA, ZAXPY, ZCOPY, ZGEMV, ZGGRQF,
-     $                   ZTRMV,
-     $                   ZTRTRS, ZUNMQR, ZUNMRQ
+      EXTERNAL           XERBLA, ZAXPY, ZCOPY, ZDSCAL, ZGEMV, ZGGRQF,
+     $                   ZLASCL, ZTRMV, ZTRTRS, ZUNMQR, ZUNMRQ
 *     ..
 *     .. External Functions ..
       INTEGER            ILAENV
-      EXTERNAL           ILAENV
+      DOUBLE PRECISION   DLAMCH, ZLANGE
+      EXTERNAL           DLAMCH, ILAENV, ZLANGE
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          INT, MAX, MIN
+      INTRINSIC          EXPONENT, HUGE, INT, MAX, MIN, SCALE
 *     ..
 *     .. Executable Statements ..
 *
@@ -276,6 +283,72 @@ SUBROUTINE ZGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK,
       IF( N.EQ.0 )
      $   RETURN
 *
+*     Get machine parameters
+*
+      SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
+      BIGNUM = ONE / SMLNUM
+      ISML = EXPONENT( SMLNUM )
+      IBIG = EXPONENT( BIGNUM ) - 1
+*
+*     Scale A, B, c and d by powers of two so that their largest
+*     entries lie in [SMLNUM,BIGNUM), the values whose EXPONENT lies
+*     in [ISML,IBIG].  A and B are scaled independently, by 2**KA and
+*     2**KB.  c is scaled by 2**(KA+KT) and d by 2**(KB+KT), which
+*     keeps the constraint and the residual consistent and scales x by
+*     2**KT; KT brings the larger of the two right-hand sides into
+*     range.  Scaling by a power of two is exact.  A norm that is zero,
+*     infinite or NaN takes no part.  ZLASCL is called with both
+*     endpoints at or above one so that forming the factor raises no
+*     underflow.
+*
+      ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK )
+      KA = 0
+      IF( ANRM.GT.ZERO .AND. ANRM.LE.HUGE( ZERO ) ) THEN
+         IA = EXPONENT( ANRM )
+         IF( IA.LT.ISML ) THEN
+            KA = ISML - IA
+         ELSE IF( IA.GT.IBIG ) THEN
+            KA = IBIG - IA
+         END IF
+      END IF
+      IF( KA.NE.0 )
+     $   CALL ZLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KA, 0 ) ),
+     $                SCALE( ONE, MAX( KA, 0 ) ), M, N, A, LDA, INFO )
+*
+      BNRM = ZLANGE( 'M', P, N, B, LDB, RWORK )
+      KB = 0
+      IF( BNRM.GT.ZERO .AND. BNRM.LE.HUGE( ZERO ) ) THEN
+         IB = EXPONENT( BNRM )
+         IF( IB.LT.ISML ) THEN
+            KB = ISML - IB
+         ELSE IF( IB.GT.IBIG ) THEN
+            KB = IBIG - IB
+         END IF
+      END IF
+      IF( KB.NE.0 )
+     $   CALL ZLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KB, 0 ) ),
+     $                SCALE( ONE, MAX( KB, 0 ) ), P, N, B, LDB, INFO )
+*
+      CNRM = ZLANGE( 'M', M, 1, C, MAX( 1, M ), RWORK )
+      DNRM = ZLANGE( 'M', P, 1, D, MAX( 1, P ), RWORK )
+      IC = -HUGE( 0 )
+      IF( CNRM.GT.ZERO .AND. CNRM.LE.HUGE( ZERO ) )
+     $   IC = EXPONENT( CNRM ) + KA
+      ID = -HUGE( 0 )
+      IF( DNRM.GT.ZERO .AND. DNRM.LE.HUGE( ZERO ) )
+     $   ID = EXPONENT( DNRM ) + KB
+      IU = MAX( IC, ID )
+      KT = 0
+      IF( IU.GT.IBIG ) THEN
+         KT = IBIG - IU
+      ELSE IF( IU.LT.ISML .AND. IU.NE.-HUGE( 0 ) ) THEN
+         KT = ISML - IU
+      END IF
+      IF( KA+KT.NE.0 )
+     $   CALL ZDSCAL( M, SCALE( ONE, KA+KT ), C, 1 )
+      IF( KB+KT.NE.0 )
+     $   CALL ZDSCAL( P, SCALE( ONE, KB+KT ), D, 1 )
+*
 *     Compute the GRQ factorization of matrices B and A:
 *
 *            B*Q**H = (  0  T12 ) P   Z**H*A*Q**H = ( R11 R12 ) N-P
@@ -356,6 +429,15 @@ SUBROUTINE ZGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK,
 *
       CALL ZUNMRQ( 'Left', 'Conjugate Transpose', N, 1, P, B, LDB,
      $             WORK( 1 ), X, N, WORK( P+MN+1 ), LWORK-P-MN, INFO )
+*
+*     Undo scaling: x carries 2**KT and the residual in c(N-P+1:M)
+*     carries 2**(KA+KT)
+*
+      IF( KT.NE.0 )
+     $   CALL ZDSCAL( N, SCALE( ONE, -KT ), X, 1 )
+      IF( KA+KT.NE.0 .AND. M+P.GT.N )
+     $   CALL ZDSCAL( M+P-N, SCALE( ONE, -KA-KT ), C( N-P+1 ), 1 )
+*
       WORK( 1 ) = P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) )
 *
       RETURN
diff --git a/TESTING/EIG/cglmts.f b/TESTING/EIG/cglmts.f
index 9bd27abc1..8445bcb32 100644
--- a/TESTING/EIG/cglmts.f
+++ b/TESTING/EIG/cglmts.f
@@ -28,7 +28,7 @@
 *> \verbatim
 *>
 *> CGLMTS tests CGGGLM - a subroutine for solving the generalized
-*> linear model problem.
+*> linear model problem, including independent scaling of its inputs.
 *> \endverbatim
 *
 *  Arguments:
@@ -168,22 +168,28 @@ SUBROUTINE CGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF,
 *     .. Parameters ..
       REAL               ZERO
       PARAMETER          ( ZERO = 0.0E+0 )
+      REAL              ONE
+      PARAMETER          ( ONE = 1.0E+0 )
+      INTEGER            MAXEXP
+      PARAMETER          ( MAXEXP = MAXEXPONENT( ZERO ) - 2 )
       COMPLEX            CONE
       PARAMETER          ( CONE = 1.0E+0 )
 *     ..
 *     .. Local Scalars ..
-      INTEGER            INFO
+      INTEGER            INFO, ISCALE, J
       REAL               ANORM, BNORM, EPS, XNORM, YNORM, DNORM, UNFL
+      REAL               ASCL, BSCL, DSCL, SCL, TNRM
 *     ..
 *     .. External Functions ..
       REAL               SCASUM, SLAMCH, CLANGE
-      EXTERNAL           SCASUM, SLAMCH, CLANGE
+      LOGICAL            SISNAN
+      EXTERNAL           SISNAN, SCASUM, SLAMCH, CLANGE
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           CCOPY, CGEMV, CGGGLM, CLACPY
+      EXTERNAL           CCOPY, CGEMV, CGGGLM, CLACPY, CSSCAL
 *
 *     .. Intrinsic Functions ..
-      INTRINSIC          MAX
+      INTRINSIC          EXPONENT, MAX, MAXEXPONENT, SCALE
 *     ..
 *     .. Executable Statements ..
 *
@@ -226,6 +232,68 @@ SUBROUTINE CGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF,
       ELSE
          RESULT =  ( ( DNORM / YNORM ) / XNORM ) /EPS
       END IF
+*
+*     A, B and d may be scaled independently: for factors a, b, d,
+*     the solutions become (d/a)*x and (d/b)*u.  Test large and
+*     tiny inputs, then undo solution scaling before the residual.
+*
+      TNRM = MAX( CLANGE( 'M', N, M, A, LDA, RWORK ),
+     $           CLANGE( 'M', N, P, B, LDB, RWORK ),
+     $           CLANGE( 'M', N, 1, D, N, RWORK ) )
+      IF( TNRM.GT.ZERO .AND. TNRM.LE.SLAMCH( 'Overflow' ) ) THEN
+*
+*        Cases: common large, common tiny, tiny d, tiny (A,d),
+*        and tiny (B,d).
+*
+         DO 30 ISCALE = 1, 5
+            IF( ISCALE.EQ.1 ) THEN
+               SCL = SCALE( ONE, MAXEXP-EXPONENT( TNRM ) )
+            ELSE
+               SCL = SLAMCH( 'Safe minimum' ) /
+     $               SLAMCH( 'Precision' )
+               SCL = SCALE( ONE, EXPONENT( SCL )-4-
+     $                      EXPONENT( TNRM ) )
+            END IF
+            ASCL = SCL
+            BSCL = SCL
+            DSCL = SCL
+            IF( ISCALE.EQ.3 .OR. ISCALE.EQ.5 ) ASCL = ONE
+            IF( ISCALE.EQ.3 .OR. ISCALE.EQ.4 ) BSCL = ONE
+            CALL CLACPY( 'Full', N, M, A, LDA, AF, LDA )
+            CALL CLACPY( 'Full', N, P, B, LDB, BF, LDB )
+            CALL CCOPY( N, D, 1, DF, 1 )
+            DO 10 J = 1, M
+               CALL CSSCAL( N, ASCL, AF( 1, J ), 1 )
+   10       CONTINUE
+            DO 20 J = 1, P
+               CALL CSSCAL( N, BSCL, BF( 1, J ), 1 )
+   20       CONTINUE
+            CALL CSSCAL( N, DSCL, DF, 1 )
+*
+            CALL CGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK,
+     $                   LWORK, INFO )
+            IF( INFO.NE.0 ) THEN
+               RESULT = ONE / EPS
+               GO TO 30
+            END IF
+            CALL CSSCAL( M, ASCL / DSCL, X, 1 )
+            CALL CSSCAL( P, BSCL / DSCL, U, 1 )
+*
+            CALL CCOPY( N, D, 1, DF, 1 )
+            CALL CGEMV( 'No transpose', N, M, -CONE, A, LDA, X, 1, CONE,
+     $                  DF, 1 )
+            CALL CGEMV( 'No transpose', N, P, -CONE, B, LDB, U, 1, CONE,
+     $                  DF, 1 )
+            DNORM = SCASUM( N, DF, 1 )
+            XNORM = SCASUM( M, X, 1 ) + SCASUM( P, U, 1 )
+            IF( SISNAN( DNORM ) .OR. SISNAN( XNORM ) ) THEN
+               RESULT = ONE / EPS
+            ELSE IF( XNORM.GT.ZERO ) THEN
+               RESULT = MAX( RESULT,
+     $                       ( ( DNORM / YNORM ) / XNORM ) / EPS )
+            END IF
+   30    CONTINUE
+      END IF
 *
       RETURN
 *
diff --git a/TESTING/EIG/clsets.f b/TESTING/EIG/clsets.f
index 7e940230e..fad4a0169 100644
--- a/TESTING/EIG/clsets.f
+++ b/TESTING/EIG/clsets.f
@@ -27,7 +27,7 @@
 *> \verbatim
 *>
 *> CLSETS tests CGGLSE - a subroutine for solving linear equality
-*> constrained least square problem (LSE).
+*> constrained least square problem (LSE), including scaled inputs.
 *> \endverbatim
 *
 *  Arguments:
@@ -171,10 +171,25 @@ SUBROUTINE CLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF,
 *
 *     ..
 *     .. Local Scalars ..
-      INTEGER            INFO
+      INTEGER            INFO, ISCALE, J
+      REAL              ASCL, BSCL, CSCL, DSCL, RESID, SCL, TNRM
+*     ..
+*     .. Parameters ..
+      REAL              ZERO, ONE
+      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
+      INTEGER            MAXEXP
+      PARAMETER          ( MAXEXP = MAXEXPONENT( ZERO ) - 2 )
+*     ..
+*     .. External Functions ..
+      LOGICAL            SISNAN
+      REAL              CLANGE, SLAMCH
+      EXTERNAL           SISNAN, CLANGE, SLAMCH
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           CCOPY, CGGLSE, CLACPY, CGET02
+      EXTERNAL           CCOPY, CGGLSE, CLACPY, CGET02, CSSCAL
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          EXPONENT, MAX, MAXEXPONENT, SCALE
 *     ..
 *     .. Executable Statements ..
 *
@@ -204,6 +219,76 @@ SUBROUTINE CLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF,
 *
       CALL CGET02( 'No transpose', P, N, 1, B, LDB, X, N, DF, P,
      $             RWORK, RESULT( 2 ) )
+*
+*     Scaling (A,c) and (B,d) independently leaves x unchanged.
+*     Scaling only (c,d) scales x by the same factor.  Exercise
+*     both ends of the range and check residuals at the input scale.
+*
+      TNRM = MAX( CLANGE( 'M', M, N, A, LDA, RWORK ),
+     $           CLANGE( 'M', P, N, B, LDB, RWORK ),
+     $           CLANGE( 'M', M, 1, C, M, RWORK ),
+     $           CLANGE( 'M', P, 1, D, P, RWORK ) )
+      IF( TNRM.GT.ZERO .AND. TNRM.LE.SLAMCH( 'Overflow' ) ) THEN
+*
+*        Cases: common large, common tiny, tiny (c,d), tiny (A,c),
+*        and tiny (B,d).
+*
+         DO 30 ISCALE = 1, 5
+            IF( ISCALE.EQ.1 ) THEN
+               SCL = SCALE( ONE, MAXEXP-EXPONENT( TNRM ) )
+            ELSE
+               SCL = SLAMCH( 'Safe minimum' ) /
+     $               SLAMCH( 'Precision' )
+               SCL = SCALE( ONE, EXPONENT( SCL )-4-
+     $                      EXPONENT( TNRM ) )
+            END IF
+            ASCL = SCL
+            BSCL = SCL
+            CSCL = SCL
+            DSCL = SCL
+            IF( ISCALE.EQ.3 .OR. ISCALE.EQ.5 ) ASCL = ONE
+            IF( ISCALE.EQ.3 .OR. ISCALE.EQ.4 ) BSCL = ONE
+            IF( ISCALE.EQ.4 ) DSCL = ONE
+            IF( ISCALE.EQ.5 ) CSCL = ONE
+            CALL CLACPY( 'Full', M, N, A, LDA, AF, LDA )
+            CALL CLACPY( 'Full', P, N, B, LDB, BF, LDB )
+            CALL CCOPY( M, C, 1, CF, 1 )
+            CALL CCOPY( P, D, 1, DF, 1 )
+            DO 10 J = 1, N
+               CALL CSSCAL( M, ASCL, AF( 1, J ), 1 )
+               CALL CSSCAL( P, BSCL, BF( 1, J ), 1 )
+   10       CONTINUE
+            CALL CSSCAL( M, CSCL, CF, 1 )
+            CALL CSSCAL( P, DSCL, DF, 1 )
+*
+            CALL CGGLSE( M, N, P, AF, LDA, BF, LDB, CF, DF, X, WORK,
+     $                   LWORK, INFO )
+            IF( INFO.NE.0 ) THEN
+               RESULT( 1 ) = ONE / SLAMCH( 'Epsilon' )
+               RESULT( 2 ) = RESULT( 1 )
+               GO TO 30
+            END IF
+            IF( ISCALE.EQ.3 )
+     $         CALL CSSCAL( N, ONE / SCL, X, 1 )
+*
+            CALL CCOPY( M, C, 1, CF, 1 )
+            CALL CCOPY( P, D, 1, DF, 1 )
+            CALL CGET02( 'No transpose', M, N, 1, A, LDA, X, N, CF, M,
+     $                   RWORK, RESID )
+            IF( SISNAN( RESID ) ) THEN
+               RESULT( 1 ) = ONE / SLAMCH( 'Epsilon' )
+            ELSE
+               RESULT( 1 ) = MAX( RESULT( 1 ), RESID )
+            END IF
+            CALL CGET02( 'No transpose', P, N, 1, B, LDB, X, N, DF, P,
+     $                   RWORK, RESID )
+            IF( SISNAN( RESID ) ) THEN
+               RESULT( 2 ) = ONE / SLAMCH( 'Epsilon' )
+            ELSE
+               RESULT( 2 ) = MAX( RESULT( 2 ), RESID )
+            END IF
+   30    CONTINUE
+      END IF
 *
       RETURN
 *
diff --git a/TESTING/EIG/dglmts.f b/TESTING/EIG/dglmts.f
index d2d24fea3..03be4917a 100644
--- a/TESTING/EIG/dglmts.f
+++ b/TESTING/EIG/dglmts.f
@@ -24,7 +24,7 @@
 *> \verbatim
 *>
 *> DGLMTS tests DGGGLM - a subroutine for solving the generalized
-*> linear model problem.
+*> linear model problem, including independent scaling of its inputs.
 *> \endverbatim
 *
 *  Arguments:
@@ -164,21 +164,25 @@ SUBROUTINE DGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U,
 *     .. Parameters ..
       DOUBLE PRECISION   ZERO, ONE
       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
+      INTEGER            MAXEXP
+      PARAMETER          ( MAXEXP = MAXEXPONENT( ZERO ) - 2 )
 *     ..
 *     .. Local Scalars ..
-      INTEGER            INFO
+      INTEGER            INFO, ISCALE, J
       DOUBLE PRECISION   ANORM, BNORM, DNORM, EPS, UNFL, XNORM, YNORM
+      DOUBLE PRECISION   ASCL, BSCL, DSCL, SCL, TNRM
 *     ..
 *     .. External Functions ..
       DOUBLE PRECISION   DASUM, DLAMCH, DLANGE
-      EXTERNAL           DASUM, DLAMCH, DLANGE
+      LOGICAL            DISNAN
+      EXTERNAL           DISNAN, DASUM, DLAMCH, DLANGE
 *     ..
 *     .. External Subroutines ..
 *
-      EXTERNAL           DCOPY, DGEMV, DGGGLM, DLACPY
+      EXTERNAL           DCOPY, DGEMV, DGGGLM, DLACPY, DSCAL
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          MAX
+      INTRINSIC          EXPONENT, MAX, MAXEXPONENT, SCALE
 *     ..
 *     .. Executable Statements ..
 *
@@ -219,6 +223,68 @@ SUBROUTINE DGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U,
       ELSE
          RESULT = ( ( DNORM / YNORM ) / XNORM ) / EPS
       END IF
+*
+*     A, B and d may be scaled independently: for factors a, b, d,
+*     the solutions become (d/a)*x and (d/b)*u.  Test large and
+*     tiny inputs, then undo solution scaling before the residual.
+*
+      TNRM = MAX( DLANGE( 'M', N, M, A, LDA, RWORK ),
+     $           DLANGE( 'M', N, P, B, LDB, RWORK ),
+     $           DLANGE( 'M', N, 1, D, N, RWORK ) )
+      IF( TNRM.GT.ZERO .AND. TNRM.LE.DLAMCH( 'Overflow' ) ) THEN
+*
+*        Cases: common large, common tiny, tiny d, tiny (A,d),
+*        and tiny (B,d).
+*
+         DO 30 ISCALE = 1, 5
+            IF( ISCALE.EQ.1 ) THEN
+               SCL = SCALE( ONE, MAXEXP-EXPONENT( TNRM ) )
+            ELSE
+               SCL = DLAMCH( 'Safe minimum' ) /
+     $               DLAMCH( 'Precision' )
+               SCL = SCALE( ONE, EXPONENT( SCL )-4-
+     $                      EXPONENT( TNRM ) )
+            END IF
+            ASCL = SCL
+            BSCL = SCL
+            DSCL = SCL
+            IF( ISCALE.EQ.3 .OR. ISCALE.EQ.5 ) ASCL = ONE
+            IF( ISCALE.EQ.3 .OR. ISCALE.EQ.4 ) BSCL = ONE
+            CALL DLACPY( 'Full', N, M, A, LDA, AF, LDA )
+            CALL DLACPY( 'Full', N, P, B, LDB, BF, LDB )
+            CALL DCOPY( N, D, 1, DF, 1 )
+            DO 10 J = 1, M
+               CALL DSCAL( N, ASCL, AF( 1, J ), 1 )
+   10       CONTINUE
+            DO 20 J = 1, P
+               CALL DSCAL( N, BSCL, BF( 1, J ), 1 )
+   20       CONTINUE
+            CALL DSCAL( N, DSCL, DF, 1 )
+*
+            CALL DGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK,
+     $                   LWORK, INFO )
+            IF( INFO.NE.0 ) THEN
+               RESULT = ONE / EPS
+               GO TO 30
+            END IF
+            CALL DSCAL( M, ASCL / DSCL, X, 1 )
+            CALL DSCAL( P, BSCL / DSCL, U, 1 )
+*
+            CALL DCOPY( N, D, 1, DF, 1 )
+            CALL DGEMV( 'No transpose', N, M, -ONE, A, LDA, X, 1, ONE,
+     $                  DF, 1 )
+            CALL DGEMV( 'No transpose', N, P, -ONE, B, LDB, U, 1, ONE,
+     $                  DF, 1 )
+            DNORM = DASUM( N, DF, 1 )
+            XNORM = DASUM( M, X, 1 ) + DASUM( P, U, 1 )
+            IF( DISNAN( DNORM ) .OR. DISNAN( XNORM ) ) THEN
+               RESULT = ONE / EPS
+            ELSE IF( XNORM.GT.ZERO ) THEN
+               RESULT = MAX( RESULT,
+     $                       ( ( DNORM / YNORM ) / XNORM ) / EPS )
+            END IF
+   30    CONTINUE
+      END IF
 *
       RETURN
 *
diff --git a/TESTING/EIG/dlsets.f b/TESTING/EIG/dlsets.f
index afb46a727..4ec6926b8 100644
--- a/TESTING/EIG/dlsets.f
+++ b/TESTING/EIG/dlsets.f
@@ -23,7 +23,7 @@
 *> \verbatim
 *>
 *> DLSETS tests DGGLSE - a subroutine for solving linear equality
-*> constrained least square problem (LSE).
+*> constrained least square problem (LSE), including scaled inputs.
 *> \endverbatim
 *
 *  Arguments:
@@ -166,10 +166,25 @@ SUBROUTINE DLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF, D, DF,
      $                   RESULT( 2 ), RWORK( * ), WORK( LWORK ), X( * )
 *     ..
 *     .. Local Scalars ..
-      INTEGER            INFO
+      INTEGER            INFO, ISCALE, J
+      DOUBLE PRECISION   ASCL, BSCL, CSCL, DSCL, RESID, SCL, TNRM
+*     ..
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
+      INTEGER            MAXEXP
+      PARAMETER          ( MAXEXP = MAXEXPONENT( ZERO ) - 2 )
+*     ..
+*     .. External Functions ..
+      LOGICAL            DISNAN
+      DOUBLE PRECISION   DLANGE, DLAMCH
+      EXTERNAL           DISNAN, DLANGE, DLAMCH
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           DCOPY, DGET02, DGGLSE, DLACPY
+      EXTERNAL           DCOPY, DGET02, DGGLSE, DLACPY, DSCAL
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          EXPONENT, MAX, MAXEXPONENT, SCALE
 *     ..
 *     .. Executable Statements ..
 *
@@ -199,6 +214,76 @@ SUBROUTINE DLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF, D, DF,
 *
       CALL DGET02( 'No transpose', P, N, 1, B, LDB, X, N, DF, P, RWORK,
      $             RESULT( 2 ) )
+*
+*     Scaling (A,c) and (B,d) independently leaves x unchanged.
+*     Scaling only (c,d) scales x by the same factor.  Exercise
+*     both ends of the range and check residuals at the input scale.
+*
+      TNRM = MAX( DLANGE( 'M', M, N, A, LDA, RWORK ),
+     $           DLANGE( 'M', P, N, B, LDB, RWORK ),
+     $           DLANGE( 'M', M, 1, C, M, RWORK ),
+     $           DLANGE( 'M', P, 1, D, P, RWORK ) )
+      IF( TNRM.GT.ZERO .AND. TNRM.LE.DLAMCH( 'Overflow' ) ) THEN
+*
+*        Cases: common large, common tiny, tiny (c,d), tiny (A,c),
+*        and tiny (B,d).
+*
+         DO 30 ISCALE = 1, 5
+            IF( ISCALE.EQ.1 ) THEN
+               SCL = SCALE( ONE, MAXEXP-EXPONENT( TNRM ) )
+            ELSE
+               SCL = DLAMCH( 'Safe minimum' ) /
+     $               DLAMCH( 'Precision' )
+               SCL = SCALE( ONE, EXPONENT( SCL )-4-
+     $                      EXPONENT( TNRM ) )
+            END IF
+            ASCL = SCL
+            BSCL = SCL
+            CSCL = SCL
+            DSCL = SCL
+            IF( ISCALE.EQ.3 .OR. ISCALE.EQ.5 ) ASCL = ONE
+            IF( ISCALE.EQ.3 .OR. ISCALE.EQ.4 ) BSCL = ONE
+            IF( ISCALE.EQ.4 ) DSCL = ONE
+            IF( ISCALE.EQ.5 ) CSCL = ONE
+            CALL DLACPY( 'Full', M, N, A, LDA, AF, LDA )
+            CALL DLACPY( 'Full', P, N, B, LDB, BF, LDB )
+            CALL DCOPY( M, C, 1, CF, 1 )
+            CALL DCOPY( P, D, 1, DF, 1 )
+            DO 10 J = 1, N
+               CALL DSCAL( M, ASCL, AF( 1, J ), 1 )
+               CALL DSCAL( P, BSCL, BF( 1, J ), 1 )
+   10       CONTINUE
+            CALL DSCAL( M, CSCL, CF, 1 )
+            CALL DSCAL( P, DSCL, DF, 1 )
+*
+            CALL DGGLSE( M, N, P, AF, LDA, BF, LDB, CF, DF, X, WORK,
+     $                   LWORK, INFO )
+            IF( INFO.NE.0 ) THEN
+               RESULT( 1 ) = ONE / DLAMCH( 'Epsilon' )
+               RESULT( 2 ) = RESULT( 1 )
+               GO TO 30
+            END IF
+            IF( ISCALE.EQ.3 )
+     $         CALL DSCAL( N, ONE / SCL, X, 1 )
+*
+            CALL DCOPY( M, C, 1, CF, 1 )
+            CALL DCOPY( P, D, 1, DF, 1 )
+            CALL DGET02( 'No transpose', M, N, 1, A, LDA, X, N, CF, M,
+     $                   RWORK, RESID )
+            IF( DISNAN( RESID ) ) THEN
+               RESULT( 1 ) = ONE / DLAMCH( 'Epsilon' )
+            ELSE
+               RESULT( 1 ) = MAX( RESULT( 1 ), RESID )
+            END IF
+            CALL DGET02( 'No transpose', P, N, 1, B, LDB, X, N, DF, P,
+     $                   RWORK, RESID )
+            IF( DISNAN( RESID ) ) THEN
+               RESULT( 2 ) = ONE / DLAMCH( 'Epsilon' )
+            ELSE
+               RESULT( 2 ) = MAX( RESULT( 2 ), RESID )
+            END IF
+   30    CONTINUE
+      END IF
 *
       RETURN
 *
diff --git a/TESTING/EIG/sglmts.f b/TESTING/EIG/sglmts.f
index b4b3e228d..176674530 100644
--- a/TESTING/EIG/sglmts.f
+++ b/TESTING/EIG/sglmts.f
@@ -27,7 +27,7 @@
 *> \verbatim
 *>
 *> SGLMTS tests SGGGLM - a subroutine for solving the generalized
-*> linear model problem.
+*> linear model problem, including independent scaling of its inputs.
 *> \endverbatim
 *
 *  Arguments:
@@ -166,20 +166,24 @@ SUBROUTINE SGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF,
 *     .. Parameters ..
       REAL               ZERO, ONE
       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
+      INTEGER            MAXEXP
+      PARAMETER          ( MAXEXP = MAXEXPONENT( ZERO ) - 2 )
 *     ..
 *     .. Local Scalars ..
-      INTEGER            INFO
+      INTEGER            INFO, ISCALE, J
       REAL               ANORM, BNORM, EPS, XNORM, YNORM, DNORM, UNFL
+      REAL               ASCL, BSCL, DSCL, SCL, TNRM
 *     ..
 *     .. External Functions ..
       REAL               SASUM, SLAMCH, SLANGE
-      EXTERNAL           SASUM, SLAMCH, SLANGE
+      LOGICAL            SISNAN
+      EXTERNAL           SISNAN, SASUM, SLAMCH, SLANGE
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           SCOPY, SGEMV, SGGGLM, SLACPY
+      EXTERNAL           SCOPY, SGEMV, SGGGLM, SLACPY, SSCAL
 *
 *     .. Intrinsic Functions ..
-      INTRINSIC          MAX
+      INTRINSIC          EXPONENT, MAX, MAXEXPONENT, SCALE
 *     ..
 *     .. Executable Statements ..
 *
@@ -222,6 +226,68 @@ SUBROUTINE SGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF,
       ELSE
          RESULT =  ( ( DNORM / YNORM ) / XNORM ) /EPS
       END IF
+*
+*     A, B and d may be scaled independently: for factors a, b, d,
+*     the solutions become (d/a)*x and (d/b)*u.  Test large and
+*     tiny inputs, then undo solution scaling before the residual.
+*
+      TNRM = MAX( SLANGE( 'M', N, M, A, LDA, RWORK ),
+     $           SLANGE( 'M', N, P, B, LDB, RWORK ),
+     $           SLANGE( 'M', N, 1, D, N, RWORK ) )
+      IF( TNRM.GT.ZERO .AND. TNRM.LE.SLAMCH( 'Overflow' ) ) THEN
+*
+*        Cases: common large, common tiny, tiny d, tiny (A,d),
+*        and tiny (B,d).
+*
+         DO 30 ISCALE = 1, 5
+            IF( ISCALE.EQ.1 ) THEN
+               SCL = SCALE( ONE, MAXEXP-EXPONENT( TNRM ) )
+            ELSE
+               SCL = SLAMCH( 'Safe minimum' ) /
+     $               SLAMCH( 'Precision' )
+               SCL = SCALE( ONE, EXPONENT( SCL )-4-
+     $                      EXPONENT( TNRM ) )
+            END IF
+            ASCL = SCL
+            BSCL = SCL
+            DSCL = SCL
+            IF( ISCALE.EQ.3 .OR. ISCALE.EQ.5 ) ASCL = ONE
+            IF( ISCALE.EQ.3 .OR. ISCALE.EQ.4 ) BSCL = ONE
+            CALL SLACPY( 'Full', N, M, A, LDA, AF, LDA )
+            CALL SLACPY( 'Full', N, P, B, LDB, BF, LDB )
+            CALL SCOPY( N, D, 1, DF, 1 )
+            DO 10 J = 1, M
+               CALL SSCAL( N, ASCL, AF( 1, J ), 1 )
+   10       CONTINUE
+            DO 20 J = 1, P
+               CALL SSCAL( N, BSCL, BF( 1, J ), 1 )
+   20       CONTINUE
+            CALL SSCAL( N, DSCL, DF, 1 )
+*
+            CALL SGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK,
+     $                   LWORK, INFO )
+            IF( INFO.NE.0 ) THEN
+               RESULT = ONE / EPS
+               GO TO 30
+            END IF
+            CALL SSCAL( M, ASCL / DSCL, X, 1 )
+            CALL SSCAL( P, BSCL / DSCL, U, 1 )
+*
+            CALL SCOPY( N, D, 1, DF, 1 )
+            CALL SGEMV( 'No transpose', N, M, -ONE, A, LDA, X, 1, ONE,
+     $                  DF, 1 )
+            CALL SGEMV( 'No transpose', N, P, -ONE, B, LDB, U, 1, ONE,
+     $                  DF, 1 )
+            DNORM = SASUM( N, DF, 1 )
+            XNORM = SASUM( M, X, 1 ) + SASUM( P, U, 1 )
+            IF( SISNAN( DNORM ) .OR. SISNAN( XNORM ) ) THEN
+               RESULT = ONE / EPS
+            ELSE IF( XNORM.GT.ZERO ) THEN
+               RESULT = MAX( RESULT,
+     $                       ( ( DNORM / YNORM ) / XNORM ) / EPS )
+            END IF
+   30    CONTINUE
+      END IF
 *
       RETURN
 *
diff --git a/TESTING/EIG/slsets.f b/TESTING/EIG/slsets.f
index 53fc03af7..6dffa7504 100644
--- a/TESTING/EIG/slsets.f
+++ b/TESTING/EIG/slsets.f
@@ -27,7 +27,7 @@
 *> \verbatim
 *>
 *> SLSETS tests SGGLSE - a subroutine for solving linear equality
-*> constrained least square problem (LSE).
+*> constrained least square problem (LSE), including scaled inputs.
 *> \endverbatim
 *
 *  Arguments:
@@ -171,10 +171,25 @@ SUBROUTINE SLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF,
 *
 *     ..
 *     .. Local Scalars ..
-      INTEGER            INFO
+      INTEGER            INFO, ISCALE, J
+      REAL              ASCL, BSCL, CSCL, DSCL, RESID, SCL, TNRM
+*     ..
+*     .. Parameters ..
+      REAL              ZERO, ONE
+      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
+      INTEGER            MAXEXP
+      PARAMETER          ( MAXEXP = MAXEXPONENT( ZERO ) - 2 )
+*     ..
+*     .. External Functions ..
+      LOGICAL            SISNAN
+      REAL              SLANGE, SLAMCH
+      EXTERNAL           SISNAN, SLANGE, SLAMCH
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           SCOPY, SGGLSE, SLACPY, SGET02
+      EXTERNAL           SCOPY, SGGLSE, SLACPY, SGET02, SSCAL
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          EXPONENT, MAX, MAXEXPONENT, SCALE
 *     ..
 *     .. Executable Statements ..
 *
@@ -204,6 +219,76 @@ SUBROUTINE SLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF,
 *
       CALL SGET02( 'No transpose', P, N, 1, B, LDB, X, N, DF, P,
      $             RWORK, RESULT( 2 ) )
+*
+*     Scaling (A,c) and (B,d) independently leaves x unchanged.
+*     Scaling only (c,d) scales x by the same factor.  Exercise
+*     both ends of the range and check residuals at the input scale.
+*
+      TNRM = MAX( SLANGE( 'M', M, N, A, LDA, RWORK ),
+     $           SLANGE( 'M', P, N, B, LDB, RWORK ),
+     $           SLANGE( 'M', M, 1, C, M, RWORK ),
+     $           SLANGE( 'M', P, 1, D, P, RWORK ) )
+      IF( TNRM.GT.ZERO .AND. TNRM.LE.SLAMCH( 'Overflow' ) ) THEN
+*
+*        Cases: common large, common tiny, tiny (c,d), tiny (A,c),
+*        and tiny (B,d).
+*
+         DO 30 ISCALE = 1, 5
+            IF( ISCALE.EQ.1 ) THEN
+               SCL = SCALE( ONE, MAXEXP-EXPONENT( TNRM ) )
+            ELSE
+               SCL = SLAMCH( 'Safe minimum' ) /
+     $               SLAMCH( 'Precision' )
+               SCL = SCALE( ONE, EXPONENT( SCL )-4-
+     $                      EXPONENT( TNRM ) )
+            END IF
+            ASCL = SCL
+            BSCL = SCL
+            CSCL = SCL
+            DSCL = SCL
+            IF( ISCALE.EQ.3 .OR. ISCALE.EQ.5 ) ASCL = ONE
+            IF( ISCALE.EQ.3 .OR. ISCALE.EQ.4 ) BSCL = ONE
+            IF( ISCALE.EQ.4 ) DSCL = ONE
+            IF( ISCALE.EQ.5 ) CSCL = ONE
+            CALL SLACPY( 'Full', M, N, A, LDA, AF, LDA )
+            CALL SLACPY( 'Full', P, N, B, LDB, BF, LDB )
+            CALL SCOPY( M, C, 1, CF, 1 )
+            CALL SCOPY( P, D, 1, DF, 1 )
+            DO 10 J = 1, N
+               CALL SSCAL( M, ASCL, AF( 1, J ), 1 )
+               CALL SSCAL( P, BSCL, BF( 1, J ), 1 )
+   10       CONTINUE
+            CALL SSCAL( M, CSCL, CF, 1 )
+            CALL SSCAL( P, DSCL, DF, 1 )
+*
+            CALL SGGLSE( M, N, P, AF, LDA, BF, LDB, CF, DF, X, WORK,
+     $                   LWORK, INFO )
+            IF( INFO.NE.0 ) THEN
+               RESULT( 1 ) = ONE / SLAMCH( 'Epsilon' )
+               RESULT( 2 ) = RESULT( 1 )
+               GO TO 30
+            END IF
+            IF( ISCALE.EQ.3 )
+     $         CALL SSCAL( N, ONE / SCL, X, 1 )
+*
+            CALL SCOPY( M, C, 1, CF, 1 )
+            CALL SCOPY( P, D, 1, DF, 1 )
+            CALL SGET02( 'No transpose', M, N, 1, A, LDA, X, N, CF, M,
+     $                   RWORK, RESID )
+            IF( SISNAN( RESID ) ) THEN
+               RESULT( 1 ) = ONE / SLAMCH( 'Epsilon' )
+            ELSE
+               RESULT( 1 ) = MAX( RESULT( 1 ), RESID )
+            END IF
+            CALL SGET02( 'No transpose', P, N, 1, B, LDB, X, N, DF, P,
+     $                   RWORK, RESID )
+            IF( SISNAN( RESID ) ) THEN
+               RESULT( 2 ) = ONE / SLAMCH( 'Epsilon' )
+            ELSE
+               RESULT( 2 ) = MAX( RESULT( 2 ), RESID )
+            END IF
+   30    CONTINUE
+      END IF
 *
       RETURN
 *
diff --git a/TESTING/EIG/zglmts.f b/TESTING/EIG/zglmts.f
index 54ceca051..764980869 100644
--- a/TESTING/EIG/zglmts.f
+++ b/TESTING/EIG/zglmts.f
@@ -24,7 +24,7 @@
 *> \verbatim
 *>
 *> ZGLMTS tests ZGGGLM - a subroutine for solving the generalized
-*> linear model problem.
+*> linear model problem, including independent scaling of its inputs.
 *> \endverbatim
 *
 *  Arguments:
@@ -165,23 +165,29 @@ SUBROUTINE ZGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U,
 *     .. Parameters ..
       DOUBLE PRECISION   ZERO
       PARAMETER          ( ZERO = 0.0D+0 )
+      DOUBLE PRECISION   ONE
+      PARAMETER          ( ONE = 1.0D+0 )
+      INTEGER            MAXEXP
+      PARAMETER          ( MAXEXP = MAXEXPONENT( ZERO ) - 2 )
       COMPLEX*16         CONE
       PARAMETER          ( CONE = 1.0D+0 )
 *     ..
 *     .. Local Scalars ..
-      INTEGER            INFO
+      INTEGER            INFO, ISCALE, J
       DOUBLE PRECISION   ANORM, BNORM, DNORM, EPS, UNFL, XNORM, YNORM
+      DOUBLE PRECISION   ASCL, BSCL, DSCL, SCL, TNRM
 *     ..
 *     .. External Functions ..
       DOUBLE PRECISION   DLAMCH, DZASUM, ZLANGE
-      EXTERNAL           DLAMCH, DZASUM, ZLANGE
+      LOGICAL            DISNAN
+      EXTERNAL           DISNAN, DLAMCH, DZASUM, ZLANGE
 *     ..
 *     .. External Subroutines ..
 *
-      EXTERNAL           ZCOPY, ZGEMV, ZGGGLM, ZLACPY
+      EXTERNAL           ZCOPY, ZGEMV, ZGGGLM, ZLACPY, ZDSCAL
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          MAX
+      INTRINSIC          EXPONENT, MAX, MAXEXPONENT, SCALE
 *     ..
 *     .. Executable Statements ..
 *
@@ -224,6 +230,68 @@ SUBROUTINE ZGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U,
       ELSE
          RESULT = ( ( DNORM / YNORM ) / XNORM ) / EPS
       END IF
+*
+*     A, B and d may be scaled independently: for factors a, b, d,
+*     the solutions become (d/a)*x and (d/b)*u.  Test large and
+*     tiny inputs, then undo solution scaling before the residual.
+*
+      TNRM = MAX( ZLANGE( 'M', N, M, A, LDA, RWORK ),
+     $           ZLANGE( 'M', N, P, B, LDB, RWORK ),
+     $           ZLANGE( 'M', N, 1, D, N, RWORK ) )
+      IF( TNRM.GT.ZERO .AND. TNRM.LE.DLAMCH( 'Overflow' ) ) THEN
+*
+*        Cases: common large, common tiny, tiny d, tiny (A,d),
+*        and tiny (B,d).
+*
+         DO 30 ISCALE = 1, 5
+            IF( ISCALE.EQ.1 ) THEN
+               SCL = SCALE( ONE, MAXEXP-EXPONENT( TNRM ) )
+            ELSE
+               SCL = DLAMCH( 'Safe minimum' ) /
+     $               DLAMCH( 'Precision' )
+               SCL = SCALE( ONE, EXPONENT( SCL )-4-
+     $                      EXPONENT( TNRM ) )
+            END IF
+            ASCL = SCL
+            BSCL = SCL
+            DSCL = SCL
+            IF( ISCALE.EQ.3 .OR. ISCALE.EQ.5 ) ASCL = ONE
+            IF( ISCALE.EQ.3 .OR. ISCALE.EQ.4 ) BSCL = ONE
+            CALL ZLACPY( 'Full', N, M, A, LDA, AF, LDA )
+            CALL ZLACPY( 'Full', N, P, B, LDB, BF, LDB )
+            CALL ZCOPY( N, D, 1, DF, 1 )
+            DO 10 J = 1, M
+               CALL ZDSCAL( N, ASCL, AF( 1, J ), 1 )
+   10       CONTINUE
+            DO 20 J = 1, P
+               CALL ZDSCAL( N, BSCL, BF( 1, J ), 1 )
+   20       CONTINUE
+            CALL ZDSCAL( N, DSCL, DF, 1 )
+*
+            CALL ZGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK,
+     $                   LWORK, INFO )
+            IF( INFO.NE.0 ) THEN
+               RESULT = ONE / EPS
+               GO TO 30
+            END IF
+            CALL ZDSCAL( M, ASCL / DSCL, X, 1 )
+            CALL ZDSCAL( P, BSCL / DSCL, U, 1 )
+*
+            CALL ZCOPY( N, D, 1, DF, 1 )
+            CALL ZGEMV( 'No transpose', N, M, -CONE, A, LDA, X, 1, CONE,
+     $                  DF, 1 )
+            CALL ZGEMV( 'No transpose', N, P, -CONE, B, LDB, U, 1, CONE,
+     $                  DF, 1 )
+            DNORM = DZASUM( N, DF, 1 )
+            XNORM = DZASUM( M, X, 1 ) + DZASUM( P, U, 1 )
+            IF( DISNAN( DNORM ) .OR. DISNAN( XNORM ) ) THEN
+               RESULT = ONE / EPS
+            ELSE IF( XNORM.GT.ZERO ) THEN
+               RESULT = MAX( RESULT,
+     $                       ( ( DNORM / YNORM ) / XNORM ) / EPS )
+            END IF
+   30    CONTINUE
+      END IF
 *
       RETURN
 *
diff --git a/TESTING/EIG/zlsets.f b/TESTING/EIG/zlsets.f
index ede45d7af..c9d00ab4e 100644
--- a/TESTING/EIG/zlsets.f
+++ b/TESTING/EIG/zlsets.f
@@ -23,7 +23,7 @@
 *> \verbatim
 *>
 *> ZLSETS tests ZGGLSE - a subroutine for solving linear equality
-*> constrained least square problem (LSE).
+*> constrained least square problem (LSE), including scaled inputs.
 *> \endverbatim
 *
 *  Arguments:
@@ -167,10 +167,25 @@ SUBROUTINE ZLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF, D, DF,
      $                   WORK( LWORK ), X( * )
 *     ..
 *     .. Local Scalars ..
-      INTEGER            INFO
+      INTEGER            INFO, ISCALE, J
+      DOUBLE PRECISION   ASCL, BSCL, CSCL, DSCL, RESID, SCL, TNRM
+*     ..
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
+      INTEGER            MAXEXP
+      PARAMETER          ( MAXEXP = MAXEXPONENT( ZERO ) - 2 )
+*     ..
+*     .. External Functions ..
+      LOGICAL            DISNAN
+      DOUBLE PRECISION   ZLANGE, DLAMCH
+      EXTERNAL           DISNAN, ZLANGE, DLAMCH
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           ZCOPY, ZGET02, ZGGLSE, ZLACPY
+      EXTERNAL           ZCOPY, ZGET02, ZGGLSE, ZLACPY, ZDSCAL
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          EXPONENT, MAX, MAXEXPONENT, SCALE
 *     ..
 *     .. Executable Statements ..
 *
@@ -200,6 +215,76 @@ SUBROUTINE ZLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF, D, DF,
 *
       CALL ZGET02( 'No transpose', P, N, 1, B, LDB, X, N, DF, P, RWORK,
      $             RESULT( 2 ) )
+*
+*     Scaling (A,c) and (B,d) independently leaves x unchanged.
+*     Scaling only (c,d) scales x by the same factor.  Exercise
+*     both ends of the range and check residuals at the input scale.
+*
+      TNRM = MAX( ZLANGE( 'M', M, N, A, LDA, RWORK ),
+     $           ZLANGE( 'M', P, N, B, LDB, RWORK ),
+     $           ZLANGE( 'M', M, 1, C, M, RWORK ),
+     $           ZLANGE( 'M', P, 1, D, P, RWORK ) )
+      IF( TNRM.GT.ZERO .AND. TNRM.LE.DLAMCH( 'Overflow' ) ) THEN
+*
+*        Cases: common large, common tiny, tiny (c,d), tiny (A,c),
+*        and tiny (B,d).
+*
+         DO 30 ISCALE = 1, 5
+            IF( ISCALE.EQ.1 ) THEN
+               SCL = SCALE( ONE, MAXEXP-EXPONENT( TNRM ) )
+            ELSE
+               SCL = DLAMCH( 'Safe minimum' ) /
+     $               DLAMCH( 'Precision' )
+               SCL = SCALE( ONE, EXPONENT( SCL )-4-
+     $                      EXPONENT( TNRM ) )
+            END IF
+            ASCL = SCL
+            BSCL = SCL
+            CSCL = SCL
+            DSCL = SCL
+            IF( ISCALE.EQ.3 .OR. ISCALE.EQ.5 ) ASCL = ONE
+            IF( ISCALE.EQ.3 .OR. ISCALE.EQ.4 ) BSCL = ONE
+            IF( ISCALE.EQ.4 ) DSCL = ONE
+            IF( ISCALE.EQ.5 ) CSCL = ONE
+            CALL ZLACPY( 'Full', M, N, A, LDA, AF, LDA )
+            CALL ZLACPY( 'Full', P, N, B, LDB, BF, LDB )
+            CALL ZCOPY( M, C, 1, CF, 1 )
+            CALL ZCOPY( P, D, 1, DF, 1 )
+            DO 10 J = 1, N
+               CALL ZDSCAL( M, ASCL, AF( 1, J ), 1 )
+               CALL ZDSCAL( P, BSCL, BF( 1, J ), 1 )
+   10       CONTINUE
+            CALL ZDSCAL( M, CSCL, CF, 1 )
+            CALL ZDSCAL( P, DSCL, DF, 1 )
+*
+            CALL ZGGLSE( M, N, P, AF, LDA, BF, LDB, CF, DF, X, WORK,
+     $                   LWORK, INFO )
+            IF( INFO.NE.0 ) THEN
+               RESULT( 1 ) = ONE / DLAMCH( 'Epsilon' )
+               RESULT( 2 ) = RESULT( 1 )
+               GO TO 30
+            END IF
+            IF( ISCALE.EQ.3 )
+     $         CALL ZDSCAL( N, ONE / SCL, X, 1 )
+*
+            CALL ZCOPY( M, C, 1, CF, 1 )
+            CALL ZCOPY( P, D, 1, DF, 1 )
+            CALL ZGET02( 'No transpose', M, N, 1, A, LDA, X, N, CF, M,
+     $                   RWORK, RESID )
+            IF( DISNAN( RESID ) ) THEN
+               RESULT( 1 ) = ONE / DLAMCH( 'Epsilon' )
+            ELSE
+               RESULT( 1 ) = MAX( RESULT( 1 ), RESID )
+            END IF
+            CALL ZGET02( 'No transpose', P, N, 1, B, LDB, X, N, DF, P,
+     $                   RWORK, RESID )
+            IF( DISNAN( RESID ) ) THEN
+               RESULT( 2 ) = ONE / DLAMCH( 'Epsilon' )
+            ELSE
+               RESULT( 2 ) = MAX( RESULT( 2 ), RESID )
+            END IF
+   30    CONTINUE
+      END IF
 *
       RETURN
 *