Loading...
Vector dot product is a fundamental operation in linear algebra, appearing in everything from neural networks to physics simulations. Computing the sum of element-wise products seems trivial, but achieving high performance requires sophisticated reduction techniques that exploit GPU hardware efficiently. A naive reduction can be 100x slower than optimized code. The challenge is combining results from thousands of threads while avoiding serialization bottlenecks. Modern GPUs provide warp-level primitives like __shfl_down_sync that enable lock-free reduction within warps, eliminating shared memory overhead. The optimization journey from naive atomic additions to hierarchical warp shuffles teaches critical parallel programming patterns: divergence avoidance, bank conflict elimination, and leveraging hardware primitives. These techniques apply to all reduction operations: sum, max, min, and custom associative operators.
Use __shfl_down_sync to perform reduction within warps without shared memory. This achieves maximum bandwidth and zero memory overhead for warp-level reductions.
__inline__ __device__
float warp_reduce_sum(float val) {
unsigned mask = 0xffffffff;
for (int offset = 16; offset > 0; offset /= 2) {
val += __shfl_down_sync(mask, val, offset);
}
return val;
}
__global__ void dot_product_shuffle(float* a, float* b, float* out, int n) {
int tid = blockIdx.x * blockDim.x + threadIdx.x;
int lane = threadIdx.x % 32;
float sum = 0.0f;
for (int i = tid; i < n; i += blockDim.x * gridDim.x) {
sum += a[i] * b[i];
}
// Reduce within warp
sum = warp_reduce_sum(sum);
// First thread in warp writes to output
if (lane == 0) {
atomicAdd(out, sum);
}
}Use shared memory for block-level reduction with sequential addressing to avoid divergence. Combine with warp shuffles for hierarchical reduction.
__global__ void dot_product_shared(float* a, float* b, float* out, int n) {
__shared__ float sdata[256];
int tid = threadIdx.x;
int idx = blockIdx.x * blockDim.x + threadIdx.x;
// Grid-stride loop for large vectors
float sum = 0.0f;
for (int i = idx; i < n; i += blockDim.x * gridDim.x) {
sum += a[i] * b[i];
}
sdata[tid] = sum;
__syncthreads();
// Sequential addressing (no divergence)
for (int s = blockDim.x / 2; s > 32; s >>= 1) {
if (tid < s) {
sdata[tid] += sdata[tid + s];
}
__syncthreads();
}
// Warp shuffle for final reduction
if (tid < 32) {
float warp_sum = sdata[tid];
warp_sum = warp_reduce_sum(warp_sum);
if (tid == 0) atomicAdd(out, warp_sum);
}
}First phase reduces to per-block values in global memory, second phase reduces blocks on CPU or with small kernel. Avoids atomic contention.
__global__ void dot_product_phase1(float* a, float* b,
float* block_sums, int n) {
__shared__ float sdata[256];
int tid = threadIdx.x;
int idx = blockIdx.x * blockDim.x + threadIdx.x;
float sum = 0.0f;
for (int i = idx; i < n; i += blockDim.x * gridDim.x) {
sum += a[i] * b[i];
}
// Block-level reduction
sdata[tid] = sum;
__syncthreads();
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) block_sums[blockIdx.x] = sdata[0];
}
__global__ void dot_product_phase2(float* block_sums,
float* result, int num_blocks) {
// Small kernel reduces block sums (single block)
__shared__ float sdata[256];
int tid = threadIdx.x;
sdata[tid] = (tid < num_blocks) ? block_sums[tid] : 0.0f;
__syncthreads();
for (int s = blockDim.x / 2; s > 0; s >>= 1) {
if (tid < s) sdata[tid] += sdata[tid + s];
__syncthreads();
}
if (tid == 0) *result = sdata[0];
}Use float4/double2 to load multiple elements per transaction, improving memory bandwidth utilization and reducing instruction count.
__global__ void dot_product_vectorized(float4* a, float4* b,
float* out, int n) {
int tid = blockIdx.x * blockDim.x + threadIdx.x;
float sum = 0.0f;
for (int i = tid; i < n / 4; i += blockDim.x * gridDim.x) {
float4 va = a[i];
float4 vb = b[i];
// Unrolled dot product of 4-element vectors
sum += va.x * vb.x;
sum += va.y * vb.y;
sum += va.z * vb.z;
sum += va.w * vb.w;
}
// Reduction using warp shuffle
sum = warp_reduce_sum(sum);
if ((threadIdx.x % 32) == 0) {
atomicAdd(out, sum);
}
}Naive atomic reduction causes severe contention. Thousands of threads serialize on a single memory location, achieving <1% of peak performance.
__global__ void dot_product_naive(float* a, float* b,
float* result, int n) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < n) {
float product = a[idx] * b[idx];
// Every thread atomically adds - massive contention!
atomicAdd(result, product);
}
}
void compute_dot_naive(float* d_a, float* d_b, float* d_result, int n) {
// Initialize result
cudaMemset(d_result, 0, sizeof(float));
int threads = 256;
int blocks = (n + threads - 1) / threads;
// All threads contend on single atomic
dot_product_naive<<<blocks, threads>>>(d_a, d_b, d_result, n);
cudaDeviceSynchronize();
}
// Performance: ~5 GFLOPS for 1M elements
// Problem: Serialized by atomic contentionOptimized reduction uses warp shuffles, vectorized loads, and two-phase reduction to eliminate atomic contention and maximize throughput.
__inline__ __device__ float warp_reduce_sum(float val) {
unsigned mask = 0xffffffff;
#pragma unroll
for (int offset = 16; offset > 0; offset >>= 1) {
val += __shfl_down_sync(mask, val, offset);
}
return val;
}
__global__ void dot_product_optimized(float4* a, float4* b,
float* block_results, int n) {
__shared__ float warp_sums[8]; // 256 threads = 8 warps
int tid = threadIdx.x;
int idx = blockIdx.x * blockDim.x + threadIdx.x;
int lane = tid % 32;
int warp_id = tid / 32;
// Grid-stride loop with vectorized access
float sum = 0.0f;
for (int i = idx; i < n / 4; i += blockDim.x * gridDim.x) {
float4 va = a[i];
float4 vb = b[i];
sum += va.x * vb.x + va.y * vb.y + va.z * vb.z + va.w * vb.w;
}
// Warp-level reduction
sum = warp_reduce_sum(sum);
// First thread in each warp writes to shared memory
if (lane == 0) {
warp_sums[warp_id] = sum;
}
__syncthreads();
// Final reduction by first warp
if (warp_id == 0) {
float warp_sum = (lane < 8) ? warp_sums[lane] : 0.0f;
warp_sum = warp_reduce_sum(warp_sum);
if (lane == 0) {
block_results[blockIdx.x] = warp_sum;
}
}
}
__global__ void reduce_blocks(float* block_results, float* final_result,
int num_blocks) {
__shared__ float sdata[256];
int tid = threadIdx.x;
float sum = 0.0f;
for (int i = tid; i < num_blocks; i += blockDim.x) {
sum += block_results[i];
}
sdata[tid] = sum;
__syncthreads();
// Reduce within block
for (int s = blockDim.x / 2; s > 32; s >>= 1) {
if (tid < s) sdata[tid] += sdata[tid + s];
__syncthreads();
}
if (tid < 32) {
float val = sdata[tid];
val = warp_reduce_sum(val);
if (tid == 0) *final_result = val;
}
}
void compute_dot_optimized(float* d_a, float* d_b, float* d_result, int n) {
int threads = 256;
int blocks = min(1024, (n / 4 + threads - 1) / threads);
float* d_block_results;
cudaMalloc(&d_block_results, blocks * sizeof(float));
// Phase 1: Reduce to per-block sums
dot_product_optimized<<<blocks, threads>>>((float4*)d_a, (float4*)d_b,
d_block_results, n);
// Phase 2: Reduce block sums
reduce_blocks<<<1, 256>>>(d_block_results, d_result, blocks);
cudaFree(d_block_results);
}
// Performance: 500+ GFLOPS for 1M elements
// 100x faster than naive| Metric | Naive | Optimized | Improvement |
|---|---|---|---|
| Throughput (1M elements) | 5.2 GFLOPS | 520 GFLOPS | 100x faster |
| Memory bandwidth utilization | 8% (atomic bottleneck) | 85% (coalesced) | 10.6x higher |
| Latency (1M elements) | 380 μs | 3.8 μs | 100x faster |
| vs cuBLAS cublasSdot | 1.2% of cuBLAS | 95% of cuBLAS | 79x closer |
Warp shuffles (__shfl_down_sync) exchange data between threads in a warp without memory accesses. This is 5-10x faster than shared memory (register vs memory latency) and requires zero memory. Shared memory is still needed for reductions larger than warp size (32 threads).
Sequential addressing means active threads have consecutive indices. In reduction, instead of having threads 0,2,4,... active (interleaved), use threads 0,1,2,... (sequential). This avoids warp divergence where half the warp is idle, improving efficiency from 50% to 100%.
Use two-phase for accuracy and performance: it avoids non-deterministic atomic ordering and contention. Use atomic only when you need the simplicity and have few blocks (<32) so contention is minimal. For production code, always use two-phase reduction.
Yes! The same patterns apply to any associative operator. Replace += with max(), min(), or custom operators. Warp shuffles work with any data type that fits in 32 bits. For larger types (double, int64), use __shfl_down_sync with appropriate casting or multiple shuffles.
Element-wise operation before reduction
Core pattern used in dot product
Multiple dot products in parallel
Ready to optimize your CUDA code? Download RightNow AI and get real-time performance analysis for your kernels.