【Udacity并行计算课程笔记】- Lesson 4 Fundamental GPU Algorithms (Applications of Sort and Scan)

2022-09-29 15:10:14

I. Scan应用——Compact

在介绍这节之前，首先给定一个情景方便理解，就是因为某种原因我们需要从扑克牌中选出方块的牌。

更formal一点的说法如下，输入是 \(s_0,s_1,...\), 我们提前预设条件来得到 Predicate，即每个元素都会根据条件输出True或False。然后我们根据Predicate(比如做与运算)就可以输出我们想要的值。

但是如下图示，我们的输出Output有两种表达形式:

第一种是 Sparse，即 \(s_0, - , s_2 , -, ...\);
第二种是Dense，即 \(s_0,s_2,s_4,...\)

以上两种表达方式实际计算过程如下图示：

对于Sparse而言，每次都需要判断当前卡片是不是方块的，如果是则执行computercard()操作。

而对于Dense而言，则需要先将方块卡片进行压缩compact(Cards, isDiamond())，然后再通过map操作对压缩后的方块卡片进行下一步的计算。

下面是两种方法所需要启动的线程数的示意图。

视频中老师对上图的说法是：

咋一看，Sparse需要启动的线程数更少，他总共用了52个线程，而Dense用了52+13个线程。但是需要注意的是其实Sparse有39个线程是闲置的，但是因为在log step 中空闲和非空闲线程都在运行，所以Sparse 其实付出了4倍于Dense的多线程成本。

我没太懂,疑问是Sparse不是有个if判断语句了吗，如果判断不是方块，这个线程不就结束了吗，而且Dense不还得压缩一下吗，而且也用了52个线程，为什么还说付出了更多线程呢？：

ex1:When to use Compact

上题想表达的就是Compact适用于有很多元素而且每个元素要执行的运算非常expensive。的时候，

ex2:Core Algorithm to Compact

假设我们有一组Predicate，我们希望输出这样一组数据，即输出True所属第几个，例如第一个T输出0，第二个是F，则输出 —，遍历到第二个T输出1，以此类推。

我们可以用什么运算方法实现呢？

思考几秒钟。

顶顶顶顶。。。是Scan。

ex3:Steps to Compact

下图介绍了Compact的步骤：

对输入的每个元素做出判定，得到Predicate。
创建一个大小与输入一致的的扫入数组(Scan-IN Array),对于每个元素，若判断其为真，则对应位置的扫如数组置为1，反之为0。
对扫入数组进行Exclusive-Sum-Scan，这个输出结果是压缩数组(Compacted Array)的分散地址(Scatter Address)。
根据得到的Address将输入scatter到输出。

下面是习题：

假设现在有数据是从1到一百万，有两个操作，一个是选择出能被17整除的数，另一个是选择出不能被31整除的数。问题是如何针对Predicate、Scan、Scatter三个运算而言，上面两个操作哪一个更快或者使用时间相同。

predicate, scan很好理解，因为都需要对所有元素遍历一遍，所以两个操作消耗的时间是一样的。
scatter其实也好理解，因为能被17整除的数远比不能被31整除的数少，所以显然前者消耗的时间更少。

ex4:Allocate

Compact is an operation that allocates exactly 1 item in the output for each true input and 0 items in the output for each false inpout.

下面举个栗子：假设有如下一组三角形

然后呢，有些三角形会出现在屏幕上，而有些不会。所以在我们把三角形传递到管道进行后续处理之前我们通常测试来判断每个三角形是否可见。这是一项压缩操作：我们把三角形输入流（包括可见的和不可见的）压缩到更小的输出流。其中输出流中的所有三角形都可见，即除了3,5外其他三角形都可见。

但是现在有个更复杂的问题就是假如有如下图的两个绿色的三角形怎么办？它们都只有一部分可见。为了将它们转化成三角形，我们需要结合边界对它们进行切割从而得到三角形。

以右边的绿色三角形为例，很简单，只需要用边界进行切割即可得到一个三角形。

下面的绿色三角形略微复杂一点点，其在屏幕内的边和边界构成了一个四边形，为了得到三角形可以如下图示的方法进行切割(当然也可以从左上角连接到右下角)。

上面两个绿色三角形在切割后可以分别得到1个和2个三角形，那么很自然地我们想知道最多可以得到几个三角形呢？

答案是5个，原因见下图。

possible allocate strategy

接着上面的内容可以知道一个可行的分配的方法是为每一个元素分配一个最大空间，以上面的三角形为例，此时最大空间为5。也就是说每一个输入三角形都需要创建一个长度为5个输出三角形的输出数组(中间数组)，接着对结果进行压缩即可。

