# Multi-threaded inverse of a matrix

**URL:** <https://discourse.julialang.org/t/multi-threaded-inverse-of-a-matrix/79543>\
**Category:** General Usage\
**Tags:** blas, multithreading, linearalgebra\
**Created:** [April 15, 2022, 5:18pm UTC](https://discourse.julialang.org/t/multi-threaded-inverse-of-a-matrix/79543 "2022-04-15T17:18:32Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![fipelle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fipelle/32/4772_2.png) [@fipelle](https://discourse.julialang.org/u/fipelle)\
**Post date:** [April 15, 2022, 5:18pm UTC](https://discourse.julialang.org/t/multi-threaded-inverse-of-a-matrix/79543/1 "2022-04-15T17:18:32Z")

</div>

Hi there,

I thought that [`BLAS.set_num_threads(n)`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.BLAS.set_num_threads) was enough to run all `LinearAlgebra` functions in a multi-threaded fashion. However, I came to realise that this actually applies mostly to matrix products. Inverses seem to be excluded from it - I suppose because they are based on LAPACK.

I was wondering if there’s a way to compute inverses in a multi-threaded fashion either within `LinearAlgebra` or via external dependencies. Right now, I would prefer using CPU rather than GPU.

I was also curious to know if there are caps in terms of the number of threads than `LinearAlgebra` can support. As highlighted in a [different post](https://discourse.julialang.org/t/ideal-number-of-blas-threads/79197/3), I have also noticed that running `BLAS.set_num_threads(64)` does not actually imply that 64 threads will be used or initialised (in my experiment the cluster set 32 threads and it used about 10% of the total CPU - comprising 64 cores).

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [April 15, 2022, 5:27pm UTC](https://discourse.julialang.org/t/multi-threaded-inverse-of-a-matrix/79543/2 "2022-04-15T17:27:39Z")

</div>

> [@fipelle](#):
>
> Inverses seem to be excluded from it - I suppose because they are based on LAPACK.

Why do you think that? My test

```julia
using Random
using LinearAlgebra
using BenchmarkTools

Random.seed!(42)

A = rand(2000, 2000)
b = rand(2000)

BLAS.set_num_threads(1)
@btime $A \ $b

BLAS.set_num_threads(6)
@btime $A \ $b

BLAS.set_num_threads(1)
@btime inv($A)

BLAS.set_num_threads(6)
@btime inv($A)

```

shows speedup

```julia
  133.715 ms (4 allocations: 30.55 MiB) # A \ b, 1 thread 
  43.437 ms (4 allocations: 30.55 MiB) # A \ b, 6 threads
  420.827 ms (5 allocations: 31.51 MiB) # inv(A), 1 thread
  173.065 ms (5 allocations: 31.51 MiB) # inv(A), 6 threads

```

---

<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:** [April 15, 2022, 7:05pm UTC](https://discourse.julialang.org/t/multi-threaded-inverse-of-a-matrix/79543/3 "2022-04-15T19:05:16Z")

</div>

> [@fipelle](#):
>
> Inverses seem to be excluded from it - I suppose because they are based on LAPACK.

Inverses aren’t excluded, _because_ they’re based on BLAS/LAPACK (just as matmuls).

> [@fipelle](#):
>
> I was also curious to know if there are caps in terms of the number of threads than `LinearAlgebra` can support. As highlighted in a [different post](https://discourse.julialang.org/t/ideal-number-of-blas-threads/79197/3), I have also noticed that running `BLAS.set_num_threads(64)` does not actually imply that 64 threads will be used or initialised (in my experiment the cluster set 32 threads and it used about 10% of the total CPU - comprising 64 cores).

Yes, IIRC the max number of BLAS threads is set to 32 by Julia. This has changed in 1.8 where you can get, I believe, up to 512 BLAS threads. See [BLAS performance testing for Julia 1.8](https://discourse.julialang.org/t/blas-performance-testing-for-julia-1-8/69520)

---

<div class="post-metadata">

**Author:** ![fipelle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fipelle/32/4772_2.png) [@fipelle](https://discourse.julialang.org/u/fipelle)\
**Post date:** [April 15, 2022, 7:06pm UTC](https://discourse.julialang.org/t/multi-threaded-inverse-of-a-matrix/79543/4 "2022-04-15T19:06:22Z")

</div>

> [@carstenbauer](#):
>
> Yes, IIRC the max number of BLAS threads is set to 32 by Julia. This has changed in 1.8 where you can get, I believe, up to 512 BLAS threads.

Thanks, that must be it then. Is there any way for increasing this number in the LTS version?

---

<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:** [April 15, 2022, 7:08pm UTC](https://discourse.julialang.org/t/multi-threaded-inverse-of-a-matrix/79543/5 "2022-04-15T19:08:39Z")

</div>

No, I don’t think so (apart from changing the number in the source code and compiling from source).

(@viralbshah or do you know a trick here ?)

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [April 15, 2022, 9:18pm UTC](https://discourse.julialang.org/t/multi-threaded-inverse-of-a-matrix/79543/6 "2022-04-15T21:18:57Z")

</div>

> [@carstenbauer](#):
>
> No, I don’t think so (apart from changing the number in the source code and compiling from source).

I do understand @fipelle’s concern. In other words, is there any reason not to change this in an upcoming release of LTS?

---

<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:** [April 16, 2022, 6:07am UTC](https://discourse.julialang.org/t/multi-threaded-inverse-of-a-matrix/79543/7 "2022-04-16T06:07:57Z")

</div>

Don’t get me wrong, I also understand his concern and personally think that this is a footgun. I’m just answering the questions 🙂

(BTW, I think that defaulting to something else than 1 for ‘OPENBLAS\_NUM\_THREADS’ when running Julia in multithreaded mode is a footgun as well…)

> [@goerch](#):
>
> In other words, is there any reason not to change this in an upcoming release of LTS?

Don’t know but I’d guess no. The thing is that patch releases are generally for fixing bugs. Question is whether we can (deliberately) label this as a bug and then fix it. But even if we would, it might not be easy to technically fix this in the LTS because it involves OpenBLAS(\_jll.jl), i.e. an external dependency. Others are more qualified to speak to this. Perhaps filling an issue and at least discussing this would make sense.