GPU Fundamentals: SIMT & Memory Hierarchy

A GPU runs thousands of threads, but it doesn’t schedule them one by one. It groups them into warps of 32 and issues one instruction to the whole warp. That model is called SIMT (single instruction, multiple threads). Every program below ran on a Tesla T4 with CUDA 13.0, built with plain nvcc.

The hardware you are programming

Before writing a kernel, ask the GPU what it has. cudaGetDeviceProperties reports the sizes of every level of the hierarchy.

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

int main() {
  cudaDeviceProp p;
  cudaGetDeviceProperties(&p, 0);
  printf("GPU:                  %s (compute %d.%d)\n", p.name, p.major, p.minor);
  printf("SMs:                  %d\n", p.multiProcessorCount);
  printf("warp size:            %d threads\n", p.warpSize);
  printf("max threads per SM:   %d (%d warps)\n", p.maxThreadsPerMultiProcessor,
         p.maxThreadsPerMultiProcessor / p.warpSize);
  printf("registers per SM:     %d x 32-bit\n", p.regsPerMultiprocessor);
  printf("shared memory per SM: %zu KB\n", p.sharedMemPerMultiprocessor / 1024);
  printf("L2 cache:             %d KB\n", p.l2CacheSize / 1024);
  printf("global memory (DRAM): %zu MB\n", p.totalGlobalMem / (1024 * 1024));
  return 0;
}

$ nvcc -o hierarchy hierarchy.cu && ./hierarchy
GPU:                  Tesla T4 (compute 7.5)
SMs:                  40
warp size:            32 threads
max threads per SM:   1024 (32 warps)
registers per SM:     65536 x 32-bit
shared memory per SM: 64 KB
L2 cache:             4096 KB
global memory (DRAM): 14912 MB

Read it as a hierarchy, fastest and smallest first:

• Registers: 65,536 per SM, shared out among the threads on that SM. If all 1,024 threads are resident, each gets 64. Use more and fewer threads fit.
• Shared memory: 64 KB per SM, split between the blocks running there. It’s a scratchpad you fill and read yourself.
• L2 cache: 4 MB for the whole chip, managed by the hardware.
• DRAM (global memory): about 15 GB. It’s huge, and it’s the slowest level by far.

The 40 SMs are the GPU’s “cores” in the CPU sense. Each holds up to 32 warps at once and switches between them whenever one waits on memory.

The global index and the bounds guard

Every thread works out which element it owns from its block and thread number. The grid is rounded up to whole blocks, so some threads have nothing to do.

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

__global__ void scale(const float* x, float* y, int n, float a) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i == 255 || i == 256 || i == 999 || i == 1000 || i == 1023)
    printf("block %d, thread %d -> global index %d: %s\n", blockIdx.x, threadIdx.x, i,
           i < n ? "does work" : "stopped by the guard");
  if (i < n) y[i] = a * x[i];
}

int main() {
  const int n = 1000, threads = 256;
  int blocks = (n + threads - 1) / threads;   // round up: 4 blocks
  float *x, *y;
  cudaMallocManaged(&x, n * sizeof(float));
  cudaMallocManaged(&y, n * sizeof(float));
  for (int i = 0; i < n; i++) x[i] = i;
  scale<<<blocks, threads>>>(x, y, n, 2.0f);
  cudaDeviceSynchronize();
  printf("%d blocks x %d threads: %d threads for %d elements; y[999] is %.0f\n",
         blocks, threads, blocks * threads, n, y[999]);
  return 0;
}

$ nvcc -o index index.cu && ./index
block 1, thread 0 -> global index 256: does work
block 0, thread 255 -> global index 255: does work
block 3, thread 231 -> global index 999: does work
block 3, thread 232 -> global index 1000: stopped by the guard
block 3, thread 255 -> global index 1023: stopped by the guard
4 blocks x 256 threads: 1024 threads for 1000 elements; y[999] is 1998

The formula blockIdx.x * blockDim.x + threadIdx.x turns thread 0 of block 1 into element 256. Four blocks of 256 make 1,024 threads for 1,000 elements, so the last 24 threads (1000 to 1023) exist but must not touch memory. Without if (i < n) they would read past the end of x and write past the end of y.

Block 1’s line printed before block 0’s. Blocks run in whatever order the hardware picks, and no order is promised between them. Never write a kernel that relies on block order.

Seeing divergence in the active mask

__activemask() returns a 32-bit mask of the lanes executing alongside the current thread. Inside a branch, it shows exactly which threads of the warp are on that path.

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

