在前一课(第10课)中,我们学习了Reduction归约操作——将多个数值聚合为一个。今天,我们将探讨另一个同样关键的并行算法:Prefix Sum(前缀和)。该算法堪称GPU上诸多复杂算子的基础,许多常见操作底层都依赖它:例如stream compaction(流压缩)、radix sort(基数排序)、histogram(直方图生成)、稀疏矩阵CSR格式构建、token offset计算、batch packing(批次打包)、并行过滤、并行分配输出位置等。这些看似无关的任务,实际上都需要高效的前缀和算法支撑。因此,深入理解本章内容非常值得。
1. 本课目标
我们将围绕六个相互关联的核心问题展开,每一环都至关重要:
- 理解inclusive scan(包含扫描)与exclusive scan(排他扫描)的定义及区别,到底有何不同?
- 为什么scan这一看似简单的操作,无法像vector add那样直接并行化?
- Blelloch Scan中的up-sweep(向上规约)和down-sweep(向下传播)思想,如何巧妙绕开致命的依赖链?
- 如何利用shared memory高效实现block-level scan(块级扫描)。
- 当数组规模超出单个block处理能力时,如何通过多轮block scan完成大数组的前缀和计算。
- 最后,用实际性能数据让CPU与GPU正面交锋,验证并行优势。
2. 什么是 Prefix Sum?
先从一个具体例子入手。假设我们有以下输入数组:
input = [3, 1, 7, 0, 4, 1, 6, 3]
那么,前缀和是什么?它有两种最常见的变体。
2.1 Inclusive Scan
“包含自身”的前缀和,指每个位置的输出都是从开头累加至当前位置(包括自身)。计算结果为:
inclusive = [3, 4, 11, 11, 15, 16, 22, 25]
其数学定义如下:
out[i] = input[0] + input[1] + ... + input[i]
2.2 Exclusive Scan
“不包含自身”的前缀和(又称排他性扫描)。每个位置的输出是从开头累加到前一个元素,不含自身:
exclusive = [0, 3, 4, 11, 11, 15, 16, 22]
数学定义为:
out[i] = input[0] + input[1] + ... + input[i-1]
在CUDA实际应用中,exclusive scan更为常见。为什么?举个例子你就明白了。假设有一个flag数组,用于标记哪些元素是有效的:
flag = [1, 0, 1, 1, 0, 1]
对该flag数组执行exclusive scan后,得到:
pos = [0, 1, 1, 2, 3, 3]
可以看到,这个pos数组直接告诉了我们:每个有效元素应被写入输出数组的哪个位置。这是一种极其优雅的“地址计算”方式,是stream compaction等算法的核心。
3. 为什么 Prefix Sum 难以直接并行?
用CPU串行实现前缀和,代码写起来非常简单:
out[0] = 0;
for (int i = 1; i < n; ++i) {
out[i] = out[i - 1] + input[i - 1];
}
但这背后的问题很明显:每一个out[i]都直接依赖于out[i-1],而out[i-1]又依赖于out[i-2],依此类推,形成一道锁链般的依赖关系。GPU上那种“每个线程独立计算一个元素”的简单并行策略在此完全失效。要想在GPU上高效实现前缀和,必须从计算结构上彻底重构。
4. Blelloch Scan 核心思想
Blelloch Scan正是破解该难题的利器。其核心思想分为两个阶段:
- Up-sweep(向上规约):构建一棵部分和树。与我们之前学习的Reduction类似,两两相加,再四四相加,直至得到整个数组的总和。
- Down-sweep(向下传播):将总和置为0,然后沿着这棵“和树”将前缀和逐步向下分发,最终得到exclusive scan的结果。
仍然以刚才的8个元素为例演示这一过程:
输入: [3, 1, 7, 0, 4, 1, 6, 3]
Up-sweep阶段,先两两求和,再四个一组求和,最后得到总和25。接着Down-sweep阶段,将树根(总和位置)置0,然后一层层向下传播,把前面累加的部分和“推”下去。
最终得到exclusive scan的结果:
exclusive scan = [0, 3, 4, 11, 11, 15, 16, 22]
该算法的精妙之处在于:串行复杂度为O(n),而并行scan只需O(log n)层同步。在一个block内,我们可以利用shared memory高效实现这个树形结构。

