GPU Memory Hierarchy and Occupancy: Why a Tiled Kernel Is Just Cache Blocking by Hand

A GPU is the memory-hierarchy lesson of the matrix-multiplication post with part of the cache taken out of the hardware's hands and placed in yours. The CPU hides latency with large caches and speculation; the GPU has caches too, but relies mainly on keeping thousands of threads resident and switching between them while loads are in flight. The consequences are that the layout of memory accesses across a group of threads matters as much as the layout within one thread, that there is a small, fast memory you must fill explicitly, and that how many threads can be resident at once (occupancy) is a resource budget you can compute by hand.

WARNING

This laptop has no NVIDIA GPU and no CUDA toolkit, so the kernels in this post were written and reasoned about but not compiled or run here. There are no timings. Every quantity I use is either an architectural constant from the CUDA C Programming Guide (warp size, bank count, sector size) or a worked example whose assumptions are stated next to it. Treat the code as a careful draft and run it through nvcc and compute-sanitizer before trusting it.

The execution model

A kernel launch names a grid of blocks, each block a 1-, 2- or 3-D array of threads. The programmer sees threads and blocks; the hardware sees warps of 32 threads that execute one instruction at a time in lockstep (SIMT), and streaming multiprocessors (SMs) onto which whole blocks are placed.

Three rules fall out of this picture and drive everything else in the post:

• Threads in a warp share an instruction stream. A branch that splits a warp is executed both ways with masks, so divergence inside a warp costs real time; divergence between warps costs nothing.
• Threads can synchronise and share fast memory only within a block, because only a block is guaranteed to be co-resident on one SM.
• An SM holds several blocks at once and interleaves their warps. A warp that is waiting on a load is simply not scheduled; a different warp runs. Latency is hidden mainly by having enough warps, not by caches.

Five memories, and who manages them

The important asymmetry is between the two fast memories of an SM. L1, like the CPU's caches, is filled by hardware on a miss and evicts by its own policy; shared memory is filled by instructions you write. That is a loss of convenience and a gain of control: a tile you stage in shared memory stays there until you overwrite it, with no replacement policy to second-guess.

Coalescing: the warp is the unit of memory access

When a warp executes a load, the hardware does not issue 32 independent requests. It looks at the 32 addresses, groups them into 32-byte sectors (four sectors make a 128-byte L2 line) and fetches each distinct sector once. If thread $t$ reads a 4-byte element at base + 4t and base is 32-byte aligned, the 32 addresses fall into exactly four sectors and every fetched byte is used; shift base by one element and the same 128 bytes straddle five sectors, 160 bytes moved for 128 used. If thread $t$ reads at base + 32t, every address lands in its own sector and the warp pulls 1024 bytes to use 128 of them.

A useful way to score an access pattern is its efficiency,

which is $1$ for the aligned contiguous case, $128/160=0.8$ for the misaligned one and $128/1024=1/8$ for a stride of 32 bytes. For a row-major matrix this means threadIdx.x must index the column (the fastest-varying dimension), exactly the stride-1 rule that made the i-k-j loop order win on the CPU. A warp walking down a column of a row-major matrix is the GPU form of the stride-$N$ walk down B in the textbook i-j-k loop.

Shared memory and bank conflicts

Shared memory is split into 32 banks, each 4 bytes wide, and word $w$ lives in bank $wmod32$. One bank can serve one address per cycle (a broadcast of the same address to several threads is free). A warp reading tile[ty][tx] along tx touches 32 consecutive words, one per bank, and finishes in one pass. A warp reading tile[tx][ty] along tx with a 32-wide tile touches words $32t$, which all map to bank 0: a 32-way conflict, serialised into 32 passes.

The classic fix is to declare the tile as tile[32][33]. The row stride becomes 33 words, $33tmod32=t$, and the column read is conflict-free again. It costs one wasted word per row and is the same move as padding per-thread counters onto separate cache lines in the false-sharing post: change the layout so that the hardware's hashing function spreads the accesses.

A tiled matrix multiplication

The naive kernel assigns one thread per element of $C$ and reads $A$ and $B$ straight from global memory. Within a warp the A[row*N + k] read is a broadcast (same address) and the B[k*N + col] read is coalesced, so the access pattern is fine, but every element of $A$ and $B$ is requested $N$ times across the grid. Broadcasts and the L1 and L2 caches absorb some of those requests; nothing in the kernel makes the reuse explicit, and nothing guarantees it.

• Naive
• Tiled

• __global__ void matmul_naive(const float *A, const float *B, float *C, int N)
  {
      int row = blockIdx.y * blockDim.y + threadIdx.y;
      int col = blockIdx.x * blockDim.x + threadIdx.x;
      if (row >= N || col >= N) return;
      float acc = 0.f;
      for (int k = 0; k < N; k++) acc += A[row * N + k] * B[k * N + col];
      C[row * N + col] = acc;
  }
  

  Two 4-byte loads per fused multiply-add: $2$ flops per $8$ bytes requested from the memory system, an arithmetic intensity of $0.25$ flop/byte before any cache help.

