CUDA Memory Hierarchy & Optimization

Most CUDA kernels spend more time waiting for memory than doing arithmetic. The levels, from fastest to slowest:

This lesson measures each one instead of taking it on trust. Every program ran on a Tesla T4 with CUDA 13.0, whose memory peaks at 320 GB/s. Timed kernels are averaged over 10 launches after a warm-up launch, and the numbers moved between runs, so the text gives the range across all the runs.

Coalesced vs Scattered Reads

Both kernels copy every element of a 4096 × 4096 matrix once. The only difference is the order in which a warp’s 32 threads read.

// coalesce.cu - copy a 4096 x 4096 matrix reading along rows, then down columns
#include <cstdio>
#include <cuda_runtime.h>

__global__ void init(float* A) { int t = blockIdx.x * blockDim.x + threadIdx.x; A[t] = (float)t; }

__global__ void copy_rows(const float* A, float* out, int N) {
  int t = blockIdx.x * blockDim.x + threadIdx.x;
  out[t] = A[t];                                  // a warp reads 32 neighbouring floats
}

__global__ void copy_cols(const float* A, float* out, int N) {
  int t = blockIdx.x * blockDim.x + threadIdx.x;
  out[t] = A[(t % N) * N + t / N];                // a warp reads 32 floats, each N apart
}

template <typename K> float avg_ms(K kernel, const float* A, float* out, int N) {
  cudaEvent_t a, b;
  cudaEventCreate(&a); cudaEventCreate(&b);
  kernel<<<N * N / 256, 256>>>(A, out, N);        // warm-up
  cudaEventRecord(a);
  for (int r = 0; r < 10; r++) kernel<<<N * N / 256, 256>>>(A, out, N);
  cudaEventRecord(b); cudaEventSynchronize(b);
  float ms; cudaEventElapsedTime(&ms, a, b);
  return ms / 10;
}

int main() {
  const int N = 4096;
  const size_t bytes = (size_t)N * N * sizeof(float);   // 64 MiB per matrix
  float *A, *out, h;
  cudaMalloc(&A, bytes); cudaMalloc(&out, bytes);
  init<<<N * N / 256, 256>>>(A);

  float r = avg_ms(copy_rows, A, out, N);
  float c = avg_ms(copy_cols, A, out, N);
  cudaMemcpy(&h, out + 1, sizeof(float), cudaMemcpyDeviceToHost);
  double gb = 2.0 * bytes / 1e9;                         // read A + write out
  printf("copy_rows: %.3f ms  %3.0f GB/s\n", r, gb / (r / 1e3));
  printf("copy_cols: %.3f ms  %3.0f GB/s  (%.1fx slower)\n", c, gb / (c / 1e3), c / r);
  printf("copy_cols out[1] = %.0f (A[1][0])\n", h);
  cudaFree(A); cudaFree(out);
  return 0;
}

$ nvcc -O2 -arch=sm_75 -o coalesce coalesce.cu && ./coalesce
copy_rows: 0.533 ms  252 GB/s
copy_cols: 1.322 ms  102 GB/s  (2.5x slower)
copy_cols out[1] = 4096 (A[1][0])

Across four runs, copy_rows took 0.533 to 0.547 ms (245 to 252 GB/s, up to 79% of the T4’s peak) and copy_cols took 1.13 to 1.32 ms (102 to 118 GB/s). The column version was 2.1 to 2.5 times slower, for exactly the same bytes and arithmetic.

In copy_rows, a warp’s 32 reads are 128 contiguous bytes, which the hardware serves with four 32-byte sectors. In copy_cols, the 32 threads read one float from each of 32 rows, so each float sits in a different sector and most of every sector is waste. Counting reads and writes, that’s 36 sectors per warp against 8, so a naive estimate is about 4.5 times slower. It’s 2.1 to 2.5 times slower, most likely because the caches keep some of those sectors long enough for later threads to use the rest. out[1] = 4096 confirms the column version really did read A[1][0].

A Stencil: Shared Memory vs the Cache

