GPU Architecture and CUDA Programming Model: Warps, Memory Hierarchy, and Thread Divergence

GPU Architecture and CUDA Programming Model: Warps, Memory Hierarchy, and Thread Divergence
GPU Architecture & HPC · Professor Pixel

GPU Architecture and CUDA Programming Model: Warps, Memory Hierarchy, and Thread Divergence

GPUs power every large language model, every modern game renderer, and every deep learning training run. But the difference between a naive GPU kernel and an optimized one can be 50× in performance — not because of algorithmic changes, but because of how you interact with the hardware. Understanding warps, shared memory bank conflicts, and memory coalescing is the difference between using a GPU and using it well.


1. GPU vs CPU: Two Fundamentally Different Design Philosophies

1.1 Latency-Optimized vs Throughput-Optimized

A modern CPU is a latency-optimized processor. It dedicates most of its die area to mechanisms that minimize the time to complete a single operation: large multi-level caches (often 30–50% of die area), sophisticated out-of-order execution engines, branch predictors, speculative execution units, and deep pipeline stages. A high-end CPU core can execute a single instruction in 1–4 nanoseconds, and it achieves this by executing instructions out of their program order, hiding cache miss latency speculatively, and prefetching data it thinks you'll need next. This makes CPUs excellent at latency-sensitive, sequential workloads — operating system scheduling, single-threaded business logic, random data access patterns.

A GPU is a throughput-optimized processor. Instead of minimizing latency per operation, it maximizes the total number of operations completed per second across thousands of concurrent threads. A modern NVIDIA H100 GPU has 16,896 CUDA cores (FP32), roughly 100× more than a high-end CPU. These cores are not individually sophisticated — they have tiny caches, simple in-order pipelines, and no branch prediction. They achieve high throughput not by making each operation fast, but by having so many concurrent operations in flight that memory latency is always hidden: while thousands of threads are stalled waiting for memory, other thousands are actively executing. This latency hiding through massive parallelism is the core GPU design principle.

The practical implication: CPUs win when you have complex control flow, random memory access patterns, or sequential dependencies between operations. GPUs win when you have large data-parallel workloads where the same operation is applied independently to millions of elements — exactly the structure of neural network forward and backward passes, graphics shading, physics simulations, and cryptographic operations.

Developer Pitfall — GPU Overhead for Small Workloads:

Launching a CUDA kernel has a fixed overhead of 5–15 microseconds (kernel launch latency, scheduling, driver synchronization). For a workload that takes 2 microseconds to compute, this overhead makes GPU execution 5–7× slower than CPU. GPUs are not universally faster — they require large enough work to amortize the launch overhead. A common mistake is moving small matrix operations (e.g., 16×16 matrix multiply) to the GPU to "use the GPU." The correct threshold: your kernel should execute for at least hundreds of microseconds to benefit from GPU execution. Batch small operations together or use libraries like cuBLAS that amortize launch overhead across many queued operations.


2. Streaming Multiprocessors: The GPU's Core Unit of Execution

2.1 SM Architecture: What's Inside a Streaming Multiprocessor

A GPU is organized as a collection of Streaming Multiprocessors (SMs). An NVIDIA A100 has 108 SMs; an H100 (SXM5) has 132 SMs. Each SM is an independent execution engine containing: a set of CUDA cores (FP32/INT32 execution units), Tensor Cores (specialized matrix multiplication units), a warp scheduler (typically 4 per SM), a register file (shared among all active threads on the SM), L1 cache / shared memory (a fast SRAM bank, configurable split), load/store units, and special function units (SFUs) for transcendental operations (sin, cos, exp).

The register file is a critical resource: on an A100 SM, each SM has 65,536 32-bit registers. These registers are divided among all threads actively resident on the SM. If a kernel uses 64 registers per thread and runs 1024 threads per SM, it uses exactly 65,536 registers — the maximum. Adding one more active thread would require more registers than available, forcing the compiler to spill registers to slower local memory (global memory). Register pressure is a direct constraint on how many threads can simultaneously reside on an SM, and therefore on the GPU's ability to hide latency.

2.2 The CUDA Thread Hierarchy

CUDA organizes threads into a three-level hierarchy that maps directly to the GPU hardware:

flowchart TB subgraph Grid["CUDA Grid (entire kernel launch)"] subgraph B1["Thread Block 0\n(up to 1024 threads)"] W1["Warp 0\n(threads 0-31)"] W2["Warp 1\n(threads 32-63)"] W3["Warp N\n(threads ...)"] end subgraph B2["Thread Block 1\n(up to 1024 threads)"] W4["Warp 0\n(threads 0-31)"] W5["Warp 1\n(threads 32-63)"] end subgraph BN["Thread Block M\n(...)"] W6["..."] end end subgraph HW["Hardware Mapping"] SM0["SM 0\n(executes Block 0 + Block 1)"] SM1["SM 1\n(executes Block 2 + Block 3)"] SMN["SM N\n(executes Block M...)"] end B1 -->|"Scheduled to"| SM0 B2 -->|"Scheduled to"| SM0 BN -->|"Scheduled to"| SMN style Grid fill:#fef3c7,stroke:#f59e0b style HW fill:#dbeafe,stroke:#3b82f6

