Shared Memory & Warp-Level Programming in CUDA C++

Global memory is large and slow; shared memory is a small on-chip scratchpad that every thread in a block can read and write. This lesson uses it to speed up a matrix multiply, shows the two classic ways to get it wrong, then moves values between threads with warp shuffles, which skip memory altogether. Every program ran on a Tesla T4 with CUDA 13.0.

The naive multiply

Each thread computes one element of C = A x B by walking a row of A and a column of B. The inputs are small integers so the result is exact and can be checked against the CPU.

// naive.cu - matrix multiply, every operand read from global memory (N = 1024)
#include <cstdio>
#include <cuda_runtime.h>
#define T 16

__global__ void matmul(const float* A, const float* B, float* C, int n) {
  int r = blockIdx.y * T + threadIdx.y, c = blockIdx.x * T + threadIdx.x;
  float s = 0.0f;
  for (int k = 0; k < n; k++) s += A[r * n + k] * B[k * n + c];  // 2 global reads per multiply-add
  C[r * n + c] = s;
}

int main() {
  const int n = 1024; float *A, *B, *C;
  cudaMallocManaged(&A, n * n * 4); cudaMallocManaged(&B, n * n * 4); cudaMallocManaged(&C, n * n * 4);
  for (int i = 0; i < n * n; i++) { A[i] = i % 7; B[i] = i % 5; }
  dim3 block(T, T), grid(n / T, n / T);
  matmul<<<grid, block>>>(A, B, C, n);  // warm-up: moves the managed pages to the GPU
  cudaEvent_t a, b; cudaEventCreate(&a); cudaEventCreate(&b);
  cudaEventRecord(a); for (int r = 0; r < 10; r++) matmul<<<grid, block>>>(A, B, C, n);
  cudaEventRecord(b); cudaEventSynchronize(b);
  float ms; cudaEventElapsedTime(&ms, a, b); ms /= 10;  // average of 10 launches
  double ref = 0; for (int k = 0; k < n; k++) ref += A[5 * n + k] * B[k * n + 7];
  printf("naive: %.2f ms, %.0f GFLOP/s, C[5][7] = %.0f, CPU says %.0f\n",
         ms, 2.0 * n * n * n / ms / 1e6, C[5 * n + 7], ref);
}

$ nvcc -o naive naive.cu && ./naive
naive: 6.33 ms, 339 GFLOP/s, C[5][7] = 6144, CPU says 6144

The answer is right. Each of the million threads does 1024 multiply-adds, and every one of them reads both operands from global memory: over two billion loads in total, although neighbouring threads read many of the same values. Across four runs the average time per launch was 4.2 to 6.3 ms. Even the best of those is about 516 GFLOP/s, around 6% of the T4’s 8.1 TFLOP/s FP32 peak.

Tiles in shared memory

The block now copies a 16 x 16 tile of A and a tile of B into __shared__ arrays, then every thread does its 16 multiply-adds from there before moving on to the next pair of tiles.

// tiled.cu - the same multiply, staging 16 x 16 tiles in shared memory
#include <cstdio>
#include <cuda_runtime.h>
#define T 16

__global__ void matmul(const float* A, const float* B, float* C, int n) {
  __shared__ float As[T][T], Bs[T][T];
  int tx = threadIdx.x, ty = threadIdx.y, r = blockIdx.y * T + ty, c = blockIdx.x * T + tx;
  float s = 0.0f;
  for (int t = 0; t < n; t += T) {
    As[ty][tx] = A[r * n + t + tx];  // each thread loads one element of each tile
    Bs[ty][tx] = B[(t + ty) * n + c];
    __syncthreads();                 // the whole tile is loaded before anyone reads it
    for (int k = 0; k < T; k++) s += As[ty][k] * Bs[k][tx];
    __syncthreads();                 // everyone is done before the next tile overwrites it
  }
  C[r * n + c] = s;
}

