Applied GPU Programming

Two ways to build a processor

Why CPUs and GPUs are designed so differently, and why real programs use both.

• After this module you can: define latency and throughput, explain where a CPU and a GPU spend their transistors, and decide which processor fits a given workload.

1.1 Two measures of speed

Picture moving people across a city. A sports car takes one passenger in 10 minutes. A bus takes 60 passengers in 30 minutes. Which is faster? It depends on what you measure:

• Latency is the time to finish one task. The car wins (10 vs 30 minutes).
• Throughput is the number of tasks finished per unit of time. The bus wins: 2 passengers/minute against 0.1.

Processors face the same trade-off. A chip has a fixed budget of transistors, power and area. You can spend it on making each instruction stream finish quickly, or on running as many streams as possible at once. You cannot maximize both.

1.2 The CPU: built for latency

A CPU runs operating systems, browsers and compilers: code with lots of branches, pointer chasing and dependencies, where step n+1 needs the result of step n. For that code, the only way to go fast is to make one thread fast. So a CPU core spends most of its area on things that are not arithmetic:

• Large caches (often tens of MB) so that most memory accesses never go to slow DRAM.
• Branch prediction to guess which way an if will go and start executing before it is known.
• Out-of-order execution to find independent instructions further ahead and run them while waiting on others.
• High clock speed (4–5 GHz) and a handful of powerful cores (8–64).

1.3 The GPU: built for throughput

A GPU was first built for graphics: compute the colour of millions of pixels, each independently, 60 times per second. Nobody cares how long one pixel takes, only that all of them are done in 16 ms. So a GPU spends its area the opposite way:

• Thousands of simple arithmetic units (2,560 on a T4, 16,896 on an H100).
• Tiny control logic: one instruction decoder drives 32 threads at once (Module 2), and there is no branch prediction or out-of-order execution.
• Small caches per thread. Instead of avoiding memory stalls, the GPU hides them by always having other work ready (Module 2).
• Very high memory bandwidth (320 GB/s on a T4, about 3.35 TB/s on an H100, versus about 100 GB/s for a desktop CPU), because feeding thousands of units needs a wide pipe.

1.4 Heterogeneous computing

The textbook's central idea is that real applications contain both kinds of code: serial control logic (reading input, making decisions) and large data-parallel loops (the same operation on millions of elements). The best design uses each processor for what it is good at. In CUDA this is explicit:

• The host (CPU) runs main(), does the serial parts and decides what the GPU does.
• The device (GPU) runs kernels: functions executed by thousands of threads in parallel.
• They have separate memories connected by a PCIe link (16–64 GB/s), so data has to be copied between them. This cost appears in every assignment.

Exercises

Warps, SIMT and hiding latency

How thousands of slow threads add up to a fast processor.

• After this module you can: explain what a warp is, the difference between SIMD and SIMT, how latency hiding works (with numbers), and what warp divergence costs.

2.1 The problem: memory is slow

Reading from GPU global memory takes roughly 400–800 clock cycles. An arithmetic instruction takes about 4. If a thread computed, then waited 500 cycles for its next number, its arithmetic unit would be idle more than 99% of the time. A CPU fights this with big caches. A GPU uses a different trick.

2.2 Warps: the unit of execution

The GPU groups threads into bundles of 32 called warps. A warp is the real unit of scheduling: the hardware fetches and decodes one instruction and all 32 threads execute it, each on its own data. Threads 0–31 of a block form warp 0, threads 32–63 warp 1, and so on.

This is why GPU control logic is so cheap: one decoder serves 32 threads. It is also why block sizes should be multiples of 32. A block of 48 threads makes 2 warps, the second only half full, wasting 16 lanes on every instruction.

2.3 SIMD versus SIMT

• SIMD (single instruction, multiple data), used in CPU vector units like AVX: one thread issues an instruction that operates on a vector register of 8 or 16 values. The programmer or compiler must pack data into vectors explicitly.
• SIMT (single instruction, multiple threads), NVIDIA's model: you write code for one scalar thread, and every thread gets its own ID and registers. The hardware groups 32 threads into a warp and runs them together. You get SIMD efficiency while writing simple scalar code.

2.4 Latency hiding by switching warps

Each SM holds many warps at once (up to 32 on a T4, 64 on an H100). Every cycle, the warp scheduler picks a warp whose next instruction is ready and issues it. When warp 0 issues a memory load, it simply waits; the scheduler runs warp 1, then warp 2, and so on. By the time it cycles back, warp 0's data has arrived.

Why is switching free? On a CPU, switching threads means saving and restoring registers to memory. On a GPU, every resident warp keeps its own registers in the SM's huge register file (65,536 32-bit registers per SM). Nothing is saved or restored; the scheduler just picks a different warp next cycle.

2.5 How many warps are enough? (Little's law)

Queueing theory gives a rule of thumb: the amount of work that must be in flight equals latency × throughput. If a load takes 500 cycles and the SM can issue one instruction per cycle per scheduler, you need about 500 independent instructions in flight per scheduler to never stall. Each warp contributes only a few independent instructions before it waits, so you need many warps. That is the intuition behind occupancy (Module 8).

2.6 Warp divergence

All 32 threads of a warp share one instruction stream. What happens with if (x > 0) A(); else B(); when some threads have x > 0 and others don't? The warp runs both paths one after the other: first A with the "else" threads masked off (idle), then B with the "if" threads masked. The results are correct; the cost is time.

• Divergence only matters within a warp. Different warps taking different branches cost nothing extra.
• Bounds checks like if (i < n) diverge in at most one warp (the last), so they're essentially free.

Exercises

#include <cstdio>
__global__ void who() {
  int t = threadIdx.x;
  printf("thread %2d  warp %d  lane %2d\n", t, t / 32, t % 32);
}
int main() { who<<<1, 96>>>(); cudaDeviceSynchronize(); }

Inside the GPU: SMs and the memory hierarchy

The physical layout you are programming for, and why it shapes every optimization in this course.

• After this module you can: describe the GPU → SM → core hierarchy, list the memory levels with their scope and speed, compute core counts and bandwidth from specs, and explain what MIG means for you on the cluster.

3.1 The building block: the Streaming Multiprocessor

A GPU is mostly copies of one unit, the Streaming Multiprocessor (SM). Each SM is like a small processor of its own, containing:

• CUDA cores (FP32/INT32 units), plus FP64 units, special-function units (sin, exp, sqrt) and, on newer GPUs, Tensor Cores for matrix math.
• Warp schedulers, usually 4 per SM, each issuing for its own set of warps.
• A register file of 65,536 32-bit registers, shared by all resident threads.
• Shared memory / L1 cache: a fast on-chip memory (64 KB on a T4, up to 228 KB on an H100) that blocks on this SM can use.

When you launch a kernel, the hardware hands out thread blocks to SMs. A block runs entirely on one SM and never moves. An SM can host several blocks at once if it has enough registers, shared memory and warp slots. More SMs therefore means more blocks running at once. This is why code with many blocks automatically gets faster on bigger GPUs.

The compute capability (7.5, 9.0) identifies the instruction set. It is what you pass to the compiler: nvcc -arch=sm_75 for a T4, -arch=sm_90 for an H100. Compiling for the wrong one gives "no kernel image is available" errors, or silently runs nothing if you don't check errors.

3.2 The memory hierarchy

The closer memory is to the cores, the faster and smaller it is, and the fewer threads can see it. Every performance technique in this course is about keeping data in the faster levels.

3.3 MIG on the school cluster

The cluster's H100s are partitioned with Multi-Instance GPU (MIG): one physical card is split into isolated slices, each with its own SMs, memory and L2 share. You get a slice (for example 16 SMs and around 10 GB), not 132 SMs. Run nvidia-smi to see what you got, and compare cluster results fairly: a MIG slice can be slower than a full T4 for some kernels.

Exercises

#include <cstdio>
#include <cuda_runtime.h>
int main() {
  cudaDeviceProp p; cudaGetDeviceProperties(&p, 0);
  printf("%s  cc %d.%d\n", p.name, p.major, p.minor);
  printf("SMs %d  warp %d  maxThr/block %d  maxThr/SM %d\n",
         p.multiProcessorCount, p.warpSize, p.maxThreadsPerBlock, p.maxThreadsPerMultiProcessor);
  printf("regs/SM %d  smem/SM %zu KB  L2 %d MB  mem %.1f GB\n",
         p.regsPerMultiprocessor, p.sharedMemPerMultiprocessor/1024,
         p.l2CacheSize/(1024*1024), p.totalGlobalMem/1e9);
}

How fast can it possibly go?

Two back-of-the-envelope tools that predict performance before you write any code.

• After this module you can: apply Amdahl's law, compute arithmetic intensity for a kernel, decide whether it is memory- or compute-bound, and estimate its best possible runtime.

4.1 Amdahl's law: the serial part limits you

Split a program's original runtime into a fraction s that stays serial and a fraction 1 − s that you speed up by a factor p. The new runtime is s + (1 − s)/p of the original, so:

The serial part never gets faster, so eventually it is all that is left. With 10% serial code, an infinitely fast GPU gives at most 10×. With 1%, at most 100×.

On a GPU, host↔device copies behave like serial time, since the GPU computes nothing useful during them unless you overlap them (streams, a later lecture).

4.2 Two ceilings: compute and memory

A kernel is limited by whichever runs out first:

• Compute ceiling: the GPU can't do more than its peak FLOP/s (8.1 TFLOP/s FP32 on a T4).
• Memory ceiling: it can't move more than its peak bandwidth (320 GB/s on a T4).

Which one you hit depends on the kernel's arithmetic intensity:

The time can't be less than either bound, so a good estimate of the best possible time is:

The crossover is the machine balance = peak FLOP/s ÷ peak bandwidth (≈ 25 FLOP/byte on a T4, from Module 3). Kernels with AI below it are memory-bound: only moving fewer bytes helps. Above it they are compute-bound: only doing fewer or faster operations helps. Plotting attainable FLOP/s against AI gives the roofline model, which Nsight Compute can draw for your kernels.

Exercises

The CUDA workflow and moving data

The skeleton every CUDA program follows, and what each call really does.

• After this module you can: use the three function qualifiers correctly, write the allocate → copy → launch → copy → free sequence, add error checking, and explain which calls are synchronous.

5.1 Three questions

Every CUDA program must answer three questions:

1. Where does each function run? CPU or GPU, set by qualifiers.
2. Where does each piece of data live? Host memory or device memory, set by how you allocate it.
3. What does each thread do? Determined by the thread's index (Module 6).

5.2 Function qualifiers

5.3 Two separate memories

The CPU and GPU have separate address spaces. A pointer returned by malloc points to host memory; one from cudaMalloc points to device memory. They're both just float* to the compiler, so mixing them up compiles fine and fails at runtime. That's why everyone uses the naming convention h_ for host and d_ for device.

• Dereferencing a device pointer on the CPU (d_out[0] in main) usually crashes with a segmentation fault.
• Passing a host pointer to a kernel gives an "illegal memory access" error when the kernel runs.

5.4 The five steps

// 1. allocate on the GPU
float *d_a, *d_c;
cudaMalloc(&d_a, bytes);              // note: address of the pointer
cudaMalloc(&d_c, bytes);

// 2. copy inputs host → device       (dst, src, size, direction)
cudaMemcpy(d_a, h_a, bytes, cudaMemcpyHostToDevice);

// 3. launch the kernel (returns immediately!)
myKernel<<<blocks, threads>>>(d_a, d_c, n);

// 4. copy results device → host (waits for the kernel to finish first)
cudaMemcpy(h_c, d_c, bytes, cudaMemcpyDeviceToHost);

// 5. free
cudaFree(d_a); cudaFree(d_c);

cudaMalloc takes &d_a (a pointer to your pointer) because it needs to write the device address into your variable. Its return value is reserved for an error code.

5.5 Synchronous or asynchronous?

The GPU has a work queue called a stream. The CPU puts operations in it, and the GPU does them in order.

• Kernel launches are asynchronous: the CPU queues the kernel and continues immediately, often before the kernel even starts.
• cudaMemcpy is synchronous: it waits until everything before it in the stream is finished, then does the copy, then returns. So step 4 above automatically waits for the kernel.
• cudaDeviceSynchronize() blocks the CPU until all queued GPU work is done. You need it when nothing else waits: before stopping a timer, before checking for kernel errors, and before exiting if the kernel uses printf.

5.6 Error checking

Every CUDA API call returns a cudaError_t. If you ignore it, failures are silent and you debug garbage output for hours. Wrap every call:

#define CHECK(call) do {                                         \
    cudaError_t err = (call);                                    \
    if (err != cudaSuccess) {                                    \
        fprintf(stderr, "CUDA error %s at %s:%d\n",              \
                cudaGetErrorString(err), __FILE__, __LINE__);    \
        exit(1);                                                 \
    } } while (0)

