# Multithreading and generated code

**URL:** <https://discourse.julialang.org/t/multithreading-and-generated-code/60828>\
**Category:** Performance\
**Tags:** multithreading\
**Created:** [May 9, 2021, 4:04pm UTC](https://discourse.julialang.org/t/multithreading-and-generated-code/60828 "2021-05-09T16:04:11Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![gideonsimpson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gideonsimpson/32/1928_2.png) [@gideonsimpson](https://discourse.julialang.org/u/gideonsimpson)\
**Post date:** [May 9, 2021, 4:04pm UTC](https://discourse.julialang.org/t/multithreading-and-generated-code/60828/1 "2021-05-09T16:04:11Z")

</div>

A while back, I solved this problem, [Performance and Allocation with Arrays of Functions](https://discourse.julialang.org/t/performance-and-allocation-issue-with-arrays-of-functions-v1-5/45013/7). More recently, I wanted a multithreaded loop inside of this, so that the code now reads:

```julia
function update(x, Δt)
    return x + 0.5 * Δt * x
end

@generated function compute_threaded(x0_vals, Δt, nΔt, f::Tuple{Vararg{<:Any,K}}) where {K}
    quote
        x_vals = copy(x0_vals);
        
        nx = length(x0_vals);
        f_vals = zeros($K, nx);
        f_avg_vals = zeros($K, nΔt)
        
        for n in 1:nΔt
            # @. x_vals = update(x_vals, Δt); switch to this, and it works fine
            Threads.@threads for j in 1:nx
                x_vals[j] = update(x_vals[j], Δt)
            end
            Base.Cartesian.@nexprs $K k -> f_vals[k,:] .= (f[k]).(x_vals);
            for j in 1:nx
                Base.Cartesian.@nexprs $K k -> f_avg_vals[k,n] += f_vals[k,j]/nx
            end
        end
        return f_avg_vals
    end
end

Random.seed!(100);
x0 = randn(10);
x₀ = [1.0];
Δt = 0.5;
nΔt = 10^2;

f1(x) = sin(x)
f2(x) = x^2

compute_threaded(x0, Δt, nΔt, (f1,))

```

I get the error:

```julia
The function body AST defined by this @generated function is not pure. This likely means it contains a closure or comprehension.

```

The motivation here is that the **real** `update` function will be very expensive and will benefit from multithreading. I am open to moving away from `@generated` code, but the motivation from the earlier post still remains: I want to evaluate some collection of observable functions, on the fly, on a time series, as it is generated, and I may have a variable number of such observables for different problems.

---

<div class="post-metadata">

**Author:** ![simeonschaub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simeonschaub/32/216566_2.png) [@simeonschaub](https://discourse.julialang.org/u/simeonschaub)\
**Post date:** [May 9, 2021, 4:28pm UTC](https://discourse.julialang.org/t/multithreading-and-generated-code/60828/2 "2021-05-09T16:28:35Z")

</div>

There is really no reason to use `@nexprs` together with a generated function here. You can just use `ntuple`, e.g.

```julia
ntuple(k -> f_vals[k,:] .= (f[k]).(x_vals), K)

```

instead of

```julia
Base.Cartesian.@nexprs $K k -> f_vals[k,:] .= (f[k]).(x_vals)

```

so you don’t need a generated function.

---

<div class="post-metadata">

**Author:** ![gideonsimpson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gideonsimpson/32/1928_2.png) [@gideonsimpson](https://discourse.julialang.org/u/gideonsimpson)\
**Post date:** [May 9, 2021, 4:54pm UTC](https://discourse.julialang.org/t/multithreading-and-generated-code/60828/3 "2021-05-09T16:54:05Z")

</div>

There’s a performance hit with that. Consider the following simplified example:

```julia
function update(x, Δt)
    return x + 0.5 * Δt * x
end

@generated function compute_values(x₀, Δt, nΔt, f::Tuple{Vararg{<:Any,K}}) where {K}
    quote
        x = x₀;
        f_vals = zeros($K, nΔt)
        for n in 1:nΔt
            x =update(x,Δt)
            Base.Cartesian.@nexprs $K k -> f_vals[k,n] = (f[k])(x);
        end
        return f_vals
    end
end

function compute_values_ntuple(x₀, Δt, nΔt, f::Tuple{Vararg{<:Any,K}}) where {K}
    x = x₀;
    f_vals = zeros(K, nΔt)
    for n in 1:nΔt
        x =update(x,Δt)
        ntuple(k -> f_vals[k,n] = (f[k])(x), K)
    end
    return f_vals
end

x₀ = 1.0;
Δt = 0.5;
nΔt = 10^3;

f1(x) = sin(x)
f2(x) = x^2
@btime compute_values($x₀, $Δt, $nΔt, $(f1,))
@btime compute_values($x₀, $Δt, $nΔt, $(f2,))
@btime compute_values($x₀, $Δt, $nΔt, $(f1,f2))

  36.936 μs (1 allocation: 7.94 KiB)
  3.044 μs (1 allocation: 7.94 KiB)
  37.929 μs (1 allocation: 15.75 KiB)

@btime compute_values_ntuple($x₀, $Δt, $nΔt, $(f1,))
@btime compute_values_ntuple($x₀, $Δt, $nΔt, $(f2,))
@btime compute_values_ntuple($x₀, $Δt, $nΔt, $(f1,f2))

 230.513 μs (3003 allocations: 54.84 KiB)
  189.271 μs (3003 allocations: 54.84 KiB)
  458.585 μs (4003 allocations: 78.28 KiB)

```

It’s about an order of magnitude slower, and the allocations now scale linearly with `nΔt`

---

<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:** [May 9, 2021, 5:08pm UTC](https://discourse.julialang.org/t/multithreading-and-generated-code/60828/4 "2021-05-09T17:08:48Z")

</div>

Try recursion.

```julia
@inline recursive_assign!(f::Tuple{}, f_vals, k, n, x) = nothing
@inline function recursive_assign!(f::Tuple{F,Vararg}, f_vals, x, k, n) where {F}
    f_vals[k,n] = first(f)(x)
    recursive_assign!(Base.tail(f), f_vals, x, k+1, n)
end

function compute_values_recurse(x₀, Δt, nΔt, f::Tuple{Vararg{<:Any,K}}) where {K}
    x = x₀;
    f_vals = zeros(K, nΔt)
    for n in 1:nΔt
        x =update(x,Δt)
        recursive_assign!(f, f_vals, x, 1, n) 
    end
    return f_vals
end

```

I get:

```julia
julia> @btime compute_values_recurse($x₀, $Δt, $nΔt, $(f1,))
  24.771 μs (1 allocation: 7.94 KiB)
1×1000 Matrix{Float64}:
 0.948985 0.999966 0.927798 0.64436 0.0897141 -0.623416 -0.998433 … -0.395455 -0.376427 -0.298363 -0.929138 0.982109 -0.998104 0.891708

julia> @btime compute_values_recurse($x₀, $Δt, $nΔt, $(f2,))
  2.277 μs (1 allocation: 7.94 KiB)
1×1000 Matrix{Float64}:
 1.5625 2.44141 3.8147 5.96046 9.31323 14.5519 22.7374 35.5271 55.5112 … 1.10853e193 1.73207e193 2.70636e193 4.22869e193 6.60733e193

julia> @btime compute_values_recurse($x₀, $Δt, $nΔt, $(f1, f2,))
  25.134 μs (1 allocation: 15.75 KiB)
2×1000 Matrix{Float64}:
 0.948985 0.999966 0.927798 0.64436 0.0897141 -0.623416 -0.998433 -0.317148 … -0.298363 -0.929138 0.982109 -0.998104 0.891708
 1.5625 2.44141 3.8147 5.96046 9.31323 14.5519 22.7374 35.5271 1.10853e193 1.73207e193 2.70636e193 4.22869e193 6.60733e193

julia> @btime compute_values($x₀, $Δt, $nΔt, $(f1,))
  24.643 μs (1 allocation: 7.94 KiB)
1×1000 Matrix{Float64}:
 0.948985 0.999966 0.927798 0.64436 0.0897141 -0.623416 -0.998433 … -0.395455 -0.376427 -0.298363 -0.929138 0.982109 -0.998104 0.891708

julia> @btime compute_values($x₀, $Δt, $nΔt, $(f2,))
  2.276 μs (1 allocation: 7.94 KiB)
1×1000 Matrix{Float64}:
 1.5625 2.44141 3.8147 5.96046 9.31323 14.5519 22.7374 35.5271 55.5112 … 1.10853e193 1.73207e193 2.70636e193 4.22869e193 6.60733e193

julia> @btime compute_values($x₀, $Δt, $nΔt, $(f1,f2))
  25.263 μs (1 allocation: 15.75 KiB)
2×1000 Matrix{Float64}:
 0.948985 0.999966 0.927798 0.64436 0.0897141 -0.623416 -0.998433 -0.317148 … -0.298363 -0.929138 0.982109 -0.998104 0.891708
 1.5625 2.44141 3.8147 5.96046 9.31323 14.5519 22.7374 35.5271 1.10853e193 1.73207e193 2.70636e193 4.22869e193 6.60733e193

julia> @btime compute_values_ntuple($x₀, $Δt, $nΔt, $(f1,))
  156.152 μs (3003 allocations: 54.84 KiB)
1×1000 Matrix{Float64}:
 0.948985 0.999966 0.927798 0.64436 0.0897141 -0.623416 -0.998433 … -0.395455 -0.376427 -0.298363 -0.929138 0.982109 -0.998104 0.891708

julia> @btime compute_values_ntuple($x₀, $Δt, $nΔt, $(f2,))
  131.799 μs (3003 allocations: 54.84 KiB)
1×1000 Matrix{Float64}:
 1.5625 2.44141 3.8147 5.96046 9.31323 14.5519 22.7374 35.5271 55.5112 … 1.10853e193 1.73207e193 2.70636e193 4.22869e193 6.60733e193

julia> @btime compute_values_ntuple($x₀, $Δt, $nΔt, $(f1,f2))
  298.434 μs (4003 allocations: 78.28 KiB)
2×1000 Matrix{Float64}:
 0.948985 0.999966 0.927798 0.64436 0.0897141 -0.623416 -0.998433 -0.317148 … -0.298363 -0.929138 0.982109 -0.998104 0.891708
 1.5625 2.44141 3.8147 5.96046 9.31323 14.5519 22.7374 35.5271 1.10853e193 1.73207e193 2.70636e193 4.22869e193 6.60733e193

```

---

<div class="post-metadata">

**Author:** ![simeonschaub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simeonschaub/32/216566_2.png) [@simeonschaub](https://discourse.julialang.org/u/simeonschaub)\
**Post date:** [May 9, 2021, 5:24pm UTC](https://discourse.julialang.org/t/multithreading-and-generated-code/60828/5 "2021-05-09T17:24:10Z")

</div>

> [@gideonsimpson](#):
>
> There’s a performance hit with that. Consider the following simplified example:

Ah, that’s because `x` is assigned to multiple times so it gets boxed in the closure. This can be circumvented by using a `let` block:

```julia
function compute_values_ntuple(x₀, Δt, nΔt, f::NTuple{K}) where {K}
    x = x₀;
    f_vals = zeros(K, nΔt)
    for n in 1:nΔt
        x =update(x,Δt)
        let x=x
            ntuple(k -> f_vals[k,n] = (f[k])(x), K)
        end
    end
    return f_vals
end

```

With this change, both perform identically for me.

---

<div class="post-metadata">

**Author:** ![gideonsimpson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gideonsimpson/32/1928_2.png) [@gideonsimpson](https://discourse.julialang.org/u/gideonsimpson)\
**Post date:** [May 9, 2021, 5:42pm UTC](https://discourse.julialang.org/t/multithreading-and-generated-code/60828/6 "2021-05-09T17:42:50Z")

</div>

Yea, that works well. And for my original problem (without the threading) this runs pretty well:

```julia
function update(x, Δt)
    return x + 0.5 * Δt * x
end

function compute_avg_values_mean_ntuple(x0_vals, Δt, nΔt, f::Tuple{Vararg{<:Any,K}}) where {K}
    x_vals = copy(x0_vals);

    nx = length(x0_vals);
    f_avg_vals = zeros(K, nΔt)

    for n in 1:nΔt
        @. x_vals = update(x_vals, Δt);
        for j in nx
            ntuple(k->f_avg_vals[k,n] += (f[k])(x_vals[j])/nx, K)
        end
    end
    return f_avg_vals
end

Random.seed!(100);
x0 = randn(100);
x₀ = [1.0];
Δt = 0.5;
nΔt = 10^3;

f1(x) = sin(x)
f2(x) = x^2
@btime compute_avg_values_mean_ntuple($x0, $Δt, $nΔt, $(f1,))
@btime compute_avg_values_mean_ntuple($x0, $Δt, $nΔt, $(f2,))
@btime compute_avg_values_mean_ntuple($x0, $Δt, $nΔt, $(f1,f2))

  66.077 μs (2 allocations: 8.81 KiB)
  26.393 μs (2 allocations: 8.81 KiB)
  67.609 μs (2 allocations: 16.62 KiB)

```

I’ve got some separate threading issue though that seems to be independent of the `tuple` of unctions issue, and I’m gonna start a new thread.