但是！！！这样做有如下缺点：

浪费空间
需要扫描整个中间数组，而这个数组过于庞大。

那么怎么改进呢？

我们需要做的是为每个输入元素设计一张分配请求的列表，接着取回一个位置来写入请求。

举个栗子，假如输入表示三角形的个数，分别是[1,0,1,2,1,0,3,0],为了节省空间每个输出元素不再是有最大空间的数组，而是根据输入的大小动态分配空间。

例如第一个输入元素值是1，那么他的输出值则为0，表示index=0的空间需要分配给第一个输出元素。

第二个输入元素值是0，则其输出值为1

第三个输入元素值为1，因为第二个输入为0，所以其输出也为1

第四个输入元素值为2，同理它的输出需要往后挪一位，即为2

第五个输入元素值为1，而又因为前面的输入元素值为2，所以他需要挪两位，所以输出值为4

其他同理不再赘述。

仔细想想这是不是就是之前内容中的Exclusive Scan呢？

不仅是这个例子，Scan在GPU运算中还有很多应用，例如GPU快速排序中也许要用到Scan运算，所以Scan非常的重要。

Ex: Segmented Scan

在一些应用中，我们可能需要进行过很多次小扫描，而不是一次大扫描。当我们在GPU启动一个内核，我们一般希望在该内核中进行很多工作。所以就每个小扫描单独启动一个扫描内核不是很有意义，而且也浪费资源。相反我们可以把这些小扫描作为segment打包进一个大数组。然后利用一个特殊扫描运算符来单独扫描每一分段。

通常为了指出分段在数组中的起始位置，我们用第二个数组，该数组以1代表段首，0代表非段首。举个栗子：

下面我们需要对这个数组[1 2 3 4 5 6 7 8]进行Exclusive sum scan可以得到[0 1 3 6 10 15 21 28]。

那么segmented scan的方法是首先将原来的数组分段，假设这样分 [1 2 | 3 4 5 | 6 7 8]。

然后使用一个中间数组来记录分段的起始位置，即[1 0 1 0 0 1 0 0]

最后对每段在进行Exclusive sum scan运算，可得[0 1 | 0 3 7 | 0 6 13]

还是以上面的例子为例，如果现在需要做Inclusive Sum Scan,那此时的结果是什么呢？

结果如下图示：

SpMv (Sparse Matrix vector)

什么是稀疏矩阵

在实际应用中我们经常会遇到各种稀疏矩阵，即矩阵中很多元素值为0。其中最著名的稀疏矩阵当属PageRank。

PageRank是一个系数矩阵，用来统计全世界的网页之间的关联性。如下图示矩阵的行和列分别表示某一个网页。如果网页R和网页C有链接，那么对应位置上的值不为0。但是很显然实际上很多网页之间是没有关联的，所以最终这是一个非常庞大的稀疏矩阵。

压缩稀疏行, CSR

表示系数矩阵的传统方式叫做压缩稀疏行(Compressed Sparse Row, CSR)

下面举个栗子：
假设需要对如下两个矩阵做乘法运算，显然左边矩阵中包含有3个0元素，我们希望能够对该矩阵进行压缩从而达到节省空间和提高计算效率的目的。（为了方便说明所以只选取了简单的矩阵，所以可能压缩效果不太明显，但是当矩阵变得特别大的时候会非常有效。）

在CSR格式中，我们需要设置三个向量对矩阵进行压缩，分别是：

Value Vector(值向量)：用于存储非0值，左边的矩阵展开后得到向量 [a b c d e f]
Column Vector（列向量）:用来指示每个元素处于哪一列，即 [0 2 0 1 2 2]
Rowptr（行指针）: 注意这里不再为每一个元素标明所在行号了，仔细观察可以知道上面的列向量中的元素值是从小到大排列，如果后面一个元素值小于前面一个说明是新的一行了。所以只需要指明从哪一个元素开始表示新的一行即可。例如第一行开始元素是a，它在value vector的位置是0。第二行开始的元素是c，它在value vector的索引是2，同理f索引是5。所以最终的Rowptr向量是 [0 2 5]

下面做个练习题看看你做对了吗：

如何应用CSR？

有了CSR格式向量后，如何应用到矩阵相乘呢？

下图给出了详细的步骤

