GPU Fundamentals: Coalescing & Occupancy

A GPU kernel that does almost no arithmetic is limited by how fast it can move bytes. Two things decide that: whether a warp’s memory accesses can be coalesced into a few wide transactions, and whether enough warps are resident (occupancy) to keep memory busy. Every program below ran on a Tesla T4 with CUDA 13.0. Its peak DRAM bandwidth is 320 GB/s.

Coalesced versus strided reads

This kernel copies 16M floats, but each thread reads from stride floats past its neighbour. Every stride moves exactly the same 128 MB; only the address pattern changes.

// coalescing.cu
#include <cstdio>
#include <cuda_runtime.h>

__global__ void gather(const float* in, float* out, int n, int stride) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  int rows = n / stride;            // neighbouring threads read stride floats apart
  if (i < n) out[i] = in[(i % rows) * stride + i / rows];
}

int main() {
  const int n = 1 << 24;            // 16M floats: 64 MB per buffer
  float *in, *out; cudaMalloc(&in, n * 4); cudaMalloc(&out, n * 4);
  cudaEvent_t a, b; cudaEventCreate(&a); cudaEventCreate(&b);
  gather<<<n / 256, 256>>>(in, out, n, 1);   // warm-up
  for (int stride = 1; stride <= 32; stride *= 2) {
    cudaEventRecord(a);
    gather<<<n / 256, 256>>>(in, out, n, stride);
    cudaEventRecord(b); cudaEventSynchronize(b);
    float ms; cudaEventElapsedTime(&ms, a, b);
    printf("stride %2d: %6.2f ms, %5.1f GB/s\n", stride, ms, 2.0 * n * 4 / ms / 1e6);
  }
  return 0;
}

$ nvcc -o coalescing coalescing.cu && ./coalescing
stride  1:   0.55 ms, 244.1 GB/s
stride  2:   0.78 ms, 172.2 GB/s
stride  4:   1.25 ms, 107.8 GB/s
stride  8:   2.25 ms,  59.6 GB/s
stride 16:   4.44 ms,  30.2 GB/s
stride 32:   4.91 ms,  27.4 GB/s

With stride 1, a warp’s 32 reads cover 128 contiguous bytes: four 32-byte sectors, all of them used. That reaches 244 GB/s, 76% of the T4’s 320 GB/s peak.

Up to stride 16, each doubling of the stride roughly halves the bandwidth. The memory system fetches whole chunks, not single floats, so as the threads’ addresses spread out, more of each fetched chunk goes unused. From 16 to 32 the drop nearly stops. That fits DRAM being read in 64-byte pieces: 16 floats is 64 bytes, so from stride 16 on each thread already gets a piece to itself, and spreading further adds little waste. By stride 32 the copy runs about 9x slower, and nothing changed except the indexing.

The writes to out[i] stay contiguous at every stride. Only the reads change, and they decide the speed. Across three runs the first stride-1 measurement once came in at 198 GB/s; the rest held within a few percent. Expect your numbers to wobble, and compare ratios.

Rows, columns and a row-major matrix

C and CUDA store a 2D matrix row by row. Whether threadIdx.x (the index that changes fastest inside a warp) walks along a row or down a column decides whether the warp’s accesses coalesce.

// rows_cols.cu
#include <cstdio>
#include <cuda_runtime.h>

__global__ void copy(const float* M, float* out, int n, bool x_walks_cols) {
  int a = blockIdx.x * blockDim.x + threadIdx.x;   // changes fastest inside a warp
  int b = blockIdx.y * blockDim.y + threadIdx.y;
  int row = x_walks_cols ? b : a, col = x_walks_cols ? a : b;
  out[row * n + col] = M[row * n + col];           // row-major storage
}

float time_ms(const float* M, float* out, int n, bool x_walks_cols) {
  cudaEvent_t s, e; cudaEventCreate(&s); cudaEventCreate(&e);
  dim3 block(32, 8), grid(n / 32, n / 8);
  cudaEventRecord(s); copy<<<grid, block>>>(M, out, n, x_walks_cols);
  cudaEventRecord(e); cudaEventSynchronize(e);
  float ms; cudaEventElapsedTime(&ms, s, e); return ms;
}

