From c7a7e38ca1bbf048085055144a8edc1335a39ccf Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Mon, 7 Sep 2026 21:38:00 -0700 Subject: [PATCH 1/2] Scale A, B and the right-hand sides of xGGLSE and xGGGLM independently, by exact powers of two Review of the previous revision found three regressions against the unscaled drivers. A and B shared one factor from their combined max-norm, so A = 2^1023, B = 2^-1023, d = 2^-1023 (x = 1) had B rounded to zero and returned INFO = 1. The two undo scalings of x cancelled through an intermediate that could flush, so x = (1, 2^-1050) came back as (1, 0). All of c was rescaled although only c(N-P+1:M) holds the residual, and the intermediate left in c(1:N-P) could overflow. Scale each operand by a power of two chosen from the EXPONENT of its max-norm instead. A and B get independent factors 2^KA and 2^KB. In xGGLSE, c carries 2^(KA+KT) and d carries 2^(KB+KT), which keeps the constraint and the residual consistent and scales x by 2^KT, with KT chosen to bring the larger right-hand side into range; on exit x is rescaled by 2^-KT and the residual in c(N-P+1:M) by 2^-(KA+KT). In xGGGLM, A, B and d are scaled independently and x and y are rescaled by 2^(KA-KD) and 2^(KB-KD). The factors are exact, so a problem the unscaled computation handles returns the same solution to the last bit, and each undo is one exact operation. xLASCL is called with both endpoints at or above one so that forming the factor raises no underflow; vectors are scaled with xSCAL. A zero, infinite or NaN norm takes no part. The three review cases and two neighbouring corners agree with the parent bit for bit and raise the same IEEE flags. Over the exponent sweep of 2860 cases the branch fails in none (the parent in 160) and, in the 1480 scaled cases with normal-range inputs, returns the solution of the unscaled twin bit for bit in 1388; the rest have entries beyond the thresholds at which xNRM2 switches accumulators and agree with the twin to rounding. xGLMTS and xLSETS solve their problem a second time with every operand scaled by one power of two, so that the largest entry sits just below the overflow threshold. The problem is exactly invariant under that scaling, so the residual of the second solution against the original data has to match the first; a NaN counts as a failure, which the comparison with the threshold would otherwise pass over. On the parent commit the GLM twin fails one of the 48 ratios in every precision. The LSE twin passes on both, because the right-hand side is the largest operand there and scaling it to the top of the range leaves the matrix below the exponent at which the reflector generator overflows; it covers the new scaling path rather than the defect. The full LAPACK test suite passes with the same totals as the parent. xGLMTS and xLSETS declare the new xSCAL calls EXTERNAL: the extended-API build renames only the routines a file declares, so without the declaration the xeigtst*_64 executables failed to link against the 64-bit BLAS. Co-Authored-By: Claude Fable 5.1 --- SRC/cggglm.f | 81 +++++++++++++++++++++++++++++++++++--- SRC/cgglse.f | 92 ++++++++++++++++++++++++++++++++++++++++--- SRC/dggglm.f | 79 ++++++++++++++++++++++++++++++++++--- SRC/dgglse.f | 94 ++++++++++++++++++++++++++++++++++++++++---- SRC/sggglm.f | 79 ++++++++++++++++++++++++++++++++++--- SRC/sgglse.f | 94 ++++++++++++++++++++++++++++++++++++++++---- SRC/zggglm.f | 81 +++++++++++++++++++++++++++++++++++--- SRC/zgglse.f | 92 ++++++++++++++++++++++++++++++++++++++++--- TESTING/EIG/cglmts.f | 51 ++++++++++++++++++++++-- TESTING/EIG/clsets.f | 61 +++++++++++++++++++++++++++- TESTING/EIG/dglmts.f | 49 +++++++++++++++++++++-- TESTING/EIG/dlsets.f | 61 +++++++++++++++++++++++++++- TESTING/EIG/sglmts.f | 49 +++++++++++++++++++++-- TESTING/EIG/slsets.f | 61 +++++++++++++++++++++++++++- TESTING/EIG/zglmts.f | 51 ++++++++++++++++++++++-- TESTING/EIG/zlsets.f | 61 +++++++++++++++++++++++++++- 16 files changed, 1068 insertions(+), 68 deletions(-) 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..08aee9eac 100644 --- a/TESTING/EIG/cglmts.f +++ b/TESTING/EIG/cglmts.f @@ -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, J REAL ANORM, BNORM, EPS, XNORM, YNORM, DNORM, UNFL + REAL SCL * .. * .. 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,43 @@ SUBROUTINE CGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, ELSE RESULT = ( ( DNORM / YNORM ) / XNORM ) /EPS END IF +* +* The problem is exactly invariant under scaling A, B and d by one +* power of two, so solving it again with the largest entry near the +* overflow threshold has to give the same residual. +* + SCL = MAX( CLANGE( 'M', N, M, A, LDA, RWORK ), + $ CLANGE( 'M', N, P, B, LDB, RWORK ), + $ CLANGE( 'M', N, 1, D, N, RWORK ) ) + IF( SCL.GT.ZERO .AND. SCL.LE.SLAMCH( 'Overflow' ) ) THEN + SCL = SCALE( ONE, MAXEXP-EXPONENT( SCL ) ) + 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, SCL, AF( 1, J ), 1 ) + 10 CONTINUE + DO 20 J = 1, P + CALL CSSCAL( N, SCL, BF( 1, J ), 1 ) + 20 CONTINUE + CALL CSSCAL( N, SCL, DF, 1 ) +* + CALL CGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, LWORK, + $ INFO ) +* + 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 + END IF * RETURN * diff --git a/TESTING/EIG/clsets.f b/TESTING/EIG/clsets.f index 7e940230e..8bc2c4cdb 100644 --- a/TESTING/EIG/clsets.f +++ b/TESTING/EIG/clsets.f @@ -171,10 +171,25 @@ SUBROUTINE CLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF, * * .. * .. Local Scalars .. - INTEGER INFO + INTEGER INFO, J + REAL RESID, SCL +* .. +* .. 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,48 @@ 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 ) ) +* +* The problem is exactly invariant under scaling A, B, c and d by +* one power of two, so solving it again with the largest entry near +* the overflow threshold has to give the same residuals. +* + SCL = 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( SCL.GT.ZERO .AND. SCL.LE.SLAMCH( 'Overflow' ) ) THEN + SCL = SCALE( ONE, MAXEXP-EXPONENT( SCL ) ) + 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, SCL, AF( 1, J ), 1 ) + CALL CSSCAL( P, SCL, BF( 1, J ), 1 ) + 10 CONTINUE + CALL CSSCAL( M, SCL, CF, 1 ) + CALL CSSCAL( P, SCL, DF, 1 ) +* + CALL CGGLSE( M, N, P, AF, LDA, BF, LDB, CF, DF, X, WORK, + $ LWORK, INFO ) +* + 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 + END IF * RETURN * diff --git a/TESTING/EIG/dglmts.f b/TESTING/EIG/dglmts.f index d2d24fea3..3adf38b51 100644 --- a/TESTING/EIG/dglmts.f +++ b/TESTING/EIG/dglmts.f @@ -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, J DOUBLE PRECISION ANORM, BNORM, DNORM, EPS, UNFL, XNORM, YNORM + DOUBLE PRECISION SCL * .. * .. 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,43 @@ SUBROUTINE DGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U, ELSE RESULT = ( ( DNORM / YNORM ) / XNORM ) / EPS END IF +* +* The problem is exactly invariant under scaling A, B and d by one +* power of two, so solving it again with the largest entry near the +* overflow threshold has to give the same residual. +* + SCL = MAX( DLANGE( 'M', N, M, A, LDA, RWORK ), + $ DLANGE( 'M', N, P, B, LDB, RWORK ), + $ DLANGE( 'M', N, 1, D, N, RWORK ) ) + IF( SCL.GT.ZERO .AND. SCL.LE.DLAMCH( 'Overflow' ) ) THEN + SCL = SCALE( ONE, MAXEXP-EXPONENT( SCL ) ) + 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, SCL, AF( 1, J ), 1 ) + 10 CONTINUE + DO 20 J = 1, P + CALL DSCAL( N, SCL, BF( 1, J ), 1 ) + 20 CONTINUE + CALL DSCAL( N, SCL, DF, 1 ) +* + CALL DGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, LWORK, + $ INFO ) +* + 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 + END IF * RETURN * diff --git a/TESTING/EIG/dlsets.f b/TESTING/EIG/dlsets.f index afb46a727..a42b60777 100644 --- a/TESTING/EIG/dlsets.f +++ b/TESTING/EIG/dlsets.f @@ -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, J + DOUBLE PRECISION RESID, SCL +* .. +* .. 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,48 @@ 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 ) ) +* +* The problem is exactly invariant under scaling A, B, c and d by +* one power of two, so solving it again with the largest entry near +* the overflow threshold has to give the same residuals. +* + SCL = 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( SCL.GT.ZERO .AND. SCL.LE.DLAMCH( 'Overflow' ) ) THEN + SCL = SCALE( ONE, MAXEXP-EXPONENT( SCL ) ) + 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, SCL, AF( 1, J ), 1 ) + CALL DSCAL( P, SCL, BF( 1, J ), 1 ) + 10 CONTINUE + CALL DSCAL( M, SCL, CF, 1 ) + CALL DSCAL( P, SCL, DF, 1 ) +* + CALL DGGLSE( M, N, P, AF, LDA, BF, LDB, CF, DF, X, WORK, + $ LWORK, INFO ) +* + 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 + END IF * RETURN * diff --git a/TESTING/EIG/sglmts.f b/TESTING/EIG/sglmts.f index b4b3e228d..5bf1255e0 100644 --- a/TESTING/EIG/sglmts.f +++ b/TESTING/EIG/sglmts.f @@ -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, J REAL ANORM, BNORM, EPS, XNORM, YNORM, DNORM, UNFL + REAL SCL * .. * .. 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,43 @@ SUBROUTINE SGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, ELSE RESULT = ( ( DNORM / YNORM ) / XNORM ) /EPS END IF +* +* The problem is exactly invariant under scaling A, B and d by one +* power of two, so solving it again with the largest entry near the +* overflow threshold has to give the same residual. +* + SCL = MAX( SLANGE( 'M', N, M, A, LDA, RWORK ), + $ SLANGE( 'M', N, P, B, LDB, RWORK ), + $ SLANGE( 'M', N, 1, D, N, RWORK ) ) + IF( SCL.GT.ZERO .AND. SCL.LE.SLAMCH( 'Overflow' ) ) THEN + SCL = SCALE( ONE, MAXEXP-EXPONENT( SCL ) ) + 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, SCL, AF( 1, J ), 1 ) + 10 CONTINUE + DO 20 J = 1, P + CALL SSCAL( N, SCL, BF( 1, J ), 1 ) + 20 CONTINUE + CALL SSCAL( N, SCL, DF, 1 ) +* + CALL SGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, LWORK, + $ INFO ) +* + 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 + END IF * RETURN * diff --git a/TESTING/EIG/slsets.f b/TESTING/EIG/slsets.f index 53fc03af7..0c0c9a997 100644 --- a/TESTING/EIG/slsets.f +++ b/TESTING/EIG/slsets.f @@ -171,10 +171,25 @@ SUBROUTINE SLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF, * * .. * .. Local Scalars .. - INTEGER INFO + INTEGER INFO, J + REAL RESID, SCL +* .. +* .. 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,48 @@ 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 ) ) +* +* The problem is exactly invariant under scaling A, B, c and d by +* one power of two, so solving it again with the largest entry near +* the overflow threshold has to give the same residuals. +* + SCL = 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( SCL.GT.ZERO .AND. SCL.LE.SLAMCH( 'Overflow' ) ) THEN + SCL = SCALE( ONE, MAXEXP-EXPONENT( SCL ) ) + 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, SCL, AF( 1, J ), 1 ) + CALL SSCAL( P, SCL, BF( 1, J ), 1 ) + 10 CONTINUE + CALL SSCAL( M, SCL, CF, 1 ) + CALL SSCAL( P, SCL, DF, 1 ) +* + CALL SGGLSE( M, N, P, AF, LDA, BF, LDB, CF, DF, X, WORK, + $ LWORK, INFO ) +* + 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 + END IF * RETURN * diff --git a/TESTING/EIG/zglmts.f b/TESTING/EIG/zglmts.f index 54ceca051..463cdf4bc 100644 --- a/TESTING/EIG/zglmts.f +++ b/TESTING/EIG/zglmts.f @@ -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, J DOUBLE PRECISION ANORM, BNORM, DNORM, EPS, UNFL, XNORM, YNORM + DOUBLE PRECISION SCL * .. * .. 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,43 @@ SUBROUTINE ZGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U, ELSE RESULT = ( ( DNORM / YNORM ) / XNORM ) / EPS END IF +* +* The problem is exactly invariant under scaling A, B and d by one +* power of two, so solving it again with the largest entry near the +* overflow threshold has to give the same residual. +* + SCL = MAX( ZLANGE( 'M', N, M, A, LDA, RWORK ), + $ ZLANGE( 'M', N, P, B, LDB, RWORK ), + $ ZLANGE( 'M', N, 1, D, N, RWORK ) ) + IF( SCL.GT.ZERO .AND. SCL.LE.DLAMCH( 'Overflow' ) ) THEN + SCL = SCALE( ONE, MAXEXP-EXPONENT( SCL ) ) + 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, SCL, AF( 1, J ), 1 ) + 10 CONTINUE + DO 20 J = 1, P + CALL ZDSCAL( N, SCL, BF( 1, J ), 1 ) + 20 CONTINUE + CALL ZDSCAL( N, SCL, DF, 1 ) +* + CALL ZGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, LWORK, + $ INFO ) +* + 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 + END IF * RETURN * diff --git a/TESTING/EIG/zlsets.f b/TESTING/EIG/zlsets.f index ede45d7af..5fd3ddaab 100644 --- a/TESTING/EIG/zlsets.f +++ b/TESTING/EIG/zlsets.f @@ -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, J + DOUBLE PRECISION RESID, SCL +* .. +* .. 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,48 @@ 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 ) ) +* +* The problem is exactly invariant under scaling A, B, c and d by +* one power of two, so solving it again with the largest entry near +* the overflow threshold has to give the same residuals. +* + SCL = 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( SCL.GT.ZERO .AND. SCL.LE.DLAMCH( 'Overflow' ) ) THEN + SCL = SCALE( ONE, MAXEXP-EXPONENT( SCL ) ) + 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, SCL, AF( 1, J ), 1 ) + CALL ZDSCAL( P, SCL, BF( 1, J ), 1 ) + 10 CONTINUE + CALL ZDSCAL( M, SCL, CF, 1 ) + CALL ZDSCAL( P, SCL, DF, 1 ) +* + CALL ZGGLSE( M, N, P, AF, LDA, BF, LDB, CF, DF, X, WORK, + $ LWORK, INFO ) +* + 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 + END IF * RETURN * From 6b929b3960068eda2448180c58025868bfd3954a Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Sat, 12 Sep 2026 12:20:52 -0700 Subject: [PATCH 2/2] TESTING: Cover tiny inputs in xGGLSE and xGGGLM Extend the existing large-input checks with common and independent small-input scaling cases. Undo the corresponding solution scaling and check residuals against the original problem in all four precisions. Check INFO before using each scaled solve's result. All eight GLM/LSE suites pass with GNU 13. Local gcov confirms execution of all 24 previously uncovered upward-scaling assignments in the eight drivers. --- TESTING/EIG/cglmts.f | 95 ++++++++++++++++++++++++-------------- TESTING/EIG/clsets.f | 106 +++++++++++++++++++++++++++---------------- TESTING/EIG/dglmts.f | 95 ++++++++++++++++++++++++-------------- TESTING/EIG/dlsets.f | 106 +++++++++++++++++++++++++++---------------- TESTING/EIG/sglmts.f | 95 ++++++++++++++++++++++++-------------- TESTING/EIG/slsets.f | 106 +++++++++++++++++++++++++++---------------- TESTING/EIG/zglmts.f | 95 ++++++++++++++++++++++++-------------- TESTING/EIG/zlsets.f | 106 +++++++++++++++++++++++++++---------------- 8 files changed, 508 insertions(+), 296 deletions(-) diff --git a/TESTING/EIG/cglmts.f b/TESTING/EIG/cglmts.f index 08aee9eac..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: @@ -176,9 +176,9 @@ SUBROUTINE CGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, PARAMETER ( CONE = 1.0E+0 ) * .. * .. Local Scalars .. - INTEGER INFO, J + INTEGER INFO, ISCALE, J REAL ANORM, BNORM, EPS, XNORM, YNORM, DNORM, UNFL - REAL SCL + REAL ASCL, BSCL, DSCL, SCL, TNRM * .. * .. External Functions .. REAL SCASUM, SLAMCH, CLANGE @@ -233,41 +233,66 @@ SUBROUTINE CGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, RESULT = ( ( DNORM / YNORM ) / XNORM ) /EPS END IF * -* The problem is exactly invariant under scaling A, B and d by one -* power of two, so solving it again with the largest entry near the -* overflow threshold has to give the same residual. +* 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. * - SCL = MAX( CLANGE( 'M', N, M, A, LDA, RWORK ), + TNRM = MAX( CLANGE( 'M', N, M, A, LDA, RWORK ), $ CLANGE( 'M', N, P, B, LDB, RWORK ), $ CLANGE( 'M', N, 1, D, N, RWORK ) ) - IF( SCL.GT.ZERO .AND. SCL.LE.SLAMCH( 'Overflow' ) ) THEN - SCL = SCALE( ONE, MAXEXP-EXPONENT( SCL ) ) - 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, SCL, AF( 1, J ), 1 ) - 10 CONTINUE - DO 20 J = 1, P - CALL CSSCAL( N, SCL, BF( 1, J ), 1 ) - 20 CONTINUE - CALL CSSCAL( N, SCL, DF, 1 ) -* - CALL CGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, LWORK, - $ INFO ) -* - 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 + 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 8bc2c4cdb..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,8 +171,8 @@ SUBROUTINE CLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF, * * .. * .. Local Scalars .. - INTEGER INFO, J - REAL RESID, SCL + INTEGER INFO, ISCALE, J + REAL ASCL, BSCL, CSCL, DSCL, RESID, SCL, TNRM * .. * .. Parameters .. REAL ZERO, ONE @@ -220,46 +220,74 @@ 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 ) ) * -* The problem is exactly invariant under scaling A, B, c and d by -* one power of two, so solving it again with the largest entry near -* the overflow threshold has to give the same residuals. +* 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. * - SCL = MAX( CLANGE( 'M', M, N, A, LDA, RWORK ), + 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( SCL.GT.ZERO .AND. SCL.LE.SLAMCH( 'Overflow' ) ) THEN - SCL = SCALE( ONE, MAXEXP-EXPONENT( SCL ) ) - 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, SCL, AF( 1, J ), 1 ) - CALL CSSCAL( P, SCL, BF( 1, J ), 1 ) - 10 CONTINUE - CALL CSSCAL( M, SCL, CF, 1 ) - CALL CSSCAL( P, SCL, DF, 1 ) -* - CALL CGGLSE( M, N, P, AF, LDA, BF, LDB, CF, DF, X, WORK, - $ LWORK, INFO ) -* - 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 + 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 3adf38b51..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: @@ -168,9 +168,9 @@ SUBROUTINE DGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U, PARAMETER ( MAXEXP = MAXEXPONENT( ZERO ) - 2 ) * .. * .. Local Scalars .. - INTEGER INFO, J + INTEGER INFO, ISCALE, J DOUBLE PRECISION ANORM, BNORM, DNORM, EPS, UNFL, XNORM, YNORM - DOUBLE PRECISION SCL + DOUBLE PRECISION ASCL, BSCL, DSCL, SCL, TNRM * .. * .. External Functions .. DOUBLE PRECISION DASUM, DLAMCH, DLANGE @@ -224,41 +224,66 @@ SUBROUTINE DGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U, RESULT = ( ( DNORM / YNORM ) / XNORM ) / EPS END IF * -* The problem is exactly invariant under scaling A, B and d by one -* power of two, so solving it again with the largest entry near the -* overflow threshold has to give the same residual. +* 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. * - SCL = MAX( DLANGE( 'M', N, M, A, LDA, RWORK ), + TNRM = MAX( DLANGE( 'M', N, M, A, LDA, RWORK ), $ DLANGE( 'M', N, P, B, LDB, RWORK ), $ DLANGE( 'M', N, 1, D, N, RWORK ) ) - IF( SCL.GT.ZERO .AND. SCL.LE.DLAMCH( 'Overflow' ) ) THEN - SCL = SCALE( ONE, MAXEXP-EXPONENT( SCL ) ) - 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, SCL, AF( 1, J ), 1 ) - 10 CONTINUE - DO 20 J = 1, P - CALL DSCAL( N, SCL, BF( 1, J ), 1 ) - 20 CONTINUE - CALL DSCAL( N, SCL, DF, 1 ) -* - CALL DGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, LWORK, - $ INFO ) -* - 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 + 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 a42b60777..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,8 +166,8 @@ 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, J - DOUBLE PRECISION RESID, SCL + INTEGER INFO, ISCALE, J + DOUBLE PRECISION ASCL, BSCL, CSCL, DSCL, RESID, SCL, TNRM * .. * .. Parameters .. DOUBLE PRECISION ZERO, ONE @@ -215,46 +215,74 @@ 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 ) ) * -* The problem is exactly invariant under scaling A, B, c and d by -* one power of two, so solving it again with the largest entry near -* the overflow threshold has to give the same residuals. +* 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. * - SCL = MAX( DLANGE( 'M', M, N, A, LDA, RWORK ), + 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( SCL.GT.ZERO .AND. SCL.LE.DLAMCH( 'Overflow' ) ) THEN - SCL = SCALE( ONE, MAXEXP-EXPONENT( SCL ) ) - 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, SCL, AF( 1, J ), 1 ) - CALL DSCAL( P, SCL, BF( 1, J ), 1 ) - 10 CONTINUE - CALL DSCAL( M, SCL, CF, 1 ) - CALL DSCAL( P, SCL, DF, 1 ) -* - CALL DGGLSE( M, N, P, AF, LDA, BF, LDB, CF, DF, X, WORK, - $ LWORK, INFO ) -* - 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 + 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 5bf1255e0..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: @@ -170,9 +170,9 @@ SUBROUTINE SGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, PARAMETER ( MAXEXP = MAXEXPONENT( ZERO ) - 2 ) * .. * .. Local Scalars .. - INTEGER INFO, J + INTEGER INFO, ISCALE, J REAL ANORM, BNORM, EPS, XNORM, YNORM, DNORM, UNFL - REAL SCL + REAL ASCL, BSCL, DSCL, SCL, TNRM * .. * .. External Functions .. REAL SASUM, SLAMCH, SLANGE @@ -227,41 +227,66 @@ SUBROUTINE SGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, RESULT = ( ( DNORM / YNORM ) / XNORM ) /EPS END IF * -* The problem is exactly invariant under scaling A, B and d by one -* power of two, so solving it again with the largest entry near the -* overflow threshold has to give the same residual. +* 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. * - SCL = MAX( SLANGE( 'M', N, M, A, LDA, RWORK ), + TNRM = MAX( SLANGE( 'M', N, M, A, LDA, RWORK ), $ SLANGE( 'M', N, P, B, LDB, RWORK ), $ SLANGE( 'M', N, 1, D, N, RWORK ) ) - IF( SCL.GT.ZERO .AND. SCL.LE.SLAMCH( 'Overflow' ) ) THEN - SCL = SCALE( ONE, MAXEXP-EXPONENT( SCL ) ) - 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, SCL, AF( 1, J ), 1 ) - 10 CONTINUE - DO 20 J = 1, P - CALL SSCAL( N, SCL, BF( 1, J ), 1 ) - 20 CONTINUE - CALL SSCAL( N, SCL, DF, 1 ) -* - CALL SGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, LWORK, - $ INFO ) -* - 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 + 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 0c0c9a997..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,8 +171,8 @@ SUBROUTINE SLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF, * * .. * .. Local Scalars .. - INTEGER INFO, J - REAL RESID, SCL + INTEGER INFO, ISCALE, J + REAL ASCL, BSCL, CSCL, DSCL, RESID, SCL, TNRM * .. * .. Parameters .. REAL ZERO, ONE @@ -220,46 +220,74 @@ 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 ) ) * -* The problem is exactly invariant under scaling A, B, c and d by -* one power of two, so solving it again with the largest entry near -* the overflow threshold has to give the same residuals. +* 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. * - SCL = MAX( SLANGE( 'M', M, N, A, LDA, RWORK ), + 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( SCL.GT.ZERO .AND. SCL.LE.SLAMCH( 'Overflow' ) ) THEN - SCL = SCALE( ONE, MAXEXP-EXPONENT( SCL ) ) - 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, SCL, AF( 1, J ), 1 ) - CALL SSCAL( P, SCL, BF( 1, J ), 1 ) - 10 CONTINUE - CALL SSCAL( M, SCL, CF, 1 ) - CALL SSCAL( P, SCL, DF, 1 ) -* - CALL SGGLSE( M, N, P, AF, LDA, BF, LDB, CF, DF, X, WORK, - $ LWORK, INFO ) -* - 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 + 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 463cdf4bc..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: @@ -173,9 +173,9 @@ SUBROUTINE ZGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U, PARAMETER ( CONE = 1.0D+0 ) * .. * .. Local Scalars .. - INTEGER INFO, J + INTEGER INFO, ISCALE, J DOUBLE PRECISION ANORM, BNORM, DNORM, EPS, UNFL, XNORM, YNORM - DOUBLE PRECISION SCL + DOUBLE PRECISION ASCL, BSCL, DSCL, SCL, TNRM * .. * .. External Functions .. DOUBLE PRECISION DLAMCH, DZASUM, ZLANGE @@ -231,41 +231,66 @@ SUBROUTINE ZGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U, RESULT = ( ( DNORM / YNORM ) / XNORM ) / EPS END IF * -* The problem is exactly invariant under scaling A, B and d by one -* power of two, so solving it again with the largest entry near the -* overflow threshold has to give the same residual. +* 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. * - SCL = MAX( ZLANGE( 'M', N, M, A, LDA, RWORK ), + TNRM = MAX( ZLANGE( 'M', N, M, A, LDA, RWORK ), $ ZLANGE( 'M', N, P, B, LDB, RWORK ), $ ZLANGE( 'M', N, 1, D, N, RWORK ) ) - IF( SCL.GT.ZERO .AND. SCL.LE.DLAMCH( 'Overflow' ) ) THEN - SCL = SCALE( ONE, MAXEXP-EXPONENT( SCL ) ) - 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, SCL, AF( 1, J ), 1 ) - 10 CONTINUE - DO 20 J = 1, P - CALL ZDSCAL( N, SCL, BF( 1, J ), 1 ) - 20 CONTINUE - CALL ZDSCAL( N, SCL, DF, 1 ) -* - CALL ZGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, LWORK, - $ INFO ) -* - 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 + 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 5fd3ddaab..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,8 +167,8 @@ SUBROUTINE ZLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF, D, DF, $ WORK( LWORK ), X( * ) * .. * .. Local Scalars .. - INTEGER INFO, J - DOUBLE PRECISION RESID, SCL + INTEGER INFO, ISCALE, J + DOUBLE PRECISION ASCL, BSCL, CSCL, DSCL, RESID, SCL, TNRM * .. * .. Parameters .. DOUBLE PRECISION ZERO, ONE @@ -216,46 +216,74 @@ 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 ) ) * -* The problem is exactly invariant under scaling A, B, c and d by -* one power of two, so solving it again with the largest entry near -* the overflow threshold has to give the same residuals. +* 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. * - SCL = MAX( ZLANGE( 'M', M, N, A, LDA, RWORK ), + 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( SCL.GT.ZERO .AND. SCL.LE.DLAMCH( 'Overflow' ) ) THEN - SCL = SCALE( ONE, MAXEXP-EXPONENT( SCL ) ) - 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, SCL, AF( 1, J ), 1 ) - CALL ZDSCAL( P, SCL, BF( 1, J ), 1 ) - 10 CONTINUE - CALL ZDSCAL( M, SCL, CF, 1 ) - CALL ZDSCAL( P, SCL, DF, 1 ) -* - CALL ZGGLSE( M, N, P, AF, LDA, BF, LDB, CF, DF, X, WORK, - $ LWORK, INFO ) -* - 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 + 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