1.首先要将值向量和行指针向量共同创建一个值向量的分段表示，也就是说每一段表示稀疏矩阵的一行，即得到 [ a b | c d e | f ]
2.结合列向量索引值得到需要相乘的向量的索引。例如a的列索引是0，那么对应的与之相乘的元素的行索引也应为0，这样就可以找到是x。同理列索引为2的b对应行索引为2的z。其他同理，不再赘述。因此可得到如下向量值 [ x z | x y z | y ]
3.对应元素相乘： [ ax bz | cx dy ez | fy ]
4.最后使用前面介绍的Exclusive sum segmented scan便可求出结果。

II.Sort

排序在GPU应用中有不少挑战，大多数的算法都是串行的，或者说通常以串行方式体现。很多我们在学校学到的算法在此系列课程中可能并不适用，这在以后的内容中会体现出来。所以接下来我们会找一些高效地并行算法，这些算法有如下特性：

keep hardware busy
limit branch divergence
prefer coalesced memory access

1. 冒泡排序

下面举个栗子：

对 [5 1 4 2 3]使用冒泡排序：

我们都知道串行方式的冒泡排序是每次都需要比较相邻的元素。如果第一个比第二个大，就交换他们两个。重复这个过程直到完成排序，所以时间复杂度是O(n^2)。

那么如果以并行方式的话是怎么做呢？下图给出了示例：

第一次遍历：每两个元素组成一组进行比较，如果前者比后者大，则两者交换位置；例如 5 1组合， 4 2组合，3无法组合，可以暂时不管。最终比较置换后得到 [ 1 5 2 4 3]
第二次遍历：此时不能再像上一次遍历一样组队了，此时组队需要往后挪一位。也就是说5 2组合， 4 3组合，1暂时不管。同理得到[ 1 2 5 3 4 ]
第三次遍历：同理，此时组队方式需要往前挪一位，即1 2组合，5 3组合，4暂时不管。最终得到 [ 1 2 3 5 4]
...
最终得到 [ 1 2 3 4 5 ]

奇偶排序(odd and even sort)

那么以并行方式运算的冒泡排序的效率如何呢？

其实上图也可以称为奇偶排序。例如第一行从0开始配对，叫做偶数排序。同理第二行叫做奇数排序。

之前介绍过并行计算评估标准有Step和Work，所以下面计算这两个标准复杂度。

Step: 由上图中的红色斜线可以看到，数字5从最左依次移动到最右，也就是说遍历了整个数组，所以该冒泡并行排序的Step复杂度是n。
Work：step复杂度确定之后，work复杂度就好理解了。可以看到每一个step都需要n/2的work（作比较），然后step复杂度又等于n/2,所以work复杂度为n^2。

2. 归并排序(merge sort)

1) 方法回顾

下图展示了传统的归并排序：

2) 并行方法复杂度

上图中的需要排序的元素数量很少，如果数量达到一百万会怎么样呢？由下图可以看到如要使用归并排序，首先需要将1M的数据分成两半，即500K+500K。然后再重复分半，直到最后得到1M单独的元素。很明显Step复杂度为O(log(n)),work复杂度为O(nlog(n))。

仔细观察，上面的计算其实可以划分成如下3个阶段：

第一阶段就是下图中蓝色字体部分，该部分有大量运算量较小的任务组成，每个任务分配一个thread。
第二阶段是黑色字体部分，该部分的任务计算量中等，每个任务分配一个Block用于计算
第三阶段就是绿色字体部分。此时只有一个计算量很大的任务。

3) 如何比较置换？

再仔细回顾一下归并排序，其主要思想是以大化小，然后拼接排序。如下图示，假设已经得到了两个有序数组。接下来要做的是就是分别比较两个数组中的第一个元素，然后输出较小的元素。迭代这一过程直到完成排序。

但是上面的方式并不适用于GPU并行计算，所以怎么办呢？此时需要借鉴上面的内容:

如下图示，通过scater运算可以得到每个元素指定的输出索引。例如输入数字5的输出索引为3,21的输出索引为5。

上面的是scater这一思想可以借鉴过来解决归并排序中的比较排序问题。如下图示有两个已经排好序的输入数组，分别是

[1 3 12 28]
[2 10 15 21]

我们可以通过得到每个元素对应的输出索引，再通过索引值只需要进行拼接即可完成排序。

但是你是不是不禁想问这些索引值又是怎么来的呢？为方便说明，以输入数组1中的数字12为例。

首先数字12可以通过它的thread id很容易地知道自己在input list1中的位置，没错是2(索引从0开始，你也可以理解成其前面有几个元素)。

但是怎么知道在input list2中的位置呢？因为两个list都是有序的，所以可以使用二叉树排序方法(复杂度为O(log(n)))求出它在list2中的相对索引位置，即2。所以最终12的输出索引值应该是2+2=4

4) SM闲置问题

