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

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â€™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.