基于CUDA的高效排序算法（Sort on GPU）

一、前 言

排序（sort），是计算机内经常进行的一种操作，其目的是将一组“无序”的记录序列调整为“有序”的记录序列，是各领域中广泛应用的核心算法之一。受限于算法本身的限制（所有基于比较的排序算法复杂度下限为O(n )，对于一组超大量数据单纯的cpu排序所需时间过长，通常1000万的数据排序就需要1秒左右，更大量的数据排序则需要分钟级的运行时间。因此，可以通过并行排序操作来进一步加速排序算法。

二、基于cpu的排序算法

经过数十年的发展，目前的排序算法已有基数排序、冒泡排序、选择排序、插入排序、希尔排序、堆排序、归并排序、快速排序等。显然，我们希望在已有的高效排序算法上进行进一步的优化，因些我们忽略冒泡排序、选择排序、希尔排序这一类指数级复杂度的算法，也忽略基数排序这一类不适于于浮点型变量的算法，而在堆排序、归并排序、快速排序这些O(n )算法中，直观的我们能发现归并排序、快速排序都是基于二分的策略，所以可以自然地同时对左右两半同时进行排序；而归并的过程也可以通过直接查找每个元素在合并数组中的位置来契合并行化，但快速排序中划分的过程无法并行实现（至少我没有想到，欢迎各路大神提出好的思路）；因此归并排序是最适于并行化的排序算法。

然而，这些已有的高效排序算法通常只在数据量较大时才有着更高的效率，在小数据时反而常数更小的选择排序、插入排序等有着更出色的性能；如 中的std::sort中在长度小于16或32时不再继续划分，而是调用插入排序来实现；因此在基于cuda的排序中也需要在小数据上选择另一种算法以达到更出色的性能。在此，我们选择并不常见、但复杂度O(n n)相对较低的双调排序来实现，自然是因为它能很好的并行化。

因此我们选择双调排序和归并排序来实现CUDA上的排序操作，如需详细了解这两个算法请查阅相关资料。

三、基于cuda的排序算法的实现

1. 双调排序的并行化

我们暂定一个线程块的线程数为 256（tile_size），因为双调排序中每次需要进行一次比较来判断是否需要交换，因此每一个线程块处理连续的tile_size x 2个数据（以下称为一个数据段），每个线程负责一次比较操作。因此该核函数线程块大小为tile_size，网格大小为总共的数据段数量。

#define tile_size 256u
#define log_tile_size 8u