CHECK(cudaMalloc(&d_a, bytes));
myKernel<<<g, b>>>(…);
CHECK(cudaGetLastError());        // catches bad launch configurations
CHECK(cudaDeviceSynchronize());   // catches errors during kernel execution

Kernel launches return nothing, so there are two checks: cudaGetLastError() right after the launch catches configuration errors (too many threads per block), and the synchronize catches errors while the kernel ran (out-of-bounds accesses).

int main() {
  float h_x[1000000];                                     // A
  float *d_x;
  cudaMalloc(d_x, 1000000 * sizeof(float));               // B
  cudaMemcpy(h_x, d_x, 1000000*sizeof(float), cudaMemcpyHostToDevice);  // C
  square<<<1000000/256, 256>>>(d_x);                       // D
  printf("%f\n", d_x[0]);                                 // E
}

Exercises

Launching kernels and finding your thread

How a launch creates threads, and how each thread works out which data is its own.

• After this module you can: describe the grid/block/thread hierarchy, compute global indices in 1D and 2D, size grids for any problem, respect hardware limits, and use vector types.

6.1 The thread hierarchy

A launch kernel<<<G, B>>> creates a grid of G blocks, each with B threads: G × B threads in total, all running the same kernel code. Why two levels instead of one flat list of threads?

• Blocks are independent and can run in any order on any SM. This is what lets the same code scale from 40 SMs to 132.
• Threads within a block can cooperate: they run on the same SM, can share __shared__ memory and can synchronize with __syncthreads(). Threads in different blocks can't.

Inside the kernel, each thread reads four built-in variables (each has .x .y .z):

6.2 The 1D global index

Think of blocks as rows of seats in a cinema. Your seat number in the whole cinema is (row × seats per row) + your seat in the row:

6.3 Sizing the grid: round up, then guard

N is rarely a multiple of the block size. Integer division rounds down, so N / TPB blocks would leave the tail unprocessed. Always round up, and then make the extra threads in the last block do nothing:

__global__ void vecAdd(const float *a, const float *b, float *c, int n) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i >= n) return;              // guard: the last block may have spare threads
  c[i] = a[i] + b[i];
}
int TPB = 256;
int blocks = (n + TPB - 1) / TPB;   // ceiling division
vecAdd<<<blocks, TPB>>>(d_a, d_b, d_c, n);

Why does (n + TPB − 1) / TPB round up? Adding TPB − 1 pushes any remainder over the next multiple. For n = 1000, TPB = 256: (1000 + 255)/256 = 1255/256 = 4, while 1000/256 = 3.

Choosing TPB: use a multiple of 32 (Module 2), usually 128 or 256. Too small (32) and you hit the per-SM block limit before filling the warp slots; 1024 leaves less flexibility for the scheduler. 256 is a good default.

6.4 Hardware limits

• At most 1024 threads per block in total (x·y·z ≤ 1024). A block of (32, 32) is allowed; (64, 32) is not.
• Grid size: x up to 2³¹ − 1 blocks, y and z up to 65,535.
• Exceeding a limit doesn't crash. The launch silently fails with "invalid configuration argument", which you only see by checking cudaGetLastError().

6.5 2D grids for 2D data

For matrices and images, use 2D blocks and grids with the dim3 type. Map x to columns and y to rows:

__global__ void scale2D(float *m, int W, int H, float s) {
  int c = blockIdx.x * blockDim.x + threadIdx.x;   // column  0..W-1
  int r = blockIdx.y * blockDim.y + threadIdx.y;   // row     0..H-1
  if (c >= W || r >= H) return;                    // guard BOTH dimensions
  m[r * W + c] *= s;                                // row-major: row × width + column
}
dim3 block(16, 16);                                 // 256 threads
dim3 grid((W + 15) / 16, (H + 15) / 16);            // round up each dimension
scale2D<<<grid, block>>>(d_m, W, H, 2.0f);

Why x → column? Threads in a warp are numbered x first: threads (0,0), (1,0), (2,0) … are consecutive lanes. With x as the column, the 32 lanes of a warp access 32 adjacent floats, which the hardware combines into a few large memory transactions. This is coalescing. Swapping x and y makes the lanes access addresses W floats apart, which can be many times slower.

6.6 Vector types

CUDA provides small structs such as int2, int3, float2, float3, float4, double2 with fields .x .y .z .w. Build them with make_float3(a, b, c). dim3 is a vector of unsigned ints where unspecified fields default to 1, so dim3 b(256) means (256, 1, 1). float4 is aligned to 16 bytes and can be loaded in one instruction, a common trick in memory-bound kernels.

Exercises

__global__ void show(int W, int H) {
  int c = blockIdx.x * blockDim.x + threadIdx.x;
  int r = blockIdx.y * blockDim.y + threadIdx.y;
  if (c >= W || r >= H) return;
  printf("block(%d,%d) thread(%d,%d) -> r=%d c=%d idx=%d\n",
         blockIdx.x, blockIdx.y, threadIdx.x, threadIdx.y, r, c, r * W + c);
}
int main() {
  int W = 7, H = 5;
  dim3 block(4, 4), grid((W + 3) / 4, (H + 3) / 4);   // grid = (2, 2)
  show<<<grid, block>>>(W, H);
  cudaDeviceSynchronize();
}

Walkthrough: the distance kernel

Read a complete CUDA program line by line, then find what's wrong with it.

• After this module you can: explain every line of the lecture's first program, predict its output, find its latent bugs, and explain why GPU and CPU floating-point results can differ.

7.1 What it computes

For N = 64 points evenly spaced on [0, 1], x_i = i/(N − 1), compute the distance to a reference point ref = 0.5: |x_i − 0.5|. One thread per point.

#define N 64
#define TPB 32

__device__ float scale(int i, int n) {          // GPU helper: index → position in [0,1]
  return ((float)i) / (n - 1);
}
__device__ float distance(float x1, float x2) {  // GPU helper: |x2 - x1| written as sqrt of square
  return sqrt((x2 - x1) * (x2 - x1));
}
__global__ void distanceKernel(float *d_out, float ref, int len) {
  const int i = blockIdx.x * blockDim.x + threadIdx.x;   // which point am I?
  const float x = scale(i, len);
  d_out[i] = distance(x, ref);                            // write my result
}

int main() {
  const float ref = 0.5f;
  float *h_out = (float *)malloc(N * sizeof(float));     // host buffer for results
  float *d_out = 0;
  cudaMalloc(&d_out, N * sizeof(float));                 // device buffer
  distanceKernel<<<N / TPB, TPB>>>(d_out, ref, N);       // 2 blocks × 32 threads
  cudaMemcpy(h_out, d_out, N * sizeof(float), cudaMemcpyDeviceToHost);  // waits for kernel
  for (int i = 0; i < N; i++) std::cout << i << " = " << h_out[i] << std::endl;
  cudaFree(d_out); free(h_out);
}

Note what's not there: no input copy to the GPU. The kernel computes its inputs from its index, so only the output is copied back.

7.2 The output

The values go 0.5, 0.484, … down to 0.0079 at i = 31 and back up to 0.5 at i = 63, the V shape of |x − 0.5|. There's no exact zero because 0.5 lies between x₃₁ = 31/63 and x₃₂ = 32/63.

Look closely: 31 = 0.00793651 but 32 = 0.00793654. Mathematically both are 1/126. They differ in the last digits because 32/63 and 31/63 round differently in 32-bit float, and those rounding errors are visible once you subtract 0.5 (subtracting nearly equal numbers amplifies relative error).

7.3 Why floats differ between CPU and GPU

A float stores about 7 significant decimal digits, so nearly every operation rounds. IEEE-754 guarantees that a single add, multiply, divide or sqrt is correctly rounded on both CPU and GPU. Results can still differ because:

• FMA: the GPU compiler fuses a*b + c into one instruction with one rounding instead of two. More accurate, but different.
• Operation order: floating-point addition isn't associative: (a + b) + c can differ from a + (b + c). Parallel reductions sum in a different order from a CPU loop.
• Fast math functions: with --use_fast_math, or with intrinsics like __sinf, results are approximate.
• CPU precision: the CPU compiler may keep intermediates in higher precision or vectorize differently.

So validate GPU results with a tolerance: fabs(gpu − cpu) <= 1e-6 * fabs(cpu) or similar, never with ==.

Exercises

Profiling: Nsight Systems, Nsight Compute and occupancy

Measure instead of guessing, and know which tool answers which question.

• After this module you can: choose between nsys and ncu, read their key tables, add NVTX ranges, time code correctly, and compute theoretical occupancy by hand.

8.1 Two tools, two questions

Rule of thumb: start with nsys to find the biggest cost. Only if it's a kernel do you dig in with ncu. ncu replays each kernel several times to collect counters (the "3 passes" in Tutorial 1), so it's slow and its timings aren't meaningful.

8.2 Reading nsys output

• cuda_api_sum: CPU-side time in each CUDA call. Remember context creation inflates the first call (Module 5).
• cuda_gpu_kern_sum: GPU time per kernel: total, instances, average, min, max. You'll use this table for Assignment 1 Q3 (QuEST's top kernels).
• cuda_gpu_mem_time_sum / _size_sum: time and bytes per copy direction. Divide bytes by time to get your effective PCIe bandwidth.

