Unnecessary add operations performed when multiplying in-place by diagonal matrix

[tristanmontoya]

Left-multiplication of a vector of length N by an N-by-N diagonal matrix should require a total of N multiplications. However, there are certain cases when (at least according to GFlops.jl) Julia appears to perform an additional N additions, effectively doubling the cost of such an operation. For example, the following code, which uses an out-of-place matrix-vector multiplication,

using LinearAlgebra, GFlops
A = Diagonal(rand(10))
x = rand(10)
@count_ops y = A * x

produces the following output:

Flop Counter: 10 flop
┌─────┬─────────┐
│     │ Float64 │
├─────┼─────────┤
│ mul │      10 │
└─────┴─────────┘

However, the equivalent in-place operation,

y = similar(x)
@count_ops mul!(y, A, x)

returns the following:

Flop Counter: 20 flop
┌─────┬─────────┐
│     │ Float64 │
├─────┼─────────┤
│ add │      10 │
│ mul │      10 │
└─────┴─────────┘

It seems like these additions are not necessary; the same operation could be done without any addition as follows:

function LinearAlgebra.mul!(y::AbstractVector, 
    A::Diagonal, x::AbstractVector)
    y[:] = A.diag .* x
    return y
end

Is this an oversight or is there a reason for such behaviour?

[mcabbott]

What version is this on? GFlops.jl appears not to work on 1.9.

This might have been fixed by #44651, as the present code seems to check β and not add if it is zero:

julia/stdlib/LinearAlgebra/src/diagonal.jl

Line 274 in 25e3809

if iszero(beta)

[tristanmontoya]

@mcabbott - you are right, this was fixed by that PR. My tests were on 1.7 but it's all fine on 1.8. Thanks!