CUDA Atomics & Parallel Reductions

When thousands of threads update the same value, ordinary arithmetic breaks: updates are lost, and nothing reports an error. Atomics make an update indivisible, but they serialize when threads compete for one address. The way out is a reduction: combine values inside each block or warp first, and only touch the shared result rarely. Every program ran on a Tesla T4 with CUDA 13.0.

A race on one counter

Every thread in this program adds 1 to the same integer: once with +=, once with atomicAdd.

// race.cu - a million threads add 1 to the same counter
#include <cstdio>
#include <cuda_runtime.h>

__global__ void plain(int* c) { *c += 1; }           // read, add, write: three separate steps
__global__ void atomic(int* c) { atomicAdd(c, 1); }  // one indivisible read-modify-write

int main() {
  int* c; cudaMallocManaged(&c, sizeof(int));
  *c = 0; plain<<<4096, 256>>>(c); cudaDeviceSynchronize();
  printf("*c += 1:      %d\n", *c);
  *c = 0; atomic<<<4096, 256>>>(c); cudaDeviceSynchronize();
  printf("atomicAdd:    %d\n", *c);
  printf("threads:      %d\n", 4096 * 256);
  cudaFree(c);
}

$ nvcc -o race race.cu && ./race
*c += 1:      75
atomicAdd:    1048576
threads:      1048576

*c += 1 is three steps: load the value, add 1, store it. Thousands of threads load the same old value at about the same time, each stores its own +1, and they overwrite each other. Over a million increments survived as 75; three more runs gave 85, 85 and 70. A race gives wrong answers that change from run to run, and no error.

atomicAdd does the load, add and store as one indivisible operation, so no update is lost and the count is exactly 4096 x 256 = 1,048,576.

Contention: when every thread wants the same address

A histogram is the classic use of atomics: each thread adds 1 to the bin its value falls in. The same 16 million values are counted twice, first spread over 256 bins, then with every value in one bin.

// hist.cu - a 256-bin histogram with global atomics: spread out, then all in one bin
#include <cstdio>
#include <cuda_runtime.h>

__global__ void histo(const unsigned* d, unsigned* h, int n) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i < n) atomicAdd(&h[d[i] % 256], 1u);
}

float run(const unsigned* d, unsigned* h, unsigned* out, int n) {
  cudaEvent_t a, b; cudaEventCreate(&a); cudaEventCreate(&b);
  cudaMemset(h, 0, 256 * 4);
  cudaEventRecord(a); histo<<<n / 256, 256>>>(d, h, n); cudaEventRecord(b);
  cudaMemcpy(out, h, 256 * 4, cudaMemcpyDeviceToHost);
  float ms; cudaEventElapsedTime(&ms, a, b); return ms;
}

int main() {
  const int n = 1 << 24; unsigned *x = new unsigned[n], out[256], *d, *h, cpu7 = 0;
  for (int i = 0; i < n; i++) { x[i] = (i * 2654435761u) >> 8; cpu7 += x[i] % 256 == 7; }
  cudaMalloc(&d, n * 4); cudaMalloc(&h, 256 * 4);
  cudaMemcpy(d, x, n * 4, cudaMemcpyHostToDevice);
  run(d, h, out, n);  // warm-up
  float spread = run(d, h, out, n);
  printf("spread over 256 bins: %6.2f ms, bin 7 = %u, CPU says %u\n", spread, out[7], cpu7);
  cudaMemset(d, 0, n * 4);  // now every value is 0: every thread hits bin 0
  float same = run(d, h, out, n);
  printf("all in one bin:       %6.2f ms, bin 0 = %u\n", same, out[0]);
}

$ nvcc -o hist hist.cu && ./hist
spread over 256 bins:   9.72 ms, bin 7 = 65536, CPU says 65536
all in one bin:        21.55 ms, bin 0 = 16777216

Both counts are correct: bin 7 holds 65,536 values, matching the CPU count, and all 16,777,216 values reach bin 0 in the second run. But the work is identical and the time more than doubled, from 9.7 to 21.6 ms (four runs: 8.3 to 11.1 ms, and 18.4 to 24.6 ms). Atomic updates to the same address are handled one at a time, so 16 million of them to a single bin form one long queue.

Summing 16 million numbers with atomics

The same problem appears in a sum: here every thread adds its element straight into one 64-bit total.

// reduce_atomic.cu - sum 16M integers with one global atomicAdd per element
#include <cstdio>
#include <cuda_runtime.h>

__global__ void sum(const unsigned* x, unsigned long long* total, int n) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i < n) atomicAdd(total, (unsigned long long)x[i]);  // 16M updates of one address
}

