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/CMakeLists.txt b/TESTING/LIN/CMakeLists.txt index 2313fa0c4..ec9bb134a 100644 --- a/TESTING/LIN/CMakeLists.txt +++ b/TESTING/LIN/CMakeLists.txt @@ -257,7 +257,7 @@ set(ZLINTSTRFP chkxer.f xerbla.f alaerh.f aladhd.f alahd.f alasvm.f) function(add_lin_executable name) - cmake_parse_arguments(PARSE_ARGV 1 LIN "" "" "SOURCES;DEFAULT_API;EXT_API") + cmake_parse_arguments(PARSE_ARGV 1 LIN "TEST" "" "SOURCES;DEFAULT_API;EXT_API") set(sources ${LIN_SOURCES} ${LIN_DEFAULT_API}) set(ext_sources ${LIN_SOURCES} ${LIN_EXT_API}) @@ -265,6 +265,9 @@ function(add_lin_executable name) add_executable(${name} ${sources}) target_link_libraries(${name} PRIVATE ${TMGLIB} ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES}) lapack_add_coverage(${name}) + if(LIN_TEST) + add_test(NAME LAPACK-${name} COMMAND $) + endif() endif() if(BUILD_INDEX64_EXT_API) @@ -276,6 +279,9 @@ function(add_lin_executable name) target_compile_options(${name}_64 PRIVATE ${FOPT_ILP64}) target_link_libraries(${name}_64 PRIVATE ${TMGLIB} ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES}) lapack_add_coverage(${name}_64) + if(LIN_TEST) + add_test(NAME LAPACK-${name}_64 COMMAND $) + endif() # Add depedency to the global codegen target. Since we reuse the same # generated source file for multiple tests, the generation could be @@ -310,6 +316,8 @@ if(BUILD_COMPLEX) DEFAULT_API ${CXLINTST} EXT_API ${CLINTST_NO_XBLAS}) add_lin_executable(xlintstrfc SOURCES ${CLINTSTRFP}) + add_lin_executable(xrq_helpers_c TEST + SOURCES ctest_rq_helpers.f90 cgerqs.f cget02.f) endif() if(BUILD_COMPLEX16) @@ -318,6 +326,8 @@ if(BUILD_COMPLEX16) DEFAULT_API ${ZXLINTST} EXT_API ${ZLINTST_NO_XBLAS}) add_lin_executable(xlintstrfz SOURCES ${ZLINTSTRFP}) + add_lin_executable(xrq_helpers_z TEST + SOURCES ztest_rq_helpers.f90 zgerqs.f zget02.f) endif() if(BUILD_COMPLEX AND BUILD_COMPLEX16) 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 ac181df8c..a46fe949c 100644 --- a/TESTING/LIN/cchkaa.F +++ b/TESTING/LIN/cchkaa.F @@ -1038,7 +1038,7 @@ PROGRAM CCHKAA * * QR: QR factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -1055,7 +1055,7 @@ PROGRAM CCHKAA * * LQ: LQ factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -1072,7 +1072,7 @@ PROGRAM CCHKAA * * QL: QL factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -1089,7 +1089,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/cgerqs.f b/TESTING/LIN/cgerqs.f index ac28a3643..0a48f6252 100644 --- a/TESTING/LIN/cgerqs.f +++ b/TESTING/LIN/cgerqs.f @@ -29,7 +29,8 @@ *> min || A*X - B || *> using the RQ factorization *> A = R*Q -*> computed by CGERQF. +*> computed by CGERQF. The triangular solve uses scaled +*> complex division to avoid intermediate overflow. *> \endverbatim * * Arguments: @@ -136,12 +137,18 @@ SUBROUTINE CGERQS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, * ===================================================================== * * .. Parameters .. - COMPLEX CZERO, CONE - PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), - $ CONE = ( 1.0E+0, 0.0E+0 ) ) + COMPLEX CZERO + PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) ) +* .. +* .. Local Scalars .. + INTEGER I, J +* .. +* .. External Functions .. + COMPLEX CLADIV + EXTERNAL CLADIV * .. * .. External Subroutines .. - EXTERNAL CLASET, CTRSM, CUNMRQ, XERBLA + EXTERNAL CAXPY, CLASET, CUNMRQ, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -177,8 +184,17 @@ SUBROUTINE CGERQS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, * * Solve R*X = B(n-m+1:n,:) * - CALL CTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', M, NRHS, - $ CONE, A( 1, N-M+1 ), LDA, B( N-M+1, 1 ), LDB ) +* LADIV avoids overflow in intrinsic complex division. Scaling +* all of B instead could erase small solution components. +* + DO 20 J = 1, NRHS + DO 10 I = M, 1, -1 + B( N-M+I, J ) = CLADIV( B( N-M+I, J ), + $ A( I, N-M+I ) ) + CALL CAXPY( I-1, -B( N-M+I, J ), A( 1, N-M+I ), + $ 1, B( N-M+1, J ), 1 ) + 10 CONTINUE + 20 CONTINUE * * Set B(1:n-m,:) to zero * diff --git a/TESTING/LIN/cget02.f b/TESTING/LIN/cget02.f index cb92c3553..3e4c7e16a 100644 --- a/TESTING/LIN/cget02.f +++ b/TESTING/LIN/cget02.f @@ -168,7 +168,7 @@ SUBROUTINE CGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, EXTERNAL CGEMM * .. * .. Intrinsic Functions .. - INTRINSIC MAX + INTRINSIC HUGE, MAX * .. * .. Executable Statements .. * @@ -187,7 +187,7 @@ SUBROUTINE CGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, N2 = N END IF * -* Exit with RESID = 1/EPS if ANORM = 0. +* Exit with RESID = 1/EPS if ANORM is zero or nonfinite. * EPS = SLAMCH( 'Epsilon' ) IF( LSAME( TRANS, 'N' ) ) THEN @@ -195,7 +195,8 @@ SUBROUTINE CGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, ELSE ANORM = CLANGE( 'I', M, N, A, LDA, RWORK ) END IF - IF( ANORM.LE.ZERO ) THEN + IF( .NOT.( ANORM.GT.ZERO .AND. + $ ANORM.LE.HUGE( ANORM ) ) ) THEN RESID = ONE / EPS RETURN END IF @@ -212,8 +213,12 @@ SUBROUTINE CGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, DO 10 J = 1, NRHS BNORM = SCASUM( N1, B( 1, J ), 1 ) XNORM = SCASUM( N2, X( 1, J ), 1 ) - IF( XNORM.LE.ZERO ) THEN +* Reject nonfinite norms before MAX can hide a NaN ratio. + IF( .NOT.( XNORM.GT.ZERO .AND. + $ XNORM.LE.HUGE( XNORM ) .AND. + $ BNORM.LE.HUGE( BNORM ) ) ) THEN RESID = ONE / EPS + RETURN ELSE RESID = MAX( RESID, ( ( BNORM/ANORM )/XNORM )/EPS ) END IF diff --git a/TESTING/LIN/ctest_rq_helpers.f90 b/TESTING/LIN/ctest_rq_helpers.f90 new file mode 100644 index 000000000..ca77488cb --- /dev/null +++ b/TESTING/LIN/ctest_rq_helpers.f90 @@ -0,0 +1,68 @@ +! SPDX-FileCopyrightText: The LAPACK Authors +! SPDX-License-Identifier: BSD-3-Clause +! +! Regression tests for the RQ solve and residual test helpers. +program ctest_rq_helpers + use, intrinsic :: ieee_arithmetic + implicit none + integer, parameter :: wp = kind(1.0) + complex(wp) :: a(2,3), b(3,2), x(3,2), tau(2), work(128) + real(wp) :: large, small, resid, rwork(3), bad(2) + integer :: info, i, j, k + external :: cgerqs, cget02 + character :: trans(3) = ['N', 'T', 'C'] + + ! R has both ends of the normal range in one right-hand side. + ! Uniform scaling of B would lose the first solution component. + large = 0.75_wp*huge(1.0_wp) + small = tiny(1.0_wp) + a = 0 + a(1,2) = small + a(2,3) = large + tau = 0 + b = 0 + b(2,1) = small + b(3,1) = cmplx(-0.95_wp*large, 0.95_wp*large, wp) + b(2,2) = 2*small + b(3,2) = 0.5_wp*large + x = 0 + x(2,1) = 1 + x(3,1) = cmplx(-0.95_wp, 0.95_wp, wp) + x(2,2) = 2 + x(3,2) = 0.5_wp + call cgerqs(2, 3, 2, a, 2, tau, b, 3, work, 128, info) + if (info /= 0) stop 1 + if (.not. all(abs(b-x) <= 8*epsilon(1.0_wp))) stop 2 + + ! Check ordinary magnitudes as well. + a(1,2) = 2 + a(2,3) = 4 + b = 0 + b(2,:) = 2*x(2,:) + b(3,:) = 4*x(3,:) + call cgerqs(2, 3, 2, a, 2, tau, b, 3, work, 128, info) + if (info /= 0) stop 3 + if (.not. all(abs(b-x) <= 8*epsilon(1.0_wp))) stop 4 + + bad = [ieee_value(0.0_wp, ieee_quiet_nan), & + ieee_value(0.0_wp, ieee_positive_inf)] + do k = 1, 3 + do j = 1, 3 + do i = 1, 2 + a = 1 + x = 1 + b = 1 + if (j == 1) a(1,1) = bad(i) + if (j == 2) x(1,1) = bad(i) + if (j == 3) b(1,1) = bad(i) + call cget02(trans(k), 1, 1, 1, a, 2, x, 3, b, 3, rwork, resid) + if (.not. (resid > 30)) stop 5 + end do + end do + a = 1 + x = 1 + b = 1 + call cget02(trans(k), 1, 1, 1, a, 2, x, 3, b, 3, rwork, resid) + if (resid /= 0) stop 6 + end do +end program diff --git a/TESTING/LIN/dchkaa.F b/TESTING/LIN/dchkaa.F index d0f7dbcb4..1411d5633 100644 --- a/TESTING/LIN/dchkaa.F +++ b/TESTING/LIN/dchkaa.F @@ -880,7 +880,7 @@ PROGRAM DCHKAA * * QR: QR factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -897,7 +897,7 @@ PROGRAM DCHKAA * * LQ: LQ factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -914,7 +914,7 @@ PROGRAM DCHKAA * * QL: QL factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -931,7 +931,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 dd3f8de4b..d7a0390f5 100644 --- a/TESTING/LIN/schkaa.F +++ b/TESTING/LIN/schkaa.F @@ -874,7 +874,7 @@ PROGRAM SCHKAA * * QR: QR factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -891,7 +891,7 @@ PROGRAM SCHKAA * * LQ: LQ factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -908,7 +908,7 @@ PROGRAM SCHKAA * * QL: QL factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -925,7 +925,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 fb1b5cd59..b64e3f31d 100644 --- a/TESTING/LIN/zchkaa.F +++ b/TESTING/LIN/zchkaa.F @@ -1037,7 +1037,7 @@ PROGRAM ZCHKAA * * QR: QR factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -1054,7 +1054,7 @@ PROGRAM ZCHKAA * * LQ: LQ factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -1071,7 +1071,7 @@ PROGRAM ZCHKAA * * QL: QL factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -1088,7 +1088,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/LIN/zgerqs.f b/TESTING/LIN/zgerqs.f index ce3021284..2ff05346e 100644 --- a/TESTING/LIN/zgerqs.f +++ b/TESTING/LIN/zgerqs.f @@ -29,7 +29,8 @@ *> min || A*X - B || *> using the RQ factorization *> A = R*Q -*> computed by ZGERQF. +*> computed by ZGERQF. The triangular solve uses scaled +*> complex division to avoid intermediate overflow. *> \endverbatim * * Arguments: @@ -136,12 +137,18 @@ SUBROUTINE ZGERQS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, * ===================================================================== * * .. Parameters .. - COMPLEX*16 CZERO, CONE - PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), - $ CONE = ( 1.0D+0, 0.0D+0 ) ) + COMPLEX*16 CZERO + PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) ) +* .. +* .. Local Scalars .. + INTEGER I, J +* .. +* .. External Functions .. + COMPLEX*16 ZLADIV + EXTERNAL ZLADIV * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZLASET, ZTRSM, ZUNMRQ + EXTERNAL XERBLA, ZAXPY, ZLASET, ZUNMRQ * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -177,8 +184,17 @@ SUBROUTINE ZGERQS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, * * Solve R*X = B(n-m+1:n,:) * - CALL ZTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', M, NRHS, - $ CONE, A( 1, N-M+1 ), LDA, B( N-M+1, 1 ), LDB ) +* LADIV avoids overflow in intrinsic complex division. Scaling +* all of B instead could erase small solution components. +* + DO 20 J = 1, NRHS + DO 10 I = M, 1, -1 + B( N-M+I, J ) = ZLADIV( B( N-M+I, J ), + $ A( I, N-M+I ) ) + CALL ZAXPY( I-1, -B( N-M+I, J ), A( 1, N-M+I ), + $ 1, B( N-M+1, J ), 1 ) + 10 CONTINUE + 20 CONTINUE * * Set B(1:n-m,:) to zero * diff --git a/TESTING/LIN/zget02.f b/TESTING/LIN/zget02.f index 8f45177b5..ffe41aba0 100644 --- a/TESTING/LIN/zget02.f +++ b/TESTING/LIN/zget02.f @@ -168,7 +168,7 @@ SUBROUTINE ZGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, EXTERNAL ZGEMM * .. * .. Intrinsic Functions .. - INTRINSIC MAX + INTRINSIC HUGE, MAX * .. * .. Executable Statements .. * @@ -187,7 +187,7 @@ SUBROUTINE ZGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, N2 = N END IF * -* Exit with RESID = 1/EPS if ANORM = 0. +* Exit with RESID = 1/EPS if ANORM is zero or nonfinite. * EPS = DLAMCH( 'Epsilon' ) IF( LSAME( TRANS, 'N' ) ) THEN @@ -195,7 +195,8 @@ SUBROUTINE ZGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, ELSE ANORM = ZLANGE( 'I', M, N, A, LDA, RWORK ) END IF - IF( ANORM.LE.ZERO ) THEN + IF( .NOT.( ANORM.GT.ZERO .AND. + $ ANORM.LE.HUGE( ANORM ) ) ) THEN RESID = ONE / EPS RETURN END IF @@ -212,8 +213,12 @@ SUBROUTINE ZGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, DO 10 J = 1, NRHS BNORM = DZASUM( N1, B( 1, J ), 1 ) XNORM = DZASUM( N2, X( 1, J ), 1 ) - IF( XNORM.LE.ZERO ) THEN +* Reject nonfinite norms before MAX can hide a NaN ratio. + IF( .NOT.( XNORM.GT.ZERO .AND. + $ XNORM.LE.HUGE( XNORM ) .AND. + $ BNORM.LE.HUGE( BNORM ) ) ) THEN RESID = ONE / EPS + RETURN ELSE RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS ) END IF diff --git a/TESTING/LIN/ztest_rq_helpers.f90 b/TESTING/LIN/ztest_rq_helpers.f90 new file mode 100644 index 000000000..2030c8adb --- /dev/null +++ b/TESTING/LIN/ztest_rq_helpers.f90 @@ -0,0 +1,68 @@ +! SPDX-FileCopyrightText: The LAPACK Authors +! SPDX-License-Identifier: BSD-3-Clause +! +! Regression tests for the RQ solve and residual test helpers. +program ztest_rq_helpers + use, intrinsic :: ieee_arithmetic + implicit none + integer, parameter :: wp = kind(1.0d0) + complex(wp) :: a(2,3), b(3,2), x(3,2), tau(2), work(128) + real(wp) :: large, small, resid, rwork(3), bad(2) + integer :: info, i, j, k + external :: zgerqs, zget02 + character :: trans(3) = ['N', 'T', 'C'] + + ! R has both ends of the normal range in one right-hand side. + ! Uniform scaling of B would lose the first solution component. + large = 0.75_wp*huge(1.0_wp) + small = tiny(1.0_wp) + a = 0 + a(1,2) = small + a(2,3) = large + tau = 0 + b = 0 + b(2,1) = small + b(3,1) = cmplx(-0.95_wp*large, 0.95_wp*large, wp) + b(2,2) = 2*small + b(3,2) = 0.5_wp*large + x = 0 + x(2,1) = 1 + x(3,1) = cmplx(-0.95_wp, 0.95_wp, wp) + x(2,2) = 2 + x(3,2) = 0.5_wp + call zgerqs(2, 3, 2, a, 2, tau, b, 3, work, 128, info) + if (info /= 0) stop 1 + if (.not. all(abs(b-x) <= 8*epsilon(1.0_wp))) stop 2 + + ! Check ordinary magnitudes as well. + a(1,2) = 2 + a(2,3) = 4 + b = 0 + b(2,:) = 2*x(2,:) + b(3,:) = 4*x(3,:) + call zgerqs(2, 3, 2, a, 2, tau, b, 3, work, 128, info) + if (info /= 0) stop 3 + if (.not. all(abs(b-x) <= 8*epsilon(1.0_wp))) stop 4 + + bad = [ieee_value(0.0_wp, ieee_quiet_nan), & + ieee_value(0.0_wp, ieee_positive_inf)] + do k = 1, 3 + do j = 1, 3 + do i = 1, 2 + a = 1 + x = 1 + b = 1 + if (j == 1) a(1,1) = bad(i) + if (j == 2) x(1,1) = bad(i) + if (j == 3) b(1,1) = bad(i) + call zget02(trans(k), 1, 1, 1, a, 2, x, 3, b, 3, rwork, resid) + if (.not. (resid > 30)) stop 5 + end do + end do + a = 1 + x = 1 + b = 1 + call zget02(trans(k), 1, 1, 1, a, 2, x, 3, b, 3, rwork, resid) + if (resid /= 0) stop 6 + end do +end program 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