From 26df4b359992758bea3cf3acbb39778c8d17212e Mon Sep 17 00:00:00 2001 From: Angelika Schwarz Date: Wed, 1 Nov 2023 17:13:10 +0100 Subject: [PATCH] Improve testing of NRM2 Add tests covering arrays with extreme numbers provided in Anderson E. (2017) Algorithm 978: Safe Scaling in the Level 1 BLAS ACM Trans Math Softw 44:1--28 https://doi.org/10.1145/3061665 [D,S]NRM2 coverage improves lines 59.1% -> 100.0% branches 30.0% -> 90.0%. [DZ,SC]NRM2 coverage improves lines 67.0% -> 100.0% branches 41.7% -> 91.7% --- BLAS/TESTING/cblat1.f | 234 +++++++++++++++++++++++++++++++++++++++++- BLAS/TESTING/dblat1.f | 219 ++++++++++++++++++++++++++++++++++++++- BLAS/TESTING/sblat1.f | 220 ++++++++++++++++++++++++++++++++++++++- BLAS/TESTING/zblat1.f | 234 +++++++++++++++++++++++++++++++++++++++++- 4 files changed, 899 insertions(+), 8 deletions(-) diff --git a/BLAS/TESTING/cblat1.f b/BLAS/TESTING/cblat1.f index c6dc453b45..83b02f0cac 100644 --- a/BLAS/TESTING/cblat1.f +++ b/BLAS/TESTING/cblat1.f @@ -121,7 +121,8 @@ SUBROUTINE HEADER SUBROUTINE CHECK1(SFAC) * .. Parameters .. INTEGER NOUT - PARAMETER (NOUT=6) + REAL THRESH + PARAMETER (NOUT=6, THRESH=10.0E0) * .. Scalar Arguments .. REAL SFAC * .. Scalars in Common .. @@ -141,7 +142,7 @@ SUBROUTINE CHECK1(SFAC) INTEGER ICAMAX EXTERNAL SCASUM, SCNRM2, ICAMAX * .. External Subroutines .. - EXTERNAL CSCAL, CSSCAL, CTEST, ITEST1, STEST1 + EXTERNAL CB1NRM2, CSCAL, CSSCAL, CTEST, ITEST1, STEST1 * .. Intrinsic Functions .. INTRINSIC MAX * .. Common blocks .. @@ -256,6 +257,10 @@ SUBROUTINE CHECK1(SFAC) 20 CONTINUE IF (ICASE.EQ.6) THEN * .. SCNRM2 .. +* Test scaling when some entries are tiny or huge + CALL CB1NRM2(N,(INCX-2)*2,THRESH) + CALL CB1NRM2(N,INCX,THRESH) +* Test with hardcoded mid range entries CALL STEST1(SCNRM2(N,CX,INCX),STRUE2(NP1),STRUE2(NP1), + SFAC) ELSE IF (ICASE.EQ.7) THEN @@ -782,3 +787,228 @@ SUBROUTINE ITEST1(ICOMP,ITRUE) * End of ITEST1 * END + SUBROUTINE CB1NRM2(N,INCX,THRESH) +* Compare NRM2 with a reference computation using combinations +* of the following values: +* +* 0, very small, small, ulp, 1, 1/ulp, big, very big, infinity, NaN +* +* one of these values is used to initialize x(1) and x(2:N) is +* filled with random values from [-1,1] scaled by another of +* these values. +* +* This routine is adapted from the test suite provided by +* Anderson E. (2017) +* Algorithm 978: Safe Scaling in the Level 1 BLAS +* ACM Trans Math Softw 44:1--28 +* https://doi.org/10.1145/3061665 +* +* .. Scalar Arguments .. + INTEGER INCX, N + REAL THRESH +* +* ===================================================================== +* .. Parameters .. + INTEGER NMAX, NOUT, NV + PARAMETER (NMAX=20, NOUT=6, NV=10) + REAL HALF, ONE, THREE, TWO, ZERO + PARAMETER (HALF=0.5E+0, ONE=1.0E+0, TWO= 2.0E+0, + & THREE=3.0E+0, ZERO=0.0E+0) +* .. External Functions .. + REAL SCNRM2 + EXTERNAL SCNRM2 +* .. Intrinsic Functions .. + INTRINSIC AIMAG, ABS, CMPLX, MAX, MIN, REAL, SQRT +* .. Model parameters .. + REAL BIGNUM, SAFMAX, SAFMIN, SMLNUM, ULP + PARAMETER (BIGNUM=0.1014120480E+32, + & SAFMAX=0.8507059173E+38, + & SAFMIN=0.1175494351E-37, + & SMLNUM=0.9860761315E-31, + & ULP=0.1192092896E-06) +* .. Local Scalars .. + COMPLEX ROGUE + REAL SNRM, TRAT, V0, V1, WORKSSQ, Y1, Y2, + & YMAX, YMIN, YNRM, ZNRM + INTEGER I, IV, IW, IX, KS + LOGICAL FIRST +* .. Local Arrays .. + COMPLEX X(NMAX), Z(NMAX) + REAL VALUES(NV), WORK(NMAX) +* .. Executable Statements .. + VALUES(1) = ZERO + VALUES(2) = TWO*SAFMIN + VALUES(3) = SMLNUM + VALUES(4) = ULP + VALUES(5) = ONE + VALUES(6) = ONE / ULP + VALUES(7) = BIGNUM + VALUES(8) = SAFMAX + VALUES(9) = SXVALS(V0,2) + VALUES(10) = SXVALS(V0,3) + ROGUE = CMPLX(1234.5678E+0,-1234.5678E+0) + FIRST = .TRUE. +* +* Check that the arrays are large enough +* + IF (N*ABS(INCX).GT.NMAX) THEN + WRITE (NOUT,99) "SCNRM2", NMAX, INCX, N, N*ABS(INCX) + RETURN + END IF +* +* Zero-sized inputs are tested in STEST1. + IF (N.LE.0) THEN + RETURN + END IF +* +* Generate 2*(N-1) values in (-1,1). +* + KS = 2*(N-1) + DO I = 1, KS + CALL RANDOM_NUMBER(WORK(I)) + WORK(I) = ONE - TWO*WORK(I) + END DO +* +* Compute the sum of squares of the random values +* by an unscaled algorithm. +* + WORKSSQ = ZERO + DO I = 1, KS + WORKSSQ = WORKSSQ + WORK(I)*WORK(I) + END DO +* +* Construct the test vector with one known value +* and the rest from the random work array multiplied +* by a scaling factor. +* + DO IV = 1, NV + V0 = VALUES(IV) + IF (ABS(V0).GT.ONE) THEN + V0 = V0*HALF*HALF + END IF + Z(1) = CMPLX(V0,-THREE*V0) + DO IW = 1, NV + V1 = VALUES(IW) + IF (ABS(V1).GT.ONE) THEN + V1 = (V1*HALF) / SQRT(REAL(KS+1)) + END IF + DO I = 1, N-1 + Z(I+1) = CMPLX(V1*WORK(2*I-1),V1*WORK(2*I)) + END DO +* +* Compute the expected value of the 2-norm +* + Y1 = ABS(V0) * SQRT(10.0E0) + IF (N.GT.1) THEN + Y2 = ABS(V1)*SQRT(WORKSSQ) + ELSE + Y2 = ZERO + END IF + YMIN = MIN(Y1, Y2) + YMAX = MAX(Y1, Y2) +* +* Expected value is NaN if either is NaN. The test +* for YMIN == YMAX avoids further computation if both +* are infinity. +* + IF ((Y1.NE.Y1).OR.(Y2.NE.Y2)) THEN +* add to propagate NaN + YNRM = Y1 + Y2 + ELSE IF (YMIN == YMAX) THEN + YNRM = SQRT(TWO)*YMAX + ELSE IF (YMAX == ZERO) THEN + YNRM = ZERO + ELSE + YNRM = YMAX*SQRT(ONE + (YMIN / YMAX)**2) + END IF +* +* Fill the input array to SCNRM2 with steps of incx +* + DO I = 1, N + X(I) = ROGUE + END DO + IX = 1 + IF (INCX.LT.0) IX = 1 - (N-1)*INCX + DO I = 1, N + X(IX) = Z(I) + IX = IX + INCX + END DO +* +* Call SCNRM2 to compute the 2-norm +* + SNRM = SCNRM2(N,X,INCX) +* +* Compare SNRM and ZNRM. Roundoff error grows like O(n) +* in this implementation so we scale the test ratio accordingly. +* + IF (INCX.EQ.0) THEN + Y1 = ABS(REAL(X(1))) + Y2 = ABS(AIMAG(X(1))) + YMIN = MIN(Y1, Y2) + YMAX = MAX(Y1, Y2) + IF ((Y1.NE.Y1).OR.(Y2.NE.Y2)) THEN +* add to propagate NaN + ZNRM = Y1 + Y2 + ELSE IF (YMIN == YMAX) THEN + ZNRM = SQRT(TWO)*YMAX + ELSE IF (YMAX == ZERO) THEN + ZNRM = ZERO + ELSE + ZNRM = YMAX * SQRT(ONE + (YMIN / YMAX)**2) + END IF + ZNRM = SQRT(REAL(n)) * ZNRM + ELSE + ZNRM = YNRM + END IF +* +* The tests for NaN rely on the compiler not being overly +* aggressive and removing the statements altogether. + IF ((SNRM.NE.SNRM).OR.(ZNRM.NE.ZNRM)) THEN + IF ((SNRM.NE.SNRM).NEQV.(ZNRM.NE.ZNRM)) THEN + TRAT = ONE / ULP + ELSE + TRAT = ZERO + END IF + ELSE IF (ZNRM == ZERO) THEN + TRAT = SNRM / ULP + ELSE + TRAT = (ABS(SNRM-ZNRM) / ZNRM) / (TWO*REAL(N)*ULP) + END IF + IF ((TRAT.NE.TRAT).OR.(TRAT.GE.THRESH)) THEN + IF (FIRST) THEN + FIRST = .FALSE. + WRITE(NOUT,99999) + END IF + WRITE (NOUT,98) "SCNRM2", N, INCX, IV, IW, TRAT + END IF + END DO + END DO +99999 FORMAT (' FAIL') + 99 FORMAT ( ' Not enough space to test ', A6, ': NMAX = ',I6, + + ', INCX = ',I6,/,' N = ',I6,', must be at least ',I6 ) + 98 FORMAT( 1X, A6, ': N=', I6,', INCX=', I4, ', IV=', I2, ', IW=', + + I2, ', test=', E15.8 ) + RETURN + CONTAINS + REAL FUNCTION SXVALS(XX,K) +* .. Scalar Arguments .. + REAL XX + INTEGER K +* .. Local Scalars .. + REAL X, Y, YY, Z +* .. Intrinsic Functions .. + INTRINSIC HUGE +* .. Executable Statements .. + Y = HUGE(XX) + Z = YY + IF (K.EQ.1) THEN + X = -Z + ELSE IF (K.EQ.2) THEN + X = Z + ELSE IF (K.EQ.3) THEN + X = Z / Z + END IF + SXVALS = X + RETURN + END + END diff --git a/BLAS/TESTING/dblat1.f b/BLAS/TESTING/dblat1.f index 9a27a249b2..063ffac3d2 100644 --- a/BLAS/TESTING/dblat1.f +++ b/BLAS/TESTING/dblat1.f @@ -247,8 +247,9 @@ SUBROUTINE CHECK0(SFAC) END SUBROUTINE CHECK1(SFAC) * .. Parameters .. + DOUBLE PRECISION THRESH INTEGER NOUT - PARAMETER (NOUT=6) + PARAMETER (NOUT=6, THRESH=10.0D0) * .. Scalar Arguments .. DOUBLE PRECISION SFAC * .. Scalars in Common .. @@ -266,7 +267,7 @@ SUBROUTINE CHECK1(SFAC) INTEGER IDAMAX EXTERNAL DASUM, DNRM2, IDAMAX * .. External Subroutines .. - EXTERNAL ITEST1, DSCAL, STEST, STEST1 + EXTERNAL ITEST1, DB1NRM2, DSCAL, STEST, STEST1 * .. Intrinsic Functions .. INTRINSIC MAX * .. Common blocks .. @@ -319,6 +320,10 @@ SUBROUTINE CHECK1(SFAC) * IF (ICASE.EQ.7) THEN * .. DNRM2 .. +* Test scaling when some entries are tiny or huge + CALL DB1NRM2(N,(INCX-2)*2,THRESH) + CALL DB1NRM2(N,INCX,THRESH) +* Test with hardcoded mid range entries STEMP(1) = DTRUE1(NP1) CALL STEST1(DNRM2(N,SX,INCX),STEMP(1),STEMP,SFAC) ELSE IF (ICASE.EQ.8) THEN @@ -1129,3 +1134,213 @@ SUBROUTINE ITEST1(ICOMP,ITRUE) * End of ITEST1 * END + SUBROUTINE DB1NRM2(N,INCX,THRESH) +* Compare NRM2 with a reference computation using combinations +* of the following values: +* +* 0, very small, small, ulp, 1, 1/ulp, big, very big, infinity, NaN +* +* one of these values is used to initialize x(1) and x(2:N) is +* filled with random values from [-1,1] scaled by another of +* these values. +* +* This routine is adapted from the test suite provided by +* Anderson E. (2017) +* Algorithm 978: Safe Scaling in the Level 1 BLAS +* ACM Trans Math Softw 44:1--28 +* https://doi.org/10.1145/3061665 +* +* .. Scalar Arguments .. + INTEGER INCX, N + DOUBLE PRECISION THRESH +* +* ===================================================================== +* .. Parameters .. + INTEGER NMAX, NOUT, NV + PARAMETER (NMAX=20, NOUT=6, NV=10) + DOUBLE PRECISION HALF, ONE, TWO, ZERO + PARAMETER (HALF=0.5D+0, ONE=1.0D+0, TWO= 2.0D+0, + & ZERO=0.0D+0) +* .. External Functions .. + DOUBLE PRECISION DNRM2 + EXTERNAL DNRM2 +* .. Intrinsic Functions .. + INTRINSIC ABS, DBLE, MAX, MIN, SQRT +* .. Model parameters .. + DOUBLE PRECISION BIGNUM, SAFMAX, SAFMIN, SMLNUM, ULP + PARAMETER (BIGNUM=0.99792015476735990583D+292, + & SAFMAX=0.44942328371557897693D+308, + & SAFMIN=0.22250738585072013831D-307, + & SMLNUM=0.10020841800044863890D-291, + & ULP=0.22204460492503130808D-015) +* .. Local Scalars .. + DOUBLE PRECISION ROGUE, SNRM, TRAT, V0, V1, WORKSSQ, Y1, Y2, + & YMAX, YMIN, YNRM, ZNRM + INTEGER I, IV, IW, IX + LOGICAL FIRST +* .. Local Arrays .. + DOUBLE PRECISION VALUES(NV), WORK(NMAX), X(NMAX), Z(NMAX) +* .. Executable Statements .. + VALUES(1) = ZERO + VALUES(2) = TWO*SAFMIN + VALUES(3) = SMLNUM + VALUES(4) = ULP + VALUES(5) = ONE + VALUES(6) = ONE / ULP + VALUES(7) = BIGNUM + VALUES(8) = SAFMAX + VALUES(9) = DXVALS(V0,2) + VALUES(10) = DXVALS(V0,3) + ROGUE = -1234.5678D+0 + FIRST = .TRUE. +* +* Check that the arrays are large enough +* + IF (N*ABS(INCX).GT.NMAX) THEN + WRITE (NOUT,99) "DNRM2", NMAX, INCX, N, N*ABS(INCX) + RETURN + END IF +* +* Zero-sized inputs are tested in STEST1. + IF (N.LE.0) THEN + RETURN + END IF +* +* Generate (N-1) values in (-1,1). +* + DO I = 2, N + CALL RANDOM_NUMBER(WORK(I)) + WORK(I) = ONE - TWO*WORK(I) + END DO +* +* Compute the sum of squares of the random values +* by an unscaled algorithm. +* + WORKSSQ = ZERO + DO I = 2, N + WORKSSQ = WORKSSQ + WORK(I)*WORK(I) + END DO +* +* Construct the test vector with one known value +* and the rest from the random work array multiplied +* by a scaling factor. +* + DO IV = 1, NV + V0 = VALUES(IV) + IF (ABS(V0).GT.ONE) THEN + V0 = V0*HALF + END IF + Z(1) = V0 + DO IW = 1, NV + V1 = VALUES(IW) + IF (ABS(V1).GT.ONE) THEN + V1 = (V1*HALF) / SQRT(DBLE(N)) + END IF + DO I = 2, N + Z(I) = V1*WORK(I) + END DO +* +* Compute the expected value of the 2-norm +* + Y1 = ABS(V0) + IF (N.GT.1) THEN + Y2 = ABS(V1)*SQRT(WORKSSQ) + ELSE + Y2 = ZERO + END IF + YMIN = MIN(Y1, Y2) + YMAX = MAX(Y1, Y2) +* +* Expected value is NaN if either is NaN. The test +* for YMIN == YMAX avoids further computation if both +* are infinity. +* + IF ((Y1.NE.Y1).OR.(Y2.NE.Y2)) THEN +* Add to propagate NaN + YNRM = Y1 + Y2 + ELSE IF (YMAX == ZERO) THEN + YNRM = ZERO + ELSE IF (YMIN == YMAX) THEN + YNRM = SQRT(TWO)*YMAX + ELSE + YNRM = YMAX*SQRT(ONE + (YMIN / YMAX)**2) + END IF +* +* Fill the input array to DNRM2 with steps of incx +* + DO I = 1, N + X(I) = ROGUE + END DO + IX = 1 + IF (INCX.LT.0) IX = 1 - (N-1)*INCX + DO I = 1, N + X(IX) = Z(I) + IX = IX + INCX + END DO +* +* Call DNRM2 to compute the 2-norm +* + SNRM = DNRM2(N,X,INCX) +* +* Compare SNRM and ZNRM. Roundoff error grows like O(n) +* in this implementation so we scale the test ratio accordingly. +* + IF (INCX.EQ.0) THEN + ZNRM = SQRT(DBLE(N))*ABS(X(1)) + ELSE + ZNRM = YNRM + END IF +* +* The tests for NaN rely on the compiler not being overly +* aggressive and removing the statements altogether. + IF ((SNRM.NE.SNRM).OR.(ZNRM.NE.ZNRM)) THEN + IF ((SNRM.NE.SNRM).NEQV.(ZNRM.NE.ZNRM)) THEN + TRAT = ONE / ULP + ELSE + TRAT = ZERO + END IF + ELSE IF (SNRM == ZNRM) THEN + TRAT = ZERO + ELSE IF (ZNRM == ZERO) THEN + TRAT = SNRM / ULP + ELSE + TRAT = (ABS(SNRM-ZNRM) / ZNRM) / (DBLE(N)*ULP) + END IF + IF ((TRAT.NE.TRAT).OR.(TRAT.GE.THRESH)) THEN + IF (FIRST) THEN + FIRST = .FALSE. + WRITE(NOUT,99999) + END IF + WRITE (NOUT,98) "DNRM2", N, INCX, IV, IW, TRAT + END IF + END DO + END DO +99999 FORMAT (' FAIL') + 99 FORMAT ( ' Not enough space to test ', A6, ': NMAX = ',I6, + + ', INCX = ',I6,/,' N = ',I6,', must be at least ',I6 ) + 98 FORMAT( 1X, A6, ': N=', I6,', INCX=', I4, ', IV=', I2, ', IW=', + + I2, ', test=', E15.8 ) + RETURN + CONTAINS + DOUBLE PRECISION FUNCTION DXVALS(XX,K) +* .. Scalar Arguments .. + DOUBLE PRECISION XX + INTEGER K +* .. Local Scalars .. + DOUBLE PRECISION X, Y, YY, Z +* .. Intrinsic Functions .. + INTRINSIC HUGE +* .. Executable Statements .. + Y = HUGE(XX) + Z = YY + IF (K.EQ.1) THEN + X = -Z + ELSE IF (K.EQ.2) THEN + X = Z + ELSE IF (K.EQ.3) THEN + X = Z / Z + END IF + DXVALS = X + RETURN + END + END diff --git a/BLAS/TESTING/sblat1.f b/BLAS/TESTING/sblat1.f index d81496b6a3..4dc537e2f0 100644 --- a/BLAS/TESTING/sblat1.f +++ b/BLAS/TESTING/sblat1.f @@ -248,7 +248,8 @@ SUBROUTINE CHECK0(SFAC) SUBROUTINE CHECK1(SFAC) * .. Parameters .. INTEGER NOUT - PARAMETER (NOUT=6) + REAL THRESH + PARAMETER (NOUT=6, THRESH=10.0E0) * .. Scalar Arguments .. REAL SFAC * .. Scalars in Common .. @@ -266,7 +267,7 @@ SUBROUTINE CHECK1(SFAC) INTEGER ISAMAX EXTERNAL SASUM, SNRM2, ISAMAX * .. External Subroutines .. - EXTERNAL ITEST1, SSCAL, STEST, STEST1 + EXTERNAL ITEST1, SB1NRM2, SSCAL, STEST, STEST1 * .. Intrinsic Functions .. INTRINSIC MAX * .. Common blocks .. @@ -319,6 +320,10 @@ SUBROUTINE CHECK1(SFAC) * IF (ICASE.EQ.7) THEN * .. SNRM2 .. +* Test scaling when some entries are tiny or huge + CALL SB1NRM2(N,(INCX-2)*2,THRESH) + CALL SB1NRM2(N,INCX,THRESH) +* Test with hardcoded mid range entries STEMP(1) = DTRUE1(NP1) CALL STEST1(SNRM2(N,SX,INCX),STEMP(1),STEMP,SFAC) ELSE IF (ICASE.EQ.8) THEN @@ -1080,3 +1085,214 @@ SUBROUTINE ITEST1(ICOMP,ITRUE) * End of ITEST1 * END + SUBROUTINE SB1NRM2(N,INCX,THRESH) +* Compare NRM2 with a reference computation using combinations +* of the following values: +* +* 0, very small, small, ulp, 1, 1/ulp, big, very big, infinity, NaN +* +* one of these values is used to initialize x(1) and x(2:N) is +* filled with random values from [-1,1] scaled by another of +* these values. +* +* This routine is adapted from the test suite provided by +* Anderson E. (2017) +* Algorithm 978: Safe Scaling in the Level 1 BLAS +* ACM Trans Math Softw 44:1--28 +* https://doi.org/10.1145/3061665 +* + IMPLICIT NONE +* .. Scalar Arguments .. + INTEGER INCX, N + REAL THRESH +* +* ===================================================================== +* .. Parameters .. + INTEGER NMAX, NOUT, NV + PARAMETER (NMAX=20, NOUT=6, NV=10) + REAL HALF, ONE, TWO, ZERO + PARAMETER (HALF=0.5E+0, ONE=1.0E+0, TWO= 2.0E+0, + & ZERO=0.0E+0) +* .. External Functions .. + REAL SNRM2 + EXTERNAL SNRM2 +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX, MIN, REAL, SQRT +* .. Model parameters .. + REAL BIGNUM, SAFMAX, SAFMIN, SMLNUM, ULP + PARAMETER (BIGNUM=0.1014120480E+32, + & SAFMAX=0.8507059173E+38, + & SAFMIN=0.1175494351E-37, + & SMLNUM=0.9860761315E-31, + & ULP=0.1192092896E-06) +* .. Local Scalars .. + REAL ROGUE, SNRM, TRAT, V0, V1, WORKSSQ, Y1, Y2, + & YMAX, YMIN, YNRM, ZNRM + INTEGER I, IV, IW, IX + LOGICAL FIRST +* .. Local Arrays .. + REAL VALUES(NV), WORK(NMAX), X(NMAX), Z(NMAX) +* .. Executable Statements .. + VALUES(1) = ZERO + VALUES(2) = TWO*SAFMIN + VALUES(3) = SMLNUM + VALUES(4) = ULP + VALUES(5) = ONE + VALUES(6) = ONE / ULP + VALUES(7) = BIGNUM + VALUES(8) = SAFMAX + VALUES(9) = SXVALS(V0,2) + VALUES(10) = SXVALS(V0,3) + ROGUE = -1234.5678E+0 + FIRST = .TRUE. +* +* Check that the arrays are large enough +* + IF (N*ABS(INCX).GT.NMAX) THEN + WRITE (NOUT,99) "SNRM2", NMAX, INCX, N, N*ABS(INCX) + RETURN + END IF +* +* Zero-sized inputs are tested in STEST1. + IF (N.LE.0) THEN + RETURN + END IF +* +* Generate (N-1) values in (-1,1). +* + DO I = 2, N + CALL RANDOM_NUMBER(WORK(I)) + WORK(I) = ONE - TWO*WORK(I) + END DO +* +* Compute the sum of squares of the random values +* by an unscaled algorithm. +* + WORKSSQ = ZERO + DO I = 2, N + WORKSSQ = WORKSSQ + WORK(I)*WORK(I) + END DO +* +* Construct the test vector with one known value +* and the rest from the random work array multiplied +* by a scaling factor. +* + DO IV = 1, NV + V0 = VALUES(IV) + IF (ABS(V0).GT.ONE) THEN + V0 = V0*HALF + END IF + Z(1) = V0 + DO IW = 1, NV + V1 = VALUES(IW) + IF (ABS(V1).GT.ONE) THEN + V1 = (V1*HALF) / SQRT(REAL(N)) + END IF + DO I = 2, N + Z(I) = V1*WORK(I) + END DO +* +* Compute the expected value of the 2-norm +* + Y1 = ABS(V0) + IF (N.GT.1) THEN + Y2 = ABS(V1)*SQRT(WORKSSQ) + ELSE + Y2 = ZERO + END IF + YMIN = MIN(Y1, Y2) + YMAX = MAX(Y1, Y2) +* +* Expected value is NaN if either is NaN. The test +* for YMIN == YMAX avoids further computation if both +* are infinity. +* + IF ((Y1.NE.Y1).OR.(Y2.NE.Y2)) THEN +* add to propagate NaN + YNRM = Y1 + Y2 + ELSE IF (YMIN == YMAX) THEN + YNRM = SQRT(TWO)*YMAX + ELSE IF (YMAX == ZERO) THEN + YNRM = ZERO + ELSE + YNRM = YMAX*SQRT(ONE + (YMIN / YMAX)**2) + END IF +* +* Fill the input array to SNRM2 with steps of incx +* + DO I = 1, N + X(I) = ROGUE + END DO + IX = 1 + IF (INCX.LT.0) IX = 1 - (N-1)*INCX + DO I = 1, N + X(IX) = Z(I) + IX = IX + INCX + END DO +* +* Call SNRM2 to compute the 2-norm +* + SNRM = SNRM2(N,X,INCX) +* +* Compare SNRM and ZNRM. Roundoff error grows like O(n) +* in this implementation so we scale the test ratio accordingly. +* + IF (INCX.EQ.0) THEN + ZNRM = SQRT(REAL(N))*ABS(X(1)) + ELSE + ZNRM = YNRM + END IF +* +* The tests for NaN rely on the compiler not being overly +* aggressive and removing the statements altogether. + IF ((SNRM.NE.SNRM).OR.(ZNRM.NE.ZNRM)) THEN + IF ((SNRM.NE.SNRM).NEQV.(ZNRM.NE.ZNRM)) THEN + TRAT = ONE / ULP + ELSE + TRAT = ZERO + END IF + ELSE IF (SNRM == ZNRM) THEN + TRAT = ZERO + ELSE IF (ZNRM == ZERO) THEN + TRAT = SNRM / ULP + ELSE + TRAT = (ABS(SNRM-ZNRM) / ZNRM) / (REAL(N)*ULP) + END IF + IF ((TRAT.NE.TRAT).OR.(TRAT.GE.THRESH)) THEN + IF (FIRST) THEN + FIRST = .FALSE. + WRITE(NOUT,99999) + END IF + WRITE (NOUT,98) "SNRM2", N, INCX, IV, IW, TRAT + END IF + END DO + END DO +99999 FORMAT (' FAIL') + 99 FORMAT ( ' Not enough space to test ', A6, ': NMAX = ',I6, + + ', INCX = ',I6,/,' N = ',I6,', must be at least ',I6 ) + 98 FORMAT( 1X, A6, ': N=', I6,', INCX=', I4, ', IV=', I2, ', IW=', + + I2, ', test=', E15.8 ) + RETURN + CONTAINS + REAL FUNCTION SXVALS(XX,K) +* .. Scalar Arguments .. + REAL XX + INTEGER K +* .. Local Scalars .. + REAL X, Y, YY, Z +* .. Intrinsic Functions .. + INTRINSIC HUGE +* .. Executable Statements .. + Y = HUGE(XX) + Z = YY + IF (K.EQ.1) THEN + X = -Z + ELSE IF (K.EQ.2) THEN + X = Z + ELSE IF (K.EQ.3) THEN + X = Z / Z + END IF + SXVALS = X + RETURN + END + END diff --git a/BLAS/TESTING/zblat1.f b/BLAS/TESTING/zblat1.f index 118212d098..ef6deff921 100644 --- a/BLAS/TESTING/zblat1.f +++ b/BLAS/TESTING/zblat1.f @@ -121,7 +121,8 @@ SUBROUTINE HEADER SUBROUTINE CHECK1(SFAC) * .. Parameters .. INTEGER NOUT - PARAMETER (NOUT=6) + DOUBLE PRECISION THRESH + PARAMETER (NOUT=6, THRESH=10.0D0) * .. Scalar Arguments .. DOUBLE PRECISION SFAC * .. Scalars in Common .. @@ -141,7 +142,7 @@ SUBROUTINE CHECK1(SFAC) INTEGER IZAMAX EXTERNAL DZASUM, DZNRM2, IZAMAX * .. External Subroutines .. - EXTERNAL ZSCAL, ZDSCAL, CTEST, ITEST1, STEST1 + EXTERNAL ZB1NRM2, ZSCAL, ZDSCAL, CTEST, ITEST1, STEST1 * .. Intrinsic Functions .. INTRINSIC MAX * .. Common blocks .. @@ -256,6 +257,10 @@ SUBROUTINE CHECK1(SFAC) 20 CONTINUE IF (ICASE.EQ.6) THEN * .. DZNRM2 .. +* Test scaling when some entries are tiny or huge + CALL ZB1NRM2(N,(INCX-2)*2,THRESH) + CALL ZB1NRM2(N,INCX,THRESH) +* Test with hardcoded mid range entries CALL STEST1(DZNRM2(N,CX,INCX),STRUE2(NP1),STRUE2(NP1), + SFAC) ELSE IF (ICASE.EQ.7) THEN @@ -782,3 +787,228 @@ SUBROUTINE ITEST1(ICOMP,ITRUE) * End of ITEST1 * END + SUBROUTINE ZB1NRM2(N,INCX,THRESH) +* Compare NRM2 with a reference computation using combinations +* of the following values: +* +* 0, very small, small, ulp, 1, 1/ulp, big, very big, infinity, NaN +* +* one of these values is used to initialize x(1) and x(2:N) is +* filled with random values from [-1,1] scaled by another of +* these values. +* +* This routine is adapted from the test suite provided by +* Anderson E. (2017) +* Algorithm 978: Safe Scaling in the Level 1 BLAS +* ACM Trans Math Softw 44:1--28 +* https://doi.org/10.1145/3061665 +* +* .. Scalar Arguments .. + INTEGER INCX, N + DOUBLE PRECISION THRESH +* +* ===================================================================== +* .. Parameters .. + INTEGER NMAX, NOUT, NV + PARAMETER (NMAX=20, NOUT=6, NV=10) + DOUBLE PRECISION HALF, ONE, THREE, TWO, ZERO + PARAMETER (HALF=0.5D+0, ONE=1.0D+0, TWO= 2.0D+0, + & THREE=3.0D+0, ZERO=0.0D+0) +* .. External Functions .. + DOUBLE PRECISION DZNRM2 + EXTERNAL DZNRM2 +* .. Intrinsic Functions .. + INTRINSIC AIMAG, ABS, DCMPLX, DBLE, MAX, MIN, SQRT +* .. Model parameters .. + DOUBLE PRECISION BIGNUM, SAFMAX, SAFMIN, SMLNUM, ULP + PARAMETER (BIGNUM=0.99792015476735990583D+292, + & SAFMAX=0.44942328371557897693D+308, + & SAFMIN=0.22250738585072013831D-307, + & SMLNUM=0.10020841800044863890D-291, + & ULP=0.22204460492503130808D-015) +* .. Local Scalars .. + COMPLEX*16 ROGUE + DOUBLE PRECISION SNRM, TRAT, V0, V1, WORKSSQ, Y1, Y2, + & YMAX, YMIN, YNRM, ZNRM + INTEGER I, IV, IW, IX, KS + LOGICAL FIRST +* .. Local Arrays .. + COMPLEX*16 X(NMAX), Z(NMAX) + DOUBLE PRECISION VALUES(NV), WORK(NMAX) +* .. Executable Statements .. + VALUES(1) = ZERO + VALUES(2) = TWO*SAFMIN + VALUES(3) = SMLNUM + VALUES(4) = ULP + VALUES(5) = ONE + VALUES(6) = ONE / ULP + VALUES(7) = BIGNUM + VALUES(8) = SAFMAX + VALUES(9) = DXVALS(V0,2) + VALUES(10) = DXVALS(V0,3) + ROGUE = DCMPLX(1234.5678D+0,-1234.5678D+0) + FIRST = .TRUE. +* +* Check that the arrays are large enough +* + IF (N*ABS(INCX).GT.NMAX) THEN + WRITE (NOUT,99) "DZNRM2", NMAX, INCX, N, N*ABS(INCX) + RETURN + END IF +* +* Zero-sized inputs are tested in STEST1. + IF (N.LE.0) THEN + RETURN + END IF +* +* Generate 2*(N-1) values in (-1,1). +* + KS = 2*(N-1) + DO I = 1, KS + CALL RANDOM_NUMBER(WORK(I)) + WORK(I) = ONE - TWO*WORK(I) + END DO +* +* Compute the sum of squares of the random values +* by an unscaled algorithm. +* + WORKSSQ = ZERO + DO I = 1, KS + WORKSSQ = WORKSSQ + WORK(I)*WORK(I) + END DO +* +* Construct the test vector with one known value +* and the rest from the random work array multiplied +* by a scaling factor. +* + DO IV = 1, NV + V0 = VALUES(IV) + IF (ABS(V0).GT.ONE) THEN + V0 = V0*HALF*HALF + END IF + Z(1) = DCMPLX(V0,-THREE*V0) + DO IW = 1, NV + V1 = VALUES(IW) + IF (ABS(V1).GT.ONE) THEN + V1 = (V1*HALF) / SQRT(DBLE(KS+1)) + END IF + DO I = 1, N-1 + Z(I+1) = DCMPLX(V1*WORK(2*I-1),V1*WORK(2*I)) + END DO +* +* Compute the expected value of the 2-norm +* + Y1 = ABS(V0) * SQRT(10.0D0) + IF (N.GT.1) THEN + Y2 = ABS(V1)*SQRT(WORKSSQ) + ELSE + Y2 = ZERO + END IF + YMIN = MIN(Y1, Y2) + YMAX = MAX(Y1, Y2) +* +* Expected value is NaN if either is NaN. The test +* for YMIN == YMAX avoids further computation if both +* are infinity. +* + IF ((Y1.NE.Y1).OR.(Y2.NE.Y2)) THEN +* add to propagate NaN + YNRM = Y1 + Y2 + ELSE IF (YMIN == YMAX) THEN + YNRM = SQRT(TWO)*YMAX + ELSE IF (YMAX == ZERO) THEN + YNRM = ZERO + ELSE + YNRM = YMAX*SQRT(ONE + (YMIN / YMAX)**2) + END IF +* +* Fill the input array to DZNRM2 with steps of incx +* + DO I = 1, N + X(I) = ROGUE + END DO + IX = 1 + IF (INCX.LT.0) IX = 1 - (N-1)*INCX + DO I = 1, N + X(IX) = Z(I) + IX = IX + INCX + END DO +* +* Call DZNRM2 to compute the 2-norm +* + SNRM = DZNRM2(N,X,INCX) +* +* Compare SNRM and ZNRM. Roundoff error grows like O(n) +* in this implementation so we scale the test ratio accordingly. +* + IF (INCX.EQ.0) THEN + Y1 = ABS(DBLE(X(1))) + Y2 = ABS(AIMAG(X(1))) + YMIN = MIN(Y1, Y2) + YMAX = MAX(Y1, Y2) + IF ((Y1.NE.Y1).OR.(Y2.NE.Y2)) THEN +* add to propagate NaN + ZNRM = Y1 + Y2 + ELSE IF (YMIN == YMAX) THEN + ZNRM = SQRT(TWO)*YMAX + ELSE IF (YMAX == ZERO) THEN + ZNRM = ZERO + ELSE + ZNRM = YMAX * SQRT(ONE + (YMIN / YMAX)**2) + END IF + ZNRM = SQRT(DBLE(n)) * ZNRM + ELSE + ZNRM = YNRM + END IF +* +* The tests for NaN rely on the compiler not being overly +* aggressive and removing the statements altogether. + IF ((SNRM.NE.SNRM).OR.(ZNRM.NE.ZNRM)) THEN + IF ((SNRM.NE.SNRM).NEQV.(ZNRM.NE.ZNRM)) THEN + TRAT = ONE / ULP + ELSE + TRAT = ZERO + END IF + ELSE IF (ZNRM == ZERO) THEN + TRAT = SNRM / ULP + ELSE + TRAT = (ABS(SNRM-ZNRM) / ZNRM) / (TWO*DBLE(N)*ULP) + END IF + IF ((TRAT.NE.TRAT).OR.(TRAT.GE.THRESH)) THEN + IF (FIRST) THEN + FIRST = .FALSE. + WRITE(NOUT,99999) + END IF + WRITE (NOUT,98) "DZNRM2", N, INCX, IV, IW, TRAT + END IF + END DO + END DO +99999 FORMAT (' FAIL') + 99 FORMAT ( ' Not enough space to test ', A6, ': NMAX = ',I6, + + ', INCX = ',I6,/,' N = ',I6,', must be at least ',I6 ) + 98 FORMAT( 1X, A6, ': N=', I6,', INCX=', I4, ', IV=', I2, ', IW=', + + I2, ', test=', E15.8 ) + RETURN + CONTAINS + DOUBLE PRECISION FUNCTION DXVALS(XX,K) +* .. Scalar Arguments .. + DOUBLE PRECISION XX + INTEGER K +* .. Local Scalars .. + DOUBLE PRECISION X, Y, YY, Z +* .. Intrinsic Functions .. + INTRINSIC HUGE +* .. Executable Statements .. + Y = HUGE(XX) + Z = YY + IF (K.EQ.1) THEN + X = -Z + ELSE IF (K.EQ.2) THEN + X = Z + ELSE IF (K.EQ.3) THEN + X = Z / Z + END IF + DXVALS = X + RETURN + END + END