int main() {
  const int n = 1024; float *A, *B, *C;
  cudaMallocManaged(&A, n * n * 4); cudaMallocManaged(&B, n * n * 4); cudaMallocManaged(&C, n * n * 4);
  for (int i = 0; i < n * n; i++) { A[i] = i % 7; B[i] = i % 5; }
  dim3 block(T, T), grid(n / T, n / T);
  matmul<<<grid, block>>>(A, B, C, n);
  cudaEvent_t a, b; cudaEventCreate(&a); cudaEventCreate(&b);
  cudaEventRecord(a); for (int r = 0; r < 10; r++) matmul<<<grid, block>>>(A, B, C, n);
  cudaEventRecord(b); cudaEventSynchronize(b);
  float ms; cudaEventElapsedTime(&ms, a, b); ms /= 10;  // average of 10 launches
  double ref = 0; for (int k = 0; k < n; k++) ref += A[5 * n + k] * B[k * n + 7];
  printf("tiled: %.2f ms, %.0f GFLOP/s, C[5][7] = %.0f, CPU says %.0f\n",
         ms, 2.0 * n * n * n / ms / 1e6, C[5 * n + 7], ref);
}

$ nvcc -o tiled tiled.cu && ./tiled
tiled: 3.88 ms, 554 GFLOP/s, C[5][7] = 6144, CPU says 6144

Each element is now loaded from global memory once per block, and then read 16 times from shared memory, which is on the chip. The two __syncthreads() calls are barriers for the whole block. The first makes sure the tile is complete before anyone uses it; the second stops a fast warp from loading the next tile over data that slower warps are still reading.

Across four runs the tiled kernel took 3.2 to 3.9 ms. That looks faster than the naive one, but the two came from separate processes, and the T4’s clock speed changes from moment to moment. For a fair comparison, measure both under fixed conditions.

A fair comparison with Nsight Compute

ncu profiles one launch of each program (-s 1 skips the warm-up) and reports the kernel’s time and how many global load requests it issued.

$ for p in naive tiled; do echo "$p:"; ncu --metrics gpu__time_duration.sum,l1tex__t_requests_pipe_lsu_mem_global_op_ld.sum -s 1 -c 1 ./$p | grep -E 'duration|requests'; done
naive:
    gpu__time_duration.sum                                   ms         9.18
    l1tex__t_requests_pipe_lsu_mem_global_op_ld.sum               67,108,864
tiled:
    gpu__time_duration.sum                                   ms         5.77
    l1tex__t_requests_pipe_lsu_mem_global_op_ld.sum                4,194,304

ncu locks the GPU to its base clock while it measures, so both times are slower than in the event timings above, but they are comparable. Three runs gave 9.18 to 9.21 ms for the naive kernel and 5.76 to 5.77 ms for the tiled one: about 1.6 times faster.

The request counts explain the change, and they’re exact. The naive kernel makes 2 x 1024³ loads, and one request serves a whole warp of 32 threads: 2,147,483,648 / 32 = 67,108,864. The tiled kernel needs 16 times fewer, 4,194,304, because each element is fetched once per tile instead of once per use. The time fell by 1.6 times, not 16, so global load requests weren’t the only cost: the tiled loop adds its own shared-memory reads and two barriers per tile.

Forget one __syncthreads()

This is the tiled kernel with the first barrier commented out. It checks every 64th row of the result against the CPU.

// nosync.cu - the tiled multiply with the first __syncthreads() removed
#include <cstdio>
#include <cuda_runtime.h>
#define T 16

__global__ void matmul(const float* A, const float* B, float* C, int n) {
  __shared__ float As[T][T], Bs[T][T];
  int tx = threadIdx.x, ty = threadIdx.y, r = blockIdx.y * T + ty, c = blockIdx.x * T + tx;
  float s = 0.0f;
  for (int t = 0; t < n; t += T) {
    As[ty][tx] = A[r * n + t + tx];
    Bs[ty][tx] = B[(t + ty) * n + c];
    // missing: __syncthreads();  other warps may not have stored their part yet
    for (int k = 0; k < T; k++) s += As[ty][k] * Bs[k][tx];
    __syncthreads();
  }
  C[r * n + c] = s;
}

