Inplace-multiplication of sparse matrices with diagonal matrices is slow

[dkarrasch]

This is the same issue reported in #28934, but for the inplace multiplication. The fix there does not seem to include rmul! and lmul!, as far as I can tell.

EDIT: Now I'm confused. Are there instances of rmul!(A::SparseMatrixCSC, D::Diagonal) and lmul!(D::Diagonal, A::SparseMatrixCSC) to begin with? I can't seem to find them. In case there aren't, is that on purpose?

[dkarrasch]

This issue arose from the observation that Julia 0.6's scale!(A::SparseMatrixCSC, b::Vector) outperforms the current rmul!(A, Diagonal(b)) by (more than) an order of magnitude. Internally, scale!(A, A, b) used to be called. Aliasing the first argument with any other is prohibited in Julia 1.0's mul!. Slightly modifying the 0.6's yields quite nice performance.

using LinearAlgebra, SparseArrays, BenchmarkTools

function scale!(A::SparseMatrixCSC, D::Diagonal{T, <:Vector}) where T
    m, n = size(A)
    b    = D.diag
    (n==length(b)) || throw(DimensionMismatch())
    Anzval = A.nzval
    for col = 1:n, p = A.colptr[col]:(A.colptr[col+1]-1)
        @inbounds Anzval[p] *= b[col]
    end
    A
end

P = sprand(12000, 12000, 0.01)
Q = copy(P);
q = Diagonal(dropdims(sum(P, dims=2), dims=2))

@btime scale!($P, $q)

which should be compared to

@btime rmul!($Q, $q)

which doesn't seem to use sparsity (it is based on broadcasting A .= A .* transpose(D.diag)).

ADDENDUM: The difference to lmul! seems to be even more dramatic:

function scale!(b::Diagonal, A::SparseMatrixCSC)
    m, n = size(A)
    (m==length(b.diag)) || throw(DimensionMismatch())
    Anzval = A.nzval
    Arowval = A.rowval
    for col = 1:n, p = A.colptr[col]:(A.colptr[col+1]-1)
        @inbounds Anzval[p] *= b.diag[Arowval[p]]
    end
    A
end

@btime scale!($q, $P)
@btime lmul!($q, $Q)