Edited the code of xGBMV, xGEMV and xGEMM by removing the test for an entry of

B (or X) being zero to skip a loop so as to enable better NaN (or Inf)
propagation.

Taking DGEMM as an example, if you use notation:
	C(I,J) = C(I,J) + alpha * A(I,L) * B(L,J),
the reference BLAS DGEMM is using the JLI version of matrix matrix multiply and
is checking, for each J (from 1 to N) and for each L (from 1 to K), whether
B(L,J) is zero (or not) to save (or not) the 2M following operations.

(See the "IF (B(L,J).NE.ZERO) THEN" in the code below.)

The snippets of code is as follows

              DO 90 J = 1,N
                  DO 80 L = 1,K
                      IF (B(L,J).NE.ZERO) THEN
                          TEMP = ALPHA*B(L,J)
                          DO 70 I = 1,M
                              C(I,J) = C(I,J) + TEMP*A(I,L)
   70                     CONTINUE
                      END IF
   80             CONTINUE
   90         CONTINUE

This induces some non NaN-propagation in a pretty ad-hoc way. For better NaN
propagation, this patch removes the above IF statement.

The snippet of code now becomes

              DO 90 J = 1,N
                  DO 80 L = 1,K
                      TEMP = ALPHA*B(L,J)
                      DO 70 I = 1,M
                          C(I,J) = C(I,J) + TEMP*A(I,L)
   70                 CONTINUE
   80             CONTINUE
   90         CONTINUE

This enables correct NaN propagation for this piece of code.

Rationale: BLAS does not correctly propagate all NaNs (and Infs). We still have
no NaN propagation where for example ALPHA=0, etc. The goal of this commit is
to have correct NaN propagation no matter what the entries of the input
matrices/vectors (A, B, C, X, etc.) are.  BLAS do not correctly propagate NaNs
and Infs based on some values of the scalars (ALPHA, BETA, etc.).

See below the email from Tom Callaway from RedHat, sent on July 9th to
lapack@cs.utk.edu.

Hello LAPACK people,

Martyn & Lejeczek (on CC) reported an issue to Fedora relating to R using our
system copy of BLAS (from LAPACK).

As noted in the R administration and Installation Manual, "R relies on ISO/IEC
60559 compliance of an external BLAS. This can be broken if for example the
code assumes that terms with a zero factor are always zero and do not need to
be computed - whereas x*0 can be NaN. This is checked in the test suite."

In the stock BLAS, DGBMV, DGEMM, and DGEMV fail this. R has been patching their
bundled BLAS to resolve this issue since 2010, but Fedora now uses the system
BLAS.

Attached is a patch (from upstream R) to fix this issue in the LAPACK BLAS.
Please consider applying it.

Thanks,

~tom

diff --git a/BLAS/SRC/cgbmv.f b/BLAS/SRC/cgbmv.f
index cd597f9..62001f2 100644 (file)
--- a/BLAS/SRC/cgbmv.f
+++ b/BLAS/SRC/cgbmv.f
@@ -319,26 +319,22 @@
           JX = KX
           IF (INCY.EQ.1) THEN
               DO 60 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      K = KUP1 - J
-                      DO 50 I = MAX(1,J-KU),MIN(M,J+KL)
-                          Y(I) = Y(I) + TEMP*A(K+I,J)
-   50                 CONTINUE
-                  END IF
+                  TEMP = ALPHA*X(JX)
+                  K = KUP1 - J
+                  DO 50 I = MAX(1,J-KU),MIN(M,J+KL)
+                      Y(I) = Y(I) + TEMP*A(K+I,J)
+   50             CONTINUE
                   JX = JX + INCX
    60         CONTINUE
           ELSE
               DO 80 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      IY = KY
-                      K = KUP1 - J
-                      DO 70 I = MAX(1,J-KU),MIN(M,J+KL)
-                          Y(IY) = Y(IY) + TEMP*A(K+I,J)
-                          IY = IY + INCY
-   70                 CONTINUE
-                  END IF
+                  TEMP = ALPHA*X(JX)
+                  IY = KY
+                  K = KUP1 - J
+                  DO 70 I = MAX(1,J-KU),MIN(M,J+KL)
+                      Y(IY) = Y(IY) + TEMP*A(K+I,J)
+                      IY = IY + INCY
+   70             CONTINUE
                   JX = JX + INCX
                   IF (J.GT.KU) KY = KY + INCY
    80         CONTINUE

