Fix issue #1327

[langou]

Fix issue #1327

AXPBY(N, α, X, INCX, β, Y, INCY) calls SCAL(N, β, Y, INCY) in the case ( α.EQ.0 )

However, as pointed by @venovako in Iisue #1327, SCAL does not handle INCX ≤ 0 so the fix is that we only call SCAL when INCY > 0.

I also fixed the two typos pointed by @venovako

Also, for some reasons, AXPBY(N, α, X, INCX, β, Y, INCY) was not calling SCAL(N, β, Y, INCY) in the case ( α.EQ.0 ) .AND. ( β.EQ.0 ), but I think it is good to call SCAL in this case too so I removed the ( β.NE.0 ) condition.

I do not know why SCAL does not handle INCX ≤ 0. This is weird and maybe a better fix is to change the behavior of SCAL so that SCAL handles INCX ≤ 0. However the behavior seemed to be intentional as mentioned in the comment:

lapack/BLAS/SRC/cscal.f

Line 72 in c1c5b72

*> modified 3/93 to return if incx .le. 0.

Anyone knows why?

[langou]

Reviewing the Level 1 BLAS:

• the same quick return for INCX ≤ 0 and 3/93 comment are present in SCAL, ASUM and IAMAX.
• And then the ROT, ROTM, SWAP, COPY, AXPY, DOT, and NRM2 do not have it and so handle INCX ≤ 0 just fine

I am a little confused as to why someone decided that SCAL, ASUM and IAMAX should not handle INCX ≤ 0.

It seems that if we look at

Lawson, Hanson, Kincaid, Krogh,
Algorithm 539: Basic Linear Algebra Subprograms for Fortran Usage
ACM TOMS 5(3), September 1979, pp. 324–325.
https://www.netlib.org/toms/539.gz

the original codes for SCAL, ASUM and IAMAX were handling INCX ≤ 0 just fine.

[martin-frbg]

Hmm. I do remember that NRM2 only (re?)gained that capability with its "recent" rewrite around three years ago. But if there was any discussion about this aspect back then, I don't recall it

[martin-frbg]

That was #514, citing an Edward Anderson paper from 2017.
That paper notes that support for negative increments was inconsistent in BLAS, being supported in levels 2 and 3 more as an afterthought, while historically the routines operating on a single vector were expected to handle positive increments only, and behaved irrationally on negative inputs as a check was missing "until the 90s" (!)
(sorry, can't provide page number as the pdf doesn't appear to have any)

[langou]

Thanks @martin-frbg, the reference to Ed's paper is great and useful. This is page 12:2 (top left corner of PDF) and reference for the record is:

E. Anderson, Algorithm 978: Safe scaling in the Level 1 BLAS, ACM Trans. Math. Software, 44 (2017), pp. 1–28 (art. no. 12), https://doi.org/10.1145/3061665

Thanks!

[langou]

I would be in favor of having the code for SCAL, ASUM and IAMAX to handle INCX ≤ 0.

[langou]

We need to speak about the current implementation of AXPBY and

      IF (α .EQ. 0) THEN
         CALL SCAL(N, β, Y, INCY)

For AXPBY( α, X, β, Y ), the call to SCAL( β, Y ) when α = 0 means that,

• if α = 0, and X(1) = NaN, then AXPBY, Y ← α X + β Y, would return Y(1) = β Y(1) and so a NaN in X(1) would not propagate to Y(1)
• Also, if α = 0, and X(1) = ±∞, and Y(1) = β Y(1), and so the operation 0 * ±∞ is avoided and no NaN is created in Y(1)

I do not think we have been thinking at exception handling for AXPBY.

Current behavior for AXPBY( α, X, β, Y ), Y ← α X + β Y, is that

• when β = 0, NaNs in Y propagate
• whereas, when α = 0, NaNs in X do not propagate (strong zero)

Is that we want?

For GEMM( α, A, B, β, C ), ( C ← α AB + β C ), the behavior is

• when β = 0, NaNs in AB do not propagate
• when α = 0, NaNs in C do not propagate

For AXPY( α, X, Y ), ( Y ← α X + Y ), the behavior is

• when α = 0, NaNs in X do not propagate

For SCAL( α, X ), ( X ← α X ), the behavior is

• when α = 0, NaNs in X propagates

I think what we want is

Desired behavior for AXPBY( α, X, β, Y ), Y ← α X + β Y,

• when β = 0, NaNs in Y do not propagate to Y
• when α = 0, NaNs in X propagates to X

(Which is the reverse of what we have.)

Opinions welcome.

[martin-frbg]

I would be in favor of having the code for SCAL, ASUM and IAMAX to handle INCX ≤ 0.

It would seem to make sense (though OpenBLAS would probably have to fix assembly implementations across all architectures). I haven't checked what MKL does here but its documentation does not mention any constraints on the increment