# Why does Julia threads in the following way?

**URL:** <https://discourse.julialang.org/t/why-does-julia-threads-in-the-following-way/124953>\
**Category:** General Usage\
**Created:** [January 20, 2025, 8:25am UTC](https://discourse.julialang.org/t/why-does-julia-threads-in-the-following-way/124953 "2025-01-20T08:25:53Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![ttx002000](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ttx002000/32/208489_2.png) [@ttx002000](https://discourse.julialang.org/u/ttx002000)\
**Post date:** [January 20, 2025, 8:25am UTC](https://discourse.julialang.org/t/why-does-julia-threads-in-the-following-way/124953/1 "2025-01-20T08:25:53Z")

</div>

Here is a working example. There are two functions, the Coulomb\_matrix() function just takes in a vector of Float numbers, and gives out a matrix of float numbers. The test() is having strange behaviour. Here is the full script:

```julia
function Coulomb_matrix(kvec::Vector{Float64})::Matrix{Float64}
   
    s1=[9047.5636*exp(-norm(kvec)) for ja in 1:10, jb in 1:10]
    cmatrix=kron(ones(Float64,4,4),s1)
    
    return cmatrix
  
end
  

  
function test()
   radius=1.4
   num_points=20
    dimension=40
    kx_grid=collect(range(-radius/2, stop=+radius/2, length=num_points))
    ky_grid=collect(range(-radius/2, stop=+radius/2, length=num_points))

   
  k_set=Vector{Float64}[]
  
      
  for ja in eachindex(kx_grid),jb in eachindex(ky_grid)  
      push!(k_set,[kx_grid[ja],ky_grid[jb]]) 
  end

   l1=zeros(ComplexF64,dimension,dimension,length(k_set))
   l2=zeros(ComplexF64,dimension,dimension,length(k_set))
  
   density_matrix=ones(ComplexF64,dimension,dimension,length(k_set))

  
  
    Threads.@threads for ja in eachindex(k_set)
        for jb in eachindex(k_set)
         cmatrix=Coulomb_matrix(k_set[ja]-k_set[jb])
        l2[:,:,ja]+=(cmatrix.*density_matrix[:,:,jb])
        end
    end

    println("maxl2",sort(vec(abs.(l2)))[end])
  
    
  
    Threads.@threads for ja in eachindex(k_set)
    for jb in eachindex(k_set)
    fcmatrix=Coulomb_matrix(k_set[ja]-k_set[jb])
    l1[:,:,ja]+=(fcmatrix.*density_matrix[:,:,jb])
    end
   end
   println("maxl1",sort(abs.(vec(l1)))[end])
     
  
  
  println("l1,l2",l1==l2)
  println(sum(abs.(l1-l2)))

 
 cmatrix=Coulomb_matrix([0.0,0.0])
  l3=diagm(cmatrix*diag(dropdims(sum(density_matrix,dims=3),dims=3)))
   
 
  
    
    
  end
  

```

I am seeing strange behaviour in the function test(). Supposedly l1 and l2 should be equal to each other. They are both three-dimensional ComplexF64 arrays. The only difference is that when I compute l1 and l2, they need to use some intermediate arrays from the Coulomb\_matrix() function. In one case, I call the intermediate array cmatrix, the other case, I call it fcmatrix (below)

```julia
 Threads.@threads for ja in eachindex(k_set)
        for jb in eachindex(k_set)
         cmatrix=Coulomb_matrix(k_set[ja]-k_set[jb])
        l2[:,:,ja]+=(cmatrix.*density_matrix[:,:,jb])
        end
    end

  Threads.@threads for ja in eachindex(k_set)
    for jb in eachindex(k_set)
    fcmatrix=Coulomb_matrix(k_set[ja]-k_set[jb])
    l1[:,:,ja]+=(fcmatrix.*density_matrix[:,:,jb])
    end
   end

```

The strange thing is that the test(), as it is written, will report l1 is not equal to l2. I was guessing that this is because the threading syntax I wrote was not appropriate, i.e. I probably should not thread over l2[:,:,index], and perform in-place addition (actually why? I had the impression that this should be safe). Since in the code below, I only access one slice of l2 every time.

```julia
 Threads.@threads for ja in eachindex(k_set)
        for jb in eachindex(k_set)
         cmatrix=Coulomb_matrix(k_set[ja]-k_set[jb])
        l2[:,:,ja]+=(cmatrix.*density_matrix[:,:,jb])
        end
    end

```

But the strangest thing is that if you just delete the last two lines, the test() will just report l1 is indeed equal to l2, which confuses a lot. Why would this happen?

```julia
 cmatrix=Coulomb_matrix([0.0,0.0])
  l3=diagm(cmatrix*diag(dropdims(sum(density_matrix,dims=3),dims=3)))
   

```

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [January 20, 2025, 9:29am UTC](https://discourse.julialang.org/t/why-does-julia-threads-in-the-following-way/124953/2 "2025-01-20T09:29:26Z")

</div>

Here’s the important parts:

> [@ttx002000](#):
>
> ```julia
> Threads.@threads for ja in eachindex(k_set)
> for jb in eachindex(k_set)
> cmatrix=Coulomb_matrix(k_set[ja]-k_set[jb])
> l2[:,:,ja]+=(cmatrix.*density_matrix[:,:,jb])
> end
> end
> [...]
> cmatrix=Coulomb_matrix([0.0,0.0])
> 
> ```

What’s going wrong here is that the `cmatrix` inside the `Threads.@threads` block has the same name as the `cmatrix` outside of it, so julia is lexically capturing the variable `cmatrix` into the closure inside the `@threads` loop.

This makes it so that every iteration of the `@threads` block actually has a race condition in it when accessing `cmatrix`.

This is something I really dislike about `Threads.@threads` and how it makes it very easy to write code that causes subtle concurrency bugs.

One way you could fix this would be to write

```julia
    Threads.@threads for ja in eachindex(k_set)
        for jb in eachindex(k_set)
            local cmatrix=Coulomb_matrix(k_set[ja]-k_set[jb])
            l2[:,:,ja]+=(cmatrix.*density_matrix[:,:,jb])
        end
    end

```

and then `cmatrix` will be local to the closure and there’ll be no race condition.

---

<div class="post-metadata">

**Author:** ![eldee](https://avatars.discourse-cdn.com/v4/letter/e/b5a626/32.png) [@eldee](https://discourse.julialang.org/u/eldee)\
**Post date:** [January 20, 2025, 9:34am UTC](https://discourse.julialang.org/t/why-does-julia-threads-in-the-following-way/124953/3 "2025-01-20T09:34:24Z")

</div>

> [@ttx002000](#):
>
> The only difference is that when I compute l1 and l2, […] In one case, I call the intermediate array cmatrix, the other case, I call it fcmatrix (below)

I think the issue is due to the scoping rules. Because you (later) declare a variable `cmatrix`, the `cmatrix` within the threaded for-loop will be same variable and get reassigned and reused by all threads, leading to race conditions.

See also [Race condition caused by variable scope getting lifted from a multithreaded context · Issue #14948 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/issues/14948).

> **Simpler example**
>
> Consider these basic `sum` functions:
> 
> ```julia
> function sum1(x)
> tasks = map(1:nthreads()) do i0
> @spawn begin
> sum = 0.
> for j = i0:nthreads():length(x)
> sum += x[j]
> end
> return sum
> end
> end
> 
> tsums = fetch.(tasks)
> sum = 0
> for i = 1:nthreads()
> sum += tsums[i]
> end
> return sum
> end
> 
> function sum2(x)
> tasks = map(1:nthreads()) do i0
> @spawn begin
> tsum = 0.
> for j = i0:nthreads():length(x)
> tsum += x[j]
> end
> return tsum
> end
> end
> 
> tsums = fetch.(tasks)
> sum = 0
> for i = 1:nthreads()
> sum += tsums[i]
> end
> return sum
> end
> 
> ```
> 
> where the only difference is that in the second case we called the per-thread sum `tsum` instead of `sum`. You’ll note that the behaviour of `sum1` is non-deterministic
> 
> ```julia-repl
> julia> x = rand(1000);
> 
> julia> sum1(x)
> 505.8078846565804
> 
> julia> sum1(x)
> 502.79831723789334
> 
> julia> sum1(x)
> 494.28161467864845
> 
> julia> sum1(x)
> 494.28161467864845
> 
> julia> sum1(x)
> 496.9478723181834
> 
> ```
> 
> while `sum2` works as you’d expect
> 
> ```julia-repl
> julia> sum(x)
> 494.2816146786484
> 
> julia> sum2(x)
> 494.28161467864845
> 
> julia> sum2(x)
> 494.28161467864845
> 
> ```
> 
> In `sum1` the `sum` variable within the tasks is the same as outside of them, meaning that multiple threads will write to the same variable concurrently, leading to data race issues. You can fix this by renaming the variable, or explicitly declaring the inner `sum` variable to be `local`.

---

<div class="post-metadata">

**Author:** ![sgaure](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sgaure/32/14779_2.png) [@sgaure](https://discourse.julialang.org/u/sgaure)\
**Post date:** [January 20, 2025, 11:26am UTC](https://discourse.julialang.org/t/why-does-julia-threads-in-the-following-way/124953/4 "2025-01-20T11:26:43Z")

</div>

While we’re at it, the `l2[:, :, ja] += ...` is lowered to `l2[:, :, ja] = l2[:, :, ja] + ...`, causing an allocation of the sub matrix `l2[:, :, ja]`, as well as `density_matrix[:, : , jb]`.

You can avoid it by using views and broadcasting. Excessive allocations (and resulting garbage collections) can be devastating for performance in parallel programs. If possible, the coulomb matrices should also be precomputed.

```julia
Threads.@threads for ja in eachindex(k_set)
    for jb in eachindex(k_set)
        local cmatrix = Coulomb_matrix(k_set[ja]-k_set[jb])
        @views l2[:, :, ja] .+= cmatrix .* density_matrix[:, :, jb]
    end
end

```

Due to these unfortunate scoping rules, I declare everything `local` inside such constructs.

---

<div class="post-metadata">

**Author:** ![ttx002000](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ttx002000/32/208489_2.png) [@ttx002000](https://discourse.julialang.org/u/ttx002000)\
**Post date:** [January 20, 2025, 4:26pm UTC](https://discourse.julialang.org/t/why-does-julia-threads-in-the-following-way/124953/5 "2025-01-20T16:26:03Z")

</div>

Hi, Thanks to this reply. But I still have one confusion. So the first time I decaled variable “cmatrix” outside the Threads loop is actually after the threads loop. Will this also make julia think that the “cmatrix” inside Threads loop is a global variable?

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [January 20, 2025, 4:27pm UTC](https://discourse.julialang.org/t/why-does-julia-threads-in-the-following-way/124953/6 "2025-01-20T16:27:25Z")

</div>

Yes.

---

<div class="post-metadata">

**Author:** ![Benny](https://avatars.discourse-cdn.com/v4/letter/b/49beb7/32.png) [@Benny](https://discourse.julialang.org/u/Benny)\
**Post date:** [January 20, 2025, 7:29pm UTC](https://discourse.julialang.org/t/why-does-julia-threads-in-the-following-way/124953/7 "2025-01-20T19:29:20Z")

</div>

> [@ttx002000](#):
>
> “cmatrix” inside Threads loop is a global variable

Not a global variable, those originate in the global scope where your function definitions live. But yes, a local variable originates in the outermost local scope it is declared or assigned, and inner local scopes will use it by default even if they were written in preceding lines; you opt out and create a new variable with explicit `local` declarations. Where local variables originate are determined when the entire function scope is parsed from the global scope, so the order of the subsequent execution doesn’t matter.