CUDA自学笔记01—Reduction规约求和

备注

这个笔记是在自学过程的一些记录，中间包含了很多个人的理解以及问AI的回答，因为找不到很系统的教程，官方文档有时候不会讲的特别细，比方说某个指标的具体计算过程网上也没找到什么资料，所以在这里记录下，后续自己有了新的理解或者发现写的不对的地方都会在这不断更新，如有不对的地方，还请指正。

Reduction规约求和

Reduction就是对一组数据经过一些操作变换成一个结果，这里以求和为例进行记录。

每个数加上它右边间隔s的数，开始s=1，每迭代一轮s*2，迭代log₂ⁿ轮，n是数据量，这里假设n是2的幂次。迭代完后第一个数就是求和结果。

测试过程数据是4000000个int的求和，然后默认blocksize=512，显卡是RTX 4070 Super。

Reduce朴素版本

kernel代码如下：

__global__voidreduce_naive(int*g_idata,int*g_odata){inttid=threadIdx.x;inti=blockIdx.x*blockDim.x+threadIdx.x;for(ints=1;s<blockDim.x;s*=2){if(tid%(2*s)==0){g_idata[i]+=g_idata[i+s];}__syncthreads();}if(tid==0)g_odata[blockIdx.x]=g_idata[blockIdx.x*blockDim.x];}

这段代码相当于在每个block的512个thread里面做规约，g_odata里面是每个block里512个数据的和，最后把g_odata求和就得到最终结果。

Reduce1——使用共享内存

我们都知道，共享内存能够帮助我们减少内存搬运耗时，因为它是片上内存，我们可以先把数据搬运到共享内存上再去计算。

kernel代码如下：

__global__voidreduce1(int*g_idata,int*g_odata){extern__shared__intsdata[];inttid=threadIdx.x;inti=blockIdx.x*blockDim.x+threadIdx.x;sdata[tid]=g_idata[i];__syncthreads();for(ints=1;s<blockDim.x;s*=2){if(tid%(2*s)==0){sdata[tid]+=sdata[tid+s];}__syncthreads();}if(tid==0)g_odata[blockIdx.x]=sdata[0];}

这里我用Nsight Compute分析对比上面两个kernel时发现耗时没有明显降低，然后上网搜索了一下说是Naive的全局显存数据请求到数据后会暂存在L2 cache中复用，所以相对于共享内存版本耗时差异不大。从下表中可以看到Naive版本的L2 cache throughput很高。（PS：这里其实没有完全理解，在这里Mark一下。我当前的猜测是因为大部分操作数据在内存上都是相邻的，一次事务可以取32个4字节的int数，所以内存事务复用率较高，所以和共享内存版本差异不大）

Reduce2——去除分支发散和除法操作

Reduce1有两个问题，一个是分支发散比较严重，另一个是’%'操作比较耗时，需要优化一下。

__global__voidreduce2(int*g_idata,int*g_odata){extern__shared__intsdata[];inttid=threadIdx.x;inti=blockIdx.x*blockDim.x+threadIdx.x;sdata[tid]=g_idata[i];__syncthreads();for(ints=1;s<blockDim.x;s*=2){intidx=2*s*tid;if(idx<blockDim.x){sdata[idx]+=sdata[idx+s];}__syncthreads();}if(tid==0)g_odata[blockIdx.x]=sdata[0];}

从下表可以看出我们实现了规约操作，且减少了分支发散和除法操作，但是这里带来了一个新的问题就是引入了Bank conflict。例如第一轮迭代tid=0的线程会读取idx=0的数据，tid=16的线程会读取idx=32的数据，而这两个地址的数据在同一Bank0（idx%32=0）里面，这就相当于不同线程在读取同一Bank的不同地址的数据，这就是Bank conflict，导致并行变成串行，增加耗时。不过相比于前面的kernel，还是能看到耗时明显降低，说明之前的kernel分支发散和除法操作比bank conflict更耗时。

bank conflict计算

先按一个512个int、blocksize=512的demo计算bank conflicts

指标定义参考Nvidia Nsight Compute官方文档

对于512个线程，每个线程处理一个int数据，可以分为512/32=16个warp

Share Store：
1、对于sdata[tid] = g_idata[i];这是在给共享内存写入，所以对应Share Store，总共16个warp，对应16个Instructions和Requests；这里由于都是连续写入，所以不存在bank conflicts，wavefronts=16

2、对于sdata[idx] += sdata[idx + s]，第一次迭代也就是s=1的时候是只在idx——0~255的索引下写入，所以只有8个warp，总共9次迭代依次是8、4、2、1、1、1、1、1、1个warp在做共享内存写入，所以对应8+4+2+1+1+1+1+1+1=20个Instructions和Requests；

这里我举了s=1和s=16时的索引，可以看到s=1时，idx和idx+s都是间隔2的索引，就是一个2-way的bank conflict，一个warp就需要2次transactions，依次类推s=2时4-way的bank conflict需要4次transactions，而到了s=16时，注意tid=16时，idx已经到了512，索引已经超出了，超出的这部分是无效的，虽然这个warp是32-way的，但是我只需要发送16次transactions就满足了，后续迭代以此类推。最终得到1次load或1次store的wavefronts=16+16+16+16+16+8+4+2+1=95