有个问题在前面没有提到，就是归并排序的第三阶段许多SM会开始闲置，因为此时我们只有一项任务。

所以为了解决这一问题，我们希望设计一个算法，可以允许很多SM进行单个归并工作。怎么做呢？还是看下面的栗子：

假设有两个大且长的数组需要归并，为了让多个SM同时工作，我们可以每隔n个元素进行切割，这样就可以得到若干个数组子集，我们把这些我们挑选的元素叫做分解器(splitters)。示意图如下：

上图中的A B C D ...就是分解器，它们是单个元素，每个分解器之间相隔一定数量的元素。可以看到上图中画了很多虚线箭头，它们表示各个分解器在另一个数组的相对位置。通过相对位置，我们能够把分解器归并成一个有序数组:[ E A B F C G D H ]。

有了这个有序数组，接下来要做的是就是将分解器之间的元素进行排序，因为A B之间的元素和E F之间的元素的相对大小不确定。 这句话什么意思呢？看个例子就明白了：

假设

A B段数据为: [ 5 8 12 14 17 ], 即A=5,B=17.
E F段数据为: [ 1 4 9 15 20 ],即E=1, F=20.

我们可以通过二叉树排序确定A和B在E F段的位置，但是A B之间的数据，即8 12 14的还需要和9 15进行归并排序。当然也有可能是C D的那种情况，即C位于F G之间，而D位于G H之间，所以不同分解器之间的归并排序计算量不一样,也就意味着不同的SM的任务量也不一样。但是总的来说，通过这样分解我们能把很多空闲的SM都利用起来了。

3. 双调排序(Bitonic Sort)

不同于上面的排序方法，双调排序是一种与数据无关的排序方法。该算法特别适用于GPU并行计算。

在介绍双调排序之间需要先介绍什么是双调序列。双调序列是指先单调递增后单调递减或先单调递减后单调递增的序列。例如[ 1 5 9 6 4]和[9 6 4 5 8]都是双调序列，而[1 4 2 3]则不是双调序列。所以双调序列其实就是两个方向不同的单调(monotonic)序列的组合。

双调序列的严格定义如下：

定义:一个序列a1,a2,…,an是双调序列(Bitonic Sequence)，如果：
(1)存在一个ak(1≤k≤n), 使得a1≥…≥ak≤…≤an成立；或者
(2)序列能够循环移位满足条件(1)

双调归并网络是基于Batcher定理而构建的。Batcher定理是说

将任意一个长为2n的双调序列A分为等长的两半X和Y，将X中的元素与Y中的元素一一按原序比较，即a[i]与ai+n比较，将较大者放入MAX序列，较小者放入MIN序列。则得到的MAX和MIN序列仍然是双调序列，并且MAX序列中的任意一个元素不小于MIN序列中的任意一个元素。

说了这么多看下面例子：

可以看到序列 [ 4 2 1 3 ]是一个双调序列，绿色连线表示比较器，通过若干次比较后得到了有序序列。

更多的细节可以阅读双调排序Bitonic Sort，适合并行计算的排序算法。

4. Radix Sort(基数排序)

基数排序（英语：Radix sort）是一种非比较型整数排序算法，其原理是将整数按位数切割成不同的数字，然后按每个位数分别比较。它是这样实现的：将所有待比较数值（正整数）统一为同样的数字长度，数字较短的数前面补零。然后，从最低位开始，依次进行一次排序。这样从最低位排序一直到最高位排序完成以后，数列就变成一个有序序列。

基数排序的方式可以采用LSD（Least significant digital）或MSD（Most significant digital），LSD的排序方式由键值的最右边开始，而MSD则相反，由键值的最左边开始。

基数排序的时间复杂度是 \({\displaystyle O(k\cdot n)}\)，其中 {\displaystyle n} n是排序元素个数， k是数字位数。注意这不是说这个时间复杂度一定优于 \({\displaystyle O\left(n\cdot \log \left(n\right)\right)}\)， k的大小取决于数字位的选择（比如比特位数），和待排序数据所属数据类型的全集的大小；k决定了进行多少轮处理，而n是每轮处理的操作数目。

看一下示例就明白了(下面例子参考基数排序（Radix Sort）)：

以LSD为例，假设原来有一串数值如下所示：

73, 22, 93, 43, 55, 14, 28, 65, 39, 81

首先根据个位数的数值，在走访数值时将它们分配至编号0到9的桶子中：

分配过程：

接下来将这些桶子中的数值重新串接起来，成为以下的数列：
收集过程：

81, 22, 73, 93, 43, 14, 55, 65, 28, 39

接着再进行一次分配，这次是根据十位数来分配：
分配过程：