# Is there an easy way to parallelise matrix multiplication?

**URL:** <https://discourse.julialang.org/t/is-there-an-easy-way-to-parallelise-matrix-multiplication/23266>\
**Category:** Performance\
**Created:** [April 18, 2019, 9:03am UTC](https://discourse.julialang.org/t/is-there-an-easy-way-to-parallelise-matrix-multiplication/23266 "2019-04-18T09:03:55Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![boutor2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/boutor2/32/27869_2.png) [@boutor2](https://discourse.julialang.org/u/boutor2)\
**Post date:** [April 18, 2019, 9:03am UTC](https://discourse.julialang.org/t/is-there-an-easy-way-to-parallelise-matrix-multiplication/23266/1 "2019-04-18T09:03:55Z")

</div>

Hi,  
I am working with large sparse matrices (10^5 x 10^5). I would like to make full use of the cores / threads of my CPU to accelerate the multiplication process. As of now, I am using mul!:

```julia
using SparseArrays,LinearAlgebra
n = 10000
L = sprand(Float64,n,n,0.3)
u0 = sprand(Float64,n,0.3)
u1 = similar(u0)
mul!(u1,L,u0);

```

However, this only uses one thread.  
Is there a package or a trick to parallelise the multiplication? I know I could always do something like:

```julia
function parallel_mul!(u1,L,u0)
 Threads.@threads for i in 1:size(u1)[1]
  for j in 1:size(u1)[2]
   mul!(u1[i,j],L[i,:],u0[:,j])
  end
 end
end

```

But maybe there is something more straightforward? Else, may be it is worth making `parallel_mul!` more robust and implement it in some package?

---

<div class="post-metadata">

**Author:** ![Per](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/per/32/10387_2.png) [@Per](https://discourse.julialang.org/u/Per)\
**Post date:** [April 18, 2019, 9:46am UTC](https://discourse.julialang.org/t/is-there-an-easy-way-to-parallelise-matrix-multiplication/23266/2 "2019-04-18T09:46:40Z")

</div>

I think BLAS will run multi-threaded by default in Julia, but otherwise you can control it via `BLAS.set_num_threads`.

However, a matrix-vector multiplication, as in your example, is limited by memory bandwidth, not CPU, so there’s no point running it multi-threaded.

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [April 18, 2019, 9:53am UTC](https://discourse.julialang.org/t/is-there-an-easy-way-to-parallelise-matrix-multiplication/23266/3 "2019-04-18T09:53:58Z")

</div>

> [@Per](#):
>
> I think BLAS will run multi-threaded by default in Julia, but otherwise you can control it via `BLAS.set_num_threads` .

The matrix is sparse in the example

> [@Per](#):
>
> However, a matrix-vector multiplication, as in your example, is limited by memory bandwidth, not CPU, so there’s no point running it multi-threaded.

This is not true (dense matrix, `n = 50000`)

```julia
julia> BLAS.set_num_threads(1)

julia> @time mul!(u1,L,u0);
  1.429343 seconds (4 allocations: 160 bytes)

julia> BLAS.set_num_threads(8)

julia> @time mul!(u1,L,u0);
  0.577491 seconds (4 allocations: 160 bytes)

```

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [April 18, 2019, 9:57am UTC](https://discourse.julialang.org/t/is-there-an-easy-way-to-parallelise-matrix-multiplication/23266/4 "2019-04-18T09:57:40Z")

</div>

Using a dense RHS will improve performance by a lot:

```julia
julia> @time mul!(u1,L,u0);
  0.584000 seconds (4 allocations: 160 bytes)

julia> u0_d = rand(Float64,n)

julia> u1_d = similar(u0_d)

julia> @time mul!(u1_d,L,u0_d);
  0.022939 seconds (4 allocations: 160 bytes)

```

There is often no advantage of having the vector RHS be sparse unless you have hyper sparse matrices (but then you shouldn’t use CSC sparse format anyway)

For threaded sparse matvec (at least in the transpose case), see [https://github.com/JuliaLang/julia/pull/29525](https://github.com/JuliaLang/julia/pull/29525). It cannot really be merged yet because the threading infrastructure in Base is not there yet.

---

<div class="post-metadata">

**Author:** ![Per](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/per/32/10387_2.png) [@Per](https://discourse.julialang.org/u/Per)\
**Post date:** [April 18, 2019, 6:20pm UTC](https://discourse.julialang.org/t/is-there-an-easy-way-to-parallelise-matrix-multiplication/23266/5 "2019-04-18T18:20:18Z")

</div>

> [@kristoffer.carlsson](#):
>
> This is not true (dense matrix, `n = 50000` )

That’s interesting. On my system, there’s no difference in timing when changing the number of threads.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 18, 2019, 7:48pm UTC](https://discourse.julialang.org/t/is-there-an-easy-way-to-parallelise-matrix-multiplication/23266/6 "2019-04-18T19:48:00Z")

</div>

> [@boutor2](#):
>
> `n = 10000; L = sprand(Float64,n,n,0.3)`

Note that this matrix is neither large nor sparse by the standards of realistic applications, so benchmarks of this case may be misleading.

30% nonzero is not usually sparse enough to be worthwhile using sparse algorithms for (unless the sparsity has a _very_ special pattern). Real sparse matrices typically have only small number of nonzero elements per row, e.g. ≲ 10 is typical in FEM-like situations. A _moderate_ sized problem of this sort (I wouldn’t call it “large”) would be 10⁵×10⁵ with 10⋅10⁵ nonzeros, or 0.01% (`1e-4`) sparsity. (For example `sprand(10^5, 10^5, 1e-4)`, though sparsity patterns are rarely random in practice.)

---

<div class="post-metadata">

**Author:** ![boutor2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/boutor2/32/27869_2.png) [@boutor2](https://discourse.julialang.org/u/boutor2)\
**Post date:** [April 19, 2019, 8:36am UTC](https://discourse.julialang.org/t/is-there-an-easy-way-to-parallelise-matrix-multiplication/23266/7 "2019-04-19T08:36:30Z")

</div>

Thanks all, now it makes sense. I knew about `BLAS.set_num_threads()` but it did not change anything in my case as computation time was limited - as said by @Per - by memory bandwidth.

> [@stevengj](#):
>
> Note that this matrix is neither large nor sparse by the standards of realistic applications, so benchmarks of this case may be misleading.

Right - the actual matrices I am playing with have `1e-3` sparsity, with a block pattern, where in this case it really makes sense to use sparse algorithms.

> [@kristoffer.carlsson](#):
>
> Using a dense RHS will improve performance by a lot:
> 
> ```julia
> 
> ```

Wow, yes it changes a lot! I don’t really understand why but, it is for sure the way to go.

Now using dense matrices, BLAS runs indeed multi-threaded and setting multiple threads improve a lot the operation.

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [April 19, 2019, 8:44am UTC](https://discourse.julialang.org/t/is-there-an-easy-way-to-parallelise-matrix-multiplication/23266/8 "2019-04-19T08:44:18Z")

</div>

> [@boutor2](#):
>
> Thanks all, now it makes sense. I knew about `BLAS.set_num_threads()` but it did not change anything in my case as computation time was limited - as said by @Per - by memory bandwidth.

You where using a sparse matrix which is not using BLAS but julias native sparse matrix multiplier, so how many threads BLAS uses has no effect.

---

<div class="post-metadata">

**Author:** ![boutor2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/boutor2/32/27869_2.png) [@boutor2](https://discourse.julialang.org/u/boutor2)\
**Post date:** [April 19, 2019, 8:53am UTC](https://discourse.julialang.org/t/is-there-an-easy-way-to-parallelise-matrix-multiplication/23266/9 "2019-04-19T08:53:34Z")

</div>

Oh, right… Got it !