• #define TILE 32
  __global__ void matmul_tiled(const float *A, const float *B, float *C, int N)
  {
      __shared__ float As[TILE][TILE];
      __shared__ float Bs[TILE][TILE + 1];      // +1 pad keeps column reads conflict-free
  
      int tx = threadIdx.x, ty = threadIdx.y;
      int row = blockIdx.y * TILE + ty;
      int col = blockIdx.x * TILE + tx;
      float acc = 0.f;
  
      for (int k0 = 0; k0 < N; k0 += TILE) {
          // cooperative, coalesced loads: tx is the fastest-varying index for both
          As[ty][tx] = (row < N && k0 + tx < N) ? A[row * N + k0 + tx] : 0.f;
          Bs[ty][tx] = (k0 + ty < N && col < N) ? B[(k0 + ty) * N + col] : 0.f;
          __syncthreads();                       // the whole tile is in before anyone reads it
  
  #pragma unroll
          for (int k = 0; k < TILE; k++)
              acc += As[ty][k] * Bs[k][tx];      // As: broadcast across the warp; Bs: stride 1
          __syncthreads();                       // everyone is done before the tile is overwritten
      }
      if (row < N && col < N) C[row * N + col] = acc;
  }
  

  Each thread now loads two floats per tile step and performs TILE multiply-adds out of shared memory: $2\cdot 32=64$ flops per $8$ bytes, an intensity of $8$ flop/byte against global memory, 32 times the naive kernel.

In roofline terms (the first post of this series), and counting the loads the source issues, tiling moves the kernel TILE times further to the right:

Whether 8 flop/byte is past the ridge depends on the GPU: the ridge is that device's own ratio of peak arithmetic to memory bandwidth, and on a part whose single-precision peak in flop/s is more than eight times its bandwidth in bytes per second, the tiled kernel is still on the bandwidth side. What tiling guarantees is the factor of 32, not the crossing.

The two __syncthreads() calls are load-bearing. The first guarantees that the tile is complete before any thread reads it; the second guarantees that no thread starts loading the next tile while a slower warp is still reading the current one. Dropping the second is the classic race that produces results that are usually right. The pad on Bs is not needed for the read pattern in this kernel (Bs[k][tx] is already a row read) but it makes the kernel safe to transpose later without rediscovering bank conflicts.

dim3 block(TILE, TILE);                               // 1024 threads = 32 warps
dim3 grid((N + TILE - 1) / TILE, (N + TILE - 1) / TILE);

cudaEvent_t t0, t1;
cudaEventCreate(&t0); cudaEventCreate(&t1);
cudaEventRecord(t0);
matmul_tiled<<<grid, block>>>(dA, dB, dC, N);
cudaEventRecord(t1);
cudaEventSynchronize(t1);
float ms; cudaEventElapsedTime(&ms, t0, t1);

int blocksPerSM;
cudaOccupancyMaxActiveBlocksPerMultiprocessor(&blocksPerSM, matmul_tiled, TILE * TILE, 0);
cudaDeviceProp prop; cudaGetDeviceProperties(&prop, 0);
printf("%d blocks/SM, %d of %d warps resident\n", blocksPerSM,
       blocksPerSM * TILE * TILE / 32, prop.maxThreadsPerMultiProcessor / 32);

A block reduction

Summing $n$ numbers is the other kernel every GPU course writes, because it is where divergence, bank conflicts and synchronisation all show up in a few lines. The structure that avoids the known traps:

1. Each thread starts by summing many elements in a grid-stride loop, so the number of blocks is decoupled from $n$ and the per-thread work amortises launch and reduction overhead.
2. The block then reduces in shared memory with sequential addressing: at step $s$, thread $t<s$ adds sdata[t + s] into sdata[t]. Active threads are contiguous, so whole warps retire together and no warp diverges until fewer than 32 threads remain, and the addresses are consecutive, so there are no bank conflicts.
3. The last 32 values are reduced inside one warp with __shfl_down_sync, which exchanges registers directly without shared memory or a barrier.

__inline__ __device__ float warp_reduce(float v)
{
    for (int offset = 16; offset > 0; offset >>= 1)
        v += __shfl_down_sync(0xffffffffu, v, offset);
    return v;
}