diff --git a/BLAS/SRC/cgemm.f b/BLAS/SRC/cgemm.f
index ecfe12e..8d92e24 100644 (file)
--- a/BLAS/SRC/cgemm.f
+++ b/BLAS/SRC/cgemm.f
@@ -317,12 +317,10 @@
    60                 CONTINUE
                   END IF
                   DO 80 L = 1,K
-                      IF (B(L,J).NE.ZERO) THEN
-                          TEMP = ALPHA*B(L,J)
-                          DO 70 I = 1,M
-                              C(I,J) = C(I,J) + TEMP*A(I,L)
-   70                     CONTINUE
-                      END IF
+                      TEMP = ALPHA*B(L,J)
+                      DO 70 I = 1,M
+                          C(I,J) = C(I,J) + TEMP*A(I,L)
+   70                 CONTINUE
    80             CONTINUE
    90         CONTINUE
           ELSE IF (CONJA) THEN
@@ -376,17 +374,15 @@
   170                 CONTINUE
                   END IF
                   DO 190 L = 1,K
-                      IF (B(J,L).NE.ZERO) THEN
-                          TEMP = ALPHA*CONJG(B(J,L))
-                          DO 180 I = 1,M
-                              C(I,J) = C(I,J) + TEMP*A(I,L)
-  180                     CONTINUE
-                      END IF
+                      TEMP = ALPHA*CONJG(B(J,L))
+                      DO 180 I = 1,M
+                          C(I,J) = C(I,J) + TEMP*A(I,L)
+  180                 CONTINUE
   190             CONTINUE
   200         CONTINUE
           ELSE
 *
-*           Form  C := alpha*A*B**T          + beta*C
+*           Form  C := alpha*A*B**T + beta*C
 *
               DO 250 J = 1,N
                   IF (BETA.EQ.ZERO) THEN
@@ -399,12 +395,10 @@
   220                 CONTINUE
                   END IF
                   DO 240 L = 1,K
-                      IF (B(J,L).NE.ZERO) THEN
-                          TEMP = ALPHA*B(J,L)
-                          DO 230 I = 1,M
-                              C(I,J) = C(I,J) + TEMP*A(I,L)
-  230                     CONTINUE
-                      END IF
+                      TEMP = ALPHA*B(J,L)
+                      DO 230 I = 1,M
+                          C(I,J) = C(I,J) + TEMP*A(I,L)
+  230                 CONTINUE
   240             CONTINUE
   250         CONTINUE
           END IF

diff --git a/BLAS/SRC/cgemv.f b/BLAS/SRC/cgemv.f
index 507d19e..5447086 100644 (file)
--- a/BLAS/SRC/cgemv.f
+++ b/BLAS/SRC/cgemv.f
@@ -285,24 +285,20 @@
           JX = KX
           IF (INCY.EQ.1) THEN
               DO 60 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      DO 50 I = 1,M
-                          Y(I) = Y(I) + TEMP*A(I,J)
-   50                 CONTINUE
-                  END IF
+                  TEMP = ALPHA*X(JX)
+                  DO 50 I = 1,M
+                      Y(I) = Y(I) + TEMP*A(I,J)
+   50             CONTINUE
                   JX = JX + INCX
    60         CONTINUE
           ELSE
               DO 80 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      IY = KY
-                      DO 70 I = 1,M
-                          Y(IY) = Y(IY) + TEMP*A(I,J)
-                          IY = IY + INCY
-   70                 CONTINUE
-                  END IF
+                  TEMP = ALPHA*X(JX)
+                  IY = KY
+                  DO 70 I = 1,M
+                      Y(IY) = Y(IY) + TEMP*A(I,J)
+                      IY = IY + INCY
+   70             CONTINUE
                   JX = JX + INCX
    80         CONTINUE
           END IF

diff --git a/BLAS/SRC/dgbmv.f b/BLAS/SRC/dgbmv.f
index 4a608bd..b20516f 100644 (file)
--- a/BLAS/SRC/dgbmv.f
+++ b/BLAS/SRC/dgbmv.f
@@ -312,26 +312,22 @@
           JX = KX
           IF (INCY.EQ.1) THEN
               DO 60 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      K = KUP1 - J
-                      DO 50 I = MAX(1,J-KU),MIN(M,J+KL)
-                          Y(I) = Y(I) + TEMP*A(K+I,J)
-   50                 CONTINUE
-                  END IF
+                  TEMP = ALPHA*X(JX)
+                  K = KUP1 - J
+                  DO 50 I = MAX(1,J-KU),MIN(M,J+KL)
+                      Y(I) = Y(I) + TEMP*A(K+I,J)
+   50             CONTINUE
                   JX = JX + INCX
    60         CONTINUE
           ELSE
               DO 80 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      IY = KY