Diagram 1: CUDA Thread Hierarchy. A Grid is subdivided into Thread Blocks; each Block is subdivided into Warps of 32 threads. Thread Blocks are scheduled to SMs by the GPU thread block scheduler. Warps are the actual unit of execution within an SM.

Developer Pitfall — Thread Block Size Not a Multiple of 32:

Warps always contain exactly 32 threads. If your thread block size is not a multiple of 32, the last warp is padded with inactive threads (predicated off). A thread block of 100 threads creates 4 warps: warps 0-2 have 32 active threads each, warp 3 has only 4 active threads and 28 inactive. The inactive threads still occupy warp scheduler slots and consume register file space — they participate in execution but with their results discarded. This wastes 28/32 = 87.5% of the last warp's compute capacity. Always size thread blocks as multiples of 32. Common choices: 128, 256, or 512 threads per block, selected to maximize occupancy (covered in Section 10).


3. Warps and SIMT: Single Instruction, Multiple Threads

3.1 What a Warp Actually Is

A warp is a group of 32 threads that execute the same instruction simultaneously in a lockstep fashion on the SM's CUDA cores. This is the GPU's fundamental execution unit — not a single thread, not a thread block, but a warp of 32. The warp scheduler selects a ready warp (one where threads are not stalled waiting for memory or synchronization) and issues one instruction from that warp to the execution units. All 32 threads execute that instruction in the same cycle, operating on their own private register data. This model is called SIMT — Single Instruction, Multiple Threads.

Why exactly 32? This is an architectural decision by NVIDIA, driven by the tradeoff between scheduling flexibility (more warps = more opportunities to hide latency) and hardware cost (warp context storage, scheduler complexity). The number 32 has remained constant across every NVIDIA GPU architecture since Tesla (2006). AMD's equivalent unit on RDNA/CDNA GPUs is called a wavefront and contains 64 threads (though RDNA 2+ supports a 32-thread subgroup mode). The differences in warp/wavefront size affect how you tune kernel code for each vendor's hardware.

An SM on the A100 can have up to 64 warps resident simultaneously (2048 threads total). The warp scheduler runs every cycle and can select among all resident warps. A warp that issues a global memory load will stall for 200–800 cycles waiting for the data. During those 800 cycles, the scheduler switches to other ready warps — hiding the memory latency completely, as long as there are enough other warps to fill those cycles. This is how GPUs "tolerate" high memory latency: not by making memory faster, but by having enough concurrent work to keep the execution units busy while memory requests are outstanding.

// CUDA kernel: each thread processes one element
__global__ void addVectors(float* A, float* B, float* C, int N) {
int tid = blockIdx.x * blockDim.x + threadIdx.x;
if (tid < N) {
C[tid] = A[tid] + B[tid]; // All 32 warp threads execute this simultaneously
}
}
 
// Launch: 256 threads per block, enough blocks for N=1,048,576 elements
int blocks = (N + 256 - 1) / 256; // ceil(N / 256)
addVectors<<<blocks, 256>>>(d_A, d_B, d_C, N);
// Each SM receives thread blocks and executes their warps in parallel
// A100 (108 SMs): processes all 4096 blocks across 108 SMs simultaneously

4. Thread Divergence: The Silent Performance Killer

4.1 What Divergence Is and Why It Hurts

Since all 32 threads in a warp execute the same instruction simultaneously, they must all follow the same code path. When a branch instruction produces different outcomes for different threads in the same warp — some threads take the if branch and others take the else branch — the warp diverges. The GPU hardware handles this by serializing the two paths: it executes the if branch with the "true" threads active and the "false" threads predicated off (their results are discarded), then executes the else branch with the "false" threads active and the "true" threads predicated off. The warp reconverges after both branches complete.

The performance cost is proportional to the number of divergent paths. A warp that splits into two paths takes 2× as long as a non-divergent warp. A warp with 32 fully divergent threads (each taking a unique path) takes up to 32× as long. In practice, most divergence is binary (if/else), and the actual slowdown depends on the relative work in each branch. But in the worst case — image processing where some pixels are in a bright region and others are in a dark region, each requiring different processing — divergence can cut effective GPU throughput by 50%.

// DIVERGENT kernel — different threads take different branches
__global__ void processData(float* data, int N) {
int tid = blockIdx.x * blockDim.x + threadIdx.x;
if (tid < N) {
if (data[tid] > 0.5f) { // ← DIVERGENCE POINT
data[tid] = sqrtf(data[tid]); // threads 0,3,7,12... execute
} else {
data[tid] = data[tid] * data[tid]; // threads 1,2,4,5... execute
} // ← RECONVERGENCE
}
}
 
