# CellListMap: How do I loop through its neighbour list in a thread-safe way?

**URL:** <https://discourse.julialang.org/t/celllistmap-how-do-i-loop-through-its-neighbour-list-in-a-thread-safe-way/92143>\
**Category:** Modelling & Simulations\
**Created:** [December 26, 2022, 12:10pm UTC](https://discourse.julialang.org/t/celllistmap-how-do-i-loop-through-its-neighbour-list-in-a-thread-safe-way/92143 "2022-12-26T12:10:20Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [December 26, 2022, 12:10pm UTC](https://discourse.julialang.org/t/celllistmap-how-do-i-loop-through-its-neighbour-list-in-a-thread-safe-way/92143/1 "2022-12-26T12:10:20Z")

</div>

Hello!

In the following piece of code, I noticed that adding `Threads.@threads` would lead to the wrong result. My limited understanding is that it is because my arrays are not “thread-safe” - I want to learn how to make them so. The code is given below:

```julia
using CellListMap
using StaticArrays

N = 100
points = rand(SVector{3,Float64},N)
system = InPlaceNeighborList(x=points, cutoff=0.5, parallel=false)
list = neighborlist!(system)
# List is an vector of tuples in the format of; (index of point i, index of point j, euclidean distance between i and j)

result_single_thread = zeros(N)
for L in list
    i = L[1]; j = L[2]; d = L[3]

    result_single_thread[i] += d
    result_single_thread[j] += d
end

result_multi_thread = zeros(N)
Threads.@threads for L in list
    i = L[1]; j = L[2]; d = L[3]

    result_multi_thread[i] += d
    result_multi_thread[j] += d
end

deviance_between_single_and_multi_thread = sum(result_single_thread .- result_multi_thread)

println("Deviance between single and multi thread")
println(deviance_between_single_and_multi_thread)
println("Sum SINGLE: $(sum(result_single_thread))")
println("Sum SINGLE: $(sum(result_multi_thread))")

```

Which returns:

```julia
Deviance between single and multi thread
17.944560782803904
Sum SINGLE: 1050.8378293270782
Sum MULTI: 1032.8932685442744

```

Which means that my result using the multi-threaded approach is wrong.

How would I fix this issue?

Kind regards

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [December 26, 2022, 12:11pm UTC](https://discourse.julialang.org/t/celllistmap-how-do-i-loop-through-its-neighbour-list-in-a-thread-safe-way/92143/2 "2022-12-26T12:11:25Z")

</div>

@lmiq pinging you as you told me it was fine for me to do so 🙂

I don’t think it has anything to do with CellListMap as such, I just included it in the title if someone else in the future might have a similar struggle as mine / gap in knowledge.

Kind regards

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [December 26, 2022, 1:27pm UTC](https://discourse.julialang.org/t/celllistmap-how-do-i-loop-through-its-neighbour-list-in-a-thread-safe-way/92143/3 "2022-12-26T13:27:12Z")

</div>

You need to allocate one array of results for each thread and add them at the end. Packages as `Floops.jl` do that for you.

I use this package in many cases: [GitHub - m3g/ChunkSplitters.jl: Simple chunk splitters for parallel loop executions](https://github.com/m3g/ChunkSplitters.jl)

With this pattern:

```julia
using ChunkSplitters
function chunk_splitters(N, list; nchunks = Threads.nthreads())
    results = [zeros(N) for _ in 1:nchunks]
    Threads.@threads for (list_range, ichunk) in chunks(list, nchunks)
        for i in list_range
            i, j, d = list[i] 
            results[ichunk][i] += d
            results[ichunk][j] += d
        end 
    end
    return sum(results)
end

```

Note, however, that with such a simple operation, you will probably be memory-bounded (see [this thread](https://discourse.julialang.org/t/how-to-perform-parallel-vector-addition/92122/2) for some discussion about this).

I can get a good speedup from multi-threading for a list of pairs of about 3 millions pairs.

Here is a complete running example, in which I preallocate the arrrays used for multi-threading:

> **Code**
>
> ```julia
> using Test
> using CellListMap
> using StaticArrays
> using BenchmarkTools
> using ChunkSplitters
> 
> function single_threaded(result, list)
> result .= 0.0
> for i in eachindex(list)
> i, j, d = list[i]
> result[i] += d
> result[j] += d
> end
> return result
> end
> 
> function chunk_splitters(result, list, nchunks, results_threaded)
> # reset
> result .= 0.0
> for i in eachindex(results_threaded)
> results_threaded[i] .= 0.0
> end
> # compute for each chunk
> Threads.@threads for (Lrange, ichunk) in chunks(list, nchunks)
> for i in Lrange
> i, j, d = list[i] 
> results_threaded[ichunk][i] += d
> results_threaded[ichunk][j] += d
> end 
> end
> # reduce
> for i in eachindex(results_threaded)
> result .+= results_threaded[i]
> end
> return result
> end
> 
> function main()
> 
> N = 10^4
> # generates points with the atomic density of water and a realistic cutoff
> # for molecular systems
> points, box = CellListMap.xatomic(N)
> 
> system = InPlaceNeighborList(x=points, cutoff=box.cutoff)
> list = neighborlist!(system)
> result = zeros(N)
> 
> nchunks = Threads.nthreads()
> results_threaded = [copy(result) for _ in 1:nchunks]
> 
> println("Number of pairs in list: $(length(list))")
> 
> @test single_threaded(result, list) ≈ chunk_splitters(result, list, nchunks, results_threaded)
> @btime single_threaded($result, $list)
> @btime chunk_splitters($result, $list, $nchunks, $results_threaded)
> 
> nothing
> end
> 
> ```

With which I get:

```julia
julia> main()
Number of pairs in list: 2689668
  7.570 ms (0 allocations: 0 bytes)
  2.188 ms (49 allocations: 4.36 KiB)

```

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [December 26, 2022, 3:48pm UTC](https://discourse.julialang.org/t/celllistmap-how-do-i-loop-through-its-neighbour-list-in-a-thread-safe-way/92143/4 "2022-12-26T15:48:32Z")

</div>

Thank you very much!

I gave it a go, but it broke my code, so I will have to try re-implementing it at some point step-by-step. I could confirm though that the error in result is gone, so the splitting of array results is working as intended.

As you mention it will probably not speed up a simulation with a low N, but I like learning how to do both single and multi-threaded versions.

Kind regards

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [December 29, 2022, 3:44pm UTC](https://discourse.julialang.org/t/celllistmap-how-do-i-loop-through-its-neighbour-list-in-a-thread-safe-way/92143/5 "2022-12-29T15:44:31Z")

</div>

Just a short update.

I got it to work now, no issues at all, so that is great!

As you mentioned at this low resolution I am memory bound, when I increase number of particles then I get “execution” bound - so at the end it was not worth it to do a threaded version for the CPU code 🙂

Thank you for me showing me how, I learned quite a bit 🙂

Kind regards

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [December 29, 2022, 3:48pm UTC](https://discourse.julialang.org/t/celllistmap-how-do-i-loop-through-its-neighbour-list-in-a-thread-safe-way/92143/6 "2022-12-29T15:48:02Z")

</div>

> [@Ahmed\_Salih](#):
>
> I increase number of particles then I get “execution” bound

What you mean?

By the way, since you are updating the lists at every step, the best is to not compute the neighbor lists at all, and just map the force computation directly using CellListMap, and that will run already in parallel, by default.

For large systems the memory used to store the neighbor list can be limiting, and computing the list is almost as expensive as computing forces for the same set of particles.

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [December 29, 2022, 3:52pm UTC](https://discourse.julialang.org/t/celllistmap-how-do-i-loop-through-its-neighbour-list-in-a-thread-safe-way/92143/7 "2022-12-29T15:52:17Z")

</div>

I don’t know the right word, but what I meant is something like this:

When I have few particles, the threaded approach is slower  
When I have a lot of particles, the threaded approach is faster, but still too slow for my purposes

So when I am using the code “to learn” it is better for me to use fewer particles and a single-threaded approach

And when I have finally learned enough, I would try to use CellListMap with GPU to speed up the process, since the force calculation becomes very slow on CPU 🙂

Kind regards