-                      K = KUP1 - J
-                      DO 70 I = MAX(1,J-KU),MIN(M,J+KL)
-                          Y(IY) = Y(IY) + TEMP*A(K+I,J)
-                          IY = IY + INCY
-   70                 CONTINUE
-                  END IF
+                  TEMP = ALPHA*X(JX)
+                  IY = KY
+                  K = KUP1 - J
+                  DO 70 I = MAX(1,J-KU),MIN(M,J+KL)
+                      Y(IY) = Y(IY) + TEMP*A(K+I,J)
+                      IY = IY + INCY
+   70             CONTINUE
                   JX = JX + INCX
                   IF (J.GT.KU) KY = KY + INCY
    80         CONTINUE

diff --git a/BLAS/SRC/dgemm.f b/BLAS/SRC/dgemm.f
index 45d001b..2136740 100644 (file)
--- a/BLAS/SRC/dgemm.f
+++ b/BLAS/SRC/dgemm.f
@@ -311,12 +311,10 @@
    60                 CONTINUE
                   END IF
                   DO 80 L = 1,K
-                      IF (B(L,J).NE.ZERO) THEN
-                          TEMP = ALPHA*B(L,J)
-                          DO 70 I = 1,M
-                              C(I,J) = C(I,J) + TEMP*A(I,L)
-   70                     CONTINUE
-                      END IF
+                      TEMP = ALPHA*B(L,J)
+                      DO 70 I = 1,M
+                          C(I,J) = C(I,J) + TEMP*A(I,L)
+   70                 CONTINUE
    80             CONTINUE
    90         CONTINUE
           ELSE
@@ -353,12 +351,10 @@
   140                 CONTINUE
                   END IF
                   DO 160 L = 1,K
-                      IF (B(J,L).NE.ZERO) THEN
-                          TEMP = ALPHA*B(J,L)
-                          DO 150 I = 1,M
-                              C(I,J) = C(I,J) + TEMP*A(I,L)
-  150                     CONTINUE
-                      END IF
+                      TEMP = ALPHA*B(J,L)
+                      DO 150 I = 1,M
+                          C(I,J) = C(I,J) + TEMP*A(I,L)
+  150                 CONTINUE
   160             CONTINUE
   170         CONTINUE
           ELSE

diff --git a/BLAS/SRC/dgemv.f b/BLAS/SRC/dgemv.f
index 675257f..6c29338 100644 (file)
--- a/BLAS/SRC/dgemv.f
+++ b/BLAS/SRC/dgemv.f
@@ -278,24 +278,20 @@
           JX = KX
           IF (INCY.EQ.1) THEN
               DO 60 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      DO 50 I = 1,M
-                          Y(I) = Y(I) + TEMP*A(I,J)
-   50                 CONTINUE
-                  END IF
+                  TEMP = ALPHA*X(JX)
+                  DO 50 I = 1,M
+                      Y(I) = Y(I) + TEMP*A(I,J)
+   50             CONTINUE
                   JX = JX + INCX
    60         CONTINUE
           ELSE
               DO 80 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      IY = KY
-                      DO 70 I = 1,M
-                          Y(IY) = Y(IY) + TEMP*A(I,J)
-                          IY = IY + INCY
-   70                 CONTINUE
-                  END IF
+                  TEMP = ALPHA*X(JX)
+                  IY = KY
+                  DO 70 I = 1,M
+                      Y(IY) = Y(IY) + TEMP*A(I,J)
+                      IY = IY + INCY
+   70             CONTINUE
                   JX = JX + INCX
    80         CONTINUE
           END IF

diff --git a/BLAS/SRC/sgbmv.f b/BLAS/SRC/sgbmv.f
index 797ac7f..26e5116 100644 (file)
--- a/BLAS/SRC/sgbmv.f
+++ b/BLAS/SRC/sgbmv.f
@@ -312,26 +312,22 @@
           JX = KX
           IF (INCY.EQ.1) THEN
               DO 60 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      K = KUP1 - J
-                      DO 50 I = MAX(1,J-KU),MIN(M,J+KL)
-                          Y(I) = Y(I) + TEMP*A(K+I,J)
-   50                 CONTINUE
-                  END IF
+                  TEMP = ALPHA*X(JX)
+                  K = KUP1 - J
+                  DO 50 I = MAX(1,J-KU),MIN(M,J+KL)
+                      Y(I) = Y(I) + TEMP*A(K+I,J)
+   50             CONTINUE
                   JX = JX + INCX
    60         CONTINUE
           ELSE
               DO 80 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      IY = KY