// Warp execution timeline (32 threads, assume ~50% take each branch):
// Cycle 1-5: threads[0,3,7,12,...] execute sqrtf; threads[1,2,4,5,...] IDLE
// Cycle 6-8: threads[1,2,4,5,...] execute multiply; threads[0,3,7,12,...] IDLE
// Effective throughput: ~50% of non-divergent warp
 
// AVOIDING DIVERGENCE — sort/partition data so warps are coherent:
// If all threads in a warp have data[tid] > 0.5 OR all have data[tid] <= 0.5,
// the branch is taken by all 32 threads → no divergence.

4.2 Strategies to Reduce Divergence

Data layout reorganization is the most effective strategy: sort or partition your data so that threads within the same warp tend to take the same branch. In a BVH tree traversal (ray tracing), rays from nearby pixels tend to diverge when they hit different geometry — grouping rays by direction before dispatching them to warps significantly reduces divergence. In machine learning inference with variable-length sequences, padding all sequences in a batch to the same length (or grouping sequences of similar length together) eliminates the divergence that would otherwise occur when some threads have finished processing their sequence while others haven't.

Predication instead of branching is another approach: when the branch body is very short (1–3 instructions), using conditional moves (the compiler may do this automatically) avoids the branch entirely. Both paths execute, but the result of the "wrong" path is discarded via a predicate mask. This eliminates the serialization overhead of a divergent branch when the branch body itself is shorter than the divergence penalty. For complex branches with many instructions in each arm, predication is worse than divergence because it executes both paths unconditionally — choose based on the branch body length.

Developer Pitfall — Divergence on threadIdx Within Warps Across Blocks:

A subtle but common divergence mistake: branching on threadIdx.x % 2 == 0 (odd vs even threads) within a block. Since threads 0–31 form warp 0 and threads 32–63 form warp 1, the odd/even split occurs WITHIN each warp — half the threads in every warp take one branch and half take the other. This maximally diverges every warp in the kernel. Compare to branching on blockIdx.x % 2 == 0 (odd vs even blocks): entire blocks take the same branch, so entire warps take the same branch — zero divergence. Always think about divergence at the warp granularity (groups of 32 consecutive threads), not at the thread granularity.


5. The GPU Memory Hierarchy: Where Your Data Lives Determines Your Speed

5.1 Six Memory Spaces, Six Different Performance Profiles

GPU performance is fundamentally about memory. The arithmetic peak of a GPU (FLOPs) is almost always higher than what most kernels achieve — the bottleneck is moving data to and from the compute units. Understanding where each memory space lives, how fast it is, and how to use it correctly is the single most impactful GPU optimization skill. The following chart shows the approximate bandwidth of each memory space on an NVIDIA A100 (a commonly used reference GPU for production ML):

Chart 1: Approximate read bandwidth for each CUDA memory space on NVIDIA A100. Registers are effectively infinite bandwidth (direct SM pipeline access); each subsequent level involves increasing latency and decreasing bandwidth. The ~80× difference between registers and global memory is why kernel optimization is almost entirely about data placement.

5.2 CUDA Memory Types Reference

Memory Type Location Scope Latency Size (A100) Best For
Registers On-SM Thread-private 1 cycle 64K/SM Local variables, loop counters
Shared Memory On-SM (SRAM) Block-shared ~5 cycles 192KB/SM Inter-thread data sharing, tiling
L1 Cache On-SM (unified with shared) Automatic ~30 cycles Configurable Read-only data, scatter-gather
L2 Cache On-chip Device-wide ~200 cycles 40MB Frequently reused global data
Global Memory HBM2e (off-chip) Device-wide ~500–800 cycles 80GB Main data storage (arrays, tensors)
Constant Memory HBM + on-SM cache Device-wide, read-only 1 cycle (if cached) 64KB Kernel parameters, lookup tables

Developer Pitfall — Register Spilling to Local Memory:

When a kernel uses more registers than the SM's register file can hold for the number of resident threads, the compiler automatically spills excess register values to local memory — which is physically the same as global memory (HBM), with 500–800 cycle latency. Register spilling silently destroys performance. Check for spilling with nvcc --ptxas-options=-v, which prints "spills" in the compilation output. Reduce register usage by: splitting large kernels into multiple smaller kernels, reducing the number of simultaneously live variables, or using __launch_bounds__(maxThreadsPerBlock, minBlocksPerSM) to tell the compiler to optimize register allocation for a specific occupancy target.


6. Shared Memory Bank Conflicts: The Hidden Bottleneck

6.1 How Shared Memory Is Organized Into Banks

Shared memory is organized into 32 banks (matching the warp size of 32). Each bank services one 32-bit (4-byte) read or write per cycle. Successive 32-bit words are interleaved across banks: word 0 goes to bank 0, word 1 to bank 1, ..., word 31 to bank 31, word 32 back to bank 0, and so on. When 32 threads in a warp each access a different bank, all 32 accesses complete in one cycle (perfect parallelism). When two or more threads in a warp access the same bank (but different addresses), a bank conflict occurs — the accesses are serialized, reducing effective bandwidth.