int main() {
  const int n = 4096; float *M, *out;
  cudaMalloc(&M, n * n * 4); cudaMalloc(&out, n * n * 4);
  time_ms(M, out, n, true);                        // warm-up
  printf("threadIdx.x walks along a row:   %5.2f ms\n", time_ms(M, out, n, true));
  printf("threadIdx.x walks down a column: %5.2f ms\n", time_ms(M, out, n, false));
  return 0;
}

$ nvcc -o rows_cols rows_cols.cu && ./rows_cols
threadIdx.x walks along a row:    0.58 ms
threadIdx.x walks down a column:  2.05 ms

Both versions copy every element of a 4096 x 4096 matrix exactly once. When threadIdx.x picks the column, a warp touches 32 neighbouring floats in one row: one coalesced 128-byte run. When it picks the row, the warp’s 32 threads land in 32 different rows, 16 KB apart. That’s 32 separate transactions, and the copy is about 3.5x slower.

Over three runs the column walk took between 2.05 and 2.41 ms; the row walk held at 0.55 to 0.58 ms. The rule for 2D kernels: map threadIdx.x to the index that is contiguous in memory, which for row-major storage is the column.

Asking the runtime about occupancy

cudaOccupancyMaxActiveBlocksPerMultiprocessor reports how many blocks of a kernel fit on one SM at a given block size. From that you get warps per SM, and occupancy.

// occupancy.cu
#include <cstdio>
#include <cuda_runtime.h>

__global__ void work(const float* x, float* y, int n) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i < n) y[i] = x[i] * 2.0f;
}

int main() {
  cudaDeviceProp p; cudaGetDeviceProperties(&p, 0);
  int maxWarps = p.maxThreadsPerMultiProcessor / 32;
  int minGrid, best;
  cudaOccupancyMaxPotentialBlockSize(&minGrid, &best, work, 0, 0);
  printf("suggested block size: %d threads (grid of at least %d blocks)\n", best, minGrid);
  for (int threads = 32; threads <= 1024; threads *= 2) {
    int blocks;
    cudaOccupancyMaxActiveBlocksPerMultiprocessor(&blocks, work, threads, 0);
    int warps = blocks * threads / 32;
    printf("%4d threads/block: %2d blocks per SM, %2d of %d warps, %3d%% occupancy\n",
           threads, blocks, warps, maxWarps, 100 * warps / maxWarps);
  }
  return 0;
}

$ nvcc -o occupancy occupancy.cu && ./occupancy
suggested block size: 1024 threads (grid of at least 40 blocks)
  32 threads/block: 16 blocks per SM, 16 of 32 warps,  50% occupancy
  64 threads/block: 16 blocks per SM, 32 of 32 warps, 100% occupancy
 128 threads/block:  8 blocks per SM, 32 of 32 warps, 100% occupancy
 256 threads/block:  4 blocks per SM, 32 of 32 warps, 100% occupancy
 512 threads/block:  2 blocks per SM, 32 of 32 warps, 100% occupancy
1024 threads/block:  1 blocks per SM, 32 of 32 warps, 100% occupancy

A T4 SM holds at most 32 warps and at most 16 blocks. With 32-thread blocks, the 16-block limit is reached when only 16 warps are resident, so half the SM’s warp slots stay empty. From 64 threads up, the warp limit is reached first, and every size fills the SM.

The suggested size of 1024 is just the largest block that achieves full occupancy for this tiny kernel. “Grid of at least 40 blocks” is one block per SM on the T4’s 40 SMs. It’s a starting point, not a measurement.

Shared memory limits occupancy too

Registers and shared memory are split among the blocks on an SM. Ask for more per block and fewer blocks fit. Here the same kernel requests more and more dynamic shared memory at launch.

// smem_limit.cu
#include <cstdio>
#include <cuda_runtime.h>

__global__ void reverse(float* y) {
  extern __shared__ float buf[];                  // size set at launch
  buf[threadIdx.x] = threadIdx.x;
  __syncthreads();
  y[blockIdx.x * blockDim.x + threadIdx.x] = buf[blockDim.x - 1 - threadIdx.x];
}