In Lecture 2's profile, the distance kernel took 2,464 ns and the copy back 1,920 ns, but cudaMalloc took 183 ms. Almost all of the program's time is start-up overhead, which tells you the problem is far too small for a GPU.

8.3 NVTX: labelling your code

#include <nvtx3/nvToolsExt.h>
nvtxRangePushA("H2D copy");
cudaMemcpy(d_a, h_a, bytes, cudaMemcpyHostToDevice);
nvtxRangePop();

Profile with nsys profile --trace=cuda,nvtx --stats=true ./app. The nvtx_sum table and the timeline then show your named phases. Ranges are measured on the CPU, so a range around a kernel launch alone measures only the launch, the same trap as with timers.

8.4 Timing correctly

For the assignments' stacked bar charts you need CPU timers around each phase. Two rules: synchronize before reading the timer, and exclude one-off start-up costs.

double now_ms() { using namespace std::chrono;
  return duration<double, std::milli>(steady_clock::now().time_since_epoch()).count(); }

cudaFree(0);                         // create the context before timing
double t0 = now_ms();
cudaMemcpy(d_a, h_a, bytes, cudaMemcpyHostToDevice);
double t1 = now_ms();
kernel<<<g, b>>>(…);
cudaDeviceSynchronize();             // without this you time only the launch
double t2 = now_ms();

For more reliable numbers, run each size a few times and report the median; the first launch of a kernel is often slower (module loading).

8.5 Occupancy

Theoretical occupancy is the best the resource limits allow. An SM's resident blocks are limited by whichever runs out first:

1. Warp slots: max warps per SM ÷ warps per block.
2. Block slots: max blocks per SM (16 on a T4, 32 on an H100).
3. Registers: 65,536 ÷ (registers per thread × threads per block).
4. Shared memory: shared memory per SM ÷ shared memory per block.

Achieved occupancy is what was measured while the kernel ran. It's lower than theoretical when there aren't enough blocks to fill the GPU, when blocks finish unevenly, or when the kernel is so short that start-up and tail dominate.

Exercises

First parallel patterns: map, stencil, matrix multiply

Three shapes of computation that cover most of what you'll write in this course.

• After this module you can: recognise map, stencil and matmul patterns, write each as a kernel with correct boundaries, explain why stencils need two buffers, and count the reads that later optimizations remove.

9.1 Map: each output depends on one input

Vector add, scaling and the distance kernel are all maps: output i depends only on input i. One thread per element, no communication. Maps are always memory-bound unless the per-element work is heavy.

9.2 Stencil: each output depends on its neighbours

The 1D heat equation models heat diffusing along a rod. Each step, a point moves toward the average of its neighbours:

__global__ void heatStep(const float *u, float *u_new, int n, float alpha) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i >= n) return;
  if (i == 0 || i == n - 1) { u_new[i] = u[i]; return; }    // fixed boundary values
  u_new[i] = u[i] + alpha * (u[i-1] - 2.0f * u[i] + u[i+1]);
}
// host: iterate, swapping buffers between steps
for (int t = 0; t < steps; t++) {
  heatStep<<<blocks, TPB>>>(d_u, d_unew, n, 0.25f);
  float *tmp = d_u; d_u = d_unew; d_unew = tmp;              // swap pointers, not data
}

Three things matter:

• Boundaries: thread 0 has no u[-1]. Decide a rule (fixed values here) and guard it.
• Two buffers: if a thread wrote into u while its neighbour still needs the old value, the result would depend on which thread runs first, and threads have no guaranteed order. Read from one array, write to another.
• Synchronizing between steps: threads in different blocks can't sync with each other inside a kernel, so each time step is a separate launch. Launches in the same stream run in order, so step t+1 always sees step t's results.

The same pattern is a 1D convolution with mask (α, 1−2α, α). General convolution uses a mask of any width.

9.3 2D Jacobi: the 5-point stencil

__global__ void jacobi(const float *in, float *out, int W, int H) {
  int c = blockIdx.x * blockDim.x + threadIdx.x;
  int r = blockIdx.y * blockDim.y + threadIdx.y;
  if (c < 1 || r < 1 || c >= W - 1 || r >= H - 1) return;   // skip boundary
  out[r*W + c] = 0.25f * (in[(r-1)*W + c] + in[(r+1)*W + c]
                        + in[r*W + c - 1] + in[r*W + c + 1]);
}

Each input value is read by 4 different threads (its 4 neighbours). That reuse is wasted when every read goes to global memory. Loading a tile into shared memory once and reading it 4 times from there is the classic optimization in a later lecture.

9.4 Matrix multiply

C (M×P) = A (M×N) · B (N×P). Element C[r][c] is the dot product of row r of A and column c of B. One thread per element of C, with a 2D grid as in Module 6:

__global__ void matmul(const float *A, const float *B, float *C, int M, int N, int P) {
  int c = blockIdx.x * blockDim.x + threadIdx.x;
  int r = blockIdx.y * blockDim.y + threadIdx.y;
  if (r >= M || c >= P) return;
  float sum = 0.0f;
  for (int k = 0; k < N; k++)
    sum += A[r * N + k] * B[k * P + c];
  C[r * P + c] = sum;
}

• FLOPs: each of the M·P outputs does N multiplies and N adds: 2·M·N·P in total.
• Global reads: each output reads N elements of A and N of B: 2·M·N·P reads, plus M·P writes.
• Reuse: each element of A is read P times (by every thread in its row of C) and each element of B M times. A only has M·N distinct elements, so almost all reads are repeats. Tiling in shared memory removes them.
• Access patterns: within a warp, B[k*P + c] is coalesced (c varies by lane), and A[r*N + k] is the same address for all lanes in a row, a broadcast. Both are efficient; the issue is the volume of reads, not their pattern.

Exercises

int n = 1024, steps = 1000, TPB = 256, blocks = (n + TPB - 1) / TPB;
std::vector<float> h(n, 0.0f);
for (int i = n / 2 - 5; i < n / 2 + 5; i++) h[i] = 1.0f;
float *d_u, *d_unew;
cudaMalloc(&d_u, n * 4); cudaMalloc(&d_unew, n * 4);
cudaMemcpy(d_u, h.data(), n * 4, cudaMemcpyHostToDevice);
cudaMemcpy(d_unew, h.data(), n * 4, cudaMemcpyHostToDevice);   // boundaries too
for (int t = 1; t <= steps; t++) {
  heatStep<<<blocks, TPB>>>(d_u, d_unew, n, 0.25f);
  std::swap(d_u, d_unew);
  if (t % 250 == 0) {
    cudaMemcpy(h.data(), d_u, n * 4, cudaMemcpyDeviceToHost);
    float mx = 0, s = 0; for (float v : h) { mx = fmaxf(mx, v); s += v; }
    printf("step %4d  max %.4f  sum %.2f\n", t, mx, s);
  }
}

Why parallel, what speed-ups are realistic, and why it's hard

The textbook's big-picture argument, which the course assumes you know.

• After this module you can: explain why the industry moved to parallel hardware, use the book's Amdahl examples, describe the "peach" model of an application, name the four challenges of parallel programming, and compare CUDA with MPI, OpenMP, OpenACC and OpenCL.

10.1 The end of the free lunch (§1.1)

Until about 2003, each new CPU ran existing programs faster because clock speeds kept rising. Then power and heat stopped that: a faster clock costs disproportionately more energy. Since then, vendors add more cores instead. A sequential program no longer gets faster on new hardware; only parallel programs do. The book calls this the concurrency revolution, and it describes two design paths:

• Multicore (CPUs): a few cores, each keeping the full speed of a sequential processor.
• Many-thread (GPUs): huge numbers of simple threads, optimized for total throughput (Module 1).

The book's main motivation for learning CUDA is scalability: write code that keeps getting faster as future GPUs add more SMs.

10.2 Why more speed? (§1.3)

The book argues that tomorrow's mass-market applications are today's "super-applications": molecular simulation in biology and medicine, high-resolution video processing, natural user interfaces, physics-based games. They all simulate a concurrent physical world with huge amounts of data, so they're naturally parallel. Its two rules of thumb:

• An application with data parallelism can often get a 10× speed-up with a few hours of work.
• A good, carefully optimized GPU implementation can reach more than 100× over one CPU core.

This is directly relevant to your programme: molecular and medical simulation is one of the book's headline examples.

10.3 Speeding up real applications (§1.4)

The book makes Amdahl's law (Module 4) concrete with two examples. Check them yourself:

Three further points from this section:

• 100×+ needs more than 99.9% of the time in parallel code, and that usually requires improving the algorithm, not just porting it.
• Straightforward ports usually saturate DRAM bandwidth and stop at around 10×. Going further means using on-chip memories (shared memory, registers) to cut DRAM accesses, then working around their limited capacity. This is the arc of the whole course.
• Give the CPU a fair chance. Some code runs very well on CPUs. A GPU speed-up measured against a poor, single-threaded CPU baseline is misleading. Keep this in mind when you report speed-ups in your project.

10.4 Four challenges of parallel programming (§1.5)

"If you don't care about performance, parallel programming is very easy." The hard part is speed, and the book names four reasons:

