__syncthreads(), warp-level primitives, atomic operations, race conditions, and parallel reduction patterns.
This tutorial explores how CUDA threads coordinate their work. We cover barriers, warp-level intrinsics, atomic operations, and put them all together with parallel reduction and histogram patterns.
Tutorials 01–04 (GPU Architecture, First Kernel, Thread Hierarchy, Memory Model). You should be comfortable writing kernels that use shared memory and understand the thread/block/grid hierarchy.
CUDA launches thousands of threads that execute concurrently. Without coordination, threads reading and writing shared data produce data races — the results depend on unpredictable execution order.
Consider two threads both trying to increment a shared counter:
Multiple threads read/write the same location without ordering. The final value is undefined.
Lost updates, partial writes, or stale reads lead to wrong answers that may look plausible.
Results change between runs. The bug may appear intermittently, making debugging extremely difficult.
| Mechanism | Scope | Cost |
|---|---|---|
__syncthreads() |
All threads in a block | Low (hardware barrier) |
| Warp-level primitives | 32 threads in a warp | Very low (single instruction) |
| Atomic operations | Global or shared memory | Medium (serialises at address) |
| Cooperative Groups | Flexible (warp to grid) | Varies by scope |
Synchronisation always has a cost. The art of CUDA programming is choosing the narrowest scope that guarantees correctness — warp-level if possible, block-level if needed, grid-level as a last resort.
__syncthreads() is a barrier: every thread in the block must reach this point before any thread is allowed to proceed past it. It also acts as a memory fence — all shared memory writes before the barrier are visible to all threads after it.
__shared__ float tile[BLOCK_SIZE];
// Phase 1: Each thread loads one element
tile[threadIdx.x] = input[blockIdx.x * blockDim.x + threadIdx.x];
// BARRIER: ensure all loads complete before any reads
__syncthreads();
// Phase 2: Now safe to read any element in the tile
float left = (threadIdx.x > 0) ? tile[threadIdx.x - 1] : 0.0f;
float right = (threadIdx.x < BLOCK_SIZE - 1) ? tile[threadIdx.x + 1] : 0.0f;
output[idx] = 0.25f * left + 0.5f * tile[threadIdx.x] + 0.25f * right;If __syncthreads() appears inside a conditional branch, all threads in the block must still reach it. If some threads take the branch and others don't, the block deadlocks — threads that reached the barrier wait forever for threads that will never arrive.
if (threadIdx.x < 32) {
smem[threadIdx.x] = data;
__syncthreads(); // Threads 32+ never reach here!
}if (threadIdx.x < 32) {
smem[threadIdx.x] = data;
}
__syncthreads(); // All threads reach this__syncthreads() must be reached by every thread in the block, or by no thread. Never place it inside an if that excludes some threads. Also note: it only synchronises within a single block — there is no built-in way to synchronise across blocks (until Cooperative Groups).
Since all 32 threads in a warp execute in lockstep (SIMT), CUDA provides warp-level intrinsics that let threads communicate directly through registers — no shared memory needed, no barrier needed.
All warp primitives take a mask parameter — a 32-bit integer where each bit indicates which threads participate. Use 0xFFFFFFFF for all 32 threads, or a custom mask for a subset.
Directly read a register value from another thread in the warp. Four variants:
| Intrinsic | Description | Use Case |
|---|---|---|
__shfl_sync(mask, val, src) |
Read val from lane src |
Broadcast a value to all lanes |
__shfl_up_sync(mask, val, delta) |
Read from lane id - delta |
Prefix sums (scan) |
__shfl_down_sync(mask, val, delta) |
Read from lane id + delta |
Reductions |
__shfl_xor_sync(mask, val, mask) |
Read from lane id ^ mask |
Butterfly reductions |
__device__ float warpReduceSum(float val) {
for (int offset = 16; offset > 0; offset >>= 1) {
val += __shfl_down_sync(0xFFFFFFFF, val, offset);
}
return val; // Result valid in lane 0 only
}This completes a full 32-element reduction in 5 instructions, entirely in registers:
| Intrinsic | Returns | Use Case |
|---|---|---|
__ballot_sync(mask, predicate) |
Bitmask of threads where predicate is true | Count matching threads, compact |
__all_sync(mask, predicate) |
1 if all threads' predicate is true | Early exit optimisation |
__any_sync(mask, predicate) |
1 if any thread's predicate is true | Detect if work remains |
__popc(__ballot_sync(...)) |
Population count of ballot result | Count how many threads match |
Warp-level primitives are the fastest synchronisation mechanism in CUDA. Each compiles to a single hardware instruction. Always prefer shuffle-based reductions over shared memory when operating within a single warp.
Atomic operations perform a read-modify-write as a single, indivisible operation. No other thread can see the intermediate state. They solve the data race problem from Slide 01 — but at a performance cost.
| Function | Operation | Types |
|---|---|---|
atomicAdd(&addr, val) |
*addr += val |
int, unsigned, float, double* |
atomicSub(&addr, val) |
*addr -= val |
int, unsigned |
atomicMax(&addr, val) |
*addr = max(*addr, val) |
int, unsigned |
atomicMin(&addr, val) |
*addr = min(*addr, val) |
int, unsigned |
atomicExch(&addr, val) |
*addr = val (returns old) |
int, unsigned, float |
atomicCAS(&addr, compare, val) |
If *addr == compare, set to val |
int, unsigned, unsigned long long |
* atomicAdd for double requires compute capability 6.0+
Compare-And-Swap is the most powerful atomic: it can implement any atomic operation. Here's how to build an atomic double-precision add on older hardware:
__device__ double atomicAddDouble(double* addr, double val) {
unsigned long long* addr_as_ull = (unsigned long long*)addr;
unsigned long long old = *addr_as_ull, assumed;
do {
assumed = old;
old = atomicCAS(addr_as_ull, assumed,
__double_as_longlong(__longlong_as_double(assumed) + val));
} while (assumed != old);
return __longlong_as_double(old);
}Atomics serialise at the address level. If 1000 threads all atomicAdd to the same global address, they execute one at a time. Solution: reduce locally first (within warp, then block), then perform a single atomic per block to global memory.
Reduction is the canonical example for synchronisation: combine N values into one (sum, max, min, etc.). Let's build it from naive to optimal.
O(N) steps, 1 thread. No parallelism at all.
O(log N) steps, N/2 threads. Each step halves active threads.
Tree in shared memory, then warp shuffle for final 32 elements. Fastest.
Eight values reduced in 3 steps (log₂8 = 3):
#include <cstdio>
#include <cstdlib>
#include <cuda_runtime.h>
#define BLOCK_SIZE 256
// ── Warp-level reduction using shuffle ──
__device__ float warpReduceSum(float val) {
for (int offset = 16; offset > 0; offset >>= 1)
val += __shfl_down_sync(0xFFFFFFFF, val, offset);
return val;
}
// ── Block-level reduction: shared memory tree + warp shuffle finish ──
__device__ float blockReduceSum(float val) {
__shared__ float warp_sums[32]; // max 32 warps per block (1024/32)
int lane = threadIdx.x & 31; // thread index within warp (0-31)
int wid = threadIdx.x >> 5; // warp index within block
// Step 1: reduce within each warp
val = warpReduceSum(val);
// Step 2: lane 0 of each warp writes its result to shared memory
if (lane == 0) warp_sums[wid] = val;
__syncthreads();
// Step 3: first warp reduces the warp_sums
int num_warps = blockDim.x >> 5;
val = (threadIdx.x < num_warps) ? warp_sums[threadIdx.x] : 0.0f;
if (wid == 0) val = warpReduceSum(val);
return val;
}
// ── Kernel: each block reduces its portion, atomicAdd to global result ──
__global__ void reduceSum(const float* input, float* output, int N) {
float sum = 0.0f;
// Grid-stride loop: each thread accumulates multiple elements
for (int i = blockIdx.x * blockDim.x + threadIdx.x;
i < N;
i += blockDim.x * gridDim.x) {
sum += input[i];
}
// Reduce within block
sum = blockReduceSum(sum);
// Thread 0 of each block adds to global result
if (threadIdx.x == 0)
atomicAdd(output, sum);
}
int main() {
const int N = 1 << 20; // 1M elements
size_t bytes = N * sizeof(float);
// Host allocation and initialisation
float* h_input = (float*)malloc(bytes);
float h_output = 0.0f;
float expected = 0.0f;
for (int i = 0; i < N; i++) {
h_input[i] = 1.0f; // all ones: sum should be N
expected += h_input[i];
}
// Device allocation
float *d_input, *d_output;
cudaMalloc(&d_input, bytes);
cudaMalloc(&d_output, sizeof(float));
// Copy input and zero the output
cudaMemcpy(d_input, h_input, bytes, cudaMemcpyHostToDevice);
cudaMemset(d_output, 0, sizeof(float));
// Launch: use enough blocks to keep the GPU busy
int numBlocks = (N + BLOCK_SIZE - 1) / BLOCK_SIZE;
if (numBlocks > 1024) numBlocks = 1024; // cap for grid-stride loop
reduceSum<<<numBlocks, BLOCK_SIZE>>>(d_input, d_output, N);
// Copy result back
cudaMemcpy(&h_output, d_output, sizeof(float), cudaMemcpyDeviceToHost);
printf("Sum = %.0f (expected %.0f) — %s\n",
h_output, expected,
(fabsf(h_output - expected) < 1.0f) ? "PASS" : "FAIL");
// Cleanup
cudaFree(d_input);
cudaFree(d_output);
free(h_input);
return 0;
}Grid-stride loop → each thread accumulates a partial sum. Warp shuffle → reduces 32 partial sums within each warp. Shared memory + __syncthreads → reduces warp results within each block. atomicAdd → combines block results into a single global sum.
Histograms are a classic use case for atomics: many threads categorise data into bins concurrently. The naive approach (global atomics for every element) is slow. The optimised approach uses shared memory per-block histograms merged into global memory at the end.
#include <cstdio>
#include <cstdlib>
#include <cuda_runtime.h>
#define NUM_BINS 256
#define BLOCK_SIZE 256
// ── Phase 1: per-block shared memory histogram ──
// ── Phase 2: merge per-block histogram into global ──
__global__ void histogram(const unsigned char* data, int N,
unsigned int* globalHist) {
// Shared memory histogram for this block
__shared__ unsigned int smemHist[NUM_BINS];
// Cooperatively zero the shared histogram
for (int i = threadIdx.x; i < NUM_BINS; i += blockDim.x)
smemHist[i] = 0;
__syncthreads();
// Phase 1: each thread processes elements via grid-stride loop
for (int i = blockIdx.x * blockDim.x + threadIdx.x;
i < N;
i += blockDim.x * gridDim.x) {
atomicAdd(&smemHist[data[i]], 1);
}
__syncthreads();
// Phase 2: merge shared histogram into global histogram
for (int i = threadIdx.x; i < NUM_BINS; i += blockDim.x) {
if (smemHist[i] > 0)
atomicAdd(&globalHist[i], smemHist[i]);
}
}
int main() {
const int N = 1 << 22; // 4M bytes
size_t dataBytes = N * sizeof(unsigned char);
size_t histBytes = NUM_BINS * sizeof(unsigned int);
// Host allocations
unsigned char* h_data = (unsigned char*)malloc(dataBytes);
unsigned int* h_hist = (unsigned int*)calloc(NUM_BINS, sizeof(unsigned int));
unsigned int* h_ref = (unsigned int*)calloc(NUM_BINS, sizeof(unsigned int));
// Generate random data and compute CPU reference
srand(42);
for (int i = 0; i < N; i++) {
h_data[i] = (unsigned char)(rand() % NUM_BINS);
h_ref[h_data[i]]++;
}
// Device allocations
unsigned char* d_data;
unsigned int* d_hist;
cudaMalloc(&d_data, dataBytes);
cudaMalloc(&d_hist, histBytes);
cudaMemcpy(d_data, h_data, dataBytes, cudaMemcpyHostToDevice);
cudaMemset(d_hist, 0, histBytes);
// Launch kernel
int numBlocks = (N + BLOCK_SIZE - 1) / BLOCK_SIZE;
if (numBlocks > 1024) numBlocks = 1024;
histogram<<<numBlocks, BLOCK_SIZE>>>(d_data, N, d_hist);
// Copy results back
cudaMemcpy(h_hist, d_hist, histBytes, cudaMemcpyDeviceToHost);
// Verify
int errors = 0;
for (int i = 0; i < NUM_BINS; i++) {
if (h_hist[i] != h_ref[i]) {
printf("Mismatch at bin %d: GPU=%u CPU=%u\n", i, h_hist[i], h_ref[i]);
errors++;
}
}
printf("Histogram: %d/%d bins correct — %s\n",
NUM_BINS - errors, NUM_BINS,
(errors == 0) ? "PASS" : "FAIL");
// Print first 10 bins
printf("First 10 bins: ");
for (int i = 0; i < 10; i++)
printf("[%d]=%u ", i, h_hist[i]);
printf("\n");
// Cleanup
cudaFree(d_data);
cudaFree(d_hist);
free(h_data);
free(h_hist);
free(h_ref);
return 0;
}With 1024 blocks, a naive global-only histogram has thousands of threads contending on each bin. The shared-memory approach reduces contention to at most BLOCK_SIZE threads per bin, and the merge phase only has numBlocks atomic operations per bin on global memory.
Before CUDA 9, synchronisation was limited to __syncthreads() (block-level) and atomics. Cooperative Groups provides a flexible, composable API for synchronising at any granularity — from sub-warp tiles to entire grids.
| Group Type | Scope | How to Create |
|---|---|---|
thread_block |
All threads in a block | this_thread_block() |
tiled_partition<N> |
Sub-warp tile (N = 1, 2, 4, 8, 16, 32) | tiled_partition<16>(block) |
coalesced_group |
Active threads only (after divergence) | coalesced_threads() |
grid_group |
All threads in the grid | this_grid() |
#include <cooperative_groups.h>
namespace cg = cooperative_groups;
__device__ float tileReduceSum(float val, cg::thread_block_tile<32> tile) {
for (int offset = tile.size() / 2; offset > 0; offset /= 2)
val += tile.shfl_down(val, offset);
return val; // valid in thread_rank() == 0
}
__global__ void kernel(float* data) {
auto block = cg::this_thread_block();
auto tile = cg::tiled_partition<32>(block);
float val = data[block.thread_rank()];
float sum = tileReduceSum(val, tile);
if (tile.thread_rank() == 0)
data[tile.meta_group_rank()] = sum;
}For the first time, all blocks in a grid can synchronise without launching a new kernel:
#include <cooperative_groups.h>
namespace cg = cooperative_groups;
__global__ void multiPhaseKernel(float* data, int N) {
auto grid = cg::this_grid();
// Phase 1: all blocks process data
int idx = grid.thread_rank();
if (idx < N) data[idx] *= 2.0f;
grid.sync(); // All blocks must reach here before continuing
// Phase 2: all blocks can safely read Phase 1 results
if (idx < N && idx > 0)
data[idx] += data[idx - 1];
}
// Must launch with cudaLaunchCooperativeKernel:
// void* args[] = { &d_data, &N };
// cudaLaunchCooperativeKernel((void*)multiPhaseKernel,
// numBlocks, blockSize, args);Grid-level sync requires that all blocks fit on the GPU simultaneously (the GPU cannot context-switch blocks). Use cudaOccupancyMaxActiveBlocksPerMultiprocessor to determine the maximum safe grid size. If you exceed it, the kernel will deadlock.
Warp shuffle > shared memory barrier > atomics > grid sync. Use the narrowest scope that guarantees correctness.
Always combine values at the warp and block level before touching global memory. This is the single most important optimisation pattern.
Profiling & Performance Analysis — use Nsight Systems, Nsight Compute, and nvprof to identify bottlenecks, measure occupancy, and systematically optimise your CUDA kernels.