int main() {
  const int n = 1 << 24; unsigned* h = new unsigned[n]; unsigned long long cpu = 0, total;
  for (int i = 0; i < n; i++) { h[i] = i % 10; cpu += h[i]; }
  unsigned* x; unsigned long long* d_total;
  cudaMalloc(&x, n * 4); cudaMalloc(&d_total, 8);
  cudaMemcpy(x, h, n * 4, cudaMemcpyHostToDevice);
  sum<<<n / 256, 256>>>(x, d_total, n);  // warm-up
  cudaEvent_t a, b; cudaEventCreate(&a); cudaEventCreate(&b);
  cudaMemset(d_total, 0, 8);
  cudaEventRecord(a); sum<<<n / 256, 256>>>(x, d_total, n); cudaEventRecord(b);
  cudaMemcpy(&total, d_total, 8, cudaMemcpyDeviceToHost);  // waits for the kernel
  float ms; cudaEventElapsedTime(&ms, a, b);
  printf("one atomic per element: %.2f ms, sum = %llu, CPU = %llu\n", ms, total, cpu);
}

$ nvcc -o reduce_atomic reduce_atomic.cu && ./reduce_atomic
one atomic per element: 22.97 ms, sum = 75497460, CPU = 75497460

The sum is exact, and it took 23 to 24.6 ms over three runs. That’s 16 million updates to one address, the worst case for contention.

Reduce inside the block first

Instead, each block sums its 256 elements in shared memory, halving the number of active threads each step, and only thread 0 sends the block’s total to global memory.

// reduce_block.cu - reduce each block in shared memory, then one atomicAdd per block
#include <cstdio>
#include <cuda_runtime.h>

__global__ void sum(const unsigned* x, unsigned long long* total, int n) {
  __shared__ unsigned s[256];
  int t = threadIdx.x, i = blockIdx.x * 256 + t;
  s[t] = i < n ? x[i] : 0;
  __syncthreads();
  for (int stride = 128; stride > 0; stride /= 2) {  // halve the active threads each step
    if (t < stride) s[t] += s[t + stride];
    __syncthreads();
  }
  if (t == 0) atomicAdd(total, (unsigned long long)s[0]);  // one atomic for 256 elements
}

int main() {
  const int n = 1 << 24; unsigned* h = new unsigned[n]; unsigned long long cpu = 0, total;
  for (int i = 0; i < n; i++) { h[i] = i % 10; cpu += h[i]; }
  unsigned* x; unsigned long long* d_total;
  cudaMalloc(&x, n * 4); cudaMalloc(&d_total, 8);
  cudaMemcpy(x, h, n * 4, cudaMemcpyHostToDevice);
  sum<<<n / 256, 256>>>(x, d_total, n);  // warm-up
  cudaEvent_t a, b; cudaEventCreate(&a); cudaEventCreate(&b);
  cudaMemset(d_total, 0, 8);
  cudaEventRecord(a); sum<<<n / 256, 256>>>(x, d_total, n); cudaEventRecord(b);
  cudaMemcpy(&total, d_total, 8, cudaMemcpyDeviceToHost);
  float ms; cudaEventElapsedTime(&ms, a, b);
  printf("one atomic per block:   %.2f ms, sum = %llu, CPU = %llu\n", ms, total, cpu);
}

$ nvcc -o reduce_block reduce_block.cu && ./reduce_block
one atomic per block:   1.10 ms, sum = 75497460, CPU = 75497460

The loop takes 8 steps: 128 threads add pairs, then 64, and so on down to 1, with a barrier after each. Only 65,536 atomics reach the total instead of 16 million, and the time dropped to 0.87 to 1.10 ms over three runs, more than 20 times faster.

Reduce inside each warp with shuffles

Warp shuffles, from the previous lesson, sum 32 values without shared memory or barriers. One lane per warp then does the atomic.

// reduce_warp.cu - sum each warp with shuffles, then one atomicAdd per warp
#include <cstdio>
#include <cuda_runtime.h>

__global__ void sum(const unsigned* x, unsigned long long* total, int n) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  unsigned v = i < n ? x[i] : 0;
  for (int offset = 16; offset > 0; offset /= 2) v += __shfl_down_sync(0xffffffff, v, offset);
  if (threadIdx.x % 32 == 0) atomicAdd(total, (unsigned long long)v);  // one atomic per 32 elements
}