综合下来Share Store Instructions和Requests = 16 + 20 = 36，wavefronts=95+16=111

Share Load：
1、对于sdata[idx] += sdata[idx + s]，分析思路和Share Store一致，不过有个区别是这里每个线程读取是有两个地址idx和idx+s，所以Share Load Instructions和Requests = 20 + 20 = 40；wavefronts参考前面的Store分析wavefronts=95*2=190（乘2是因为做了两次读取）

2、if (tid == 0) g_odata[blockIdx.x] = sdata[0];这行语句也有一个读取共享内存操作，不过只做一次； sdata[0]相当于做了一次内存读取，也就相当于一次wavefront
综合下来Share Load Instructions和Requests = 40 + 1 = 41，wavefronts=190+1=191

而bank conflicts就相当于额外的wavefronts，结果也就等于总的wavefronts减去请求的requests
Share Load Bank Conflicts = 191 - 41 = 150
Share Store Bank Conflicts = 111 - 36 = 75
注意到这里Load的冲突数刚好是Store的两倍，这是因为发生冲突的语句sdata[idx] += sdata[idx + s]刚好是两次读取一次写入。

我后面又尝试了不同尺寸的输入，发现Share Load wavefronts不能通过单个block的wavefronts * 总的blocks计算（而且每次跑可能会不一样，虽然相比于总数不会差太多），其他几个都可以。问了下AI说是Wavefronts 是流水线层面的“物理动作”。在 GPU 的硬件流水线中，当某几个 Warp 的执行由于某些外部原因（例如抢占、上下文切换、或规约中的 __syncthreads() 同步导致暂停再恢复时），LSU 阶段可能为了确保数据的危害（Hazard）检测通过，在极个别周期发生了物理重复请求（Retry）。这种极低概率的硬件自发抖动，会被物理计数器记录下来。

Reduce3——减少Bank Conflicts

我们可以将求和的间隔由小到大改成由大到小，如下图所示，这样就可以避免绝大多数的Bank conflict

可以看到Bank Conflits明显减少

Reduce4——减少空闲线程

前面的Reduce方法有比较多的线程冗余，比如说我们每个block开了512个线程，但是后面256个线程只做了一次数据搬运就没有其他操作了，第一轮迭代只有前面256个线程在做加法运算运行，后面256个线程就都空置了，后续每次迭代空置的线程也就越来越多。
为了把线程利用起来，可以在把数据往共享内存搬运的时候就做一次加法，比方说我一个block还是512个线程，但是我用这512个线程处理1024个数据，在把数据从全局内存搬运到共享内存的时候，我读取第0个数据和第512个数据，把他们相加然后写入共享内存，然后一直到第511个数据和第1023个数据，把他们相加然后写入共享内存，这样就相当于每个线程都做了一次加法操作，相比上面的Reduce把空闲线程也做了利用。

__global__voidreduce4(int*g_idata,int*g_odata,intsize){extern__shared__intsdata[];inttid=threadIdx.x;inti=blockIdx.x*(blockDim.x*2)+threadIdx.x;if(i<size)sdata[tid]=g_idata[i]+g_idata[i+blockDim.x];elsesdata[tid]=0;__syncthreads();for(ints=blockDim.x/2;s>0;s>>=1){if(tid<s){sdata[tid]+=sdata[tid+s];}__syncthreads();}if(tid==0)g_odata[blockIdx.x]=sdata[0];}

这里需要注意，由于我们每个block是用512个线程处理1024个数据，所以总的block数相比于之前的Reduce要减半。

可以看到kernel耗时显著降低

Reduce5——warp展开

我们知道一个warp里面有32个线程，且是同步执行的，所以当间隔s从大间隔降到32的时候，我们可以不用 __syncthreads()命令，直接把内部规约展开，注意kenel中for循环中止条件从s>0变成了s>32

__device__voidwrapReduce5(volatileint*sdata,inttid){sdata[tid]+=sdata[tid+32];sdata[tid]+=sdata[tid+16];sdata[tid]+=sdata[tid+8];sdata[tid]+=sdata[tid+4];sdata[tid]+=sdata[tid+2];sdata[tid]+=sdata[tid+1];}__global__voidreduce5(int*g_idata,int*g_odata,intsize){extern__shared__intsdata[];inttid=threadIdx.x;inti=blockIdx.x*(blockDim.x*2)+threadIdx.x;if(i<size)sdata[tid]=g_idata[i]+g_idata[i+blockDim.x];elsesdata[tid]=0;__syncthreads();for(ints=blockDim.x/2;s>32;s>>=1){if(tid<s){sdata[tid]+=sdata[tid+s];}__syncthreads();}if(tid<32)wrapReduce5(sdata,tid);if(tid==0)g_odata[blockIdx.x]=sdata[0];}

