summaryrefslogtreecommitdiff
path: root/BLAS/SRC/sgbmv.f
diff options
context:
space:
mode:
authorlangou <langou@users.noreply.github.com>2014-07-12 19:58:39 (GMT)
committerlangou <langou@users.noreply.github.com>2014-07-12 19:58:39 (GMT)
commit2b911746d5d7173ccee0067a4ad09159cdfce1a2 (patch)
treedbd4feb2c6a281cc6cdd9d124a49a0dcb431c1af /BLAS/SRC/sgbmv.f
parent986c48f6312031373a6fa11220cbbcb5556e6782 (diff)
downloadlapack-2b911746d5d7173ccee0067a4ad09159cdfce1a2.zip
lapack-2b911746d5d7173ccee0067a4ad09159cdfce1a2.tar.gz
lapack-2b911746d5d7173ccee0067a4ad09159cdfce1a2.tar.bz2
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
Diffstat (limited to 'BLAS/SRC/sgbmv.f')
-rw-r--r--BLAS/SRC/sgbmv.f28
1 files changed, 12 insertions, 16 deletions
diff --git a/BLAS/SRC/sgbmv.f b/BLAS/SRC/sgbmv.f
index 797ac7f..26e5116 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