5 — GPU

Thousands of threads — provided you feed them

Emmanuel Pilliat — ENSAI 3A

Goals

• A GPU is not a faster CPU: the opposite bet.
• CuArray ≠ Array: the bus can cost more than the maths.
• Feed the GPU: when it pays off, when it loses.
• Coalescing: locality on the card.

▶ Script: 4_memory_and_gpu.jl (sections 4–7) · Notebook: Memory & GPU · 🐍 PyTorch CPU vs CUDA (Colab)

⬅ 4 — Memory & batching · Synthesis ➡

The batching deck ended on “feed the hardware”. This is the hardware that starves fastest.

Part A — moving to the card

Two strategies for the same problem

Same wall: memory is far away. Opposite bets.

	CPU	GPU
threads	few, fast	thousands, slow
caches	big (24 MiB L3)	small
hides the wait by	caching	switching threads

CPU avoids the wait, GPU covers it. Too little work: the wait shows.

GPU: when one group of threads waits on memory, another group runs. With too little work there is nobody to switch to.

Same code, different hardware

Ag = cu(A)          # copy to VRAM
Bg = cu(B)

A  * B              # CPU
Ag * Bg             # GPU: the SAME `*`

• Multiple dispatch picks the GPU method.
• The algorithm didn’t change.

VRAM = the card’s own memory; the bus links RAM and VRAM. * is a generic function; its arguments are now CuArrays. That is the dispatch module’s promise, cashed in.

The array changed, even if the code didn’t

A  = rand(3, 3)     # Matrix{Float64}, RAM
Ag = cu(A)          # CuArray{Float32,…}, VRAM

• cu copies across the bus: two objects, two memories.
• Float64 → Float32 (CuArray(A) keeps Float64).
• Ag[1] is an error: one element = one bus round trip.

▶ Script: cell “4. The GPU — a CuArray is NOT a CPU array”

The type change is a copy you paid for.

Dispatch picked the GPU method because the type changed. Measured live: typeof(Ag) prints CuArray{Float32, 2, CUDACore.DeviceMemory}. Float32 is the card’s native precision; the Float64 -> Float32 narrowing is the most common “why are my results different on GPU” question. Scalar indexing throws by design (CUDA.@allowscalar overrides it — mention only if asked).

What the copy costs

Two 4096² Float32 matrices (64 MiB each):

step	time
host → device, cu(A)	8.5 ms
compute, Ag * Bg	35 ms
device → host, Array(Ag)	52 ms

• Round trip = 1.7× the computation. Back is the expensive way.
• Don’t ping-pong: upload once, bring back a scalar.

▶ Script: cell “What the copy costs”

A pipeline that copies back after every step lives on the bus, not on the GPU. The torch twin (python/4) measures the same shape: ~10 ms up, ~50 ms down, ~28 ms of compute.

Part B — feed the GPU

The big matmul

Same *, 4096², Float32:

	time
CPU, A * B	hundreds of ms (varies)
GPU, CUDA.@sync Ag * Bg	~35 ms: ~15–20×

• CUDA.@sync: GPU calls return immediately — else you time the launch.
• 1 FLOP = one float operation; this matmul = 2·4096³ ≈ 137 G.

▶ Script: cell “The matrix product, CPU vs GPU”

Big, regular job: what the GPU wants. Smaller ones?

The CPU side swings a lot on this laptop (366–942 ms seen in one session, power profile and hybrid cores); the GPU side is stable at ~35 ms. Quote the order of magnitude, not one ratio.

Bet before you look 🎲

N × N matmul, CPU vs GPU, N from 64 to 4096.

Does the GPU always win?

▶ Script: cell “5. Feed the GPU — when does it start paying off?”

When does the GPU start paying off?

Float32, CPU measured first:

N	CPU	GPU	speedup
64	0.01 ms	0.03 ms	0.25× loses
128	0.03 ms	0.03 ms	1.2×
512	1.15 ms	0.09 ms	13.0×
2048	69.8 ms	3.02 ms	23.1×
4096	hundreds of ms	~35 ms	~15–20×

• Q1. Why slower for small N? Which costs don’t depend on N?
• Q2. Speedup: grows or plateaus? What’s the ceiling?

Not faster: wider.

Narrow work doesn’t get faster on a wider machine. Q1 of the script: fixed costs (launch, sync) don’t depend on N; at N = 64 the GPU spends its time being launched. Q2: the speedup peaks (1024–2048), then levels off — both sides hit their throughput ceiling (the card goes from ~18 to ~5000+ GFLOP/s once the work fills it). The dip at 4096 is not explained here; plausibly a laptop power/thermal limit. Shape is robust, exact ratios are not. On a laptop, measure the CPU BEFORE the GPU: they share a power budget, and the CPU slows down after the GPU has worked (measured once: 75 ms vs 2.97 ms for the same gemm, throttled).

