LAPACK Archives

[Lapack] NaN propagation in BLAS


Hi all, Okey doc. I just changed this. See log below. This is for the better. 
Cheers, Julien.


Author: langou
Date: Sat Jul 12 15:58:39 2014
New Revision: 1488

URL: https://icl.cs.utk.edu/trac/lapack-dev/changeset/1488
Log:

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@Domain.Removed.

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





Modified:
   lapack/trunk/BLAS/SRC/cgbmv.f
   lapack/trunk/BLAS/SRC/cgemm.f
   lapack/trunk/BLAS/SRC/cgemv.f
   lapack/trunk/BLAS/SRC/dgbmv.f
   lapack/trunk/BLAS/SRC/dgemm.f
   lapack/trunk/BLAS/SRC/dgemv.f
   lapack/trunk/BLAS/SRC/sgbmv.f
   lapack/trunk/BLAS/SRC/sgemm.f
   lapack/trunk/BLAS/SRC/sgemv.f
   lapack/trunk/BLAS/SRC/zgbmv.f
   lapack/trunk/BLAS/SRC/zgemm.f
   lapack/trunk/BLAS/SRC/zgemv.f

Modified: lapack/trunk/BLAS/SRC/cgbmv.f
URL: 
https://icl.cs.utk.edu/trac/lapack-dev/file/lapack/trunk/BLAS/SRC/cgbmv.f?rev=1488
==============================================================================
--- lapack/trunk/BLAS/SRC/cgbmv.f (original)
+++ lapack/trunk/BLAS/SRC/cgbmv.f Sat Jul 12 15:58:39 2014
@@ -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

Modified: lapack/trunk/BLAS/SRC/cgemm.f
URL: 
https://icl.cs.utk.edu/trac/lapack-dev/file/lapack/trunk/BLAS/SRC/cgemm.f?rev=1488
==============================================================================
--- lapack/trunk/BLAS/SRC/cgemm.f (original)
+++ lapack/trunk/BLAS/SRC/cgemm.f Sat Jul 12 15:58:39 2014
@@ -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

Modified: lapack/trunk/BLAS/SRC/cgemv.f
URL: 
https://icl.cs.utk.edu/trac/lapack-dev/file/lapack/trunk/BLAS/SRC/cgemv.f?rev=1488
==============================================================================
--- lapack/trunk/BLAS/SRC/cgemv.f (original)
+++ lapack/trunk/BLAS/SRC/cgemv.f Sat Jul 12 15:58:39 2014
@@ -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

Modified: lapack/trunk/BLAS/SRC/dgbmv.f
URL: 
https://icl.cs.utk.edu/trac/lapack-dev/file/lapack/trunk/BLAS/SRC/dgbmv.f?rev=1488
==============================================================================
--- lapack/trunk/BLAS/SRC/dgbmv.f (original)
+++ lapack/trunk/BLAS/SRC/dgbmv.f Sat Jul 12 15:58:39 2014
@@ -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

Modified: lapack/trunk/BLAS/SRC/dgemm.f
URL: 
https://icl.cs.utk.edu/trac/lapack-dev/file/lapack/trunk/BLAS/SRC/dgemm.f?rev=1488
==============================================================================
--- lapack/trunk/BLAS/SRC/dgemm.f (original)
+++ lapack/trunk/BLAS/SRC/dgemm.f Sat Jul 12 15:58:39 2014
@@ -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

Modified: lapack/trunk/BLAS/SRC/dgemv.f
URL: 
https://icl.cs.utk.edu/trac/lapack-dev/file/lapack/trunk/BLAS/SRC/dgemv.f?rev=1488
==============================================================================
--- lapack/trunk/BLAS/SRC/dgemv.f (original)
+++ lapack/trunk/BLAS/SRC/dgemv.f Sat Jul 12 15:58:39 2014
@@ -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

