GPU版基数排序（radix sort）

基数排序算法的思想很简单，但如何提高并行性，使得能用GPU进行高效的排序计算呢？本文介绍适合于GPU计算版本的排序算法。

---

为了说明GPU版的排序算法，我们先复习一下最初版本的排序算法：

基数排序算法

A picture is worth than a thousand words.

因此我们上图说话。

基数排序是从最低位看起，把这一位数字值相同的按照扫描顺序放入同一个桶里面，值小的桶在前面。当所有数字都扫描完，再使用高一位，循环上述步骤，直至达到所有数字所有的最高位数，最后输出的就是排序后的答案。

CPU版本迭代基数排序算法

input: unsigned int inputVals[n]
output: unsigned int outputVals[n]

for(i=0; i< numBins; i+= numBits)    1. 计算inputVal中第i位值为0到numBins-1的各有多少个，其值分别放入binHistogram[numBins]中    2. 对binHistogram[numBins]进行scan计算，其值分别放入binScan[numBins]中，即    binScan[k]=binHistogram[0]+···+binHistogram[k-1]    3. 遍历inputVal数组，计算第j个量时，计算其值所属的桶bin的值      outputVals[binScan[bin]] = vals_src[j];      binScan[bin]++;       4. swap(inputVals,outputVals)

• 其中numBins是数字由多少位表示的，比如int就是32个bit位表示。
• numBits表示一次考虑多少位，比如可=1,2,…
• binHistogram[k] 记录了第k个桶中的数目
• binScan[k]记录了第k个桶中的应该起始的位置

GPU版本基数排序算法

GPU版本的排序算法与CPU在前两步都是类似的。

for(i=0; i< numBins; i+= numBits)    1. 同CPU    2. 同CPU    3. 决定每个数值的相对位置，比如：[0 0 1 1 0 0 1]                             -> [0 1 0 1 2 3 2]    4. 根据第二步和第三步计算的结果，决定每个每个数字的绝对位置，通过map运算从相应的位置取值

那些需要避免的GPU实现的坑

终于到了我最想提醒大家的部分了。

1. 如何决定每个数值的相对位置？

决定每个数值的相对位置其实使用的是compact方法，以input为[0 0 1 1 0 0 1]为例子，计算如下：

• 计算值为0的相对位置
  计算针对0的Predicate: [true true false false true true false]
  其中true=1, false=0: [1 1 0 0 1 1 0]
  exclusive scan: [0 1 2 2 2 3 4]
  对于predicate为true的值，留下上一步的结果[0 1 Nan Nan 2 3 Nan]
• 计算值为1的相对位置
  计算针对1的Predicate:[false false true true false false true]
  其中true=1, false=0: [0 0 1 1 0 0 1]
  exclusive scan: [0 0 0 1 2 2 2]
  对于predicate为true的值，留下上一步的结果[0 1 0 1 2 3 2]
  ……

其中scan使用Hill/Steele算法。

2. 但用以上的方法能算出真正的相对位置么？

不能。
What？不能你还写这么久？
其实也不是完全不能，当排序的数目小于block中最大threads数目时（比如我设置的是2048），是可以正确排序的。为什么呢？因为我们使用的scan方法中需要同步。
先复习Hill/Steele算法，以加法为例，对[1, 2, 3, 4, 5, 6, 7, 8]进行scan计算

平行计算中最怕的就是没能同步啦，那么再上Hill/Steele算法的代码。

// Scan__global__ void scan_cal(unsigned int *d_cdf, const unsigned int *d_bins, const size_t size){       // if the thread exceed the array size, return it    int myId = blockDim.x*blockIdx.x + threadIdx.x;    if (myId >= size)        return;    int current_value = d_bins[myId];    int current_cdf = d_bins[myId];         for (unsigned int interval = 1 ; interval < blockDim.x; interval<<= 1)    {        if (tid >= interval)            current_cdf += d_bins[myId - interval];        __syncthreads();        d_bins[myId] = current_cdf;        __syncthreads();    }    // exclusive scan    d_cdf[myId]= d_bins[myId]-current_value;    __syncthreads();} 

