本周主要内容如下：
- 如何分析GPU算法的速度和效率
- 3个新的基本算法：归约、扫描和直方图(Reduce、Scan、Histogram)
一、评估标准
首先介绍用于评估GPU计算的两个标准：
- step ：完成某特定计算所需时间--挖洞操作（Operation Hole Digging）
- work：工作总量
如下图示，第一种情况只有一个工人挖洞，他需要8小时才能完成，所以工作总量（Work）是8小时。第二种情况是有4个工人，它们2个小时就能完成挖洞任务，此时工作总量是8小时。第三种情况同理不加赘述。
另一种直观理解如下图示。该计算任务步骤复杂度是3，工作复杂度是7。而接下来的课程的目的则是学会如何优化GPU算法。
二、3个新的基本算法
2.1 Reduce
2.1.1 Reduce运算基本介绍
下图展示的是reduce运算。
reduce运算由两部分组成：输入值和 归约运算符（Reduction Operation）。
其中归约运算符需要满足下面两个条件：
- Binary 二元运算符：即对两个输入对象进行操作得到一个结果
- Associative 结合性：直白解释就是运算符需要与顺序无关，例如加号就满足这一性质，而减号不满足。
2.1.2 串行方式实现归约和并行方式实现归约
- 串行方式实现归约
reduce运算和map运算有点类似，但不同的地方在于map是可以对多个元素做单独操作的，而reduce每次计算都需要依赖上一次运算的结果。示意图如下：
- 并行方式实现归约
并行方式与串行方式的不同点在于计算的方式不同，如果以穿行方式对a,b,c,d进行加法归约运算，那么计算方式则为\(((a+b)+c)+d\)。而如果以并行方式的话，则为\((a+b)+(c+d)\)。
两种方式的复杂度如下图示：
上面的例子貌似并不能很明显的看出串行和并行的差异，下面我们通过大致计算来更直观的看出不同：
首先假设使用的并行的方式计算，当输入对象有两个，那么step=1。当有4个输入对象，那么就有step=2。以此类推，当有N个输入对象，那么step=logN。而串行的step=N-1。因此当N一定大的时候，并行方式的优势就很明显了。
下面通过具体示例来看看使用 Global Memory和 Shared Memory进行并行归约运算的区别。
如下图示，一共有1024个线程块，其中每个线程块有1024个线程需要进行归约运算。运算分成两个步骤，首先块内运算得到1个值，之后再将1024个值再做归约运算。
下面给出了使用全局变量的代码。实现思路是将线程分成两份，然后等距相加。
具体来说就是首先将后512个线程值添加对应添加到前512个线程中去，之后再折半添加，即将前面512个线程拆成前256和后256个线程，以此类推，直到最后计算得到一个值。
下面给出了使用 Shared Memory的代码，它的思路与上面类似，但是因为使用共享内存，速度得到了提升
理论上来说使用Global Memory需要进行2n次读操作，n次写操作，而Shared Memory需要进行n次读操作，1次写操作，所以后者应该比前者快3倍左右。
但是因为该示例并没有使得计算饱和，所以最终测试结果并没有3倍，具体时间效果如下图。
2.2 Scan
2.2.1 Scan基本介绍
下图给出了Scan的例子，输入为1,2,3,4，运算符是+，输出结果是当前输入值与前面值的总和。
咋看貌似并不像是并行计算，但是Scan运算对于并行计算具有很大的作用。
下图给出了Scan的在实际生活中的例子，即银行存款账户余额情况，左边表示存钱，取钱数，右边表示余额。
2.2.2 Scan运算组成
Scan运算由输入向量，二元结合运算符(与Reduce运算类似)以及标识值组成。
基本上该课程中提到的运算符都需要具有Associative(结合性)，这样更加符合并行计算的特点。
为方便说明，标识值(Identify element)用I表示，它需要满足\(I op a = a\),其中op表示二元结合运算符，a表示任意值。举例来说可能更好理解：
例如如果运算符是and，那么标识值就是1。如果运算符是or，那么标识值就是0。
2.2.3 Exclusive Scan & Inclusive Scan
Scan运算分为两种：Exclusive和 Inclusive。区别是前者的输出不包含当前的输入值，后者则包括。下图给出了具体例子，使用的运算符是add。
下面以Inclusive Scan为例来计算Step复杂度和Work复杂度。
- Step复杂度
它的Steps复杂度和Reduce运算相同都是\(n(logn)\)，不再详细解释。
- Work复杂度
Work复杂度是\(o(n^2)\)。如下图示，第一个输出值的计算量是0，第二个是1，第三个是2，以此类推，最终复杂度是\(0+1+2+...+(n-1)≈\frac{n*(n-1+0)}{2}=o(n^2)\)
2.2.4 改进Scan运算的实现方法
如果采用上面的实现方式，Work复杂度太大，在实际运用中消耗太多资源，为了改善这一问题，出现了如 Hillis Steele和 Blelloch Scan方法，下面分别进行介绍。
- 1. Hillis Steele Inclusive Scan
还是从1到8做scan运算。该方法的主要实现思路是对于step i，每个位置的输出值是当前值加上其\(2^i\)左边的值。
Step 0：每个输出值是当前值加上前一个值。（绿色）
Step 1：每个输出值是当前值加上前两个值（蓝色）。
Step 2：每个输出值是当前值加上前4 (\(=2^2\)) 个值（红色）。
总结起来计算过程就是：
step 0: out[i] = in[i] + in[i-\(2^0\)];
step 1: out[i] = in[i] + in[i-\(2^1\)];
step 2: out[i] = in[i] + in[i-\(2^2\)];
通过下面完整的计算过程，我们可以明显看到Step复杂度依旧是\(logn\).
Work复杂度不能很明显的算出来，我们可以把所有的计算过程看成一个矩阵，纵向长度是\(logn\),横向长度是\(n\),所以Work复杂度是\(nlogn\)，该算法相对于之前的\(o(n^2)\)Work复杂度有了一定的提升。
- 2. Blelloch Scan
这是一种优化Work complexity的算法，比上面的要复杂一些。主要分为reduce和downsweep两步。
该算法由两部分组成：Reduce+Downsweep。
以exclusive scan加法运算为例，如下图示。
为更好说明用下图进行解释说明(图片来源：CUDA系列学习（五）GPU基础算法: Reduce, Scan, Histogram
如图示，上部分是Reduce，不再赘述。
下部分的Downsweep其实可以理解成Reduce的镜像操作，对于每一组输入in1 & in2有两个输出，左边输出out1 = in2，右边输出out2 = in1 op in2 （这里的op就是reduce部分的op），如图：
该算法有点绕，是否明白了呢？下面做个题来看看是不是真的明白了：
运算符是Max，将框中答案补充完整。
Ans：
介绍了这么多，那该方法的复杂度如何呢？
因为我们已经知道Reduce的Step复杂度是\(O(logn)\),Work复杂度是\(O(n)\)。
而Downsweep其实可以理解成Reduce的镜像运算，所以复杂度与之相同。
所以该算法整体的Step复杂度是\(O(2logn)\),Work复杂度是\(O(n)\)。这里对于Work复杂度存疑，不明白为什么不是\(O(2n)\)，但是网上好几篇博文都说是\(O(n)\)）
2.2.5 如何选择合适的算法
上面已经分别介绍了 Hillis Steele和Blelloch Scan算法，下面将二者的复杂度总结如下：
| Method | Step complexity | Work Complexity | type |
|---|---|---|---|
| Hillis Steele | \(O(logn)\) | \(O(nlogn)\) | SE (Step Effient) |
| Blelloch | \(O(2logn)\) | \(O(n)\) | WE (Work Efficient) |
在实际中可能会遇到如下图示的情况，随着工作量的变化我们需要动态改变算法策略。
例如最开始可能任务中地Work要远多于已有的处理器数量，这种情况执行速度受到了处理器数量的限制，此时我们则需要选择 更具高效工作性的算法(如Blelloch)。
随着工作量不断完成，渐渐地可用处理器数量比工作更多，此时我们则可以选择 更高效步骤的算法实现(如Hillis Steele算法)。
下面来做个题看是否已经理解了如何选择合适的算法。
| 情况 | Serial(串行) | Hillis Steele | Blelloch |
|---|---|---|---|
| 128k的vector， 1个processor | √ | ||
| 512个元素的vector， 512个processor | √ | ||
| 一百万的vector， 512个processor | √ |
解析：
- 第一种情况只有一个处理器，难不成还能选择并行算法？
- 第二种情况处理器数量刚好合适，工作量也不大，所以可以选择Hillis Steele
- 第三种情况工作量远大于处理器数量，所以选择Blelloch算法来提高工作效率。
2.3 Histogram
2.3.1 Histogram是啥
顾名思义，统计直方图就是将一个统计量在直方图中显示出来。
2.3.2 Histogram 的 Serial 实现
分两部分：1. 初始化，2. 统计
for(i = 0; i < bin.count; i++)
res[i] = 0;
for(i = 0; i<nElements; i++)
res[computeBin(i)] ++; // computeBin(i)函数是用来判断第i个元素属于哪一类2.3.3 Histogram 的 Parallel 实现
直接实现：
__global__ void naive_histo(int* d_bins, const int* d_in, const in BIN_COUNT){
int myID = threadIdx.x + blockDim.x * blockIdx.x;
int myItem = d_in[myID];
int myBin = myItem % BIN_COUNT;
d_bins[myBin]++;
}仔细看上面的代码，我们很容易看出d_bins[myBin]++;语句在并行计算过程中会出现read-modify-write（从全局内存中读取数据，修改数据，写回数据） 冲突问题，进而造成结果出错。
而serial implementation不会有这个问题，那么想实现parallel histogram计算有什么方法呢？
1. accumulate using atomics
下图给出了read-modify-write介绍，可以看出读取，修改，写入是3个独立的原子运算，但是如果我们将这3个操作整合成1个原子操作那么就可以很好地解决上述问题。而且现如今的GPU能够锁定特定的内存地址，因此其他的线程就无法访问该地址。
具体实现代码只需要将d_bins[myBin]++;修改成atomicAdd(&(d_bins[myBin]), 1);即可。
__global__ void simple_histo(int* d_bins, const int* d_in, const in BIN_COUNT){
int myID = threadIdx.x + blockDim.x * blockIdx.x;
int myItem = d_in[myID];
int myBin = myItem % BIN_COUNT;
atomicAdd(&(d_bins[myBin]), 1);
}但是对于atomics的方法而言，不管GPU多好，并行线程数都被限制到histogram个数N，也就是最多只有N个线程并行。
2. local memory + reduce
思路原理：设置n个并行线程，每个线程都有自己的local histogram（一个长为bin数的vector）；即每个local histogram都被一个thread顺序访问，所以这样没有shared memory，即便没有用atomics也不会出现read-modify-write问题。然后，我们将这n个histogram进行合并（即加和），可以通过reduce实现。
举例：
如下图示，假设有128个item，有8个线程可以使用，然后需要分成3类。
这样每个线程理论上需要处理16个item，即将这16个item分成3类。
通过使用并行线程，我们无需使用上面那个方法将read-modify-write原子化，因为这8个线程使用的是自己的local memory，而不是shared memory。如此一来我们只需要在所有线程计算结束后，再使用Reduce运算将8个线程的计算结果对应相加即可。
3. sort then reduce by key
该方法的主要思路是采用键值对(类似Python中的字典)的方式来记录数据，之后对键值对按键的大小排序得到下图中蓝色键值对。
如此一来键相同的挨在一起，之后便可用Reduce运算将相同类别的值相加即可。
三、作业应用
Tone Mapping(色调映射) 是转换一组颜色到另一组颜色的过程。之所以要色调映射是因为真实世界中的光谱强度是非常广的，例如白天太阳的光亮程度到晚上月亮的光亮程度，这之间的范围十分大。而我们的设备，如电脑显示屏，手机所能显示的光亮程度远小于真实世界，所以我们需要进行色调映射。我们手机上拍照功能中的HDR模式就是色调映射的作用。
该过程如果没有处理好，那么得到的效果则要么变得过于暗淡，要么过于明亮。
下图很好地展示了这两种情况。
下图左曝光较低，很多细节丢失。下图右虽然细节保留更多，但是窗户部分因为曝光过多基本模糊掉了。两图片下面的直方图也很直观的呈现了像素分布情况。
而色调映射就是为了解决上述情况，以期望达到如下图的效果。
本次作业需要结合所学的三种运算Reduce，Scan，Histogram。
MARSGGBO原创
2018-7-9