unsigned int b = n + (tile_size <> log_tile_size + 1, s = tile_size > 1>>>(data, n);

首先先将需要的数据载入共享内存中，由于每一个线程块处理tile_size x 2个数据，所以一个线程块负责的数据段起始为“线程块索引 x tile_size x 2”；每一个线程将相隔tile_size的两个数据载入共享内存中，超过数据规模的部分直接填充浮点数的最大值，这样一个线程束访问一段连续的地址以达到更高的访存效率；还需添加一个线程块内同步函数，确保所有数据都已载入共享内存后再开始排序。

__shared__  double buffer[tile_size <
unsigned int d = blockIdx.x buffer[i] = i buffer[i + blockDim.x] = i + blockDim.x __syncthreads();

然后需要将随机序列转换为双调序列，在这个过程中需要交换操作，因此需要先定义一个设备端的交换函数，该交换函数可定义为模板函数以支持不同类型的数据交换。

    
// auxiliary functions
templatetype>__device__ void swap(type& v1, type& v2)
{
 type t = v1;
 v1 = v2;
 v2 = t;
}

unsigned int s, t;
for (unsigned int k = 2; k <= blockDim.x; k <<= 1)
 for (unsigned int p = k >> 1; p; p >>= 1)
 {
  s = (i & -p) <  t = s | p;

  if (s & k ? buffer[s]  buffer[t])
   swap(buffer[s], buffer[t]);
  __syncthreads();
 }

该过程转化为代码的操作与上述合并过程中的内循环完全一致。

for (unsigned int p = blockDim.x; p; p >>= 1)
{
 s = (i & -p) < t = s | p;
 if (buffer[s] > buffer[t])
  swap(buffer[s], buffer[t]);
 __syncthreads();
}

最后将该数据段的有效排序结果写回全局内存。

if (i  data[i] = buffer[i];
if (i + blockDim.x  data[i + blockDim.x] = buffer[i + blockDim.x];

最终分块双调排序的代码如下：

#include 

// binotic sort
__global__ void binoticsort_gpu(double* data, unsigned int n)
{
 __shared__  double buffer[tile_size <
 unsigned int d = blockIdx.x < data += d; n -= d;
 buffer[i] = i  buffer[i + blockDim.x] = i + blockDim.x  __syncthreads();

 unsigned int s, t;
 for (unsigned int k = 2; k <= blockDim.x; k <<= 1)
  for (unsigned int p = k >> 1; p; p >>= 1)
  {
   s = (i & -p) <   t = s | p;

   if (s & k ? buffer[s] > buffer[t] : buffer[s]     swap(buffer[s], buffer[t]);
   __syncthreads();
  }
 for (unsigned int p = blockDim.x; p; p >>= 1)
 {
  s = (i & -p) <  t = s | p;
  if (buffer[s] > buffer[t])
   swap(buffer[s], buffer[t]);
  __syncthreads();
 }

 if (i   data[i] = buffer[i];
 if (i + blockDim.x   data[i + blockDim.x] = buffer[i + blockDim.x];
}

经过该分块双调排序后，规模为n的数据被划分为多个长为 tile_size x 2的数据段，每个数据段内的 tile_size x 2个数据已经被成功排序，下面只需要逐步将这些有序数据段合并即可。

2. 归并操作的并行化

有序数据段合并操作即像归并排序中的归并一样，但是需要对其做一定的修改以更好地适用于并行化。通常的合并操作需从头逐个比较两个数组中的元素，取较小者置于结果数组中——但是这个循环是具有数据依赖性的，故这种相互之间不独立的操作无法进行并行化处理。因此我们需要用另一种方法来进行合并操作，我们直接计算两个待合并数组中每个元素在结果数组中应当存在的位置，然后把该元素置于该位置即可。

显然对于其中一待合并数组中比某元素小的元素个数必然是该元素的索引，因此只需找到另一待合并数组中小于（或小于等于）该元素的元素个数即可，注意这儿必须要一个数组用小于而另一个数组用小于等于，否则两个相同的数字就会落在同一位置。在此直接引用 中的std::lower_bound和std::upper_bound中的算法来实现，只是需要将其编译为设备端函数。

// auxiliary functions
templatetype>__device__ unsigned int lower_bound(const type* data, unsigned int n, const type& value)
{
 unsigned int len = n, half;
 const type* start = data, * mid;
 while (len)
 {
  half = len >> 1;
  mid = data + half;
  if (*mid   {
   data = mid + 1;
   len = len - half - 1;
  }
  else
   len = half;
 }
 return data - start;
}
templatetype>__device__ unsigned int upper_bound(const type* data, unsigned int n, const type& value)
{
 unsigned int len = n, half;
 const type* start = data, *mid;
 while (len)
 {
  half = len >> 1;
  mid = data + half;
  if (value    len = half;
  else
  {
   data = mid + 1;
   len = len - half - 1;
  }
 }
 return data - start;
}

mergedirect_gpu<<> log_tile_size, std::min(b, 32768u), b + 32767 >> 15), tile_size>>>(data, n, s, tmp);

// merge directly
__global__ void mergedirect_gpu(const double* data, unsigned int n, unsigned int s, double* result)
{
 unsigned int d1 = (blockIdx.z * gridDim.y + blockIdx.y) * s < if (i   result[d2  if (s + i   result[i + upper_bound(data + d1, s, data[s + i])] = data[s + i];
}

然后将这两个函数整合成一个排序函数即可，该函数接受两个参数：指向数据首地址的指针和数据规模，返回指向已排序数据首地址的指针，其中的指针均指向设备端。该函数首先需要申请一块与输入数据等长的临时数组用于存放合并的结果，然后先进行一次双调排序将输入数据分块排序为多个长为 tile_size x 2的有序数据段，再来一个循环将这些有序数据段逐轮合并直至数据段长度大于数据规模，最后释放临时数组空间并返回结果。

double* sortdirect_gpu(double* data, unsigned int n)
{
 double* tmp;
 cudaMalloc(&tmp, sizeof(double) * n);
 unsigned int b = n + (tile_size <> log_tile_size + 1, s = tile_size < binoticsort_gpu<<> 1>>>(data, n);
 for (b = b + 1 >> 1; s > 1)
 {
  mergedirect_gpu<<> log_tile_size, std::min(b, 32768u), b + 32767 >> 15), tile_size>>>(data, n, s, tmp);
  std::swap(data, tmp);
 }
 cudaDeviceSynchronize();
 cudaFree(tmp);

 return data;
}

至此，我们已初步实现了基于CUDA的排序算法，我们选择了不同的数据规模在RTX A6000上对它的性能进行测试，测试代码如下：

// test.cu
#include 
#include 
#include 
#include 
#include 
#include 

int main(int argc, char** argv)
{
 cudaSetDevice(0);
 cudaFree(0);

 unsigned int test[] = { 1025u, 65537u, 1048577u, 16777217u, 134217729u, 1073741825u };
 const unsigned int runtime = 3;

 std::chrono::high_resolution_clock::time_point start, end;
 unsigned long long elasped[3] = { 0ull };

 for (unsigned int t = 0; t test) / sizeof(unsigned int); ++t)
 {
  unsigned int n = test[t];
  double* original = new double[n];
  double* data = new double[n];

  elasped[0] = elasped[1] = elasped[2] = 0ull;
  srand(time(nullptr));
  for (int i = 0; i   {
   for (unsigned int i = 0; i     original[i] = data[i] = (rand() * 2.0 / RAND_MAX - 1.0) * rand();

   // cpu sort 1
   start = std::chrono::high_resolution_clock::now();
   qsort(data, n, sizeof(double), [](const void* left, const void* right)->int {double tmp = *(double*)left - *(double*)right; if (tmp > 0) return 1; else if (tmp return -1; else return 0; });
   end = std::chrono::high_resolution_clock::now();
   elasped[0] += std::chrono::duration_cast<:chrono::nanoseconds>(end - start).count();

   // cpu sort 2
   memcpy(data, original, sizeof(double) * n);
   start = std::chrono::high_resolution_clock::now();
   std::sort(data, data + n);
   end = std::chrono::high_resolution_clock::now();
   elasped[1] += std::chrono::duration_cast<:chrono::nanoseconds>(end - start).count();

   // gpu sort
   double* pd_data;
   cudaMalloc(&pd_data, sizeof(double) * n);
   cudaMemcpy(pd_data, original, sizeof(double) * n, cudaMemcpyHostToDevice);
   start = std::chrono::high_resolution_clock::now();
   pd_data = sortdirect_gpu(pd_data, n);
   end = std::chrono::high_resolution_clock::now();
   elasped[2] += std::chrono::duration_cast<:chrono::nanoseconds>(end - start).count();
   cudaMemcpy(data, pd_data, sizeof(double) * n, cudaMemcpyDeviceToHost);
   cudaFree(pd_data);
  }

  printf("data number: %u\n", n);
  printf("cpu merge sort cost %.5f ms\n", elasped[0] / 1e6 / runtime);
  printf("cpu quick sort cost %.5f ms\n", elasped[1] / 1e6 / runtime);
  printf("gpu sort cost %.5f ms\n", elasped[2] / 1e6 / runtime);
  printf("-------------------------------------------\n");

  delete[] data;
  delete[] original;
 }

 return 0;
}

我们先不加任何优化进行编译

nvcc -o ./benchmark -g ./test.cu && ./benchmark