# Multithreaded broadcast?

**URL:** <https://discourse.julialang.org/t/multithreaded-broadcast/26786>\
**Category:** Internals & Design\
**Tags:** multithreading, broadcast\
**Created:** [July 25, 2019, 9:08am UTC](https://discourse.julialang.org/t/multithreaded-broadcast/26786 "2019-07-25T09:08:08Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [July 25, 2019, 9:08am UTC](https://discourse.julialang.org/t/multithreaded-broadcast/26786/1 "2019-07-25T09:08:08Z")

</div>

So, with the new improvements to threading, I’m wondering how that can be exposed easily to users. One of the most common parallelizable pattern is broadcasting (discussed in [https://github.com/JuliaLang/julia/issues/19777](https://github.com/JuliaLang/julia/issues/19777)). The conservative approach is to make `@threads` work for these, so we can do `@threads a.= b.+ c`. This is fine, but adding annotations gets annoying really fast (“so, what combination of `@threads @simd @inbounds @.` do I need this time?”), especially in matlab/numpy vector code with a lot of such operations. The other non-breaking approach is to have a MultithreadedArray type that implements a custom broadcast; this is also annoying. The disrupting approach is to do that automatically, which would really be an awesome feature. As far as I know matlab also does this (for a subset of builtin functions).

There are two aspects to whether that is doable and desirable: correctness and performance.

Correctness is an issue in general, because there is no way to know whether the function calls will be thread safe. However, it doesn’t seem that bad in practice. I did a `for i in $(find . -name *.jl) ; do grep "\.(" $i ; done` in my `.julia/packages`, and most of the dot calls are of simple functions: exp/log/sin/cos/abs/min, type conversion and so on. There are a couple of more complex functions, but all I’ve checked look OK. Of course some existing code will likely break, but that feels like the kind of major change for which it’s worth it to break stuff.

Regarding performance, I seem to be getting a constant overhead of about 5 microseconds. For a simple `x .= x .+ y`, the crossover point is about `N=50_000` with two threads on my laptop. For a more complex `x .= abs.(cos.(sqrt.(x)))` it was about `N=3_000`. That’s tricky, because broadcast is commonly used both for small and large arrays, so some kind of thresholding would have to be implemented for generic code to make use of multithreading.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [July 25, 2019, 9:19am UTC](https://discourse.julialang.org/t/multithreaded-broadcast/26786/2 "2019-07-25T09:19:14Z")

</div>

I’m not really in favor of automatic multithreading if it is not 100% certain that it will make things faster. Having it opt-in would allow me to choose for myself. How about

```julia
@threads begin
    a .= b .+c
    ...
    more code ...
end

```

?

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [July 25, 2019, 9:41am UTC](https://discourse.julialang.org/t/multithreaded-broadcast/26786/3 "2019-07-25T09:41:24Z")

</div>

There are really two questions here : should threading for broadcast be automatic or should it be opt-in, and should it threshold. Even if threading is opt-in, there still is a dilemma for package writers (eg diff eq, optim, etc) : without some kind of thresholding, adding threads speeds up on large arrays, but slows down on small ones.

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [July 25, 2019, 12:48pm UTC](https://discourse.julialang.org/t/multithreaded-broadcast/26786/4 "2019-07-25T12:48:44Z")

</div>

I’m hoping it’ll be eventually both automatic and opt-in. We have quite a bit of information about the expression that we can use to determine whether it’s safe and/or big enough. But that will necessarily be over-conservative.

On the safety front, it’s much worse than function calls that immediately appear within a broadcasted expression. There’s the [BitArray problem](https://github.com/JuliaLang/julia/issues/27482) — arrays can expose multiple “independent” elements that are actually backed by the same byte.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [July 25, 2019, 1:38pm UTC](https://discourse.julialang.org/t/multithreaded-broadcast/26786/5 "2019-07-25T13:38:34Z")

</div>

FWIW, I just tested the overhead of C/OpenMP threading; for a simple sum the crossover point appears to be at N=3\_000, which would seems to suggest that the overhead in julia is ~20x larger than in OpenMP.

---

<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:** [July 28, 2019, 11:24pm UTC](https://discourse.julialang.org/t/multithreaded-broadcast/26786/6 "2019-07-28T23:24:23Z")

</div>

> [@mbauman](#):
>
> I’m hoping it’ll be eventually both automatic and opt-in. We have quite a bit of information about the expression that we can use to determine whether it’s safe and/or big enough. But that will necessarily be over-conservative.

I’m guessing that it’d be a challenging project as you’d need to know if the function is pure-ish (I mean not `@pure`-pure). For example, multi-threading `xs .+ rand.(Ref(rng))` is a bad idea. Maybe the only way to automate “thread safety” detection would be to implement a bunch of function traits like `ispure(::Type{Tuple{typeof(+),Float64,Float64}}) = true`?

By the way, speaking of RNG, I suppose a proper way to use it with threaded broadcasting is to define thread-able RNG where you create different states for each thread using the jump function. I think the threaded `copyto!` needs to be aware of this thread-able RNG so that appropriate state can be used for each thread.

Given the generality of broadcasting and heuristics needed to solve the problem, I think it makes sense for “thread-ability hinting” mechanism to evolve outside the `Base`. Thanks to the very flexible broadcasting mechanism, it’s pretty simple to implement a basic version based on `BroadcastStyle`:

```julia
module ThreadHints

export threaded

using Base.Broadcast: Broadcasted, BroadcastStyle, throwdm, preprocess

function threaded end

struct ThreadedStyle <: BroadcastStyle end

Broadcast.broadcasted(f::typeof(threaded), x) =
    Broadcasted{ThreadedStyle}(identity, (x,))

@inline function Base.copyto!(dest::AbstractArray, bc::Broadcasted{ThreadedStyle})
    axes(dest) == axes(bc) || throwdm(axes(dest), axes(bc))
    @assert bc.f === identity
    @assert bc.args isa Tuple{Any}
    bc′ = preprocess(dest, bc.args[1])
    Threads.@threads for I in eachindex(bc′)
        @inbounds dest[I] = bc′[I]
    end
    return dest
end

end # module
threaded = ThreadHints.threaded

using BenchmarkTools
N = 50_000
xs = randn(N)
ys = randn(N)
zs = randn(N)
@btime $xs .= $threaded.($ys .+ $zs)
# 26.005 μs (13 allocations: 928 bytes) # 1 thread
# 15.832 μs (34 allocations: 3.13 KiB) # 4 threads

```

---

<div class="post-metadata">

**Author:** ![juthohaegeman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juthohaegeman/32/8620_2.png) [@juthohaegeman](https://discourse.julialang.org/u/juthohaegeman)\
**Post date:** [July 29, 2019, 8:29am UTC](https://discourse.julialang.org/t/multithreaded-broadcast/26786/7 "2019-07-29T08:29:28Z")

</div>

For regular arrays, or, more generally, arrays with a strided memory layout, you can use the `@strided` macro from the [Strided.jl](https://github.com/Jutho/Strided.jl) package.

It will make broadcasting multithreaded, and do a few other things, like automatic views (if they are still strided, that is, if they are indexed by ranges, not by arbitrary vectors) and a view-behaviour for `permutedims` and certain `reshape` operations.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [July 29, 2019, 9:22am UTC](https://discourse.julialang.org/t/multithreaded-broadcast/26786/8 "2019-07-29T09:22:56Z")

</div>

Awesome! Your examples are for quite large arrays, do you have plans to deal with small arrays by disabling multithreading?

---

<div class="post-metadata">

**Author:** ![juthohaegeman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juthohaegeman/32/8620_2.png) [@juthohaegeman](https://discourse.julialang.org/u/juthohaegeman)\
**Post date:** [July 29, 2019, 10:11am UTC](https://discourse.julialang.org/t/multithreaded-broadcast/26786/9 "2019-07-29T10:11:45Z")

</div>

There is a minimal block size which is just handled by a single thread, so it won’t use multiple threads for small arrays. But I am sure all the infrastructure to analyse the operations involved introduce some overhead which can dominate the actual computation for small arrays.

It will be interesting to see if using the new `@spawn` in Julia 1.3 allows to simplify/improve the implementation.

---

<div class="post-metadata">

**Author:** ![gdkrmr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdkrmr/32/2791_2.png) [@gdkrmr](https://discourse.julialang.org/u/gdkrmr)\
**Post date:** [February 10, 2021, 3:21pm UTC](https://discourse.julialang.org/t/multithreaded-broadcast/26786/10 "2021-02-10T15:21:41Z")

</div>

Is there any news on this?

---

<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:** [March 19, 2021, 1:09pm UTC](https://discourse.julialang.org/t/multithreaded-broadcast/26786/11 "2021-03-19T13:09:53Z")

</div>

Nope

---

<div class="post-metadata">

**Author:** ![dilumaluthge](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dilumaluthge/32/29283_2.png) [@dilumaluthge](https://discourse.julialang.org/u/dilumaluthge)\
**Post date:** [March 20, 2021, 3:53pm UTC](https://discourse.julialang.org/t/multithreaded-broadcast/26786/12 "2021-03-20T15:53:53Z")

</div>

I think @chriselrod is planning on adding multithreaded broadcasting to LoopVectorization.jl.

---

<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 21, 2021, 3:43am UTC](https://discourse.julialang.org/t/multithreaded-broadcast/26786/13 "2021-03-21T03:43:24Z")

</div>

Now in LoopVectorization v0.12.5:

```julia
julia> x = rand(3000); y = similar(x);

julia> using LoopVectorization

julia> @benchmark @avxt $y .= abs.(cos.(sqrt.($x))) # multithreaded
BenchmarkTools.Trial:
  memory estimate: 0 bytes
  allocs estimate: 0
  --------------
  minimum time: 988.300 ns (0.00% GC)
  median time: 1.016 μs (0.00% GC)
  mean time: 1.021 μs (0.00% GC)
  maximum time: 13.133 μs (0.00% GC)
  --------------
  samples: 10000
  evals/sample: 10

julia> y ≈ abs.(cos.(sqrt.(x)))
true

julia> @benchmark @avx $y .= abs.(cos.(sqrt.($x))) # single-threaded
BenchmarkTools.Trial:
  memory estimate: 0 bytes
  allocs estimate: 0
  --------------
  minimum time: 2.870 μs (0.00% GC)
  median time: 3.088 μs (0.00% GC)
  mean time: 3.108 μs (0.00% GC)
  maximum time: 7.177 μs (0.00% GC)
  --------------
  samples: 10000
  evals/sample: 9

julia> @benchmark $y .= abs.(cos.(sqrt.($x))) # base
BenchmarkTools.Trial:
  memory estimate: 0 bytes
  allocs estimate: 0
  --------------
  minimum time: 14.779 μs (0.00% GC)
  median time: 15.066 μs (0.00% GC)
  mean time: 15.167 μs (0.00% GC)
  maximum time: 50.096 μs (0.00% GC)
  --------------
  samples: 10000
  evals/sample: 1

julia> versioninfo()
Julia Version 1.7.0-DEV.713
Commit ae26fc6d5f* (2021-03-16 06:50 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, tigerlake)
Environment:
  JULIA_NUM_THREADS = 8

```

At the cost of composability, overhead is substantially smaller than 5 microseconds, which is why on a 4-core CPU it can speed up the `abs.(cos.(sqrt.(x)))` example to 1 microsecond.  
15x faster than `Base` broadcasting.

---

<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:** [March 21, 2021, 3:53am UTC](https://discourse.julialang.org/t/multithreaded-broadcast/26786/14 "2021-03-21T03:53:04Z")

</div>

I would like you to know that you personally have made it so I have to spend as much time precompiling 1.6 as I did 1.5 since the 12.x release of `LoopVectorization` has been getting incredible releases multiple times a week.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [March 21, 2021, 6:55am UTC](https://discourse.julialang.org/t/multithreaded-broadcast/26786/15 "2021-03-21T06:55:33Z")

</div>

Awesome! I’m marking this as a solution until the day we have fast composable automatic threading of broadcasts 🙂

---

<div class="post-metadata">

**Author:** ![jzr](https://avatars.discourse-cdn.com/v4/letter/j/eb9ed0/32.png) [@jzr](https://discourse.julialang.org/u/jzr)\
**Post date:** [March 21, 2021, 6:06pm UTC](https://discourse.julialang.org/t/multithreaded-broadcast/26786/16 "2021-03-21T18:06:47Z")

</div>

> [@Elrod](#):
>
> At the cost of composability

Is noncomposability a temporary issue or a necessary consequence of LoopVectorization’s design?

---

<div class="post-metadata">

**Author:** ![longemen3000](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/longemen3000/32/7298_2.png) [@longemen3000](https://discourse.julialang.org/u/longemen3000)\
**Post date:** [March 22, 2021, 12:03am UTC](https://discourse.julialang.org/t/multithreaded-broadcast/26786/17 "2021-03-22T00:03:28Z")

</div>

> [@jzr](#):
>
> Is noncomposability a temporary issue or a necessary consequence of LoopVectorization’s design?

as far as i understand, LoopVectorization only works with cpu numbers (floats, ints). base broadcast can receive anything as input

---

<div class="post-metadata">

**Author:** ![jzr](https://avatars.discourse-cdn.com/v4/letter/j/eb9ed0/32.png) [@jzr](https://discourse.julialang.org/u/jzr)\
**Post date:** [March 22, 2021, 12:12am UTC](https://discourse.julialang.org/t/multithreaded-broadcast/26786/18 "2021-03-22T00:12:23Z")

</div>

I believe “composability” in @Elrod’s post refers to [composable multithreaded parallelism](https://julialang.org/blog/2019/07/multithreading/).

---

<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 22, 2021, 12:53am UTC](https://discourse.julialang.org/t/multithreaded-broadcast/26786/19 "2021-03-22T00:53:03Z")

</div>

> [@jzr](#):
>
> Is noncomposability a temporary issue or a necessary consequence of LoopVectorization’s design?

`@avxt` already composes with `CheapThreads.batch`.

Should be temporary. tkf made a good suggestion on how to fix it, so I’ll probably do that fairly soon.

The fix would be to make the `wait` function it uses for synchronization check if tasks were started yet. If they weren’t, it’ll take the tasks back.  
The tasks get sent to specific threads, so if that thread is running an `@spawn`ed task, it’d have to yield before LV’s thread can actually begin. If it never does, the task will never start, and the thread running the `@avxt` code would wait forever.  
So the solution is check if it’s been started. If not, remove the task and have this main thread begin. That’d mean it runs serially, but that’d be a lot better than deadlocking or waiting who-knows-how-long before the task can start.

---

<div class="post-metadata">

**Author:** ![jzr](https://avatars.discourse-cdn.com/v4/letter/j/eb9ed0/32.png) [@jzr](https://discourse.julialang.org/u/jzr)\
**Post date:** [May 5, 2021, 10:16pm UTC](https://discourse.julialang.org/t/multithreaded-broadcast/26786/20 "2021-05-05T22:16:54Z")

</div>

> [@Elrod](#):
>
> Should be temporary. tkf made a good suggestion on how to fix it, so I’ll probably do that fairly soon.

Is there a tracking issue for this somewhere?

[Next page](https://discourse.julialang.org/t/multithreaded-broadcast/26786.md?page=2)