Modified: lapack/trunk/BLAS/SRC/sgbmv.f
URL: 
https://icl.cs.utk.edu/trac/lapack-dev/file/lapack/trunk/BLAS/SRC/sgbmv.f?rev=1488
==============================================================================
--- lapack/trunk/BLAS/SRC/sgbmv.f (original)
+++ lapack/trunk/BLAS/SRC/sgbmv.f Sat Jul 12 15:58:39 2014
@@ -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

Modified: lapack/trunk/BLAS/SRC/sgemm.f
URL: 
https://icl.cs.utk.edu/trac/lapack-dev/file/lapack/trunk/BLAS/SRC/sgemm.f?rev=1488
==============================================================================
--- lapack/trunk/BLAS/SRC/sgemm.f (original)
+++ lapack/trunk/BLAS/SRC/sgemm.f Sat Jul 12 15:58:39 2014
@@ -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

Modified: lapack/trunk/BLAS/SRC/sgemv.f
URL: 
https://icl.cs.utk.edu/trac/lapack-dev/file/lapack/trunk/BLAS/SRC/sgemv.f?rev=1488
==============================================================================
--- lapack/trunk/BLAS/SRC/sgemv.f (original)
+++ lapack/trunk/BLAS/SRC/sgemv.f Sat Jul 12 15:58:39 2014
@@ -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

Modified: lapack/trunk/BLAS/SRC/zgbmv.f
URL: 
https://icl.cs.utk.edu/trac/lapack-dev/file/lapack/trunk/BLAS/SRC/zgbmv.f?rev=1488
==============================================================================
--- lapack/trunk/BLAS/SRC/zgbmv.f (original)
+++ lapack/trunk/BLAS/SRC/zgbmv.f Sat Jul 12 15:58:39 2014
@@ -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

Modified: lapack/trunk/BLAS/SRC/zgemm.f
URL: 
https://icl.cs.utk.edu/trac/lapack-dev/file/lapack/trunk/BLAS/SRC/zgemm.f?rev=1488
==============================================================================
--- lapack/trunk/BLAS/SRC/zgemm.f (original)
+++ lapack/trunk/BLAS/SRC/zgemm.f Sat Jul 12 15:58:39 2014
@@ -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

Modified: lapack/trunk/BLAS/SRC/zgemv.f
URL: 
https://icl.cs.utk.edu/trac/lapack-dev/file/lapack/trunk/BLAS/SRC/zgemv.f?rev=1488
==============================================================================
--- lapack/trunk/BLAS/SRC/zgemv.f (original)
+++ lapack/trunk/BLAS/SRC/zgemv.f Sat Jul 12 15:58:39 2014
@@ -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



On Jul 10, 2014, at 10:32 AM, James Demmel <demmel@Domain.Removed> wrote:

How about this for a careful (nonreckless) matmul (or other BLAS):
   Spend O(n^2) to check the inputs for NaNs (and perhaps infinities) ... 
cheap
   If (a NaN appears)   ... very unlikely 
        use a version of matmul that we know will propagate them correctly 
            ... could be slower matmul, or we could just fill in selected 
rows and 
            ... columns of product with NaNs depending on where NaNs appeared 
in the input
   else
        use your favorite performance tuned (reckless) version, including 
Strassen etc.
   end if

How to do this of course depends on how we decide to treat the cases alpha=0, 
beta=0.

A less portable version might use some kind of exceptions flags at the end 
(of course
the usual ones will not be set by quiet NaNs) to see if a NaN appeared, and 
recomputation
is necessary (so slow with low probability). This would of course not fix the 
problem identified
in Julien's email, where the NaN may not be touched at all.


Jim

   
On 7/10/14, 8:34 AM, Julien Langou wrote:
Hi all, 

1) Below is an email that we just received from some R & redhat folks about 
NaN propagation in the BLAS. Please read the email discussion between them 
and I. 
I summarized their request in my reply to them. (So you can just read my 
answer.)

2) Please see as well my question to them. Starting with: "A related 
question. How do you feel about the other of IF statements if DGEMM?? (I do 
not have an answer yet from them.)