1. Algorithmic complexity. Some parallel algorithms do much more total work than the sequential version (they aren't work-efficient), so they can be slower on large inputs. Example: a naive parallel prefix sum does O(n log n) work against O(n) sequentially.
2. Memory-bound execution. Many applications are limited by memory access speed rather than by instructions (Module 4), which calls for new ways of improving data access.
3. Input sensitivity. Parallel performance varies much more with the characteristics of the input: irregular sizes, sparse structure, unpredictable data rates. This causes load imbalance and warp divergence.
4. Recurrences. Problems naturally written as "step n depends on step n−1" need non-intuitive reformulations to parallelize, sometimes with redundant work.

The good news: researchers have solved most of these for common patterns (reduction, scan, stencil, histogram, sparse matrix-vector product, and so on), and the rest of the book and course teaches them pattern by pattern.

10.5 Programming models compared (§1.6)

The book also mentions unified address spaces between CPU and GPU (the idea behind CUDA's Unified Memory, one of your possible project topics), which remove much of the manual copying.

10.6 The book's three goals (§1.7)

1. High performance through "computational thinking": reasoning about problems in ways that suit parallel hardware. You need hardware knowledge, taught as part of programming techniques.
2. Correctness and reliability: fast code you can debug and support, using simple barrier synchronization, memory consistency and atomics.
3. Scalability across future hardware: the key is to regularize and localize memory accesses, minimizing use of critical resources and conflicts when updating data.

Exercises

Oral exam drill

Answer each question out loud in about a minute, with no notes. If you hesitate, go back to the module in brackets. Repeat weekly until all of them are easy.

Concepts

Your Assignment 1 code

Big picture

Memory spaces and coalesced access

Why the order in which threads touch memory can change speed by 5× or more.

• After this module you can: say where each kind of variable lives, explain DRAM bursts, decide whether a warp's access is coalesced, count the memory transactions it causes, and confirm it in Nsight Compute.

12.1 Where does each variable live?

The lecture's check question: in distanceKernel, the local variable x is in a register. Each thread has its own copy and no other thread can see it. d_out is different: the pointer is passed by value into a register in each thread, but the array it points to is in global memory, visible to all threads.

12.2 How DRAM is read: bursts

DRAM stores each bit in a tiny capacitor, and reading a cell is slow: select a bank, open a row, select a column. The DRAM core runs at only a fraction (1/4 or 1/8) of the interface speed. To compensate, every access reads a whole burst of consecutive bytes into a buffer and sends them all. On a modern GPU the unit is a 32-byte sector, grouped into 128-byte lines.

The consequence: reading 4 bytes costs about the same as reading the 32 bytes around them. If threads use those neighbouring bytes, the transfer is fully used. If they don't, most of the bandwidth is wasted.

12.3 Coalescing

All 32 threads of a warp execute a load instruction together. The hardware collects their 32 addresses and issues as few sector requests as possible. When the threads access consecutive addresses, typically array[base + threadIdx.x], the 32 floats (128 bytes) fit in 4 sectors. That's a coalesced access.

When the addresses are spread out, the warp needs many sectors and uses only 4 bytes of each. Examples: a[threadIdx.x * stride], a column of a row-major matrix, or an array of large structs.

12.4 Matrices: rows good, columns bad

For a row-major matrix, element (r, c) is at r * cols + c. Consecutive columns are adjacent in memory and consecutive rows are cols elements apart. So:

• If neighbouring threads (consecutive threadIdx.x) handle neighbouring columns, accesses coalesce.
• If they handle neighbouring rows, each thread is cols × 4 bytes from the next. In the demo that's 32 KB, a different sector for every thread.

This is the whole difference between the two kernels in 3-memory-coalescing.cu. They do the same copy; only the mapping of x and y changes. In naive matrix multiply (Module 9), threads in a warp read B[k*P + c], which is coalesced, and the same A[r*N + k], a broadcast. That's why naive matmul is still reasonable.

12.5 Book Chapter 5: two more levers

• Memory parallelism (§5.2). DRAM has several channels and banks that work in parallel. To use them all, many requests must be in flight at once, which is another reason to have many warps resident (Module 2's Little's law).
• Thread granularity (§5.5). Sometimes it pays to give each thread more work, e.g. two outputs instead of one, when neighbouring threads would otherwise load the same data twice. Fewer threads, less redundant traffic.

Exercises

cudaEvent_t s, e; cudaEventCreate(&s); cudaEventCreate(&e); float ms;
cudaEventRecord(s);
copyRowWise<<<gridRow, block>>>(d_in, d_out, ROWS, COLS);
cudaEventRecord(e); cudaEventSynchronize(e); cudaEventElapsedTime(&ms, s, e);
printf("row-wise    %.2f ms  %.0f GB/s\n", ms, matBytes / ms / 1e6);
// same for copyColumnWise

Shared memory, synchronization and atomics

A fast scratchpad that threads of a block share, and the rules that keep its use correct.

• After this module you can: declare static and dynamic shared arrays, use __syncthreads() correctly (and spot when it deadlocks), handle halo cells, explain race conditions and fix them with atomics, write a tree reduction, and explain bank conflicts.

13.1 Why shared memory?

In the 1D stencil out[i] = (in[i-1] - 2·in[i] + in[i+1]) / h², every input value is read from global memory by three threads. Shared memory removes that redundancy: the block loads its piece of the input once into on-chip memory, then all threads read their neighbours from there. The recipe from the lecture:

1. Break the global array into tiles, one per block.
2. Allocate a tile in shared memory.
3. Fill the tile collectively: each thread loads one (or a few) elements.
4. Synchronize, then compute using the tile.

13.2 Declaring shared arrays

// static: size known at compile time
__shared__ float s_in[TPB + 2 * RAD];

// dynamic: size given at launch, as the 3rd launch parameter (in bytes)
extern __shared__ float s_in[];
size_t smem = (TPB + 2 * RAD) * sizeof(float);
kernel<<<blocks, TPB, smem>>>(…);

• A block can normally use up to 48 KB. More needs an opt-in: cudaFuncSetAttribute(kernel, cudaFuncAttributeMaxDynamicSharedMemorySize, bytes), up to prop.sharedMemPerBlockOptin (64 KB on a T4, 227 KB on an H100).
• Several dynamic arrays share one extern buffer: declare one array and split it with pointers, e.g. float *sA = buf; float *sB = buf + tileA_elems;.
• Shared memory limits occupancy, like registers (Module 8). With 64 KB per SM on a T4 and 16 KB per block, at most 4 blocks fit per SM.

13.3 __syncthreads(): the block barrier

Threads of a block run in no particular order, warp by warp. After thread 5 writes s_in[5], thread 6 can't assume thread 7 has written s_in[7] yet. __syncthreads() makes every thread in the block wait until all threads of that block have reached it. Every write before the barrier is then visible to every thread after it.

13.4 Halo cells

A stencil of radius RAD needs RAD extra values on each side of the tile, called halo (or ghost) cells. The tile holds blockDim.x + 2·RAD values, and thread t stores its own value at local index s_idx = t + RAD. The first RAD threads also load the left and right halos. Here's a version that is safe at every boundary:

#define RAD 1
__global__ void ddShared(float *out, const float *in, int n, float h) {
  extern __shared__ float s[];
  int i = blockIdx.x * blockDim.x + threadIdx.x;   // global index
  int s_idx = threadIdx.x + RAD;                   // local index
  s[s_idx] = (i < n) ? in[i] : 0.0f;               // body (0 if past the end)
  if (threadIdx.x < RAD) {                         // halos
    int L = i - RAD, R = i + blockDim.x;
    s[s_idx - RAD]        = (L >= 0) ? in[L] : 0.0f;
    s[s_idx + blockDim.x] = (R < n)  ? in[R] : 0.0f;
  }
  __syncthreads();                                 // ALL threads reach this
  if (i > 0 && i < n - 1)                           // guard only the write
    out[i] = (s[s_idx - 1] - 2.f * s[s_idx] + s[s_idx + 1]) / (h * h);
}
// launch: ddShared<<<blocks, TPB, (TPB + 2*RAD) * sizeof(float)>>>(…)

13.5 Race conditions and atomics

The dot product example (main.cpp, kernel.cu) first has each block sum its own products in shared memory, then adds the block sum to one global result. With ATOMIC 0 it does *d_res += blockSum;. That single line is really three steps: read d_res, add, write. When two blocks do this at once:

One update is lost. This is a race condition: the result depends on timing, so it can differ between runs and even be right sometimes, which makes it dangerous. atomicAdd(d_res, blockSum) does the read-add-write as one indivisible operation, so no update is lost. Atomics exist for add, sub, min, max, exch, CAS and bitwise operations, on int, float and (since Pascal) double.

Cost: atomics on the same address are serialized. One atomic per block (as here) is cheap; one per thread on a single counter can be very slow. That's why we reduce within the block first.

13.6 Doing the block sum in parallel: tree reduction

In kernel.cu, thread 0 adds all 64 products alone while 63 threads wait, which is O(TPB) serial steps. A tree reduction halves the active threads each step, so it finishes in log₂(TPB) steps:

__shared__ int s[TPB];                                  // TPB a power of 2
s[tid] = (idx < n) ? a[idx] * b[idx] : 0;              // 0 for threads past n
__syncthreads();
for (int stride = blockDim.x / 2; stride > 0; stride /= 2) {
  if (tid < stride) s[tid] += s[tid + stride];
  __syncthreads();                                      // outside the if: every thread reaches it
}
if (tid == 0) atomicAdd(res, s[0]);

With stride halving from the top, active threads are always the contiguous tid < stride, so whole warps either work or idle. That avoids divergence until stride < 32.

13.7 Bank conflicts

Shared memory is divided into 32 banks, each 4 bytes wide; element s[i] (float) is in bank i % 32. If the 32 threads of a warp access 32 different banks (or the same word, which broadcasts), the access takes one cycle. If k threads hit different words in the same bank, it takes k cycles: a k-way bank conflict.

• s[threadIdx.x]: conflict-free.
• s[threadIdx.x * 2]: 2-way conflict.
• Column of a 32×32 tile, tile[threadIdx.x][c]: every element is in the same bank, a 32-way conflict. The classic fix is padding: __shared__ float tile[32][33]; shifts each row by one bank.

The 2D Shared Memory tutorial (needed for Assignment 3) explores exactly this with setRowReadRow (conflict-free) against setColReadCol (conflicts), using a 16 × 32 tile. Nsight Compute reports conflicts under Memory Workload Analysis.

Exercises

__global__ void dotKernel(int *d_res, const int *d_a, const int *d_b, int n) {
  __shared__ int s[TPB];
  int idx = blockIdx.x * blockDim.x + threadIdx.x, tid = threadIdx.x;
  s[tid] = (idx < n) ? d_a[idx] * d_b[idx] : 0;
  __syncthreads();
  for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) {
    if (tid < stride) s[tid] += s[tid + stride];
    __syncthreads();
  }
  if (tid == 0) atomicAdd(d_res, s[0]);
}

Tiling: convolution and matrix multiplication

Load data into shared memory once and reuse it many times. This is the most important optimization in the course.

• After this module you can: write basic and tiled 1D convolution, derive global read counts for both, use constant memory for a mask, explain the tiled matrix multiply phase by phase, derive its reduction in global reads, handle sizes that aren't multiples of the tile, and reason about choosing a tile size.

14.1 1D convolution

Convolution slides a small mask M of width MW = 2r + 1 over an input array and computes a weighted sum at every position:

It's used for signal and audio filtering, image blurring and edge detection, CNN layers, and stencils in simulations such as the heat equation (whose mask is [α, 1−2α, α], Module 9).

Basic kernel: one thread per output, reading MW inputs straight from global memory:

#define R 1                                     // example: mask width 3 (the assignment uses 5)
__constant__ float d_M[2 * R + 1];              // mask in constant memory
__global__ void conv_basic(const float *N, float *P, int n) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i >= n) return;                           // fine here: no __syncthreads
  float sum = 0.0f;
  for (int j = -R; j <= R; j++) {
    int k = i + j;
    if (k >= 0 && k < n) sum += d_M[j + R] * N[k];   // out-of-range inputs count as 0
  }
  P[i] = sum;
}
// host: cudaMemcpyToSymbol(d_M, h_M, sizeof(h_M));

Why constant memory for the mask? Every thread reads the same d_M[j] at the same time. Constant memory has a cache that broadcasts one value to the whole warp, so the mask costs almost no global traffic. Mention this in your report: the mask's reads are separate from the input's.

14.2 Tiled 1D convolution

Each input value is needed by MW different outputs. In the tiled kernel, a block loads a tile of TILE inputs plus r halo cells on each side into shared memory once, then each thread computes from shared memory. This is Module 13's halo pattern with the radius as a parameter:

