# Multithreading workflow - alternative for :static scheduling

**URL:** <https://discourse.julialang.org/t/multithreading-workflow-alternative-for-static-scheduling/87271>\
**Category:** General Usage\
**Tags:** multithreading\
**Created:** [September 14, 2022, 12:06pm UTC](https://discourse.julialang.org/t/multithreading-workflow-alternative-for-static-scheduling/87271 "2022-09-14T12:06:20Z")\
**Posts on this page:** 3\
**Page:** 1

<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:** [September 14, 2022, 12:06pm UTC](https://discourse.julialang.org/t/multithreading-workflow-alternative-for-static-scheduling/87271/1 "2022-09-14T12:06:20Z")

</div>

As stated in the docs [(Documentation for Threads.@threads)](https://github.com/JuliaLang/julia/blob/afb6c60d69a38e8a2442a0c7e87c47b8880ad294/base/threadingconstructs.jl#L99-L205):

“_`:static` scheduling exists for supporting transition of code written before Julia 1.3.  
In newly written library functions, `:static` scheduling is discouraged because the  
functions using this option cannot be called from arbitrary worker threads._”

To avoid racing conditions, I usually divide the work manually between workers and then each thread writes to an own column of a result matrix. With this workflow, `:static` scheduling is needed to obtain correct results.

Consider the following example:

```julia
using Base.Threads

function scheduling_comparison()
    n = 100
    x = 0.01 .* collect(1:n)
    indices = rand(1:n, n*1000)
    owned_indices = defaultdist(length(indices), nthreads())
    result_serial = serial_solution(x, indices, n)
    result_static = static_scheduling(x, indices, owned_indices, n)
    result_dynamic = dynamic_scheduling(x, indices, owned_indices, n)
    @show result_serial ≈ result_static
    @show result_serial ≈ result_dynamic
    return nothing
end

some_big_calculation(x) = sin(x) - exp(x) * x^2

function serial_solution(x, indices, n)
    result = zeros(n)
    for index in indices
        result[index] += some_big_calculation(x[index])
    end
    return result
end

function static_scheduling(x, indices, owned_indices, n)
    result_matrix = zeros(n, nthreads())
    @threads :static for _ in 1:nthreads()
        for index_pos in owned_indices[threadid()]
            index = indices[index_pos]
            result_matrix[index, threadid()] += some_big_calculation(x[index])
        end
    end
    result = zeros(n)
    @threads for i in 1:n
        result[i] = sum(@view result_matrix[i, :])
    end
    return result
end

function dynamic_scheduling(x, indices, owned_indices, n)
    result_matrix = zeros(n, nthreads())
    @threads :dynamic for _ in 1:nthreads()
        for index_pos in owned_indices[threadid()]
            index = indices[index_pos]
            result_matrix[index, threadid()] += some_big_calculation(x[index])
        end
    end
    result = zeros(n)
    @threads for i in 1:n
        result[i] = sum(@view result_matrix[i, :])
    end
    return result
end

function defaultdist(sz::Int, nc::Int)
    if sz >= nc
        chunk_size = div(sz,nc)
        remainder = rem(sz,nc)
        sidx = zeros(Int64, nc+1)
        for i = 1:(nc+1)
            sidx[i] += (i-1)*chunk_size + 1
            if i<= remainder
                sidx[i] += i-1
            else
                sidx[i] += remainder
            end
        end
        grid = fill(0:0,nc)
        for i = 1:nc
            grid[i] = sidx[i]:sidx[i+1]-1
        end
        return grid
    else
        sidx = [1:(sz+1);]
        grid = fill(0:0,nc)
        for i = 1:sz
            grid[i] = sidx[i]:sidx[i+1]-1
        end
        return grid
    end
end

```

As you can see, for `:dynamic` scheduling sometimes wrong results are determined:

```julia-repl
julia> scheduling_comparison()
result_serial ≈ result_static = true
result_serial ≈ result_dynamic = false

julia> scheduling_comparison()
result_serial ≈ result_static = true
result_serial ≈ result_dynamic = true

```

The problem with `:dynamic` scheduling is that the value of [`threadid()`](https://docs.julialang.org/en/v1/base/multi-threading/#Base.Threads.threadid) is not guranteed to be constant within one iteration and thus certain indices are not used for the calculation.

Since I am in the process of preparing my first Julia package:  
What is the official way if it is important that [`threadid()`](https://docs.julialang.org/en/v1/base/multi-threading/#Base.Threads.threadid) has a constant value per iteration?

Any suggestions, comments or help with my workflow would be greatly appreciated!

---

<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:** [September 14, 2022, 4:53pm UTC](https://discourse.julialang.org/t/multithreading-workflow-alternative-for-static-scheduling/87271/2 "2022-09-14T16:53:55Z")

</div>

```julia
@threads :dynamic for tid in 1:nthreads()
        for index_pos in owned_indices[tid]
            index = indices[index_pos]
            result_matrix[index, tid] += some_big_calculation(x[index])
        end
    end

```

---

<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:** [September 15, 2022, 9:15am UTC](https://discourse.julialang.org/t/multithreading-workflow-alternative-for-static-scheduling/87271/3 "2022-09-15T09:15:35Z")

</div>

Well, sometimes you do not see the wood for the trees…  
Thank you for your answer, @Elrod!

In case somebody is interested, here the performance differences (8 threads, Julia 1.8.1, Intel i9 on macOS Monterey):

```julia
function scheduling_benchmarks()
    n = 100000
    x = 0.01 .* collect(1:n)
    indices = rand(1:n, n*1000)
    owned_indices = defaultdist(length(indices), nthreads())
    @btime result_serial = serial_solution($x, $indices, $n)
    @btime result_static = static_scheduling($x, $indices, $owned_indices, $n)
    @btime result_dynamic = dynamic_scheduling($x, $indices, $owned_indices, $n)
    return nothing
end

```

```julia-repl
julia> scheduling_benchmarks()
  3.058 s (2 allocations: 781.30 KiB)
  532.266 ms (102 allocations: 6.87 MiB)
  540.981 ms (103 allocations: 6.87 MiB)

```

> **Complete code**
>
> ```julia
> using Base.Threads
> using BenchmarkTools
> 
> function scheduling_comparison()
> n = 100
> x = 0.01 .* collect(1:n)
> indices = rand(1:n, n*1000)
> owned_indices = defaultdist(length(indices), nthreads())
> result_serial = serial_solution(x, indices, n)
> result_static = static_scheduling(x, indices, owned_indices, n)
> result_dynamic = dynamic_scheduling(x, indices, owned_indices, n)
> @show result_serial ≈ result_static
> @show result_serial ≈ result_dynamic
> return nothing
> end
> 
> function scheduling_benchmarks()
> n = 100000
> x = 0.01 .* collect(1:n)
> indices = rand(1:n, n*1000)
> owned_indices = defaultdist(length(indices), nthreads())
> @btime result_serial = serial_solution($x, $indices, $n)
> @btime result_static = static_scheduling($x, $indices, $owned_indices, $n)
> @btime result_dynamic = dynamic_scheduling($x, $indices, $owned_indices, $n)
> return nothing
> end
> 
> some_big_calculation(x) = sin(x) - exp(x) * x^2
> 
> function serial_solution(x, indices, n)
> result = zeros(n)
> for index in indices
> result[index] += some_big_calculation(x[index])
> end
> return result
> end
> 
> function dynamic_scheduling(x, indices, owned_indices, n)
> result_matrix = zeros(n, nthreads())
> @threads :dynamic for tid in 1:nthreads()
> for index_pos in owned_indices[tid]
> index = indices[index_pos]
> result_matrix[index, tid] += some_big_calculation(x[index])
> end
> end
> result = zeros(n)
> @threads for i in 1:n
> result[i] = sum(@view result_matrix[i, :])
> end
> return result
> end
> 
> function static_scheduling(x, indices, owned_indices, n)
> result_matrix = zeros(n, nthreads())
> @threads :static for tid in 1:nthreads()
> for index_pos in owned_indices[tid]
> index = indices[index_pos]
> result_matrix[index, tid] += some_big_calculation(x[index])
> end
> end
> result = zeros(n)
> @threads for i in 1:n
> result[i] = sum(@view result_matrix[i, :])
> end
> return result
> end
> 
> function defaultdist(sz::Int, nc::Int)
> if sz >= nc
> chunk_size = div(sz,nc)
> remainder = rem(sz,nc)
> sidx = zeros(Int64, nc+1)
> for i = 1:(nc+1)
> sidx[i] += (i-1)*chunk_size + 1
> if i<= remainder
> sidx[i] += i-1
> else
> sidx[i] += remainder
> end
> end
> grid = fill(0:0,nc)
> for i = 1:nc
> grid[i] = sidx[i]:sidx[i+1]-1
> end
> return grid
> else
> sidx = [1:(sz+1);]
> grid = fill(0:0,nc)
> for i = 1:sz
> grid[i] = sidx[i]:sidx[i+1]-1
> end
> return grid
> end
> end
> 
> ```