Special case `lmul!(0, ::AbstractArray)`?

Internals & Design
,
1 
julia> using LinearAlgebra

julia> function lmul2!(s::Number, X::AbstractArray)
           iszero(s) && return fill!(X, s)
           @simd for I in eachindex(X)
               @inbounds X[I] = s*X[I]
           end
           X
       end
lmul2! (generic function with 1 method)

julia> A = rand(1000, 1000);

julia> @btime (A -> lmul!(0, A))($A);
  290.313 μs (0 allocations: 0 bytes)

julia> @btime (A -> lmul2!(0, A))($A);
  222.413 μs (0 allocations: 0 bytes)

Perhaps the compiler should transform the loop to a fill!. It does seem to handle 1 specially, transforming the loop to a no-op.

julia> @btime (A -> lmul!(1, A))($A);
  2.883 ns (0 allocations: 0 bytes)

If this is complicated, adding an iszero check may help with performance.

2

The compiler is doing the best it can in the face of NaNs and Infs, I think. 1*[2.3 NaN Inf] is indeed [2.3 NaN Inf], but 0*[2.3 NaN Inf] isn’t [0.0 0.0 0.0].

2 Likes
3

One can debate the merits, but there is precedent for strong zeros in linear algebra functions:

using LinearAlgebra

julia> A = fill(1.0,3,3); B = fill(NaN,3,3); C = fill(NaN,3,3);

julia> mul!(C,A,B,0.0,0.0)
3×3 Matrix{Float64}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0

Even if we don’t like numeric zeros acting as strong zeros, we could certainly adopt false since it already acts that way.

4

Regardless of other discussion,

iszero(s) && return fill!(X, s*one(eltype(X)))

might be appropriate (on the fence, the original looks much cleaner)