# Questions on parallel programming and Threads.@threads vs @parallel with SharedArray

**URL:** https://discourse.julialang.org/t/questions-on-parallel-programming-and-threads-threads-vs-parallel-with-sharedarray/7985
**Category:** Julia at Scale
**Created:** [December 25, 2017, 11:42am UTC](https://discourse.julialang.org/t/questions-on-parallel-programming-and-threads-threads-vs-parallel-with-sharedarray/7985 "2017-12-25T11:42:25Z")
**Posts on this page:** 8
**Page:** 1

<div class="post-metadata">

### Author: ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)
#### Post date: [December 25, 2017, 11:42am UTC](https://discourse.julialang.org/t/questions-on-parallel-programming-and-threads-threads-vs-parallel-with-sharedarray/7985/1 "2017-12-25T11:42:25Z")

</div>

Hi!

I have a few questions on parallel programming and multi-threading. If I understand correctly, the `Threads.@threads` feature is modeled after OpenMP’s and Cilk’s multi-threading schemes, where one and only one thread is running on one logical processor, but somehow the heap memory is shared between all processors/processes/threads. The stack on the other hand is split between the threads. This seems verifiable from the following code which allocates constant memory regardless of the input vector size, and the input vector is a normal array not a `SharedArray`.

```julia
using BenchmarkTools

function g(a, n)
   @inbounds Threads.@threads for i in 1:n
       a[i] = i*rand() + rand()^2 * (i - 1)^2 - 2;
   end
   return
end

n = 100000; a = zeros(Float64, n); @btime g(a, n)
# 1.573 ms (1 allocation: 32 bytes)

```

The above code is run with `Threads.nthreads()` equal 8.

1. Can someone please verify the above understanding? If my understanding is correct, then every stack-allocated variable read from a threaded block will be copied over to each thread’s stack when creating new threads, which is useful to keep in mind.

In contrast, `@parallel` requires an array to be a `SharedArray` to be visible to all logical processors/processes/threads or it needs to be sent over, so the heap and stack memories are not shared by default. The following shows a comparison of the running times and allocation using `Threads.@threads` and `@parallel`, which shows a clear advantage to the `@parallel` construct even though `@parallel` uses only 7 workers for 8 processes while `Threads.@threads` uses all the 8 threads AFAICT (correct me if I am wrong). On the other hand, `@parallel` has a larger overhead in both time and allocations which reverse the ranking when testing with simpler functions.

1. So can someone explain the difference below, when it would be favorable to use `Threads.@threads` over `@parallel`, and if the following benchmarks are expected to change in v1.0? I currently don’t have access to nightly or Linux so this is on v0.6.1 on a Windows machine.

```julia
function f(a, n)
   @inbounds @parallel for i in 1:n
       a[i] = i*rand() + rand()^2 * (i - 1)^2 - 2;
   end
   return
end

function g(a, n)
   @inbounds Threads.@threads for i in 1:n
       a[i] = i*rand() + rand()^2 * (i - 1)^2 - 2;
   end
   return
end

n = 100000; a = SharedArray{Float64}(n); @btime f(a, n)
# 399.680 μs (942 allocations: 38.25 KiB)

n = 100000; a = zeros(Float64, n); @btime g(a, n)
# 1.661 ms (1 allocation: 32 bytes)

```

1. Does the Julia backend of GPUArrays.jl use `Threads.@threads` or `@parallel` constructs?

2. Is it possible to use `@simd` with multi-threading? It gives an error anywhere I try it.

3. Don’t you think it would be nice to have an inplace `pmap!` like `map!`?

Thanks a lot in advance!

---

<div class="post-metadata">

### Author: ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)
#### Post date: [December 25, 2017, 12:04pm UTC](https://discourse.julialang.org/t/questions-on-parallel-programming-and-threads-threads-vs-parallel-with-sharedarray/7985/2 "2017-12-25T12:04:12Z")

</div>

> [@mohamed82008](#):
>
> Does the Julia backend of GPUArrays.jl use Threads.@threads or @parallel constructs?

It seems to be using coroutines/tasks/cooperative threading (@async routines) which are different from the multi-threading above in that the cooperative threads are not executed in parallel and they are all run on a single logical processor/process at different times.

---

<div class="post-metadata">

### Author: ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)
#### Post date: [December 25, 2017, 8:34pm UTC](https://discourse.julialang.org/t/questions-on-parallel-programming-and-threads-threads-vs-parallel-with-sharedarray/7985/3 "2017-12-25T20:34:58Z")

</div>

Yes, the heap is shared by threads and stacks are separate. But threads can still see each others’ stack memory, so I don’t think the extra copying you mention is needed.

Your example uses the default RNG which is not threaded, so the threads are competing for access. The simplest fix is to use a vector of RNGs, indexed by thread id. This should not only be more correct (i.e. actually produce reliably random numbers), but also avoids cache collisions.

I think you need to use `@sync @parallel` to get true execution times for the multiprocess version. On my Linux system that makes `@parallel` look less efficient than `@threads`. (A factor of 2 here, with the above-mentioned vector of RNGs for threads.)

So `@threads` is likely to be a bit more efficient and gives immediate access to all existing objects, at the cost of needing to watch out for resource contention and avoid thread-unsafe functions. Unfortunately since multithreading is still “experimental”, some causes of contention and lack of safety are not documented. (But Yuyichao has been doing great work to shrink their number!)

`@simd` should be fine in a threaded context, but you may need to build an inner SIMD loop to avoid confusing the threading macro.

---

<div class="post-metadata">

### Author: ![vchuravy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vchuravy/32/8_2.png) [@vchuravy](https://discourse.julialang.org/u/vchuravy)
#### Post date: [December 25, 2017, 8:44pm UTC](https://discourse.julialang.org/t/questions-on-parallel-programming-and-threads-threads-vs-parallel-with-sharedarray/7985/4 "2017-12-25T20:44:13Z")

</div>

> [@mohamed82008](#):
>
> Is it possible to use @simd with multi-threading? It gives an error anywhere I try it.

Yes using `@simd` together with `@threads` is possible. See [julia-parallelism](http://slides.com/valentinchuravy/julia-parallelism#/5/6) for an example.

> [@mohamed82008](#):
>
> Don’t you think it would be nice to have an inplace pmap! like map!?

I would recommend looking into `DistributedArrays.jl` which does offer that.

---

<div class="post-metadata">

### Author: ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)
#### Post date: [December 25, 2017, 9:54pm UTC](https://discourse.julialang.org/t/questions-on-parallel-programming-and-threads-threads-vs-parallel-with-sharedarray/7985/5 "2017-12-25T21:54:07Z")

</div>

> [@Ralph\_Smith](#):
>
> Yes, the heap is shared by threads and stacks are separate. But threads can still see each others’ stack memory, so I don’t think the extra copying you mention is needed.

Hmm, this bit is kind of confusing, is there a reference I can take a look at? I got rid of `rand()` for easier benchmarking, and did the following instead and the results are very much in favor of `@threads`. Thanks!

```julia
using BenchmarkTools
using StaticArrays

function g(a, t, n)
   @inbounds Threads.@threads for i in 1:n
       a[i] = i + sum(t)/length(t) - log(sum(t))/2;
   end
   return
end

function f(a, t, n)
   @inbounds @sync @parallel for i in 1:n
       a[i] = i + sum(t)/length(t) - log(sum(t))/2;
   end
   return
end

n = 10000; a = zeros(Float64, n); t = SVector(Tuple(i for i in 1:1000)); @btime g($a, $t, $n)
# 432.866 μs (1 allocation: 7.94 KiB)

n = 10000; a = SharedArray{Float64}(n); t = SVector(Tuple(i for i in 1:1000)); @btime f($a, $t, $n)
# 2.385 ms (8326 allocations: 318.86 KiB)

n = 10000; a = zeros(Float64, n); t = SVector(Tuple(i for i in 1:100)); @btime g($a, $t, $n)
# 77.310 μs (1 allocation: 896 bytes)

n = 10000; a = SharedArray{Float64}(n); t = SVector(Tuple(i for i in 1:100)); @btime f($a, $t, $n)
# 2.009 ms (2018 allocations: 84.13 KiB)

n = 10000; a = zeros(Float64, n); t = SVector(Tuple(i for i in 1:10)); @btime g($a, $t, $n)
# 30.268 μs (1 allocation: 112 bytes)

n = 10000; a = SharedArray{Float64}(n); t = SVector(Tuple(i for i in 1:10)); @btime f($a, $t, $n)
# 1.918 ms (1547 allocations: 62.64 KiB)

```

Side note: when I used a `Tuple` instead of `SVector`, I was getting an inefficiency when summing tuples of size 15 and more, apparently the tuple was moved to the heap before summing as allocations were proportional to the tuple size. This is not related to parallel programming per se but I think it is relevant here.

---

<div class="post-metadata">

### Author: ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)
#### Post date: [December 25, 2017, 9:57pm UTC](https://discourse.julialang.org/t/questions-on-parallel-programming-and-threads-threads-vs-parallel-with-sharedarray/7985/6 "2017-12-25T21:57:45Z")

</div>

That presentation is awesome, thanks a lot!

Also `DistributedArrays.jl` seems great, but I think it only supports inplace `map!` with `DArray`s, not `SharedArray`s, so I think the latter would still be interesting to have.

Edit: actually this seems trivial to do with `@parallel` and a for loop, so unless the inplace `pmap!` also activates SIMD by default, it’s probably not much value added.

---

<div class="post-metadata">

### Author: ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)
#### Post date: [December 26, 2017, 1:10am UTC](https://discourse.julialang.org/t/questions-on-parallel-programming-and-threads-threads-vs-parallel-with-sharedarray/7985/7 "2017-12-26T01:10:52Z")

</div>

The current implementation of threading uses the `libuv` interface; for Linux this is a thin wrapper around native `pthreads` and I presume the logic is made similar for other platforms. So manual pages for `pthreads` (easy to find online) and related discussions may be helpful.

With regard to stack variables, my understanding is this: `@threads` builds a function which will be called in the context of each thread. Arguments are put onto the appropriate stack for the thread where the function will execute (and local variables may be added there). Variables in earlier stack frames can be referenced by pointer and need not be in the local stack. Apparently large immutables like your `SVector` are copied (once!) to new memory on the heap, and the address of this is one such reference.

The strange handling of large tuples seems to be a trick to avoid overly expensive inference.

---

<div class="post-metadata">

### Author: ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)
#### Post date: [December 26, 2017, 4:16am UTC](https://discourse.julialang.org/t/questions-on-parallel-programming-and-threads-threads-vs-parallel-with-sharedarray/7985/8 "2017-12-26T04:16:50Z")

</div>

Great, thanks for all the comments and pointer!
