# Multithreading problem

**URL:** <https://discourse.julialang.org/t/multithreading-problem/70357>\
**Category:** Performance\
**Tags:** multithreading\
**Created:** [October 25, 2021, 10:25pm UTC](https://discourse.julialang.org/t/multithreading-problem/70357 "2021-10-25T22:25:39Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![Jake](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jake/32/46007_2.png) [@Jake](https://discourse.julialang.org/u/Jake)\
**Post date:** [October 25, 2021, 10:25pm UTC](https://discourse.julialang.org/t/multithreading-problem/70357/1 "2021-10-25T22:25:39Z")

</div>

This is my first time trying to use multithreading for real work. I have an inner loop in my code for doing sinc interpolation of a recorded time history file. A minimum working example is

```julia
# Minumum working example for threaded operation of a for loop

using Base.Threads
using Random

@inline h(x) = (cos(π/2 * x/n)) ^2 # Hanning window where n is the number of data points in each direction
@inline wsinc(x) = sinc(x) * h(x) # windowed sinc function

n = 10
ordinate_time = rand(31)
fpart, ipart = modf(13.78)
ipart = Integer(ipart)
ya = 0.0

for (i, filtervalue) in enumerate(-n+fpart : n)
    index = ipart+1-n+i
    ya += ordinate_time[index] * wsinc(filtervalue)
end

ya

```

This gives a result. When I try multithreading I replace the for loop with

```julia
@threads for (i, filtervalue) in enumerate(-n+fpart : n)
    index = ipart+1-n+i
    #ya += ordinate_time[index] * wsinc(filtervalue)
    atomic_add!(ya, ordinate_time[index] * wsinc(filtervalue))
end

```

I don’t think I need the atomic\_add, but I tried it because I really don’t know what I am doing wrong here.

The error is

```julia
    nested task error: MethodError: no method matching firstindex(::Base.Iterators.Enumerate{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}})

```

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [October 26, 2021, 4:26am UTC](https://discourse.julialang.org/t/multithreading-problem/70357/2 "2021-10-26T04:26:50Z")

</div>

This is [a known issue](https://github.com/JuliaLang/julia/pull/40772), but I don’t think it’s solved yet

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [October 26, 2021, 4:49am UTC](https://discourse.julialang.org/t/multithreading-problem/70357/3 "2021-10-26T04:49:07Z")

</div>

The problem is `enumerate`, which (like many of Julia’s `Base.Iterators` constructs) doesn’t currently support indexing. You can work around it with something like this:

```julia
nrange = -n+fpart : n
Threads.@threads for i in eachindex(nrange)
    filtervalue = nrange[i]
    [...]
end

```

External libraries like [ThreadsX](https://github.com/tkf/ThreadsX.jl) and [LoopVectorization](https://github.com/JuliaSIMD/LoopVectorization.jl) tend to have better support for generic threading constructs at the moment, e.g.

```julia
using ThreadsX

ya = ThreadsX.sum(enumerate(-n+fpart : n)) do (i, filtervalue)
    index = ipart + 1 - n + i
    return ordinate_time[index] * wsinc(filtervalue)
end

```

---

<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:** [October 26, 2021, 5:17am UTC](https://discourse.julialang.org/t/multithreading-problem/70357/4 "2021-10-26T05:17:06Z")

</div>

FYI, I wrote [A quick introduction to data parallelism in Julia](https://juliafolds.github.io/data-parallelism/tutorials/quick-introduction/) for covering patterns like this.

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [October 26, 2021, 6:33am UTC](https://discourse.julialang.org/t/multithreading-problem/70357/5 "2021-10-26T06:33:29Z")

</div>

A separate performance quibble (maybe not reflected in your actual code, but noted just in case): `h(x)` depends on the global variable `n`, which [kills performance](https://docs.julialang.org/en/v1/manual/performance-tips/#Avoid-global-variables). Better to specify it as `h(x, n)`.

You might want to try [Lanczos interpolation](http://juliamath.github.io/Interpolations.jl/latest/api/#Interpolations.Lanczos) from Interpolations.jl (it’s just a sinc function windowed with a sinc).

```julia
support = 10
scale = 4
ya = interpolate(ordinate_time, Interpolations.Lanczos{support}(scale))[13.78]

```

---

<div class="post-metadata">

**Author:** ![Jake](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jake/32/46007_2.png) [@Jake](https://discourse.julialang.org/u/Jake)\
**Post date:** [October 27, 2021, 1:58am UTC](https://discourse.julialang.org/t/multithreading-problem/70357/6 "2021-10-27T01:58:36Z")

</div>

Thank you for pointing me to the ThreadsX library, and showing the code snippet. I think this is the first time I have seen the CPU at 100% on my computer. And also for mentioning the performance quibble. I was looking at that and wondering if I should change it, now I have. And thanks to tkf and all those that have made Julia what it is. I will check out the Lanczos interpolation.

Should I also implement the fastmath library for additional speed?

---

<div class="post-metadata">

**Author:** ![Jake](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jake/32/46007_2.png) [@Jake](https://discourse.julialang.org/u/Jake)\
**Post date:** [October 27, 2021, 8:11pm UTC](https://discourse.julialang.org/t/multithreading-problem/70357/7 "2021-10-27T20:11:28Z")

</div>

Doing some more experimentation I find that any attempt to parallelize this loop results in a loss of performance. Would I be correct in assuming that this is because the number of elements to be summed is too small (20 elements)?

My minimum working example now is:

```julia
# Minumum working example for threaded operation of a for loop
using ThreadsX, Base.Threads
using Random

function mwe()
@inline h(x,n) = (cos(π/2 * x/n)) ^2 # Hanning window where n is the number of data points in each direction
@inline wsinc(x,n) = sinc(x) * h(x,n) # windowed sinc function

n = 10
for i = 1:1000000
    ordinate_time = rand(31)
    fpart, ipart = modf(13.78)
    ipart = Integer(ipart)
    #=
    ya = 0.0
    nrange = -n+fpart : n
    for i in eachindex(nrange)
        filtervalue = nrange[i]
        index = ipart+1-n+i
        ya += ordinate_time[index] * wsinc(filtervalue,n)
    end
    =#
    
    ya = ThreadsX.sum(enumerate(-n+fpart : n)) do (i, filtervalue)
        index = ipart+1-n+i
        return ordinate_time[index] * wsinc(filtervalue,n)
    end
end
end

```

where I comment out the respective blocks of code to try each technique

The results for ThreadsX are:

```julia
ya = ThreadsX.sum(enumerate(-n+fpart : n)) do (i, filtervalue)
        index = ipart+1-n+i
        return ordinate_time[index] * wsinc(filtervalue,n)
    end

julia> Threads.nthreads()
12
julia> @btime mwe()
  50.776 s (307802549 allocations: 18.50 GiB)

julia> Threads.nthreads()
1
julia> @btime mwe()
  15.994 s (193000000 allocations: 8.43 GiB)

Use sum instead of ThreadsX.sum and get the following
julia> @btime mwe()
  3.425 s (82000000 allocations: 1.52 GiB)

```

Clearly the most performant code is without any parallelizing

Now try using the @threads macro with a for loop:

```julia
    Threads.@threads for i in eachindex(nrange)
        filtervalue = nrange[i]
        index = ipart+1-n+i
        ya += ordinate_time[index] * wsinc(filtervalue,n)
    end

julia> Threads.nthreads()
12
julia> @btime mwe()
  18.827 s (145337991 allocations: 7.51 GiB)

julia> Threads.nthreads()
1
julia> @btime mwe()
  4.828 s (90000000 allocations: 2.18 GiB)

Remove the Threads.threads before the for loop
julia> @btime mwe()
  716.158 ms (1000000 allocations: 320.43 MiB)

```

Clearly the for loop outperforms the sum function and adding parallelizing reduces performance for all cases.

I have an outer for loop that I may be able to use LoopVectorization on. I will try that next.

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [October 27, 2021, 9:10pm UTC](https://discourse.julialang.org/t/multithreading-problem/70357/8 "2021-10-27T21:10:25Z")

</div>

It’s almost always fastest to apply multithreading to the outer loop. Can you illustrate your outer loop in a MWE?

---

<div class="post-metadata">

**Author:** ![Jake](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jake/32/46007_2.png) [@Jake](https://discourse.julialang.org/u/Jake)\
**Post date:** [October 27, 2021, 9:32pm UTC](https://discourse.julialang.org/t/multithreading-problem/70357/9 "2021-10-27T21:32:04Z")

</div>

Here is a MWE.

```julia
# Minimum working example for threaded operation of a for loop
# parallel summation of the inner loop is too small for this to be effective
using ThreadsX, Base.Threads
using Random

function mwe()
@inline h(x,n) = (cos(π/2 * x/n)) ^2 # Hanning window where n is the number of data points in each direction
@inline wsinc(x,n) = sinc(x) * h(x,n) # windowed sinc function

t = 1:1000000
ordinate_time = sinpi.(2*10*t) .+ cospi.(2*14.7*t)
ordinate_angle = zeros(Float64, length(ordinate_time))
n = 10
for i = 15:1000000-15 # so I don't deal with the end conditions 
    fpart, ipart = modf(13.78)
    ipart = Integer(ipart)

    ya = 0.0
    nrange = -n+fpart : n
    for i in eachindex(nrange)
        filtervalue = nrange[i]
        index = ipart+1-n+i
        ya += ordinate_time[index] * wsinc(filtervalue,n)
    end
    ordinate_angle[i] = ya
    #=
    ya = ThreadsX.sum(enumerate(-n+fpart : n)) do (i, filtervalue)
        index = ipart+1-n+i
        return ordinate_time[index] * wsinc(filtervalue,n)
    end=#
end
return ordinate_angle
end

```

Edit: I added

```julia
using LoopVectorization 

@turbo for i = 15:1000000-15 # so I don't deal with the end conditions 

...

```

and now get the error

```julia
┌ Error: Failed to revise C:\Users\jakez\.julia\dev\SignalProcessing\example\MEW parallel sum.jl
│ exception =
│ LoadError: BoundsError: attempt to access 0-element Vector{LoopVectorization.Operation} at index [1]
│ in expression starting at C:\Users\jakez\.julia\dev\SignalProcessing\example\MEW parallel sum.jl:14
└ @ Revise C:\Users\jakez\.julia\packages\Revise\3RMhb\src\packagedef.jl:716

```

I am not sure what is going on.

---

<div class="post-metadata">

**Author:** ![kfrb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kfrb/32/29557_2.png) [@kfrb](https://discourse.julialang.org/u/kfrb)\
**Post date:** [November 3, 2021, 8:58am UTC](https://discourse.julialang.org/t/multithreading-problem/70357/10 "2021-11-03T08:58:08Z")

</div>

After reworking your MWE, here’s what I get:

```julia
using Base.Threads, Random, BenchmarkTools, Test

@inline h(x::Float64,n::Int) = (cos(π/2 * x/n))^2
@inline wsinc(x::Float64,n::Int) = sinc(x) * h(x,n)

@inline function calc_ya(n::Int, ipart::Int ,fpart::Float64, ordinate_time::Vector{Float64})
    ya = 0.0
    for (j, filtervalue) in enumerate(-n+fpart : n)
        index = ipart+1-n+j
        ya += ordinate_time[index] * wsinc(filtervalue,n)
    end
    return ya
end

function mwe()
    t = 1:1000000
    n = 10
    ordinate_time = sinpi.(2*10*t) .+ cospi.(2*14.7*t)
    ordinate_angle = zeros(Float64, length(ordinate_time))
    fpart, ipart = modf(13.78)
    ipart = Integer(ipart)
    tsteps = 15:length(t)-15
    for i ∈ tsteps
        ordinate_angle[i] = calc_ya(n,ipart,fpart,ordinate_time)
    end
    return ordinate_angle
end

function mwe_threads()
    t = 1:1000000
    n = 10
    ordinate_time = sinpi.(2*10*t) .+ cospi.(2*14.7*t)
    ordinate_angle = zeros(Float64, length(ordinate_time))
    fpart, ipart = modf(13.78)
    ipart = Integer(ipart)
    tsteps = 15:length(t)-15
    @threads for i ∈ tsteps
        ordinate_angle[i] = calc_ya(n,ipart,fpart,ordinate_time)
    end
    return ordinate_angle
end

```

With 8 threads I get the following performance:

```julia-repl
julia> @btime mwe();
  588.646 ms (4 allocations: 15.26 MiB)

julia> @btime mwe_threads();
  108.995 ms (47 allocations: 15.26 MiB)

julia> @test mwe() == mwe_threads()
Test Passed

```

I hope this helps you and I have understood your MWE correctly. Note that all variables that do not depend on `i` can be calculated outside the for-loop.

---

<div class="post-metadata">

**Author:** ![Jake](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jake/32/46007_2.png) [@Jake](https://discourse.julialang.org/u/Jake)\
**Post date:** [November 4, 2021, 12:03pm UTC](https://discourse.julialang.org/t/multithreading-problem/70357/11 "2021-11-04T12:03:13Z")

</div>

Thank you kfrb for this reworked example. I was able to take your modification to the MWE and incorporate it into the main code.

The reason for some of the variables that were within the for loop were there because in the main code they did depend upon i, but thanks for the tip.

I think the main take-away from your rework is to put the inner for loop into its own function.