diff --git a/BLAS/TESTING/cblat1.f b/BLAS/TESTING/cblat1.f index f07e4e189..c968f5a5e 100644 --- a/BLAS/TESTING/cblat1.f +++ b/BLAS/TESTING/cblat1.f @@ -65,7 +65,7 @@ PROGRAM CBLAT1 * .. Executable Statements .. CALL CPU_TIME( S1 ) WRITE (NOUT,99999) - DO 20 IC = 1, 11 + DO 20 IC = 1, 13 ICASE = IC CALL HEADER * @@ -80,7 +80,10 @@ PROGRAM CBLAT1 INCX = 9999 INCY = 9999 MODE = 9999 - IF (ICASE.LE.5 .OR. ICASE.EQ.11) THEN + IF (ICASE.EQ.12) THEN + CALL CHECK0(SFAC) + ELSE IF (ICASE.LE.5 .OR. ICASE.EQ.11 .OR. + + ICASE.EQ.13) THEN CALL CHECK2(SFAC) ELSE IF (ICASE.GE.6) THEN CALL CHECK1(SFAC) @@ -112,7 +115,7 @@ SUBROUTINE HEADER INTEGER ICASE, INCX, INCY, MODE, N LOGICAL PASS * .. Local Arrays .. - CHARACTER*6 L(11) + CHARACTER*6 L(13) * .. Common blocks .. COMMON /COMBLA/ICASE, N, INCX, INCY, MODE, PASS COMMON /NAMBLA/SUBNAM @@ -128,6 +131,8 @@ SUBROUTINE HEADER DATA L(9)/'CSSCAL'/ DATA L(10)/'ICAMAX'/ DATA L(11)/'CAXPBY'/ + DATA L(12)/'CROTG '/ + DATA L(13)/'CSROT '/ * .. Executable Statements .. SUBNAM = L(ICASE) @@ -137,6 +142,149 @@ SUBROUTINE HEADER 99999 FORMAT (/' Test of subprogram number',I3,12X,A6) * * End of HEADER +* + END + SUBROUTINE CHECK0(SFAC) + IMPLICIT NONE +* .. Parameters .. + INTEGER NOUT + PARAMETER (NOUT=6) +* .. Scalar Arguments .. + REAL SFAC +* .. Scalars in Common .. + INTEGER ICASE, INCX, INCY, MODE, N + LOGICAL PASS +* .. Local Scalars .. + COMPLEX CS, CSA, CSB, CA0, CB0 + REAL SC, SGOT, RTMIN, RTMAX, SONE, SSIZ, + + SZERO + INTEGER K +* .. Local Arrays .. + COMPLEX CA1(8), CB1(8), CSTRUE(8), CATRUE(8), + + CGOT(1), CWANT(1), CSIZE(1) + REAL SWANT(1) + REAL CCTRUE(8) +* .. External Subroutines .. + EXTERNAL CROTG, CTEST, STEST1 +* .. Intrinsic Functions .. + INTRINSIC ABS, AIMAG, HUGE, REAL, SQRT, TINY +* .. Common blocks .. + COMMON /COMBLA/ICASE, N, INCX, INCY, MODE, PASS +* .. Data statements .. +* +* CROTG is exercised on the branches its argument pair selects: +* B zero, A zero with B real, imaginary or general, and the ordinary +* path with real, imaginary and negative A. Each expected value +* satisfies C*A + S*B = R and C**2 + ABS(S)**2 = 1 exactly. +* + DATA CA1/(1.0E0,0.0E0), (0.0E0,0.0E0), + + (0.0E0,0.0E0), (0.0E0,0.0E0), + + (3.0E0,0.0E0), (0.0E0,3.0E0), + + (3.0E0,4.0E0), (-3.0E0,0.0E0)/ + DATA CB1/(0.0E0,0.0E0), (3.0E0,0.0E0), + + (0.0E0,4.0E0), (3.0E0,4.0E0), + + (4.0E0,0.0E0), (0.0E0,4.0E0), + + (0.0E0,0.0E0), (4.0E0,0.0E0)/ + DATA CCTRUE/1.0E0, 0.0E0, 0.0E0, 0.0E0, + + 0.6E0, 0.6E0, 1.0E0, 0.6E0/ + DATA CSTRUE/(0.0E0,0.0E0), (1.0E0,0.0E0), + + (0.0E0,-1.0E0), (0.6E0,-0.8E0), + + (0.8E0,0.0E0), (0.8E0,0.0E0), + + (0.0E0,0.0E0), (-0.8E0,0.0E0)/ + DATA CATRUE/(1.0E0,0.0E0), (3.0E0,0.0E0), + + (4.0E0,0.0E0), (5.0E0,0.0E0), + + (5.0E0,0.0E0), (0.0E0,5.0E0), + + (3.0E0,4.0E0), (-5.0E0,0.0E0)/ +* .. Executable Statements .. +* + DO 20 K = 1, 8 +* .. Set N=K for identification in output if any .. + N = K + CSA = CA1(K) + CSB = CB1(K) + CALL CROTG(CSA,CSB,SC,CS) + CGOT(1) = CSA + CWANT(1) = CATRUE(K) + CALL CTEST(1,CGOT,CWANT,CWANT,SFAC) + CGOT(1) = CS + CWANT(1) = CSTRUE(K) + CALL CTEST(1,CGOT,CWANT,CWANT,SFAC) + CALL STEST1(SC,CCTRUE(K),CCTRUE(K),SFAC) + 20 CONTINUE +* +* The scaled arms of CROTG. Their results are awkward to tabulate, +* so check the identities that define the rotation instead: +* C*A + S*B = R, and C**2 + ABS(S)**2 = 1. RTMIN and RTMAX bracket +* the window inside which the routine uses its unscaled algorithm. +* + SONE = 1E0 + SZERO = 0E0 + RTMIN = SQRT(TINY(SONE)) + RTMAX = SQRT(HUGE(SONE)) + DO 40 K = 1, 8 + N = 8 + K + IF (K.EQ.1) THEN +* A zero and B below the window, with both parts nonzero so +* the general arm is taken rather than the real or imaginary +* shortcut. + CA0 = CMPLX(SZERO,SZERO) + CB0 = CMPLX(RTMIN*0.25E0,RTMIN*0.25E0) + ELSE IF (K.EQ.2) THEN +* Both nonzero, B far above the window. + CA0 = CMPLX(SONE,SZERO) + CB0 = CMPLX(RTMAX,SZERO) + ELSE IF (K.EQ.3) THEN +* Both inside the window, but A so small against B that the +* routine switches to its F2 < H2*SAFMIN formulation. + CA0 = CMPLX(RTMIN*2E0,SZERO) + CB0 = CMPLX(RTMAX*0.25E0,SZERO) + ELSE IF (K.EQ.4) THEN +* Both above the window and of a size, so the routine scales +* A and B by the same factor rather than one of its own. + CA0 = CMPLX(RTMAX,SZERO) + CB0 = CMPLX(RTMAX,SZERO) + ELSE IF (K.EQ.5) THEN +* The same, below the window. + CA0 = CMPLX(RTMIN*0.25E0,SZERO) + CB0 = CMPLX(RTMIN*0.25E0,SZERO) + ELSE IF (K.EQ.6) THEN +* B above the window with both parts set, so the routine +* scales, and A small enough against it that the scaled +* algorithm takes its F2 < H2*SAFMIN arm as well. + CA0 = CMPLX(SONE,SZERO) + CB0 = CMPLX(RTMAX*0.5E0,RTMAX*0.5E0) + ELSE IF (K.EQ.7) THEN +* The same, with A large enough that SQRT(F2*H2) is no longer +* known to stay in range, so S is formed from R/H2 instead. + CA0 = CMPLX(1.5E0,SZERO) + CB0 = CMPLX(RTMAX*0.5E0,RTMAX*0.5E0) + ELSE +* Both inside the window, A just above its lower edge, where +* the unscaled algorithm forms S from R/H2 for that reason. + CA0 = CMPLX(RTMIN*2E0,SZERO) + CB0 = CMPLX(SONE,SONE) + END IF + CSA = CA0 + CSB = CB0 + CALL CROTG(CSA,CSB,SC,CS) +* CSA holds R on return. C*A and S*B are each bounded by the +* larger of ABS(A) and ABS(B), so their sum is the scale the +* residual of the identity has to be measured against. Passing +* CWANT as the size would ask the imaginary part for an exact +* zero, which it cannot give: for operands near either end of +* the range it is a cancellation of two quantities of that size. + SSIZ = ABS(CA0) + ABS(CB0) + CGOT(1) = SC*CA0 + CS*CB0 + CWANT(1) = CSA + CSIZE(1) = CMPLX(SSIZ,SSIZ) + CALL CTEST(1,CGOT,CWANT,CSIZE,SFAC) + SGOT = SC*SC + REAL(CS)*REAL(CS) + AIMAG(CS)*AIMAG(CS) + SWANT(1) = SONE + CALL STEST1(SGOT,SONE,SWANT,SFAC) + 40 CONTINUE + RETURN +* +* End of CHECK0 * END SUBROUTINE CHECK1(SFAC) @@ -380,17 +528,19 @@ SUBROUTINE CHECK2(SFAC) COMPLEX CA, CB INTEGER I, J, KI, KN, KSIZE, LENX, LENY, LINCX, LINCY, + MX, MY + REAL SC, SS * .. 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), CTY0(1), CX(7), CX0(1), CX1(7), - + CY(7), CY0(1), CY1(7), CT11(7,4,4) + + CY(7), CY0(1), CY1(7), CT11(7,4,4), + + CTX(7), CTY(7) INTEGER INCXS(4), INCYS(4), LENS(4,2), NS(4) * .. External Functions .. COMPLEX CDOTC, CDOTU EXTERNAL CDOTC, CDOTU * .. External Subroutines .. - EXTERNAL CAXPY, CAXPBY, CCOPY, CSWAP, CTEST + EXTERNAL CAXPY, CAXPBY, CCOPY, CSWAP, CTEST, CSROT * .. Intrinsic Functions .. INTRINSIC ABS, MIN * .. Common blocks .. @@ -615,6 +765,65 @@ SUBROUTINE CHECK2(SFAC) + (-0.2E0,-1.27E0)/ * .. Executable Statements .. +* +* CSROT applies a real plane rotation to a pair of complex vectors. +* The expected values are formed here with explicit indexing, so the +* stride arithmetic inside the routine is what is under test; the +* negative stride exercises its IX = (-N+1)*INCX + 1 start. +* + IF (ICASE.EQ.13) THEN + SC = 0.6E0 + SS = 0.8E0 + DO 100 KI = 1, 4 + IF (KI.EQ.1) THEN + LINCX = 1 + LINCY = 1 + ELSE IF (KI.EQ.2) THEN + LINCX = 2 + LINCY = 1 + ELSE IF (KI.EQ.3) THEN + LINCX = -2 + LINCY = 1 + ELSE + LINCX = -1 + LINCY = -2 + END IF + N = 3 + DO 80 I = 1, 7 + CX(I) = CX1(I) + CY(I) = CY1(I) + CTX(I) = CX1(I) + CTY(I) = CY1(I) + 80 CONTINUE + MX = 1 + MY = 1 + IF (LINCX.LT.0) MX = (-N+1)*LINCX + 1 + IF (LINCY.LT.0) MY = (-N+1)*LINCY + 1 + DO 90 I = 1, N + CTX(MX) = SC*CX1(MX) + SS*CY1(MY) + CTY(MY) = SC*CY1(MY) - SS*CX1(MX) + MX = MX + LINCX + MY = MY + LINCY + 90 CONTINUE + CALL CSROT(N,CX,LINCX,CY,LINCY,SC,SS) + CALL CTEST(7,CX,CTX,CTX,SFAC) + CALL CTEST(7,CY,CTY,CTY,SFAC) + 100 CONTINUE +* +* N = 0 must leave both vectors untouched. +* + N = 0 + DO 110 I = 1, 7 + CX(I) = CX1(I) + CY(I) = CY1(I) + CTX(I) = CX1(I) + CTY(I) = CY1(I) + 110 CONTINUE + CALL CSROT(N,CX,1,CY,1,SC,SS) + CALL CTEST(7,CX,CTX,CTX,SFAC) + CALL CTEST(7,CY,CTY,CTY,SFAC) + RETURN + END IF DO 60 KI = 1, 4 INCX = INCXS(KI) INCY = INCYS(KI) diff --git a/BLAS/TESTING/dblat1.f b/BLAS/TESTING/dblat1.f index 43590bd23..0b6a95289 100644 --- a/BLAS/TESTING/dblat1.f +++ b/BLAS/TESTING/dblat1.f @@ -161,25 +161,27 @@ SUBROUTINE CHECK0(SFAC) DOUBLE PRECISION SA, SB, SC, SS, D12 INTEGER I, K * .. Local Arrays .. - DOUBLE PRECISION DA1(8), DATRUE(8), DB1(8), DBTRUE(8), DC1(8), - $ DS1(8), DAB(4,9), DTEMP(9), DTRUE(9,9) + DOUBLE PRECISION DA1(9), DATRUE(9), DB1(9), DBTRUE(9), DC1(9), + $ DS1(9), DAB(4,13), DTEMP(9), DTRUE(9,13) * .. External Subroutines .. EXTERNAL DROTG, DROTMG, STEST, STEST1 +* .. Intrinsic Functions .. + INTRINSIC NEAREST * .. Common blocks .. COMMON /COMBLA/ICASE, N, INCX, INCY, PASS * .. Data statements .. DATA DA1/0.3D0, 0.4D0, -0.3D0, -0.4D0, -0.3D0, 0.0D0, - + 0.0D0, 1.0D0/ + + 0.0D0, 1.0D0, 0.0D0/ DATA DB1/0.4D0, 0.3D0, 0.4D0, 0.3D0, -0.4D0, 0.0D0, - + 1.0D0, 0.0D0/ + + 1.0D0, 0.0D0, 0.0D0/ DATA DC1/0.6D0, 0.8D0, -0.6D0, 0.8D0, 0.6D0, 1.0D0, - + 0.0D0, 1.0D0/ + + 0.0D0, 1.0D0, 0.0D0/ DATA DS1/0.8D0, 0.6D0, 0.8D0, -0.6D0, 0.8D0, 0.0D0, - + 1.0D0, 0.0D0/ + + 1.0D0, 0.0D0, 0.0D0/ DATA DATRUE/0.5D0, 0.5D0, 0.5D0, -0.5D0, -0.5D0, - + 0.0D0, 1.0D0, 1.0D0/ + + 0.0D0, 1.0D0, 1.0D0, 0.0D0/ DATA DBTRUE/0.0D0, 0.6D0, 0.0D0, -0.6D0, 0.0D0, - + 0.0D0, 1.0D0, 0.0D0/ + + 0.0D0, 1.0D0, 0.0D0, 0.0D0/ * INPUT FOR MODIFIED GIVENS DATA DAB/ .1D0,.3D0,1.2D0,.2D0, A .7D0, .2D0, .6D0, 4.2D0, @@ -189,7 +191,11 @@ SUBROUTINE CHECK0(SFAC) E 4.D10, 2.D-2, 1.D-5, 10.D0, F 2.D-10, 4.D-2, 1.D5, 10.D0, G 2.D10, 4.D-2, 1.D-5, 10.D0, - H 4.D-9, 2.D-9, 2.D0, 1.D0/ + H 4.D-9, 2.D-9, 2.D0, 1.D0, + I -1.D0, 1.D0, 1.D0, 1.D0, + J 0.D0,0.D0,0.D0,0.D0, + K 0.D0,0.D0,0.D0,0.D0, + L 0.D0,0.D0,0.D0,0.D0/ * TRUE RESULTS FOR MODIFIED GIVENS DATA DTRUE/0.D0,0.D0, 1.3D0, .2D0, 0.D0,0.D0,0.D0, .5D0, 0.D0, A 0.D0,0.D0, 4.5D0, 4.2D0, 1.D0, .5D0, 0.D0,0.D0,0.D0, @@ -202,7 +208,15 @@ SUBROUTINE CHECK0(SFAC) H 0.D0,0.D0,15.D0,10.D0,-1.D0, 5.D-5, 0.D0,1.D0,0.D0, I 0.D0,0.D0, 15.D0, 10.D0, -1. D0, 5.D5, -4096.D0, J 1.D0, 4096.D-6, - K 0.D0,0.D0, 7.D0, 4.D0, 0.D0,0.D0, -.5D0, -.25D0, 0.D0/ + K 0.D0,0.D0, 7.D0, 4.D0, 0.D0,0.D0, -.5D0, -.25D0, 0.D0, + Z 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + Z 0.D0,0.D0, + Z 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + Z 0.D0,0.D0, + Z 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + Z 0.D0,0.D0, + Z 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + Z 0.D0,0.D0/ * 4096 = 2 ** 12 DATA D12 /4096.D0/ DTRUE(1,1) = 12.D0 / 130.D0 @@ -233,6 +247,81 @@ SUBROUTINE CHECK0(SFAC) DTRUE(7,9) = -1.2207031250000000D-04 DTRUE(8,9) = 6.1035156250000000D-05 DTRUE(9,9) = 2.4414062500000000D-04 +* +* The tenth case has DD1 < 0, which ?ROTMG rejects by zeroing H, +* D and DX1 and returning DFLAG = -1. DY1 is left alone. +* + DTRUE(1,10) = 0.D0 + DTRUE(2,10) = 0.D0 + DTRUE(3,10) = 0.D0 + DTRUE(4,10) = 1.D0 + DTRUE(5,10) = -1.D0 + DTRUE(6,10) = 0.D0 + DTRUE(7,10) = 0.D0 + DTRUE(8,10) = 0.D0 + DTRUE(9,10) = 0.D0 +* +* Two more cases enter ?ROTMG's scale check holding a flag that no +* existing input pairs with it: the eleventh rescales D1 having set +* DFLAG = 1, the twelfth rescales D2 having set DFLAG = 0. Each of +* those arms rewrites H into the DFLAG = -1 form before scaling. +* The inputs are powers of two, so the eleventh case is exact +* throughout; the twelfth divides by SU = 1 + 2**(-10) = 1025/1024, +* which leaves its D1 and D2 as the quotients below. +* + DAB(1,11) = 2.0D0**(-40) + DAB(2,11) = 2.0D0**(-40) + DAB(3,11) = 1.D0 + DAB(4,11) = 1.D0 + DTRUE(1,11) = 2.0D0**(-17) + DTRUE(2,11) = 2.0D0**(-17) + DTRUE(3,11) = 2.0D0**(-11) + DTRUE(4,11) = 1.D0 + DTRUE(5,11) = -1.D0 + DTRUE(6,11) = 2.0D0**(-12) + DTRUE(7,11) = -2.0D0**(-12) + DTRUE(8,11) = 2.0D0**(-12) + DTRUE(9,11) = 2.0D0**(-12) +* + DAB(1,12) = 2.0D0**(-20) + DAB(2,12) = 2.0D0**(-30) + DAB(3,12) = 1.D0 + DAB(4,12) = 1.D0 + DTRUE(1,12) = 1.D0 / 1049600.D0 + DTRUE(2,12) = 16.D0 / 1025.D0 + DTRUE(3,12) = 1025.D0 / 1024.D0 + DTRUE(4,12) = 1.D0 + DTRUE(5,12) = -1.D0 + DTRUE(6,12) = 1.D0 + DTRUE(7,12) = -2.0D0**(-12) + DTRUE(8,12) = 2.0D0**(-10) + DTRUE(9,12) = 2.0D0**(-12) +* +* A thirteenth case reaches the block ?ROTMG keeps for safety +* against rounding. ABS(SQ1) > ABS(SQ2) holds, yet the quotients +* SP2/SP1 and SY1/SX1 each round up, so that their exact product is +* at least one and SU = 1 - SH12*SH21 is not positive. The block +* then zeroes H, D and X1; without it D1/SU would be infinite and +* the scale check would never terminate. Because the product is at +* least one exactly, the outcome does not depend on whether the +* compiler fuses it. The ulp alignment behind this is specific to +* the precision, so the constants differ between the real testers, +* and NEAREST steps D2 one ulp toward zero exactly in any arithmetic. +* + DAB(1,13) = 1.D0 + DAB(3,13) = 5.D0 + DAB(4,13) = 27.D0 + DAB(2,13) = NEAREST(-(DAB(1,13)*DAB(3,13)*DAB(3,13) + + /(DAB(4,13)*DAB(4,13))), 1.D0) + DTRUE(1,13) = 0.D0 + DTRUE(2,13) = 0.D0 + DTRUE(3,13) = 0.D0 + DTRUE(4,13) = DAB(4,13) + DTRUE(5,13) = -1.D0 + DTRUE(6,13) = 0.D0 + DTRUE(7,13) = 0.D0 + DTRUE(8,13) = 0.D0 + DTRUE(9,13) = 0.D0 * .. Executable Statements .. * * Compute true values which cannot be prestored @@ -242,12 +331,24 @@ SUBROUTINE CHECK0(SFAC) DBTRUE(3) = -1.0D0/0.6D0 DBTRUE(5) = 1.0D0/0.6D0 * - DO 20 K = 1, 9 +* A ninth ?ROTG case, built here because it needs run-time +* arithmetic: A is so small against B that C = A/R underflows to +* zero, which is the only way to reach the Z = ONE branch. Both +* values are powers of two, so the expected results are exact. +* + DA1(9) = 2.0D0**(-1000) + DB1(9) = 2.0D0**100 + DATRUE(9) = 2.0D0**100 + DBTRUE(9) = 1.0D0 + DC1(9) = 0.0D0 + DS1(9) = 1.0D0 +* + DO 20 K = 1, 13 * .. Set N=K for identification in output if any .. N = K IF (ICASE.EQ.3) THEN * .. DROTG .. - IF (K.GT.8) GO TO 40 + IF (K.GT.9) GO TO 40 SA = DA1(K) SB = DB1(K) CALL DROTG(SA,SB,SC,SS) diff --git a/BLAS/TESTING/sblat1.f b/BLAS/TESTING/sblat1.f index 4ea895abb..65e250ca7 100644 --- a/BLAS/TESTING/sblat1.f +++ b/BLAS/TESTING/sblat1.f @@ -160,25 +160,27 @@ SUBROUTINE CHECK0(SFAC) REAL D12, SA, SB, SC, SS INTEGER I, K * .. Local Arrays .. - REAL DA1(8), DATRUE(8), DB1(8), DBTRUE(8), DC1(8), - + DS1(8), DAB(4,9), DTEMP(9), DTRUE(9,9) + REAL DA1(9), DATRUE(9), DB1(9), DBTRUE(9), DC1(9), + + DS1(9), DAB(4,13), DTEMP(9), DTRUE(9,13) * .. External Subroutines .. EXTERNAL SROTG, SROTMG, STEST, STEST1 +* .. Intrinsic Functions .. + INTRINSIC NEAREST * .. Common blocks .. COMMON /COMBLA/ICASE, N, INCX, INCY, PASS * .. Data statements .. DATA DA1/0.3E0, 0.4E0, -0.3E0, -0.4E0, -0.3E0, 0.0E0, - + 0.0E0, 1.0E0/ + + 0.0E0, 1.0E0, 0.0E0/ DATA DB1/0.4E0, 0.3E0, 0.4E0, 0.3E0, -0.4E0, 0.0E0, - + 1.0E0, 0.0E0/ + + 1.0E0, 0.0E0, 0.0E0/ DATA DC1/0.6E0, 0.8E0, -0.6E0, 0.8E0, 0.6E0, 1.0E0, - + 0.0E0, 1.0E0/ + + 0.0E0, 1.0E0, 0.0E0/ DATA DS1/0.8E0, 0.6E0, 0.8E0, -0.6E0, 0.8E0, 0.0E0, - + 1.0E0, 0.0E0/ + + 1.0E0, 0.0E0, 0.0E0/ DATA DATRUE/0.5E0, 0.5E0, 0.5E0, -0.5E0, -0.5E0, - + 0.0E0, 1.0E0, 1.0E0/ + + 0.0E0, 1.0E0, 1.0E0, 0.0E0/ DATA DBTRUE/0.0E0, 0.6E0, 0.0E0, -0.6E0, 0.0E0, - + 0.0E0, 1.0E0, 0.0E0/ + + 0.0E0, 1.0E0, 0.0E0, 0.0E0/ * INPUT FOR MODIFIED GIVENS DATA DAB/ .1E0,.3E0,1.2E0,.2E0, A .7E0, .2E0, .6E0, 4.2E0, @@ -188,7 +190,11 @@ SUBROUTINE CHECK0(SFAC) E 4.E10, 2.E-2, 1.E-5, 10.E0, F 2.E-10, 4.E-2, 1.E5, 10.E0, G 2.E10, 4.E-2, 1.E-5, 10.E0, - H 4.E-9, 2.E-9, 2.E0, 1.E0/ + H 4.E-9, 2.E-9, 2.E0, 1.E0, + I -1.E0, 1.E0, 1.E0, 1.E0, + J 0.E0,0.E0,0.E0,0.E0, + K 0.E0,0.E0,0.E0,0.E0, + L 0.E0,0.E0,0.E0,0.E0/ * TRUE RESULTS FOR MODIFIED GIVENS DATA DTRUE/0.E0,0.E0, 1.3E0, .2E0, 0.E0,0.E0,0.E0, .5E0, 0.E0, A 0.E0,0.E0, 4.5E0, 4.2E0, 1.E0, .5E0, 0.E0,0.E0,0.E0, @@ -201,7 +207,15 @@ SUBROUTINE CHECK0(SFAC) H 0.E0,0.E0,15.E0,10.E0,-1.E0, 5.E-5, 0.E0,1.E0,0.E0, I 0.E0,0.E0, 15.E0, 10.E0, -1. E0, 5.E5, -4096.E0, J 1.E0, 4096.E-6, - K 0.E0,0.E0, 7.E0, 4.E0, 0.E0,0.E0, -.5E0, -.25E0, 0.E0/ + K 0.E0,0.E0, 7.E0, 4.E0, 0.E0,0.E0, -.5E0, -.25E0, 0.E0, + Z 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + Z 0.E0,0.E0, + Z 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + Z 0.E0,0.E0, + Z 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + Z 0.E0,0.E0, + Z 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + Z 0.E0,0.E0/ * 4096 = 2 ** 12 DATA D12 /4096.E0/ DTRUE(1,1) = 12.E0 / 130.E0 @@ -232,6 +246,81 @@ SUBROUTINE CHECK0(SFAC) DTRUE(7,9) = -1.2207031250000000E-04 DTRUE(8,9) = 6.1035156250000000E-05 DTRUE(9,9) = 2.4414062500000000E-04 +* +* The tenth case has DD1 < 0, which ?ROTMG rejects by zeroing H, +* D and DX1 and returning DFLAG = -1. DY1 is left alone. +* + DTRUE(1,10) = 0.E0 + DTRUE(2,10) = 0.E0 + DTRUE(3,10) = 0.E0 + DTRUE(4,10) = 1.E0 + DTRUE(5,10) = -1.E0 + DTRUE(6,10) = 0.E0 + DTRUE(7,10) = 0.E0 + DTRUE(8,10) = 0.E0 + DTRUE(9,10) = 0.E0 +* +* Two more cases enter ?ROTMG's scale check holding a flag that no +* existing input pairs with it: the eleventh rescales D1 having set +* DFLAG = 1, the twelfth rescales D2 having set DFLAG = 0. Each of +* those arms rewrites H into the DFLAG = -1 form before scaling. +* The inputs are powers of two, so the eleventh case is exact +* throughout; the twelfth divides by SU = 1 + 2**(-10) = 1025/1024, +* which leaves its D1 and D2 as the quotients below. +* + DAB(1,11) = 2.0E0**(-40) + DAB(2,11) = 2.0E0**(-40) + DAB(3,11) = 1.E0 + DAB(4,11) = 1.E0 + DTRUE(1,11) = 2.0E0**(-17) + DTRUE(2,11) = 2.0E0**(-17) + DTRUE(3,11) = 2.0E0**(-11) + DTRUE(4,11) = 1.E0 + DTRUE(5,11) = -1.E0 + DTRUE(6,11) = 2.0E0**(-12) + DTRUE(7,11) = -2.0E0**(-12) + DTRUE(8,11) = 2.0E0**(-12) + DTRUE(9,11) = 2.0E0**(-12) +* + DAB(1,12) = 2.0E0**(-20) + DAB(2,12) = 2.0E0**(-30) + DAB(3,12) = 1.E0 + DAB(4,12) = 1.E0 + DTRUE(1,12) = 1.E0 / 1049600.E0 + DTRUE(2,12) = 16.E0 / 1025.E0 + DTRUE(3,12) = 1025.E0 / 1024.E0 + DTRUE(4,12) = 1.E0 + DTRUE(5,12) = -1.E0 + DTRUE(6,12) = 1.E0 + DTRUE(7,12) = -2.0E0**(-12) + DTRUE(8,12) = 2.0E0**(-10) + DTRUE(9,12) = 2.0E0**(-12) +* +* A thirteenth case reaches the block ?ROTMG keeps for safety +* against rounding. ABS(SQ1) > ABS(SQ2) holds, yet the quotients +* SP2/SP1 and SY1/SX1 each round up, so that their exact product is +* at least one and SU = 1 - SH12*SH21 is not positive. The block +* then zeroes H, D and X1; without it D1/SU would be infinite and +* the scale check would never terminate. Because the product is at +* least one exactly, the outcome does not depend on whether the +* compiler fuses it. The ulp alignment behind this is specific to +* the precision, so the constants differ between the real testers, +* and NEAREST steps D2 one ulp toward zero exactly in any arithmetic. +* + DAB(1,13) = 1.E0 + DAB(3,13) = 13.E0 + DAB(4,13) = 9.E0 + DAB(2,13) = NEAREST(-(DAB(1,13)*DAB(3,13)*DAB(3,13) + + /(DAB(4,13)*DAB(4,13))), 1.E0) + DTRUE(1,13) = 0.E0 + DTRUE(2,13) = 0.E0 + DTRUE(3,13) = 0.E0 + DTRUE(4,13) = DAB(4,13) + DTRUE(5,13) = -1.E0 + DTRUE(6,13) = 0.E0 + DTRUE(7,13) = 0.E0 + DTRUE(8,13) = 0.E0 + DTRUE(9,13) = 0.E0 * .. Executable Statements .. * * Compute true values which cannot be prestored @@ -241,12 +330,24 @@ SUBROUTINE CHECK0(SFAC) DBTRUE(3) = -1.0E0/0.6E0 DBTRUE(5) = 1.0E0/0.6E0 * - DO 20 K = 1, 9 +* A ninth ?ROTG case, built here because it needs run-time +* arithmetic: A is so small against B that C = A/R underflows to +* zero, which is the only way to reach the Z = ONE branch. Both +* values are powers of two, so the expected results are exact. +* + DA1(9) = 2.0E0**(-100) + DB1(9) = 2.0E0**60 + DATRUE(9) = 2.0E0**60 + DBTRUE(9) = 1.0E0 + DC1(9) = 0.0E0 + DS1(9) = 1.0E0 +* + DO 20 K = 1, 13 * .. Set N=K for identification in output if any .. N = K IF (ICASE.EQ.3) THEN * .. SROTG .. - IF (K.GT.8) GO TO 40 + IF (K.GT.9) GO TO 40 SA = DA1(K) SB = DB1(K) CALL SROTG(SA,SB,SC,SS) diff --git a/BLAS/TESTING/zblat1.f b/BLAS/TESTING/zblat1.f index 7742590d7..486bdeaef 100644 --- a/BLAS/TESTING/zblat1.f +++ b/BLAS/TESTING/zblat1.f @@ -65,7 +65,7 @@ PROGRAM ZBLAT1 * .. Executable Statements .. CALL CPU_TIME( S1 ) WRITE (NOUT,99999) - DO 20 IC = 1, 11 + DO 20 IC = 1, 13 ICASE = IC CALL HEADER * @@ -80,7 +80,10 @@ PROGRAM ZBLAT1 INCX = 9999 INCY = 9999 MODE = 9999 - IF (ICASE.LE.5 .OR. ICASE.EQ.11) THEN + IF (ICASE.EQ.12) THEN + CALL CHECK0(SFAC) + ELSE IF (ICASE.LE.5 .OR. ICASE.EQ.11 .OR. + + ICASE.EQ.13) THEN CALL CHECK2(SFAC) ELSE IF (ICASE.GE.6) THEN CALL CHECK1(SFAC) @@ -111,7 +114,7 @@ SUBROUTINE HEADER INTEGER ICASE, INCX, INCY, MODE, N LOGICAL PASS * .. Local Arrays .. - CHARACTER*6 L(11) + CHARACTER*6 L(13) * .. Common blocks .. COMMON /COMBLA/ICASE, N, INCX, INCY, MODE, PASS COMMON /NAMBLA/SUBNAM @@ -127,6 +130,8 @@ SUBROUTINE HEADER DATA L(9)/'ZDSCAL'/ DATA L(10)/'IZAMAX'/ DATA L(11)/'ZAXPBY'/ + DATA L(12)/'ZROTG '/ + DATA L(13)/'ZDROT '/ * .. Executable Statements .. SUBNAM = L(ICASE) @@ -136,6 +141,149 @@ SUBROUTINE HEADER 99999 FORMAT (/' Test of subprogram number',I3,12X,A6) * * End of HEADER +* + END + SUBROUTINE CHECK0(SFAC) + IMPLICIT NONE +* .. Parameters .. + INTEGER NOUT + PARAMETER (NOUT=6) +* .. Scalar Arguments .. + DOUBLE PRECISION SFAC +* .. Scalars in Common .. + INTEGER ICASE, INCX, INCY, MODE, N + LOGICAL PASS +* .. Local Scalars .. + COMPLEX*16 CS, CSA, CSB, CA0, CB0 + DOUBLE PRECISION SC, SGOT, RTMIN, RTMAX, + + SONE, SSIZ, SZERO + INTEGER K +* .. Local Arrays .. + COMPLEX*16 CA1(8), CB1(8), CSTRUE(8), CATRUE(8), + + CGOT(1), CWANT(1), CSIZE(1) + DOUBLE PRECISION SWANT(1) + DOUBLE PRECISION CCTRUE(8) +* .. External Subroutines .. + EXTERNAL ZROTG, CTEST, STEST1 +* .. Intrinsic Functions .. + INTRINSIC ABS, AIMAG, HUGE, REAL, SQRT, TINY +* .. Common blocks .. + COMMON /COMBLA/ICASE, N, INCX, INCY, MODE, PASS +* .. Data statements .. +* +* ZROTG is exercised on the branches its argument pair selects: +* B zero, A zero with B real, imaginary or general, and the ordinary +* path with real, imaginary and negative A. Each expected value +* satisfies C*A + S*B = R and C**2 + ABS(S)**2 = 1 exactly. +* + DATA CA1/(1.0D0,0.0D0), (0.0D0,0.0D0), + + (0.0D0,0.0D0), (0.0D0,0.0D0), + + (3.0D0,0.0D0), (0.0D0,3.0D0), + + (3.0D0,4.0D0), (-3.0D0,0.0D0)/ + DATA CB1/(0.0D0,0.0D0), (3.0D0,0.0D0), + + (0.0D0,4.0D0), (3.0D0,4.0D0), + + (4.0D0,0.0D0), (0.0D0,4.0D0), + + (0.0D0,0.0D0), (4.0D0,0.0D0)/ + DATA CCTRUE/1.0D0, 0.0D0, 0.0D0, 0.0D0, + + 0.6D0, 0.6D0, 1.0D0, 0.6D0/ + DATA CSTRUE/(0.0D0,0.0D0), (1.0D0,0.0D0), + + (0.0D0,-1.0D0), (0.6D0,-0.8D0), + + (0.8D0,0.0D0), (0.8D0,0.0D0), + + (0.0D0,0.0D0), (-0.8D0,0.0D0)/ + DATA CATRUE/(1.0D0,0.0D0), (3.0D0,0.0D0), + + (4.0D0,0.0D0), (5.0D0,0.0D0), + + (5.0D0,0.0D0), (0.0D0,5.0D0), + + (3.0D0,4.0D0), (-5.0D0,0.0D0)/ +* .. Executable Statements .. +* + DO 20 K = 1, 8 +* .. Set N=K for identification in output if any .. + N = K + CSA = CA1(K) + CSB = CB1(K) + CALL ZROTG(CSA,CSB,SC,CS) + CGOT(1) = CSA + CWANT(1) = CATRUE(K) + CALL CTEST(1,CGOT,CWANT,CWANT,SFAC) + CGOT(1) = CS + CWANT(1) = CSTRUE(K) + CALL CTEST(1,CGOT,CWANT,CWANT,SFAC) + CALL STEST1(SC,CCTRUE(K),CCTRUE(K),SFAC) + 20 CONTINUE +* +* The scaled arms of ZROTG. Their results are awkward to tabulate, +* so check the identities that define the rotation instead: +* C*A + S*B = R, and C**2 + ABS(S)**2 = 1. RTMIN and RTMAX bracket +* the window inside which the routine uses its unscaled algorithm. +* + SONE = 1D0 + SZERO = 0D0 + RTMIN = SQRT(TINY(SONE)) + RTMAX = SQRT(HUGE(SONE)) + DO 40 K = 1, 8 + N = 8 + K + IF (K.EQ.1) THEN +* A zero and B below the window, with both parts nonzero so +* the general arm is taken rather than the real or imaginary +* shortcut. + CA0 = DCMPLX(SZERO,SZERO) + CB0 = DCMPLX(RTMIN*0.25D0,RTMIN*0.25D0) + ELSE IF (K.EQ.2) THEN +* Both nonzero, B far above the window. + CA0 = DCMPLX(SONE,SZERO) + CB0 = DCMPLX(RTMAX,SZERO) + ELSE IF (K.EQ.3) THEN +* Both inside the window, but A so small against B that the +* routine switches to its F2 < H2*SAFMIN formulation. + CA0 = DCMPLX(RTMIN*2D0,SZERO) + CB0 = DCMPLX(RTMAX*0.25D0,SZERO) + ELSE IF (K.EQ.4) THEN +* Both above the window and of a size, so the routine scales +* A and B by the same factor rather than one of its own. + CA0 = DCMPLX(RTMAX,SZERO) + CB0 = DCMPLX(RTMAX,SZERO) + ELSE IF (K.EQ.5) THEN +* The same, below the window. + CA0 = DCMPLX(RTMIN*0.25D0,SZERO) + CB0 = DCMPLX(RTMIN*0.25D0,SZERO) + ELSE IF (K.EQ.6) THEN +* B above the window with both parts set, so the routine +* scales, and A small enough against it that the scaled +* algorithm takes its F2 < H2*SAFMIN arm as well. + CA0 = DCMPLX(SONE,SZERO) + CB0 = DCMPLX(RTMAX*0.5D0,RTMAX*0.5D0) + ELSE IF (K.EQ.7) THEN +* The same, with A large enough that SQRT(F2*H2) is no longer +* known to stay in range, so S is formed from R/H2 instead. + CA0 = DCMPLX(1.5D0,SZERO) + CB0 = DCMPLX(RTMAX*0.5D0,RTMAX*0.5D0) + ELSE +* Both inside the window, A just above its lower edge, where +* the unscaled algorithm forms S from R/H2 for that reason. + CA0 = DCMPLX(RTMIN*2D0,SZERO) + CB0 = DCMPLX(SONE,SONE) + END IF + CSA = CA0 + CSB = CB0 + CALL ZROTG(CSA,CSB,SC,CS) +* CSA holds R on return. C*A and S*B are each bounded by the +* larger of ABS(A) and ABS(B), so their sum is the scale the +* residual of the identity has to be measured against. Passing +* CWANT as the size would ask the imaginary part for an exact +* zero, which it cannot give: for operands near either end of +* the range it is a cancellation of two quantities of that size. + SSIZ = ABS(CA0) + ABS(CB0) + CGOT(1) = SC*CA0 + CS*CB0 + CWANT(1) = CSA + CSIZE(1) = DCMPLX(SSIZ,SSIZ) + CALL CTEST(1,CGOT,CWANT,CSIZE,SFAC) + SGOT = SC*SC + REAL(CS)*REAL(CS) + AIMAG(CS)*AIMAG(CS) + SWANT(1) = SONE + CALL STEST1(SGOT,SONE,SWANT,SFAC) + 40 CONTINUE + RETURN +* +* End of CHECK0 * END SUBROUTINE CHECK1(SFAC) @@ -379,17 +527,19 @@ SUBROUTINE CHECK2(SFAC) COMPLEX*16 CA, CB INTEGER I, J, KI, KN, KSIZE, LENX, LENY, LINCX, LINCY, + MX, MY + DOUBLE PRECISION SC, SS * .. 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), CTY0(1), CX(7), CX0(1), CX1(7), - + CY(7), CY0(1), CY1(7), CT11(7,4,4) + + CY(7), CY0(1), CY1(7), CT11(7,4,4), + + CTX(7), CTY(7) INTEGER INCXS(4), INCYS(4), LENS(4,2), NS(4) * .. External Functions .. COMPLEX*16 ZDOTC, ZDOTU EXTERNAL ZDOTC, ZDOTU * .. External Subroutines .. - EXTERNAL ZAXPY, ZAXPBY, ZCOPY, ZSWAP, CTEST + EXTERNAL ZAXPY, ZAXPBY, ZCOPY, ZSWAP, CTEST, ZDROT * .. Intrinsic Functions .. INTRINSIC ABS, MIN * .. Common blocks .. @@ -615,6 +765,65 @@ SUBROUTINE CHECK2(SFAC) * .. Executable Statements .. +* +* ZDROT applies a real plane rotation to a pair of complex vectors. +* The expected values are formed here with explicit indexing, so the +* stride arithmetic inside the routine is what is under test; the +* negative stride exercises its IX = (-N+1)*INCX + 1 start. +* + IF (ICASE.EQ.13) THEN + SC = 0.6D0 + SS = 0.8D0 + DO 100 KI = 1, 4 + IF (KI.EQ.1) THEN + LINCX = 1 + LINCY = 1 + ELSE IF (KI.EQ.2) THEN + LINCX = 2 + LINCY = 1 + ELSE IF (KI.EQ.3) THEN + LINCX = -2 + LINCY = 1 + ELSE + LINCX = -1 + LINCY = -2 + END IF + N = 3 + DO 80 I = 1, 7 + CX(I) = CX1(I) + CY(I) = CY1(I) + CTX(I) = CX1(I) + CTY(I) = CY1(I) + 80 CONTINUE + MX = 1 + MY = 1 + IF (LINCX.LT.0) MX = (-N+1)*LINCX + 1 + IF (LINCY.LT.0) MY = (-N+1)*LINCY + 1 + DO 90 I = 1, N + CTX(MX) = SC*CX1(MX) + SS*CY1(MY) + CTY(MY) = SC*CY1(MY) - SS*CX1(MX) + MX = MX + LINCX + MY = MY + LINCY + 90 CONTINUE + CALL ZDROT(N,CX,LINCX,CY,LINCY,SC,SS) + CALL CTEST(7,CX,CTX,CTX,SFAC) + CALL CTEST(7,CY,CTY,CTY,SFAC) + 100 CONTINUE +* +* N = 0 must leave both vectors untouched. +* + N = 0 + DO 110 I = 1, 7 + CX(I) = CX1(I) + CY(I) = CY1(I) + CTX(I) = CX1(I) + CTY(I) = CY1(I) + 110 CONTINUE + CALL ZDROT(N,CX,1,CY,1,SC,SS) + CALL CTEST(7,CX,CTX,CTX,SFAC) + CALL CTEST(7,CY,CTY,CTY,SFAC) + RETURN + END IF DO 60 KI = 1, 4 INCX = INCXS(KI) INCY = INCYS(KI)