Reduction & Scan

Classic optimization steps (Mark Harris's "Optimizing Parallel Reduction"): interleaved → sequential addressing (avoid divergence and bank conflicts), first add during load, unroll the last warp, complete unrolling with templates, multiple elements per thre…

14.1 Parallel reduction (sum)

Classic optimization steps (Mark Harris's "Optimizing Parallel Reduction"): interleaved → sequential addressing (avoid divergence and bank conflicts), first add during load, unroll the last warp, complete unrolling with templates, multiple elements per thread, warp shuffles. Use CUB cub::DeviceReduce::Sum in production.

// Block-level reduction: shared memory tree + warp shuffles; one atomic per block
__global__ void reduceSum(const float* in, float* out, int n) {
    __shared__ float warpSums[32];               // up to 1024 threads → 32 warps
    float sum = 0.0f;
    for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += blockDim.x * gridDim.x)
        sum += in[i];                             // grid-stride accumulation (sequential part)

    sum = warpReduceSum(sum);                     // reduce within each warp
    int lane = threadIdx.x & 31, wid = threadIdx.x >> 5;
    if (lane == 0) warpSums[wid] = sum;
    __syncthreads();

    if (wid == 0) {                               // first warp reduces the warp sums
        sum = (threadIdx.x < blockDim.x / 32) ? warpSums[lane] : 0.0f;
        sum = warpReduceSum(sum);
        if (lane == 0) atomicAdd(out, sum);
    }
}

14.2 Prefix sum (scan)

Inclusive scan: [3, 1, 7, 0, 4] → [3, 4, 11, 11, 15]; exclusive: [0, 3, 4, 11, 11].

Algorithms: Hillis–Steele (work-inefficient, O(n log n)), Blelloch (work-efficient up-sweep/down-sweep), decoupled look-back single-pass scan (CUB). Applications: stream compaction, radix sort, sparse matrix ops, allocation, histogram equalization.

// Warp-level inclusive scan with shuffles
__inline__ __device__ int warpInclusiveScan(int val) {
    int lane = threadIdx.x & 31;
    for (int offset = 1; offset < 32; offset <<= 1) {
        int n = __shfl_up_sync(0xffffffff, val, offset);
        if (lane >= offset) val += n;
    }
    return val;
}

// Block scan: scan each warp, scan the warp totals, add offsets
__global__ void blockScan(const int* in, int* out, int n) {
    __shared__ int warpTotals[32];
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    int v = (i < n) ? in[i] : 0;
    int lane = threadIdx.x & 31, wid = threadIdx.x >> 5;
    v = warpInclusiveScan(v);
    if (lane == 31) warpTotals[wid] = v;
    __syncthreads();
    if (wid == 0) {
        int t = (lane < (blockDim.x + 31) / 32) ? warpTotals[lane] : 0;
        warpTotals[lane] = warpInclusiveScan(t);
    }
    __syncthreads();
    if (wid > 0) v += warpTotals[wid - 1];
    if (i < n) out[i] = v;
    // Multi-block scans: scan block sums in a second pass, then add (or decoupled look-back as in CUB)
}