A 7-point moving average reads every input seven times, so it looks like a classic case for shared memory: load each input once into a tile, then read the neighbours from on-chip memory.

// stencil.cu - 7-point moving average: global reads vs a shared-memory tile
#include <cstdio>
#include <cuda_runtime.h>
#define R 3                                       // radius: 3 neighbours each side
#define BLOCK 256

__global__ void avg_global(const float* in, float* out, int N) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i < R || i >= N - R) return;
  float s = 0.0f;
  for (int k = -R; k <= R; k++) s += in[i + k];   // 7 global loads per thread
  out[i] = s / (2 * R + 1);
}

__global__ void avg_shared(const float* in, float* out, int N) {
  __shared__ float tile[BLOCK + 2 * R];           // the block's inputs plus a halo each side
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  int t = threadIdx.x + R;
  if (i < N) tile[t] = in[i];                     // one global load per thread
  if (threadIdx.x < R) {                          // the first R threads also load the halos
    if (i >= R) tile[t - R] = in[i - R];
    if (i + BLOCK < N) tile[t + BLOCK] = in[i + BLOCK];
  }
  __syncthreads();                                // wait until the whole tile is loaded
  if (i < R || i >= N - R) return;
  float s = 0.0f;
  for (int k = -R; k <= R; k++) s += tile[t + k]; // 7 shared-memory loads
  out[i] = s / (2 * R + 1);
}

template <typename K> float run(K kernel, const float* in, float* out, int N, int* wrong) {
  cudaMemset(out, 0, N * sizeof(float));
  kernel<<<N / BLOCK, BLOCK>>>(in, out, N);       // warm-up, also the result we check
  float* h = new float[N];
  cudaMemcpy(h, out, N * sizeof(float), cudaMemcpyDeviceToHost);
  *wrong = 0;
  for (int i = R; i < N - R; i++) if (h[i] != 3.0f) (*wrong)++;
  delete[] h;
  cudaEvent_t a, b; cudaEventCreate(&a); cudaEventCreate(&b);
  cudaEventRecord(a);
  for (int r = 0; r < 10; r++) kernel<<<N / BLOCK, BLOCK>>>(in, out, N);
  cudaEventRecord(b); cudaEventSynchronize(b);
  float ms; cudaEventElapsedTime(&ms, a, b);
  return ms / 10;
}

int main() {
  const int N = 1 << 24;
  float* h = new float[N];
  for (int i = 0; i < N; i++) h[i] = (float)(i % 7);   // any 7 neighbours average to 3
  float *in, *out;
  cudaMalloc(&in, N * sizeof(float)); cudaMalloc(&out, N * sizeof(float));
  cudaMemcpy(in, h, N * sizeof(float), cudaMemcpyHostToDevice);

  int wg, ws;
  float g = run(avg_global, in, out, N, &wg);
  float s = run(avg_shared, in, out, N, &ws);
  double gb = 2.0 * N * sizeof(float) / 1e9;
  printf("avg_global: %.3f ms  %3.0f GB/s  wrong=%d\n", g, gb / (g / 1e3), wg);
  printf("avg_shared: %.3f ms  %3.0f GB/s  wrong=%d\n", s, gb / (s / 1e3), ws);
  delete[] h; cudaFree(in); cudaFree(out);
  return 0;
}

$ nvcc -O2 -arch=sm_75 -o stencil stencil.cu && ./stencil
avg_global: 0.547 ms  245 GB/s  wrong=0
avg_shared: 0.588 ms  228 GB/s  wrong=0

Both versions are correct, since every window of 7 consecutive values of i % 7 averages to exactly 3. But the shared-memory version was slower in every run: 0.59 to 0.92 ms against 0.55 to 0.71 ms for the plain version. (Both slowed down over the session, probably as the GPU’s clock dropped.) The tile saves global loads, but the T4 had already been saving them.

Counting the Loads with ncu

Nsight Compute can count what each kernel actually asked for and what really came from DRAM. -k picks the kernel and -c 1 profiles one launch.

