diff --git a/SRC/clarfg.f b/SRC/clarfg.f index d9e6cf4bc..b77be5446 100644 --- a/SRC/clarfg.f +++ b/SRC/clarfg.f @@ -118,8 +118,9 @@ SUBROUTINE CLARFG( N, ALPHA, X, INCX, TAU ) * ===================================================================== * * .. Parameters .. - REAL ONE, ZERO - PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) + REAL ONE, ZERO, HALF + PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0, + $ HALF = 0.5E+0 ) * .. * .. Local Scalars .. INTEGER J, KNT @@ -131,7 +132,7 @@ SUBROUTINE CLARFG( N, ALPHA, X, INCX, TAU ) EXTERNAL SCNRM2, SLAMCH, SLAPY3, CLADIV * .. * .. Intrinsic Functions .. - INTRINSIC ABS, AIMAG, CMPLX, REAL, SIGN + INTRINSIC ABS, AIMAG, CMPLX, HUGE, REAL, SIGN * .. * .. External Subroutines .. EXTERNAL CSCAL, CSSCAL @@ -179,16 +180,31 @@ SUBROUTINE CLARFG( N, ALPHA, X, INCX, TAU ) XNORM = SCNRM2( N-1, X, INCX ) ALPHA = CMPLX( ALPHR, ALPHI ) BETA = -SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) + ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN +* +* |ALPHA| <= |BETA|, so ALPHA-BETA can overflow only when +* |BETA| > HUGE/2; scale X down and recompute them. +* + KNT = -1 + CALL CSSCAL( N-1, SAFMIN, X, INCX ) + ALPHI = ALPHI*SAFMIN + ALPHR = ALPHR*SAFMIN + XNORM = SCNRM2( N-1, X, INCX ) + ALPHA = CMPLX( ALPHR, ALPHI ) + BETA = -SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) END IF TAU = CMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA ) ALPHA = CLADIV( CMPLX( ONE ), ALPHA-BETA ) CALL CSCAL( N-1, ALPHA, X, INCX ) * -* If ALPHA is subnormal, it may lose relative accuracy +* Undo the scaling. If ALPHA is subnormal, it may lose relative +* accuracy * DO 20 J = 1, KNT BETA = BETA*SAFMIN 20 CONTINUE + IF( KNT.LT.0 ) + $ BETA = BETA*RSAFMN ALPHA = BETA END IF * diff --git a/SRC/clarfgp.f b/SRC/clarfgp.f index ca24dc8ec..e0b4d6d68 100644 --- a/SRC/clarfgp.f +++ b/SRC/clarfgp.f @@ -116,8 +116,9 @@ SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU ) * ===================================================================== * * .. Parameters .. - REAL TWO, ONE, ZERO - PARAMETER ( TWO = 2.0E+0, ONE = 1.0E+0, ZERO = 0.0E+0 ) + REAL TWO, ONE, ZERO, HALF + PARAMETER ( TWO = 2.0E+0, ONE = 1.0E+0, ZERO = 0.0E+0, + $ HALF = 0.5E+0 ) * .. * .. Local Scalars .. INTEGER J, KNT @@ -131,7 +132,7 @@ SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU ) $ CLADIV * .. * .. Intrinsic Functions .. - INTRINSIC ABS, AIMAG, CMPLX, REAL, SIGN + INTRINSIC ABS, AIMAG, CMPLX, HUGE, REAL, SIGN * .. * .. External Subroutines .. EXTERNAL CSCAL, CSSCAL @@ -193,6 +194,18 @@ SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU ) XNORM = SCNRM2( N-1, X, INCX ) ALPHA = CMPLX( ALPHR, ALPHI ) BETA = SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) + ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN +* +* |ALPHA| <= |BETA|, so ALPHA+BETA can overflow only when +* |BETA| > HUGE/2; scale X down and recompute them. +* + KNT = -1 + CALL CSSCAL( N-1, SMLNUM, X, INCX ) + ALPHI = ALPHI*SMLNUM + ALPHR = ALPHR*SMLNUM + XNORM = SCNRM2( N-1, X, INCX ) + ALPHA = CMPLX( ALPHR, ALPHI ) + BETA = SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) END IF SAVEALPHA = ALPHA ALPHA = ALPHA + BETA @@ -245,11 +258,14 @@ SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU ) * END IF * -* If BETA is subnormal, it may lose relative accuracy +* Undo the scaling. If BETA is subnormal, it may lose relative +* accuracy * DO 20 J = 1, KNT BETA = BETA*SMLNUM 20 CONTINUE + IF( KNT.LT.0 ) + $ BETA = BETA*BIGNUM ALPHA = BETA END IF * diff --git a/SRC/dlarfg.f b/SRC/dlarfg.f index 38f65af0f..03c489b89 100644 --- a/SRC/dlarfg.f +++ b/SRC/dlarfg.f @@ -118,8 +118,9 @@ SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU ) * ===================================================================== * * .. Parameters .. - DOUBLE PRECISION ONE, ZERO - PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) + DOUBLE PRECISION ONE, ZERO, HALF + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0, + $ HALF = 0.5D+0 ) * .. * .. Local Scalars .. INTEGER J, KNT @@ -130,7 +131,7 @@ SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU ) EXTERNAL DLAMCH, DLAPY2, DNRM2 * .. * .. Intrinsic Functions .. - INTRINSIC ABS, SIGN + INTRINSIC ABS, HUGE, SIGN * .. * .. External Subroutines .. EXTERNAL DSCAL @@ -173,15 +174,29 @@ SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU ) * XNORM = DNRM2( N-1, X, INCX ) BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) + ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN +* +* |ALPHA| <= |BETA|, so ALPHA-BETA can overflow only when +* |BETA| > HUGE/2; scale X down and recompute them. +* + RSAFMN = ONE / SAFMIN + KNT = -1 + CALL DSCAL( N-1, SAFMIN, X, INCX ) + ALPHA = ALPHA*SAFMIN + XNORM = DNRM2( N-1, X, INCX ) + BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) END IF TAU = ( BETA-ALPHA ) / BETA CALL DSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX ) * -* If ALPHA is subnormal, it may lose relative accuracy +* Undo the scaling. If ALPHA is subnormal, it may lose relative +* accuracy * DO 20 J = 1, KNT BETA = BETA*SAFMIN 20 CONTINUE + IF( KNT.LT.0 ) + $ BETA = BETA*RSAFMN ALPHA = BETA END IF * diff --git a/SRC/dlarfgp.f b/SRC/dlarfgp.f index fe6613997..c74daa312 100644 --- a/SRC/dlarfgp.f +++ b/SRC/dlarfgp.f @@ -116,8 +116,9 @@ SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU ) * ===================================================================== * * .. Parameters .. - DOUBLE PRECISION TWO, ONE, ZERO - PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 ) + DOUBLE PRECISION TWO, ONE, ZERO, HALF + PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0, + $ HALF = 0.5D+0 ) * .. * .. Local Scalars .. INTEGER J, KNT @@ -128,7 +129,7 @@ SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU ) EXTERNAL DLAMCH, DLAPY2, DNRM2 * .. * .. Intrinsic Functions .. - INTRINSIC ABS, SIGN + INTRINSIC ABS, HUGE, SIGN * .. * .. External Subroutines .. EXTERNAL DSCAL @@ -185,6 +186,17 @@ SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU ) * XNORM = DNRM2( N-1, X, INCX ) BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) + ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN +* +* |ALPHA| <= |BETA|, so ALPHA+BETA can overflow only when +* |BETA| > HUGE/2; scale X down and recompute them. +* + BIGNUM = ONE / SMLNUM + KNT = -1 + CALL DSCAL( N-1, SMLNUM, X, INCX ) + ALPHA = ALPHA*SMLNUM + XNORM = DNRM2( N-1, X, INCX ) + BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) END IF SAVEALPHA = ALPHA ALPHA = ALPHA + BETA @@ -229,11 +241,14 @@ SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU ) * END IF * -* If BETA is subnormal, it may lose relative accuracy +* Undo the scaling. If BETA is subnormal, it may lose relative +* accuracy * DO 20 J = 1, KNT BETA = BETA*SMLNUM 20 CONTINUE + IF( KNT.LT.0 ) + $ BETA = BETA*BIGNUM ALPHA = BETA END IF * diff --git a/SRC/slarfg.f b/SRC/slarfg.f index a0893439d..2e7163c81 100644 --- a/SRC/slarfg.f +++ b/SRC/slarfg.f @@ -118,8 +118,9 @@ SUBROUTINE SLARFG( N, ALPHA, X, INCX, TAU ) * ===================================================================== * * .. Parameters .. - REAL ONE, ZERO - PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) + REAL ONE, ZERO, HALF + PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0, + $ HALF = 0.5E+0 ) * .. * .. Local Scalars .. INTEGER J, KNT @@ -130,7 +131,7 @@ SUBROUTINE SLARFG( N, ALPHA, X, INCX, TAU ) EXTERNAL SLAMCH, SLAPY2, SNRM2 * .. * .. Intrinsic Functions .. - INTRINSIC ABS, SIGN + INTRINSIC ABS, HUGE, SIGN * .. * .. External Subroutines .. EXTERNAL SSCAL @@ -173,15 +174,29 @@ SUBROUTINE SLARFG( N, ALPHA, X, INCX, TAU ) * XNORM = SNRM2( N-1, X, INCX ) BETA = -SIGN( SLAPY2( ALPHA, XNORM ), ALPHA ) + ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN +* +* |ALPHA| <= |BETA|, so ALPHA-BETA can overflow only when +* |BETA| > HUGE/2; scale X down and recompute them. +* + RSAFMN = ONE / SAFMIN + KNT = -1 + CALL SSCAL( N-1, SAFMIN, X, INCX ) + ALPHA = ALPHA*SAFMIN + XNORM = SNRM2( N-1, X, INCX ) + BETA = -SIGN( SLAPY2( ALPHA, XNORM ), ALPHA ) END IF TAU = ( BETA-ALPHA ) / BETA CALL SSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX ) * -* If ALPHA is subnormal, it may lose relative accuracy +* Undo the scaling. If ALPHA is subnormal, it may lose relative +* accuracy * DO 20 J = 1, KNT BETA = BETA*SAFMIN 20 CONTINUE + IF( KNT.LT.0 ) + $ BETA = BETA*RSAFMN ALPHA = BETA END IF * diff --git a/SRC/slarfgp.f b/SRC/slarfgp.f index 0a3d72313..d23a74b95 100644 --- a/SRC/slarfgp.f +++ b/SRC/slarfgp.f @@ -116,8 +116,9 @@ SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU ) * ===================================================================== * * .. Parameters .. - REAL TWO, ONE, ZERO - PARAMETER ( TWO = 2.0E+0, ONE = 1.0E+0, ZERO = 0.0E+0 ) + REAL TWO, ONE, ZERO, HALF + PARAMETER ( TWO = 2.0E+0, ONE = 1.0E+0, ZERO = 0.0E+0, + $ HALF = 0.5E+0 ) * .. * .. Local Scalars .. INTEGER J, KNT @@ -128,7 +129,7 @@ SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU ) EXTERNAL SLAMCH, SLAPY2, SNRM2 * .. * .. Intrinsic Functions .. - INTRINSIC ABS, SIGN + INTRINSIC ABS, HUGE, SIGN * .. * .. External Subroutines .. EXTERNAL SSCAL @@ -185,6 +186,17 @@ SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU ) * XNORM = SNRM2( N-1, X, INCX ) BETA = SIGN( SLAPY2( ALPHA, XNORM ), ALPHA ) + ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN +* +* |ALPHA| <= |BETA|, so ALPHA+BETA can overflow only when +* |BETA| > HUGE/2; scale X down and recompute them. +* + BIGNUM = ONE / SMLNUM + KNT = -1 + CALL SSCAL( N-1, SMLNUM, X, INCX ) + ALPHA = ALPHA*SMLNUM + XNORM = SNRM2( N-1, X, INCX ) + BETA = SIGN( SLAPY2( ALPHA, XNORM ), ALPHA ) END IF SAVEALPHA = ALPHA ALPHA = ALPHA + BETA @@ -229,11 +241,14 @@ SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU ) * END IF * -* If BETA is subnormal, it may lose relative accuracy +* Undo the scaling. If BETA is subnormal, it may lose relative +* accuracy * DO 20 J = 1, KNT BETA = BETA*SMLNUM 20 CONTINUE + IF( KNT.LT.0 ) + $ BETA = BETA*BIGNUM ALPHA = BETA END IF * diff --git a/SRC/zlarfg.f b/SRC/zlarfg.f index 576d9c64d..ece8ed5c9 100644 --- a/SRC/zlarfg.f +++ b/SRC/zlarfg.f @@ -118,8 +118,9 @@ SUBROUTINE ZLARFG( N, ALPHA, X, INCX, TAU ) * ===================================================================== * * .. Parameters .. - DOUBLE PRECISION ONE, ZERO - PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) + DOUBLE PRECISION ONE, ZERO, HALF + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0, + $ HALF = 0.5D+0 ) * .. * .. Local Scalars .. INTEGER J, KNT @@ -131,7 +132,7 @@ SUBROUTINE ZLARFG( N, ALPHA, X, INCX, TAU ) EXTERNAL DLAMCH, DLAPY3, DZNRM2, ZLADIV * .. * .. Intrinsic Functions .. - INTRINSIC ABS, DBLE, DCMPLX, DIMAG, SIGN + INTRINSIC ABS, DBLE, DCMPLX, DIMAG, HUGE, SIGN * .. * .. External Subroutines .. EXTERNAL ZDSCAL, ZSCAL @@ -179,16 +180,31 @@ SUBROUTINE ZLARFG( N, ALPHA, X, INCX, TAU ) XNORM = DZNRM2( N-1, X, INCX ) ALPHA = DCMPLX( ALPHR, ALPHI ) BETA = -SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) + ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN +* +* |ALPHA| <= |BETA|, so ALPHA-BETA can overflow only when +* |BETA| > HUGE/2; scale X down and recompute them. +* + KNT = -1 + CALL ZDSCAL( N-1, SAFMIN, X, INCX ) + ALPHI = ALPHI*SAFMIN + ALPHR = ALPHR*SAFMIN + XNORM = DZNRM2( N-1, X, INCX ) + ALPHA = DCMPLX( ALPHR, ALPHI ) + BETA = -SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) END IF TAU = DCMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA ) ALPHA = ZLADIV( DCMPLX( ONE ), ALPHA-BETA ) CALL ZSCAL( N-1, ALPHA, X, INCX ) * -* If ALPHA is subnormal, it may lose relative accuracy +* Undo the scaling. If ALPHA is subnormal, it may lose relative +* accuracy * DO 20 J = 1, KNT BETA = BETA*SAFMIN 20 CONTINUE + IF( KNT.LT.0 ) + $ BETA = BETA*RSAFMN ALPHA = BETA END IF * diff --git a/SRC/zlarfgp.f b/SRC/zlarfgp.f index 29a6b824f..87d49259e 100644 --- a/SRC/zlarfgp.f +++ b/SRC/zlarfgp.f @@ -116,8 +116,9 @@ SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU ) * ===================================================================== * * .. Parameters .. - DOUBLE PRECISION TWO, ONE, ZERO - PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 ) + DOUBLE PRECISION TWO, ONE, ZERO, HALF + PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0, + $ HALF = 0.5D+0 ) * .. * .. Local Scalars .. INTEGER J, KNT @@ -131,7 +132,7 @@ SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU ) $ ZLADIV * .. * .. Intrinsic Functions .. - INTRINSIC ABS, DBLE, DCMPLX, DIMAG, SIGN + INTRINSIC ABS, DBLE, DCMPLX, DIMAG, HUGE, SIGN * .. * .. External Subroutines .. EXTERNAL ZDSCAL, ZSCAL @@ -193,6 +194,18 @@ SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU ) XNORM = DZNRM2( N-1, X, INCX ) ALPHA = DCMPLX( ALPHR, ALPHI ) BETA = SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) + ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN +* +* |ALPHA| <= |BETA|, so ALPHA+BETA can overflow only when +* |BETA| > HUGE/2; scale X down and recompute them. +* + KNT = -1 + CALL ZDSCAL( N-1, SMLNUM, X, INCX ) + ALPHI = ALPHI*SMLNUM + ALPHR = ALPHR*SMLNUM + XNORM = DZNRM2( N-1, X, INCX ) + ALPHA = DCMPLX( ALPHR, ALPHI ) + BETA = SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) END IF SAVEALPHA = ALPHA ALPHA = ALPHA + BETA @@ -245,11 +258,14 @@ SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU ) * END IF * -* If BETA is subnormal, it may lose relative accuracy +* Undo the scaling. If BETA is subnormal, it may lose relative +* accuracy * DO 20 J = 1, KNT BETA = BETA*SMLNUM 20 CONTINUE + IF( KNT.LT.0 ) + $ BETA = BETA*BIGNUM ALPHA = BETA END IF * diff --git a/TESTING/LIN/alahd.f b/TESTING/LIN/alahd.f index b04a3f796..85d216b3a 100644 --- a/TESTING/LIN/alahd.f +++ b/TESTING/LIN/alahd.f @@ -945,7 +945,8 @@ SUBROUTINE ALAHD( IOUNIT, PATH ) $ '2. Upper triangular', 16X, '6. Random, CNDNUM = 0.1/EPS', $ / 4X, '3. Lower triangular', 16X, $ '7. Scaled near underflow', / 4X, '4. Random, CNDNUM = 2', - $ 14X, '8. Scaled near overflow' ) + $ 14X, '8. Scaled near overflow', / 39X, + $ '9. Leading entry near overflow' ) * * QP matrix types * diff --git a/TESTING/LIN/cchkaa.F b/TESTING/LIN/cchkaa.F index fa4000050..0c57bef29 100644 --- a/TESTING/LIN/cchkaa.F +++ b/TESTING/LIN/cchkaa.F @@ -1017,7 +1017,7 @@ PROGRAM CCHKAA * * QR: QR factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -1034,7 +1034,7 @@ PROGRAM CCHKAA * * LQ: LQ factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -1051,7 +1051,7 @@ PROGRAM CCHKAA * * QL: QL factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -1068,7 +1068,7 @@ PROGRAM CCHKAA * * RQ: RQ factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN diff --git a/TESTING/LIN/cchklq.f b/TESTING/LIN/cchklq.f index 175fd1a84..e9648ff25 100644 --- a/TESTING/LIN/cchklq.f +++ b/TESTING/LIN/cchklq.f @@ -219,9 +219,10 @@ SUBROUTINE CCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - REAL ZERO - PARAMETER ( ZERO = 0.0E0 ) + PARAMETER ( NTYPES = 9 ) + REAL ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, + $ QUARTER = 0.25E0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -235,6 +236,11 @@ SUBROUTINE CCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) REAL RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL SISNAN + REAL SLAMCH + EXTERNAL SISNAN, SLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, CERRLQ, CGELQF, CGELS, $ CGET02, CLACPY, CLARHS, CLATB4, CLATMS, CLQT01, @@ -314,6 +320,14 @@ SUBROUTINE CCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( 1 ) = ( ONE - QUARTER )*SLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of CLQT01; other values are * used in the calls of CLQT02, and must not exceed MINMN. @@ -423,7 +437,8 @@ SUBROUTINE CCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ SISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/cchkql.f b/TESTING/LIN/cchkql.f index 8b029263b..db3c2bf94 100644 --- a/TESTING/LIN/cchkql.f +++ b/TESTING/LIN/cchkql.f @@ -219,9 +219,10 @@ SUBROUTINE CCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - REAL ZERO - PARAMETER ( ZERO = 0.0E0 ) + PARAMETER ( NTYPES = 9 ) + REAL ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, + $ QUARTER = 0.25E0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -235,6 +236,11 @@ SUBROUTINE CCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) REAL RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL SISNAN + REAL SLAMCH + EXTERNAL SISNAN, SLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, CERRQL, CGEQLS, CGET02, $ CLACPY, CLARHS, CLATB4, CLATMS, CQLT01, CQLT02, @@ -314,6 +320,15 @@ SUBROUTINE CCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( ( N-1 )*LDA+M ) = ( ONE-QUARTER )* + $ SLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of CQLT01; other values are * used in the calls of CQLT02, and must not exceed MINMN. @@ -410,7 +425,8 @@ SUBROUTINE CCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ SISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/cchkqr.f b/TESTING/LIN/cchkqr.f index 2f3e64000..556ade52d 100644 --- a/TESTING/LIN/cchkqr.f +++ b/TESTING/LIN/cchkqr.f @@ -224,9 +224,10 @@ SUBROUTINE CCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 9 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - REAL ZERO - PARAMETER ( ZERO = 0.0E0 ) + PARAMETER ( NTYPES = 9 ) + REAL ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, + $ QUARTER = 0.25E0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -241,8 +242,11 @@ SUBROUTINE CCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, REAL RESULT( NTESTS ) * .. * .. External Functions .. + LOGICAL SISNAN + REAL SLAMCH LOGICAL CGENND EXTERNAL CGENND + EXTERNAL SISNAN, SLAMCH * .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, CERRQR, CGELS, CGET02, @@ -323,6 +327,14 @@ SUBROUTINE CCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( 1 ) = ( ONE - QUARTER )*SLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of CQRT01; other values are * used in the calls of CQRT02, and must not exceed MINMN. @@ -434,7 +446,8 @@ SUBROUTINE CCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NTESTS - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ SISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/cchkrq.f b/TESTING/LIN/cchkrq.f index 4d77d192f..4fb7f6c5a 100644 --- a/TESTING/LIN/cchkrq.f +++ b/TESTING/LIN/cchkrq.f @@ -224,9 +224,10 @@ SUBROUTINE CCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - REAL ZERO - PARAMETER ( ZERO = 0.0E0 ) + PARAMETER ( NTYPES = 9 ) + REAL ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, + $ QUARTER = 0.25E0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -240,6 +241,11 @@ SUBROUTINE CCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) REAL RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL SISNAN + REAL SLAMCH + EXTERNAL SISNAN, SLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, CERRRQ, CGERQS, CGET02, $ CLACPY, CLARHS, CLATB4, CLATMS, CRQT01, CRQT02, @@ -319,6 +325,15 @@ SUBROUTINE CCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( ( N-1 )*LDA+M ) = ( ONE-QUARTER )* + $ SLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of CRQT01; other values are * used in the calls of CRQT02, and must not exceed MINMN. @@ -415,7 +430,8 @@ SUBROUTINE CCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ SISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/cgeqls.f b/TESTING/LIN/cgeqls.f index 232b5dd23..da9e95c58 100644 --- a/TESTING/LIN/cgeqls.f +++ b/TESTING/LIN/cgeqls.f @@ -29,7 +29,8 @@ *> min || A*X - B || *> using the QL factorization *> A = Q*L -*> computed by CGEQLF. +*> computed by CGEQLF. Large right hand sides are scaled before +*> applying Q, and the solution is returned at the original scale. *> \endverbatim * * Arguments: @@ -139,11 +140,19 @@ SUBROUTINE CGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, COMPLEX ONE PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) ) * .. +* .. Local Scalars .. + REAL BIGNUM, BNRM + INTEGER I, J +* .. +* .. External Functions .. + REAL SLAMCH + EXTERNAL SLAMCH +* .. * .. External Subroutines .. - EXTERNAL CTRSM, CUNMQL, XERBLA + EXTERNAL CLASCL, CTRSM, CUNMQL, XERBLA * .. * .. Intrinsic Functions .. - INTRINSIC MAX + INTRINSIC ABS, AIMAG, MAX, REAL * .. * .. Executable Statements .. * @@ -174,6 +183,23 @@ SUBROUTINE CGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, IF( N.EQ.0 .OR. NRHS.EQ.0 .OR. M.EQ.0 ) $ RETURN * +* Scale large right hand sides before applying Q: tau*w can +* overflow even when Q' * B and the solution are representable. +* Use component magnitudes since ABS of a finite complex entry +* can exceed the overflow threshold. +* + BNRM = 0.0E+0 + DO 20 J = 1, NRHS + DO 10 I = 1, M + BNRM = MAX( BNRM, ABS( REAL( B( I, J ) ) ), + $ ABS( AIMAG( B( I, J ) ) ) ) + 10 CONTINUE + 20 CONTINUE + BIGNUM = SLAMCH( 'Precision' ) / SLAMCH( 'Safe minimum' ) + IF( BNRM.GT.BIGNUM ) + $ CALL CLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, + $ INFO ) +* * B := Q' * B * CALL CUNMQL( 'Left', 'Conjugate transpose', M, NRHS, N, A, LDA, @@ -183,6 +209,10 @@ SUBROUTINE CGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, * CALL CTRSM( 'Left', 'Lower', 'No transpose', 'Non-unit', N, NRHS, $ ONE, A( M-N+1, 1 ), LDA, B( M-N+1, 1 ), LDB ) +* + IF( BNRM.GT.BIGNUM ) + $ CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, M, NRHS, B, LDB, + $ INFO ) * RETURN * diff --git a/TESTING/LIN/dchkaa.F b/TESTING/LIN/dchkaa.F index 856185c65..09a332690 100644 --- a/TESTING/LIN/dchkaa.F +++ b/TESTING/LIN/dchkaa.F @@ -859,7 +859,7 @@ PROGRAM DCHKAA * * QR: QR factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -876,7 +876,7 @@ PROGRAM DCHKAA * * LQ: LQ factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -893,7 +893,7 @@ PROGRAM DCHKAA * * QL: QL factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -910,7 +910,7 @@ PROGRAM DCHKAA * * RQ: RQ factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN diff --git a/TESTING/LIN/dchklq.f b/TESTING/LIN/dchklq.f index 8fa7c8d6f..d297af90d 100644 --- a/TESTING/LIN/dchklq.f +++ b/TESTING/LIN/dchklq.f @@ -219,9 +219,10 @@ SUBROUTINE DCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - DOUBLE PRECISION ZERO - PARAMETER ( ZERO = 0.0D0 ) + PARAMETER ( NTYPES = 9 ) + DOUBLE PRECISION ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, + $ QUARTER = 0.25D0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -235,6 +236,11 @@ SUBROUTINE DCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) DOUBLE PRECISION RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION DLAMCH + EXTERNAL DISNAN, DLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, DERRLQ, DGELQF, DGELS, $ DGET02, DLACPY, DLARHS, DLATB4, DLATMS, DLQT01, @@ -314,6 +320,14 @@ SUBROUTINE DCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( 1 ) = ( ONE - QUARTER )*DLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of DLQT01; other values are * used in the calls of DLQT02, and must not exceed MINMN. @@ -433,7 +447,8 @@ SUBROUTINE DCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ DISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/dchkql.f b/TESTING/LIN/dchkql.f index 24837f1c4..0f907b827 100644 --- a/TESTING/LIN/dchkql.f +++ b/TESTING/LIN/dchkql.f @@ -219,9 +219,10 @@ SUBROUTINE DCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - DOUBLE PRECISION ZERO - PARAMETER ( ZERO = 0.0D0 ) + PARAMETER ( NTYPES = 9 ) + DOUBLE PRECISION ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, + $ QUARTER = 0.25D0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -235,6 +236,11 @@ SUBROUTINE DCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) DOUBLE PRECISION RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION DLAMCH + EXTERNAL DISNAN, DLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, DERRQL, DGEQLS, DGET02, $ DLACPY, DLARHS, DLATB4, DLATMS, DQLT01, DQLT02, @@ -314,6 +320,15 @@ SUBROUTINE DCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( ( N-1 )*LDA+M ) = ( ONE-QUARTER )* + $ DLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of DQLT01; other values are * used in the calls of DQLT02, and must not exceed MINMN. @@ -410,7 +425,8 @@ SUBROUTINE DCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ DISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/dchkqr.f b/TESTING/LIN/dchkqr.f index 013f7676a..157e9b5cb 100644 --- a/TESTING/LIN/dchkqr.f +++ b/TESTING/LIN/dchkqr.f @@ -224,9 +224,10 @@ SUBROUTINE DCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 9 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - DOUBLE PRECISION ZERO - PARAMETER ( ZERO = 0.0D0 ) + PARAMETER ( NTYPES = 9 ) + DOUBLE PRECISION ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, + $ QUARTER = 0.25D0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -241,8 +242,11 @@ SUBROUTINE DCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, DOUBLE PRECISION RESULT( NTESTS ) * .. * .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION DLAMCH LOGICAL DGENND EXTERNAL DGENND + EXTERNAL DISNAN, DLAMCH * .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, DERRQR, DGELS, DGET02, @@ -323,6 +327,14 @@ SUBROUTINE DCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( 1 ) = ( ONE - QUARTER )*DLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of DQRT01; other values are * used in the calls of DQRT02, and must not exceed MINMN. @@ -435,7 +447,8 @@ SUBROUTINE DCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NTESTS - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ DISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/dchkrq.f b/TESTING/LIN/dchkrq.f index e85ac113b..3d2d663d5 100644 --- a/TESTING/LIN/dchkrq.f +++ b/TESTING/LIN/dchkrq.f @@ -224,9 +224,10 @@ SUBROUTINE DCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - DOUBLE PRECISION ZERO - PARAMETER ( ZERO = 0.0D0 ) + PARAMETER ( NTYPES = 9 ) + DOUBLE PRECISION ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, + $ QUARTER = 0.25D0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -240,6 +241,11 @@ SUBROUTINE DCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) DOUBLE PRECISION RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION DLAMCH + EXTERNAL DISNAN, DLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, DERRRQ, DGERQS, DGET02, $ DLACPY, DLARHS, DLATB4, DLATMS, DRQT01, DRQT02, @@ -319,6 +325,15 @@ SUBROUTINE DCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( ( N-1 )*LDA+M ) = ( ONE-QUARTER )* + $ DLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of DRQT01; other values are * used in the calls of DRQT02, and must not exceed MINMN. @@ -416,7 +431,8 @@ SUBROUTINE DCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ DISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/dgeqls.f b/TESTING/LIN/dgeqls.f index 749acb1ea..2061e14ec 100644 --- a/TESTING/LIN/dgeqls.f +++ b/TESTING/LIN/dgeqls.f @@ -29,7 +29,8 @@ *> min || A*X - B || *> using the QL factorization *> A = Q*L -*> computed by DGEQLF. +*> computed by DGEQLF. Large right hand sides are scaled before +*> applying Q, and the solution is returned at the original scale. *> \endverbatim * * Arguments: @@ -139,8 +140,18 @@ SUBROUTINE DGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D+0 ) * .. +* .. Local Scalars .. + DOUBLE PRECISION BIGNUM, BNRM +* .. +* .. Local Arrays .. + DOUBLE PRECISION RWORK( 1 ) +* .. +* .. External Functions .. + DOUBLE PRECISION DLAMCH, DLANGE + EXTERNAL DLAMCH, DLANGE +* .. * .. External Subroutines .. - EXTERNAL DORMQL, DTRSM, XERBLA + EXTERNAL DLASCL, DORMQL, DTRSM, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -174,6 +185,15 @@ SUBROUTINE DGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, IF( N.EQ.0 .OR. NRHS.EQ.0 .OR. M.EQ.0 ) $ RETURN * +* Scale large right hand sides before applying Q: tau*w can +* overflow even when Q' * B and the solution are representable. +* + BNRM = DLANGE( 'M', M, NRHS, B, LDB, RWORK ) + BIGNUM = DLAMCH( 'Precision' ) / DLAMCH( 'Safe minimum' ) + IF( BNRM.GT.BIGNUM ) + $ CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, + $ INFO ) +* * B := Q' * B * CALL DORMQL( 'Left', 'Transpose', M, NRHS, N, A, LDA, TAU, B, LDB, @@ -183,6 +203,10 @@ SUBROUTINE DGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, * CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Non-unit', N, NRHS, $ ONE, A( M-N+1, 1 ), LDA, B( M-N+1, 1 ), LDB ) +* + IF( BNRM.GT.BIGNUM ) + $ CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, M, NRHS, B, LDB, + $ INFO ) * RETURN * diff --git a/TESTING/LIN/schkaa.F b/TESTING/LIN/schkaa.F index 25c455d4f..ae1325398 100644 --- a/TESTING/LIN/schkaa.F +++ b/TESTING/LIN/schkaa.F @@ -853,7 +853,7 @@ PROGRAM SCHKAA * * QR: QR factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -870,7 +870,7 @@ PROGRAM SCHKAA * * LQ: LQ factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -887,7 +887,7 @@ PROGRAM SCHKAA * * QL: QL factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -904,7 +904,7 @@ PROGRAM SCHKAA * * RQ: RQ factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN diff --git a/TESTING/LIN/schklq.f b/TESTING/LIN/schklq.f index 5dce71ecd..ae2a2ca23 100644 --- a/TESTING/LIN/schklq.f +++ b/TESTING/LIN/schklq.f @@ -219,9 +219,10 @@ SUBROUTINE SCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - REAL ZERO - PARAMETER ( ZERO = 0.0E0 ) + PARAMETER ( NTYPES = 9 ) + REAL ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, + $ QUARTER = 0.25E0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -235,6 +236,11 @@ SUBROUTINE SCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) REAL RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL SISNAN + REAL SLAMCH + EXTERNAL SISNAN, SLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, SERRLQ, SGELQF, SGELS, $ SGET02, SLACPY, SLARHS, SLATB4, SLATMS, SLQT01, @@ -314,6 +320,14 @@ SUBROUTINE SCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( 1 ) = ( ONE - QUARTER )*SLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of SLQT01; other values are * used in the calls of SLQT02, and must not exceed MINMN. @@ -423,7 +437,8 @@ SUBROUTINE SCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ SISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/schkql.f b/TESTING/LIN/schkql.f index 20726e188..7ea38fb0b 100644 --- a/TESTING/LIN/schkql.f +++ b/TESTING/LIN/schkql.f @@ -219,9 +219,10 @@ SUBROUTINE SCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - REAL ZERO - PARAMETER ( ZERO = 0.0E0 ) + PARAMETER ( NTYPES = 9 ) + REAL ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, + $ QUARTER = 0.25E0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -235,6 +236,11 @@ SUBROUTINE SCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) REAL RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL SISNAN + REAL SLAMCH + EXTERNAL SISNAN, SLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, SERRQL, SGEQLS, SGET02, $ SLACPY, SLARHS, SLATB4, SLATMS, SQLT01, SQLT02, @@ -314,6 +320,15 @@ SUBROUTINE SCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( ( N-1 )*LDA+M ) = ( ONE-QUARTER )* + $ SLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of SQLT01; other values are * used in the calls of SQLT02, and must not exceed MINMN. @@ -410,7 +425,8 @@ SUBROUTINE SCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ SISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/schkqr.f b/TESTING/LIN/schkqr.f index 49d775c17..a680e7b84 100644 --- a/TESTING/LIN/schkqr.f +++ b/TESTING/LIN/schkqr.f @@ -224,9 +224,10 @@ SUBROUTINE SCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 9 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - REAL ZERO - PARAMETER ( ZERO = 0.0E0 ) + PARAMETER ( NTYPES = 9 ) + REAL ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, + $ QUARTER = 0.25E0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -241,8 +242,11 @@ SUBROUTINE SCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, REAL RESULT( NTESTS ) * .. * .. External Functions .. + LOGICAL SISNAN + REAL SLAMCH LOGICAL SGENND EXTERNAL SGENND + EXTERNAL SISNAN, SLAMCH * .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, SERRQR, SGELS, SGET02, @@ -323,6 +327,14 @@ SUBROUTINE SCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( 1 ) = ( ONE - QUARTER )*SLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of SQRT01; other values are * used in the calls of SQRT02, and must not exceed MINMN. @@ -434,7 +446,8 @@ SUBROUTINE SCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NTESTS - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ SISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/schkrq.f b/TESTING/LIN/schkrq.f index 1d919f19c..5bd14fdc8 100644 --- a/TESTING/LIN/schkrq.f +++ b/TESTING/LIN/schkrq.f @@ -224,9 +224,10 @@ SUBROUTINE SCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - REAL ZERO - PARAMETER ( ZERO = 0.0E0 ) + PARAMETER ( NTYPES = 9 ) + REAL ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, + $ QUARTER = 0.25E0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -240,6 +241,11 @@ SUBROUTINE SCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) REAL RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL SISNAN + REAL SLAMCH + EXTERNAL SISNAN, SLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, SERRRQ, SGERQS, SGET02, $ SLACPY, SLARHS, SLATB4, SLATMS, SRQT01, SRQT02, @@ -319,6 +325,15 @@ SUBROUTINE SCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( ( N-1 )*LDA+M ) = ( ONE-QUARTER )* + $ SLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of SRQT01; other values are * used in the calls of SRQT02, and must not exceed MINMN. @@ -415,7 +430,8 @@ SUBROUTINE SCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ SISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/sgeqls.f b/TESTING/LIN/sgeqls.f index e89fb37b7..1baf89239 100644 --- a/TESTING/LIN/sgeqls.f +++ b/TESTING/LIN/sgeqls.f @@ -29,7 +29,8 @@ *> min || A*X - B || *> using the QL factorization *> A = Q*L -*> computed by SGEQLF. +*> computed by SGEQLF. Large right hand sides are scaled before +*> applying Q, and the solution is returned at the original scale. *> \endverbatim * * Arguments: @@ -139,8 +140,18 @@ SUBROUTINE SGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, REAL ONE PARAMETER ( ONE = 1.0E+0 ) * .. +* .. Local Scalars .. + REAL BIGNUM, BNRM +* .. +* .. Local Arrays .. + REAL RWORK( 1 ) +* .. +* .. External Functions .. + REAL SLAMCH, SLANGE + EXTERNAL SLAMCH, SLANGE +* .. * .. External Subroutines .. - EXTERNAL SORMQL, STRSM, XERBLA + EXTERNAL SLASCL, SORMQL, STRSM, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -174,6 +185,15 @@ SUBROUTINE SGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, IF( N.EQ.0 .OR. NRHS.EQ.0 .OR. M.EQ.0 ) $ RETURN * +* Scale large right hand sides before applying Q: tau*w can +* overflow even when Q' * B and the solution are representable. +* + BNRM = SLANGE( 'M', M, NRHS, B, LDB, RWORK ) + BIGNUM = SLAMCH( 'Precision' ) / SLAMCH( 'Safe minimum' ) + IF( BNRM.GT.BIGNUM ) + $ CALL SLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, + $ INFO ) +* * B := Q' * B * CALL SORMQL( 'Left', 'Transpose', M, NRHS, N, A, LDA, TAU, B, LDB, @@ -183,6 +203,10 @@ SUBROUTINE SGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, * CALL STRSM( 'Left', 'Lower', 'No transpose', 'Non-unit', N, NRHS, $ ONE, A( M-N+1, 1 ), LDA, B( M-N+1, 1 ), LDB ) +* + IF( BNRM.GT.BIGNUM ) + $ CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, M, NRHS, B, LDB, + $ INFO ) * RETURN * diff --git a/TESTING/LIN/zchkaa.F b/TESTING/LIN/zchkaa.F index a86793907..e9720a807 100644 --- a/TESTING/LIN/zchkaa.F +++ b/TESTING/LIN/zchkaa.F @@ -1016,7 +1016,7 @@ PROGRAM ZCHKAA * * QR: QR factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -1033,7 +1033,7 @@ PROGRAM ZCHKAA * * LQ: LQ factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -1050,7 +1050,7 @@ PROGRAM ZCHKAA * * QL: QL factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -1067,7 +1067,7 @@ PROGRAM ZCHKAA * * RQ: RQ factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN diff --git a/TESTING/LIN/zchklq.f b/TESTING/LIN/zchklq.f index d4106fa4c..1305ba8ca 100644 --- a/TESTING/LIN/zchklq.f +++ b/TESTING/LIN/zchklq.f @@ -219,9 +219,10 @@ SUBROUTINE ZCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - DOUBLE PRECISION ZERO - PARAMETER ( ZERO = 0.0D0 ) + PARAMETER ( NTYPES = 9 ) + DOUBLE PRECISION ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, + $ QUARTER = 0.25D0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -235,6 +236,11 @@ SUBROUTINE ZCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) DOUBLE PRECISION RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION DLAMCH + EXTERNAL DISNAN, DLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, ZERRLQ, ZGELQF, ZGELS, $ ZGET02, ZLACPY, ZLARHS, ZLATB4, ZLATMS, ZLQT01, @@ -314,6 +320,14 @@ SUBROUTINE ZCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( 1 ) = ( ONE - QUARTER )*DLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of ZLQT01; other values are * used in the calls of ZLQT02, and must not exceed MINMN. @@ -423,7 +437,8 @@ SUBROUTINE ZCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ DISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/zchkql.f b/TESTING/LIN/zchkql.f index cafb2c623..3f8b005d4 100644 --- a/TESTING/LIN/zchkql.f +++ b/TESTING/LIN/zchkql.f @@ -219,9 +219,10 @@ SUBROUTINE ZCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - DOUBLE PRECISION ZERO - PARAMETER ( ZERO = 0.0D0 ) + PARAMETER ( NTYPES = 9 ) + DOUBLE PRECISION ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, + $ QUARTER = 0.25D0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -235,6 +236,11 @@ SUBROUTINE ZCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) DOUBLE PRECISION RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION DLAMCH + EXTERNAL DISNAN, DLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, ZERRQL, ZGEQLS, $ ZGET02, ZLACPY, ZLARHS, ZLATB4, ZLATMS, ZQLT01, @@ -314,6 +320,15 @@ SUBROUTINE ZCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( ( N-1 )*LDA+M ) = ( ONE-QUARTER )* + $ DLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of ZQLT01; other values are * used in the calls of ZQLT02, and must not exceed MINMN. @@ -410,7 +425,8 @@ SUBROUTINE ZCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ DISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/zchkqr.f b/TESTING/LIN/zchkqr.f index 740abf4df..6cb7d2c39 100644 --- a/TESTING/LIN/zchkqr.f +++ b/TESTING/LIN/zchkqr.f @@ -224,9 +224,10 @@ SUBROUTINE ZCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 9 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - DOUBLE PRECISION ZERO - PARAMETER ( ZERO = 0.0D0 ) + PARAMETER ( NTYPES = 9 ) + DOUBLE PRECISION ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, + $ QUARTER = 0.25D0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -241,8 +242,11 @@ SUBROUTINE ZCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, DOUBLE PRECISION RESULT( NTESTS ) * .. * .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION DLAMCH LOGICAL ZGENND EXTERNAL ZGENND + EXTERNAL DISNAN, DLAMCH * .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, ZERRQR, ZGELS, @@ -323,6 +327,14 @@ SUBROUTINE ZCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( 1 ) = ( ONE - QUARTER )*DLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of ZQRT01; other values are * used in the calls of ZQRT02, and must not exceed MINMN. @@ -434,7 +446,8 @@ SUBROUTINE ZCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NTESTS - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ DISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/zchkrq.f b/TESTING/LIN/zchkrq.f index 6bcb3d6b9..7e2ccc80a 100644 --- a/TESTING/LIN/zchkrq.f +++ b/TESTING/LIN/zchkrq.f @@ -224,9 +224,10 @@ SUBROUTINE ZCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - DOUBLE PRECISION ZERO - PARAMETER ( ZERO = 0.0D0 ) + PARAMETER ( NTYPES = 9 ) + DOUBLE PRECISION ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, + $ QUARTER = 0.25D0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -240,6 +241,11 @@ SUBROUTINE ZCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) DOUBLE PRECISION RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION DLAMCH + EXTERNAL DISNAN, DLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, ZERRRQ, ZGERQS, $ ZGET02, ZLACPY, ZLARHS, ZLATB4, ZLATMS, ZRQT01, @@ -319,6 +325,15 @@ SUBROUTINE ZCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( ( N-1 )*LDA+M ) = ( ONE-QUARTER )* + $ DLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of ZRQT01; other values are * used in the calls of ZRQT02, and must not exceed MINMN. @@ -415,7 +430,8 @@ SUBROUTINE ZCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ DISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/zgeqls.f b/TESTING/LIN/zgeqls.f index 5d0eb55e8..aa92bb044 100644 --- a/TESTING/LIN/zgeqls.f +++ b/TESTING/LIN/zgeqls.f @@ -29,7 +29,8 @@ *> min || A*X - B || *> using the QL factorization *> A = Q*L -*> computed by ZGEQLF. +*> computed by ZGEQLF. Large right hand sides are scaled before +*> applying Q, and the solution is returned at the original scale. *> \endverbatim * * Arguments: @@ -139,11 +140,19 @@ SUBROUTINE ZGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, COMPLEX*16 ONE PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) * .. +* .. Local Scalars .. + DOUBLE PRECISION BIGNUM, BNRM + INTEGER I, J +* .. +* .. External Functions .. + DOUBLE PRECISION DLAMCH + EXTERNAL DLAMCH +* .. * .. External Subroutines .. - EXTERNAL XERBLA, ZTRSM, ZUNMQL + EXTERNAL XERBLA, ZLASCL, ZTRSM, ZUNMQL * .. * .. Intrinsic Functions .. - INTRINSIC MAX + INTRINSIC ABS, DBLE, DIMAG, MAX * .. * .. Executable Statements .. * @@ -174,6 +183,23 @@ SUBROUTINE ZGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, IF( N.EQ.0 .OR. NRHS.EQ.0 .OR. M.EQ.0 ) $ RETURN * +* Scale large right hand sides before applying Q: tau*w can +* overflow even when Q' * B and the solution are representable. +* Use component magnitudes since ABS of a finite complex entry +* can exceed the overflow threshold. +* + BNRM = 0.0D+0 + DO 20 J = 1, NRHS + DO 10 I = 1, M + BNRM = MAX( BNRM, ABS( DBLE( B( I, J ) ) ), + $ ABS( DIMAG( B( I, J ) ) ) ) + 10 CONTINUE + 20 CONTINUE + BIGNUM = DLAMCH( 'Precision' ) / DLAMCH( 'Safe minimum' ) + IF( BNRM.GT.BIGNUM ) + $ CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, + $ INFO ) +* * B := Q' * B * CALL ZUNMQL( 'Left', 'Conjugate transpose', M, NRHS, N, A, LDA, @@ -183,6 +209,10 @@ SUBROUTINE ZGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, * CALL ZTRSM( 'Left', 'Lower', 'No transpose', 'Non-unit', N, NRHS, $ ONE, A( M-N+1, 1 ), LDA, B( M-N+1, 1 ), LDB ) +* + IF( BNRM.GT.BIGNUM ) + $ CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, M, NRHS, B, LDB, + $ INFO ) * RETURN * diff --git a/TESTING/ctest.in b/TESTING/ctest.in index 4e30224d7..a16b56295 100644 --- a/TESTING/ctest.in +++ b/TESTING/ctest.in @@ -37,10 +37,10 @@ CSP 11 List types on next line if 0 < NTYPES < 11 CTR 18 List types on next line if 0 < NTYPES < 18 CTP 18 List types on next line if 0 < NTYPES < 18 CTB 17 List types on next line if 0 < NTYPES < 17 -CQR 8 List types on next line if 0 < NTYPES < 8 -CRQ 8 List types on next line if 0 < NTYPES < 8 -CLQ 8 List types on next line if 0 < NTYPES < 8 -CQL 8 List types on next line if 0 < NTYPES < 8 +CQR 9 List types on next line if 0 < NTYPES < 9 +CRQ 9 List types on next line if 0 < NTYPES < 9 +CLQ 9 List types on next line if 0 < NTYPES < 9 +CQL 9 List types on next line if 0 < NTYPES < 9 CQP 6 List types on next line if 0 < NTYPES < 6 CQK 19 List types on next line if 0 < NTYPES < 19 CCX 19 LIst types on next line if 0 < NTYPES < 19 diff --git a/TESTING/dtest.in b/TESTING/dtest.in index cde62db50..2c5ce7659 100644 --- a/TESTING/dtest.in +++ b/TESTING/dtest.in @@ -31,10 +31,10 @@ DSP 10 List types on next line if 0 < NTYPES < 10 DTR 18 List types on next line if 0 < NTYPES < 18 DTP 18 List types on next line if 0 < NTYPES < 18 DTB 17 List types on next line if 0 < NTYPES < 17 -DQR 8 List types on next line if 0 < NTYPES < 8 -DRQ 8 List types on next line if 0 < NTYPES < 8 -DLQ 8 List types on next line if 0 < NTYPES < 8 -DQL 8 List types on next line if 0 < NTYPES < 8 +DQR 9 List types on next line if 0 < NTYPES < 9 +DRQ 9 List types on next line if 0 < NTYPES < 9 +DLQ 9 List types on next line if 0 < NTYPES < 9 +DQL 9 List types on next line if 0 < NTYPES < 9 DQP 6 List types on next line if 0 < NTYPES < 6 DQK 19 LIst types on next line if 0 < NTYPES < 19 DCX 19 LIst types on next line if 0 < NTYPES < 19 diff --git a/TESTING/stest.in b/TESTING/stest.in index abfd639fd..b6a859ef1 100644 --- a/TESTING/stest.in +++ b/TESTING/stest.in @@ -31,10 +31,10 @@ SSP 10 List types on next line if 0 < NTYPES < 10 STR 18 List types on next line if 0 < NTYPES < 18 STP 18 List types on next line if 0 < NTYPES < 18 STB 17 List types on next line if 0 < NTYPES < 17 -SQR 8 List types on next line if 0 < NTYPES < 8 -SRQ 8 List types on next line if 0 < NTYPES < 8 -SLQ 8 List types on next line if 0 < NTYPES < 8 -SQL 8 List types on next line if 0 < NTYPES < 8 +SQR 9 List types on next line if 0 < NTYPES < 9 +SRQ 9 List types on next line if 0 < NTYPES < 9 +SLQ 9 List types on next line if 0 < NTYPES < 9 +SQL 9 List types on next line if 0 < NTYPES < 9 SQP 6 List types on next line if 0 < NTYPES < 6 SQK 19 List types on next line if 0 < NTYPES < 19 SCX 19 LIst types on next line if 0 < NTYPES < 19 diff --git a/TESTING/ztest.in b/TESTING/ztest.in index bf4c9d100..f123e6630 100644 --- a/TESTING/ztest.in +++ b/TESTING/ztest.in @@ -37,10 +37,10 @@ ZSP 11 List types on next line if 0 < NTYPES < 11 ZTR 18 List types on next line if 0 < NTYPES < 18 ZTP 18 List types on next line if 0 < NTYPES < 18 ZTB 17 List types on next line if 0 < NTYPES < 17 -ZQR 8 List types on next line if 0 < NTYPES < 8 -ZRQ 8 List types on next line if 0 < NTYPES < 8 -ZLQ 8 List types on next line if 0 < NTYPES < 8 -ZQL 8 List types on next line if 0 < NTYPES < 8 +ZQR 9 List types on next line if 0 < NTYPES < 9 +ZRQ 9 List types on next line if 0 < NTYPES < 9 +ZLQ 9 List types on next line if 0 < NTYPES < 9 +ZQL 9 List types on next line if 0 < NTYPES < 9 ZQP 6 List types on next line if 0 < NTYPES < 6 ZQK 19 List types on next line if 0 < NTYPES < 19 ZCX 19 List types on next line if 0 < NTYPES < 19