5. 实验内容
光说不练假把式。本次实验需要动手实现以下五个部分:
- 一个CPU版本的exclusive scan,作为性能对比的基准。
- 一个基于GPU block-level的Blelloch scan。
- 一个能够处理大数组的多轮GPU scan方案。
- 对GPU结果进行校验,确保与CPU结果一致。
- 最后,进行令人期待的性能对比分析。

6. CUDA C++ 代码
直接上代码。整个实现包含了上述所有核心逻辑:
#include
#include
#include
#include
#include
#include
#include
#define CUDA_CHECK(call) \
do { \
cudaError_t err = call; \
if (err != cudaSuccess) { \
std::cerr << "CUDA Error: " << cudaGetErrorString(err) \
<< " at " << __FILE__ << ":" << __LINE__ << std::endl; \
std::exit(EXIT_FAILURE); \
} \
} while (0)
/*
* CPU exclusive scan:
*
* input : [x0, x1, x2, x3, ...]
* output : [0, x0, x0+x1, x0+x1+x2, ...]
*/
void exclusive_scan_cpu(const std::vector& input,
std::vector& output) {
if (input.empty()) return;
output[0] = 0.0f;
for (size_t i = 1; i < input.size(); ++i) {
output[i] = output[i - 1] + input[i - 1];
}
}
/*
* block-level exclusive scan.
*
* 每个 block 处理 2 * blockDim.x 个元素。
*
* 例如 blockDim.x = 256:
* 一个 block 处理 512 个 float。
*
* block_sums[blockIdx.x] 保存当前 block 的总和,
* 后续要对 block_sums 再做 scan,得到每个 block 的 offset。
*/
__global__ void block_exclusive_scan_kernel(const float* input,
float* output,
float* block_sums,
size_t n) {
extern __shared__ float temp[];
unsigned int tid = threadIdx.x;
size_t block_start = static_cast(blockIdx.x) * blockDim.x * 2;
size_t idx1 = block_start + tid;
size_t idx2 = block_start + tid + blockDim.x;
/*
* 1. 从 global memory 读入 shared memory。
*
* 每个线程最多读两个元素。
* 越界位置补 0。
*/
temp[tid] = (idx1 < n) ? input[idx1] : 0.0f;
temp[tid + blockDim.x] = (idx2 < n) ? input[idx2] : 0.0f;
__syncthreads();
/*
* 2. Up-sweep 阶段。
*
* 构造一棵求和树。
*
* 例如 512 个元素:
* 512 -> 256 -> 128 -> ... -> 1
*/
unsigned int offset = 1;
unsigned int total_elems = blockDim.x * 2;
for (unsigned int d = total_elems >> 1; d > 0; d >>= 1) {
if (tid < d) {
unsigned int ai = offset * (2 * tid + 1) - 1;
unsigned int bi = offset * (2 * tid + 2) - 1;
temp[bi] += temp[ai];
}
offset <<= 1;
__syncthreads();
}
/*
* 3. 保存 block 总和,并把根节点置 0。
*
* 置 0 是 exclusive scan 的关键。
*/
if (tid == 0) {
if (block_sums != nullptr) {
block_sums[blockIdx.x] = temp[total_elems - 1];
}
temp[total_elems - 1] = 0.0f;
}
__syncthreads();
/*
* 4. Down-sweep 阶段。
*
* 把前缀和从树根向下传播。
*/
for (unsigned int d = 1; d < total_elems; d <<= 1) {
offset >>= 1;
if (tid < d) {
unsigned int ai = offset * (2 * tid + 1) - 1;
unsigned int bi = offset * (2 * tid + 2) - 1;
float t = temp[ai];
temp[ai] = temp[bi];
temp[bi] += t;
}
__syncthreads();
}
/*
* 5. 写回 global memory。
*/
if (idx1 < n) {
output[idx1] = temp[tid];
}
if (idx2 < n) {
output[idx2] = temp[tid + blockDim.x];
}
}
/*
* 给每个 block 的局部 scan 结果加上 block offset。
*
* block_offsets 是对 block_sums 做 exclusive scan 后得到的结果。
*/
__global__ void add_block_offsets_kernel(float* data,
const float* block_offsets,
size_t n,
int elems_per_block) {
size_t idx = static_cast(blockIdx.x) * blockDim.x + threadIdx.x;
if (idx < n) {
size_t block_id = idx / elems_per_block;
data[idx] += block_offsets[block_id];
}
}
/*
* 递归式 GPU exclusive scan。
*
* 思路:
* 1. 每个 block 做局部 scan,得到 output 和 block_sums
* 2. 对 block_sums 再做 exclusive scan,得到 block_offsets
* 3. output 每个元素加上自己 block 的 offset
*/
void exclusive_scan_gpu_recursive(const float* d_input,
float* d_output,
size_t n,
int block_size) {
if (n == 0) return;
int elems_per_block = block_size * 2;
size_t num_blocks = (n + elems_per_block - 1) / elems_per_block;
float* d_block_sums = nullptr;
CUDA_CHECK(cudaMalloc(&d_block_sums, num_blocks * sizeof(float)));
size_t shared_bytes = static_cast(elems_per_block) * sizeof(float);
block_exclusive_scan_kernel<<(num_blocks),
block_size,
shared_bytes>>>(d_input,
d_output,
d_block_sums,
n);
CUDA_CHECK(cudaGetLastError());
if (num_blocks > 1) {
float* d_block_offsets = nullptr;
CUDA_CHECK(cudaMalloc(&d_block_offsets, num_blocks * sizeof(float)));
/*
* 递归地对 block_sums 做 exclusive scan。
*/
exclusive_scan_gpu_recursive(d_block_sums,
d_block_offsets,
num_blocks,
block_size);
/*
* 把 block offset 加回每个元素。
*/
int add_threads = 256;
int add_grid = static_cast((n + add_threads - 1) / add_threads);
add_block_offsets_kernel<<>>(d_output,
d_block_offsets,
n,
elems_per_block);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaFree(d_block_offsets));
}
CUDA_CHECK(cudaFree(d_block_sums));
}
/*
* 计时 GPU scan。
*
* 注意:
* 这里用 cudaEvent 测 GPU kernel 路径时间。
* 临时 cudaMalloc/cudaFree 不应作为核心性能结论。
* 工程中应该预分配临时 buffer。
*/
float time_gpu_scan(const float* d_input,
float* d_output,
size_t n,
int block_size,
int repeat) {
/*
* warmup
*/
exclusive_scan_gpu_recursive(d_input, d_output, n, block_size);
CUDA_CHECK(cudaDeviceSynchronize());
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
float total_ms = 0.0f;
for (int r = 0; r < repeat; ++r) {
CUDA_CHECK(cudaEventRecord(start));
exclusive_scan_gpu_recursive(d_input, d_output, n, block_size);
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop));
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
total_ms += ms;
}
CUDA_CHECK(cudaEventDestroy(start));
CUDA_CHECK(cudaEventDestroy(stop));
return total_ms / repeat;
}
float max_abs_diff(const std::vector& a, const std::vector& b) {
float max_diff = 0.0f;
for (size_t i = 0; i < a.size(); ++i) {
float diff = std::fabs(a[i] - b[i]);
if (diff > max_diff) {
max_diff = diff;
}
}
return max_diff;
}
int main(int argc, char** argv) {
size_t n = 1ULL << 22; // 4,194,304 floats,约 16 MB
int block_size = 256;
int repeat = 10;
if (argc >= 2) {
n = static_cast(std::atoll(argv[1]));
}
if (argc >= 3) {
block_size = std::atoi(argv[2]);
}
if (argc >= 4) {
repeat = std::atoi(argv[3]);
}
if (block_size <= 0 || block_size > 1024) {
std::cerr << "block_size must be in (0, 1024]." << std::endl;
return 1;
}
/*
* 为了让 float 结果容易校验,默认全部初始化为 1.0f。
* 当 n <= 2^24 时,前缀和整数部分可以被 float 精确表示。
*/
std::vector h_input(n, 1.0f);
std::vector h_cpu(n, 0.0f);
std::vector h_gpu(n, 0.0f);
size_t bytes = n * sizeof(float);
std::cout << "CUDA Lesson 11: Parallel Prefix Sum / Scan" << std::endl;
std::cout << "Elements : " << n << std::endl;
std::cout << "Data size : " << bytes / 1024.0 / 1024.0 << " MB" << std::endl;
std::cout << "Block size: " << block_size << std::endl;
std::cout << "Repeat : " << repeat << std::endl;
/*
* CPU timing。
*/
auto cpu_start = std::chrono::high_resolution_clock::now();
exclusive_scan_cpu(h_input, h_cpu);
auto cpu_end = std::chrono::high_resolution_clock::now();
double cpu_ms = std::chrono::duration(cpu_end - cpu_start).count();
float* d_input = nullptr;
float* d_output = nullptr;
CUDA_CHECK(cudaMalloc(&d_input, bytes));
CUDA_CHECK(cudaMalloc(&d_output, bytes));
CUDA_CHECK(cudaMemcpy(d_input, h_input.data(), bytes, cudaMemcpyHostToDevice));
float gpu_ms = time_gpu_scan(d_input, d_output, n, block_size, repeat);
CUDA_CHECK(cudaMemcpy(h_gpu.data(), d_output, bytes, cudaMemcpyDeviceToHost));
float diff = max_abs_diff(h_cpu, h_gpu);
std::cout << std::fixed << std::setprecision(4);
std::cout << "[Timing]" << std::endl;
std::cout << "CPU exclusive scan time : " << cpu_ms << " ms" << std::endl;
std::cout << "GPU exclusive scan time : " << gpu_ms << " ms" << std::endl;
std::cout << "Speedup : " << cpu_ms / gpu_ms << "x" << std::endl;
std::cout << "[Check]" << std::endl;
std::cout << "Max abs diff: " << diff << std::endl;
std::cout << "Result : " << (diff < 1e-3f ? "PASS" : "CHECK") << std::endl;
std::cout << "[Sample]" << std::endl;
std::cout << "input[0..7] : ";
for (int i = 0; i < 8 && i < static_cast(n); ++i) {
std::cout << h_input[i] << " ";
}
std::cout << std::endl;
std::cout << "gpu[0..7] : ";
for (int i = 0; i < 8 && i < static_cast(n); ++i) {
std::cout << h_gpu[i] << " ";
}
std::cout << std::endl;
CUDA_CHECK(cudaFree(d_input));
CUDA_CHECK(cudaFree(d_output));
return 0;
}
7. 编译与运行
编译指令十分简洁,以Tesla T4为例:
nvcc -O3 -arch=sm_75 lesson11_prefix_scan.cu -o lesson11_scan
直接运行默认实验:
./lesson11_scan
当然,也支持自定义参数,例如:
./lesson11_scan 16777216 256 10
这三个参数依次对应:元素数量n、block_size以及重复测试次数repeat。
8. 输出
以下是实际运行的部分输出结果:
CUDA Lesson 11: Parallel Prefix Sum / Scan
Elements : 16777216
Data size : 64 MB
Block size: 256
Repeat : 10
[Timing]
CPU exclusive scan time : 23.7002 ms
GPU exclusive scan time : 2.9987 ms
Speedup : 7.9035x
[Check]
Max abs diff: 0.0000
Result : PASS
[Sample]
input[0..7] : 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
gpu[0..7] : 0.0000 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000
可以看到,GPU的执行速度达到了CPU的近8倍,结果完全正确。这充分说明,对于这种看似串行的算法,只要找到合适的并行化结构,GPU依然能展现出巨大的计算威力。