$ ncu -k avg_global -c 1 --metrics dram__bytes_read.sum,l1tex__t_requests_pipe_lsu_mem_global_op_ld.sum ./stencil | grep -E '^  avg_|dram__|l1tex__' && ncu -k avg_shared -c 1 --metrics dram__bytes_read.sum,l1tex__t_requests_pipe_lsu_mem_global_op_ld.sum ./stencil | grep -E '^  avg_|dram__|l1tex__'
  avg_global(const float *, float *, int) (65536, 1, 1)x(256, 1, 1), Context 1, Stream 7, Device 0, CC 7.5
    dram__bytes_read.sum                                  Mbyte        78.01
    l1tex__t_requests_pipe_lsu_mem_global_op_ld.sum                3,670,016
  avg_shared(const float *, float *, int) (65536, 1, 1)x(256, 1, 1), Context 1, Stream 7, Device 0, CC 7.5
    dram__bytes_read.sum                                  Mbyte        77.08
    l1tex__t_requests_pipe_lsu_mem_global_op_ld.sum                  655,358

The request counts are exact. avg_global made 65,536 blocks × 8 warps × 7 = 3,670,016 warp-wide load requests. avg_shared made 524,288 for the tile plus 2 halo loads per block, less the missing left halo of the first block and right halo of the last: 655,358, 5.6 times fewer.

The DRAM line is what matters, though. The input array is 2^24 floats, 67.1 MB, and both kernels read about 77 to 78 MB from DRAM (four runs agreed within 0.15 MB). Without a cache the plain version would have read close to 7 × 67 MB. The L1 and L2 caches had already turned seven reads into roughly one, so the tile removed requests that were cheap anyway. Its extra halo branches and __syncthreads() are most likely why it ended up slower; this lesson didn’t time them separately.

Shared memory wins when each value is reused many times in a pattern the cache can’t hold. Tiled matrix multiplication is the standard example: each element of a tile is reused 16 or 32 times. The shared memory and warps lesson measures it on the same GPU.

What the Compiler Allocated

-Xptxas -v makes the compiler’s assembler report each kernel’s resources. -c stops after compiling, so nothing runs. Without the grep, each kernel also gets a compile-time line that changes from run to run.

$ nvcc -O2 -arch=sm_75 -Xptxas -v -c stencil.cu -o /dev/null 2>&1 | grep -v 'Compile time'
ptxas info    : 0 bytes gmem
ptxas info    : Compiling entry function '_Z10avg_sharedPKfPfi' for 'sm_75'
ptxas info    : Function properties for _Z10avg_sharedPKfPfi
    0 bytes stack frame, 0 bytes spill stores, 0 bytes spill loads
ptxas info    : Used 18 registers, used 1 barriers, 1048 bytes smem, 372 bytes cmem[0]
ptxas info    : Compiling entry function '_Z10avg_globalPKfPfi' for 'sm_75'
ptxas info    : Function properties for _Z10avg_globalPKfPfi
    0 bytes stack frame, 0 bytes spill stores, 0 bytes spill loads
ptxas info    : Used 18 registers, used 0 barriers, 372 bytes cmem[0]

Kernel names appear in their C++ mangled form: _Z10avg_sharedPKfPfi is avg_shared(const float*, float*, int). Each kernel uses 18 registers per thread with no spills. avg_shared adds 1048 bytes of shared memory, which is (256 + 2 × 3) floats × 4 bytes, and 1 barrier, its __syncthreads(). cmem[0] is constant memory holding the kernel’s parameters plus constants the driver reserves. Registers and shared memory per block are what limit how many blocks fit on an SM at once; the next lesson measures that.

Common Gotchas

• Measure before tiling. The cache already handles simple neighbour reuse. Profile first, and use shared memory where the counts show DRAM traffic far above the data size.
• A missing __syncthreads(). Without the barrier, some threads read tile slots before they’re written. Practice 1 shows the result.
• Strided access hidden in an index. t % N and t / N look harmless, but they decide whether a warp’s reads are neighbours.
• Shared memory has costs too. Every byte per block reduces how many blocks an SM can hold, and banked access can serialise. The cuda-cpp lesson covers bank conflicts.

