diff options
author | langou <langou@users.noreply.github.com> | 2014-07-12 19:58:39 +0000 |
---|---|---|
committer | langou <langou@users.noreply.github.com> | 2014-07-12 19:58:39 +0000 |
commit | 2b911746d5d7173ccee0067a4ad09159cdfce1a2 (patch) | |
tree | dbd4feb2c6a281cc6cdd9d124a49a0dcb431c1af | |
parent | 986c48f6312031373a6fa11220cbbcb5556e6782 (diff) | |
download | lapack-2b911746d5d7173ccee0067a4ad09159cdfce1a2.tar.gz lapack-2b911746d5d7173ccee0067a4ad09159cdfce1a2.tar.bz2 lapack-2b911746d5d7173ccee0067a4ad09159cdfce1a2.zip |
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
-rw-r--r-- | BLAS/SRC/cgbmv.f | 28 | ||||
-rw-r--r-- | BLAS/SRC/cgemm.f | 32 | ||||
-rw-r--r-- | BLAS/SRC/cgemv.f | 24 | ||||
-rw-r--r-- | BLAS/SRC/dgbmv.f | 28 | ||||
-rw-r--r-- | BLAS/SRC/dgemm.f | 20 | ||||
-rw-r--r-- | BLAS/SRC/dgemv.f | 24 | ||||
-rw-r--r-- | BLAS/SRC/sgbmv.f | 28 | ||||
-rw-r--r-- | BLAS/SRC/sgemm.f | 20 | ||||
-rw-r--r-- | BLAS/SRC/sgemv.f | 24 | ||||
-rw-r--r-- | BLAS/SRC/zgbmv.f | 28 | ||||
-rw-r--r-- | BLAS/SRC/zgemm.f | 32 | ||||
-rw-r--r-- | BLAS/SRC/zgemv.f | 24 |
12 files changed, 130 insertions, 182 deletions
diff --git a/BLAS/SRC/cgbmv.f b/BLAS/SRC/cgbmv.f index cd597f90..62001f2c 100644 --- 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 ecfe12e8..8d92e240 100644 --- 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 507d19e1..5447086d 100644 --- 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 4a608bd6..b20516fc 100644 --- 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 45d001b7..2136740e 100644 --- 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 675257fa..6c29338d 100644 --- 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 797ac7fc..26e51160 100644 --- 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 9a3d9e1a..ce264abe 100644 --- 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 eef133f3..4ec3f061 100644 --- 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 0e7311a0..1a7fb4b9 100644 --- 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 f4233155..9d422f40 100644 --- 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 4e174c95..6242cbb2 100644 --- 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 |