-                      K = KUP1 - J
-                      DO 70 I = MAX(1,J-KU),MIN(M,J+KL)
-                          Y(IY) = Y(IY) + TEMP*A(K+I,J)
-                          IY = IY + INCY
-   70                 CONTINUE
-                  END IF
+                  TEMP = ALPHA*X(JX)
+                  IY = KY
+                  K = KUP1 - J
+                  DO 70 I = MAX(1,J-KU),MIN(M,J+KL)
+                      Y(IY) = Y(IY) + TEMP*A(K+I,J)
+                      IY = IY + INCY
+   70             CONTINUE
                   JX = JX + INCX
                   IF (J.GT.KU) KY = KY + INCY
    80         CONTINUE

diff --git a/BLAS/SRC/sgemm.f b/BLAS/SRC/sgemm.f
index 9a3d9e1..ce264ab 100644 (file)
--- a/BLAS/SRC/sgemm.f
+++ b/BLAS/SRC/sgemm.f
@@ -311,12 +311,10 @@
    60                 CONTINUE
                   END IF
                   DO 80 L = 1,K
-                      IF (B(L,J).NE.ZERO) THEN
-                          TEMP = ALPHA*B(L,J)
-                          DO 70 I = 1,M
-                              C(I,J) = C(I,J) + TEMP*A(I,L)
-   70                     CONTINUE
-                      END IF
+                      TEMP = ALPHA*B(L,J)
+                      DO 70 I = 1,M
+                          C(I,J) = C(I,J) + TEMP*A(I,L)
+   70                 CONTINUE
    80             CONTINUE
    90         CONTINUE
           ELSE
@@ -353,12 +351,10 @@
   140                 CONTINUE
                   END IF
                   DO 160 L = 1,K
-                      IF (B(J,L).NE.ZERO) THEN
-                          TEMP = ALPHA*B(J,L)
-                          DO 150 I = 1,M
-                              C(I,J) = C(I,J) + TEMP*A(I,L)
-  150                     CONTINUE
-                      END IF
+                      TEMP = ALPHA*B(J,L)
+                      DO 150 I = 1,M
+                          C(I,J) = C(I,J) + TEMP*A(I,L)
+  150                 CONTINUE
   160             CONTINUE
   170         CONTINUE
           ELSE

diff --git a/BLAS/SRC/sgemv.f b/BLAS/SRC/sgemv.f
index eef133f..4ec3f06 100644 (file)
--- a/BLAS/SRC/sgemv.f
+++ b/BLAS/SRC/sgemv.f
@@ -278,24 +278,20 @@
           JX = KX
           IF (INCY.EQ.1) THEN
               DO 60 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      DO 50 I = 1,M
-                          Y(I) = Y(I) + TEMP*A(I,J)
-   50                 CONTINUE
-                  END IF
+                  TEMP = ALPHA*X(JX)
+                  DO 50 I = 1,M
+                      Y(I) = Y(I) + TEMP*A(I,J)
+   50             CONTINUE
                   JX = JX + INCX
    60         CONTINUE
           ELSE
               DO 80 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      IY = KY
-                      DO 70 I = 1,M
-                          Y(IY) = Y(IY) + TEMP*A(I,J)
-                          IY = IY + INCY
-   70                 CONTINUE
-                  END IF
+                  TEMP = ALPHA*X(JX)
+                  IY = KY
+                  DO 70 I = 1,M
+                      Y(IY) = Y(IY) + TEMP*A(I,J)
+                      IY = IY + INCY
+   70             CONTINUE
                   JX = JX + INCX
    80         CONTINUE
           END IF

diff --git a/BLAS/SRC/zgbmv.f b/BLAS/SRC/zgbmv.f
index 0e7311a..1a7fb4b 100644 (file)
--- a/BLAS/SRC/zgbmv.f
+++ b/BLAS/SRC/zgbmv.f
@@ -319,26 +319,22 @@
           JX = KX
           IF (INCY.EQ.1) THEN
               DO 60 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      K = KUP1 - J
-                      DO 50 I = MAX(1,J-KU),MIN(M,J+KL)
-                          Y(I) = Y(I) + TEMP*A(K+I,J)
-   50                 CONTINUE
-                  END IF
+                  TEMP = ALPHA*X(JX)
+                  K = KUP1 - J
+                  DO 50 I = MAX(1,J-KU),MIN(M,J+KL)
+                      Y(I) = Y(I) + TEMP*A(K+I,J)
+   50             CONTINUE
                   JX = JX + INCX
    60         CONTINUE
           ELSE
               DO 80 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      IY = KY