int main() {
  const int n = 1 << 24; unsigned* h = new unsigned[n]; unsigned long long cpu = 0, total;
  for (int i = 0; i < n; i++) { h[i] = i % 10; cpu += h[i]; }
  unsigned* x; unsigned long long* d_total;
  cudaMalloc(&x, n * 4); cudaMalloc(&d_total, 8);
  cudaMemcpy(x, h, n * 4, cudaMemcpyHostToDevice);
  sum<<<n / 256, 256>>>(x, d_total, n);  // warm-up
  cudaEvent_t a, b; cudaEventCreate(&a); cudaEventCreate(&b);
  cudaMemset(d_total, 0, 8);
  cudaEventRecord(a); sum<<<n / 256, 256>>>(x, d_total, n); cudaEventRecord(b);
  cudaMemcpy(&total, d_total, 8, cudaMemcpyDeviceToHost);
  float ms; cudaEventElapsedTime(&ms, a, b);
  printf("one atomic per warp:    %.2f ms, sum = %llu, CPU = %llu\n", ms, total, cpu);
}

$ nvcc -o reduce_warp reduce_warp.cu && ./reduce_warp
one atomic per warp:    0.68 ms, sum = 75497460, CPU = 75497460

This version sends 524,288 atomics to the total, eight times as many as the block version, and it was still faster: 0.68 to 0.74 ms over three runs. Fewer atomics isn’t automatically better; here the extra atomics most likely cost less than the block version’s eight rounds of shared memory and __syncthreads(), which this lesson didn’t time separately.

The three versions under Nsight Compute

Event timings from separate processes move around, so this profiles one launch of each with ncu, which runs them at the same locked clock.

$ for p in reduce_atomic reduce_block reduce_warp; do echo "$p:"; ncu --metrics gpu__time_duration.sum -s 1 -c 1 ./$p | grep duration; done
reduce_atomic:
    gpu__time_duration.sum          ms        24.58
reduce_block:
    gpu__time_duration.sum          ms         1.48
reduce_warp:
    gpu__time_duration.sum          us       816.61

The ranking holds: 24.58 ms for one atomic per element, 1.48 ms for one per block and 0.82 ms for one per warp, so the warp version is about 30 times faster than the first. Three runs gave the same values to the precision ncu prints. The block and warp versions are slower here than in their event timings, most likely because ncu holds the GPU at its base clock.

Common gotchas

• A plain += on shared data. Any location that several threads update needs an atomic, or a reduction that gives each thread its own slot.
• Contention. Atomics on one address are serialized. Combine values in a warp or block first, so few atomics reach the shared result.
• Float atomics and rounding. Every float atomicAdd rounds. Adding many small values into a large total loses precision; CUDA also doesn’t guarantee the order of atomic updates, so for mixed values the result can change between runs.
• Overflow. A large sum can overflow a 32-bit integer, even when every element is small. Sum into unsigned long long when the total can grow large.
• Missing barriers in a tree reduction. Every step of the shared-memory loop needs __syncthreads(), or a thread reads a partial sum that isn’t finished yet.

Practice

1. Make the histogram faster by giving each block its own 256-bin histogram in shared memory, merging it into the global one with one atomic per bin at the end. Use 160 blocks (4 per SM across the T4’s 40 SMs) that loop over the data.

// hist_shared.cu - each block counts into its own shared-memory histogram, then merges once
#include <cstdio>
#include <cuda_runtime.h>

__global__ void histo(const unsigned* d, unsigned* h, int n) {
  __shared__ unsigned local[256];
  local[threadIdx.x] = 0;  // launched with 256 threads: one bin each
  __syncthreads();
  for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += gridDim.x * blockDim.x)
    atomicAdd(&local[d[i] % 256], 1u);  // contention stays inside the block
  __syncthreads();
  atomicAdd(&h[threadIdx.x], local[threadIdx.x]);  // 256 global atomics per block
}

float run(const unsigned* d, unsigned* h, unsigned* out, int n) {
  cudaEvent_t a, b; cudaEventCreate(&a); cudaEventCreate(&b);
  cudaMemset(h, 0, 256 * 4);
  cudaEventRecord(a); histo<<<160, 256>>>(d, h, n); cudaEventRecord(b);
  cudaMemcpy(out, h, 256 * 4, cudaMemcpyDeviceToHost);
  float ms; cudaEventElapsedTime(&ms, a, b); return ms;
}

int main() {
  const int n = 1 << 24; unsigned *x = new unsigned[n], out[256], *d, *h, cpu7 = 0;
  for (int i = 0; i < n; i++) { x[i] = (i * 2654435761u) >> 8; cpu7 += x[i] % 256 == 7; }
  cudaMalloc(&d, n * 4); cudaMalloc(&h, 256 * 4);
  cudaMemcpy(d, x, n * 4, cudaMemcpyHostToDevice);
  run(d, h, out, n);  // warm-up
  float spread = run(d, h, out, n);
  printf("spread over 256 bins: %6.2f ms, bin 7 = %u, CPU says %u\n", spread, out[7], cpu7);
  cudaMemset(d, 0, n * 4);
  float same = run(d, h, out, n);
  printf("all in one bin:       %6.2f ms, bin 0 = %u\n", same, out[0]);
}

