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> @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)
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
L being of type
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
Diagonal) then you will see what you likely expect.
@btime $B .= $D*$L
38.936 μs (5 allocations: 151.47 KiB)
@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 .
Thanks. That did it. FWIW, here’s L
returns the negative Laplacian in two space dimensions
on n x n grid.
Unit square, homogeneous Dirichlet BC
L2d=spdiagm(lowydiag, lowxdiag, mdiag, upxdiag, upydiag);