int main() {
  const int n = 1024; float *A, *B, *C;
  cudaMallocManaged(&A, n * n * 4); cudaMallocManaged(&B, n * n * 4); cudaMallocManaged(&C, n * n * 4);
  for (int i = 0; i < n * n; i++) { A[i] = i % 7; B[i] = i % 5; }
  matmul<<<dim3(n / T, n / T), dim3(T, T)>>>(A, B, C, n);
  cudaDeviceSynchronize();
  int wrong = 0;
  for (int r = 0; r < n; r += 64) for (int c = 0; c < n; c++) {  // check every 64th row on the CPU
    double ref = 0; for (int k = 0; k < n; k++) ref += A[r * n + k] * B[k * n + c];
    wrong += C[r * n + c] != ref;
  }
  printf("wrong results in %d of %d checked elements\n", wrong, (n / 64) * n);
}

$ nvcc -o nosync nosync.cu && ./nosync
wrong results in 16082 of 16384 checked elements

Almost every checked element is wrong, and nothing reported an error. Without the barrier, a warp starts multiplying as soon as its own two loads are stored. The Bs[k][tx] values it reads come from the other seven warps of the block, and some of them still hold the previous tile. Four runs gave 16,082, 16,053, 15,958 and 15,991 wrong elements: a race gives a different answer each time, which is often the first clue that a barrier is missing.

Bank conflicts, and padding them away

Shared memory is divided into 32 banks, each 4 bytes wide. This transpose writes a 32 x 32 tile by rows and reads it back by columns; the only difference between the two versions is one extra column.

// banks.cu - a shared-memory transpose, with and without a padding column
#include <cstdio>
#include <cuda_runtime.h>

template <int PAD> __global__ void transpose(const float* in, float* out, int n) {
  __shared__ float tile[32][32 + PAD];
  int x = blockIdx.x * 32 + threadIdx.x, y = blockIdx.y * 32 + threadIdx.y;
  for (int j = 0; j < 32; j += 8) tile[threadIdx.y + j][threadIdx.x] = in[(y + j) * n + x];
  __syncthreads();
  x = blockIdx.y * 32 + threadIdx.x; y = blockIdx.x * 32 + threadIdx.y;
  for (int j = 0; j < 32; j += 8) out[(y + j) * n + x] = tile[threadIdx.x][threadIdx.y + j];  // reads a column
}

template <int PAD> float run(const float* in, float* out, int n) {
  cudaEvent_t a, b; cudaEventCreate(&a); cudaEventCreate(&b);
  transpose<PAD><<<dim3(n / 32, n / 32), dim3(32, 8)>>>(in, out, n);  // warm-up
  cudaEventRecord(a);
  for (int r = 0; r < 10; r++) transpose<PAD><<<dim3(n / 32, n / 32), dim3(32, 8)>>>(in, out, n);
  cudaEventRecord(b); cudaEventSynchronize(b);
  float ms; cudaEventElapsedTime(&ms, a, b); return ms / 10;  // average of 10 launches
}

int main() {
  const int n = 4096; float *in, *out;
  cudaMalloc(&in, n * n * 4); cudaMalloc(&out, n * n * 4);
  float t0 = run<0>(in, out, n), t1 = run<1>(in, out, n);
  printf("tile[32][32]: %.3f ms, %.0f GB/s\n", t0, 2.0 * n * n * 4 / t0 / 1e6);
  printf("tile[32][33]: %.3f ms, %.0f GB/s\n", t1, 2.0 * n * n * 4 / t1 / 1e6);
}

$ nvcc -o banks banks.cu && ./banks
tile[32][32]: 1.013 ms, 133 GB/s
tile[32][33]: 0.693 ms, 194 GB/s

In tile[32][32] each row is exactly 32 floats, so all the elements of a column sit in the same bank. When a warp reads tile[threadIdx.x][...], its 32 threads hit one bank, and the hardware serves them one after another. With tile[32][33], each row starts one bank further along, so a column is spread over all 32 banks and the warp’s read is served at once.

Across four runs the padded version took 0.68 to 0.70 ms (about 190 to 197 GB/s), and the unpadded one 0.84 to 1.22 ms. It was slower every time, though by how much varied. Practice 2 counts the conflicts directly.

