Multi-threaded Array building

General Usage Performance
1

Dear all,
I have a code of this kind to build an array:

N= 2000
A = zeros(Float64,N,N)
B # some long (>>10_000) vector of custom struct
for idx in arrayofindex #> lots (>>10_000) of elements
      i,j,v = processing(B[idx])
     A[i,j] += v
end

It is quite long as processing() take some times and I want to multithreaded it ( processing() is thread-safe).
I have read https://discourse.julialang.org/t/editing-an-array-with-multiple-threads/25364 and https://discourse.julialang.org/t/thread-safe-array-building/3275/2 but these answers (that are using of A[threadid()]) seems outdated.

What is now (with julia 1.10 or 1.11) the best way to multithread that kind of array building function?

2

I just saw that a very similar question Pool of locks? appears when I was writing this. I’ll see if the solution given there answer my problem

3

You want to do data conversion? from B to A?
if this is a deterministic procedure, why don’t you just do it once and then save A?
If this data conversion procedure the main workload in your overall procedure?

I think Threads.@threads can do what you want.

4

No it is not a conversion and (i,j) can appear multiple time.
I could have use @threads to build an array of (i,j,v) for all idx and then performing the reduction but in my actual problem I can’t really store all (i,j,v) (it is a lot of memory).

5

What about

A = [Threads.Atomic{Float64}(0) for _ = 1:N, _ = 1:N]
Threads.@threads for k = eachindex(B)
    local (i, j, v) = processing(B[k])
    Threads.atomic_add!(A[i, j], v)
end

Finally recover A = map(x -> x[], A).

1 Like
6

Is your array sparse?

7

No A is basically an image with many patterns of parameters stored in B

8

Thanks, I’ll try this and the ones given in Pool of locks? .
I will report a small benchmark as soon I have the time to perform it

9

You could

  1. Use task_local_storage to have a per-task array (rather than per-thread). e.g. No more threadid indexing? [thread-local storage] - #12 by stevengj
  2. Use ChunkSplitters.jl (ala Why has simple threading using `threadid` become so complex in v1.12...? - #28 by foobar_lv2) to create a per-chunk array. (This has the downside of imposing essentially a static parallelization schedule, which might be bad if processing(B[idx]) could take a very different amount of time depending on the argument.)
  3. Use some higher-level primitives like in OhMyThreads.jl or Transducers.jl, e.g. expressing this as a parallel reduction operation.
10

Why not use atomics? E.g.

N = 2000
buffers = [Threads.Atomic{Float64}() for _ in 1:N, _ in 1:N]
B # some long (>>10_000) vector of custom struct
OhMyThreads.@tasks for idx in arrayofindex #> lots (>>10_000) of elements
      i,j,v = processing(B[idx])
     Threads.atomic_add!(buffers[i, j], v)
end
A = getindex.(buffers)