Feeding it: the batch, measured

One layer W (4096²) × B examples, as in the batching deck:

B	GPU total	GFLOP/s
1	0.41 ms	82.8
8	0.42 ms	638
64	0.45 ms	4761
256	2.23 ms	3853

• B = 1 → 64: 64× the work, < 2× the time.
• At B = 1 the card was idle.

▶ Script: cell “Batching on the GPU”

Same GPU; we just asked more of it. Do NOT say “the GPU loses at B=1”: measured, it still beats the CPU at B=1. The honest line is “you are using ~2% of what you paid for” (82.8 vs 4761 GFLOP/s, relative to the best of this run, not the spec sheet). “It loses” invites “so use a CPU”; “~2%” invites “so batch”. B=64 total is run-dependent (0.45–0.63 ms seen). Torch twin: 0.40 → 0.46 ms, B=1 → 64.

Part C — coalescing: locality on the card

What a warp wants

Warp = 32 threads in lockstep, loading together.

• Neighboring addresses → few transactions. Coalesced.
• Scattered addresses → one per thread.

The column-major lesson again, on the card.

Coalescing, measured

idx = CuArray(Int32.(((0:N-1) .* s) .% N .+ 1))   # stride s
a[idx]                                            # 20 M reads

stride s	time
1	2.24 ms	coalesced
33	5.40 ms	2.4× slower
1023	12.3 ms	5.5× slower

▶ Script: cell “6. Coalescing — locality again, at the scale of a warp”

“Only 2.4×?” Which other bytes does a[idx] touch?

Same 20 M reads, different order: i ↦ (i·s mod N) + 1. The stride must be COPRIME with N or it is not a permutation: gcd(32, 20_000_000)=32 would reread only 625k slots — a cache demo, not a coalescing demo. 33 works; 65 and 1025 do NOT (they share the factor 5). Script question: why must the index be a permutation?

Testing that explanation 🔬

a[idx] also reads idx and writes out, in order: only 1/3 of the bytes are punishable.

Index computed in the kernel → 1/2 punishable. The ratio must rise.

i = (blockIdx().x - 1) * blockDim().x + threadIdx().x   # my thread
j = ((i - 1) * s) % n + 1     # computed, never loaded
out[i] = a[j]

stride 33 vs 1	ratio
a[idx] (1/3)	2.4×
kernel (1/2)	4.18×

▶ Script: cell “Test that explanation”

Kernel = a function run once per GPU thread; thread i handles element i. The prediction held: an explanation you can test, and might have lost.

Monte-Carlo π

Monte-Carlo π on the card

function pi_gpu(n)
    x = CUDA.rand(Float32, n)            # in VRAM
    y = CUDA.rand(Float32, n)
    inside = (x.^2 .+ y.^2) .<= 1f0
    return 4 * sum(inside) / n           # one scalar back
end

• Same count as count_in_disk; 100 million GPU threads.
• Only one number crosses the bus.
• No shared counter: sum reduces.

▶ Script: cell “7. Monte-Carlo π on the GPU”

1f0 = 1 as a Float32. Minimize sharing, from 22 threads to 100 million.

π: which CPU do we compare to?

100 million points, pi_gpu ≈ 12 ms:

CPU baseline	GPU gain
pi_cpu_vec (allocates 2 × 400 MB)	~35×
count_in_disk, scalar loop	~25×

• Fair baseline: the scalar loop. ~25×.

Q3. Threads module: 22 threads → 6.3×. Why is the card so far ahead?

Q4. What in pi_gpu does not benefit from 100 M threads?

pi_cpu_vec mirrors pi_gpu but pays the allocation traffic; allocation-bound, very load-sensitive. Threads deck: @spawn ×22 gave 6.3×. Q3: the card runs thousands of threads and never leaves VRAM. Q4 answer: the final sum reduction (and the launches) — the Amdahl part.

What we saw on the GPU

• CuArray in VRAM: same code, you pay the copy.
• Bus is expensive: bring back a scalar.
• Small jobs lose: batch.
• Coalescing: neighbor threads, neighbor addresses.

The four ideas of the course: Synthesis.

Going further

• ▶ Script: scripts/4_memory_and_gpu.jl, sections 4–7
• ▶ Notebook: Memory & GPU
• 🐍 Python: Memory & GPU (Colab)
• ✏️ Activity: Julia · 🐍 Python (Colab)
• 🎯 Not faster, wider. Feed it; keep data on the card.

⬅ Previous: 4 — Memory & batching · Next: Course overview & synthesis ➡