可以看到kernel耗时进一步降低

Reduce6——for循环进一步展开

可以进一步展开kernel中的for循环

template<unsignedintblockSize>__device__voidwrapReduce(volatileint*sdata,inttid){if(blockSize>=64)sdata[tid]+=sdata[tid+32];if(blockSize>=32)sdata[tid]+=sdata[tid+16];if(blockSize>=16)sdata[tid]+=sdata[tid+8];if(blockSize>=8)sdata[tid]+=sdata[tid+4];if(blockSize>=4)sdata[tid]+=sdata[tid+2];if(blockSize>=2)sdata[tid]+=sdata[tid+1];}template<unsignedintblockSize>__global__voidreduce6(int*g_idata,int*g_odata,intsize){extern__shared__intsdata[];inttid=threadIdx.x;inti=blockIdx.x*(blockDim.x*2)+threadIdx.x;if(i<size)sdata[tid]=g_idata[i]+g_idata[i+blockDim.x];elsesdata[tid]=0;__syncthreads();// 展开循环if(blockSize>=1024){if(tid<512){sdata[tid]+=sdata[tid+512];}__syncthreads();}if(blockSize>=512){if(tid<256){sdata[tid]+=sdata[tid+256];}__syncthreads();}if(blockSize>=256){if(tid<128){sdata[tid]+=sdata[tid+128];}__syncthreads();}if(blockSize>=128){if(tid<64){sdata[tid]+=sdata[tid+64];}__syncthreads();}// Warp 级优化if(tid<32){// 强制转换为 volatile 指针传递，确保读写可见性wrapReduce<blockSize>((volatileint*)sdata,tid);}if(tid==0)g_odata[blockIdx.x]=sdata[0];}

不过耗时并没有太大变化，但是Mark Harris当时测试时有一定提升，猜测是这些年硬件有一定的迭代优化

Reduce7——网格步幅循环加载

Mark Harris文档里最后还提了一个网格步幅循环加载（Grid-Stride Loop）的规约，原理是考虑数据量很大，然后以有限的block有限的线程去做运算，就相当于每个线程处理ceil(size/gridsize)个数据

template<unsignedintblockSize>__global__voidreduce7(int*g_idata,int*g_odata,intsize){extern__shared__intsdata[];inttid=threadIdx.x;inti=blockIdx.x*(blockSize*2)+threadIdx.x;intgridSize=blockSize*2*gridDim.x;sdata[tid]=0;while(i<size){sdata[tid]+=g_idata[i];if(i+blockSize<size){// 防止越界sdata[tid]+=g_idata[i+blockSize];}i+=gridSize;}__syncthreads();// 展开循环if(blockSize>=1024){if(tid<512){sdata[tid]+=sdata[tid+512];}__syncthreads();}if(blockSize>=512){if(tid<256){sdata[tid]+=sdata[tid+256];}__syncthreads();}if(blockSize>=256){if(tid<128){sdata[tid]+=sdata[tid+128];}__syncthreads();}if(blockSize>=128){if(tid<64){sdata[tid]+=sdata[tid+64];}__syncthreads();}// Warp 级优化if(tid<32){// 强制转换为 volatile 指针传递，确保读写可见性wrapReduce<blockSize>((volatileint*)sdata,tid);}if(tid==0)g_odata[blockIdx.x]=sdata[0];}

这里我就没去调节block数了，直接和上面的kernel用的一样的block数

Reduce8——warp shuffle

现在的CUDA有个shuffle的功能支持warp间不同线程做一些数据处理，等价于把之前的取数相加赋值改成__shfl_down_sync函数调用

template<unsignedintblockSize>__global__voidreduce8(constint*__restrict__ g_idata,int*__restrict__ g_odata,intsize){extern__shared__intsdata[];unsignedinttid=threadIdx.x;unsignedintidx=blockIdx.x*(blockSize*2)+tid;intsum=0;// 处理不能整除情况if(idx<size)sum+=g_idata[idx];if(idx+blockSize<size)sum+=g_idata[idx+blockSize];sdata[tid]=sum;__syncthreads();//--------------------------------------------------// Block Reduction//--------------------------------------------------if(blockSize>=1024){if(tid<512)sdata[tid]+=sdata[tid+512];__syncthreads();}if(blockSize>=512){if(tid<256)sdata[tid]+=sdata[tid+256];__syncthreads();}if(blockSize>=256){if(tid<128)sdata[tid]+=sdata[tid+128];__syncthreads();}if(blockSize>=128){if(tid<64)sdata[tid]+=sdata[tid+64];__syncthreads();}//--------------------------------------------------// Warp Reduction using shuffle//--------------------------------------------------if(tid<32){intval=sdata[tid];if(blockSize>=64)val+=sdata[tid+32];for(intoffset=16;offset>0;offset>>=1){val+=__shfl_down_sync(0xffffffff,val,offset);}if(tid==0)g_odata[blockIdx.x]=val;}}

不过后面几种Reduce方法耗时都比较接近了

参考文献

Optimizing Parallel Reduction in CUDA