(只此一篇便绝b能懂的)cuda reduction操作优化

Reduction操作即为求一个数组的和，按一般的cpu来说我们用树形规约算法就可以达到最优，即对于长度为n的数组，运算次数为n/2+n/4+....+1。那么在GPU中依靠CUDA，最好的算法是什么呢？

一、利用shared内存进行树形规约

template <size_t BLOCK_SIZE>
void __global__ shared_reduction_kernel(
                            const size_t n,
                            const unsigned *src_d,
                            unsigned *tmp_d)
{
    const size_t global_id = threadIdx.x + blockDim.x * blockIdx.x; //线程全局id
    unsigned __shared__ shared[BLOCK_SIZE];
    shared[threadIdx.x] = src_d[global_id];
    for (size_t offset = BLOCK_SIZE >> 1; offset > 0; offset >>= 1)
    {
        __syncthreads();
        if (threadIdx.x < offset)
            shared[threadIdx.x] += shared[threadIdx.x ^ offset]; //按位异或
    }
    if (threadIdx.x == 0)
        tmp_d[global_id / BLOCK_SIZE] = shared[threadIdx.x];
}

   kernel中我们使用模版来表示BLOCK_SIZE而不是blockDim.x，这是为了后面的for循环编译器能够进行展开优化，同时我们也可以更方便地定义shared数组。这里我们来分析这个kernel，首先我们根据线程内置变量threadIdx, blockDim, blockIdx来计算出线程的全局id。再之后我们申请shared内存，并将原数据写入shared内存。这里有些初学cuda的同学可能会有疑问：为什么使用global id与threadIdx匹配，这两个变量的范围明显不一样啊。这里就需要介绍下cuda的内部存储结构，在cuda中每个block分配到到一个SM上，并且拥有一片自己的shared内存，所以对于shared数组是每个block都会申请一个，这样就解决了变量范围不一样的冲突。

接下来的每一次for循环，offset变量为16、8、4、2、1(假定BLOCK大小为32)，通过按位异或操作，我们可以让shared[16]加到shared[0]上，shared[17]加到shared[1]上，这样就完成了规约过程而且相邻的线程(0,1)所访问的16、17也是相邻的，这样也会对程序有好处。这里一些身处第二层但以为我在第一层的同学可能要问了：所谓的shared内存就是可编程的L1cache，既然都已经是cache了，那么访问顺序连续还有什么关系吗？不好意思，我是在第五层。确实如果连续的访问顺序只是为了cache命中率，那么对于shared内存没有任何意义，但这其实是为了避免bank conflict。

众所周知，cuda中线程的执行是以warp为最小单位的，也就是32个线程执行同样的代码，而访问shared内存却是以half warp为基本单位的，也就是同一时间只有16个线程在访问共享内存。因此，共享内存的硬件就被设置成了拥有16个bank，一个线程只能从bank的入口访问shared内存，所以同一时刻最多能有16个线程同时访问shared内存。以特斯拉一块SM上的share内存有16kb为例，每一个bank就有1024字节，可以存储256个int。所以如果shared内存装的是int数据，那么其就可以看作一个256行16列的矩阵，shared[0]与shared[16]都是在第一列也就是bank1，那么如果同一个half warp内的两个线程一个要访问shared[0]一个要访问shared[16]就会产生冲突，这就叫做bank conflict。

其中的__syncthreads()用于同一个block内的线程同步，确保线程们看到的shared内存是一样的。

最后我们将每个block规约出的结果放入tmp_d，就完成了第一级规约(2级后面的时间可以忽略不计)。

二、利用Warp通信进行规约

template<size_t WARP_SIZE>
void __global__ warp_reduction_kernel(
                        const size_t n,
                        const unsigned *src_d,
                        unsigned *tmp_d)
{
    const size_t global_id = threadIdx.x + blockIdx.x * blockDim.x;
    const size_t lane_id = global_id % WARP_SIZE;
    unsigned val = global_id < n ? src_d[global_id] : 0;
    for (size_t offset = WARP_SIZE >> 1; offset > 0; offset >>= 1)
        val += __shfl_xor_sync(0xffffffff, val, offset, WARP_SIZE);
    if (lane_id == 0)
        tmp_d[global_id / WARP_SIZE] = val;
}

下面介绍另一种trick的方法，利用Warp通信的原子操作实现规约，这也是目前性能最好的规约算法。__shfl_xor_sync操作允许我们通过线程的寄存器进行通信，这样做不仅节省了shared内存，而且效率比shared内存更高！

我们看到这里有了一个新变量叫lane_id，这个变量的取值是0～31，表示线程在一个warp里的编号。for循环内的操作将一个warp内的所有线程的val变量累加起来，最后lane_id=0的线程，也就是warp内的第一个线程的val变量就是整个warp的和。

下面解释一下__shfl_xor_sync的作用，它将第一个参数是一个标志warp内有效线程的mask，由于一个warp有32个线程，如果对全部warp内的线程生效就设置成8个f，这样就是32个1了。这个函数的作用就是从当前线程的lane_id值与offset做xor，然后把结果处的val变量值取来。