Warp shuffles: no shared memory at all

The 32 threads of a warp can read each other’s registers with __shfl_down_sync. Here a warp sums the numbers 0 to 31.

// shuffle.cu - a warp sums 32 values with no shared memory and no __syncthreads()
#include <cstdio>
#include <cuda_runtime.h>

__global__ void warp_sum() {
  int lane = threadIdx.x;  // 0 to 31: this kernel is launched as one warp
  int v = lane;
  for (int offset = 16; offset > 0; offset /= 2) {
    v += __shfl_down_sync(0xffffffff, v, offset);  // add the v held by lane + offset
    if (lane == 0) printf("offset %2d: lane 0 holds %d\n", offset, v);
  }
  int total = __shfl_sync(0xffffffff, v, 0);  // every lane reads lane 0's v
  if (lane == 31) printf("lane 31 sees the total: %d\n", total);
}

int main() {
  warp_sum<<<1, 32>>>();
  cudaDeviceSynchronize();
}

$ nvcc -o shuffle shuffle.cu && ./shuffle
offset 16: lane 0 holds 16
offset  8: lane 0 holds 48
offset  4: lane 0 holds 112
offset  2: lane 0 holds 240
offset  1: lane 0 holds 496
lane 31 sees the total: 496

At offset 16, lane 0 adds lane 16’s value (0 + 16). At offset 8 it adds lane 8’s value, which is by then 8 + 24 = 32, giving 48. Each step halves the distance, and after five steps lane 0 holds 0 + 1 + … + 31 = 496. The other lanes end up with partial or meaningless values, because a shuffle from beyond lane 31 returns the caller’s own value; __shfl_sync then copies lane 0’s total to every lane, so lane 31 sees it too.

No __syncthreads() is needed, because a shuffle only moves data within the warp, and the 0xffffffff mask says all 32 lanes take part. This is the usual first stage of a fast reduction; the next lesson builds a full one.

Common gotchas

• A missing barrier. Any thread that reads shared memory written by another warp needs a __syncthreads() between the write and the read, and another before the data is overwritten.
• __syncthreads() inside a branch. Every thread of the block must reach the same barrier. If some threads skip it, the block can hang or behave unpredictably.
• Bank conflicts from column access. A float[32][32] array read by column serializes a warp’s access. Pad the row to 33 floats.
• Timing different programs. GPU clocks change, so separate runs aren’t directly comparable. Time both kernels the same way, or compare them under ncu.
• Shuffles across warps. __shfl_*_sync only reaches lanes of the same warp. Combining warps needs shared memory or atomics.

Practice

1. Change the tile size to 32 x 32 (1024 threads per block). Is it faster? Check with the events and with ncu.

// tiled32.cu - the tiled multiply with 32 x 32 tiles (1024 threads per block)
#include <cstdio>
#include <cuda_runtime.h>
#define T 32

__global__ void matmul(const float* A, const float* B, float* C, int n) {
  __shared__ float As[T][T], Bs[T][T];
  int tx = threadIdx.x, ty = threadIdx.y, r = blockIdx.y * T + ty, c = blockIdx.x * T + tx;
  float s = 0.0f;
  for (int t = 0; t < n; t += T) {
    As[ty][tx] = A[r * n + t + tx];
    Bs[ty][tx] = B[(t + ty) * n + c];
    __syncthreads();
    for (int k = 0; k < T; k++) s += As[ty][k] * Bs[k][tx];
    __syncthreads();
  }
  C[r * n + c] = s;
}

int main() {
  const int n = 1024; float *A, *B, *C;
  cudaMallocManaged(&A, n * n * 4); cudaMallocManaged(&B, n * n * 4); cudaMallocManaged(&C, n * n * 4);
  for (int i = 0; i < n * n; i++) { A[i] = i % 7; B[i] = i % 5; }
  dim3 block(T, T), grid(n / T, n / T);
  matmul<<<grid, block>>>(A, B, C, n);
  cudaEvent_t a, b; cudaEventCreate(&a); cudaEventCreate(&b);
  cudaEventRecord(a); for (int r = 0; r < 10; r++) matmul<<<grid, block>>>(A, B, C, n);
  cudaEventRecord(b); cudaEventSynchronize(b);
  float ms; cudaEventElapsedTime(&ms, a, b); ms /= 10;  // average of 10 launches
  double ref = 0; for (int k = 0; k < n; k++) ref += A[5 * n + k] * B[k * n + 7];
  printf("tiled 32: %.2f ms, %.0f GFLOP/s, C[5][7] = %.0f, CPU says %.0f\n",
         ms, 2.0 * n * n * n / ms / 1e6, C[5 * n + 7], ref);
}

