Axpby! and 0*NaN

General Usage
1

There might be some inconsistencies regarding 0*NaN: if one has p = 0*NaN, then p is also NaN (ok, treating this as julia convention), however,

using LinearAlgebra
x = rand(2,2)
y = rand(2,2)
y[1,1] = NaN

axpby!(.3, x, 0.0, y)
y == .3*x  # true

y[1,1] should be equal to 0.3*x[1,1]+0.0*NaN, which by julia convention should also be NaN but is not; this is also different from the behavior of StaticArray,

using LinearAlgebra
using StaticArrays

x = rand(2,2)
y = rand(2,2)
y[1,1] = NaN

sx = MMatrix{2,2}(x)
sy = MMatrix{2,2}(y)

axpby!(.3, sx, 0.0, sy)
sy == .3*sx # false

sy[1,1] is indeed NaN. Bug or expected? Thanks!

Note: this calls BLAS.axpby, not the generic routine
https://github.com/JuliaLang/julia/blob/e5ee8da70a646f1f54730e727932ebaeaee256d5/stdlib/LinearAlgebra/src/generic.jl#L1389

2

While not ideal, it is expected. The root cause is that OpenBLAS is not “NaN correct”, it does a check if the “add-value” is 0 and then just ignores it. Julia could chose not to use BLAS or check for NaNs when the “add-value” is 0 but it is a bit odd to have a performance regression for the case that is supposed to be faster

It is a bit similar to the discussion in https://github.com/JuliaLang/julia/issues/22733.

5 Likes
3

It’s not just OpenBLAS. The BLAS design documents (some of which predate IEEE-754) specify short-circuiting for zeros in some routines, and it’s done in the model Fortran implementations.

3 Likes