int main() {
  cudaDeviceProp p; cudaGetDeviceProperties(&p, 0);
  int maxWarps = p.maxThreadsPerMultiProcessor / 32;
  for (int kb : {0, 8, 16, 32, 48}) {
    int blocks;
    cudaOccupancyMaxActiveBlocksPerMultiprocessor(&blocks, reverse, 256, kb * 1024);
    printf("256 threads + %2d KB shared per block: %d blocks per SM, %3d%% occupancy\n",
           kb, blocks, 100 * blocks * 8 / maxWarps);
  }
  return 0;
}

$ nvcc -o smem_limit smem_limit.cu && ./smem_limit
256 threads +  0 KB shared per block: 4 blocks per SM, 100% occupancy
256 threads +  8 KB shared per block: 4 blocks per SM, 100% occupancy
256 threads + 16 KB shared per block: 4 blocks per SM, 100% occupancy
256 threads + 32 KB shared per block: 2 blocks per SM,  50% occupancy
256 threads + 48 KB shared per block: 1 blocks per SM,  25% occupancy

A 256-thread block is 8 warps, which is the * 8 in the formula. Four blocks fill the SM’s 32 warp slots, and four blocks of 16 KB need exactly the T4’s 64 KB of shared memory per SM, so occupancy holds. At 32 KB per block only two fit, and at 48 KB only one. The kernel’s code didn’t change at all; a launch parameter halved, then quartered, its occupancy.

Does occupancy actually change the speed?

Occupancy is a means, not the goal. Time the simple work kernel from occupancy.cu at every block size and see where it matters.

// block_sweep.cu
#include <cstdio>
#include <cuda_runtime.h>

__global__ void work(const float* x, float* y, int n) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i < n) y[i] = x[i] * 2.0f;
}

int main() {
  const int n = 1 << 24;
  float *x, *y; cudaMalloc(&x, n * 4); cudaMalloc(&y, n * 4);
  cudaEvent_t a, b; cudaEventCreate(&a); cudaEventCreate(&b);
  work<<<n / 256, 256>>>(x, y, n);               // warm-up
  for (int threads = 32; threads <= 1024; threads *= 2) {
    cudaEventRecord(a);
    work<<<(n + threads - 1) / threads, threads>>>(x, y, n);
    cudaEventRecord(b); cudaEventSynchronize(b);
    float ms; cudaEventElapsedTime(&ms, a, b);
    printf("%4d threads/block: %5.2f ms, %5.1f GB/s\n", threads, ms, 2.0 * n * 4 / ms / 1e6);
  }
  return 0;
}

$ nvcc -o block_sweep block_sweep.cu && ./block_sweep
  32 threads/block:  0.73 ms, 183.8 GB/s
  64 threads/block:  0.55 ms, 242.1 GB/s
 128 threads/block:  0.55 ms, 243.3 GB/s
 256 threads/block:  0.55 ms, 243.8 GB/s
 512 threads/block:  0.55 ms, 244.0 GB/s
1024 threads/block:  0.55 ms, 243.6 GB/s

Only the 32-thread row is slow, and it’s the only one below 100% occupancy. With 16 warps per SM there isn’t enough in-flight work to hide memory latency. From 64 threads up, every size runs at the same 243 GB/s: once the warps are enough to keep memory busy, extra occupancy buys nothing.

The 32-thread row varied most across four runs (136 to 184 GB/s); the others stayed within 3%. That’s why “maximise occupancy” is the wrong rule. Get enough of it, then measure.

Common gotchas

• Trusting the indexing by eye. One swapped row/col turned a 0.58 ms copy into 2 ms. Time both layouts if you’re unsure.
• Chasing 100% occupancy. Above the point where latency is hidden, it doesn’t help a memory-bound kernel.
• Forgetting that shared memory and registers cost occupancy. A bigger tile can make each block faster and the SM emptier.
• Benchmarking once. The first measurement after launch is often the outlier.

Practice

1. Extend coalescing.cu to strides 64 and 128. Do you expect bandwidth to keep halving?

2. In smem_limit.cu, switch to 512-thread blocks (16 warps each). At what shared-memory size does occupancy first fall?

3. In rows_cols.cu, change the block shape to dim3 block(8, 32), grid(n / 8, n / 32). Predict both times.

4. Why does a 32-thread block reach only 50% occupancy on the T4, when 64 threads reach 100%? Confirm the limit with a device property.