# How to speed up this simple code? Multithreading, simd, inbounds

**URL:** <https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681>\
**Category:** Performance\
**Created:** [January 16, 2019, 1:50am UTC](https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681 "2019-01-16T01:50:55Z")\
**Posts on this page:** 20\
**Page:** 2

<div class="post-metadata">

**Author:** ![Juan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juan/32/7657_2.png) [@Juan](https://discourse.julialang.org/u/Juan)\
**Post date:** [January 17, 2019, 4:56pm UTC](https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681/21 "2019-01-17T16:56:29Z")

</div>

Oh, I was playing with and without Int32 and trying more modifications, but if y the input is large you get overflow errors.  
In general I’m interested in improving the speed using all tricks. But in this case I wanted to compare it to R using the same formula/algorithm on both platforms in order to get a fair comparison.  
“Avoiding unnecessary allocations” and “removing type instabilities” was other of my aims.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 17, 2019, 5:15pm UTC](https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681/22 "2019-01-17T17:15:03Z")

</div>

> [@Juan](#):
>
> But in this case I wanted to compare it to R using the same formula/algorithm on both platforms in order to get a fair comparison.

R is an interpreted language, so even mildly idiomatic Julia code **with loops** will blow it out of the water. Eg

```R
library(tictoc)

mysum <- function(n) {
    acc <- 0
    for (i in 1:n) {
        for (j in 1:i) {
            acc <- acc + j*j
        }
    }
    acc
}

```

gives

```R
> tic()
> mysum(10000)
[1] 8.336667e+14
> toc()
2.544 sec elapsed

```

Cf

```julia
using BenchmarkTools
function mysum(n)
    acc = 0
    for i in 1:n
        for j in 1:i
            acc += j*j
        end
    end
    acc
end

```

```julia
julia> @btime mysum(10000)
  21.943 μs (0 allocations: 0 bytes)
833666708335000

```

This comparison was “fair”, but quite meaningless. This is not where R shines, and this is well known; your only hope for fast runtime in R is when you find something called in C or Fortran. No one would use R for speed, only for the large selection of excellent packaged libraries. One can optimize the Julia code further I guess, but it would not make a huge difference to this conclusion.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 17, 2019, 5:27pm UTC](https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681/23 "2019-01-17T17:27:18Z")

</div>

