啥是Parallel Reduction

why so serious

主要的阅读材料有：

这里就说Reduction了，Scan和Reduction的关系可以参考：

0. What & Why & How ?

Reduce：给一组数据，一个满足结合律的二元操作符⊕，那么reduce可以表示为：

和并行啥关系：假设一组数据[1, 2, 3, 4]，⊕就是“ + ”操作，那么顺序执行的reduce就是(((1+2)+3)+4) = 10; 再假设现在有两个线程，如何并行reduce：比如可以两个线程同时分别执行(1+2) 和 (3+4)，再由第一个线程执行之前的结果 (3+7) ，即第一个线程（序号为0）会得到最后reduce的结果。

为什么要并行reduce？这个问题倒可以反过来思考，即为什么不呢？除了一半意义的并行计算的优点，reduce的特性也十分符合并行计算的特点，而Reduce和Scan也确实都是cuda并行计算的核心基础算法，再cuda的官方sample里都有源码实例。

再看一张图片：

直观的是数据多了，20个数据，竖着看，即同时启用10个线程开始计算，最后第一个线程输出结果87。

和之前举的 [1, 2, 3, 4] 的例子有些不同，这个例子每个线程取的两个数据之间的stride是前一组结果的一半。由于结合律我们知道，最终输出的数值是没差的，但是由于gpu硬件和cuda的设计原因，在最终代码实现的时候，这二者的执行效率还是有很大区别。

串行实现就不说了，来看看cuda并行实现，并提出问题，再进行一步步优化。

整体设计：在每个线程块中，使用基于树结构的方法：

并且能够在多个线程块中执行，以满足大数据量的处理，使得每一个线程块对原始数据的部分进行reduce。

那么就有一个问题，由于每个线程块是对部分数据进行reduce，那么其输出结果如何通信，cuda没有global synchronization，咋办？kernel递归，如下图：

那么这个kernel咋写？下面来了：

1. Reduction #1: Interleaved Addressing

首先明确设计思想，如上图，对于给定的输入，存储到shared memory中，设线程id是tid，对于1d的block，由于是两两相加，显然第一次iteration执行 blockDim.x / 2次计算（即有blockDim.x / 2 个激活的线程），第二次执行 blockDim.x / 4 ... 具体实现如下：

可以看到，每个block里reduce都在其shared memory里进行。如上面讨论过的，结果最终在每个block的0号线程里得到，并且存储在其shared memory sdata[0]里，而sdata[0]的值就是这个block输出的reduce的结果。逻辑上可行，但在实现上这里存在一个Warp Divergence问题。

2. Problem #1：Warp divergence

看到for循环中：

GPU是SIMT（Single Instruction Multiple Threads）结构，warp是GPU调度的基本单元，当前一个warp有32个线程，这些线程以不同数据资源执行相同的指令。如果一个warp中的线程做不同的事情时（遇到控制流结构时），例如：

一个warp中的线程分别进入了不同的分支中，这就是warp divergence，那么同一时刻除了正在执行的分之外，其余分支都被阻塞了，也就影响了性能。如果每个warp里的线程只有一个执行路径，那效率就很高咯，但是分支语句有时又不可避免，所以解决的关键在于如何降低warp divergence。（注：不同warp之间不会有warp divergence问题）如何解决：

3. Reduction #2: Interleaved Addressing

小标题没有错哈，仍然是Interleaved Addressing，只不过在 Reduction #1 基础上为了解决 highly divergent branching 问题作了些修改，把原来的for循环改成：

首先去掉了取模运算（取模运算在GPU中开销很大），再有if分支中线程是连续的，缓解了highly divergent branching问题。

但是又有了新的问题：

4. Problem #2：Shared Memory Bank Conflicts

代码中可以看到，reduce是在申请的shared memory中进行的。而这个问题就涉及到cuda中shared memory的分配策略。相比于local memory 和 global memory，shared memory 有着极高的bandwidth和更低的延迟，为此shared memory被划分为等分的内存模块，这些模块称之为bank，bank可以被多个线程同时访问。具体描述可以查看：

简单来说，bank conflict 指的是当同一个warp里的多个线程对同一个bank发出访问请求，只能进行串行访问，极大影响了并行效率。当bank conflict发生时，硬件将该访问请求拆分成多个conflict-free的请求，拆分出来的个数n，就说这个访问请求会引起 n-way bank conflicts。随后cuda引入broadcast机制，使得对同一个bank的同一个地址进行访问时，不会产生bank conflict, 一定程度上缓解了bank conflict，但对同一个bank里的不同地址进行访问时，bank conflict 问题仍然存在。举个例子：

假设有4个banks用于shared memory， bank和所分配的地址关系为：

address:  bank:
       0      0 <-----------------------Thread 0
       1      1
       2      2 
       3      3    
       4      0 <---------Thread 1/ Thread 2
       5      1
       6      2 
       7      3 <-----------------------Thread 3
       8      0 
...

这里：
Thread0 和 Thread3 访问的不同bank，不会发生 bank conflict；

Thread1 和 Thread2 访问的同一个bank的同一个地址，由于broadcas机制的引入，也不会发生bank conflict；

Thread0 和Thread1 访问的同一个bank里的不同地址，发生了bank conflict。

再来看看这张图：

https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#shared-memory-5-x

该图中，左右两侧都没有bank conflict，中间的发生了bank conflict。Reduction#1#2中的取址方式就如同中间这个图，如果能够改进成连续读取如左侧的图，不仅提高效率且避免了bank conflict的发生。

5. Reduction #3: Sequential Addressing

按照之前的讨论，新的方案中，对shared memory中数据的访问应当是：

代码修改也很简单：

上图中，绿色是原来的循环体，改成蓝色的。

还可再改进吗？可以！首先找到问题：

6. Problem #3：Idle Threads

看到for语句：

s初值是 blockDim.x 的一半，当（tid < s）再执行后面指令，那意味着在第一层iteration中就有一半线程是闲置的，这也太浪费了吧！怎么办？看到前面的给shared memory传值的语句：

在这里做做文章：

7. Reduction #4: First Add During Load

改进的思路是，在赋值的时候，就把第一层的reduce给执行了。

前面赋值改为：

这里要注意的是分配的 block 的数量减半，kernel里体现在blockIdx.x的最大值发生了变化，而blockDim是不变的。即每个block中的线程不变，但总的线程数减少了一半。

8. Problem #4：Instruction Overhead

为了继续提升性能，再分析下瓶颈所在。reduction实际上是二元运算，计算强度是比较低的，所以相对的，指令调度的开销可能是瓶颈所在，特别是还有个 for 循环。

那就把 for 给展开吧！

9. Reduction #5: Unroll the Last Warp

之前有提到warp中的32个线程是SIMT策略，意味着，当tid < 32 的时候，并不需要等待其他线程咯，也就是不需要__syncthreads() 指令，而warp的32个线程的又是顺序分配的，步进方式也知道，即右移1位，那么可以把最后的warp展开，而且也不需要每次都 if(tid < s) ：

顺着这个思路，如果所有的迭代次数能够在编译时知道，那就可以把for都给展开咯。

10. Reduction #6: Completely Unrolled

那好就好在，对于现在的GPU，其单个 block 所能容纳的最大线程数是已知的：1024个。随着软硬件更新，这个限制也发生着变化，就比如Mark Harris大神在写《Optimizing Parallel Reduction in CUDA》时，这个数字还是512。不过这不是讨论重点，重点是这个数字已知。

所以将原来的for完全展开：

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(); 
}
if (tid < 32) {
    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]; 
}