$ nvcc -o hist_shared hist_shared.cu && ./hist_shared
spread over 256 bins:   0.31 ms, bin 7 = 65536, CPU says 65536
all in one bin:         0.79 ms, bin 0 = 16777216

2. Sum a million copies of 0.1f with one float atomicAdd each, and compare the result with a double-precision sum on the CPU.

// float_sum.cu - a million float atomicAdds into one float: the rounding error piles up
#include <cstdio>
#include <cuda_runtime.h>

__global__ void sum(const float* x, float* total, int n) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i < n) atomicAdd(total, x[i]);
}

int main() {
  const int n = 1 << 20; float *x, *total;
  cudaMallocManaged(&x, n * 4); cudaMallocManaged(&total, 4);
  double cpu = 0;
  for (int i = 0; i < n; i++) { x[i] = 0.1f; cpu += x[i]; }
  for (int run = 1; run <= 3; run++) {
    *total = 0.0f;
    sum<<<n / 256, 256>>>(x, total, n);
    cudaDeviceSynchronize();
    printf("launch %d: %.4f\n", run, *total);
  }
  printf("double on the CPU: %.4f\n", cpu);
}

$ nvcc -o float_sum float_sum.cu && ./float_sum
launch 1: 105891.8438
launch 2: 105891.8438
launch 3: 105891.8438
double on the CPU: 104857.6016

3. Fix the float sum by reducing each block in shared memory first, as in reduce_block.cu.

// float_block.cu - the same float sum, reduced in each block before one atomicAdd per block
#include <cstdio>
#include <cuda_runtime.h>

__global__ void sum(const float* x, float* total, int n) {
  __shared__ float s[256];
  int t = threadIdx.x, i = blockIdx.x * 256 + t;
  s[t] = i < n ? x[i] : 0.0f;
  __syncthreads();
  for (int stride = 128; stride > 0; stride /= 2) {
    if (t < stride) s[t] += s[t + stride];
    __syncthreads();
  }
  if (t == 0) atomicAdd(total, s[0]);  // 4096 additions of 25.6 instead of a million of 0.1
}

int main() {
  const int n = 1 << 20; float *x, *total;
  cudaMallocManaged(&x, n * 4); cudaMallocManaged(&total, 4);
  double cpu = 0;
  for (int i = 0; i < n; i++) { x[i] = 0.1f; cpu += x[i]; }
  for (int run = 1; run <= 3; run++) {
    *total = 0.0f;
    sum<<<n / 256, 256>>>(x, total, n);
    cudaDeviceSynchronize();
    printf("launch %d: %.4f\n", run, *total);
  }
  printf("double on the CPU: %.4f\n", cpu);
}

$ nvcc -o float_block float_block.cu && ./float_block
launch 1: 104861.6406
launch 2: 104861.6406
launch 3: 104861.6406
double on the CPU: 104857.6016

4. CUDA has no atomicMax for floats. Build one from atomicCAS and check it against the CPU.

// cas_max.cu - a float atomicMax built from atomicCAS
#include <cstdio>
#include <cuda_runtime.h>

__device__ void atomic_max_float(float* addr, float v) {
  int* p = (int*)addr;
  int old = *p, assumed;
  do {
    assumed = old;
    if (__int_as_float(assumed) >= v) return;                // already big enough
    old = atomicCAS(p, assumed, __float_as_int(v));          // swap only if nobody changed it
  } while (old != assumed);                                  // someone did: retry with their value
}

__global__ void kmax(const float* x, float* m, int n) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i < n) atomic_max_float(m, x[i]);
}

int main() {
  const int n = 1 << 20; float *x, *m, cpu = -1e30f;
  cudaMallocManaged(&x, n * 4); cudaMallocManaged(&m, 4);
  for (int i = 0; i < n; i++) { x[i] = (1LL * i * 7919 % 1000003) / 1000.0f; cpu = x[i] > cpu ? x[i] : cpu; }
  *m = -1e30f;
  kmax<<<n / 256, 256>>>(x, m, n);
  cudaDeviceSynchronize();
  printf("GPU max = %.3f, CPU max = %.3f\n", *m, cpu);
}

$ nvcc -o cas_max cas_max.cu && ./cas_max
GPU max = 1000.002, CPU max = 1000.002