A 32-way bank conflict (all 32 threads access the same bank) reduces shared memory throughput by 32×, making it slower than a poorly-aligned global memory access. The exception: if all threads in a warp access the same address in the same bank, the hardware broadcasts the value to all threads in a single cycle — no conflict. Bank conflicts only occur when threads access different addresses that map to the same bank.

// 32 banks × 4 bytes = 128-byte interleaving pattern
// Address → Bank: bank = (address / 4) % 32
 
// CONFLICT-FREE: stride-1 access (consecutive addresses)
__shared__ float tile[32][32];
float val = tile[threadIdx.y][threadIdx.x];
// thread 0 → tile[y][0] → bank 0
// thread 1 → tile[y][1] → bank 1 ... thread 31 → bank 31. No conflict!
 
// 32-WAY BANK CONFLICT: stride-32 access
float val = tile[threadIdx.x][0]; // ALL threads access column 0
// thread 0 → tile[0][0] → address 0 → bank 0
// thread 1 → tile[1][0] → address 128 → bank 0 (128/4=32, 32%32=0)
// thread 31 → tile[31][0] → address 3968 → bank 0. 32-way conflict!
 
// FIX: pad the shared memory array by 1 column
__shared__ float tile[32][33]; // +1 column shifts bank assignments
float val = tile[threadIdx.x][0];
// Now tile[0][0]=bank0, tile[1][0]=bank(33%32)=bank1 ... no conflict!

Developer Pitfall — Transposing a Matrix via Shared Memory and Creating Bank Conflicts:

Matrix transposition is the classic bank conflict example. When reading a matrix row-by-row into shared memory (coalesced global read) and then writing it to output column-by-column (transposed write), the write step accesses shared memory with a stride-32 pattern — a 32-way bank conflict. The fix is always to pad shared memory: declare __shared__ float tile[BLOCK][BLOCK+1] instead of __shared__ float tile[BLOCK][BLOCK]. The extra padding column shifts consecutive rows to different banks, eliminating conflicts with zero algorithmic overhead. This one-character fix can improve transposition kernel performance by 10–15×. Always check the access stride when designing shared memory layouts.


7. Memory Coalescing: How to Get Full Global Memory Bandwidth

7.1 What Coalescing Means and Why It Matters

Global memory (HBM) is accessed in cache lines of 128 bytes. When a warp issues a load instruction, the memory controller services all 32 threads' addresses in as few transactions as possible. If the 32 threads access 32 consecutive floats (128 bytes = one cache line), all accesses are served in a single transaction — this is coalesced access. If the 32 threads access 32 scattered addresses (e.g., stride-2 or random), each address may fall on a different cache line, requiring up to 32 separate transactions — a 32× bandwidth reduction from the hardware's peak.

Coalescing has a dramatic real-world impact. A kernel with perfectly coalesced global memory access saturates the theoretical HBM bandwidth (2 TB/s on A100). The same kernel with stride-2 access (every other float) achieves only ~1 TB/s. With fully random access (e.g., indirect array indexing like B[perm[i]]), effective bandwidth can drop to 10–50 GB/s — below even PCIe bandwidth — making the GPU no faster than the CPU for that operation.

// COALESCED: thread i accesses element i (consecutive, 128-byte aligned)
float val = A[blockDim.x * blockIdx.x + threadIdx.x];
// Warp accesses: A[0],A[1],...,A[31] → 1 memory transaction
 
// UN-COALESCED: threads access with stride (array-of-structures pattern)
struct Particle { float x, y, z, w; }; // 16 bytes per struct
Particle* particles = ...;
float x = particles[tid].x; // stride-4 across structs → un-coalesced
// Thread 0 reads bytes 0-3, Thread 1 reads bytes 16-19... 32 transactions!
 
// FIX: Structure-of-Arrays (SoA) instead of Array-of-Structures (AoS)
float* x_arr; float* y_arr; float* z_arr; float* w_arr;
float x = x_arr[tid]; // Coalesced! x[0],x[1],...,x[31] in one transaction

Developer Pitfall — Column-Major vs Row-Major Access in 2D Arrays:

C arrays are row-major: element [i][j] is at address base + i*cols + j. If thread tid accesses column j using row index tid (i.e., A[tid][j]), consecutive threads access consecutive rows but the same column — stride-cols access, un-coalesced. If thread tid accesses A[i][tid] (same row, different columns), consecutive threads access consecutive addresses — coalesced. For matrix operations, assign the innermost loop variable (the one that varies fastest between threads) to correspond to column indices (j), and ensure j = threadIdx.x for row-major matrices. This is why cuBLAS requires column-major matrices — Fortran convention — it naturally coalesces when thread IDs map to row indices, which are consecutive in column-major storage.


8. Step-by-Step: Optimizing Matrix Multiplication from Naïve to Tiled

