diff --git a/test/cblat1.f b/test/cblat1.f index ecf2a44cb1..83b02f0cac 100644 --- a/test/cblat1.f +++ b/test/cblat1.f @@ -30,17 +30,14 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \date April 2012 -* *> \ingroup complex_blas_testing * * ===================================================================== PROGRAM CBLAT1 * -* -- Reference BLAS test routine (version 3.7.0) -- +* -- Reference BLAS test routine -- * -- Reference BLAS is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* April 2012 * * ===================================================================== * @@ -86,6 +83,9 @@ PROGRAM CBLAT1 * 99999 FORMAT (' Complex BLAS Test Program Results',/1X) 99998 FORMAT (' ----- PASS -----') +* +* End of CBLAT1 +* END SUBROUTINE HEADER * .. Parameters .. @@ -114,11 +114,15 @@ SUBROUTINE HEADER RETURN * 99999 FORMAT (/' Test of subprogram number',I3,12X,A6) +* +* End of HEADER +* END SUBROUTINE CHECK1(SFAC) * .. Parameters .. INTEGER NOUT - PARAMETER (NOUT=6) + REAL THRESH + PARAMETER (NOUT=6, THRESH=10.0E0) * .. Scalar Arguments .. REAL SFAC * .. Scalars in Common .. @@ -127,18 +131,18 @@ SUBROUTINE CHECK1(SFAC) * .. Local Scalars .. COMPLEX CA REAL SA - INTEGER I, J, LEN, NP1 + INTEGER I, IX, J, LEN, NP1 * .. Local Arrays .. - COMPLEX CTRUE5(8,5,2), CTRUE6(8,5,2), CV(8,5,2), CX(8), - + MWPCS(5), MWPCT(5) + COMPLEX CTRUE5(8,5,2), CTRUE6(8,5,2), CV(8,5,2), CVR(8), + + CX(8), CXR(15), MWPCS(5), MWPCT(5) REAL STRUE2(5), STRUE4(5) - INTEGER ITRUE3(5) + INTEGER ITRUE3(5), ITRUEC(5) * .. External Functions .. REAL SCASUM, SCNRM2 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 .. @@ -173,6 +177,9 @@ SUBROUTINE CHECK1(SFAC) + (7.0E0,2.0E0), (0.3E0,0.1E0), (5.0E0,8.0E0), + (0.5E0,0.0E0), (6.0E0,9.0E0), (0.0E0,0.5E0), + (8.0E0,3.0E0), (0.0E0,0.2E0), (9.0E0,4.0E0)/ + DATA CVR/(8.0E0,8.0E0), (-7.0E0,-7.0E0), + + (9.0E0,9.0E0), (5.0E0,5.0E0), (9.0E0,9.0E0), + + (8.0E0,8.0E0), (7.0E0,7.0E0), (7.0E0,7.0E0)/ DATA STRUE2/0.0E0, 0.5E0, 0.6E0, 0.7E0, 0.8E0/ DATA STRUE4/0.0E0, 0.7E0, 1.0E0, 1.3E0, 1.6E0/ DATA ((CTRUE5(I,J,1),I=1,8),J=1,5)/(0.1E0,0.1E0), @@ -238,6 +245,7 @@ SUBROUTINE CHECK1(SFAC) + (0.15E0,0.00E0), (6.0E0,9.0E0), (0.00E0,0.15E0), + (8.0E0,3.0E0), (0.00E0,0.06E0), (9.0E0,4.0E0)/ DATA ITRUE3/0, 1, 2, 2, 2/ + DATA ITRUEC/0, 1, 1, 1, 1/ * .. Executable Statements .. DO 60 INCX = 1, 2 DO 40 NP1 = 1, 5 @@ -249,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 @@ -268,12 +280,25 @@ SUBROUTINE CHECK1(SFAC) ELSE IF (ICASE.EQ.10) THEN * .. ICAMAX .. CALL ITEST1(ICAMAX(N,CX,INCX),ITRUE3(NP1)) + DO 160 I = 1, LEN + CX(I) = (42.0E0,43.0E0) + 160 CONTINUE + CALL ITEST1(ICAMAX(N,CX,INCX),ITRUEC(NP1)) ELSE WRITE (NOUT,*) ' Shouldn''t be here in CHECK1' STOP END IF * 40 CONTINUE + IF (ICASE.EQ.10) THEN + N = 8 + IX = 1 + DO 180 I = 1, N + CXR(IX) = CVR(I) + IX = IX + INCX + 180 CONTINUE + CALL ITEST1(ICAMAX(N,CXR,INCX),3) + END IF 60 CONTINUE * INCX = 1 @@ -315,6 +340,9 @@ SUBROUTINE CHECK1(SFAC) CALL CTEST(5,CX,MWPCT,MWPCS,SFAC) END IF RETURN +* +* End of CHECK1 +* END SUBROUTINE CHECK2(SFAC) * .. Parameters .. @@ -327,11 +355,13 @@ SUBROUTINE CHECK2(SFAC) LOGICAL PASS * .. Local Scalars .. COMPLEX CA - INTEGER I, J, KI, KN, KSIZE, LENX, LENY, MX, MY + INTEGER I, J, KI, KN, KSIZE, LENX, LENY, LINCX, LINCY, + + MX, MY * .. Local Arrays .. COMPLEX CDOT(1), CSIZE1(4), CSIZE2(7,2), CSIZE3(14), + CT10X(7,4,4), CT10Y(7,4,4), CT6(4,4), CT7(4,4), - + CT8(7,4,4), CX(7), CX1(7), CY(7), CY1(7) + + CT8(7,4,4), CTY0(1), CX(7), CX0(1), CX1(7), + + CY(7), CY0(1), CY1(7) INTEGER INCXS(4), INCYS(4), LENS(4,2), NS(4) * .. External Functions .. COMPLEX CDOTC, CDOTU @@ -546,6 +576,23 @@ SUBROUTINE CHECK2(SFAC) * .. CCOPY .. CALL CCOPY(N,CX,INCX,CY,INCY) CALL CTEST(LENY,CY,CT10Y(1,KN,KI),CSIZE3,1.0E0) + IF (KI.EQ.1) THEN + CX0(1) = (42.0E0,43.0E0) + CY0(1) = (44.0E0,45.0E0) + IF (N.EQ.0) THEN + CTY0(1) = CY0(1) + ELSE + CTY0(1) = CX0(1) + END IF + LINCX = INCX + INCX = 0 + LINCY = INCY + INCY = 0 + CALL CCOPY(N,CX0,INCX,CY0,INCY) + CALL CTEST(1,CY0,CTY0,CSIZE3,1.0E0) + INCX = LINCX + INCY = LINCY + END IF ELSE IF (ICASE.EQ.5) THEN * .. CSWAP .. CALL CSWAP(N,CX,INCX,CY,INCY) @@ -559,6 +606,9 @@ SUBROUTINE CHECK2(SFAC) 40 CONTINUE 60 CONTINUE RETURN +* +* End of CHECK2 +* END SUBROUTINE STEST(LEN,SCOMP,STRUE,SSIZE,SFAC) * ********************************* STEST ************************** @@ -615,6 +665,9 @@ SUBROUTINE STEST(LEN,SCOMP,STRUE,SSIZE,SFAC) + ' COMP(I) TRUE(I) DIFFERENCE', + ' SIZE(I)',/1X) 99997 FORMAT (1X,I4,I3,3I5,I3,2E36.8,2E12.4) +* +* End of STEST +* END SUBROUTINE STEST1(SCOMP1,STRUE1,SSIZE,SFAC) * ************************* STEST1 ***************************** @@ -640,6 +693,9 @@ SUBROUTINE STEST1(SCOMP1,STRUE1,SSIZE,SFAC) CALL STEST(1,SCOMP,STRUE,SSIZE,SFAC) * RETURN +* +* End of STEST1 +* END REAL FUNCTION SDIFF(SA,SB) * ********************************* SDIFF ************************** @@ -650,6 +706,9 @@ REAL FUNCTION SDIFF(SA,SB) * .. Executable Statements .. SDIFF = SA - SB RETURN +* +* End of SDIFF +* END SUBROUTINE CTEST(LEN,CCOMP,CTRUE,CSIZE,SFAC) * **************************** CTEST ***************************** @@ -681,6 +740,9 @@ SUBROUTINE CTEST(LEN,CCOMP,CTRUE,CSIZE,SFAC) * CALL STEST(2*LEN,SCOMP,STRUE,SSIZE,SFAC) RETURN +* +* End of CTEST +* END SUBROUTINE ITEST1(ICOMP,ITRUE) * ********************************* ITEST1 ************************* @@ -721,4 +783,232 @@ SUBROUTINE ITEST1(ICOMP,ITRUE) + ' COMP TRUE DIFFERENCE', + /1X) 99997 FORMAT (1X,I4,I3,3I5,2I36,I12) +* +* 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/test/dblat1.f b/test/dblat1.f index 28af121cd7..063ffac3d2 100644 --- a/test/dblat1.f +++ b/test/dblat1.f @@ -30,17 +30,14 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \date April 2012 -* *> \ingroup double_blas_testing * * ===================================================================== PROGRAM DBLAT1 * -* -- Reference BLAS test routine (version 3.8.0) -- +* -- Reference BLAS test routine -- * -- Reference BLAS is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* April 2012 * * ===================================================================== * @@ -91,6 +88,9 @@ PROGRAM DBLAT1 * 99999 FORMAT (' Real BLAS Test Program Results',/1X) 99998 FORMAT (' ----- PASS -----') +* +* End of DBLAT1 +* END SUBROUTINE HEADER * .. Parameters .. @@ -122,6 +122,9 @@ SUBROUTINE HEADER RETURN * 99999 FORMAT (/' Test of subprogram number',I3,12X,A6) +* +* End of HEADER +* END SUBROUTINE CHECK0(SFAC) * .. Parameters .. @@ -238,28 +241,33 @@ SUBROUTINE CHECK0(SFAC) END IF 20 CONTINUE 40 RETURN +* +* End of CHECK0 +* 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 .. INTEGER ICASE, INCX, INCY, N LOGICAL PASS * .. Local Scalars .. - INTEGER I, LEN, NP1 + INTEGER I, IX, LEN, NP1 * .. Local Arrays .. DOUBLE PRECISION DTRUE1(5), DTRUE3(5), DTRUE5(8,5,2), DV(8,5,2), - + SA(10), STEMP(1), STRUE(8), SX(8) - INTEGER ITRUE2(5) + + DVR(8), SA(10), STEMP(1), STRUE(8), SX(8), + + SXR(15) + INTEGER ITRUE2(5), ITRUEC(5) * .. External Functions .. DOUBLE PRECISION DASUM, DNRM2 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 .. @@ -280,6 +288,8 @@ SUBROUTINE CHECK1(SFAC) + 0.2D0, 3.0D0, -0.6D0, 5.0D0, 0.3D0, 2.0D0, + 2.0D0, 2.0D0, 0.1D0, 4.0D0, -0.3D0, 6.0D0, + -0.5D0, 7.0D0, -0.1D0, 3.0D0/ + DATA DVR/8.0D0, -7.0D0, 9.0D0, 5.0D0, 9.0D0, 8.0D0, + + 7.0D0, 7.0D0/ DATA DTRUE1/0.0D0, 0.3D0, 0.5D0, 0.7D0, 0.6D0/ DATA DTRUE3/0.0D0, 0.3D0, 0.7D0, 1.1D0, 1.0D0/ DATA DTRUE5/0.10D0, 2.0D0, 2.0D0, 2.0D0, 2.0D0, @@ -297,6 +307,7 @@ SUBROUTINE CHECK1(SFAC) + 0.03D0, 4.0D0, -0.09D0, 6.0D0, -0.15D0, 7.0D0, + -0.03D0, 3.0D0/ DATA ITRUE2/0, 1, 2, 2, 3/ + DATA ITRUEC/0, 1, 1, 1, 1/ * .. Executable Statements .. DO 80 INCX = 1, 2 DO 60 NP1 = 1, 5 @@ -309,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 @@ -325,13 +340,29 @@ SUBROUTINE CHECK1(SFAC) ELSE IF (ICASE.EQ.10) THEN * .. IDAMAX .. CALL ITEST1(IDAMAX(N,SX,INCX),ITRUE2(NP1)) + DO 100 I = 1, LEN + SX(I) = 42.0D0 + 100 CONTINUE + CALL ITEST1(IDAMAX(N,SX,INCX),ITRUEC(NP1)) ELSE WRITE (NOUT,*) ' Shouldn''t be here in CHECK1' STOP END IF 60 CONTINUE + IF (ICASE.EQ.10) THEN + N = 8 + IX = 1 + DO 120 I = 1, N + SXR(IX) = DVR(I) + IX = IX + INCX + 120 CONTINUE + CALL ITEST1(IDAMAX(N,SXR,INCX),3) + END IF 80 CONTINUE RETURN +* +* End of CHECK1 +* END SUBROUTINE CHECK2(SFAC) * .. Parameters .. @@ -345,7 +376,7 @@ SUBROUTINE CHECK2(SFAC) * .. Local Scalars .. DOUBLE PRECISION SA INTEGER I, J, KI, KN, KNI, KPAR, KSIZE, LENX, LENY, - $ MX, MY + $ LINCX, LINCY, MX, MY * .. Local Arrays .. DOUBLE PRECISION DT10X(7,4,4), DT10Y(7,4,4), DT7(4,4), $ DT8(7,4,4), DX1(7), @@ -354,7 +385,8 @@ SUBROUTINE CHECK2(SFAC) $ DPAR(5,4), DT19X(7,4,16),DT19XA(7,4,4), $ DT19XB(7,4,4), DT19XC(7,4,4),DT19XD(7,4,4), $ DT19Y(7,4,16), DT19YA(7,4,4),DT19YB(7,4,4), - $ DT19YC(7,4,4), DT19YD(7,4,4), DTEMP(5) + $ DT19YC(7,4,4), DT19YD(7,4,4), DTEMP(5), + $ STY0(1), SX0(1), SY0(1) INTEGER INCXS(4), INCYS(4), LENS(4,2), NS(4) * .. External Functions .. DOUBLE PRECISION DDOT, DSDOT @@ -628,6 +660,23 @@ SUBROUTINE CHECK2(SFAC) 60 CONTINUE CALL DCOPY(N,SX,INCX,SY,INCY) CALL STEST(LENY,SY,STY,SSIZE2(1,1),1.0D0) + IF (KI.EQ.1) THEN + SX0(1) = 42.0D0 + SY0(1) = 43.0D0 + IF (N.EQ.0) THEN + STY0(1) = SY0(1) + ELSE + STY0(1) = SX0(1) + END IF + LINCX = INCX + INCX = 0 + LINCY = INCY + INCY = 0 + CALL DCOPY(N,SX0,INCX,SY0,INCY) + CALL STEST(1,SY0,STY0,SSIZE2(1,1),1.0D0) + INCX = LINCX + INCY = LINCY + END IF ELSE IF (ICASE.EQ.6) THEN * .. DSWAP .. CALL DSWAP(N,SX,INCX,SY,INCY) @@ -677,6 +726,9 @@ SUBROUTINE CHECK2(SFAC) 100 CONTINUE 120 CONTINUE RETURN +* +* End of CHECK2 +* END SUBROUTINE CHECK3(SFAC) * .. Parameters .. @@ -883,6 +935,9 @@ SUBROUTINE CHECK3(SFAC) CALL STEST(5,COPYY,MWPSTY,MWPSTY,SFAC) 200 CONTINUE RETURN +* +* End of CHECK3 +* END SUBROUTINE STEST(LEN,SCOMP,STRUE,SSIZE,SFAC) * ********************************* STEST ************************** @@ -939,6 +994,9 @@ SUBROUTINE STEST(LEN,SCOMP,STRUE,SSIZE,SFAC) + ' COMP(I) TRUE(I) DIFFERENCE', + ' SIZE(I)',/1X) 99997 FORMAT (1X,I4,I3,2I5,I3,2D36.8,2D12.4) +* +* End of STEST +* END SUBROUTINE TESTDSDOT(SCOMP,STRUE,SSIZE,SFAC) * ********************************* STEST ************************** @@ -987,6 +1045,9 @@ SUBROUTINE TESTDSDOT(SCOMP,STRUE,SSIZE,SFAC) + ' COMP(I) TRUE(I) DIFFERENCE', + ' SIZE(I)',/1X) 99997 FORMAT (1X,I4,I3,1I5,I3,2E36.8,2E12.4) +* +* End of TESTDSDOT +* END SUBROUTINE STEST1(SCOMP1,STRUE1,SSIZE,SFAC) * ************************* STEST1 ***************************** @@ -1012,6 +1073,9 @@ SUBROUTINE STEST1(SCOMP1,STRUE1,SSIZE,SFAC) CALL STEST(1,SCOMP,STRUE,SSIZE,SFAC) * RETURN +* +* End of STEST1 +* END DOUBLE PRECISION FUNCTION SDIFF(SA,SB) * ********************************* SDIFF ************************** @@ -1022,6 +1086,9 @@ DOUBLE PRECISION FUNCTION SDIFF(SA,SB) * .. Executable Statements .. SDIFF = SA - SB RETURN +* +* End of SDIFF +* END SUBROUTINE ITEST1(ICOMP,ITRUE) * ********************************* ITEST1 ************************* @@ -1063,4 +1130,217 @@ SUBROUTINE ITEST1(ICOMP,ITRUE) + ' COMP TRUE DIFFERENCE', + /1X) 99997 FORMAT (1X,I4,I3,2I5,2I36,I12) +* +* 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/test/sblat1.f b/test/sblat1.f index fe05bbe870..4dc537e2f0 100644 --- a/test/sblat1.f +++ b/test/sblat1.f @@ -30,17 +30,14 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \date April 2012 -* *> \ingroup single_blas_testing * * ===================================================================== PROGRAM SBLAT1 * -* -- Reference BLAS test routine (version 3.8.0) -- +* -- Reference BLAS test routine -- * -- Reference BLAS is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* April 2012 * * ===================================================================== * @@ -91,6 +88,9 @@ PROGRAM SBLAT1 * 99999 FORMAT (' Real BLAS Test Program Results',/1X) 99998 FORMAT (' ----- PASS -----') +* +* End of SBLAT1 +* END SUBROUTINE HEADER * .. Parameters .. @@ -122,6 +122,9 @@ SUBROUTINE HEADER RETURN * 99999 FORMAT (/' Test of subprogram number',I3,12X,A6) +* +* End of HEADER +* END SUBROUTINE CHECK0(SFAC) * .. Parameters .. @@ -238,28 +241,33 @@ SUBROUTINE CHECK0(SFAC) END IF 20 CONTINUE 40 RETURN +* +* End of CHECK0 +* END SUBROUTINE CHECK1(SFAC) * .. Parameters .. INTEGER NOUT - PARAMETER (NOUT=6) + REAL THRESH + PARAMETER (NOUT=6, THRESH=10.0E0) * .. Scalar Arguments .. REAL SFAC * .. Scalars in Common .. INTEGER ICASE, INCX, INCY, N LOGICAL PASS * .. Local Scalars .. - INTEGER I, LEN, NP1 + INTEGER I, IX, LEN, NP1 * .. Local Arrays .. REAL DTRUE1(5), DTRUE3(5), DTRUE5(8,5,2), DV(8,5,2), - + SA(10), STEMP(1), STRUE(8), SX(8) - INTEGER ITRUE2(5) + + DVR(8), SA(10), STEMP(1), STRUE(8), SX(8), + + SXR(15) + INTEGER ITRUE2(5), ITRUEC(5) * .. External Functions .. REAL SASUM, SNRM2 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 .. @@ -280,6 +288,8 @@ SUBROUTINE CHECK1(SFAC) + 0.2E0, 3.0E0, -0.6E0, 5.0E0, 0.3E0, 2.0E0, + 2.0E0, 2.0E0, 0.1E0, 4.0E0, -0.3E0, 6.0E0, + -0.5E0, 7.0E0, -0.1E0, 3.0E0/ + DATA DVR/8.0E0, -7.0E0, 9.0E0, 5.0E0, 9.0E0, 8.0E0, + + 7.0E0, 7.0E0/ DATA DTRUE1/0.0E0, 0.3E0, 0.5E0, 0.7E0, 0.6E0/ DATA DTRUE3/0.0E0, 0.3E0, 0.7E0, 1.1E0, 1.0E0/ DATA DTRUE5/0.10E0, 2.0E0, 2.0E0, 2.0E0, 2.0E0, @@ -297,6 +307,7 @@ SUBROUTINE CHECK1(SFAC) + 0.03E0, 4.0E0, -0.09E0, 6.0E0, -0.15E0, 7.0E0, + -0.03E0, 3.0E0/ DATA ITRUE2/0, 1, 2, 2, 3/ + DATA ITRUEC/0, 1, 1, 1, 1/ * .. Executable Statements .. DO 80 INCX = 1, 2 DO 60 NP1 = 1, 5 @@ -309,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 @@ -325,13 +340,29 @@ SUBROUTINE CHECK1(SFAC) ELSE IF (ICASE.EQ.10) THEN * .. ISAMAX .. CALL ITEST1(ISAMAX(N,SX,INCX),ITRUE2(NP1)) + DO 100 I = 1, LEN + SX(I) = 42.0E0 + 100 CONTINUE + CALL ITEST1(ISAMAX(N,SX,INCX),ITRUEC(NP1)) ELSE WRITE (NOUT,*) ' Shouldn''t be here in CHECK1' STOP END IF 60 CONTINUE + IF (ICASE.EQ.10) THEN + N = 8 + IX = 1 + DO 120 I = 1, N + SXR(IX) = DVR(I) + IX = IX + INCX + 120 CONTINUE + CALL ITEST1(ISAMAX(N,SXR,INCX),3) + END IF 80 CONTINUE RETURN +* +* End of CHECK1 +* END SUBROUTINE CHECK2(SFAC) * .. Parameters .. @@ -345,7 +376,7 @@ SUBROUTINE CHECK2(SFAC) * .. Local Scalars .. REAL SA INTEGER I, J, KI, KN, KNI, KPAR, KSIZE, LENX, LENY, - $ MX, MY + $ LINCX, LINCY, MX, MY * .. Local Arrays .. REAL DT10X(7,4,4), DT10Y(7,4,4), DT7(4,4), $ DT8(7,4,4), DX1(7), @@ -355,7 +386,7 @@ SUBROUTINE CHECK2(SFAC) $ DT19XB(7,4,4), DT19XC(7,4,4),DT19XD(7,4,4), $ DT19Y(7,4,16), DT19YA(7,4,4),DT19YB(7,4,4), $ DT19YC(7,4,4), DT19YD(7,4,4), DTEMP(5), - $ ST7B(4,4) + $ ST7B(4,4), STY0(1), SX0(1), SY0(1) INTEGER INCXS(4), INCYS(4), LENS(4,2), NS(4) * .. External Functions .. REAL SDOT, SDSDOT @@ -631,6 +662,23 @@ SUBROUTINE CHECK2(SFAC) 60 CONTINUE CALL SCOPY(N,SX,INCX,SY,INCY) CALL STEST(LENY,SY,STY,SSIZE2(1,1),1.0E0) + IF (KI.EQ.1) THEN + SX0(1) = 42.0E0 + SY0(1) = 43.0E0 + IF (N.EQ.0) THEN + STY0(1) = SY0(1) + ELSE + STY0(1) = SX0(1) + END IF + LINCX = INCX + INCX = 0 + LINCY = INCY + INCY = 0 + CALL SCOPY(N,SX0,INCX,SY0,INCY) + CALL STEST(1,SY0,STY0,SSIZE2(1,1),1.0E0) + INCX = LINCX + INCY = LINCY + END IF ELSE IF (ICASE.EQ.6) THEN * .. SSWAP .. CALL SSWAP(N,SX,INCX,SY,INCY) @@ -680,6 +728,9 @@ SUBROUTINE CHECK2(SFAC) 100 CONTINUE 120 CONTINUE RETURN +* +* End of CHECK2 +* END SUBROUTINE CHECK3(SFAC) * .. Parameters .. @@ -886,6 +937,9 @@ SUBROUTINE CHECK3(SFAC) CALL STEST(5,COPYY,MWPSTY,MWPSTY,SFAC) 200 CONTINUE RETURN +* +* End of CHECK3 +* END SUBROUTINE STEST(LEN,SCOMP,STRUE,SSIZE,SFAC) * ********************************* STEST ************************** @@ -942,6 +996,9 @@ SUBROUTINE STEST(LEN,SCOMP,STRUE,SSIZE,SFAC) + ' COMP(I) TRUE(I) DIFFERENCE', + ' SIZE(I)',/1X) 99997 FORMAT (1X,I4,I3,2I5,I3,2E36.8,2E12.4) +* +* End of STEST +* END SUBROUTINE STEST1(SCOMP1,STRUE1,SSIZE,SFAC) * ************************* STEST1 ***************************** @@ -967,6 +1024,9 @@ SUBROUTINE STEST1(SCOMP1,STRUE1,SSIZE,SFAC) CALL STEST(1,SCOMP,STRUE,SSIZE,SFAC) * RETURN +* +* End of STEST1 +* END REAL FUNCTION SDIFF(SA,SB) * ********************************* SDIFF ************************** @@ -977,6 +1037,9 @@ REAL FUNCTION SDIFF(SA,SB) * .. Executable Statements .. SDIFF = SA - SB RETURN +* +* End of SDIFF +* END SUBROUTINE ITEST1(ICOMP,ITRUE) * ********************************* ITEST1 ************************* @@ -1018,4 +1081,218 @@ SUBROUTINE ITEST1(ICOMP,ITRUE) + ' COMP TRUE DIFFERENCE', + /1X) 99997 FORMAT (1X,I4,I3,2I5,2I36,I12) +* +* 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/test/zblat1.f b/test/zblat1.f index 2d7b884902..ef6deff921 100644 --- a/test/zblat1.f +++ b/test/zblat1.f @@ -30,17 +30,14 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \date April 2012 -* *> \ingroup complex16_blas_testing * * ===================================================================== PROGRAM ZBLAT1 * -* -- Reference BLAS test routine (version 3.7.0) -- +* -- Reference BLAS test routine -- * -- Reference BLAS is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* April 2012 * * ===================================================================== * @@ -86,6 +83,9 @@ PROGRAM ZBLAT1 * 99999 FORMAT (' Complex BLAS Test Program Results',/1X) 99998 FORMAT (' ----- PASS -----') +* +* End of ZBLAT1 +* END SUBROUTINE HEADER * .. Parameters .. @@ -114,11 +114,15 @@ SUBROUTINE HEADER RETURN * 99999 FORMAT (/' Test of subprogram number',I3,12X,A6) +* +* End of HEADER +* END 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 .. @@ -127,18 +131,18 @@ SUBROUTINE CHECK1(SFAC) * .. Local Scalars .. COMPLEX*16 CA DOUBLE PRECISION SA - INTEGER I, J, LEN, NP1 + INTEGER I, IX, J, LEN, NP1 * .. Local Arrays .. - COMPLEX*16 CTRUE5(8,5,2), CTRUE6(8,5,2), CV(8,5,2), CX(8), - + MWPCS(5), MWPCT(5) + COMPLEX*16 CTRUE5(8,5,2), CTRUE6(8,5,2), CV(8,5,2), CVR(8), + + CX(8), CXR(15), MWPCS(5), MWPCT(5) DOUBLE PRECISION STRUE2(5), STRUE4(5) - INTEGER ITRUE3(5) + INTEGER ITRUE3(5), ITRUEC(5) * .. External Functions .. DOUBLE PRECISION DZASUM, DZNRM2 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 .. @@ -173,6 +177,9 @@ SUBROUTINE CHECK1(SFAC) + (7.0D0,2.0D0), (0.3D0,0.1D0), (5.0D0,8.0D0), + (0.5D0,0.0D0), (6.0D0,9.0D0), (0.0D0,0.5D0), + (8.0D0,3.0D0), (0.0D0,0.2D0), (9.0D0,4.0D0)/ + DATA CVR/(8.0D0,8.0D0), (-7.0D0,-7.0D0), + + (9.0D0,9.0D0), (5.0D0,5.0D0), (9.0D0,9.0D0), + + (8.0D0,8.0D0), (7.0D0,7.0D0), (7.0D0,7.0D0)/ DATA STRUE2/0.0D0, 0.5D0, 0.6D0, 0.7D0, 0.8D0/ DATA STRUE4/0.0D0, 0.7D0, 1.0D0, 1.3D0, 1.6D0/ DATA ((CTRUE5(I,J,1),I=1,8),J=1,5)/(0.1D0,0.1D0), @@ -238,6 +245,7 @@ SUBROUTINE CHECK1(SFAC) + (0.15D0,0.00D0), (6.0D0,9.0D0), (0.00D0,0.15D0), + (8.0D0,3.0D0), (0.00D0,0.06D0), (9.0D0,4.0D0)/ DATA ITRUE3/0, 1, 2, 2, 2/ + DATA ITRUEC/0, 1, 1, 1, 1/ * .. Executable Statements .. DO 60 INCX = 1, 2 DO 40 NP1 = 1, 5 @@ -249,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 @@ -268,12 +280,25 @@ SUBROUTINE CHECK1(SFAC) ELSE IF (ICASE.EQ.10) THEN * .. IZAMAX .. CALL ITEST1(IZAMAX(N,CX,INCX),ITRUE3(NP1)) + DO 160 I = 1, LEN + CX(I) = (42.0D0,43.0D0) + 160 CONTINUE + CALL ITEST1(IZAMAX(N,CX,INCX),ITRUEC(NP1)) ELSE WRITE (NOUT,*) ' Shouldn''t be here in CHECK1' STOP END IF * 40 CONTINUE + IF (ICASE.EQ.10) THEN + N = 8 + IX = 1 + DO 180 I = 1, N + CXR(IX) = CVR(I) + IX = IX + INCX + 180 CONTINUE + CALL ITEST1(IZAMAX(N,CXR,INCX),3) + END IF 60 CONTINUE * INCX = 1 @@ -315,6 +340,9 @@ SUBROUTINE CHECK1(SFAC) CALL CTEST(5,CX,MWPCT,MWPCS,SFAC) END IF RETURN +* +* End of CHECK1 +* END SUBROUTINE CHECK2(SFAC) * .. Parameters .. @@ -327,11 +355,13 @@ SUBROUTINE CHECK2(SFAC) LOGICAL PASS * .. Local Scalars .. COMPLEX*16 CA - INTEGER I, J, KI, KN, KSIZE, LENX, LENY, MX, MY + INTEGER I, J, KI, KN, KSIZE, LENX, LENY, LINCX, LINCY, + + MX, MY * .. Local Arrays .. COMPLEX*16 CDOT(1), CSIZE1(4), CSIZE2(7,2), CSIZE3(14), + CT10X(7,4,4), CT10Y(7,4,4), CT6(4,4), CT7(4,4), - + CT8(7,4,4), CX(7), CX1(7), CY(7), CY1(7) + + CT8(7,4,4), CTY0(1), CX(7), CX0(1), CX1(7), + + CY(7), CY0(1), CY1(7) INTEGER INCXS(4), INCYS(4), LENS(4,2), NS(4) * .. External Functions .. COMPLEX*16 ZDOTC, ZDOTU @@ -546,6 +576,23 @@ SUBROUTINE CHECK2(SFAC) * .. ZCOPY .. CALL ZCOPY(N,CX,INCX,CY,INCY) CALL CTEST(LENY,CY,CT10Y(1,KN,KI),CSIZE3,1.0D0) + IF (KI.EQ.1) THEN + CX0(1) = (42.0D0,43.0D0) + CY0(1) = (44.0D0,45.0D0) + IF (N.EQ.0) THEN + CTY0(1) = CY0(1) + ELSE + CTY0(1) = CX0(1) + END IF + LINCX = INCX + INCX = 0 + LINCY = INCY + INCY = 0 + CALL ZCOPY(N,CX0,INCX,CY0,INCY) + CALL CTEST(1,CY0,CTY0,CSIZE3,1.0D0) + INCX = LINCX + INCY = LINCY + END IF ELSE IF (ICASE.EQ.5) THEN * .. ZSWAP .. CALL ZSWAP(N,CX,INCX,CY,INCY) @@ -559,6 +606,9 @@ SUBROUTINE CHECK2(SFAC) 40 CONTINUE 60 CONTINUE RETURN +* +* End of CHECK2 +* END SUBROUTINE STEST(LEN,SCOMP,STRUE,SSIZE,SFAC) * ********************************* STEST ************************** @@ -615,6 +665,9 @@ SUBROUTINE STEST(LEN,SCOMP,STRUE,SSIZE,SFAC) + ' COMP(I) TRUE(I) DIFFERENCE', + ' SIZE(I)',/1X) 99997 FORMAT (1X,I4,I3,3I5,I3,2D36.8,2D12.4) +* +* End of STEST +* END SUBROUTINE STEST1(SCOMP1,STRUE1,SSIZE,SFAC) * ************************* STEST1 ***************************** @@ -640,6 +693,9 @@ SUBROUTINE STEST1(SCOMP1,STRUE1,SSIZE,SFAC) CALL STEST(1,SCOMP,STRUE,SSIZE,SFAC) * RETURN +* +* End of STEST1 +* END DOUBLE PRECISION FUNCTION SDIFF(SA,SB) * ********************************* SDIFF ************************** @@ -650,6 +706,9 @@ DOUBLE PRECISION FUNCTION SDIFF(SA,SB) * .. Executable Statements .. SDIFF = SA - SB RETURN +* +* End of SDIFF +* END SUBROUTINE CTEST(LEN,CCOMP,CTRUE,CSIZE,SFAC) * **************************** CTEST ***************************** @@ -681,6 +740,9 @@ SUBROUTINE CTEST(LEN,CCOMP,CTRUE,CSIZE,SFAC) * CALL STEST(2*LEN,SCOMP,STRUE,SSIZE,SFAC) RETURN +* +* End of CTEST +* END SUBROUTINE ITEST1(ICOMP,ITRUE) * ********************************* ITEST1 ************************* @@ -721,4 +783,232 @@ SUBROUTINE ITEST1(ICOMP,ITRUE) + ' COMP TRUE DIFFERENCE', + /1X) 99997 FORMAT (1X,I4,I3,3I5,2I36,I12) +* +* 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