Practice

1. Delete the __syncthreads() from avg_shared and run the program three times. What happens to wrong?

$ nvcc -O2 -arch=sm_75 -o p1_nosync p1_nosync.cu && ./p1_nosync
avg_global: 0.588 ms  228 GB/s  wrong=0
avg_shared: 0.693 ms  194 GB/s  wrong=228782

2. Make the stencil radius 8 (17 reads per output, with inputs i % 17, which average to 8). Does shared memory win now that the reuse is larger?

$ nvcc -O2 -arch=sm_75 -o p3_radius8 p3_radius8.cu && ./p3_radius8
avg_global: 0.858 ms  156 GB/s  wrong=0
avg_shared: 1.053 ms  127 GB/s  wrong=0

3. Read the same resource numbers at run time with cudaFuncGetAttributes, and check that they match ptxas.

// p4_attrs.cu - read a kernel's registers and shared memory at run time
#include <cstdio>
#include <cuda_runtime.h>
#define R 3
#define BLOCK 256

__global__ void avg_global(const float* in, float* out, int N) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i < R || i >= N - R) return;
  float s = 0.0f;
  for (int k = -R; k <= R; k++) s += in[i + k];
  out[i] = s / (2 * R + 1);
}

__global__ void avg_shared(const float* in, float* out, int N) {
  __shared__ float tile[BLOCK + 2 * R];
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  int t = threadIdx.x + R;
  if (i < N) tile[t] = in[i];
  if (threadIdx.x < R) {
    if (i >= R) tile[t - R] = in[i - R];
    if (i + BLOCK < N) tile[t + BLOCK] = in[i + BLOCK];
  }
  __syncthreads();
  if (i < R || i >= N - R) return;
  float s = 0.0f;
  for (int k = -R; k <= R; k++) s += tile[t + k];
  out[i] = s / (2 * R + 1);
}

void show(const char* name, const void* kernel) {
  cudaFuncAttributes a;
  cudaFuncGetAttributes(&a, kernel);
  printf("%s: %d registers, %zu bytes shared, max %d threads per block\n",
         name, a.numRegs, a.sharedSizeBytes, a.maxThreadsPerBlock);
}

int main() {
  show("avg_global", (const void*)avg_global);
  show("avg_shared", (const void*)avg_shared);
  return 0;
}

$ nvcc -O2 -arch=sm_75 -o p4_attrs p4_attrs.cu && ./p4_attrs
avg_global: 18 registers, 0 bytes shared, max 1024 threads per block
avg_shared: 18 registers, 1048 bytes shared, max 1024 threads per block

4. Before profiling the radius-8 stencil, predict its two load-request counts. Then run ncu on it and look at the DRAM bytes too.

$ ncu -k avg_global -c 1 --metrics dram__bytes_read.sum,l1tex__t_requests_pipe_lsu_mem_global_op_ld.sum ./p3_radius8 | grep -E '^  avg_|dram__|l1tex__' && ncu -k avg_shared -c 1 --metrics dram__bytes_read.sum,l1tex__t_requests_pipe_lsu_mem_global_op_ld.sum ./p3_radius8 | grep -E '^  avg_|dram__|l1tex__'
  avg_global(const float *, float *, int) (65536, 1, 1)x(256, 1, 1), Context 1, Stream 7, Device 0, CC 7.5
    dram__bytes_read.sum                                  Mbyte        88.49
    l1tex__t_requests_pipe_lsu_mem_global_op_ld.sum                8,912,896
  avg_shared(const float *, float *, int) (65536, 1, 1)x(256, 1, 1), Context 1, Stream 7, Device 0, CC 7.5
    dram__bytes_read.sum                                  Mbyte        81.76
    l1tex__t_requests_pipe_lsu_mem_global_op_ld.sum                  655,358