发现问题了么？在图中每一行计算完后，需要等整个一行同步后才能进行下一次的计算。而__syncthreads();函数是block-level的同步，即在不同的block中的计算并不会同步后再进行，因而一些block没有算完的时候，其他block已经取了之前的值接着下面的计算了。因此如果只有一个block，即排序的数目小于block中最大threads数目时，还是可以正确排序的。

3. 那，怎么办呢？

这里提供一个可以排序2²⁰的方法。我们对之前的GPU算法进行小的修改：
blockSize是一个block里面有多少个线程（blockDim.x×blockDim.y×blockDim.z），对于cuda计算能力1.x的，blockSize最大值为1024；1.x+的值为2048。

for(i=0; i< numBins; i+= numBits)    1. 计算inputVal的每blockSize个sub-array中，第i位值为0到numBins-1的各有多少个，其值分别放入binHistogram[numBins×gridDim]中    2. 对binHistogram[numBins×gridDim]进行scan计算，其值分别放入binScan[numBins×gridDim]中，即    binScan[k]=binHistogram[0]+···+binHistogram[k-1]    3. block内决定每个数值的相对位置，比如：[0 0 1 1 0 0 1]                                    -> [0 1 0 1 2 3 2]    4. （不变）根据第二步和第三步计算的结果，决定每个每个数字的绝对位置，通过map运算从相应的位置取值

即分组compact。还是以[0 0 1 1 0 0 1]为例，此时我们设置blockSize=4，numBins=2（0和1），gridDim=ceil(7/blockSize)=2：

1. binHistogram=[ 2 2 | 2 1]
2. binScan=[0 2 | 4 6]，其中：
   
   • 第一个0代表第一组中值为0的从地址0开始
   • 第二个2代表第二组中值为0的从地址2开始
   • 第三个4代表第一组中值为1的从地址4开始
   • 第四个6代表第二组中值为1的从地址6开始
3. 1. 计算值为0的相对位置
      计算针对0的Predicate: [true true false false | true true false]
      其中true=1, false=0: [1 1 0 0 | 1 1 0]
      exclusive scan: [0 1 2 2 | 0 1 2]
      对于predicate为true的值，留下上一步的结果[0 1 Nan Nan | 0 1 Nan]
   2. 计算值为1的相对位置
      计算针对1的Predicate:[false false true true | false false true]
      其中true=1, false=0: [0 0 1 1 | 0 0 1]
      exclusive scan: [0 0 0 1 | 0 0 0]
      对于predicate为true的值，留下上一步的结果[0 1 0 1 | 0 1 1]
4. 绝对位置=相对位置+binScan[bin的位置]
   绝对位置=[ (0+0) (1+0) (0+4) (1+4) | (0+2) (1+2) (1+6) ]=[0 1 4 5 | 2 3 7]

求绝对地址的代码如下：

__global__ void specific_address(unsigned int *d_out_value, const unsigned int *const d_in_value,    const unsigned int *d_cdfs, const unsigned int *d_relative_pos,    const size_t numBins, const unsigned int pos, const size_t size){    int myId = blockDim.x*blockIdx.x + threadIdx.x;    if (myId >= size)        return;    unsigned int mask = (numBins - 1) << pos;    int bin = (d_in_value[myId] & mask) >> pos;    int real_bin_addr = bin*gridDim.x + blockIdx.x;    unsigned int address = d_cdfs[real_bin_addr] + d_relative_pos[myId];    d_out_value[address] = d_in_value[myId];    __syncthreads();}

还要解决的问题

之前提到的，上述算法对最多2²⁰数字进行排序（对于计算能力1.x以上的显卡可以为2²²）。因为还是scan同步那里进行了限制，因而还有修改的空间。因为我所做的工作计算量小于2²⁰因此就没接着继续研究了，有兴趣的同学欢迎分享对更高数据集排序实现的idea，我们一起学习~