Parallel operations over arrays in Julia

Question

What is the parallel (over multiple CPUs) version of this code in Julia?

V = zeros(3)
for i = 1:100000
    cc = rand(1:3)
    V[cc] += 1
end

I'm not certain that loop is actually thread safe - since two processes could grab the same index and try to update it at the same time. You may find this section of the manual and the section it's part of helpful. — kevbonham
@kevbonham It is potentially not safe but other languages have provisions for this, like atomic operations, or making copies of variables for each process and then reducing them safely. — Lado Samushia

Bogumił Kamiński Bogumił Kamiński · Accepted Answer · 2019-10-25T19:24:43

This is a direct rewrite of your loop that is thread, tread-safe and avoiding false sharing:

using Random
using Base.Threads

V = let
    mt = Tuple([MersenneTwister() for _ in 1:nthreads()])
    Vv = Tuple([zeros(3) for _ in 1:nthreads()])
    @threads for i = 1:100000
        @inbounds cc = rand(mt[threadid()], 1:3)
        @inbounds Vv[threadid()][cc] += 1
    end
    reduce(+, Vv)
end

However, in general for such a small job probably using threading will not give you much benefit. Also if you really need performance probably the code should be restructured a bit e.g. like this:

function worker(iters, rng)
    v = zeros(3)
    for i = 1:iters
        cc = rand(rng, 1:3)
        v[cc] += 1
    end
    v
end

V = let
    mt = Tuple([MersenneTwister() for _ in 1:nthreads()])
    Vv = [zeros(3) for _ in 1:nthreads()]
    jobs_per_thread = fill(div(100000, nthreads()),nthreads())
    for i in 1:100000-sum(jobs_per_thread)
        jobs_per_thread[i] += 1
    end
    @assert sum(jobs_per_thread) == 100000
    @threads for i = 1:nthreads()
        Vv[threadid()] = worker(jobs_per_thread[threadid()], mt[threadid()])
    end
    reduce(+, Vv)
end

Also under Julia 1.3 you will not have to do manual MersenneTwister management, as Julia will create separate PRNG per thread.

Parallel operations over arrays in Julia

1 Answers