# Is the linear system solver \\ also multi threaded in Julia as in Matlab? And how to “multithread” it in Julia?

**URL:** <https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404>\
**Category:** Performance\
**Tags:** linearalgebra, sparse\
**Created:** [September 28, 2020, 1:24pm UTC](https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404 "2020-09-28T13:24:38Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![nico](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nico/32/211667_2.png) [@nico](https://discourse.julialang.org/u/nico)\
**Post date:** [September 28, 2020, 1:24pm UTC](https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404/1 "2020-09-28T13:24:38Z")

</div>

I am trying to compare speed and performance between Matlab and Julia. I am looking at a code that does topology optimization of a continuum structure subjected to a given load. The code I am looking at is the public code topopt88.m: [Efficient topology optimization in MATLAB using 88 lines of code - TopOpt](https://www.topopt.mek.dtu.dk/Apps-and-software/Efficient-topology-optimization-in-MATLAB)

Essentially it is an iterative algorithm where in every iteration a Ax=b system is solved (x=A\b), where the sparse matrix A depends on the structural design (it is the finite element stiffness matrix) and it is updated in every iteration.

In Julia the same code runs slower than Matlab. I have done some code optimization in Julia, declaring types in function definitions, using functions as much as possible, avoiding global variables, and implementing other tips I found in the internet. But Julia is still slower than the same Matlab code (same in the sense of conceptual steps).

My question: since Matlab system solve "" is multi threaded by default, is it true the same for Julia? If not, how to multi thread Julia’s \ operator, or to get speed-ups from parallelization similarly?

P.S.  
The following is the current best solution that I found and that gives the best performance to me. Outside Julia in a terminal in Windows type `set JULIA_NUM_THREADS=4` . Start Julia and use the Linear Algebra module this way: `using LinearAlgebra``BLAS.set_num_threads(1)`

In this way there will be 4 threads, and by setting the BLAS threads to 1 it basically will not use its own pool that seems to conflict with Julia threads pool, or at least to slow the performance for me. I got some hint in this regard here: [Julia Threads vs BLAS threads - #15 by alkorang](https://discourse.julialang.org/t/julia-threads-vs-blas-threads/8914/15) My current BLAS setting in Julia 1.5.1 is `julia> BLAS.vendor()``:openblas64`

Using Pardiso.jl with the default MKL sparse solver didn’t help. I still need to figure out if there is a way to increase performance there too and how.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [September 28, 2020, 1:28pm UTC](https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404/2 "2020-09-28T13:28:37Z")

</div>

Any chance you can post the code you’re using in Julia? Functions like `\` have dozens of different possible methods that will be called depending on the type of arguments they receive. This it is much more predictive to talk about the performance of a specific piece of code than an algorithm in general.

---

<div class="post-metadata">

**Author:** ![nico](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nico/32/211667_2.png) [@nico](https://discourse.julialang.org/u/nico)\
**Post date:** [September 28, 2020, 1:36pm UTC](https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404/3 "2020-09-28T13:36:06Z")

</div>

The full code is a bit too long I feel, but essentially the heart of the code is a while loop until convergence of the design update. In every while loop iteration the following system is solved:

```julia
#FE-ANALYSIS
sK = reshape(KE[:]*(Emin.+xPhys[:]'.^penal*(E0-Emin)),64*nelx*nely)
K = sparse(iK,jK,sK)
K = (K+K')/2
U[freedofs] = K[freedofs,freedofs]\F 

```

where freedofs are indices of rows and columns to be selected for the analysis. K is sparse, symmetric, and positive definite hence invertible; F an U are normal arrays.

Does this help?

---

<div class="post-metadata">

**Author:** ![mzilhao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzilhao/32/9473_2.png) [@mzilhao](https://discourse.julialang.org/u/mzilhao)\
**Post date:** [September 28, 2020, 2:24pm UTC](https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404/4 "2020-09-28T14:24:56Z")

</div>

You should probably try benchmarking your code, but it could be that the way you are using the sparse matrix is non-optimal… See here: [Sparse Arrays · The Julia Language](https://docs.julialang.org/en/v1/stdlib/SparseArrays/)  
Note specially the following:

> Indexing operations, especially assignment, are expensive, when carried out one element at a time. In many cases it may be better to convert the sparse matrix into `(I,J,V)` format using [`findnz`](https://docs.julialang.org/en/v1/stdlib/SparseArrays/#SparseArrays.findnz), manipulate the values or the structure in the dense vectors `(I,J,V)` , and then reconstruct the sparse matrix.

---

<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:** [September 28, 2020, 2:41pm UTC](https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404/5 "2020-09-28T14:41:48Z")

</div>

> [@nico](#):
>
> ```julia
> U[freedofs] = K[freedofs,freedofs]\F 
> 
> ```
> 
> where freedofs are indices of rows and columns to be selected for the analysis. K is sparse, symmetric, and positive definite hence invertible; F an U are normal arrays.

Try

```julia
Symmetric(K[freedofs, freedofs]) \ F

```

even though the matrix is analytically symmetric, it may not be so numerically. You can check with `issymetric(K)`.

---

<div class="post-metadata">

**Author:** ![nico](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nico/32/211667_2.png) [@nico](https://discourse.julialang.org/u/nico)\
**Post date:** [September 28, 2020, 2:54pm UTC](https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404/6 "2020-09-28T14:54:12Z")

</div>

Thank you @kristoffer.carlsson. The symmetry is ensured numerically (and not only analytically) by doing:

```julia
K = (K+K')/2

```

Do you think that Julia would prefer the Symmetric command? Doesn’t this slow down things more? I will try anyway, thanks.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [September 28, 2020, 2:59pm UTC](https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404/7 "2020-09-28T14:59:29Z")

</div>

Even if a matrix is numerically symmetric, telling Julia that often gives a 2x speedup.

---

<div class="post-metadata">

**Author:** ![nico](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nico/32/211667_2.png) [@nico](https://discourse.julialang.org/u/nico)\
**Post date:** [September 28, 2020, 3:12pm UTC](https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404/8 "2020-09-28T15:12:54Z")

</div>

I removed  
`K = (K+K')/2` ,  
and replaced  
`U[freedofs] = K[freedofs,freedofs]\F `  
with  
`U[freedofs] = Symmetric(K[freedofs,freedofs])\F`  
and gained in performance. Thanks!  
I still think that there might be more to it. I am setting the number of threads outside Julia as 4  
`set JULIA_NUM_THREADS=4`,  
and in the script I set  
`using LinearAlgebra`  
`BLAS.set_num_threads(1)`  
and in this way I get the best performance using openblas.  
I feel that mutithreading and sparse matrix definition could be pushed more to gain more in computational speed.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [September 28, 2020, 3:22pm UTC](https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404/9 "2020-09-28T15:22:57Z")

</div>

Does your matrix have banded it blocked-banded structure? If so, using that could easily give you a factor of 10 speedup. Also, those structures are much easier to parallelize than CSC

---

<div class="post-metadata">

**Author:** ![nico](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nico/32/211667_2.png) [@nico](https://discourse.julialang.org/u/nico)\
**Post date:** [September 28, 2020, 3:39pm UTC](https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404/10 "2020-09-28T15:39:12Z")

</div>

The benchmark:

```julia
BenchmarkTools.Trial:
  memory estimate: 16.82 GiB
  allocs estimate: 29092
  --------------
  minimum time: 18.681 s (5.05% GC)
  median time: 18.681 s (5.05% GC)
  mean time: 18.681 s (5.05% GC)
  maximum time: 18.681 s (5.05% GC)
  --------------
  samples: 1
  evals/sample: 1
```

---

<div class="post-metadata">

**Author:** ![mzilhao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzilhao/32/9473_2.png) [@mzilhao](https://discourse.julialang.org/u/mzilhao)\
**Post date:** [September 28, 2020, 3:48pm UTC](https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404/11 "2020-09-28T15:48:45Z")

</div>

I was thinking more of benchmarking each individual call in your sample code, like the cost of

```julia
K = sparse(iK,jK,sK)

```

```julia
K = (K+K')/2

```

and

```julia
K[freedofs,freedofs]

```

and then compare that with the actual time for solving the system.

---

<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:** [September 28, 2020, 3:49pm UTC](https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404/12 "2020-09-28T15:49:09Z")

</div>

For large symmetric positive definite matrix you are also likely better of to use an interative solver like conjugate gradient. Would be a good idea to try that out:

[https://juliamath.github.io/IterativeSolvers.jl/dev/linear\_systems/cg/#IterativeSolvers.cg](https://juliamath.github.io/IterativeSolvers.jl/dev/linear_systems/cg/#IterativeSolvers.cg)

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [September 28, 2020, 3:50pm UTC](https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404/13 "2020-09-28T15:50:44Z")

</div>

How big is the matrix?

---

<div class="post-metadata">

**Author:** ![nico](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nico/32/211667_2.png) [@nico](https://discourse.julialang.org/u/nico)\
**Post date:** [September 28, 2020, 3:53pm UTC](https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404/14 "2020-09-28T15:53:14Z")

</div>

I have a mesh of 300x100 elements in 2D, with 2 degrees of freedom per node, which gives around 60800 dofs.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [September 28, 2020, 3:54pm UTC](https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404/15 "2020-09-28T15:54:13Z")

</div>

Then the matrix must be extremely dense, for some reason. I solve 2D heat conduction models with 1 million dofs in ~5seconds. Try to find the number of nonzeros, total and per row.

---

<div class="post-metadata">

**Author:** ![nico](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nico/32/211667_2.png) [@nico](https://discourse.julialang.org/u/nico)\
**Post date:** [September 28, 2020, 3:55pm UTC](https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404/16 "2020-09-28T15:55:29Z")

</div>

Is there a better way to impose boundary conditions in Julia than writing `K[freedofs,freedofs]`? Maybe this eats some time…

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [September 28, 2020, 3:56pm UTC](https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404/17 "2020-09-28T15:56:11Z")

</div>

I do precisely this in Elfel.jl. No problem, it only takes about 0.1 sec.

---

<div class="post-metadata">

**Author:** ![nico](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nico/32/211667_2.png) [@nico](https://discourse.julialang.org/u/nico)\
**Post date:** [September 28, 2020, 3:56pm UTC](https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404/18 "2020-09-28T15:56:26Z")

</div>

OK thank you very much

---

<div class="post-metadata">

**Author:** ![nico](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nico/32/211667_2.png) [@nico](https://discourse.julialang.org/u/nico)\
**Post date:** [September 28, 2020, 3:56pm UTC](https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404/19 "2020-09-28T15:56:58Z")

</div>

I will try

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [September 28, 2020, 3:58pm UTC](https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404/20 "2020-09-28T15:58:09Z")

</div>

Hey, just in case: check the type of the matrix when you get to solve. To see if it did not get converted to a dense matrix at some point…

[Next page](https://discourse.julialang.org/t/is-the-linear-system-solver-also-multi-threaded-in-julia-as-in-matlab-and-how-to-multithread-it-in-julia/47404.md?page=2)