__global__ void reduce_sum(const float *in, float *partial, int n)
{
    extern __shared__ float sdata[];
    unsigned tid = threadIdx.x;
    unsigned i = blockIdx.x * (blockDim.x * 2) + tid;
    unsigned stride = gridDim.x * blockDim.x * 2;

    float v = 0.f;
    for (; i < (unsigned)n; i += stride) {          // grid-stride: many elements per thread
        v += in[i];
        if (i + blockDim.x < (unsigned)n) v += in[i + blockDim.x];
    }
    sdata[tid] = v;
    __syncthreads();

    for (unsigned s = blockDim.x / 2; s > 32; s >>= 1) {   // sequential addressing
        if (tid < s) sdata[tid] += sdata[tid + s];
        __syncthreads();
    }
    if (tid < 32) {                                        // final warp in registers
        v = sdata[tid] + sdata[tid + 32];
        v = warp_reduce(v);
        if (tid == 0) partial[blockIdx.x] = v;
    }
}

The kernel produces one partial sum per block; a second launch with a single block (or a final atomicAdd) finishes. It assumes a one-dimensional block whose size is a power of two from 64 up to the device's limit, and blockDim.x * sizeof(float) bytes of dynamic shared memory, as in reduce_sum<<<blocks, 256, 256 * sizeof(float)>>>(in, partial, n); the finishing launch needs the same kind of block, however few partial sums remain. Outside that range the tree breaks: with 32 threads sdata[tid + 32] reads past the end of the array, and with 96 threads some values are added twice, so 96 ones sum to 112. The interleaved-addressing variant that textbooks show first, if (tid % (2*s) == 0), is correct but leaves most lanes of every warp idle from the first step. Its usual repair, index = 2*s*tid, keeps the active threads contiguous but strides through shared memory, which brings bank conflicts: two-way at the first step and worse after. Sequential addressing avoids both. Mark Harris's reduction slides, linked below, walk through seven versions of this kernel; the code above is close to the last one without its template unrolling, and with a warp shuffle, which the slides predate, in place of their unrolled last warp.

Occupancy as a budget

Occupancy is the fraction of an SM's warp slots that are filled: resident warps divided by the maximum. Blocks are placed whole, so the budget is counted in blocks, and the binding resource is whichever runs out first:

Here $R_{\text{SM}}$ is the register file size in 32-bit registers, $R_{\text{thread}}$ the registers the compiler assigned per thread, $T_{\text{block}}$ and $w_{\text{block}}$ the threads and warps per block, $S$ shared memory, and $W_{\text{max}}$, $B_{\text{max}}$ the architectural limits on resident warps and blocks. (Real hardware allocates registers and shared memory in granules, so the exact numbers come from the occupancy API or Nsight Compute; the formula gives the shape.)

A worked example with stated, illustrative assumptions: an SM with $R_{\text{SM}}=\mathrm{65,536}$ registers, $W_{\text{max}}=64$, $B_{\text{max}}=32$ and $S_{\text{SM}}=100$ KiB, running the tiled kernel above with 1024-thread blocks ($w_{\text{block}}=32$) that use $S_{\text{block}}=32\cdot 32\cdot 4+32\cdot 33\cdot 4=8320$ bytes of shared memory.

• If the compiler needs $R_{\text{thread}}=40$, one block needs $1024\cdot 40=\mathrm{40,960}$ registers and only $\lfloor 65536/40960\rfloor =1$ fits: one block, 32 warps, 50 % occupancy, with 24,576 registers left unused because half a block cannot run.
• If $R_{\text{thread}}\le 32$, registers allow $\lfloor 65536/(1024\cdot 32)\rfloor =2$ blocks, shared memory $\lfloor 102400/8320\rfloor =12$, block slots 32 and warp slots $64/32=2$: two blocks, 64 warps, 100 %.

Eight registers per thread decide whether the SM is half empty. This is why --ptxas-options=-v, __launch_bounds__ and the -maxrregcount flag exist, and why a block of 1024 threads is often a worse choice than 256: smaller blocks pack into the leftover space more flexibly.

TIP

High occupancy is a means, not the goal. It buys latency hiding, but a kernel can also hide latency with more independent work per thread (instruction-level parallelism) and fewer resident warps, and the register-heavy version is often faster because the extra registers hold reused tiles. Vasily Volkov's talk linked below makes that case with measurements; the matching intuition from the CPU side is that keeping data in registers beats re-reading it from any cache.

Mapping back to the CPU posts

The vocabulary is different and the numbers are different, but the questions a profiler answers are the same: how many bytes moved for the useful work done, and how much of the machine was waiting.

Further reading

• CUDA C Programming Guide: the execution model, the memory hierarchy and the compute-capability tables that give the exact limits for each architecture.
• CUDA C Best Practices Guide: coalescing, shared-memory and occupancy guidance with the official terminology.
• Mark Harris, Optimizing Parallel Reduction in CUDA: seven versions of the reduction kernel and why each is faster than the last.
• How to Access Global Memory Efficiently in CUDA C/C++ Kernels and Using Shared Memory in CUDA C/C++, two short posts with measured stride and bank-conflict experiments.
• CUDA Pro Tip: Occupancy API Simplifies Launch Configuration, the API used in the host code above.
• Vasily Volkov, Better Performance at Lower Occupancy, GTC 2010.