从CPU到GPU：手把手教你用CUDA实现高性能并行前缀和（附完整代码与性能对比）

从CPU到GPU手把手教你用CUDA实现高性能并行前缀和附完整代码与性能对比在数据处理和科学计算领域前缀和Prefix Sum是一个基础但极其重要的算法。它不仅是排序、流压缩等高级算法的基础构件更在图像处理、物理模拟等场景中扮演关键角色。传统CPU实现虽然直观易懂但在处理大规模数据时性能捉襟见肘。本文将带您从零开始通过CUDA实现四种不同版本的前缀和算法并深入分析每种实现的性能特点与优化技巧。1. 基础概念与CPU实现前缀和又称扫描Scan是指对一个输入序列[x₀, x₁,..., xₙ₋₁]计算输出序列[y₀, y₁,..., yₙ₋₁]其中每个yᵢ都是x₀到xᵢ的和。根据是否包含当前元素可分为**包含型inclusive和排除型exclusive**两种形式包含型yᵢ x₀ x₁ ... xᵢ 排除型yᵢ x₀ x₁ ... xᵢ₋₁CPU上的串行实现非常简单直接// 包含型前缀和 void cpu_inclusive_scan(float* input, float* output, int n) { output[0] input[0]; for (int i 1; i n; i) { output[i] output[i-1] input[i]; } } // 排除型前缀和 void cpu_exclusive_scan(float* input, float* output, int n) { output[0] 0; for (int i 1; i n; i) { output[i] output[i-1] input[i-1]; } }这种实现的时间复杂度为O(n)在小数据量时表现尚可。但当数据规模达到百万甚至千万级别时串行计算的局限性就暴露无遗。下表展示了不同数据规模下CPU实现的运行时间数据规模执行时间(ms)1K0.021M2.110M21.8100M218.5测试环境Intel i7-10700K 3.8GHz单线程执行2. GPU并行计算基础与CUDA实现2.1 CUDA编程模型简介CUDA是NVIDIA推出的通用并行计算架构其核心概念包括网格Grid最高层次的并行组织线程块Block共享内存和同步的基本单位线程Thread最小执行单元GPU的优势在于其大规模并行处理能力。以NVIDIA RTX 3090为例拥有10496个CUDA核心可以同时执行大量线程。2.2 朴素并行实现最直观的GPU实现方式是让每个线程独立计算自己的前缀和__global__ void naive_scan(float* input, float* output, int n) { int idx blockIdx.x * blockDim.x threadIdx.x; if (idx n) return; float sum 0; for (int i 0; i idx; i) { sum input[i]; } output[idx] sum; }这种实现虽然简单但存在严重问题每个线程重复计算前面的元素计算复杂度高达O(n²)全局内存访问模式不佳无法利用缓存线程间没有协作完全丧失了并行优势实测性能甚至比CPU串行实现还要慢这提醒我们简单的并行化不一定带来性能提升。3. 高效并行扫描算法3.1 Hillis-Steele算法Hillis-Steele算法是一种工作高效的并行扫描算法其核心思想是通过多轮迭代逐步构建前缀和for s from 1 to log2(n): for all k in parallel: if k s: x[k] x[k-s]CUDA实现需要考虑共享内存和线程同步__global__ void hillis_steele_scan(float* input, float* output, int n) { extern __shared__ float temp[]; int tid threadIdx.x; int idx blockIdx.x * blockDim.x tid; // 加载数据到共享内存 temp[tid] (idx n) ? input[idx] : 0; __syncthreads(); // 并行扫描 for (int s 1; s blockDim.x; s * 2) { if (tid s) { temp[tid] temp[tid - s]; } __syncthreads(); } // 写回结果 if (idx n) { output[idx] temp[tid]; } }该算法的时间复杂度为O(log n)但存在bank冲突问题。当多个线程同时访问共享内存中属于同一个bank的不同地址时会导致串行访问降低性能。3.2 双缓冲优化为解决同步和bank冲突问题可以采用双缓冲技术__global__ void hillis_steele_double_buffer(float* input, float* output, int n) { extern __shared__ float temp[]; float* buf1 temp; float* buf2 temp blockDim.x; int tid threadIdx.x; int idx blockIdx.x * blockDim.x tid; // 初始化双缓冲 buf1[tid] (idx n) ? input[idx] : 0; __syncthreads(); for (int s 1; s blockDim.x; s * 2) { // 交替使用缓冲区 if (tid s) { buf2[tid] buf1[tid] buf1[tid - s]; } else { buf2[tid] buf1[tid]; } __syncthreads(); // 交换缓冲区指针 float* tmp buf1; buf1 buf2; buf2 tmp; } if (idx n) { output[idx] buf1[tid]; } }3.3 Blelloch算法Blelloch算法采用不同的策略分为上扫reduce和下扫down-sweep两个阶段__global__ void blelloch_scan(float* input, float* output, int n) { extern __shared__ float temp[]; int tid threadIdx.x; int idx 2 * blockIdx.x * blockDim.x tid; // 加载两倍数据 temp[tid] (idx n) ? input[idx] : 0; temp[tid blockDim.x] (idx blockDim.x n) ? input[idx blockDim.x] : 0; __syncthreads(); // 上扫阶段 for (int s 1; s blockDim.x; s * 2) { int index 2 * s * (tid 1) - 1; if (index 2 * blockDim.x) { temp[index] temp[index - s]; } __syncthreads(); } // 清零最后一个元素排除型扫描 if (tid 0) { temp[2 * blockDim.x - 1] 0; } __syncthreads(); // 下扫阶段 for (int s blockDim.x; s 1; s / 2) { int index 2 * s * (tid 1) - 1; if (index 2 * blockDim.x) { float t temp[index]; temp[index] temp[index - s]; temp[index - s] t; } __syncthreads(); } // 写回结果 if (idx n) output[idx] temp[tid]; if (idx blockDim.x n) output[idx blockDim.x] temp[tid blockDim.x]; }3.4 消除Bank冲突的优化通过填充padding技术可以消除共享内存的bank冲突#define PADDING_SIZE 1 // 每32个元素填充1个 __global__ void blelloch_padded_scan(float* input, float* output, int n) { extern __shared__ float temp[]; int tid threadIdx.x; int idx 2 * blockIdx.x * blockDim.x tid; // 计算带padding的索引 auto padded_idx [](int i) { return i i / 32; }; // 加载数据到带padding的共享内存 temp[padded_idx(tid)] (idx n) ? input[idx] : 0; temp[padded_idx(tid blockDim.x)] (idx blockDim.x n) ? input[idx blockDim.x] : 0; __syncthreads(); // 上扫阶段与标准Blelloch相同但使用padded索引 // ... // 下扫阶段与标准Blelloch相同但使用padded索引 // ... // 写回结果 if (idx n) output[idx] temp[padded_idx(tid)]; if (idx blockDim.x n) output[idx blockDim.x] temp[padded_idx(tid blockDim.x)]; }4. 任意长度数据扫描上述算法都假设数据可以完全放入一个线程块处理。对于更大规模数据需要采用分而治之的策略将数据划分为多个块对每个块独立扫描收集每个块的总和并扫描将块总和扫描结果加到原始块上实现分为三个核心函数// 块内扫描 __global__ void block_scan(float* input, float* output, float* block_sums, int n) { // 实现块内扫描并保存块总和 // ... } // 块总和扫描 void scan_block_sums(float* block_sums, int num_blocks) { // 递归调用扫描算法 // ... } // 分散相加 __global__ void add_block_sums(float* output, float* block_sums, int n) { // 将块总和加到对应块上 // ... }5. 性能对比与优化建议我们对五种实现进行了性能测试NVIDIA RTX 3090数据规模10M算法版本执行时间(ms)加速比(相对CPU)CPU串行21.81xGPU朴素45.20.48xHillis-Steele2.110.4xBlelloch1.812.1xBlellochpadding1.316.8x任意长度Blelloch1.514.5x优化建议选择合适的算法小数据量用Hillis-Steele大数据量用Blelloch最大化并行度确保每个SM有足够的线程块优化内存访问使用共享内存减少全局内存访问合并全局内存访问避免bank冲突隐藏延迟通过足够的线程块和线程实现指令级并行6. 完整代码实现以下是经过优化的完整Blelloch扫描实现#include stdio.h #include stdlib.h #include cuda_runtime.h #define BLOCK_SIZE 256 #define PADDING(idx) ((idx) (idx)/32) void check_cuda_error(cudaError_t err, const char* msg) { if (err ! cudaSuccess) { fprintf(stderr, CUDA Error: %s: %s\n, msg, cudaGetErrorString(err)); exit(EXIT_FAILURE); } } __global__ void blelloch_scan_kernel(float* input, float* output, float* block_sums, int n) { extern __shared__ float temp[]; int tid threadIdx.x; int bid blockIdx.x; int idx 2 * bid * blockDim.x tid; // 加载数据到共享内存带padding temp[PADDING(tid)] (idx n) ? input[idx] : 0; temp[PADDING(tid blockDim.x)] (idx blockDim.x n) ? input[idx blockDim.x] : 0; __syncthreads(); // 上扫阶段 for (int s 1; s blockDim.x; s * 2) { int index 2 * s * (tid 1) - 1; if (index 2 * blockDim.x) { temp[PADDING(index)] temp[PADDING(index - s)]; } __syncthreads(); } // 排除型扫描清零最后一个元素 if (tid 0) { temp[PADDING(2 * blockDim.x - 1)] 0; } __syncthreads(); // 下扫阶段 for (int s blockDim.x; s 1; s / 2) { int index 2 * s * (tid 1) - 1; if (index 2 * blockDim.x) { float t temp[PADDING(index)]; temp[PADDING(index)] temp[PADDING(index - s)]; temp[PADDING(index - s)] t; } __syncthreads(); } // 写回结果 if (idx n) output[idx] temp[PADDING(tid)]; if (idx blockDim.x n) output[idx blockDim.x] temp[PADDING(tid blockDim.x)]; // 保存块总和 if (block_sums ! NULL tid 0) { block_sums[bid] temp[PADDING(2 * blockDim.x - 1)] ((2 * bid * blockDim.x 2 * blockDim.x - 1 n) ? input[2 * bid * blockDim.x 2 * blockDim.x - 1] : 0); } } void scan(float* d_input, float* d_output, int n) { float *d_block_sums NULL; int grid_size (n 2 * BLOCK_SIZE - 1) / (2 * BLOCK_SIZE); if (grid_size 1) { check_cuda_error(cudaMalloc(d_block_sums, grid_size * sizeof(float)), Allocate block sums); // 第一遍扫描计算块内前缀和和块总和 blelloch_scan_kernelgrid_size, BLOCK_SIZE, (2 * BLOCK_SIZE (2 * BLOCK_SIZE)/32) * sizeof(float)( d_input, d_output, d_block_sums, n); check_cuda_error(cudaGetLastError(), First scan kernel); // 递归扫描块总和 scan(d_block_sums, d_block_sums, grid_size); // 第二遍扫描将块总和加到对应块上 dim3 add_grid((n BLOCK_SIZE - 1) / BLOCK_SIZE); add_block_sumsadd_grid, BLOCK_SIZE(d_output, d_block_sums, n); check_cuda_error(cudaGetLastError(), Add block sums kernel); cudaFree(d_block_sums); } else { // 单个块可以直接处理 blelloch_scan_kernel1, BLOCK_SIZE, (2 * BLOCK_SIZE (2 * BLOCK_SIZE)/32) * sizeof(float)( d_input, d_output, NULL, n); check_cuda_error(cudaGetLastError(), Single block scan kernel); } } __global__ void add_block_sums(float* output, float* block_sums, int n) { int tid threadIdx.x; int bid blockIdx.x; int idx bid * blockDim.x tid; if (idx n) { int block_num idx / (2 * BLOCK_SIZE); if (block_num 0) { output[idx] block_sums[block_num - 1]; } } }