Parallelism & Concurrency

Parallelism & Concurrency#

Your laptop has somewhere between 4 and 16 cores. Your server has more. Your GPU has thousands. Most code uses one of them.

Julia gives you four distinct ways to change that, and the first thing to get right is knowing which one your problem needs — because they solve genuinely different problems and picking wrong wastes a lot of effort.

You have

Use

Chapter section

Waiting on I/O, network, disk

Tasks (async)

below

CPU-bound work, shared memory, one machine

Threads

below

Work that doesn’t fit on one machine

Distributed

below

Massively parallel numerical work

GPU

below

Tasks: concurrency without parallelism#

A Task is a lightweight coroutine. Tasks let a program do something else while waiting — they don’t make computation faster, they stop it from idling.

t = @async begin
    sleep(2)
    "done waiting"
end

println("this prints immediately")
result = fetch(t)      # blocks until the task finishes

The classic use is I/O. Fetching ten URLs one at a time takes ten round trips; fetching them concurrently takes about one.

using Downloads

urls = ["https://example.com/a", "https://example.com/b", "https://example.com/c"]

tasks = [@async Downloads.download(u) for u in urls]
files = fetch.(tasks)

Julia’s I/O is asynchronous throughout, so this works without any special async-aware library. There’s no async/await colouring problem: an ordinary function called inside a task just works.

Channel lets tasks communicate, which is how you’d build a producer/consumer pipeline:

ch = Channel{Int}(32)          # buffered channel

producer = @async begin
    for i in 1:100
        put!(ch, i)
    end
    close(ch)
end

consumer = @async for x in ch
    process(x)
end

Threads: real parallelism on one machine#

Threads run genuinely simultaneously on separate cores. You must start Julia with threads enabled:

julia -t auto        # use all available cores
julia -t 8           # use exactly 8

julia> Threads.nthreads()
8

Note

Threading defaults changed in 1.12 As of Julia 1.12, -t auto (or any multi-thread request) also gives you an extra interactive thread in its own pool, so a long-running @threads loop can’t lock up your REPL. Asking explicitly for -t1 still gives you exactly one thread. Julia also now respects CPU affinity settings, so it won’t grab cores your container isn’t entitled to.

On the 1.10 LTS you request the interactive pool manually: -t8,1.

@threads#

The simplest way to parallelize a loop:

using Base.Threads

function parallel_norms(rows)
    out = Vector{Float64}(undef, size(rows, 1))
    @threads for i in axes(rows, 1)
        r = @view rows[i, :]
        out[i] = sqrt(sum(abs2, r))
    end
    return out
end

The loop body must be independent across iterations — that’s the condition. Writing out[i] is fine because each iteration touches a different element.

@spawn: dynamic scheduling#

@threads splits the range into equal chunks, which is wrong when iterations take wildly different amounts of time. @spawn creates tasks the scheduler distributes as workers free up:

tasks = [Threads.@spawn expensive(x) for x in data]
results = fetch.(tasks)

Use @threads for uniform work, @spawn for uneven work.

Race conditions#

This is the part that bites.

# BROKEN
total = 0
@threads for i in 1:1_000_000
    total += i          # multiple threads writing the same variable
end

That silently gives the wrong answer, and — worse — it gives a different wrong answer each run. Two threads read total, both add, both write, one update is lost.

Three correct approaches:

# 1. Per-thread accumulators, combined at the end (usually fastest)
function threaded_sum(xs)
    partials = zeros(eltype(xs), nthreads())
    @threads for i in eachindex(xs)
        partials[threadid()] += xs[i]
    end
    return sum(partials)
end

# 2. Atomics — correct, but contention hurts at high thread counts
total = Threads.Atomic{Int}(0)
@threads for i in 1:1_000_000
    Threads.atomic_add!(total, i)
end

# 3. A lock — for anything more complex than a single number
lk = ReentrantLock()
@threads for x in data
    result = compute(x)
    lock(lk) do
        push!(shared_results, result)
    end
end

Warning

The rules of shared mutable state

  1. Two threads writing the same memory without synchronization is a bug, even if it seems to work.

  2. It will pass your tests and fail in production, because timing.

  3. Prefer designs where each thread writes only to its own slot.

