cuda reduction&reduce

概念

reduction

是一种并行计算中的操作概念或技术，指的是将一组数据通过某种特定的操作（如加法、乘法、求最大值、求最小值等）进行聚合，最终得到一个或几个汇总结果的过程。它强调的是这种数据处理的模式或任务类型，是一个较为宽泛的概念，在各种并行计算场景和框架中都可能存在。

reduce

通常是指 CUDA 中实现 reduction 操作的具体函数或方法。在 CUDA 的编程接口中，“reduce” 是用于执行归约操作的具体函数名或函数模板，是实现 “reduction” 概念的具体工具，是一个具体的编程元素。

代码1-初始

假设处理N个元素，划分M个block，然后每个线程处理一个元素，每个block一共M个thread

__global__ void reduce0(float *d_in,float *d_out){__shared__ float sdata[THREAD_PER_BLOCK];//申请共享内存，大小存放的是一个block里的thread数量//each thread loads one element from global memory to shared memunsigned int i=blockIdx.x*blockDim.x+threadIdx.x;//当前线程要索引的数据对应全局数组位置unsigned int tid=threadIdx.x;//当前线程在block里面的位置sdata[tid]=d_in[i];//这个过程是global到share__syncthreads();//同步，因为有依赖关系// do reduction in shared mem ，这里就是简单的二分合并，share对于一个block的所有线程可见，所以可以这样在share上合并去做。for(unsigned int s=1; s<blockDim.x; s*=2){if(tid%(2*s) == 0){sdata[tid]+=sdata[tid+s];}__syncthreads();//计算有依赖所以同步}// write result for this block to global memif(tid==0)d_out[blockIdx.x]=sdata[tid];//最后合并的值放在数组的最开始
}

假设32M数据，每个thread处理一个数据，一个block256个thread，
1M=1024K，所以一共32*1024/256=128个block

代码2-减少warp内分支

对于cuda来说，假设代码存在分支AB，因为是simt，所以会先执行A，再执行B，而不是并行执行AB，所以对于代码1的优化先是减少分支

__global__ void reduce1(float *d_in,float *d_out){__shared__ float sdata[THREAD_PER_BLOCK];//each thread loads one element from global memory to shared memunsigned int i=blockIdx.x*blockDim.x+threadIdx.x;unsigned int tid=threadIdx.x;sdata[tid]=d_in[i];__syncthreads();// do reduction in shared memfor(unsigned int s=1; s<blockDim.x; s*=2){int index = 2*s*tid;if(index < blockDim.x){sdata[index]+=sdata[index+s];}__syncthreads();}// write result for this block to global memif(tid==0)d_out[blockIdx.x]=sdata[tid];
}

这里的思路是256个数据有八个warp，为了不让warp出现较少分支，比如第一个循环，只让warp0~3的线程去做累加计算，4·7的warp执行另外一个分支，这样虽然空转的线程数没有减少，但是warp内没有分支，性能提升。

代码3-减少bank冲突

__global__ void reduce2(float *d_in,float *d_out){__shared__ float sdata[THREAD_PER_BLOCK];//each thread loads one element from global memory to shared memunsigned int i=blockIdx.x*blockDim.x+threadIdx.x;unsigned int tid=threadIdx.x;sdata[tid]=d_in[i];__syncthreads();// do reduction in shared memfor(unsigned int s=blockDim.x/2; s>0; s>>=1){if(tid < s){sdata[tid]+=sdata[tid+s];}__syncthreads();}// write result for this block to global memif(tid==0)d_out[blockIdx.x]=sdata[tid];
}

需要以warp为单位去避免bank

当朴素的将数据搬运到share的时候，一个数据占一个bank，那么第一轮循环中，code2的thread0占用为了索取01数据占据01bank，thread16为了索引数据3233也占据了bank01，这两个thread此时就发生了两路bank冲突。
如果不改变数据存放的方式，就需要改变索引。
对于朴素存放，数据在share中如下
0 ~ 31
32 ~ 63
64 ~ 95
96 ~127
128~159
160~191
192~223
224~255

第一轮循环：
对于修改代码，第一轮循环thread索引0和128，而这两个元素都在bank0上，对于底层代码来说，这部分计算设计两个load，第一个load，warp0的每个线程都去load0~31，这里不会冲突，并且底层会优化为128bit的大load。第二个load同理，只不过换了不同行。
第一个循环将128-255的值加到了0-128位置上
第二轮循环：
thead0 ： 0+64：同理，也是不同行
第三轮循环：
thread 0：0+32：同理
第四轮循环
thread 0：0+16，但此时因为tid<16才计算，所以也没有冲突
第五轮循环
0+8
六
0+4
七
0+2
八
0+1
得到结果存在0

标准写法

#include <iostream>
#include <cuda_runtime.h>// CUDA 内核函数：规约加法
__global__ void reductionAdd(float *input, float *output, int n) {__shared__ float partialSum[256];// 线程索引int tid = threadIdx.x;int idx = blockIdx.x * blockDim.x + threadIdx.x;// 加载数据到共享内存partialSum[tid] = (idx < n) ? input[idx] : 0.0f;__syncthreads();// 规约操作for (int s = blockDim.x / 2; s > 0; s >>= 1) {if (tid < s) {partialSum[tid] += partialSum[tid + s];}__syncthreads();}// 将每个块的部分和写入全局内存if (tid == 0) {output[blockIdx.x] = partialSum[0];}
}// 主机端函数：计算最终结果
float finalSum(float *output, int numBlocks) {float sum = 0.0f;for (int i = 0; i < numBlocks; i++) {sum += output[i];}return sum;
}int main() {const int n = 1024;float *h_input = new float[n];float *d_input, *d_output;// 初始化输入数据for (int i = 0; i < n; i++) {h_input[i] = static_cast<float>(i + 1);}// 分配设备内存cudaMalloc((void**)&d_input, n * sizeof(float));cudaMalloc((void**)&d_output, (n / 256) * sizeof(float));// 将数据从主机复制到设备cudaMemcpy(d_input, h_input, n * sizeof(float), cudaMemcpyHostToDevice);// 定义线程块和网格的维度dim3 blockSize(256);dim3 gridSize((n + blockSize.x - 1) / blockSize.x);// 调用 CUDA 内核函数reductionAdd<<<gridSize, blockSize>>>(d_input, d_output, n);// 同步设备cudaDeviceSynchronize();// 将部分和从设备复制到主机float *h_output = new float[gridSize.x];cudaMemcpy(h_output, d_output, gridSize.x * sizeof(float), cudaMemcpyDeviceToHost);// 计算最终结果float result = finalSum(h_output, gridSize.x);// 输出结果std::cout << "The sum of the array is: " << result << std::endl;// 释放内存delete[] h_input;delete[] h_output;cudaFree(d_input);cudaFree(d_output);return 0;
}