8.1 The Naïve Kernel: Each Thread Reads Entire Rows and Columns

Matrix multiplication C = A × B (M×K × K×N → M×N) is the canonical CUDA optimization problem because it starts with terrible memory behavior and can be systematically improved to near-peak bandwidth. In the naïve version, each thread computes one element of C by iterating over K elements from a row of A and a column of B:

// NAÏVE matrix multiply — each thread computes one C[row][col]
__global__ void matmulNaive(float* A, float* B, float* C, int M, int N, int K) {
int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;
if (row < M && col < N) {
float sum = 0.0f;
for (int k = 0; k < K; k++) {
sum += A[row * K + k] * B[k * N + col]; // ← PROBLEM HERE
}
C[row * N + col] = sum;
}
}
 
// Problem analysis (warp processes 32 threads, varying col):
// A[row * K + k]: all threads read same row, same k → BROADCAST ✓ (ok)
// B[k * N + col]: all threads read same k, different col → COALESCED ✓ (ok)
// But: for K=1024, each thread makes 2048 global memory reads.
// A 32×32 block of threads reads A's rows K times each → massive redundancy!
// 32 threads in a warp all read the same A[row][k] element K times.

8.2 The Tiled Kernel: Using Shared Memory to Reduce Global Reads

The fundamental inefficiency in the naïve kernel: 32 threads in a warp all read the same element of A (since they share the same row) — they're reading the same data from global memory 32 times. The fix is to use shared memory to cache tiles of A and B that are reused by multiple threads in the same block, loading each element exactly once from global memory and then reusing it from the fast shared memory:

#define TILE_SIZE 32
 
__global__ void matmulTiled(float* A, float* B, float* C, int M, int N, int K) {
__shared__ float tileA[TILE_SIZE][TILE_SIZE]; // 32×32×4 = 4KB
__shared__ float tileB[TILE_SIZE][TILE_SIZE]; // 4KB (total 8KB shared mem)
 
int row = blockIdx.y * TILE_SIZE + threadIdx.y;
int col = blockIdx.x * TILE_SIZE + threadIdx.x;
float sum = 0.0f;
 
for (int t = 0; t < (K + TILE_SIZE - 1) / TILE_SIZE; t++) {
// Phase 1: Cooperatively load tiles from global → shared memory
tileA[threadIdx.y][threadIdx.x] = (row < M && t*TILE_SIZE+threadIdx.x < K)
? A[row * K + t * TILE_SIZE + threadIdx.x] : 0.0f;
tileB[threadIdx.y][threadIdx.x] = (col < N && t*TILE_SIZE+threadIdx.y < K)
? B[(t * TILE_SIZE + threadIdx.y) * N + col] : 0.0f;
 
__syncthreads(); // Wait for all threads to finish loading tiles
 
// Phase 2: Compute partial dot product from shared memory tiles
for (int k = 0; k < TILE_SIZE; k++) {
sum += tileA[threadIdx.y][k] * tileB[k][threadIdx.x];
} // All reads from shared memory — ~5 cycle latency
 
__syncthreads(); // Wait before loading next tile (don't overwrite live data)
}
if (row < M && col < N) C[row * N + col] = sum;
}
 
// Memory traffic reduction:
// Naïve: 2*M*N*K global reads (each element of A/B read M or N times)
// Tiled: 2*M*N*K/TILE_SIZE global reads (each tile element read once → reused TILE_SIZE times)
// With TILE_SIZE=32: 32× reduction in global memory bandwidth. Real-world speedup: 10–20×

Developer Pitfall — Forgetting __syncthreads() Between Tile Phases:

The two __syncthreads() barriers in the tiled kernel are not optional. The first ensures all threads have finished loading their tile elements into shared memory before any thread starts the dot product phase (otherwise fast threads start computing with partially-loaded tiles). The second ensures all threads have finished consuming the current tile before any thread overwrites it with the next tile. Missing either barrier produces silently incorrect results — the kernel will execute, produce wrong output, and there is no runtime error. This is one of the most common CUDA correctness bugs. Use tools like cuda-memcheck --tool racecheck to detect shared memory races during development.


9. Occupancy and the Roofline Model: Diagnosing What Limits Your Kernel

9.1 GPU Occupancy: Filling the SM to Hide Latency

Occupancy is the ratio of active warps on an SM to the maximum number of warps the SM supports. An A100 SM supports up to 64 concurrent warps (2048 threads). If your kernel launches 1024 threads per SM (32 warps), occupancy is 32/64 = 50%. Higher occupancy generally gives the warp scheduler more warps to switch to when one stalls on memory — more opportunities to hide latency. But occupancy is not the only metric that matters, and blindly maximizing it can hurt performance (more active threads → fewer registers per thread → register spilling → slower).

