CUDA Atomics & Parallel Reductions

Atomics: Theory (What They Guarantee)

Atomics perform read-modify-write operations as a single, indivisible step. In CUDA, atomics can help implement:

• counters (increment)
• histograms (add to a bin)
• reductions (sum/min/max)
• lock-free algorithms (advanced)

Trade-off:

• Under contention (many threads updating the same address), atomics serialize and reduce throughput.

Code Example 1 — Atomic Counter Histogram

#include <cuda_runtime.h>
#include <cstdio>

__global__ void histogram_atomic(const unsigned* data, unsigned* hist, int N, unsigned bins) {
  int idx = blockIdx.x * blockDim.x + threadIdx.x;
  if (idx >= N) return;

  unsigned value = data[idx];
  unsigned bin = value % bins;

  // Atomic add ensures correctness when many threads update the same bin
  atomicAdd(&hist[bin], 1u);
}

int main() {
  const int N = 1 << 20;
  const unsigned bins = 256;

  unsigned *h_data = (unsigned*)malloc(N * sizeof(unsigned));
  for (int i = 0; i < N; i++) h_data[i] = i * 2654435761u;

  unsigned *d_data, *d_hist;
  cudaMalloc(&d_data, N * sizeof(unsigned));
  cudaMalloc(&d_hist, bins * sizeof(unsigned));

  cudaMemcpy(d_data, h_data, N * sizeof(unsigned), cudaMemcpyHostToDevice);
  cudaMemset(d_hist, 0, bins * sizeof(unsigned));

  int threads = 256;
  int blocks  = (N + threads - 1) / threads;

  histogram_atomic<<<blocks, threads>>>(d_data, d_hist, N, bins);

  unsigned *h_hist = (unsigned*)malloc(bins * sizeof(unsigned));
  cudaMemcpy(h_hist, d_hist, bins * sizeof(unsigned), cudaMemcpyDeviceToHost);

  printf("hist[0]=%u  hist[1]=%u  hist[2]=%u\n", h_hist[0], h_hist[1], h_hist[2]);

  cudaFree(d_data);
  cudaFree(d_hist);
  free(h_data);
  free(h_hist);
  return 0;
}

Code Example 2 — Two-Stage Parallel Reduction (Block + Global)

Efficient reductions typically:

1. Reduce within each block using shared memory
2. Write partial sums to global memory
3. Reduce partial sums (recursively or with another kernel)

Kernel: block-level reduction

#include <cuda_runtime.h>

__global__ void reduce_block_sum(const float* x, float* partial, int N) {
  extern __shared__ float sdata[]; // size = blockDim.x

  int tid = threadIdx.x;
  int i = blockIdx.x * blockDim.x * 2 + threadIdx.x; // load two elements per thread

  float sum = 0.0f;
  if (i < N) sum += x[i];
  if (i + blockDim.x < N) sum += x[i + blockDim.x];

  sdata[tid] = sum;
  __syncthreads();

  // Reduce in shared memory
  for (int s = blockDim.x / 2; s > 0; s >>= 1) {
    if (tid < s) sdata[tid] += sdata[tid + s];
    __syncthreads();
  }

  // Write block result
  if (tid == 0) partial[blockIdx.x] = sdata[0];
}

Host-side sketch for 2-stage reduction

#include <cuda_runtime.h>
#include <vector>
#include <cstdio>

void reduce_two_stage(const float* d_x, int N) {
  int threads = 256;
  int blocks  = (N + threads * 2 - 1) / (threads * 2);

  float* d_partial;
  cudaMalloc(&d_partial, blocks * sizeof(float));

  // shared memory size per block = threads * sizeof(float)
  size_t shmem = threads * sizeof(float);

  reduce_block_sum<<<blocks, threads, shmem>>>(d_x, d_partial, N);

  // Second stage: reduce partial array.
  // For illustration, you can reduce again (recursively) until size == 1.
  // In production, implement a loop until you get a single value.
  printf("Partial results computed for %d blocks\n", blocks);

  cudaFree(d_partial);
}

Common Gotchas

• High contention atomics: if many threads hit the same address, atomics can dominate runtime.
• Incorrect reduction indexing: reductions must handle N not being a power of two.
• Unnecessary global atomics: prefer block-level reduction then one final reduction stage.
• Race conditions in shared memory: always place __syncthreads() where needed.

Quick Checklist

• Use atomics for correct single-location updates; benchmark under expected contention.
• Use hierarchical reductions: block → partial → final.
• Avoid global serialization (lots of global atomics).
• Validate correctness on small arrays before scaling up.