# Mul! vs SparseMatrix

**URL:** <https://discourse.julialang.org/t/mul-vs-sparsematrix/51720>\
**Category:** Performance\
**Created:** [December 12, 2020, 2:44pm UTC](https://discourse.julialang.org/t/mul-vs-sparsematrix/51720 "2020-12-12T14:44:05Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [December 12, 2020, 2:44pm UTC](https://discourse.julialang.org/t/mul-vs-sparsematrix/51720/1 "2020-12-12T14:44:05Z")

</div>

I’m trying to use `mul!` to compute the product of a sparse matrix and a sparse diagonal matrix. My problem is that `mul!` allocates very litte, but takes a lot more time than a matrix multiply.

Is this what’s supposed to happen?

Herewith an example. The starting matrix L is the discrete Laplacian in two space dimension on a square with zero boundary conditions. The grid has n \times n interior nodes with n=31

```julia
julia> typeof(L)
SparseMatrixCSC{Float64,Int64}

julia> u=rand(n*n,);

julia> D=spdiagm(0=>u);

julia> B=copy(L);

julia> @btime $B .= $D*$L;
  34.406 μs (8 allocations: 171.64 KiB)

julia> @btime mul!($B, $D, $L);
  690.153 ms (6 allocations: 336 bytes)

```

---

<div class="post-metadata">

**Author:** ![moeddel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moeddel/32/18641_2.png) [@moeddel](https://discourse.julialang.org/u/moeddel)\
**Post date:** [December 12, 2020, 10:12pm UTC](https://discourse.julialang.org/t/mul-vs-sparsematrix/51720/2 "2020-12-12T22:12:59Z")

</div>

The reason is that `D*L` is computed by the `spmatmul` method from the `SparseArrays` package [here](https://github.com/JuliaLang/julia/blob/v1.5.3/stdlib/SparseArrays/src/linalg.jl), whereas `mul!(B,D,L)` does not invoke the specialized sparse matrix-matrix multiplication, but the `generic_matmatmul!` from the `LinearAlgebra` subpackage [here](https://github.com/JuliaLang/julia/blob/v1.5.3/stdlib/LinearAlgebra/src/matmul.jl) 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.

```julia
n=31
L=sprandn(n*n,n*n,0.01)
u=randn(n*n)
D=Diagonal(u)
B=copy(L)

```

gives

```julia
@btime $B .= $D*$L
38.936 μs (5 allocations: 151.47 KiB)

```

and

```julia
@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 😉.

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [December 12, 2020, 10:22pm UTC](https://discourse.julialang.org/t/mul-vs-sparsematrix/51720/3 "2020-12-12T22:22:55Z")

</div>

Thanks. That did it. FWIW, here’s L

Do

```julia
julia> L=Lap2d(63);

```

where

“”"

```julia
Lap2d(n)

returns the negative Laplacian in two space dimensions

on n x n grid.

Unit square, homogeneous Dirichlet BC

"""

function Lap2d(n)

h=1/(n+1)

maindiag=4*ones(n^2,)/(h*h)

mdiag=Pair(0,maindiag);

sxdiag=-ones(n^2-1,)/(h*h);

for iz=n:n:n^2-1

sxdiag[iz]=0.0;

end

upxdiag=Pair(1,sxdiag);

lowxdiag=Pair(-1,sxdiag);

sydiag=-ones(n^2-n,)/(h*h);

upydiag=Pair(n,sydiag);

lowydiag=Pair(-n,sydiag);

L2d=spdiagm(lowydiag, lowxdiag, mdiag, upxdiag, upydiag);

return L2d

end

```