Threads.@threads doesn’t check any of this for you. Julia gives you real threads and real responsibility.

Distributed: more than one machine#

When one machine isn’t enough, Distributed runs separate Julia processes — potentially on other computers — with no shared memory.

using Distributed
addprocs(4)                   # 4 local worker processes
# addprocs([("server.example.com", 8)])   # or remote, over SSH

@everywhere using Statistics  # load on all workers

results = pmap(expensive_analysis, dataset_chunks)

pmap distributes work across processes and collects the results. @distributed handles reductions:

total = @distributed (+) for i in 1:1_000_000
    sqrt(i)
end

Because processes don’t share memory, every argument and result must be serialized and sent. That’s a real cost. Distributed computing pays off when each work unit is substantial relative to the data it needs — big computation, small payload. If you’re shipping a gigabyte to do a millisecond of work, you’ve made things worse.

GPU#

For array-shaped numerical work, this is where the biggest wins are.

using CUDA        # NVIDIA. Also: AMDGPU.jl, Metal.jl (Apple), oneAPI.jl (Intel)

x_cpu = rand(Float32, 10_000_000)
x_gpu = CuArray(x_cpu)          # move to GPU memory

y = x_gpu .^ 2 .+ 1.0f0         # runs on the GPU
result = Array(y)               # bring it back

Look at that middle line: it’s the same broadcasting syntax from chapter seven. CuArray is an AbstractArray, so every generic function written against the array interface works on it. Nothing was rewritten for the GPU.

This is multiple dispatch paying off at the hardware level. A package author writes a numerical algorithm against AbstractArray; a different author writes a GPU array type; the algorithm runs on the GPU without either of them coordinating. KernelAbstractions.jl takes this further, letting you write one kernel that compiles for NVIDIA, AMD, Intel, and Apple hardware.

For anything beyond elementwise work you can write a custom kernel:

function gpu_add!(c, a, b)
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    if i <= length(c)
        @inbounds c[i] = a[i] + b[i]
    end
    return nothing
end

@cuda threads=256 blocks=cld(n, 256) gpu_add!(c, a, b)

That’s Julia code, compiled to GPU machine code, with no CUDA C anywhere. For me — someone who spent a pandemic unable to get research code onto accelerators — this is close to the whole point of the language.

GPU rules of thumb: use Float32 (consumer GPUs are dramatically slower at Float64); minimize transfers between host and device; and only bother when the arrays are large, because transfer overhead swamps small problems.

Choosing, in practice#

Ask two questions.

Am I waiting, or computing? Waiting → tasks. Computing → threads.

Does it fit on one machine? Yes → threads. No → distributed.

And before any of it: is the single-threaded version already type stable and non-allocating? Parallelizing badly-written code gives you a badly-written program running on eight cores. Fix the serial version first; sometimes that’s enough and you’re done. Amdahl’s law is unsentimental — if 20% of your runtime is serial, infinite cores only buy you 5x.

Try it yourself#

  1. Start Julia with -t auto and check Threads.nthreads().

  2. Write a single-threaded sum of sqrt(i) over 10 million values, then a @threads version with per-thread accumulators. Compare with @btime.

  3. Write the deliberately broken racy version, run it five times, and observe that the answer changes. This is worth doing once so the failure mode is real to you.

  4. Use @async to sleep for 1 second three times concurrently, and confirm the total is about 1 second, not 3.

Solutions#

# 2.
using Base.Threads, BenchmarkTools

serial_sum(n) = sum(sqrt(i) for i in 1:n)

function threaded_sum(n)
    partials = zeros(Float64, nthreads())
    @threads for i in 1:n
        partials[threadid()] += sqrt(i)
    end
    return sum(partials)
end

@btime serial_sum(10_000_000)
@btime threaded_sum(10_000_000)

# 4.
@time begin
    ts = [@async sleep(1) for _ in 1:3]
    wait.(ts)
end          # ~1 second

You can now use every core in the machine. Next, a short chapter on what’s actually happening when Julia compiles your code — because now you have the context to appreciate it.