__global__ void conv_tiled(const float *N, float *P, int n) {
  extern __shared__ float s[];                      // TILE + 2R floats
  int i = blockIdx.x * blockDim.x + threadIdx.x, t = threadIdx.x;
  s[t + R] = (i < n) ? N[i] : 0.0f;
  if (t < R) {                                       // left and right halos
    int L = i - R, Rt = i + blockDim.x;
    s[t]                  = (L >= 0) ? N[L]  : 0.0f;
    s[t + blockDim.x + R] = (Rt < n) ? N[Rt] : 0.0f;
  }
  __syncthreads();
  if (i < n) {
    float sum = 0.0f;
    for (int j = 0; j < 2 * R + 1; j++) sum += d_M[j] * s[t + j];
    P[i] = sum;
  }
}
// launch: conv_tiled<<<(n+TILE-1)/TILE, TILE, (TILE + 2*R)*sizeof(float)>>>(d_N, d_P, n);

Counting global memory accesses

Ignoring the mask (constant cache) and the edges of the whole array:

The reduction factor is about MW / (1 + 2r/TILE). It approaches MW for large tiles, because the halo overhead 2r/TILE shrinks. At the array ends, inputs outside the array aren't read at all (they're 0), so exact counts are slightly lower. A careful report states whether it counts those.

Choosing the tile size

• TILE = threads per block here, so it must be ≤ 1024 and should be a multiple of 32.
• Bigger tiles mean less halo overhead (2r/TILE) but fewer blocks. For n = 1024 and TILE = 1024 there's only one block, using 1 of 40 SMs.
• Shared memory per block is (TILE + 2r) × 4 bytes, tiny for 1D, so it doesn't limit occupancy here. That matters for the assignment's question about whether more shared memory helps.
• Decide by experiment: sweep TILE ∈ {32, 64, 128, 256, 512, 1024}, time each with events, and explain the best one with occupancy and block count.

14.3 Tiled matrix multiplication

Naive matmul (Module 9) reads every element of A P times and every element of B M times from global memory. Tiling exploits that a block of threads computing a T × T tile of C needs the same rows of A and columns of B.

The kernel runs in phases. In phase p, each thread of the block loads one element of the A tile and one of the B tile into shared memory, the block synchronizes, every thread does T multiply-adds from shared memory, and the block synchronizes again before the tiles are overwritten:

#define T 16                                        // square tile, fixed at compile time
__global__ void matmul_tiled(const float *A, const float *B, float *C, int M, int N, int P) {
  __shared__ float sA[T][T], sB[T][T];
  int tx = threadIdx.x, ty = threadIdx.y;
  int col = blockIdx.x * T + tx, row = blockIdx.y * T + ty;
  float sum = 0.0f;
  for (int p = 0; p < (N + T - 1) / T; p++) {       // phases along the shared dimension
    int aCol = p * T + tx, bRow = p * T + ty;
    sA[ty][tx] = (row < M && aCol < N) ? A[row * N + aCol] : 0.0f;   // 0-padding at edges
    sB[ty][tx] = (bRow < N && col < P) ? B[bRow * P + col] : 0.0f;
    __syncthreads();                                // tiles fully loaded
    for (int k = 0; k < T; k++) sum += sA[ty][k] * sB[k][tx];
    __syncthreads();                                // everyone done before overwrite
  }
  if (row < M && col < P) C[row * P + col] = sum;
}
dim3 block(T, T), grid((P + T - 1) / T, (M + T - 1) / T);

Points to understand line by line:

• Coalesced loads (corner turning). Both loads index the contiguous dimension with tx: A[row*N + p*T + tx] and B[(p*T+ty)*P + col]. Consecutive threads read consecutive addresses. The lecture calls this corner turning: the load runs along rows even though the computation uses a column of B.
• Padding with zeros handles M, N, P that aren't multiples of T. A 0 contributes nothing to the sum. There's no early return, because of the barriers.
• Two barriers per phase: the first so no one reads a tile before it's complete, the second so no one overwrites it while others still read.
• Shared-memory reads: sA[ty][k] is the same for all threads in a row (broadcast); sB[k][tx] spans consecutive banks. No conflicts.

Counting global reads

Each element of C still needs N products, but each thread now loads only one A and one B element per phase, and there are N/T phases. So each thread loads 2N/T elements instead of 2N:

With T = 16 that's 16× fewer global reads, and the arithmetic intensity grows by the same factor (Module 4). For non-square tiles, work out separately how many times each A and B element is loaded: each A element is loaded once per block column of the grid, each B element once per block row. That reasoning gives you the formula for a tileX × tileY kernel.

Tile size: the trade-offs

Bigger tiles aren't automatically faster. Once memory traffic is low enough, the limit moves to shared-memory bandwidth, barrier stalls and occupancy, and the best T has to be measured.

Runtime tile sizes: the assignment's tiled_gemm() takes tileX and tileY as parameters. A static __shared__ float sA[T][T] needs T at compile time, so you'll need dynamic shared memory: one extern buffer split into two arrays, indexed by hand, with the size passed at launch (Module 13.2). Think carefully about which dimension of the A tile and of the B tile each of tileX and tileY controls, and whether the shared dimension's tile width must match between them.

Exercises

Porting a CPU hotspot to the GPU

The workflow you'll use in Q3, the bonus and your llm.c project: measure, port one function correctly, find the new bottleneck, fix it, measure again.

• After this module you can: find and count a hotspot, port a matrix-vector product with a coalesced design, verify a port inside a larger program, explain why per-call allocation and copying kill performance, choose an optimization from an nsys profile, and plan the port of a particle mover (structs, deep copies, validation).

15.1 The workflow

1. Measure the baseline with the program's own metric (tokens/s, total simulation time) and a profiler.
2. Find and count the hotspot. Add a static counter to the function or use perf/gprof. Knowing the call count and per-call size tells you whether launch overhead will matter.
3. Port naively but correctly, with data allocated and copied inside the function. It's slow, but it isolates correctness.
4. Verify against the CPU version.
5. Profile with nsys: is time in kernels, in cudaMemcpy, in cudaMalloc, or in gaps between launches?
6. Fix the biggest cost, one change at a time, and re-profile to show which metrics changed.

15.2 Matrix-vector products in an LLM

In llama2.c, matmul(xout, x, w, n, d) computes xout = W · x, where W is d × n (row-major) and x has n entries. Generating each token calls it several times per transformer layer (the attention projections and the feed-forward layers) plus once for the output classifier. Count it with a static counter, don't guess. Unlike matrix-matrix multiply, a matrix-vector product does just 2 FLOPs per weight element loaded, so it's memory-bound (AI ≈ 0.5 FLOP/byte). Its speed is set by how fast you can stream W.

Design 1: one thread per output row

__global__ void matvec_row(const float *W, const float *x, float *y, int n, int d) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;          // output row
  if (i >= d) return;
  float s = 0.0f;
  for (int j = 0; j < n; j++) s += W[(size_t)i * n + j] * x[j];
  y[i] = s;
}

Simple and correct, but check the access to W with Module 12: at each step j, consecutive threads read W[i*n + j], which are n floats apart, so the access is uncoalesced. (x[j] is a broadcast, which is fine.) For d = 288 rows there are also only 288 threads, which won't fill a GPU.

Design 2: one warp per output row

__global__ void matvec_warp(const float *W, const float *x, float *y, int n, int d) {
  int warp = (blockIdx.x * blockDim.x + threadIdx.x) / 32, lane = threadIdx.x % 32;
  if (warp >= d) return;                                  // whole warp exits together
  float s = 0.0f;
  for (int j = lane; j < n; j += 32)                      // lanes read consecutive W entries
    s += W[(size_t)warp * n + j] * x[j];
  for (int off = 16; off > 0; off /= 2)                   // warp reduction in registers
    s += __shfl_down_sync(0xffffffff, s, off);
  if (lane == 0) y[warp] = s;
}
// launch: d warps → blocks = (d * 32 + 255) / 256, 256 threads

Now the 32 lanes read 32 consecutive floats of a row (coalesced), and there are 32× more threads. __shfl_down_sync lets threads of a warp read each other's registers, so the warp sums its 32 partial results in 5 steps with no shared memory and no barrier. Comparing designs 1 and 2 with ncu (sectors per request, DRAM throughput) makes a strong Module 12 story for your report.

A third option is the library: cublasSgemv. cuBLAS assumes column-major storage, so a row-major W (d × n) looks like a column-major n × d matrix, and you pass CUBLAS_OP_T. Library calls are a legitimate optimization and a project topic.

15.3 Why the naive port is slow

A naive port does cudaMalloc → copy W → copy x → kernel → copy y back → cudaFree inside every call. Per call:

• W is the same every token (model weights), yet it's re-copied each call. Copying d·n·4 bytes over PCIe takes longer than the GPU needs to read it.
• cudaMalloc/cudaFree are slow API calls (tens to hundreds of µs) that also synchronize.
• Each call is small, so fixed costs (launch ~5 µs, copy latency ~10 µs) dominate.

So the GPU version can easily be slower than the CPU at first. That's expected and is the point of the exercise. nsys's cuda_api_sum and cuda_gpu_mem_time_sum tables show exactly which of these costs dominates.

Optimizations to choose from (the assignment asks for one)

Verifying inside a larger program