Occupancy is constrained by three resources: (1) the number of thread blocks and threads per block, (2) register usage per thread (higher registers → fewer active threads), and (3) shared memory per block (more shared memory → fewer concurrent blocks → fewer warps). The CUDA Occupancy Calculator (available in NVIDIA NSight Systems and as an Excel spreadsheet from NVIDIA) computes the binding constraint given your kernel's resource usage. The target is not 100% occupancy — it's enough occupancy to hide the memory latency of your specific kernel, which depends on the memory access pattern and arithmetic intensity.

9.2 The Roofline Model: Compute-Bound vs Memory-Bound

The Roofline Model is a simple but powerful framework for characterizing GPU kernel performance. It plots kernel performance (GFLOPs/s) as a function of arithmetic intensity (FLOPs per byte of data moved from global memory). The model has two limits: the memory bandwidth roof (horizontal line at peak bandwidth × arithmetic intensity) and the compute roof (horizontal line at peak FLOPs). A kernel's performance is bounded by whichever roof it hits first.

For an NVIDIA A100: peak FP32 compute = 19.5 TFLOPS, peak HBM bandwidth = 2 TB/s. The ridge point (where compute and bandwidth limits intersect) is at arithmetic intensity = 19,500 / 2,000 = 9.75 FLOPs/byte. A kernel with arithmetic intensity below 9.75 is memory-bound — adding more compute units won't help; you need to reduce data movement. Above 9.75 FLOPs/byte is compute-bound — adding more bandwidth won't help; you need more FP32 units or better instruction-level parallelism. The tiled matrix multiply above achieves ~15–20 FLOPs/byte (compute-bound for large matrices), while the naïve version achieves ~2 FLOPs/byte (memory-bound). The roofline predicts the tiled kernel should be 7–10× faster before you even benchmark — and it is.

Developer Pitfall — Using NSight Compute's "Compute Throughput" Without Memory Throughput:

NVIDIA NSight Compute reports "Compute (SM) Throughput" and "Memory (DRAM) Throughput" as percentages of peak. A kernel with 90% compute throughput and 5% memory throughput is compute-bound — excellent! But a kernel with 50% compute throughput and 50% memory throughput is at the ridge point and may or may not have room to improve. A common mistake is seeing 90% memory throughput and 10% compute throughput and concluding "we're efficient." You're not — you're memory-bound, and 90% of your SM's compute capacity is idle waiting for memory. The fix in this case is to increase arithmetic intensity (do more math per byte loaded), not to optimize the memory access pattern (which is already saturating bandwidth).


10. Advanced: CUDA Streams and Concurrent Execution

10.1 What CUDA Streams Are and Why They Matter

By default, all CUDA operations (kernel launches, memory copies) execute sequentially in the default stream. CUDA Streams are sequences of GPU operations that execute in order within the stream, but different streams can execute concurrently — kernel execution in stream A can overlap with a memory copy in stream B. This is particularly valuable for hiding the latency of CPU-to-GPU data transfers (which are typically the production bottleneck, not the kernel itself).

The classic pattern is double buffering with streams: while the GPU is processing batch N from GPU memory, the CPU is simultaneously copying batch N+1 from host memory to GPU memory. With a single stream (no overlap), the pipeline is: copy → compute → copy → compute. With double buffering: copy₀ + compute₋₁ overlap, then copy₁ + compute₀ overlap — effectively hiding the copy latency behind compute. For large models where data loading is the bottleneck (common in training data pipelines), this overlap can reduce wall-clock time by 20–40%.

// Double buffering with two CUDA streams
cudaStream_t streams[2];
cudaStreamCreate(&streams[0]);
cudaStreamCreate(&streams[1]);
 
// Pre-allocate pinned host memory (required for async copies)
float *h_A0, *h_A1, *d_A0, *d_A1;
cudaMallocHost(&h_A0, BATCH_SIZE * sizeof(float)); // pinned memory
cudaMallocHost(&h_A1, BATCH_SIZE * sizeof(float));
cudaMalloc(&d_A0, BATCH_SIZE * sizeof(float));
cudaMalloc(&d_A1, BATCH_SIZE * sizeof(float));
 
for (int i = 0; i < numBatches; i++) {
int s = i % 2; // Alternate between streams 0 and 1
float* h_buf = (s == 0) ? h_A0 : h_A1;
float* d_buf = (s == 0) ? d_A0 : d_A1;
 
loadBatchToPinnedMemory(h_buf, i); // CPU fills pinned buffer
cudaMemcpyAsync(d_buf, h_buf, BATCH_SIZE * sizeof(float),
cudaMemcpyHostToDevice, streams[s]); // Async H2D copy
myKernel<<<grid, block, 0, streams[s]>>>(d_buf); // Kernel in same stream
// While this stream's kernel runs, other stream's copy overlaps!
}
cudaDeviceSynchronize(); // Wait for all streams to finish

Developer Pitfall — Pageable Host Memory Prevents Async Copies:

cudaMemcpyAsync only performs a truly asynchronous copy (overlapping with kernel execution) when the source/destination host memory is pinned (page-locked, allocated with cudaMallocHost() or cudaHostAlloc()). If you call cudaMemcpyAsync with regular malloc()-allocated host memory, CUDA silently falls back to synchronous behavior — the copy blocks the CPU until complete, and no overlap occurs with other streams. Always use cudaMallocHost() for host buffers involved in async transfers. The tradeoff: pinned memory cannot be paged out by the OS, reducing available virtual memory for other processes. Only pin what you actually need for transfers.


11. Frequently Asked Questions

Q1: When should I use cuBLAS instead of writing my own matrix multiplication kernel?

Almost always — use cuBLAS. NVIDIA's cuBLAS library implements GEMM (General Matrix Multiplication) with architecture-specific optimizations: register blocking, vectorized loads (LDG.128), asynchronous memory copies (LDGSTS on Ampere+), and Tensor Core utilization (for FP16/BF16/TF32). A well-tuned cuBLAS SGEMM achieves 90–95% of the A100's theoretical FP32 peak — far beyond what a hand-written tiled kernel achieves (typically 50–70% of peak). Writing a custom matmul kernel only makes sense when your matrix has a special sparsity structure, fused operations (bias add + activation in one pass), or non-standard precision that cuBLAS doesn't support. For everything else, cuBLAS is the professional standard, and using it correctly (choosing the right algorithm via cublasGemmEx with CUBLAS_GEMM_DEFAULT_TENSOR_OP) gives you immediate access to years of NVIDIA engineering effort.

Q2: What are Tensor Cores and how do they differ from CUDA cores?

