# Matrix multiplication is slower when multithreading in Julia

**URL:** <https://discourse.julialang.org/t/matrix-multiplication-is-slower-when-multithreading-in-julia/56227>\
**Category:** Performance\
**Tags:** question, multithreading, linearalgebra\
**Created:** [March 1, 2021, 3:43am UTC](https://discourse.julialang.org/t/matrix-multiplication-is-slower-when-multithreading-in-julia/56227 "2021-03-01T03:43:52Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![pablof300](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pablof300/32/22410_2.png) [@pablof300](https://discourse.julialang.org/u/pablof300)\
**Post date:** [March 1, 2021, 3:43am UTC](https://discourse.julialang.org/t/matrix-multiplication-is-slower-when-multithreading-in-julia/56227/1 "2021-03-01T03:43:53Z")

</div>

I am working with big matrices (size of 30k rows and ~100 columns). I am doing some matrix multiplication and the process would take around 20 seconds. This is my code:

```julia
    @time begin
        result = -1
        data = -1
    
        for i=1:size
            first_matrix = @view data[i * split,:]
    
            for j=1:size
                second_matrix = @view Qg[j * split,:]
                matrix_multiplication = first_matrix * second_matrix'
                current_sum = sum(matrix_multiplication)
                global result
    
                if current_sum > result
                    result = current_sum
                    data = matrix_multiplication[1,1]
                end
            end
        end
    end

```

Trying to optimize this a little more, I tried to use multi-threading (julia --thread 4) to get better performance.

```julia
    @time begin
        global result = -1
        global data = -1
        lock = ReentrantLock()
    
        for i=1:size
            first_matrix = @view data[i * split,:]
    
            Threads.@threads for j=1:size
                second_matrix = @view Qg[j * split,:]
                matrix_multiplication = first_matrix * second_matrix'
                current_sum = sum(matrix_multiplication)
                global result
    
                if current_sum > result
                    lock(lock)
                    result = current_sum
                    data = matrix_multiplication[1,1]
                    unlock(lock)
                end
            end
        end
    end

```

By adding multi-threading I thought I would get an increase in performance, but the performance got worse (~40 seconds). I removed the lock to see if that was the issue, but still got the same performance. I am running this on a Dual-Core Intel Core i5 (MacBook pro). Does anyone know why my multi-threading code doesn’t work?

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [March 1, 2021, 6:32am UTC](https://discourse.julialang.org/t/matrix-multiplication-is-slower-when-multithreading-in-julia/56227/2 "2021-03-01T06:32:13Z")

</div>

Three things first

- Avoid working in global scope like this, and definitely avoid global variables. Wrap all of this in a function, and make sure you don’t use globals.
- Arrays in Julia are column major, while you are working along rows. This will hurt performance.
- Could you make a self-contained, running code example with input data (random), so others can try running your code?

Don’t think about multithreading until your single threaded performance is good.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [March 1, 2021, 7:03am UTC](https://discourse.julialang.org/t/matrix-multiplication-is-slower-when-multithreading-in-julia/56227/3 "2021-03-01T07:03:32Z")

</div>

a) Julia’s BLAS is already multithreaded. Try `using LinearAlgebra; BLAS.set_num_threads(1)` to make it single threaded when using `Threads.@threads` for your multithreading.  
b) BLAS is fastest with 1 thread per core, so use just 2 threads on a dual core CPU.

---

<div class="post-metadata">

**Author:** ![tkf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkf/32/17635_2.png) [@tkf](https://discourse.julialang.org/u/tkf)\
**Post date:** [March 1, 2021, 7:18am UTC](https://discourse.julialang.org/t/matrix-multiplication-is-slower-when-multithreading-in-julia/56227/4 "2021-03-01T07:18:30Z")

</div>

> [@pablof300](#):
>
> ```julia
> if current_sum > result
> lock(lock)
> result = current_sum
> 
> ```

You are reading `result` outside the `lock`. It’s probably an undefined behavior. I think you’d want to check if your program is correct, before asking if it is fast.

Unless you are very aware of the pros and cons of what you are doing, I’d suggest _avoiding locks as much as possible_, if you want efficient _computation_. It is very likely there is a better way to phrase your parallel computation. FYI, my introduction of data parallelism in Julia discusses various high-level data-parallel patterns: [A quick introduction to data parallelism in Julia](https://juliafolds.github.io/data-parallelism/tutorials/quick-introduction/)

> [@pablof300](#):
>
> ```julia
> matrix_multiplication = first_matrix * second_matrix'
> current_sum = sum(matrix_multiplication)
> 
> ```

This and the following code `data = matrix_multiplication[1,1]` indicates that you don’t actually need to materialize `matrix_multiplication` as a matrix. It’d be faster (probably even with a single thread) to compute the sum and `(1, 1)`-th element directly (in one sweep).

---

<div class="post-metadata">

**Author:** ![pablof300](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pablof300/32/22410_2.png) [@pablof300](https://discourse.julialang.org/u/pablof300)\
**Post date:** [March 1, 2021, 3:21pm UTC](https://discourse.julialang.org/t/matrix-multiplication-is-slower-when-multithreading-in-julia/56227/5 "2021-03-01T15:21:37Z")

</div>

How can I avoid the lock if I need to redefine a shared variable among threads?

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [March 1, 2021, 3:41pm UTC](https://discourse.julialang.org/t/matrix-multiplication-is-slower-when-multithreading-in-julia/56227/6 "2021-03-01T15:41:18Z")

</div>

I really meant it when I said you should fix your single threaded code first. You will get much bigger gains just from fixing that, and the parallelization will work better too.

You are starting at the wrong end.

---

<div class="post-metadata">

**Author:** ![tkf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkf/32/17635_2.png) [@tkf](https://discourse.julialang.org/u/tkf)\
**Post date:** [March 1, 2021, 8:59pm UTC](https://discourse.julialang.org/t/matrix-multiplication-is-slower-when-multithreading-in-julia/56227/7 "2021-03-01T20:59:03Z")

</div>

I agree with DNF’s suggestion that first thing to try is to improve the serial code by making more idiomatic Julia.

> [@pablof300](#):
>
> if I need to redefine a shared variable among threads

Do you really _need_ to update the global? I’d suggest avoid mixing how to compute things and what to compute. If you just need to find the maximum and a corresponding value (= what), you don’t need to use globals (= how). My tutorial includes a subsection for an approach to do this “findmax” type of computations that is applicable for serial, threaded, distributed, and GPU computations.

---

<div class="post-metadata">

**Author:** ![e3c6](https://avatars.discourse-cdn.com/v4/letter/e/e79b87/32.png) [@e3c6](https://discourse.julialang.org/u/e3c6)\
**Post date:** [January 20, 2022, 10:53am UTC](https://discourse.julialang.org/t/matrix-multiplication-is-slower-when-multithreading-in-julia/56227/8 "2022-01-20T10:53:50Z")

</div>

> [@Elrod](#):
>
> a) Julia’s BLAS is already multithreaded. Try `using LinearAlgebra; BLAS.set_num_threads(1)` to make it single threaded when using `Threads.@threads` for your multithreading.

If I start Julia with `julia - t 4` (say), each of those 4 Julia threads gets _its own_ pool of BLAS threads?

If the answer is no, is there a way to make each julia thread spawn its own BLAS threads?

---

<div class="post-metadata">

**Author:** ![jpsamaroo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jpsamaroo/32/46804_2.png) [@jpsamaroo](https://discourse.julialang.org/u/jpsamaroo)\
**Post date:** [January 20, 2022, 2:53pm UTC](https://discourse.julialang.org/t/matrix-multiplication-is-slower-when-multithreading-in-julia/56227/9 "2022-01-20T14:53:20Z")

</div>

> [@e3c6](#):
>
> If I start Julia with `julia - t 4` (say), each of those 4 Julia threads gets _its own_ pool of BLAS threads?

No, changing the number of Julia threads does not affect the number of BLAS threads.

> [@e3c6](#):
>
> If the answer is no, is there a way to make each julia thread spawn its own BLAS threads?

`BLAS.set_num_threads(1)` doesn’t make all BLAS calls run on just a single thread (which would suck); it just prevents usage of the BLAS threadpool (IIUC), instead causing BLAS code to be directly executed on the Julia thread invoking it. This still allows parallelism when executing many BLAS operations concurrently, if you execute enough of them.

---

<div class="post-metadata">

**Author:** ![e3c6](https://avatars.discourse-cdn.com/v4/letter/e/e79b87/32.png) [@e3c6](https://discourse.julialang.org/u/e3c6)\
**Post date:** [January 20, 2022, 4:47pm UTC](https://discourse.julialang.org/t/matrix-multiplication-is-slower-when-multithreading-in-julia/56227/10 "2022-01-20T16:47:59Z")

</div>

> [@jpsamaroo](#):
>
> `BLAS.set_num_threads(1)` doesn’t make all BLAS calls run on just a single thread (which would suck); it just prevents usage of the BLAS threadpool (IIUC), instead causing BLAS code to be directly executed on the Julia thread invoking it. This still allows parallelism when executing many BLAS operations concurrently, if you execute enough of them.

I should call `BLAS.set_num_threads(1)` within each thread? Or just one global call?

---

<div class="post-metadata">

**Author:** ![jpsamaroo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jpsamaroo/32/46804_2.png) [@jpsamaroo](https://discourse.julialang.org/u/jpsamaroo)\
**Post date:** [January 20, 2022, 5:35pm UTC](https://discourse.julialang.org/t/matrix-multiplication-is-slower-when-multithreading-in-julia/56227/11 "2022-01-20T17:35:02Z")

</div>

Just once, the setting is global to the Julia session.

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [January 20, 2022, 5:45pm UTC](https://discourse.julialang.org/t/matrix-multiplication-is-slower-when-multithreading-in-julia/56227/12 "2022-01-20T17:45:10Z")

</div>

> [@jpsamaroo](#):
>
> `BLAS.set_num_threads(1)` doesn’t make all BLAS calls run on just a single thread (which would suck); it just prevents usage of the BLAS threadpool (IIUC), instead causing BLAS code to be directly executed on the Julia thread invoking it. This still allows parallelism when executing many BLAS operations concurrently, if you execute enough of them.

While this is true (at least I think so as well), I guess the interesting part is what happens for `BLAS.set_num_threads(N)` where `N>1`. And based on a simple test I just ran, OpenBLAS (default) and MKL behave very differently here!

Test function:

```julia
function f()
    X = rand(1000,1000);
    Y = zeros(Threads.nthreads())
    Threads.@threads for i in 1:Threads.nthreads()
        Y[i] = sum(X * X)
    end
    return Y
end

```

**MKL:**

```julia
# 4 Julia threads

julia> BLAS.set_num_threads(1)

julia> @btime f(); # ~400% CPU usage
  32.890 ms (39 allocations: 45.78 MiB)

julia> BLAS.set_num_threads(2)

julia> @btime f(); # ~800% CPU usage
  18.541 ms (39 allocations: 45.78 MiB)

julia> BLAS.get_config()
LinearAlgebra.BLAS.LBTConfig
Libraries:
└ [ILP64] libmkl_rt.so

```

**OpenBLAS:**

```julia
# 4 Julia threads

julia> BLAS.set_num_threads(1)

julia> @btime f(); # ~400% CPU usage
  36.527 ms (39 allocations: 45.78 MiB)

julia> BLAS.set_num_threads(2)

julia> @btime f(); # ~240% CPU usage
  90.159 ms (39 allocations: 45.78 MiB)

julia> BLAS.get_config()
LinearAlgebra.BLAS.LBTConfig
Libraries:
└ [ILP64] libopenblas64_.so

```

For for OpenBLAS, `BLAS.set_num_threads(2)` seems to control the total number of threads that BLAS uses. So increasing `N` in `BLAS.set_num_threads(N)` from `N=1` to `N=2` might actually lead to the computation running on fewer threads. IMHO, that’s very counterintuitive.

For MKL on the other hand, `BLAS.set_num_threads(2)` makes the BLAS computation use `2 * (# of Julia threads)` many threads.

From the above I’d conclude the following. Let’s say a machine has `C` cores (no hyperthreading) and we have `J` Julia threads. What do we have to choose for `N` in `BLAS.set_num_threads(N)` to utilize all cores? For OpenBLAS, the answer seems to be `N=C` whereas for MKL you’d want to set `N=C/J`.

I could nicely confirm this for MKL where on a `C=40` core machine I did set `J=5` and `N=8` to obtain

```julia
julia> @btime f(); # ~4000% CPU usage
  8.076 ms (39 allocations: 45.78 MiB) 

```

Note that `BLAS.set_num_threads(C)`, which strangely enough seems to be the default (why?!?), is understandably suboptimal:

```julia
julia> @btime f();
  19.382 ms (39 allocations: 45.78 MiB)

julia> BLAS.get_num_threads()
40

```

For OpenBLAS, it’s not as clean since `N=C=40` only gave me

```julia
julia> @btime f(); # ~3000% CPU usage
  17.151 ms (39 allocations: 45.78 MiB)

```

(Note that going to `N=45` didn’t improve the CPU usage.)

I guess the conclusion is: Set `BLAS.set_num_threads(1)` unless you know what you / your BLAS library is doing!

PS: Instead of running simple experiments it would of course make sense to read the OpenBLAS / MKL documentation on this matter. But experiments are just more fun 😃

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [January 21, 2022, 10:10pm UTC](https://discourse.julialang.org/t/matrix-multiplication-is-slower-when-multithreading-in-julia/56227/13 "2022-01-21T22:10:11Z")

</div>

Citing a post by myself from a follow-up discussion on Slack:

> Tried to find some documentation of `OPENBLAS_NUM_THREADS` but couldn’t find too much. There is
> 
> - [Faq · OpenMathLib/OpenBLAS Wiki · GitHub](https://github.com/xianyi/OpenBLAS/wiki/Faq#how-can-i-use-openblas-in-multi-threaded-applications) which essentially just says “always use `OPENBLAS_NUM_THREADS=1` if your application is multithreaded”, and
> - [GitHub - OpenMathLib/OpenBLAS: OpenBLAS is an optimized BLAS library based on GotoBLAS2 1.13 BSD version.](https://github.com/xianyi/OpenBLAS#setting-the-number-of-threads-using-environment-variables) which states that `OPENBLAS_NUM_THREADS` indeed specifies the “ **maximum** number of threads”
> 
> So I guess it at least behaves according to documentation.

It appears that OpenBLAS just starts `OPENBLAS_NUM_THREADS` many threads and then schedules work on them (unless `OPENBLAS_NUM_THREADS == 1` which seems to be special-cased as @jpsamaroo mentioned above).

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [January 21, 2022, 10:26pm UTC](https://discourse.julialang.org/t/matrix-multiplication-is-slower-when-multithreading-in-julia/56227/14 "2022-01-21T22:26:46Z")

</div>

And since we’re already reading documentations, the MKL behavior can be understood e.g. from the [documentation of `mkl_set_num_threads_local`](https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/support-functions/threading-control/mkl-set-num-threads-local.html) (which allows one to set the number of MKL threads per executing thread (i.e. Julia thread). It says

> If the thread-local number is not set or if this number is set to zero in a call to this function, Intel® oneAPI Math Kernel Library functions use the global number of threads.

I read this as `# of local MKL threads == # of global MKL threads` for each application thread (i.e. Julia thread).