$ nvcc -o tiled32 tiled32.cu && ./tiled32
tiled 32: 4.38 ms, 490 GFLOP/s, C[5][7] = 6144, CPU says 6144

$ ncu --metrics gpu__time_duration.sum,l1tex__t_requests_pipe_lsu_mem_global_op_ld.sum -s 1 -c 1 ./tiled32 | grep -E 'duration|requests'
    gpu__time_duration.sum                                   ms         5.27
    l1tex__t_requests_pipe_lsu_mem_global_op_ld.sum                2,097,152

2. Use Nsight Compute to count the shared-memory bank conflicts in both transpose kernels. The -s 11 skips the warm-up and the ten timed launches of the first kernel.

$ ncu --metrics l1tex__data_bank_conflicts_pipe_lsu_mem_shared_op_ld.sum -c 1 ./banks | grep -E 'void transpose|bank_conflicts' && ncu --metrics l1tex__data_bank_conflicts_pipe_lsu_mem_shared_op_ld.sum -s 11 -c 1 ./banks | grep -E 'void transpose|bank_conflicts'
  void transpose<0>(const float *, float *, int) (128, 128, 1)x(32, 8, 1), Context 1, Stream 7, Device 0, CC 7.5
    l1tex__data_bank_conflicts_pipe_lsu_mem_shared_op_ld.sum               16,252,928
  void transpose<1>(const float *, float *, int) (128, 128, 1)x(32, 8, 1), Context 1, Stream 7, Device 0, CC 7.5
    l1tex__data_bank_conflicts_pipe_lsu_mem_shared_op_ld.sum                        0

3. Find the maximum of 32 values with __shfl_xor_sync instead, so that every lane ends up with the answer, with no broadcast step.

// xormax.cu - a butterfly with __shfl_xor_sync leaves the maximum in every lane
#include <cstdio>
#include <cuda_runtime.h>

__global__ void warp_max() {
  int lane = threadIdx.x;
  int v = (lane * 13) % 29;  // 32 scattered values between 0 and 28
  for (int mask = 16; mask > 0; mask /= 2)
    v = max(v, __shfl_xor_sync(0xffffffff, v, mask));  // swap with lane ^ mask
  if (lane == 0 || lane == 17 || lane == 31) printf("lane %2d: max = %d\n", lane, v);
}

int main() {
  warp_max<<<1, 32>>>();
  cudaDeviceSynchronize();
}

$ nvcc -o xormax xormax.cu && ./xormax
lane  0: max = 28
lane 17: max = 28
lane 31: max = 28

4. Launch the shuffle sum with 64 threads, giving each thread its own threadIdx.x as the value. Do you get the sum of 0 to 63?

// twowarps.cu - the shuffle sum launched with 64 threads: two warps, two totals
#include <cstdio>
#include <cuda_runtime.h>

__global__ void warp_sum() {
  int lane = threadIdx.x % 32, warp = threadIdx.x / 32;
  int v = threadIdx.x;  // 0 to 63
  for (int offset = 16; offset > 0; offset /= 2)
    v += __shfl_down_sync(0xffffffff, v, offset);
  if (lane == 0) printf("warp %d: lane 0 holds %d\n", warp, v);
}

int main() {
  warp_sum<<<1, 64>>>();
  cudaDeviceSynchronize();
  printf("sum of 0..63 is %d\n", 63 * 64 / 2);
}

$ nvcc -o twowarps twowarps.cu && ./twowarps
warp 1: lane 0 holds 1520
warp 0: lane 0 holds 496
sum of 0..63 is 2016