# Multithreaded parallel prefix scan not faster

**URL:** <https://discourse.julialang.org/t/multithreaded-parallel-prefix-scan-not-faster/42591>\
**Category:** Performance\
**Tags:** question, performance, multithreading\
**Created:** [July 6, 2020, 7:06am UTC](https://discourse.julialang.org/t/multithreaded-parallel-prefix-scan-not-faster/42591 "2020-07-06T07:06:49Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ellipse0934](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ellipse0934/32/4012_2.png) [@Ellipse0934](https://discourse.julialang.org/u/Ellipse0934)\
**Post date:** [July 6, 2020, 7:06am UTC](https://discourse.julialang.org/t/multithreaded-parallel-prefix-scan-not-faster/42591/1 "2020-07-06T07:06:49Z")

</div>

I am trying to write the parallel prefix [scan](https://en.wikipedia.org/wiki/Prefix_sum) as a multithreaded CPU function. Since the number of processors is small relative to the length of the array It wouldn’t be fruitful to use Hillis-Steele or Belloch scan.

Instead for `n` processors my algorithm divides the array into `n` parts and then spawns a task to each thread. Storing each part’s sum to a different array and then performing a prefix scan on this array. Finally broadcasting back these sums to `n - 1` parts which need the sum of previous parts.

So, in theory for 2 processors the new algorithm should take \frac{t}{2} + \frac{t}{4} =\frac{3t}{4}.  
And for 6 processors (my case), \frac{t}{6} + \frac{5t}{36} = \frac{11t}{36}

```julia
function scan_naive!(op, output::Vector{T}, input::Vector{T}) where T <: Number
    output[1] = input[1]
    for i = 2:length(input)
        @inbounds output[i] = op(output[i - 1], input[i])
    end
    output
end

```

```julia
using Plots
sizes = 1 .<< collect(10:29)

function timer(f, op, len)
    A = ones(Int32, len)
    B = similar(A)
    x = Base.@elapsed f(op, B, A)
    y = Base.@elapsed f(op, B, A)
    z = Base.@elapsed f(op, B, A)
    min(x, y, z) # Unstable results without multiple trials for multithreaded version
end
t = [timer(scan_naive!, +, i) for i in sizes]
plot(sizes, t, title="Benchmark for prefix scan", scale=:log10, marker = :circle, legend=:bottomright)

```

![image](https://global.discourse-cdn.com/julialang/original/3X/2/7/27889d3017913202ac1ceea388ffd6635afad877.png)

```julia
function threaded_scan!(op, output::Vector{T}, input::Vector{T}) where T <: Number
    segment_length = length(input) ÷ nthreads()
    sums = Array{T}(undef, nthreads())
    @threads for i = 1:nthreads()
        low = 1 + (threadid() - 1)*segment_length
        high = threadid()*segment_length
        threadid() == nthreads() && (high = length(input))
        
        @inbounds output[low] = input[low]
        for j in (low + 1):high
            @inbounds output[j] = op(input[j], output[j - 1])
        end
        @inbounds sums[threadid()] = output[high]
    end
    scan_naive!(op, sums, sums)
    
    @threads for i in (segment_length + 1):length(input)
        segment = (i - 1) ÷ segment_length 
        i >= nthreads()*segment_length && (segment = nthreads() - 1)
        @inbounds output[i] = op(sums[segment], output[i])
    end
    output
end

```

```julia
t2 = [timer(threaded_scan!, +, i) for i in sizes]
plot(sizes, [t, t2], title="Benchmark for prefix scan", scale=:log10, marker = :circle, legend = :bottomright)

```

![image](https://global.discourse-cdn.com/julialang/original/3X/5/3/53620c8e99d479d3cf12d665f14a9830d50511d6.png)

This result was surprising to me as I expected the multithreaded version to overtake the single threaded version by the end which is an array of size 2^{28} \approx 5\cdot10^8.

Then I made a small optimization which is that I use only `n - 1` processors in the latter stage to broadcast partial sums which will remove the overhead of computing the segment and kick in SIMD instructions.

```julia
function threaded_scan2!(op, output::Vector{T}, input::Vector{T}) where T <: Number
    segment_length = length(input) ÷ nthreads()
    sums = Array{T}(undef, nthreads())
    @threads for i = 1:nthreads()
        low = 1 + (threadid() - 1)*segment_length
        high = threadid()*segment_length
        threadid() == nthreads() && (high = length(input))
        
        @inbounds output[low] = input[low]
        for j in (low + 1):high
            @inbounds output[j] = input[j] + output[j - 1]
        end
        @inbounds sums[threadid()] = output[high]
    end
    scan_naive!(op, sums, sums)
    
    tasks = []
    
    for i in 2:nthreads()
        low = 1 + (i- 1)*segment_length
        high = i*segment_length
        v = view(output, low:high)
        tsk = Threads.@spawn @inbounds v .+= sums[i - 1] # Only does addition, `op` is useless here 
        push!(tasks, tsk)
    end
    wait.(tasks)
    output
end

```

![image](https://global.discourse-cdn.com/julialang/original/3X/a/a/aae951c0172cbbdf8f8ceffcad7371d43055bfe4.png)

This appears to be better but still not beating the single threaded version. Since this is a memory intensive kernel maybe that’s the bottleneck that’s slowing this down, is there a way to confirm this suspicion ?

---

<div class="post-metadata">

**Author:** ![bernhard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bernhard/32/2619_2.png) [@bernhard](https://discourse.julialang.org/u/bernhard)\
**Post date:** [July 6, 2020, 10:39am UTC](https://discourse.julialang.org/t/multithreaded-parallel-prefix-scan-not-faster/42591/2 "2020-07-06T10:39:31Z")

</div>

You can take a look at the two macros `@btime and @benchmark`. From this package [GitHub - JuliaCI/BenchmarkTools.jl: A benchmarking framework for the Julia language](https://github.com/JuliaCI/BenchmarkTools.jl)  
The allocations (MB and count) and gc time should give you an idea about memory consumption (and time for garbage collection)

---

<div class="post-metadata">

**Author:** ![Ellipse0934](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ellipse0934/32/4012_2.png) [@Ellipse0934](https://discourse.julialang.org/u/Ellipse0934)\
**Post date:** [July 6, 2020, 11:16am UTC](https://discourse.julialang.org/t/multithreaded-parallel-prefix-scan-not-faster/42591/3 "2020-07-06T11:16:12Z")

</div>

I did that before posting, just didn’t mention it here. (2^{29} elements)

```julia
# threaded_scan2!
BenchmarkTools.Trial: 
  memory estimate: 8.83 KiB
  allocs estimate: 82
  --------------
  minimum time: 298.860 ms (0.00% GC)
  median time: 299.207 ms (0.00% GC)
  mean time: 299.275 ms (0.00% GC)
  maximum time: 300.185 ms (0.00% GC)
  --------------
  samples: 17
  evals/sample: 1

```

```julia
# naive_scan
BenchmarkTools.Trial: 
  memory estimate: 0 bytes
  allocs estimate: 0
  --------------
  minimum time: 192.783 ms (0.00% GC)
  median time: 193.240 ms (0.00% GC)
  mean time: 193.359 ms (0.00% GC)
  maximum time: 195.136 ms (0.00% GC)
  --------------
  samples: 26
  evals/sample: 1

```

I assumed that the threading is doing the allocations as internally I only declare 2 small arrays which shouldn’t be in the order of kilobytes. I even went ahead and changed `tasks = []` to an Array of undef Task with required size. This didn’t have much effect on the benchmark, only the allocs dropped to 76, 8.33 kB which isn’t significant.

---

<div class="post-metadata">

**Author:** ![eaubanel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eaubanel/32/13327_2.png) [@eaubanel](https://discourse.julialang.org/u/eaubanel)\
**Post date:** [July 7, 2020, 1:57pm UTC](https://discourse.julialang.org/t/multithreaded-parallel-prefix-scan-not-faster/42591/4 "2020-07-07T13:57:59Z")

</div>

> [@Ellipse0934](#):
>
> ```julia
> t2 = [timer(threaded_scan!, +, i) for i in sizes]
> plot(sizes, [t, t2], title="Benchmark for prefix scan", scale=:log10, marker = :circle, legend = :bottomright)
> 
> ```

What are you running on? On a 24-core Xeon Platinum 8260M I get the expected speedup for your first version with 16 threads:  
:  
 ![scanspeedupn16](https://global.discourse-cdn.com/julialang/original/3X/e/9/e9567041bdc246267d4292820fcff37809088880.png)

---

<div class="post-metadata">

**Author:** ![Ellipse0934](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ellipse0934/32/4012_2.png) [@Ellipse0934](https://discourse.julialang.org/u/Ellipse0934)\
**Post date:** [July 7, 2020, 4:55pm UTC](https://discourse.julialang.org/t/multithreaded-parallel-prefix-scan-not-faster/42591/5 "2020-07-07T16:55:58Z")

</div>

i5 8600k [link](https://ark.intel.com/content/www/us/en/ark/products/126685/intel-core-i5-8600k-processor-9m-cache-up-to-4-30-ghz.html). My benchmark was a bit noisy altering a bit from run to run. I did be a little dishonest and post on of the better runs here. I changed using the `@time` to using `@benchmark` along with a setup and things turned out much better. Visually I couldn’t notice any difference between different runs which was apparent earlier when I took the best time out of three runs.

But, even after that the multithreaded version was slower. I pulled out my system monitor and it appeared as if CPU cores would alternate between each other’s 100% usage. Since this is a **very** memory bound benchmark the OS’s policies on scheduling and swapping should play a significant part in the benchmark. I switched my OS to windows and it appears that in Windows multithreaded performance is better (sigh) here. Of course, the linux kernel version and the distro’s policies will com into play but this is interesting nonetheless.

On windows I get (with very few changes: @benchmark + using channel instead of Tasks)  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/f/6/f693d58a5e27c518661bd10e3f2812ab5083454e.png)

---

<div class="post-metadata">

**Author:** ![eaubanel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eaubanel/32/13327_2.png) [@eaubanel](https://discourse.julialang.org/u/eaubanel)\
**Post date:** [July 8, 2020, 1:12pm UTC](https://discourse.julialang.org/t/multithreaded-parallel-prefix-scan-not-faster/42591/6 "2020-07-08T13:12:45Z")

</div>

Oops, I see that I posted results with 16 threads above! Here are my results with 6:  
 ![scanspeedupn6](https://global.discourse-cdn.com/julialang/original/3X/e/2/e22dffa206bba58c61c18737d14a1e7198dba11f.png)  
So I only get a modest speedup for your first version, but the improved version gets the expected speedup. As you mention, the first version suffers from extra overhead in the second loop. In your case it’s likely memory pressure since you’re using all 6 cores, whereas I’m only using 6 out of 24.