> [@Tamas\_Papp](#):
>
> This comparison was “fair”, but quite meaningless. This is not where R shines, and this is well known

Just because it is widely known does not make it meaningless. It is still quite meaningful to quantify why you _can’t_ write performance-critical inner loops in R, and why you _can_ in Julia.

(Many people have a vague notion that “loops in R are slow” but have no idea how slow, or why, or whether loops written by mere mortals—not just the godlike beings who write the “built-in” libraries—can be fast in a different dynamic programming environment.)

But this is a terrible example to study multi-threading with for the reason @bennedich pointed out above: this kind of summation (over only 1000 elements!) is too cheap to parallelize effectively.

---

<div class="post-metadata">

**Author:** ![Juan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juan/32/7657_2.png) [@Juan](https://discourse.julialang.org/u/Juan)\
**Post date:** [January 17, 2019, 6:18pm UTC](https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681/24 "2019-01-17T18:18:56Z")

</div>

In fact I’m using the library data.table and sometimes is faster than Julia, for example to calculate the mean by groups or to read csv files.  
But yes, I know Julia is fast, that’s why I want to move to it.

---

<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:** [January 17, 2019, 6:41pm UTC](https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681/25 "2019-01-17T18:41:52Z")

</div>

> [@Juan](#):
>
> In fact I’m using the library data.table and sometimes is faster than Julia, for example to calculate the mean by groups or to read csv files.

See:

> [@Tamas\_Papp](#):
>
> your only hope for fast runtime in R is when you find something called in C or Fortran.

---

<div class="post-metadata">

**Author:** ![bennedich](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bennedich/32/4894_2.png) [@bennedich](https://discourse.julialang.org/u/bennedich)\
**Post date:** [January 17, 2019, 9:17pm UTC](https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681/26 "2019-01-17T21:17:22Z")

</div>

As @foobar_lv2 pointed out above, once you use a for loop or generator for the inner summation, the benchmark is flawed since the compiler will figure out that it can replace the loop with a formula. You can see this by either looking at the native code, or estimate if the elapsed time makes sense at all for the number of operations done, or simply grow the problem by a factor 10 and see how the time is affected. For this problem, it should be roughly quadratic, so if you see a linear growth, something is wrong.

To prevent this optimization from happening, you can instead use an array of numbers:

```julia
function my_sum(v)
    acc = zero(eltype(v))
    for i = 1:length(v), j = 1:i
        acc += v[j] * v[j]
    end
    acc
end

julia> Random.seed!(0); v = rand(100);

julia> @btime my_sum($v);
  5.776 μs (0 allocations: 0 bytes)

julia> Random.seed!(0); v = rand(1000);

julia> @btime my_sum($v);
  557.419 μs (0 allocations: 0 bytes)

```

Is this benchmark more accurate? Let’s quickly check:

- A 10x larger problem resulted in a ~100x run-time. Quadratic growth; this is what we’d expect.
- A 1000 element input vector should result in 1000\*1001/2 = 500500 additions. A time of 557 μs means ~1.1 ns per addition, or ~3.2 clock cycles on my system. This is very reasonable.

Now let’s see if we can parallelize this.

To reduce the overhead of threading, it usually makes sense to apply multithreading on as coarse of a level as possible. For this problem, we can split the input into one chunk per thread. So if the input size is 1000, and you have 8 threads available, split the input into 8 chunks of 125, let each thread finish its work, then sum up these 8 chunk-sums in the end.

To make all chunks do about the same amount of work in this problem, we can interleave the `i` indexes. This is important, since otherwise some threads will finish much faster than others, and will just idle while the remaining threads finish their work.

Below is a sample implementation.

```julia
function my_multithreaded_sum(v; T = Threads.nthreads())
    acc = zeros(eltype(v), T)
    Threads.@threads for t = 1:T
        s = zero(eltype(v))
        for i = t:T:length(v), j = 1:i
            s += v[j] * v[j]
        end
        acc[t] = s
    end
    return sum(acc)
end

```

Let’s see how efficient it is:

```julia
julia> Random.seed!(0); v = rand(1000);

julia> @btime my_sum($v) # original version
  557.422 μs (0 allocations: 0 bytes)
168854.26023687015
  
julia> @btime my_multithreaded_sum($v)
  75.522 μs (2 allocations: 192 bytes)
168854.2602369164

julia> Threads.nthreads()
8

julia> 557.422 / 75.522
7.380922115410079

```

(Note: slightly different sums due to rounding errors, this is expected)

---

<div class="post-metadata">

**Author:** ![davidbp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/davidbp/32/463_2.png) [@davidbp](https://discourse.julialang.org/u/davidbp)\
**Post date:** [January 19, 2019, 11:31pm UTC](https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681/27 "2019-01-19T23:31:33Z")

</div>

Thank you for your amazing answers!

In this case it seems that a single thread + simd is a perfect match

```julia
function my_simd_sum(v)
    s = zero(eltype(v))
    for i = 1:length(v)
        @simd for j = 1:i
        @inbounds s += v[j] * v[j]
        end
    end
    return s
end

```

```julia
@btime my_simd_sum($v); 40.130 μs
@btime my_multithreaded_sum($v; T=8); 58.376 μs (2 allocations: 192 bytes)

```

- Do you have any advice on how to implement a solution with Threads and SIMD?

Probably there is not enough computation in this example to make it worth it but, for the sake of knowledge, I would like to know how to approach a very simple problem mixing both @threads and @simd.

This is what I tried

```julia
function my_multithreaded_simd_sum(v; T = Threads.nthreads())
    partial_result_per_thread = zeros(eltype(v), T)
    Threads.@threads for t = 1:T
        s = zero(eltype(v))
        len = div(length(v), T)
        domain_per_thread = ((t-1)*len +1):t*len
        for i in domain_per_thread
           @simd for j in 1:i
            @inbounds s += v[j] * v[j]
            end
        end
        partial_result_per_thread[t] = s
    end
    return sum(partial_result_per_thread)
end

```

```julia
@btime my_multithreaded_simd_sum($v);
9.196 μs (2 allocations: 192 bytes)

```

---

<div class="post-metadata">

**Author:** ![bennedich](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bennedich/32/4894_2.png) [@bennedich](https://discourse.julialang.org/u/bennedich)\
**Post date:** [January 20, 2019, 12:22am UTC](https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681/28 "2019-01-20T00:22:00Z")

</div>

That looks alright to me, but see my note above about interleaving the `i` indexes. The point of this is to let each thread do the same amount of work. Let’s say you have 8 threads and a 1000 element vector. With your approach, thread 1 will be doing `125*126/2 = 7875` additions, while thread 8 does `(1000*1001-875*876)/2 = 117250` additions. This means that thread 1 (and 2, 3, …) will finish long before thread 8 and just sit and idle. On the contrary, by interleaving the `i` indices, all threads will do approximately the same number of additions (this is very problem specific though). On my system, this doubles the performance (also with SIMD).

Btw, the implementation above will also not work correctly if the vector size is not a multiple of the number of threads, e.g.:

```julia
julia> my_multithreaded_simd_sum([1 2 3])
0

```

To fix that, you could do something like this:

```julia
n = length(v)
domain_per_thread = 1+((t-1)*n÷T):t*n÷T

```

---

<div class="post-metadata">

**Author:** ![davidbp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/davidbp/32/463_2.png) [@davidbp](https://discourse.julialang.org/u/davidbp)\
**Post date:** [January 20, 2019, 4:31pm UTC](https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681/29 "2019-01-20T16:31:06Z")

</div>

> [@bennedich](#):
>
> On the contrary, by interleaving the `i` indices  
> …
> 
> With your approach, thread 1 will be doing `125*126/2 = 7875` additions, while thread 8 does `(1000*1001-875*876)/2 = 117250` additions.

I am not sure what you did. Could you post your function with the change you mention that is twice faster?

This type of exercices are helpful for learning a little bit how to use Threads.

---

<div class="post-metadata">

**Author:** ![bennedich](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bennedich/32/4894_2.png) [@bennedich](https://discourse.julialang.org/u/bennedich)\
**Post date:** [January 20, 2019, 4:38pm UTC](https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681/30 "2019-01-20T16:38:56Z")

</div>

Of course, it’s simply the code I posted in my previous post, but with `@simd` and `@inbounds`:

```julia
function my_multithreaded_sum_simd(v; T = Threads.nthreads())
    acc = zeros(eltype(v), T)
    Threads.@threads for t = 1:T
        s = zero(eltype(v))
        for i = t:T:length(v) # this is the "interleaving"
            @simd for j = 1:i
                @inbounds s += v[j] * v[j]
            end
        end
        acc[t] = s
    end
    return sum(acc)
end

```

---

<div class="post-metadata">

**Author:** ![KZiemian](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kziemian/32/9020_2.png) [@KZiemian](https://discourse.julialang.org/u/KZiemian)\
**Post date:** [January 28, 2019, 2:39pm UTC](https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681/31 "2019-01-28T14:39:19Z")

</div>

> [@foobar\_lv2](#):
>
> What happens is that llvm manages to optimize the inner loop into an explicit expression: `sum(j^p for j=1:N)` actually is a polynomial in `N` , llvm figures this out and then computes this in O(1).

Can you make this more clear to me? How compiler make something that is polynomial in `N` computable in O(1)? I probably misunderstood some concepts here.

---

<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:** [January 28, 2019, 2:51pm UTC](https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681/32 "2019-01-28T14:51:42Z")

</div>

Because of math.

> **[Summation | Powers and logarithm of arithmetic progressions](https://en.wikipedia.org/wiki/Summation#Powers_and_logarithm_of_arithmetic_progressions)**
>
> More generally, one has Faulhaber's formula for 
>   
>     
>       
> p
> \>
> 1
>       
>     
> {\\displaystyle p\>1}

---

<div class="post-metadata">

**Author:** ![KZiemian](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kziemian/32/9020_2.png) [@KZiemian](https://discourse.julialang.org/u/KZiemian)\
**Post date:** [January 28, 2019, 3:42pm UTC](https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681/33 "2019-01-28T15:42:39Z")

</div>

Excuse me, I need ask more about that. LLVM is able to recognize mathematical code, that have known formula and just use this formula? I’m ignorant in the world of compilers and after my C and C++ expirience, this looks unbelivebly good, so its hard to belive. To be sure, do I understand this correctly?

I don’t mean that compilers of C and C++ are bad, I just never noticing them doing something that. Maybe I just don’t look very cearfully.

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [January 28, 2019, 4:04pm UTC](https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681/34 "2019-01-28T16:04:59Z")

</div>

LLVM and compilers in general do know about basic integer arithmetic identities and can transform e.g. sum of a consecutive range of integers into the corresponding polynomial expression for the sum. It doesn’t go too much beyond that and doesn’t work for floating-point since floating-point isn’t purely associative like integers are.

---

<div class="post-metadata">

**Author:** ![yuyichao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yuyichao/32/20_2.png) [@yuyichao](https://discourse.julialang.org/u/yuyichao)\
**Post date:** [January 28, 2019, 7:40pm UTC](https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681/35 "2019-01-28T19:40:08Z")

</div>

> [@KZiemian](#):
>
> I just never noticing them doing something that.

They do. If you write a C function that add a integer sequence, they should do this optimization automatically there too (both clang ~~and gcc~~ (edit: aparently not for gcc) at least).

I actually kind of doubt how useful this is though. Unless the pattern is somehow really hiden from the programmer this should be a fairly simple transform the programmer can do. Anyway, most optimization are trivial by themselves anyway so I’m certainly not complaining that they could do this…

---

<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:** [January 28, 2019, 7:57pm UTC](https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681/36 "2019-01-28T19:57:14Z")

</div>

Clang definitely does, but [here](https://godbolt.org/z/26bJve) gcc seems not to.

---

<div class="post-metadata">

**Author:** ![Orbots](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/orbots/32/3392_2.png) [@Orbots](https://discourse.julialang.org/u/Orbots)\
**Post date:** [January 29, 2019, 12:55am UTC](https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681/37 "2019-01-29T00:55:52Z")

</div>

This looks nice. When I was using clojure I would rely on reducers quite heavily for performant code. This looks similar. It used java’s fork/join framework under the hood. I’d love to see a reducers.jl package.

[https://clojure.org/reference/reducers](https://clojure.org/reference/reducers)

---

<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:** [January 29, 2019, 2:30am UTC](https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681/38 "2019-01-29T02:30:47Z")

</div>

My [Transducers.jl](https://github.com/tkf/Transducers.jl) package (which obviously is inspired by Clojure) also does threading with `mapreduce`. But it’s very simple and not optimized ATM:

```julia
using Transducers
mapreduce(Map(i -> mapfoldl(MapSplat(*), +, zip(1:i, 1:i); simd=true)),
          +, 1:N; init=0)

```

---

<div class="post-metadata">

**Author:** ![KZiemian](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kziemian/32/9020_2.png) [@KZiemian](https://discourse.julialang.org/u/KZiemian)\
**Post date:** [January 29, 2019, 8:21pm UTC](https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681/39 "2019-01-29T20:21:39Z")

</div>

Thank you for that information. Now it is clear.

---

<div class="post-metadata">

**Author:** ![KZiemian](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kziemian/32/9020_2.png) [@KZiemian](https://discourse.julialang.org/u/KZiemian)\
**Post date:** [January 29, 2019, 8:22pm UTC](https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681/40 "2019-01-29T20:22:20Z")

</div>

I only used gcc, so this is maybe part of reason. Thank you for that information.

[Previous page](https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681.md?page=1)