The reason is that D*L
is computed by the spmatmul
method from the SparseArrays
package here, whereas mul!(B,D,L)
does not invoke the specialized sparse matrix-matrix multiplication, but the generic_matmatmul!
from the LinearAlgebra
subpackage here which does not exploit sparsity.
There is simply no specialized method mul!(B,D,L)
for B
, D
and L
being of type SparseMatrixCSC
(spdiagm(0=>u)
is also of type SparseMatrixCSC
since there is no specialized type for a sparse diagonal matrix in SparseArrays
). However, if you use a dense diagonal matrix D=Diagonal(u)
(type Diagonal
) then you will see what you likely expect.
n=31
L=sprandn(n*n,n*n,0.01)
u=randn(n*n)
D=Diagonal(u)
B=copy(L)
gives
@btime $B .= $D*$L
38.936 μs (5 allocations: 151.47 KiB)
and
@btime mul!($B, $D, $L)
17.935 μs (0 allocations: 0 bytes)
As a side note your minimal working example was not working, since you forgot to provide L
, which is likely why it took so long for someone to respond .