Mul! vs SparseMatrix

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 :wink:.

3 Likes