游乐游手机版
首页/AI教程/文章详情

CUDA编程并行前缀和扫描详解

时间:2026-07-07 15:39
前缀和(PrefixSum)有包含自身和排除自身两种形式,其中排他扫描常用于地址计算。由于存在数据依赖,简单并行失效。BlellochScan通过向上规约构建部分和树、向下传播获得结果,并行复杂度为O(logn)。实验实现了GPU多轮扫描并与CPU进行性能对比。

在前一课(第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正是破解该难题的利器。其核心思想分为两个阶段:

  1. Up-sweep(向上规约):构建一棵部分和树。与我们之前学习的Reduction类似,两两相加,再四四相加,直至得到整个数组的总和。
  2. 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. 实验内容

光说不练假把式。本次实验需要动手实现以下五个部分:

  1. 一个CPU版本的exclusive scan,作为性能对比的基准。
  2. 一个基于GPU block-level的Blelloch scan。
  3. 一个能够处理大数组的多轮GPU scan方案。
  4. 对GPU结果进行校验,确保与CPU结果一致。
  5. 最后,进行令人期待的性能对比分析。

在这里插入图片描述

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依然能展现出巨大的计算威力。

来源:https://cloud.tencent.com.cn/developer/article/2703887
上一篇RAG与微调区别别再背概念面试解析 下一篇阿里云百炼Token Plan支持OpenAI和Anthropic及主流AI编程助手
本站内容用于信息整理与展示,如有侵权或内容问题请及时联系处理。

相关推荐

补充同频道和同主题内容,方便继续浏览更多相关内容。

同类最新

继续查看同栏目最近更新的文章。

更多
科研人员必读:多肽、蛋白质、重组蛋白区别及定制指南
AI教程 · 2026-07-07

科研人员必读:多肽、蛋白质、重组蛋白区别及定制指南

Section 01 多肽 VS 蛋白质 VS 重组蛋白 多肽、蛋白质和重组蛋白,本质上是同宗同源的东西——都是氨基酸串起来的生物大分子。三者的核心区别,说到底无非是三个维度:分子大小、折叠形态,以及生产方式。 接下来是一张清晰的对比图,帮你快速建立直觉: ![对比图1](https:

知识图谱与本体语义建模的核心区别解析
AI教程 · 2026-07-07

知识图谱与本体语义建模的核心区别解析

谈到人工智能如何“理解”知识,有两个概念常被放在一起讨论:知识图谱与本体语义建模。不少人以为它们是同一事物,或者认为后者是前者的进化版。实际上,两者的分工完全不同——打个比方,一个是“记事的本子”,另一个是“写本子之前先定好的规矩”。 1 本体语义建模:先绘制一张“通用分类蓝图” 设想一下,你要整

强烈推荐工作搭子WorkBuddy
AI教程 · 2026-07-07

强烈推荐工作搭子WorkBuddy

一次偶然的机会,从朋友那里了解到WorkBuddy这个工具。说实话,在AI产品扎堆的今天,能遇到一个下载即用的助手,确实值得推荐给每一个被日常琐事缠身的人。 安装过程没什么难度,双击安装包默认安装即可。需要留意的是,如果在Windows7上折腾了半天没反应,别慌——这工具在高版本Windows下运行

跨境电商系统自动化测试与CI/CD流水线构建指南
AI教程 · 2026-07-07

跨境电商系统自动化测试与CI/CD流水线构建指南

技术方向:自动化测试与DevOps实践关键词:日本代购、一站式日淘、雅虎代拍系统、煤炉自动代拍 一、测试分层策略详解 不少人刚开始就想直接搞E2E测试,觉得跑通完整流程才够“真实”。然而,测试金字塔这么多年仍不过时,原因很简单——不同层级的测试各有分工,缺少任何一层都会不稳。来看看这张金字塔图: ┌

中小企业AI营销矩阵工具推荐:赛诺贝斯智域蒲公英
AI教程 · 2026-07-07

中小企业AI营销矩阵工具推荐:赛诺贝斯智域蒲公英

天天刷着别人的爆款内容,自己却“有心无力”——这才是2026年绝大多数中小企业运营社交媒体的真实写照。说白了,社交媒体如今早已不是“要不要做”的选择题,而是“怎么做才能真正见效”的生存考验。现实情况是,团队人力就那么几个,预算也紧巴巴,却要同时运营抖音、小红书、知乎、头条、百家号等多个阵地……文案、