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

For the record, I have experimented with it, also inspired by a talk I have watched (JuliaCon 2016 | Iterative Methods for Sparse Linear Systems in Julia | Lars Ruthotto - YouTube).
Honestly, in my experience it performed more poorly than just A\b, with A sparse with dimensions (60802, 60802) and the normal LinearAlgebra package (julia> BLAS.vendor() :openblas64).

@time K = sparse(iK,jK,sK)
@time Ksolve = Symmetric(K[freedofs,freedofs])
@time Usolve = Ksolve\F
@time U[freedofs] = Usolve
@time U[freedofs] = Symmetric(K[freedofs,freedofs])\F

gives:

0.024518 seconds (17 allocations: 47.242 MiB)
0.008867 seconds (8 allocations: 17.448 MiB)
0.141280 seconds (50 allocations: 128.413 MiB)
0.000167 seconds
0.149922 seconds (58 allocations: 145.861 MiB)

We can thus say that:

K = sparse(iK,jK,sK)
Ksolve = Symmetric(K[freedofs,freedofs])
Usolve = Ksolve\F
U[freedofs] = Usolve

takes 0.174832 s, and:

      U[freedofs] = Symmetric(K[freedofs,freedofs])\F

takes 0.149922 s.

Note: the size of K is (60802, 60802).

Perhaps a minor thing, but don’t do this. It creates an unnecessary copy. If you need to shape it into an 1D vector write

vec(KE) 

instead.

2 Likes

Are you doing the benchmarking in global scope (you shouldn’t)? And are you making sure not to include compilation time?

Also, you cannot add up runtimes like this, in fact, you should be aware that the @time macro is not well suited to doing microbenchmarks. Instead wrap your code in a function before benchmarking, and use BenchmarkTools.jl.

I am confused: didn’t you previously provide timing of 18 seconds?
0.15 seconds is nothing to complain about for a matrix this size…

18 sec for the full optimization analysis, with about 61 iterations

I see, I didn’t get that. I thought it was just the solve. Ignore my advice.

I think 0.15 sec for a 60k x 60k matrix is reasonable. What did you expect?

It all started from the fact that I am getting the same computational time for a “mirrored” (conceptually) code in Matlab. I was trying to prove to my self that Julia is faster. Instead I managed with some fine tuning to get to Matlab performance, and not better…

Let’s make it faster then. Post your code … :slight_smile:

But did you wrap the code in a function, and do the benchmarking properly?

It would be nice to see both the code and the benchmarking code.

Here we go: https://github.com/pollinico/topopt_jl