-                      K = KUP1 - J
-                      DO 70 I = MAX(1,J-KU),MIN(M,J+KL)
-                          Y(IY) = Y(IY) + TEMP*A(K+I,J)
-                          IY = IY + INCY
-   70                 CONTINUE
-                  END IF
+                  TEMP = ALPHA*X(JX)
+                  IY = KY
+                  K = KUP1 - J
+                  DO 70 I = MAX(1,J-KU),MIN(M,J+KL)
+                      Y(IY) = Y(IY) + TEMP*A(K+I,J)
+                      IY = IY + INCY
+   70             CONTINUE
                   JX = JX + INCX
                   IF (J.GT.KU) KY = KY + INCY
    80         CONTINUE

diff --git a/BLAS/SRC/zgemm.f b/BLAS/SRC/zgemm.f
index f423315..9d422f4 100644 (file)
--- a/BLAS/SRC/zgemm.f
+++ b/BLAS/SRC/zgemm.f
@@ -317,12 +317,10 @@
    60                 CONTINUE
                   END IF
                   DO 80 L = 1,K
-                      IF (B(L,J).NE.ZERO) THEN
-                          TEMP = ALPHA*B(L,J)
-                          DO 70 I = 1,M
-                              C(I,J) = C(I,J) + TEMP*A(I,L)
-   70                     CONTINUE
-                      END IF
+                      TEMP = ALPHA*B(L,J)
+                      DO 70 I = 1,M
+                          C(I,J) = C(I,J) + TEMP*A(I,L)
+   70                 CONTINUE
    80             CONTINUE
    90         CONTINUE
           ELSE IF (CONJA) THEN
@@ -376,17 +374,15 @@
   170                 CONTINUE
                   END IF
                   DO 190 L = 1,K
-                      IF (B(J,L).NE.ZERO) THEN
-                          TEMP = ALPHA*DCONJG(B(J,L))
-                          DO 180 I = 1,M
-                              C(I,J) = C(I,J) + TEMP*A(I,L)
-  180                     CONTINUE
-                      END IF
+                      TEMP = ALPHA*DCONJG(B(J,L))
+                      DO 180 I = 1,M
+                          C(I,J) = C(I,J) + TEMP*A(I,L)
+  180                 CONTINUE
   190             CONTINUE
   200         CONTINUE
           ELSE
 *
-*           Form  C := alpha*A*B**T          + beta*C
+*           Form  C := alpha*A*B**T + beta*C
 *
               DO 250 J = 1,N
                   IF (BETA.EQ.ZERO) THEN
@@ -399,12 +395,10 @@
   220                 CONTINUE
                   END IF
                   DO 240 L = 1,K
-                      IF (B(J,L).NE.ZERO) THEN
-                          TEMP = ALPHA*B(J,L)
-                          DO 230 I = 1,M
-                              C(I,J) = C(I,J) + TEMP*A(I,L)
-  230                     CONTINUE
-                      END IF
+                      TEMP = ALPHA*B(J,L)
+                      DO 230 I = 1,M
+                          C(I,J) = C(I,J) + TEMP*A(I,L)
+  230                 CONTINUE
   240             CONTINUE
   250         CONTINUE
           END IF

diff --git a/BLAS/SRC/zgemv.f b/BLAS/SRC/zgemv.f
index 4e174c9..6242cbb 100644 (file)
--- a/BLAS/SRC/zgemv.f
+++ b/BLAS/SRC/zgemv.f
@@ -285,24 +285,20 @@
           JX = KX
           IF (INCY.EQ.1) THEN
               DO 60 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      DO 50 I = 1,M
-                          Y(I) = Y(I) + TEMP*A(I,J)
-   50                 CONTINUE
-                  END IF
+                  TEMP = ALPHA*X(JX)
+                  DO 50 I = 1,M
+                      Y(I) = Y(I) + TEMP*A(I,J)
+   50             CONTINUE
                   JX = JX + INCX
    60         CONTINUE
           ELSE
               DO 80 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      IY = KY
-                      DO 70 I = 1,M
-                          Y(IY) = Y(IY) + TEMP*A(I,J)
-                          IY = IY + INCY
-   70                 CONTINUE
-                  END IF
+                  TEMP = ALPHA*X(JX)
+                  IY = KY
+                  DO 70 I = 1,M
+                      Y(IY) = Y(IY) + TEMP*A(I,J)
+                      IY = IY + INCY
+   70             CONTINUE
                   JX = JX + INCX
    80         CONTINUE
           END IF