一、前言
使用 CUDA 并行计算可以加速数组所有元素之和，由原来的逐一累加，变为多区间并行累加，并将多区间的和再累加，在数组较大时，有明显的加速效果。
该小结还使用了共享内存，提高访问效率。
二、CPU 计算
这里直接放源码：
int _sum_cpu(int *ptr, int count)
{
int sum = 0;
for (int i = 0; i < count; i++)
{
sum += ptr[i];
}
return sum;
}这是一个非常简单的 for 循环数组求和，CPU 需要依次读取计算，属于单线程依赖计算。如果改用多线程，就需要分区间进行，防止数据竞争和重复计算。
三、GPU 并行计算
对于较大的数组，会将每个所有元素均分给每个线程。例如，数组总数为 10000000，grid 大小为 65536，也就是线程数，那么每个线程负责 ceil(10000000 / 65536) = 153 个元素的计算。这样做的目的是实现动态数组大小的适配，不依赖线程数。具体做法为：
每个线程从自己的起始位置开始，以整个网格的线程总数为步长（即65536），跳跃地读取数据并进行累加，核心目的是将超出线程数量的原始数据，压缩到与线程数量相等的局部和数据（10000000 -> 65536），每个线程处理元素数（for循环次数）：ceil(10000000 / 65536) = 153。
下面是 CUDA 核函数的源码，有很多注释辅助理解：
__global__ void _sum_gpu(int *input, int count, int *output) {
__shared__ int sum_per_block[BLOCK_SIZE];
int temp = 0;
// 每个线程从自己的起始位置开始，以整个网格的线程总数为步长（即65536），跳跃地读取数据并进行累加
// 核心目的是将超出线程数量的原始数据，压缩到与线程数量相等的局部和数据（10000000 -> 65536）
// 每个线程处理元素数（for循环次数）：ceil(10000000 / 65536) = 153
for (int idx = threadIdx.x + blockDim.x * blockIdx.x;
idx < count;
idx += gridDim.x * blockDim.x)
{
temp += input[idx];
}
// 将各线程的局部和存入共享内存
// 这256个block是独立在SM上调度的，每个block中的256个线程共享同一个共享内存空间，但不同block的共享内存是完全隔离的
sum_per_block[threadIdx.x] = temp;
// 确保当前block中所有线程完成写入后再继续
__syncthreads();
//**********shared memory summation stage***********
// 每次迭代将数据规模减半（类似折半求和）
for (int length = BLOCK_SIZE / 2; length >= 1; length /= 2)
{
int double_kill = -1;
// 仅低索引线程参与计算（每次减半）
if (threadIdx.x < length)
{
// 例：1 2 3 4，这里将1和3相加，2和4相加
double_kill = sum_per_block[threadIdx.x] + sum_per_block[threadIdx.x + length];
}
// 确保所有线程完成对共享内存的读取（避免部分线程写入前被读取）
__syncthreads();
if (threadIdx.x < length)
{
// 结果回写到共享内存的前半段（覆盖旧值）
sum_per_block[threadIdx.x] = double_kill;
}
// 确保写入操作完成后再进入下一轮归约（防止脏数据）
__syncthreads();
} // the per-block partial sum is sum_per_block[0]
// 跳过超出数组范围的线程块（如网格过大时），避免无效线程块的冗余原子操作
if (blockDim.x * blockIdx.x < count)
{
// the final reduction performed by atomicAdd()
// 每个block的 0 号线程将块部分和累加到全局输出
// atomicAdd确保跨线程块的原子操作
if (threadIdx.x == 0)
atomicAdd(output, sum_per_block[0]);
}
}需要注意的是，代码中多次调用了线程同步函数：
__syncthreads();
每一步的作用不一样，读者可以结合注释理解。
四、源文件代码
sum.cu
#include<stdio.h>
#include<stdint.h>
#include"error.cuh"
#include<time.h>
#include<stdlib.h> //for srand()/rand()
#include<sys/time.h> //for gettimeofday()/struct timeval
#define N 10000000
#define BLOCK_SIZE 256
//#define BLOCKS ((N + BLOCK_SIZE - 1) / BLOCK_SIZE)
#define BLOCKS 256
// CUDA 的统一内存（Unified Memory）是通过 cudaMallocManaged 函数和 __managed__ 关键字实现的主机与设备的透明化内存访问
// 其核心原理是将物理存储位置抽象为统一的虚拟地址空间，当 CPU 或 GPU 访问数据时，系统自动完成数据迁移（按需分页迁移），开发者无需手动调用 cudaMemcpy
// 统一内存消除了通过 cudaMemcpy*() 例程进行显式数据移动的需要，而不会因将所有数据放入零拷贝内存而导致性能损失
__managed__ int source[N]; //input data
__managed__ int final_result[1] = {0}; //scalar output
void _init(int *ptr, int count)
{
uint32_t seed = (uint32_t)time(NULL);
srand(seed); //reseeding the random generator
//filling the buffer with random data
for (int i = 0; i < count; i++) ptr[i] = rand();
}
double get_time()
{
struct timeval tv;
gettimeofday(&tv, NULL);
return ((double)tv.tv_usec * 0.000001 + tv.tv_sec);
}
__global__ void _sum_gpu(int *input, int count, int *output) {
__shared__ int sum_per_block[BLOCK_SIZE];
int temp = 0;
// 每个线程从自己的起始位置开始，以整个网格的线程总数为步长（即65536），跳跃地读取数据并进行累加
// 核心目的是将超出线程数量的原始数据，压缩到与线程数量相等的局部和数据（10000000 -> 65536）
// 每个线程处理元素数（for循环次数）：ceil(10000000 / 65536) = 153
for (int idx = threadIdx.x + blockDim.x * blockIdx.x;
idx < count;
idx += gridDim.x * blockDim.x)
{
temp += input[idx];
}
// 将各线程的局部和存入共享内存
// 这256个block是独立在SM上调度的，每个block中的256个线程共享同一个共享内存空间，但不同block的共享内存是完全隔离的
sum_per_block[threadIdx.x] = temp;
// 确保当前block中所有线程完成写入后再继续
__syncthreads();
//**********shared memory summation stage***********
// 每次迭代将数据规模减半（类似折半求和）
for (int length = BLOCK_SIZE / 2; length >= 1; length /= 2)
{
int double_kill = -1;
// 仅低索引线程参与计算（每次减半）
if (threadIdx.x < length)
{
// 例：1 2 3 4，这里将1和3相加，2和4相加
double_kill = sum_per_block[threadIdx.x] + sum_per_block[threadIdx.x + length];
}
// 确保所有线程完成对共享内存的读取（避免部分线程写入前被读取）
__syncthreads();
if (threadIdx.x < length)
{
// 结果回写到共享内存的前半段（覆盖旧值）
sum_per_block[threadIdx.x] = double_kill;
}
// 确保写入操作完成后再进入下一轮归约（防止脏数据）
__syncthreads();
} // the per-block partial sum is sum_per_block[0]
// 跳过超出数组范围的线程块（如网格过大时），避免无效线程块的冗余原子操作
if (blockDim.x * blockIdx.x < count)
{
// the final reduction performed by atomicAdd()
// 每个block的 0 号线程将块部分和累加到全局输出
// atomicAdd确保跨线程块的原子操作
if (threadIdx.x == 0)
atomicAdd(output, sum_per_block[0]);
}
}
int _sum_cpu(int *ptr, int count)
{
int sum = 0;
for (int i = 0; i < count; i++)
{
sum += ptr[i];
}
return sum;
}
int main() {
fprintf(stderr, "filling the buffer with %d random elements...\n", N);
_init(source, N);
cudaDeviceSynchronize();
fprintf(stderr, "Running on GPU...\n");
double t0 = get_time();
_sum_gpu<<<BLOCKS, BLOCK_SIZE>>>(source, N, final_result);
CHECK(cudaGetLastError()); // checking for launch failures
CHECK(cudaDeviceSynchronize()); // checking for run-time failurs
double t1 = get_time();
int A = final_result[0];
fprintf(stderr, "GPU sum: %u\n", A);
fprintf(stderr, "Running on CPU...\n");
double t2 = get_time();
int B = _sum_cpu(source, N);
double t3 = get_time();
fprintf(stderr, "CPU sum: %u\n", B);
//******The last judgement**********
if (A == B)
{
fprintf(stderr, "Test Passed!\n");
}
else
{
fprintf(stderr, "Test failed!\n");
exit(-1);
}
//****and some timing details*******
fprintf(stderr, "GPU time %.3f ms\n", (t1 - t0) * 1000.0);
fprintf(stderr, "CPU time %.3f ms\n", (t3 - t2) * 1000.0);
return 0;
}运行结果：
