From 11d59375b7e3ed8e80015ef9b0e5812e76a140bb Mon Sep 17 00:00:00 2001 From: Simon Maertens Date: Mon, 14 Sep 2026 17:22:11 +0100 Subject: [PATCH] BLAS TESTING: exercise the Level 1 unrolled loops and the ?NRM2 threshold ?DOT, ?COPY, ?SCAL and ?ASUM return early below their unroll factors of 5, 7, 5 and 6, so no case at N <= 4 reached the unrolled body: CHECK1 and CHECK2 gain a pass at N = 16. ?B1NRM2 gains a value at SQRT(SAFMIN), where ?NRM2 switches to its small accumulator, and one just below it, so that the small sum can outweigh the mid-range one. Co-Authored-By: Claude Opus 5 (1M context) --- BLAS/TESTING/cblat1.f | 10 ++- BLAS/TESTING/dblat1.f | 149 ++++++++++++++++++++++++++++++++++++++++-- BLAS/TESTING/sblat1.f | 149 ++++++++++++++++++++++++++++++++++++++++-- BLAS/TESTING/zblat1.f | 10 ++- 4 files changed, 304 insertions(+), 14 deletions(-) diff --git a/BLAS/TESTING/cblat1.f b/BLAS/TESTING/cblat1.f index f07e4e189..fc38a4b6b 100644 --- a/BLAS/TESTING/cblat1.f +++ b/BLAS/TESTING/cblat1.f @@ -899,7 +899,7 @@ SUBROUTINE CB1NRM2(N,INCX,THRESH) * ===================================================================== * .. Parameters .. INTEGER NMAX, NOUT, NV - PARAMETER (NMAX=20, NOUT=6, NV=10) + PARAMETER (NMAX=20, NOUT=6, NV=12) 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) @@ -939,6 +939,14 @@ SUBROUTINE CB1NRM2(N,INCX,THRESH) VALUES(8) = SAFMAX VALUES(9) = SXVALS(V0,2) VALUES(10) = SXVALS(V0,3) +* SQRT(SAFMIN) is exactly the threshold at which ?NRM2 switches to +* its small accumulator: SAFMIN is RADIX**(MINEXPONENT-1) and the +* threshold is that exponent halved. One value at the threshold +* and one just below it let the small part of a vector outweigh +* its mid-range part, which is what the arm ordering the two sums +* needs. A wider pair leaves them too far apart for that. + VALUES(11) = SQRT(SAFMIN) + VALUES(12) = SQRT(SAFMIN)*(ONE-ULP) ROGUE = CMPLX(1234.5678E+0,-1234.5678E+0) FIRST = .TRUE. * diff --git a/BLAS/TESTING/dblat1.f b/BLAS/TESTING/dblat1.f index 43590bd23..92539d373 100644 --- a/BLAS/TESTING/dblat1.f +++ b/BLAS/TESTING/dblat1.f @@ -287,10 +287,13 @@ SUBROUTINE CHECK1(SFAC) LOGICAL PASS * .. Local Scalars .. INTEGER I, IX, LEN, NP1 + INTEGER NLONG + PARAMETER (NLONG=16) + DOUBLE PRECISION SLSUM, SLSIZE(1) * .. Local Arrays .. DOUBLE PRECISION DTRUE1(5), DTRUE3(5), DTRUE5(8,5,2), DV(8,5,2), + DVR(8), SA(10), STEMP(1), STRUE(8), SX(8), - + SXR(15) + + SXR(15), SXL(NLONG), STL(NLONG) INTEGER ITRUE2(5), ITRUEC(5) * .. External Functions .. DOUBLE PRECISION DASUM, DNRM2 @@ -299,7 +302,7 @@ SUBROUTINE CHECK1(SFAC) * .. External Subroutines .. EXTERNAL ITEST1, DB1NRM2, DSCAL, STEST, STEST1 * .. Intrinsic Functions .. - INTRINSIC MAX + INTRINSIC ABS, DBLE, MAX, MOD, SQRT * .. Common blocks .. COMMON /COMBLA/ICASE, N, INCX, INCY, PASS * .. Data statements .. @@ -389,6 +392,52 @@ SUBROUTINE CHECK1(SFAC) CALL ITEST1(IDAMAX(N,SXR,INCX),3) END IF 80 CONTINUE +* +* The Level 1 routines process their vectors in unrolled blocks of +* up to seven elements once N reaches the unroll factor, and handle +* only the remainder below it. Every case above stops at N = 4, so +* the unrolled bodies never ran. Repeat the tests at N = 16, which +* leaves a nonzero remainder for every unroll factor in use. The +* data are small integers, so the unrolled and sequential sums agree +* exactly and the expected values can be formed here. +* + IF (ICASE.GE.7 .AND. ICASE.LE.10) THEN + INCX = 1 + N = NLONG + DO 140 I = 1, N + SXL(I) = DBLE(MOD(I,7) - 3) + 140 CONTINUE +* Give the search a unique maximum away from either end. + SXL(11) = -9D0 + IF (ICASE.EQ.7) THEN +* .. ?NRM2 .. + SLSUM = 0D0 + DO 150 I = 1, N + SLSUM = SLSUM + SXL(I)*SXL(I) + 150 CONTINUE + SLSUM = SQRT(SLSUM) + SLSIZE(1) = SLSUM + CALL STEST1(DNRM2(N,SXL,INCX),SLSUM,SLSIZE,SFAC) + ELSE IF (ICASE.EQ.8) THEN +* .. ?ASUM .. + SLSUM = 0D0 + DO 160 I = 1, N + SLSUM = SLSUM + ABS(SXL(I)) + 160 CONTINUE + SLSIZE(1) = SLSUM + CALL STEST1(DASUM(N,SXL,INCX),SLSUM,SLSIZE,SFAC) + ELSE IF (ICASE.EQ.9) THEN +* .. ?SCAL .. + DO 170 I = 1, N + STL(I) = SA(1)*SXL(I) + 170 CONTINUE + CALL DSCAL(N,SA(1),SXL,INCX) + CALL STEST(N,SXL,STL,STL,SFAC) + ELSE +* .. I?AMAX .. + CALL ITEST1(IDAMAX(N,SXL,INCX),11) + END IF + END IF RETURN * * End of CHECK1 @@ -405,7 +454,9 @@ SUBROUTINE CHECK2(SFAC) INTEGER ICASE, INCX, INCY, N LOGICAL PASS * .. Local Scalars .. - DOUBLE PRECISION SA, SB + DOUBLE PRECISION SA, SB, SLSUM + INTEGER NLONG + PARAMETER (NLONG=16) INTEGER I, J, KI, KN, KNI, KPAR, KSIZE, LENX, LENY, $ LINCX, LINCY, MX, MY * .. Local Arrays .. @@ -417,7 +468,9 @@ 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), - $ STY0(1), SX0(1), SY0(1), DT20(7,4,4) + $ STY0(1), SX0(1), SY0(1), DT20(7,4,4), + $ SXL(NLONG), SYL(NLONG), STXL(NLONG), + $ STYL(NLONG), SLSIZE(1) INTEGER INCXS(4), INCYS(4), LENS(4,2), NS(4) * .. External Functions .. DOUBLE PRECISION DDOT, DSDOT @@ -426,7 +479,7 @@ SUBROUTINE CHECK2(SFAC) EXTERNAL DAXPY, DAXPBY, DCOPY, DROTM, DSWAP, STEST, $ STEST1, TESTDSDOT * .. Intrinsic Functions .. - INTRINSIC ABS, MIN + INTRINSIC ABS, DBLE, MIN, MOD * .. Common blocks .. COMMON /COMBLA/ICASE, N, INCX, INCY, PASS * .. Data statements .. @@ -786,6 +839,82 @@ SUBROUTINE CHECK2(SFAC) END IF 100 CONTINUE 120 CONTINUE +* +* As in CHECK1, the two-vector routines unroll once N reaches the +* unroll factor and the cases above stop at N = 4. Repeat them at +* N = 16 on small integer data, so the expected values are exact. +* + IF (ICASE.EQ.1 .OR. ICASE.EQ.2 .OR. ICASE.EQ.5 .OR. + + ICASE.EQ.6 .OR. ICASE.EQ.12 .OR. ICASE.EQ.13 .OR. + + ICASE.EQ.14) THEN + INCX = 1 + INCY = 1 + N = NLONG + DO 200 I = 1, N + SXL(I) = DBLE(MOD(I,7) - 3) + SYL(I) = DBLE(MOD(I,5) - 2) + 200 CONTINUE + IF (ICASE.EQ.1 .OR. ICASE.EQ.13) THEN +* .. ?DOT / ?SDOT .. + SLSUM = 0D0 + DO 210 I = 1, N + SLSUM = SLSUM + SXL(I)*SYL(I) + 210 CONTINUE + SLSIZE(1) = SLSUM + IF (ICASE.EQ.1) THEN + CALL STEST1(DDOT(N,SXL,INCX,SYL,INCY),SLSUM,SLSIZE, + + SFAC) + ELSE + CALL TESTDSDOT(REAL(DSDOT(N,REAL(SXL),INCX,REAL(SYL), + + INCY)),REAL(SLSUM),REAL(SLSUM),.3125E-1) + END IF + ELSE + DO 220 I = 1, N + STXL(I) = SXL(I) + STYL(I) = SYL(I) + 220 CONTINUE + IF (ICASE.EQ.2) THEN +* .. ?AXPY .. + DO 230 I = 1, N + STYL(I) = SYL(I) + SA*SXL(I) + 230 CONTINUE + CALL DAXPY(N,SA,SXL,INCX,SYL,INCY) + ELSE IF (ICASE.EQ.5) THEN +* .. ?COPY .. + DO 240 I = 1, N + STYL(I) = SXL(I) + 240 CONTINUE + CALL DCOPY(N,SXL,INCX,SYL,INCY) + ELSE IF (ICASE.EQ.6) THEN +* .. ?SWAP .. + DO 250 I = 1, N + STXL(I) = SYL(I) + STYL(I) = SXL(I) + 250 CONTINUE + CALL DSWAP(N,SXL,INCX,SYL,INCY) + ELSE IF (ICASE.EQ.12) THEN +* .. ?ROTM, with DFLAG = -1 so the full H is applied .. + DTEMP(1) = -1D0 + DTEMP(2) = 2D0 + DTEMP(3) = 3D0 + DTEMP(4) = 4D0 + DTEMP(5) = 5D0 + DO 260 I = 1, N + STXL(I) = SXL(I)*DTEMP(2) + SYL(I)*DTEMP(4) + STYL(I) = SXL(I)*DTEMP(3) + SYL(I)*DTEMP(5) + 260 CONTINUE + CALL DROTM(N,SXL,INCX,SYL,INCY,DTEMP) + ELSE +* .. ?AXPBY .. + DO 270 I = 1, N + STYL(I) = SA*SXL(I) + SB*SYL(I) + 270 CONTINUE + CALL DAXPBY(N,SA,SXL,INCX,SB,SYL,INCY) + END IF + CALL STEST(N,SXL,STXL,STXL,SFAC) + CALL STEST(N,SYL,STYL,STYL,SFAC) + END IF + END IF RETURN * * End of CHECK2 @@ -1236,7 +1365,7 @@ SUBROUTINE DB1NRM2(N,INCX,THRESH) * ===================================================================== * .. Parameters .. INTEGER NMAX, NOUT, NV - PARAMETER (NMAX=20, NOUT=6, NV=10) + PARAMETER (NMAX=20, NOUT=6, NV=12) DOUBLE PRECISION HALF, ONE, TWO, ZERO PARAMETER (HALF=0.5D+0, ONE=1.0D+0, TWO= 2.0D+0, & ZERO=0.0D+0) @@ -1274,6 +1403,14 @@ SUBROUTINE DB1NRM2(N,INCX,THRESH) VALUES(8) = SAFMAX VALUES(9) = DXVALS(V0,2) VALUES(10) = DXVALS(V0,3) +* SQRT(SAFMIN) is exactly the threshold at which ?NRM2 switches to +* its small accumulator: SAFMIN is RADIX**(MINEXPONENT-1) and the +* threshold is that exponent halved. One value at the threshold +* and one just below it let the small part of a vector outweigh +* its mid-range part, which is what the arm ordering the two sums +* needs. A wider pair leaves them too far apart for that. + VALUES(11) = SQRT(SAFMIN) + VALUES(12) = SQRT(SAFMIN)*(ONE-ULP) ROGUE = -1234.5678D+0 FIRST = .TRUE. * diff --git a/BLAS/TESTING/sblat1.f b/BLAS/TESTING/sblat1.f index 4ea895abb..d226f72d9 100644 --- a/BLAS/TESTING/sblat1.f +++ b/BLAS/TESTING/sblat1.f @@ -285,10 +285,13 @@ SUBROUTINE CHECK1(SFAC) LOGICAL PASS * .. Local Scalars .. INTEGER I, IX, LEN, NP1 + INTEGER NLONG + PARAMETER (NLONG=16) + REAL SLSUM, SLSIZE(1) * .. Local Arrays .. REAL DTRUE1(5), DTRUE3(5), DTRUE5(8,5,2), DV(8,5,2), + DVR(8), SA(10), STEMP(1), STRUE(8), SX(8), - + SXR(15) + + SXR(15), SXL(NLONG), STL(NLONG) INTEGER ITRUE2(5), ITRUEC(5) * .. External Functions .. REAL SASUM, SNRM2 @@ -297,7 +300,7 @@ SUBROUTINE CHECK1(SFAC) * .. External Subroutines .. EXTERNAL ITEST1, SB1NRM2, SSCAL, STEST, STEST1 * .. Intrinsic Functions .. - INTRINSIC MAX + INTRINSIC ABS, MAX, MOD, REAL, SQRT * .. Common blocks .. COMMON /COMBLA/ICASE, N, INCX, INCY, PASS * .. Data statements .. @@ -387,6 +390,52 @@ SUBROUTINE CHECK1(SFAC) CALL ITEST1(ISAMAX(N,SXR,INCX),3) END IF 80 CONTINUE +* +* The Level 1 routines process their vectors in unrolled blocks of +* up to seven elements once N reaches the unroll factor, and handle +* only the remainder below it. Every case above stops at N = 4, so +* the unrolled bodies never ran. Repeat the tests at N = 16, which +* leaves a nonzero remainder for every unroll factor in use. The +* data are small integers, so the unrolled and sequential sums agree +* exactly and the expected values can be formed here. +* + IF (ICASE.GE.7 .AND. ICASE.LE.10) THEN + INCX = 1 + N = NLONG + DO 140 I = 1, N + SXL(I) = REAL(MOD(I,7) - 3) + 140 CONTINUE +* Give the search a unique maximum away from either end. + SXL(11) = -9E0 + IF (ICASE.EQ.7) THEN +* .. ?NRM2 .. + SLSUM = 0E0 + DO 150 I = 1, N + SLSUM = SLSUM + SXL(I)*SXL(I) + 150 CONTINUE + SLSUM = SQRT(SLSUM) + SLSIZE(1) = SLSUM + CALL STEST1(SNRM2(N,SXL,INCX),SLSUM,SLSIZE,SFAC) + ELSE IF (ICASE.EQ.8) THEN +* .. ?ASUM .. + SLSUM = 0E0 + DO 160 I = 1, N + SLSUM = SLSUM + ABS(SXL(I)) + 160 CONTINUE + SLSIZE(1) = SLSUM + CALL STEST1(SASUM(N,SXL,INCX),SLSUM,SLSIZE,SFAC) + ELSE IF (ICASE.EQ.9) THEN +* .. ?SCAL .. + DO 170 I = 1, N + STL(I) = SA(1)*SXL(I) + 170 CONTINUE + CALL SSCAL(N,SA(1),SXL,INCX) + CALL STEST(N,SXL,STL,STL,SFAC) + ELSE +* .. I?AMAX .. + CALL ITEST1(ISAMAX(N,SXL,INCX),11) + END IF + END IF RETURN * * End of CHECK1 @@ -402,7 +451,9 @@ SUBROUTINE CHECK2(SFAC) INTEGER ICASE, INCX, INCY, N LOGICAL PASS * .. Local Scalars .. - REAL SA,SB + REAL SA,SB, SLSUM + INTEGER NLONG + PARAMETER (NLONG=16) INTEGER I, J, KI, KN, KNI, KPAR, KSIZE, LENX, LENY, $ LINCX, LINCY, MX, MY * .. Local Arrays .. @@ -414,7 +465,9 @@ 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), STY0(1), SX0(1), SY0(1), DT20(7,4,4) + $ ST7B(4,4), STY0(1), SX0(1), SY0(1), DT20(7,4,4), + $ SXL(NLONG), SYL(NLONG), STXL(NLONG), + $ STYL(NLONG), SLSIZE(1) INTEGER INCXS(4), INCYS(4), LENS(4,2), NS(4) * .. External Functions .. REAL SDOT, SDSDOT @@ -422,7 +475,7 @@ SUBROUTINE CHECK2(SFAC) * .. External Subroutines .. EXTERNAL SAXPY, SAXPBY,SCOPY, SROTM, SSWAP, STEST, STEST1 * .. Intrinsic Functions .. - INTRINSIC ABS, MIN + INTRINSIC ABS, MIN, MOD, REAL * .. Common blocks .. COMMON /COMBLA/ICASE, N, INCX, INCY, PASS * .. Data statements .. @@ -785,6 +838,82 @@ SUBROUTINE CHECK2(SFAC) END IF 100 CONTINUE 120 CONTINUE +* +* As in CHECK1, the two-vector routines unroll once N reaches the +* unroll factor and the cases above stop at N = 4. Repeat them at +* N = 16 on small integer data, so the expected values are exact. +* + IF (ICASE.EQ.1 .OR. ICASE.EQ.2 .OR. ICASE.EQ.5 .OR. + + ICASE.EQ.6 .OR. ICASE.EQ.12 .OR. ICASE.EQ.13 .OR. + + ICASE.EQ.14) THEN + INCX = 1 + INCY = 1 + N = NLONG + DO 200 I = 1, N + SXL(I) = REAL(MOD(I,7) - 3) + SYL(I) = REAL(MOD(I,5) - 2) + 200 CONTINUE + IF (ICASE.EQ.1 .OR. ICASE.EQ.13) THEN +* .. ?DOT / ?SDOT .. + SLSUM = 0E0 + DO 210 I = 1, N + SLSUM = SLSUM + SXL(I)*SYL(I) + 210 CONTINUE + SLSIZE(1) = SLSUM + IF (ICASE.EQ.1) THEN + CALL STEST1(SDOT(N,SXL,INCX,SYL,INCY),SLSUM,SLSIZE, + + SFAC) + ELSE + CALL STEST1(SDSDOT(N,0.0E0,SXL,INCX,SYL,INCY),SLSUM, + + SLSIZE,SFAC) + END IF + ELSE + DO 220 I = 1, N + STXL(I) = SXL(I) + STYL(I) = SYL(I) + 220 CONTINUE + IF (ICASE.EQ.2) THEN +* .. ?AXPY .. + DO 230 I = 1, N + STYL(I) = SYL(I) + SA*SXL(I) + 230 CONTINUE + CALL SAXPY(N,SA,SXL,INCX,SYL,INCY) + ELSE IF (ICASE.EQ.5) THEN +* .. ?COPY .. + DO 240 I = 1, N + STYL(I) = SXL(I) + 240 CONTINUE + CALL SCOPY(N,SXL,INCX,SYL,INCY) + ELSE IF (ICASE.EQ.6) THEN +* .. ?SWAP .. + DO 250 I = 1, N + STXL(I) = SYL(I) + STYL(I) = SXL(I) + 250 CONTINUE + CALL SSWAP(N,SXL,INCX,SYL,INCY) + ELSE IF (ICASE.EQ.12) THEN +* .. ?ROTM, with DFLAG = -1 so the full H is applied .. + DTEMP(1) = -1E0 + DTEMP(2) = 2E0 + DTEMP(3) = 3E0 + DTEMP(4) = 4E0 + DTEMP(5) = 5E0 + DO 260 I = 1, N + STXL(I) = SXL(I)*DTEMP(2) + SYL(I)*DTEMP(4) + STYL(I) = SXL(I)*DTEMP(3) + SYL(I)*DTEMP(5) + 260 CONTINUE + CALL SROTM(N,SXL,INCX,SYL,INCY,DTEMP) + ELSE +* .. ?AXPBY .. + DO 270 I = 1, N + STYL(I) = SA*SXL(I) + SB*SYL(I) + 270 CONTINUE + CALL SAXPBY(N,SA,SXL,INCX,SB,SYL,INCY) + END IF + CALL STEST(N,SXL,STXL,STXL,SFAC) + CALL STEST(N,SYL,STYL,STYL,SFAC) + END IF + END IF RETURN * * End of CHECK2 @@ -1175,7 +1304,7 @@ SUBROUTINE SB1NRM2(N,INCX,THRESH) * ===================================================================== * .. Parameters .. INTEGER NMAX, NOUT, NV - PARAMETER (NMAX=20, NOUT=6, NV=10) + PARAMETER (NMAX=20, NOUT=6, NV=12) REAL HALF, ONE, TWO, ZERO PARAMETER (HALF=0.5E+0, ONE=1.0E+0, TWO= 2.0E+0, & ZERO=0.0E+0) @@ -1213,6 +1342,14 @@ SUBROUTINE SB1NRM2(N,INCX,THRESH) VALUES(8) = SAFMAX VALUES(9) = SXVALS(V0,2) VALUES(10) = SXVALS(V0,3) +* SQRT(SAFMIN) is exactly the threshold at which ?NRM2 switches to +* its small accumulator: SAFMIN is RADIX**(MINEXPONENT-1) and the +* threshold is that exponent halved. One value at the threshold +* and one just below it let the small part of a vector outweigh +* its mid-range part, which is what the arm ordering the two sums +* needs. A wider pair leaves them too far apart for that. + VALUES(11) = SQRT(SAFMIN) + VALUES(12) = SQRT(SAFMIN)*(ONE-ULP) ROGUE = -1234.5678E+0 FIRST = .TRUE. * diff --git a/BLAS/TESTING/zblat1.f b/BLAS/TESTING/zblat1.f index 7742590d7..fa1a47b45 100644 --- a/BLAS/TESTING/zblat1.f +++ b/BLAS/TESTING/zblat1.f @@ -899,7 +899,7 @@ SUBROUTINE ZB1NRM2(N,INCX,THRESH) * ===================================================================== * .. Parameters .. INTEGER NMAX, NOUT, NV - PARAMETER (NMAX=20, NOUT=6, NV=10) + PARAMETER (NMAX=20, NOUT=6, NV=12) 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) @@ -939,6 +939,14 @@ SUBROUTINE ZB1NRM2(N,INCX,THRESH) VALUES(8) = SAFMAX VALUES(9) = DXVALS(V0,2) VALUES(10) = DXVALS(V0,3) +* SQRT(SAFMIN) is exactly the threshold at which ?NRM2 switches to +* its small accumulator: SAFMIN is RADIX**(MINEXPONENT-1) and the +* threshold is that exponent halved. One value at the threshold +* and one just below it let the small part of a vector outweigh +* its mid-range part, which is what the arm ordering the two sums +* needs. A wider pair leaves them too far apart for that. + VALUES(11) = SQRT(SAFMIN) + VALUES(12) = SQRT(SAFMIN)*(ONE-ULP) ROGUE = DCMPLX(1234.5678D+0,-1234.5678D+0) FIRST = .TRUE. *