# 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.

```{code-block} julia
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.

```{code-block} julia
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:

```{code-block} julia
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:

```{code-block} bash
julia -t auto        # use all available cores
julia -t 8           # use exactly 8
```

```{code-block} julia
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:

```{code-block} julia
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:

```{code-block} julia
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.

```{code-block} julia
# 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:

```{code-block} julia
# 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.

```{code-block} julia
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:

```{code-block} julia
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.

```{code-block} julia
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:

```{code-block} julia
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

```{code-block} julia
# 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.