• Temporarily run both versions each call and compare xout elementwise with a relative tolerance (reductions sum in a different order, Module 7).
• Run with temperature 0 (greedy sampling; llama2.c's -t 0) on CPU and GPU builds. The generated text should be identical, or diverge only after many tokens because of float rounding. Explain this in the report.

15.4 Bonus preparation: porting a particle mover

The iPIC3D-mini bonus asks you to port mover_PC(): for every particle, interpolate the electromagnetic field at its position and push it for several sub-cycles. The structure suits a GPU well: particles are independent, so use one thread per particle, keeping its whole sub-cycling loop inside the thread. What makes it harder than vector add:

• Structs with pointers. Particle and field data live in structs containing pointers to arrays. Copying the struct to the GPU copies only the host pointers, which are invalid there. You need a deep copy: cudaMalloc each array, copy its contents, and pass device pointers (in a device-side struct or as kernel arguments).
• Multi-dimensional arrays. Fields may be stored as double*** on the CPU. On the GPU, use the flat contiguous array behind them (check the code for a "flat" version) and compute the index by hand.
• Data layout. Separate arrays per component (x[], y[], …), a structure of arrays, coalesce. An array of structs doesn't (Module 12 exercise).
• Double precision. On a T4, FP64 runs at 1/32 of FP32, so don't expect a huge speed-up there; the cluster's H100 is much better at FP64. Mention which GPU you used and why it matters.
• Validation. Run CPU and GPU versions on the same input, compare the output files (positions, velocities, or the code's diagnostic outputs) with a tolerance, and report the maximum relative difference.
• Timing. Report the mover's time with and without allocation and transfers, since the bonus allows them inside mover_PC_gpu(). The difference shows how much a fuller port, keeping data resident, would gain.

Exercises

Streams: overlapping copies and computation

Keep the PCIe link and the GPU busy at the same time, instead of taking turns.

• After this module you can: explain what a stream is and how the default stream behaves, use pinned memory and cudaMemcpyAsync, pipeline a computation in segments, synchronize streams correctly, time with CUDA events, and predict the speed-up from overlap.

16.1 The problem: everything takes turns

So far, each program copies A in, copies B in, computes, then copies C out, one after another. While data crosses PCIe the GPU is idle, and while the GPU computes the PCIe link is idle. For vector add (Module 4) the copies are about 95% of the time, so the GPU sits idle almost always.

The fix is pipelining: split the arrays into segments. While segment 1 is being computed, segment 2 is being copied in and segment 0's result copied out. That's the slides' A.0/A.1/A.2 picture.

16.2 Streams and the default stream

A stream is a queue of GPU operations (copies and kernels). Operations in the same stream run in order. Operations in different streams may run at the same time. This is task parallelism on top of the data parallelism inside each kernel.

If you don't name a stream, you use the default (null) stream, which is special:

• No operation in the default stream starts until all previously issued operations in any stream have finished.
• An operation in the default stream must finish before any later operation in any stream starts.

So one default-stream call in the middle of your pipeline serializes everything. Since CUDA 7 you can compile with --default-stream per-thread, which gives each host thread its own default stream that behaves like a normal stream.

16.3 The API

cudaStream_t s[NS];
for (int k = 0; k < NS; k++) cudaStreamCreate(&s[k]);

// pinned (page-locked) host memory: REQUIRED for copies to overlap
float *h_a; cudaMallocHost(&h_a, bytes);       // free with cudaFreeHost

cudaMemcpyAsync(d_a, h_a, bytes, cudaMemcpyHostToDevice, s[k]);   // 5th arg: stream
kernel<<<grid, block, 0, s[k]>>>(…);                              // 4th launch arg: stream
cudaMemcpyAsync(h_c, d_c, bytes, cudaMemcpyDeviceToHost, s[k]);

cudaStreamSynchronize(s[k]);   // wait for stream k only
cudaDeviceSynchronize();       // wait for all streams
cudaStreamDestroy(s[k]);

The two requirements from the slides: the kernel and the copy to be overlapped must be in different non-default streams, and the host memory must be pinned. Hardware matters too: most GPUs have separate copy engines for each direction, so H2D, D2H and a kernel can all run at once. deviceQuery reports this as asyncEngineCount.

16.4 Pipelining vector add

int seg = (n + NS - 1) / NS;                    // could also use many more segments than streams
for (int k = 0; k < NS; k++) {
  int off = k * seg, len = min(seg, n - off);
  size_t b = len * sizeof(float);
  cudaMemcpyAsync(d_a + off, h_a + off, b, cudaMemcpyHostToDevice, s[k]);
  cudaMemcpyAsync(d_b + off, h_b + off, b, cudaMemcpyHostToDevice, s[k]);
  vecAdd<<<(len + 255) / 256, 256, 0, s[k]>>>(d_a + off, d_b + off, d_c + off, len);
  cudaMemcpyAsync(h_c + off, d_c + off, b, cudaMemcpyDeviceToHost, s[k]);
}
cudaDeviceSynchronize();

Note the pointer offsets (d_a + off): each stream works on its own slice, so there's no conflict. Issue order matters on some GPUs. Interleaving per segment (copy, kernel, copy for segment 0, then segment 1 …), as above, works well on modern hardware; older GPUs sometimes needed all H2D copies issued first. The nsys timeline shows which you got.

16.5 Timing with CUDA events

CPU timers measure when the host saw things happen. With asynchronous work, events are better: they're recorded in a stream and timestamped by the GPU when they're reached.

cudaEvent_t start, stop; cudaEventCreate(&start); cudaEventCreate(&stop);
cudaEventRecord(start, s[0]);
/* … work … */
cudaEventRecord(stop, s[0]);
cudaEventSynchronize(stop);
float ms; cudaEventElapsedTime(&ms, start, stop);

Events also create dependencies between streams: cudaStreamWaitEvent(s2, ev) makes stream 2 wait until ev (recorded in stream 1) is reached, without blocking the CPU.

Exercises

Scheduling, synchronization and CUDA Graphs

What the hardware does with your blocks, which kinds of synchronization exist, and how to stop launch overhead from dominating.

• After this module you can: describe how blocks and warps are scheduled, choose the right synchronization level (warp, block, grid, stream, device), recognize a launch-overhead-bound program from an nsys profile, and use CUDA Graphs to remove that overhead.

17.1 How work is scheduled

• Blocks → SMs. The block scheduler hands blocks to SMs as resources free up (Module 8's limits). When a block finishes, a waiting block takes its place. There's no guaranteed order, and a block never migrates.
• Waves. If a grid has more blocks than fit at once, it runs in waves. With 40 SMs × 4 blocks per SM = 160 resident blocks, a 170-block grid runs one full wave and then a second wave of 10 blocks with most of the GPU idle: the tail effect. Grids much larger than one wave make tails negligible.
• Warps → schedulers. Inside an SM, each of the 4 warp schedulers picks a ready warp every cycle (Module 2). Nsight Compute's Scheduler Statistics and Warp State sections show why warps aren't ready (waiting on memory, at a barrier, …).

17.2 The synchronization ladder

Rule of thumb: use the narrowest synchronization that is correct. Every wider one stops more work.

17.3 Launch overhead

Every kernel launch costs a few microseconds of CPU time (cudaLaunchKernel), and back-to-back kernels leave small gaps on the GPU. That's negligible for a 10 ms kernel and dominant for a 2 µs one.

17.4 CUDA Graphs

A CUDA Graph records a sequence of GPU operations (kernels, copies, dependencies) once, then replays the whole sequence with one launch. The driver validates and prepares everything at creation time, so each replay costs about as much as a single launch, and the GPU runs the kernels back to back. The easiest way to build one is stream capture:

cudaStream_t s; cudaStreamCreate(&s);
cudaGraph_t graph; cudaGraphExec_t exec;

cudaStreamBeginCapture(s, cudaStreamCaptureModeGlobal);
for (int k = 0; k < 100; k++) {                     // record 100 steps (an even number,
  heatStep<<<blocks, threads, 0, s>>>(d_u, d_unew, N, r);   // so the swaps end where they began)
  std::swap(d_u, d_unew);                           // host-side swap is baked into the captured args
}
cudaStreamEndCapture(s, &graph);
cudaGraphInstantiate(&exec, graph, 0);              // CUDA 12+ signature

for (int rep = 0; rep < steps / 100; rep++)
  cudaGraphLaunch(exec, s);                         // one launch = 100 kernels
cudaStreamSynchronize(s);

• Captured operations aren't executed during capture, only recorded. Kernel arguments are captured by value, so the pointer swaps are frozen into the graph, which is why an even number of steps per graph keeps d_u pointing at the latest result.
• Graphs suit work that repeats with the same structure: time-stepping, every training step of a neural network, every token in LLM inference. llm.c and PyTorch both use them.
• If parameters change between replays, update the graph with cudaGraphExecKernelNodeSetParams or re-capture.

Exercises

Unified Memory

One pointer that works on both CPU and GPU, what it really does underneath, and how to make it fast.

• After this module you can: allocate managed memory, explain on-demand page migration and why the first kernel is slow, use prefetching and memory advice, explain oversubscription, and design a fair Unified Memory experiment for the project.

18.1 One address space

Unified Memory gives you a single pointer valid on both the CPU and the GPU. Replace malloc + cudaMalloc + cudaMemcpy with one call:

float *x, *y;
cudaMallocManaged(&x, n * sizeof(float));
cudaMallocManaged(&y, n * sizeof(float));
for (int i = 0; i < n; i++) { x[i] = 1.0f; y[i] = 2.0f; }    // CPU writes
add<<<blocks, 256>>>(n, x, y);                                 // GPU reads/writes the same pointers
cudaDeviceSynchronize();                                      // REQUIRED before the CPU reads y
printf("%f\n", y[0]);                                          // CPU reads
cudaFree(x); cudaFree(y);

The synchronize is still essential: the launch is asynchronous, so without it the CPU might read y before the kernel has written it.

18.2 How it works (Pascal and newer, including T4 and H100)

• cudaMallocManaged reserves address space but may not allocate any physical memory yet. Pages are created on first touch.
• When the GPU accesses a page that lives in CPU memory, it takes a page fault. The driver migrates the page (typically 4 KB–2 MB chunks) over PCIe, updates the page tables, and the warp continues. The same happens in reverse when the CPU touches GPU-resident pages.
• Faults are handled in batches, but each batch costs tens of microseconds. A kernel touching gigabytes of non-resident data spends most of its time faulting, so the first kernel that touches managed data is often much slower than the same kernel with explicit copies.
• Pre-Pascal GPUs worked differently: all managed data was migrated to the GPU at every launch, with no faults. That's why the slides distinguish the generations.
• Oversubscription: on Pascal and newer, managed allocations can exceed GPU memory. Pages are evicted to the host as needed. It's slow, but it lets programs run that wouldn't fit otherwise.

18.3 Making it fast: prefetch and advise

int dev; cudaGetDevice(&dev);
cudaMemPrefetchAsync(x, bytes, dev, stream);             // move pages to the GPU before the kernel
kernel<<<g, b, 0, stream>>>(x, …);
cudaMemPrefetchAsync(y, bytes, cudaCpuDeviceId, stream); // bring results back before the CPU reads

cudaMemAdvise(w, bytes, cudaMemAdviseSetReadMostly, dev);  // read-only data: keep copies on both sides
cudaMemAdvise(x, bytes, cudaMemAdviseSetPreferredLocation, dev);

Prefetching turns many small fault-driven migrations into one bulk transfer, often matching explicit cudaMemcpy performance. SetReadMostly suits model weights or lookup tables: each processor keeps its own read-only copy. In CUDA 13 the prefetch and advise functions take a cudaMemLocation struct instead of a plain device ID. Check the signature in your toolkit's headers if a call won't compile.

18.4 When to use it

• Good for: complex data structures with pointers (linked structures, the particle structs of Module 15), where deep copies are painful; fast porting; data that doesn't fit in GPU memory; CPU and GPU taking turns on the same data.
• Costs: without prefetching, page faults can make it several times slower than explicit copies; behaviour is less predictable; it's harder to overlap with computation.
• On systems with a fast CPU–GPU link and shared page tables (such as Grace Hopper), the costs are much lower. Results depend on the platform, which is worth stating in the project.

Exercises

Sparse data on GPUs

Most values are zero, so store and compute only the rest, and deal with the irregularity that brings.

• After this module you can: store a matrix in CSR, write scalar and vector SpMV kernels and explain their coalescing and load balance, describe SpMM and embedding lookups, and explain 2:4 structured sparsity on Tensor Cores.

19.1 Why sparsity is different

A sparse matrix (finite-element meshes, graphs, recommendation data, pruned neural networks) is mostly zeros. Storing it densely wastes memory and computing on it wastes FLOPs. Sparse formats store only the non-zeros (nnz), but the access pattern becomes irregular: the book's challenge of "input sensitivity" (Module 10). Rows have different lengths, so threads get different amounts of work.

19.2 CSR: compressed sparse row

Three arrays: val[nnz] (non-zero values, row by row), col[nnz] (their column indices), and rowPtr[rows+1] (where each row starts in val and col).

// matrix        val    = [5, 1, 2, 3, 4]
// [5 0 0 1]     col    = [0, 3, 1, 2, 3]
// [0 2 0 0]     rowPtr = [0, 2, 3, 3, 5]     (row 2 is empty: rowPtr[2] == rowPtr[3])
// [0 0 0 0]
// [0 0 3 4]

Other formats trade differently: COO (row, col, value triples: simple, good for building), ELL (pad every row to the same length: regular and coalesced, but wasteful if lengths vary), and hybrids. cuSPARSE supports all of them.

19.3 SpMV: sparse matrix × dense vector

Scalar kernel, one thread per row:

__global__ void spmv_scalar(int rows, const int *rowPtr, const int *col,
                            const float *val, const float *x, float *y) {
  int r = blockIdx.x * blockDim.x + threadIdx.x;
  if (r >= rows) return;
  float s = 0.0f;
  for (int k = rowPtr[r]; k < rowPtr[r + 1]; k++) s += val[k] * x[col[k]];
  y[r] = s;
}

• Uncoalesced: neighbouring threads read val from different rows, far apart (the same problem as Module 15's thread-per-row matvec).
• Load imbalance and divergence: a thread with a 1,000-element row keeps its whole warp busy while threads with 3-element rows idle.
• Gather: x[col[k]] is an indirect, irregular read. Caches help if columns are clustered.

Vector kernel, one warp per row: lanes stride through the row's non-zeros (k = rowPtr[r] + lane; k < rowPtr[r+1]; k += 32), then do a shuffle reduction exactly as in Module 15. Reads of val and col are now coalesced. It's best when rows have ≳ 32 non-zeros; for very short rows most lanes idle, so libraries pick strategies based on the row-length distribution.

SpMV has very low arithmetic intensity: 2 FLOPs per non-zero for ~12 bytes read (value, index, x). It's always memory-bound, and its performance is judged in GB/s, not GFLOP/s.

19.4 SpMM and embedding lookups

• SpMM (sparse matrix × dense matrix): every non-zero is used for a whole row of the dense matrix, so each loaded index is reused many times. Intensity is higher than SpMV and coalescing is easier along the dense dimension. It's used in graph neural networks and pruned models.
• Embedding lookup (recommender systems, LLM token embeddings): gather rows of a large table by index, out[i] = table[ids[i]], often followed by a sum (EmbeddingBag). This is mathematically an SpMM with a one-hot sparse matrix, purely memory-bound. Assign a warp or block per looked-up row so that each row's elements are copied with coalesced accesses.

19.5 Structured sparsity on Tensor Cores

Unstructured sparsity is hard for hardware to exploit. Ampere and newer Tensor Cores support 2:4 structured sparsity: in every group of 4 consecutive weights, at most 2 are non-zero. The hardware stores only the 2 values plus 2-bit indices, and skips the zeros, up to 2× the Tensor Core throughput with half the weight storage.

• Networks are pruned to the 2:4 pattern and fine-tuned to recover accuracy, typically with little loss.
• Used through cuSPARSELt or framework support. It needs Ampere or newer: the H100 has it, the T4 (Turing) doesn't.
• It's the "structured sparsity" listed in project topic 2. Pair it with Module 20 (Tensor Cores).

Exercises

__global__ void spmv_vector(int rows, const int *rowPtr, const int *col,
                            const float *val, const float *x, float *y) {
  int r = (blockIdx.x * blockDim.x + threadIdx.x) / 32, lane = threadIdx.x % 32;
  if (r >= rows) return;
  float s = 0.0f;
  for (int k = rowPtr[r] + lane; k < rowPtr[r + 1]; k += 32) s += val[k] * x[col[k]];
  for (int off = 16; off > 0; off /= 2) s += __shfl_down_sync(0xffffffff, s, off);
  if (lane == 0) y[r] = s;
}

Tensor Cores and matrix cores

Hardware that does a whole small matrix multiply per instruction, and what it costs in precision.

• After this module you can: explain what a Tensor Core computes and why it's so fast, compute peak throughput from unit counts, program one with WMMA, summarize the course reading's findings on performance and precision (including its residual technique), and compare NVIDIA Tensor Cores with AMD matrix cores.

20.1 What a Tensor Core does

A normal CUDA core does one fused multiply-add (FMA) per cycle: d = a·b + c on scalars. Volta (V100, 2017) added Tensor Cores to each SM. A Tensor Core does

with A and B in half precision (FP16) and the accumulation C, D in FP16 or FP32. That's mixed precision: cheap 16-bit inputs, but products summed in 32 bits so rounding errors don't pile up. Each FMA rounds once instead of twice (Module 7), which also helps accuracy.

Why so fast? Matrix multiply has huge data reuse (Module 4). One instruction that fetches 32 inputs and does 64 FMAs needs far less control logic, register traffic and energy per FLOP than 64 separate FMA instructions. The silicon goes to arithmetic.

20.2 Peak throughput, from the slides

The V100 has 80 SMs; each SM has 4 processing blocks with 2 Tensor Cores each, so 8 per SM and 640 in total.

20.3 Three ways to program them (from the reading)

1. WMMA API (warp matrix multiply-accumulate): CUDA C++ calls where a whole warp cooperatively loads 16×16 tiles into opaque "fragments" and multiplies them.
2. CUTLASS: NVIDIA's C++ template library of GEMM building blocks (Module 22).
3. cuBLAS: call GEMM with Tensor Core math enabled (Module 22).

#include <mma.h>
using namespace nvcuda;
// one warp computes a 16x16 tile of C = A (16xK) * B (Kx16); half inputs, float accumulate
__global__ void wmma_tile(const half *A, const half *B, float *C, int K, int lda, int ldb, int ldc) {
  wmma::fragment<wmma::matrix_a, 16, 16, 16, half, wmma::row_major> a;
  wmma::fragment<wmma::matrix_b, 16, 16, 16, half, wmma::row_major> b;
  wmma::fragment<wmma::accumulator, 16, 16, 16, float> acc;
  wmma::fill_fragment(acc, 0.0f);
  for (int k = 0; k < K; k += 16) {
    wmma::load_matrix_sync(a, A + k, lda);          // whole warp loads the tile
    wmma::load_matrix_sync(b, B + k * ldb, ldb);
    wmma::mma_sync(acc, a, b, acc);                  // Tensor Core multiply-accumulate
  }
  wmma::store_matrix_sync(C, acc, ldc, wmma::mem_row_major);
}
// launch one warp: wmma_tile<<<1, 32>>>(…); a full GEMM tiles C across many warps and blocks

Note the shift in thinking: the unit of work is a warp, not a thread. Fragments are distributed across the 32 lanes in an undocumented layout, and the _sync functions must be called by all lanes. Dimensions must be multiples of 16, so pad with zeros otherwise.

20.4 What the reading found (Markidis et al., 2018, on a V100)

• Performance: up to 83 TFLOP/s in mixed precision with cuBLAS, about 7× single-precision and 3× half-precision (non-Tensor-Core) GEMM, but still only two thirds of the 125 TFLOP/s peak. A WMMA batched GEMM for many small matrices reached about 4 TFLOP/s. Tensor Cores sped up GEMM about 6× for large matrices and 2.5–12× for many small matrices in parallel.
• Programmability: WMMA gives direct control but is low-level; CUTLASS builds on it; cuBLAS is easiest and fastest for standard GEMMs.
• Precision: rounding FP32 inputs to FP16 loses accuracy, and the error grows with matrix size. That might be critical for HPC applications.
• Their fix, residual refinement: split each input into a half-precision part and the leftover, R_A = A_single − A_half (also stored in FP16). Then compute
  A·B ≈ A_h·B_h + A_h·R_B + R_A·B_h (dropping the tiny R_A·R_B)
  Each term runs on Tensor Cores. Three GEMMs instead of one buys back most of the precision, while still beating FP32 CUDA cores when Tensor Cores are ≥ 3× faster. You can refine one or both matrices depending on the accuracy you need.
• Outlook: codes built on GEMM benefit directly; others (e.g. particle codes) would need reformulating in matrix form.

20.5 AMD: MI100 and matrix cores

• AMD's data-centre line (CDNA architecture, Instinct MI100/MI200/MI300) powers Frontier, El Capitan, LUMI and KTH's own Dardel.
• Vocabulary: AMD's compute unit (CU) ≈ NVIDIA's SM; a wavefront is 64 threads (vs a 32-thread warp); MI100 has 120 CUs with HBM2.
• From the slides: MI100 offers 11.5 TFLOP/s FP64 and 23.1 TFLOP/s FP32, against the A100's 9.7 and 19.5 (non-Tensor-Core). It's competitive in raw performance, while NVIDIA's software ecosystem (CUDA, libraries, tools) is more mature.
• Matrix Core Engines are AMD's equivalent of Tensor Cores, programmed through rocBLAS, rocWMMA or compiler intrinsics. The course reading (Schieffer, De Medeiros, Faj, Marathe and Peng, ISPASS 2024) studies their performance, power efficiency and programmability.
• Portability: HIP is nearly source-compatible with CUDA (cudaMalloc → hipMalloc), and hipify tools translate code automatically. Everything in this workbook transfers, but warp-size assumptions (32 vs 64) must be checked.

Exercises

Low precision and numerics

What you lose and gain with 16-bit and 8-bit numbers, and how training stays stable.

• After this module you can: describe FP32, TF32, FP16, BF16, FP8 and INT8 by their bits, range and precision, predict overflow and underflow, explain mixed-precision training (master weights, loss scaling), and design an accuracy-vs-speed experiment.

21.1 Floating-point formats

A float is sign × mantissa × 2^exponent. The exponent bits set the range (how big and small numbers can get); the mantissa bits set the precision (how many significant digits).

• FP16 vs BF16: same 16 bits, different trade-off. FP16 has more precision but a tiny range: values above 65,504 overflow to infinity, and gradients smaller than about 6×10⁻⁸ underflow to zero. BF16 keeps FP32's range with less precision, which is why it became the standard for training: no overflow or underflow, and you can convert from FP32 by just dropping bits.
• Machine epsilon (relative gap between neighbouring numbers): FP32 ≈ 1.2×10⁻⁷, FP16 ≈ 9.8×10⁻⁴, BF16 ≈ 7.8×10⁻³. In BF16, 256 + 1 = 256.
• Fewer bytes also halves or quarters memory traffic. For memory-bound kernels that's a speed-up even without Tensor Cores (Module 4).

21.2 Mixed-precision training

Training in pure 16-bit fails: weight updates (learning rate × gradient) are often smaller than the weight's precision, so w + Δw == w and the model stops learning. The standard recipe:

1. Keep FP32 master weights.
2. Each step, cast to BF16/FP16 for the forward and backward passes; GEMMs run on Tensor Cores with FP32 accumulation.
3. Apply the gradient update to the FP32 master copy.
4. With FP16 only, use loss scaling: multiply the loss by a factor S (e.g. 1024) so small gradients don't underflow, then divide the gradients by S before the update. Dynamic loss scaling lowers S when infinities appear. BF16 usually doesn't need it, thanks to its range.
5. Keep numerically sensitive operations in FP32: softmax, layernorm statistics, the loss, and large reductions.

llm.c's GPU training follows this pattern with BF16. That's what your project topic 3 measures.

21.3 Quantization for inference

For inference only, weights can be stored as INT8 (or 4-bit) with a scale: w ≈ scale × q with q ∈ [−127, 127], per tensor, per channel or per block of weights. Memory and bandwidth drop 4× relative to FP32, which directly speeds up memory-bound inference (the llama2.c matvec in Module 15 is exactly this case). Accuracy depends on how outliers are handled, so always measure perplexity or task accuracy, not just speed.

Exercises

CUDA libraries: cuBLAS, cuSPARSE, cuFFT, cuRAND, CUTLASS

Expert-tuned building blocks, and the handful of conventions you have to get right to use them.

• After this module you can: follow the library workflow, call cuBLAS GEMM correctly from row-major C code (the column-major trick), enable Tensor Cores through cuBLAS, generate random numbers with cuRAND in kernels, and say what CUTLASS, cuSPARSE, cuFFT and Thrust are for.

22.1 The common workflow (from the tutorial)

1. Create a handle (cublasCreate): an opaque object holding the library's context. Optionally set its stream (cublasSetStream).
2. Allocate device memory.
3. Convert inputs to a supported format (layout, sparse format).
4. Copy inputs to the device (cudaMemcpy or cublasSetMatrix/SetVector).
5. Configure (descriptors, transposes, algorithms).
6. Execute.
7. Retrieve results.
8. Convert back to your format.
9. Release resources, and reuse handles and buffers across calls, since creating them is expensive (the same lesson as Module 15).

22.2 cuBLAS and column-major order

cuBLAS follows Fortran BLAS: column-major storage, but zero-based pointers. Levels: 1 = vector ops (axpy, dot), 2 = matrix-vector (gemv), 3 = matrix-matrix (gemm). The leading dimension (lda) is the stride between columns, normally the number of rows.

The trick for row-major C code. A row-major M × N array, read as column-major, is its transpose (N × M). Since (A·B)ᵀ = Bᵀ·Aᵀ, ask cuBLAS for Cᵀ = Bᵀ·Aᵀ by simply swapping the operands. No copies, no transposes:

// row-major: C (M x P) = A (M x N) * B (N x P)
cublasHandle_t h; cublasCreate(&h);
const float alpha = 1.0f, beta = 0.0f;
cublasSgemm(h, CUBLAS_OP_N, CUBLAS_OP_N,
            P, M, N,              // dimensions of the column-major problem C^T = B^T * A^T
            &alpha,
            d_B, P,               // "first" matrix is B, leading dim P
            d_A, N,               // "second" matrix is A, leading dim N
            &beta,
            d_C, P);              // C, leading dim P
cublasDestroy(h);

Getting this right first time is rare. Always check against your CPU reference. The tutorial's cublasSgemv example works on genuinely column-major data, which is why it passes CUBLAS_OP_N with lda = M.

Tensor Cores through cuBLAS

• cublasGemmEx takes separate data types for inputs, output and compute: e.g. A, B in CUDA_R_16F or CUDA_R_16BF, C in CUDA_R_32F, compute type CUBLAS_COMPUTE_32F. cuBLAS picks Tensor Core kernels automatically when the types and alignment allow.
• For FP32 data on Ampere+, CUBLAS_COMPUTE_32F_FAST_TF32 (or cublasSetMathMode(h, CUBLAS_TF32_TENSOR_OP_MATH)) runs FP32 GEMMs on Tensor Cores in TF32.
• cuBLASLt is the more flexible interface (epilogues such as bias + GELU fused into the GEMM, FP8). llm.c uses it for its matmuls.

22.3 cuRAND: random numbers inside kernels (from the tutorial)

Random number generators are sequential, so each thread needs its own state, initialized so that threads produce different, independent streams:

#include <curand_kernel.h>
__global__ void init(curandState *st, unsigned long long seed) {
  int id = blockIdx.x * blockDim.x + threadIdx.x;
  curand_init(seed, id, 0, &st[id]);         // same seed, subsequence = thread id
}
__global__ void monteCarloPi(curandState *st, int trials, unsigned long long *hits) {
  int id = blockIdx.x * blockDim.x + threadIdx.x;
  curandState s = st[id];                    // copy state to registers (faster)
  unsigned int local = 0;
  for (int t = 0; t < trials; t++) {
    float x = curand_uniform(&s), y = curand_uniform(&s);
    local += (x * x + y * y <= 1.0f);
  }
  st[id] = s;                                // save state if you'll generate more later
  atomicAdd(hits, (unsigned long long)local); // one atomic per thread (Module 13)
}
// pi ≈ 4 * hits / (threads * trials)

• The tutorial seeds with seed = id. The documented best practice is the same seed with a different subsequence per thread (as above), which guarantees non-overlapping streams.
• curand_init is slow (it skips ahead in the sequence), so do it once in a separate kernel and reuse the states.
• For filling large arrays, the host API (curandCreateGenerator, curandGenerateUniform) is simpler and faster.

22.4 The rest of the ecosystem

CUTLASS exposes the same hierarchy you built by hand in Module 14, one level per layer: a threadblock tile staged through shared memory, divided into warp tiles, each computed with Tensor Core MMA instructions (the WMMA idea of Module 20), plus an epilogue that can fuse bias, activation or scaling before writing C. You choose tile sizes and data types as template parameters. Use it when cuBLAS doesn't offer the fused operation you need; for a plain GEMM, cuBLAS is usually as fast and far simpler. (The course reading used an early CUTLASS built on WMMA.)

Exercises

Triton, and programming several GPUs

A higher-level way to write GPU kernels, and how to spread work over more than one GPU.

• After this module you can: explain Triton's block-level programming model and write a simple Triton kernel, use cudaSetDevice and streams to drive several GPUs, enable and use peer-to-peer copies, and explain data-parallel training with all-reduce.

23.1 Triton: programming blocks, not threads

Triton (from OpenAI, used heavily inside PyTorch) is a Python-embedded language for GPU kernels. In CUDA you write what one thread does; in Triton you write what one program instance does to a whole block of data, as vector operations. The compiler then decides how to map that onto threads, warps, shared memory and Tensor Cores: the coalescing, tiling and memory staging you did by hand in Part 2.

import torch, triton, triton.language as tl

@triton.jit
def add_kernel(x_ptr, y_ptr, out_ptr, n, BLOCK: tl.constexpr):
    pid = tl.program_id(0)                       # like blockIdx.x
    offs = pid * BLOCK + tl.arange(0, BLOCK)     # a whole vector of indices
    mask = offs < n                              # replaces the "if (i < n)" bounds check
    x = tl.load(x_ptr + offs, mask=mask)
    y = tl.load(y_ptr + offs, mask=mask)
    tl.store(out_ptr + offs, x + y, mask=mask)

n = 1 << 24
x = torch.rand(n, device="cuda"); y = torch.rand(n, device="cuda"); out = torch.empty_like(x)
grid = (triton.cdiv(n, 1024),)                   # like (n + BLOCK - 1) / BLOCK
add_kernel[grid](x, y, out, n, BLOCK=1024)

• Mapping to CUDA: program_id ≈ blockIdx; tl.arange ≈ the threadIdx range; mask ≈ the bounds check; there's no explicit shared memory and no __syncthreads.
• For matrix multiply, tl.dot(a_block, b_block) on 2D blocks compiles to Tensor Core instructions, and @triton.autotune searches over block sizes automatically.
• Trade-off: much faster to write and often close to cuBLAS for GEMM-like and fused element-wise kernels (FlashAttention has a Triton version), but less control than CUDA, and debugging generated code is harder. It's the "productivity" learning outcome of the course.

23.2 Several GPUs in one node (from the slides)

Multi-GPU systems are connected two ways: GPUs in one node over PCIe or NVLink (much faster, hundreds of GB/s), and nodes over a network (InfiniBand). The lecture covers one node. The workflow:

1. Count devices: cudaGetDeviceCount(&ngpus).
2. For each device: cudaSetDevice(d), then create its stream and allocate its memory. Every CUDA call after cudaSetDevice targets that device.
3. Launch work on each GPU through its stream. Calls are asynchronous, so you can switch device while earlier GPUs are still running, and all GPUs work at once.
4. Wait on each stream.

for (int d = 0; d < ngpus; d++) {
  cudaSetDevice(d);
  int off = d * chunk;
  cudaMemcpyAsync(d_a[d], h_a + off, chunk * 4, cudaMemcpyHostToDevice, s[d]);   // pinned h_a
  vecAdd<<<(chunk + 255) / 256, 256, 0, s[d]>>>(d_a[d], d_b[d], d_c[d], chunk);
  cudaMemcpyAsync(h_c + off, d_c[d], chunk * 4, cudaMemcpyDeviceToHost, s[d]);
}
for (int d = 0; d < ngpus; d++) { cudaSetDevice(d); cudaStreamSynchronize(s[d]); }

Peer-to-peer (P2P)

int can; cudaDeviceCanAccessPeer(&can, 0, 1);          // can GPU 0 access GPU 1?
if (can) { cudaSetDevice(0); cudaDeviceEnablePeerAccess(1, 0); }   // one direction only
cudaMemcpyPeerAsync(d_dst_on1, 1, d_src_on0, 0, bytes, s0);   // GPU0 → GPU1 without the host

Access is enabled one direction at a time. With P2P, kernels on GPU 0 can even dereference pointers into GPU 1's memory. Without it, copies go through host memory and are much slower. NVIDIA's p2pBandwidthLatencyTest sample measures what your node supports.

23.3 Data-parallel training and NCCL

The most common multi-GPU pattern in ML is data parallelism: every GPU holds a full copy of the model and processes a different slice of the batch. After the backward pass, the gradients are averaged across GPUs with an all-reduce, so every replica applies the same update. NCCL implements all-reduce efficiently over NVLink, PCIe and networks (ring and tree algorithms). For a ring all-reduce over G GPUs, each GPU sends and receives about 2·(G−1)/G times the gradient size, nearly constant in G, which is why it scales. llm.c supports multi-GPU training this way (NCCL, with MPI to launch processes). Larger models also split the model itself (tensor or pipeline parallelism).

Exercises

Final drill: bonus oral exam and project defence

Answer each out loud in about a minute, with no notes. Brackets point to the module to reread. Do this weekly from December, together with Module 11.

Memory and kernels

Whole program

Project (llm.c)