CUDA cores are general-purpose execution units that perform one floating-point multiply-add (FMA) per cycle per core. Tensor Cores are specialized matrix units that perform a 4×4×4 (or 16×16×16 in later generations) matrix multiply-accumulate in a single operation — equivalent to 256 FMAs. On an A100, Tensor Cores deliver 312 TFLOPS (FP16) vs 19.5 TFLOPS for CUDA cores — a 16× increase. The catch: Tensor Cores require inputs in specific reduced-precision formats (FP16, BF16, TF32, INT8) and matrix dimensions that are multiples of 16 (for Ampere). Deep learning training exploits this by using mixed-precision (FP16 forward/backward, FP32 accumulation) to access Tensor Core throughput while maintaining FP32 numerical stability for the optimizer. Access Tensor Cores via cuBLAS with the TENSOR_OP flag, or through CUTLASS (NVIDIA's open-source template library for composable GEMM kernels).

Q3: What is the difference between CUDA and OpenCL for GPU programming?

CUDA is NVIDIA-specific, tightly integrated with the NVIDIA toolchain (nvcc compiler, NVIDIA NSight profiler, cuBLAS/cuDNN libraries), and consistently offers the best performance on NVIDIA hardware because NVIDIA optimizes both the runtime and the hardware together. OpenCL is vendor-neutral and runs on GPUs from NVIDIA, AMD, Intel, and Apple Silicon, plus FPGAs and CPUs. The tradeoff: CUDA provides a better developer experience (more comprehensive documentation, more tutorials, better debugging tools), while OpenCL provides portability at the cost of higher boilerplate code and lower peak performance on any given hardware. For deep learning workloads, CUDA dominates because PyTorch, TensorFlow, and JAX all use CUDA/cuDNN as their primary backend. AMD offers HIP (Heterogeneous Integrated Platform) as a CUDA-compatible API that allows porting CUDA code to AMD GPUs with minimal changes — a better portability option than OpenCL for code already written in CUDA.

Q4: How does CUDA handle error checking, and why is it often skipped in examples?

Every CUDA API call returns a cudaError_t status, and kernel launches set an error state accessible via cudaGetLastError(). In production code, all calls must be checked. Most tutorial code omits error checks to keep examples short, which is a bad habit — CUDA errors propagate silently: a failed memory allocation returns a null pointer that causes an out-of-bounds write that produces garbage output with no runtime error on the GPU. The minimum production wrapper: #define CUDA_CHECK(call) { cudaError_t e = (call); if(e != cudaSuccess) { fprintf(stderr, "CUDA error %s:%d: %s\n", __FILE__, __LINE__, cudaGetErrorString(e)); exit(1); } }. Kernel errors are asynchronous — use CUDA_CHECK(cudaDeviceSynchronize()) after kernels during development to surface them immediately rather than discovering them when the next synchronous call fails.

Q5: What is flash attention and why does it matter for GPU memory efficiency?

Flash Attention (Dao et al., 2022) is an I/O-aware implementation of the attention mechanism in Transformers that dramatically reduces global memory reads/writes. Standard attention for sequence length N requires O(N²) global memory I/O — it materializes the full N×N attention matrix in HBM. For N=8192, this is 8192² × 4 bytes = 256GB of HBM reads/writes — completely memory-bound and the main bottleneck in long-context transformer inference. Flash Attention reformulates the softmax computation to use online normalization, allowing the algorithm to compute attention in tiles that fit in shared memory — never materializing the full attention matrix in HBM. The result: O(N) HBM I/O instead of O(N²), making attention compute-bound instead of memory-bound. Flash Attention 2 and 3 further optimize for Tensor Core utilization and asynchronous memory pipelines, enabling 2–3× faster attention vs standard implementations on A100.

Q6: How does thread block size affect performance beyond occupancy?

Thread block size affects performance through three mechanisms beyond raw occupancy: (1) Warp-level efficiency — block sizes that are exact multiples of 32 prevent inactive partial-warp threads; (2) Shared memory usage — larger blocks can cooperatively share more data in shared memory, reducing global reads per thread (the tiled matmul uses a 32×32 = 1024 thread block specifically to load 32×32 tiles); and (3) L1 cache hit rate — smaller blocks with more concurrent blocks on the SM compete for L1 cache lines, potentially reducing hit rate for streaming access patterns. The empirical rule: use 256 threads per block as a default, measure with NSight, then experiment with 128 or 512 if metrics show the kernel is register-pressure-limited or shared-memory-limited respectively. Never pick block size based on intuition alone — benchmark with NVIDIA NSight Compute's occupancy analysis tool.

Q7: What are CUDA graphs and when should I use them?

CUDA Graphs (introduced in CUDA 10) capture a sequence of operations (kernel launches, memory copies) as a single executable graph that the GPU driver can execute with minimal CPU overhead. Normally, each kernel launch and memcpy requires a CPU-side API call (~5–15µs overhead each) before the GPU starts. For iterative workloads (training loops, iterative solvers) that repeat the same sequence of kernels, each iteration pays this overhead. A CUDA graph records the sequence once and then replays it each iteration with a single CPU call — eliminating the per-launch overhead entirely. This provides 2–20× CPU overhead reduction for workloads dominated by many small kernels. The limitation: the graph structure (kernel parameters, memory addresses, block dimensions) is fixed at capture time. Any dynamic change (different batch size, different input shape) requires re-capturing the graph. PyTorch 2.0's torch.compile() uses CUDA graphs internally for static-shape workloads.

Q8: What is unified memory and when should I use it?

CUDA Unified Memory (cudaMallocManaged()) creates a single pointer accessible from both CPU and GPU, with the driver migrating pages automatically between host and device memory as needed. This dramatically simplifies memory management — no explicit cudaMemcpy calls, no separate host and device pointers. However, the migration mechanism has significant overhead: the first GPU access to a page triggers a page fault, the page is migrated from CPU RAM to GPU HBM (hundreds of microseconds per page), and the same happens in reverse for CPU access. For workloads with good data locality (all GPU accesses, then all CPU accesses), unified memory performs comparably to explicit copies after the initial migration. For workloads with frequent CPU↔GPU data sharing, the page fault overhead makes unified memory slower. Use unified memory for development (easier debugging, no memory management bugs) and switch to explicit copies for production performance-critical paths.

Q9: How does GPU memory differ between training and inference for deep learning?

Training requires significantly more GPU memory than inference for the same model size. Inference: you store model weights (FP16: 2 bytes/parameter) and the activation tensors for the current batch — for a 7B parameter model, weights need ~14GB. Training additionally stores: (1) optimizer state — AdamW stores first and second moment estimates per parameter (8 bytes/parameter for FP32 moments = 56GB for 7B params), (2) gradients (same size as weights, 14GB), and (3) activations for all layers that must be retained for backpropagation (can exceed the model size for long sequences). Gradient checkpointing (recomputing activations during backward pass instead of storing them) trades compute for memory — typically 30–40% memory reduction at 30% compute overhead. Mixed-precision training (BF16 compute, FP32 optimizer state) with DeepSpeed ZeRO-3 sharding across multiple GPUs is the standard approach for training large models on limited GPU budgets.

Q10: What profiling tools should I use to identify GPU performance bottlenecks?

NVIDIA provides two primary profiling tools: NSight Systems and NSight Compute. NSight Systems is a system-level profiler that shows CPU↔GPU interaction, kernel execution timelines, memory transfer overlaps, and stream concurrency — essential for identifying whether your CPU-to-GPU pipeline has idle gaps (copy-compute overlap issues, kernel launch delays, synchronization stalls). NSight Compute is a kernel-level profiler that drills into individual kernel performance: SM occupancy, warp stall reasons, memory throughput (global/shared), cache hit rates, instruction-level statistics, and roofline model positioning. The workflow: start with NSight Systems to identify which kernels are the hotspots and whether the pipeline structure is correct, then use NSight Compute on the hotspot kernels to diagnose the specific bottleneck (memory-bound, compute-bound, latency-bound) and measure the impact of optimizations. PyTorch also provides torch.profiler which integrates with NSight Systems and adds Python-level stack traces to GPU timeline events.


Written by Professor Pixel · CodingPancake · GPU Architecture & HPC Series

Post a Comment

Previous Post Next Post