__global__ void branches() {
  int warp = threadIdx.x / 32, lane = threadIdx.x % 32;
  if (lane % 2 == 0) {                    // splits every warp in half
    unsigned m = __activemask();
    if (lane == 0) printf("warp %d, even lanes:   active mask 0x%08x\n", warp, m);
  } else {
    unsigned m = __activemask();
    if (lane == 1) printf("warp %d, odd lanes:    active mask 0x%08x\n", warp, m);
  }
  if (warp == 0) {                        // a whole warp takes the same side
    unsigned m = __activemask();
    if (lane == 0) printf("warp %d, uniform test: active mask 0x%08x\n", warp, m);
  }
}

int main() {
  branches<<<1, 64>>>();
  cudaDeviceSynchronize();
  return 0;
}

$ nvcc -o divergence divergence.cu && ./divergence
warp 0, odd lanes:    active mask 0xaaaaaaaa
warp 1, odd lanes:    active mask 0xaaaaaaaa
warp 0, even lanes:   active mask 0x55555555
warp 1, even lanes:   active mask 0x55555555
warp 0, uniform test: active mask 0xffffffff

0x55555555 is binary 0101...0101: only the even lanes are active, and the odd half of the warp sits idle. The odd path runs with 0xaaaaaaaa, the other half. The warp went through both sides of the if, one after the other. Here the else side happened to run first; the compiler and hardware choose the order, and your code must not depend on it.

The third test splits on warp, not lane, so all 32 lanes of warp 0 agree and the mask is full: 0xffffffff. Warp 1 skipped that block entirely and paid nothing for it. Divergence only exists inside a warp.

What divergence costs

Give the two paths real work and time them with CUDA events. The same kernel runs twice: once with the branch splitting every warp, once with each warp on a single side.

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

__global__ void work(float* out, int split) {  // split 1: per thread, 32: per warp
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  float v = i * 1e-6f;
  if ((threadIdx.x / split) % 2 == 0) { for (int k = 0; k < 1000; k++) v = sinf(v) + 0.5f; }
  else                                { for (int k = 0; k < 1000; k++) v = sqrtf(v) + 0.5f; }
  out[i] = v;
}

float time_ms(float* out, int split) {
  cudaEvent_t a, b; cudaEventCreate(&a); cudaEventCreate(&b);
  cudaEventRecord(a);
  work<<<4096, 256>>>(out, split);
  cudaEventRecord(b); cudaEventSynchronize(b);
  float ms; cudaEventElapsedTime(&ms, a, b); return ms;
}

int main() {
  float* out; cudaMalloc(&out, 4096 * 256 * sizeof(float));
  time_ms(out, 32);                                      // warm-up launch
  printf("branch splits every warp:   %.2f ms\n", time_ms(out, 1));
  printf("branch is uniform per warp: %.2f ms\n", time_ms(out, 32));
  return 0;
}

$ nvcc -o divergence_cost divergence_cost.cu && ./divergence_cost
branch splits every warp:   7.51 ms
branch is uniform per warp: 3.90 ms

Both runs do the same arithmetic: half the threads run the sinf loop and half the sqrtf loop. The only change is which threads. With split = 1, neighbouring lanes disagree, so every warp runs both loops back to back and the kernel takes about twice as long. With split = 32, each warp runs just one loop.

Your times will differ. Across repeated runs on the same T4 they ranged from 6.5 to 14.7 ms for the divergent version, because the GPU’s clock speed moves with temperature and load. The ratio held steady at about 1.9x. Compare ratios, not single timings.

The first launch is a warm-up. It pays one-off setup costs that would otherwise land on the first measurement.

Common gotchas

• No bounds guard. The grid rounds up, so the last block nearly always has spare threads.
• Expecting an order between blocks. The block order in index.cu’s output isn’t guaranteed, and it can change from run to run or GPU to GPU.
• Branching on threadIdx.x % 2-style conditions in hot loops. Branch on something constant across a warp (threadIdx.x / 32, blockIdx.x) when you have the choice.
• Timing the first launch. It includes one-time setup. Warm up first.

Practice

1. In index.cu, change threads to 128. Predict the number of blocks and spare threads before you run it.

2. Remove the if (i < n) guard in index.cu and run it under compute-sanitizer ./index. What does the tool report?

3. In divergence.cu, change the first test to lane < 16. What masks do you expect?

4. Run divergence_cost with split set to 16 instead of 1. Is it fast or slow, and why?