The Memory Hierarchy in Anger: Matrix Transpose
The Memory Hierarchy in Anger: Matrix Transpose
Writing fast GPU kernels is fundamentally about choreographing data movement through the memory hierarchy. While modern GPUs boast tens of teraflops of arithmetic power, every floating-point operation is useless if data cannot reach the functional units.
In this chapter, we dissect the GPU memory hierarchy, master global memory coalescing, resolve shared memory bank conflicts, and walk through the canonical optimization benchmark: Matrix Transpose.
1. The GPU Memory Hierarchy
A modern GPU features a steeply tiered memory architecture. Each level represents a compromise between capacity, latency, and aggregate bandwidth:
+--------------------------------------------------------------+
| Registers: 64K x 32-bit registers per SM (~0-1 cycle latency)|
| Bandwidth: >30 TB/s aggregate across all SMs |
+--------------------------------------------------------------+
|
+--------------------------------------------------------------+
| Shared Memory / L1 Data Cache: 64 - 228 KB per SM (~19-32 cyc)|
| Bandwidth: ~19 TB/s aggregate (Software-managed SRAM) |
+--------------------------------------------------------------+
|
+--------------------------------------------------------------+
| L2 Cache: 40 - 50 MB on-chip shared across all SMs (~200 cyc)|
| Bandwidth: ~5.2 TB/s aggregate |
+--------------------------------------------------------------+
|
+--------------------------------------------------------------+
| Global Memory (HBM3 / GDDR6): 16 - 96 GB off-chip (~400 cyc)|
| Bandwidth: 1.5 - 3.35 TB/s |
+--------------------------------------------------------------+Memory Space Properties
| Memory Space | Scope | Lifetime | Physical Location | Speed |
|---|---|---|---|---|
| Registers | Single Thread | Thread duration | On-chip Register File | Single-cycle (~0 ns) |
| Shared Memory | Single Thread Block | Block duration | On-chip SRAM | ~15–30 cycles (~1.5 ns) |
| L1 Cache | Streaming Multiprocessor | Hardware managed | On-chip SRAM | ~15–30 cycles |
| L2 Cache | All SMs on GPU | Hardware managed | On-chip Crossbar | ~150–200 cycles (~15 ns) |
| Global Memory | All threads + Host | Application duration | Off-chip DRAM (HBM/GDDR) | ~400+ cycles (~100 ns) |
2. Global Memory Coalescing: The 10x Lever
When threads in a warp issue a read or write instruction to global memory, the memory controller does not process 32 separate individual 4-byte requests. Instead, it aggregates adjacent addresses into coalesced 32-byte or 128-byte memory transactions.
The Golden Rule of Coalescing
Consecutive threads in a warp (
threadIdx.xfrom 0 to 31) must access consecutive memory addresses (address + 0,address + 4, …,address + 124).
Coalesced Access (Consecutive):
Threads: T0 T1 T2 T3 ... T31
| | | | |
Memory Addr: 0x00 0x04 0x08 0x0C ... 0x7C
Result: ONE single 128-byte transaction! (100% bus utilization)
Uncoalesced Strided Access (Column-wise):
Threads: T0 T1 T2 ... T31
| | | |
Memory Addr: 0x0000 0x1000 0x2000 ... 0x1F000 (Stride = N elements)
Result: 32 SEPARATE 32-byte transactions! (3.1% bus utilization!)In the uncoalesced case, the GPU fetches \(32 \times 32 \text{ Bytes} = 1024 \text{ Bytes}\) over the physical memory bus to retrieve only \(32 \times 4 \text{ Bytes} = 128 \text{ Bytes}\) of useful data. 90% of your memory bandwidth is wasted.
4. Case Study: Step-by-Step Optimization of Matrix Transpose
Matrix transpose (\(B_{j,i} = A_{i,j}\)) is the ultimate memory benchmark: if you read \(A\) row-wise (coalesced), you must write \(B\) column-wise (uncoalesced), or vice-versa.
Achieving near-peak memory bandwidth requires staging tiles through shared memory.
Step 1: Naive Transpose (Uncoalesced Writes)
__global__ void transpose_naive(const float* __restrict__ in,
float* __restrict__ out,
int width, int height) {
int x = blockIdx.x * blockDim.x + threadIdx.x; // column
int y = blockIdx.y * blockDim.y + threadIdx.y; // row
if (x < width && y < height) {
// Coalesced Read: Adjacent threads have adjacent x (consecutive in memory)
float val = in[y * width + x];
// UNCOALESCED Write: Adjacent threads have adjacent y, strided by height!
out[x * height + y] = val;
}
}Benchmark Result on A100: Reaches only ~18% of peak global memory bandwidth due to serialized uncoalesced writes.
Full Production Code
#include <cuda_runtime.h>
#include <stdio.h>
#define TILE_DIM 32
#define BLOCK_ROWS 8
__global__ void transpose_optimal(const float* __restrict__ in,
float* __restrict__ out,
int width, int height) {
// 32 x 33 shared memory array: +1 padding eliminates all bank conflicts
__shared__ float tile[TILE_DIM][TILE_DIM + 1];
int x = blockIdx.x * TILE_DIM + threadIdx.x;
int y = blockIdx.y * TILE_DIM + threadIdx.y;
// 1. Coalesced Global Read -> Conflict-Free Shared Write
#pragma unroll
for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS) {
if (x < width && (y + j) < height) {
tile[threadIdx.y + j][threadIdx.x] = in[(y + j) * width + x];
}
}
__syncthreads();
// Transpose tile destination mapping
x = blockIdx.y * TILE_DIM + threadIdx.x;
y = blockIdx.x * TILE_DIM + threadIdx.y;
// 2. Conflict-Free Shared Read -> Coalesced Global Write
#pragma unroll
for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS) {
if (x < height && (y + j) < width) {
out[(y + j) * height + x] = tile[threadIdx.x][threadIdx.y + j];
}
}
}5. Performance Comparison
Benchmarked on an NVIDIA A100 GPU (\(4096 \times 4096\) FP32 Matrix, Peak HBM Bandwidth = \(2,039 \text{ GB/s}\)):
| Implementation | Effective Bandwidth | % of Peak HBM Bandwidth | Bottleneck / Limiter |
|---|---|---|---|
| Naive Transpose | 358 GB/s | 17.5% | Severe uncoalesced global write fragmentation |
| Coalesced with Bank Conflicts | 1,220 GB/s | 59.8% | 32-way shared memory bank conflict serialization |
| Padded Conflict-Free Transpose | 1,845 GB/s | 90.5% | Near-hardware speed of light saturation |
To verify bank conflicts in Nsight Compute:
ncu --metrics l1tex__data_bank_conflicts_pipe_lsu_mem_shared.sum ./transpose_testFor the padded kernel, the conflict count drops from hundreds of thousands to 0.