3) Please see as well: 
http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=5&t=4499

4) Please see as well email from Jim (Demmel) to Penny Andrerson (MathWorks) 
"[LAPACK-DEV] NaNs and consistent exception handling" from DEC-14-2012. 

So basically this is the recurring of issue of NaN propagation.

1) I think we will apply the patch of these guys. Reason: I am OK with 
keeping the non-propagation of NaNs when ALPHA and BETA have special values 
(for now), but this non-propagation here is really weird. It is specific to 
the JLI variant of DGEMM, it applies when an entry of B is ZERO, pfff . . . 
Oh my. I just have a headache to think what the output would look like. 
There would be NaN at some places in C, no NaN at other places. A mess. I am 
OK with keeping the non-propagation of NaNs when ALPHA and BETA have special 
values (for now), but not based on values of the entries of the matrices A 
and B. 

(Does this make sense?)

2) We can rethink the NaN propagation in DGEMM (and BLAS in general) for 
values of ALPHA and BETA. (???) So bottom line: do we want the reference 
DGEMM to just be three loops and propagate everything? I think this is more 
a question for users, super-users, etc. Jim might have some answers with his 
related projects. 

3) Note: Strassen version of DGEMM will not propagate NaN correctly. Too 
many NaNs will be created. (Just a random comment, I was just thinking about 
this.)

Cheers, 
Julien.



Begin forwarded message:

From: Julien Langou <julien.langou@Domain.Removed>
Subject: Re: [Lapack] Issue with BLAS causing regression test failure in R
Date: July 10, 2014 at 9:12:10 AM MDT
To: Tom Callaway <tcallawa@Domain.Removed>
Cc: Martyn Plummer <plummerm@Domain.Removed>, "peljasz@Domain.Removed" 
<peljasz@Domain.Removed>, julie langou <julie@Domain.Removed>


Hi Tom, Hi Martyn, Hi Lejeczek,

Thanks for pointing this, just to let you know that we (LAPACK folks) are 
discussing this. 

So just 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 each L (from 1 to K), whether 
B(L,J) is zero (or not) to save (or not) the 2M following operations.

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

You suggest to remove the IF statement so that NaNs are correctly 
propagated. 

So as I said, we (LAPACK folks) are discussing this.

A related question. How do you feel about the other of IF statements if 
DGEMM? 

For example: (always on DGEMM)

                (1) if BETA = 0, we set C to ZERO, then we follow with 
computation, etc. So if BETA = 0 and NaN in C, the NaN in C are removed. 
Are you OK with this? 

                (2) if ALPHA = 0 and BETA = 1, we return directly. So if 
ALPHA = 0 and BETA = 1, and NaNs in A or in B, these NaNs are not 
propagated to output C. Are you OK with this? 
What I am trying to say is that there are plenty of other places where NaN 
are not necessarily propagated within DGEMM.
(Pretty much each time there is an IF statement in the DGEMM subroutine is 
a case of this.)

So my question is how do you feel about the other places where NaN in input 
(A, B, C, or alpha, or beta) do not propagate?

Cheers, 
Julien.


On Jul 7, 2014, at 9:58 AM, Tom Callaway <tcallawa@Domain.Removed> wrote:

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

==
?.???`?.??`?.??.???`?.?><(((?> OSAS @ Red Hat
University Outreach || Fedora Special Projects || Fedora Legal
<lapack-3.5.0-R-blas-fixes.patch><tcallawa.vcf>_______________________________________________
Lapack mailing list
Lapack@Domain.Removed
http://lists.eecs.utk.edu/mailman/listinfo/lapack




-------------- next part --------------
An HTML attachment was scrubbed...
URL: 
<http://lists.eecs.utk.edu/mailman/private/lapack/attachments/20140712/ae98fb10/attachment-0001.html>

<Prev in Thread] Current Thread [Next in Thread>


For additional information you may use the LAPACK/ScaLAPACK Forum.
Or one of the mailing lists, or