简介
排序是任何应用的基本构造块的关键算法之一。有许多排序算法已经被广泛研究。，常见的排序算法时间和空间复杂度如下：
一些排序算法属于分治算法的范畴。这些算法适用于并行性，并适合 GPU 等架构，在这些架构中，要排序的数据可以被划分进行排序。快速排序就是这样一种算法。如前所述，Quicksort 属于“分而治之”的范畴。这是一个三步走的方法，如下所示:
- 从需要排序的数组中选择一个元素。该元素充当枢轴元素。
- 第二步是划分所有元素的位置。小于枢轴的所有元素都向左移动，大于或等于枢轴的所有元素都向右移动。这一步也称为分区。
- 递归地执行步骤 1 和 2，直到所有的子数组都被排序。
3 种快排枢轴元素选择方法：
- 随机（rand函数）：几乎不会出现最坏情况的复杂性
- 固定（队首、队尾）
- 三数取中（队首、队中和队尾的中间数）
与其他排序算法(如需要额外存储的合并排序)相比，快速排序的内存负载和要求更低。快速排序的更实际的实现使用随机化版本。随机化版本的预期时间复杂度为 O ( n l o g n ) O(nlogn) O(nlogn)。在随机化的版本中，最坏情况的复杂性也是可能的，但是对于特定的模式(如排序数组)不会出现这种情况，随机化的快速排序在实践中效果很好。
cpu版本的快速排序实现：
int partition(std::vector<int> &v, int low, int high) {
int pivot = v[low]; //记录第一个枢轴
while (low < high) {
while (low < high && v[high] >= pivot) {
high--; // 找到第一个比pivot小的
}
v[low] = v[high]; // high指示小于pivot，交换
while (low < high && v[low] <= pivot) {
low++; // 找到第一个比pivot大的
}
v[high] = v[low]; // low指示大于pivot，交换
}
// 记录新的枢轴,一次partition后 pivot左边比pivot小。后边比pivot大
v[low] = pivot;
return low;
}
void quick_sort(std::vector<int> &v, int low, int high) {
if (low < high) {
int i = partition(v, low, high);
quick_sort(v, low, i - 1);
quick_sort(v, i + 1, high);
}
}
在本文中，我们将利用动态并行（从 CUDA 6.0 和具有 3.5 架构的图形处理器开始引入的），在 GPU 上高效地实现快速排序。
CUDA动态并行和快速排序
CUDA动态并行
过去，在需要递归等特性的 GPU 上高效实现快速排序等算法非常困难。随着 GPU 架构 3.5 和 CUDA 5.0 的发展，引入了一个新特性，称为动态并行。
以前的快速排序算法要求递归启动内核。通过 CPU 调用内核一次，内核完成执行后，我们返回到 CPU 线程，然后重新启动它。这样做的结果是将控制权交还给中央处理器，还可能导致中央处理器和图形处理器之间的数据传输，这是一项昂贵的操作。
动态并行允许内核中的线程从图形处理器启动新内核，而无需将控制权交还给中央处理器。动态这个词来自于它是基于运行时数据的事实。线程可以同时在GPU设备端启动多个GPU内核，可以让递归算法更加清晰易懂，也更加容易理解。下图简化了解释:
如果我们将这个概念转化为快速排序的执行方式，它看起来会是这样的:
深度 0 是来自 CPU 的调用。对于每个子阵列，我们启动两个内核:一个用于左阵列，一个用于右阵列。达到内核的最大深度或元素数量小于 32(即扭曲大小)后，递归停止。为了使内核的启动处于非零流和异步状态，以便子阵列内核独立启动，我们需要在每次内核启动之前创建一个流:
cudaStream_t s;
cudaStreamCreateWithFlags( &s, cudaStreamNonBlocking );
cdp_simple_quicksort<<< 1, 1, 0, s >>>(data, left, nright, depth+1);
cudaStreamDestroy( s );
这是非常重要的一步，因为，否则，内核的启动可能会被序列化。
CUDA快速排序
对于我们的快速排序实现，我们将利用动态并行性递归启动 GPU 内核。实施快速排序的主要步骤如下:
- CPU 启动第一个内核:内核一个块一个线程启动。左边的元素是数组的开始，而右边是数组的最后一个元素(基本上是整个数组):
int main(int argc, char **argv)
{ ...
cdp_simple_quicksort<<< 1, 1 >>>(data, left, right, 0);
}
- 极限检查:在从内核内部启动内核之前，检查两个标准。首先，检查我们是否已经达到硬件允许的最大深度限制。其次，我们需要检查子数组中要排序的元素数量是否小于扭曲大小(32)。如果其中一个是真的，那么我们必须按顺序进行简单选择排序，而不是启动新的内核:
__global__ void cdp_simple_quicksort( unsigned int *data, int left, int right, int depth )
{ ...
if( depth >= MAX_DEPTH || right-left <= INSERTION_SORT )
{
selection_sort( data, left, right );
return;
}
- 划分:如果满足前面的条件，那么将数组划分为两个子数组，并启动两个新的内核，一个用于左数组，另一个用于右数组。并从从内核内部启动一个内核:
__global__ void cdp_simple_quicksort( unsigned int *data, int left, int right, int depth ) {
...
while(lptr <= rptr)
{
// Launch a new block to sort the left part.
if(left < (rptr-data))
{ // Create a new stream for the eft sub array
cdp_simple_quicksort<<< 1, 1, 0, s >>>(data, left, nright, depth+1);
}
// Launch a new block to sort the right part.
if((lptr-data) < right)
{//Create stream for the right sub array
cdp_simple_quicksort<<< 1, 1, 0, s1 >>>(data, nleft, right, depth+1);
}
}
- 执行代码:使用以下命令用
nvcc编译器编译您的应用:
$nvcc -o quick_sort --gpu-architecture=sm_70 -rdc=true quick_sort.cu
nvcc编译分成device部分编译和host部分编译
- host部分直接调用平台编译器进行编译Linux使用gcc，window使用cl.exe，
device部分编译，此部分编译分两个阶段：
- 第一阶段将源文件.cu文件的device部分编译成ptx文本指令
- 第二阶段将ptx文本指令编译成在真实架构上运行的二进制指令，第二阶段可能发生在生成可执行程序的过程中，也可能发生在运行可执行程序的过程中（just-in-time compilation）。
在生成可执行程序的过程中可以根据nvcc选项选择是否将ptx文本指令(x.ptx中间文件中)、二进制指令（x.cubin中间文件）嵌入到可执行程序中，一般有3种嵌入方式：
- 只嵌入x.ptx（第二阶段被忽略，全部依赖just-in-time compilation）；
- 只嵌入x.cubin（无法进行just-in-time compilation）；
- 两者都嵌入（运行过程中driver找到合适二进制指令镜像则加载之，否则进行just-in-time compilation再加载之）。
我们在编译中添加了两个标志:
-- gpu-architecture=sm_70:这个标志告诉nvcc为 Volta GPU 编译生成二进制ptx。如果没有添加这个标志，编译器会尝试编译从sm_20（费米架构）兼容的代码，直到新的架构，也就是sm_70，也就是Volta gpu架构。老一代卡不支持动态并行，编译将失败，动态并行从sm_35开始支持的-rdc=true:这是在 GPU 上实现动态并行的一个关键参数。
完整代码如下：
#define MAX_DEPTH 16
#define INSERTION_SORT 32
__device__ void selection_sort(unsigned int *data, int left, int right) {
for (int i = left; i <= right; i++) {
unsigned int min_val = data[i];
int min_idx = i;
for (j = i + 1; j <= right; j++) {
unsigned int val_j = data[j];
if (val_j < min_val) {
min_idx = j;
min_val = val_j;
}
if (i != min_idx) {
data[min_idx] = data[i];
data[i] = min_val;
}
}
}
}
__global void cdp_simple_quicksort(unsigned int *data, int left, int right,
int depth) {
// 当递归的深度大于设定的MAX_DEPTH或者待排序的数组长度小于设定的阈值，直接调用简单选择排序
if (depth >= MAX_DEPTH || right - left <= INSERTION_SORT) {
selection_sort(data, left, right);
return;
}
unsigned int *left_ptr = data + left;
unsigned int *right_ptr = data + right;
unsigned int pivot = data[(left + right) / 2];
// partition
while (left_ptr <= right_ptr) {
unsigned int left_val = *left_ptr;
unsigned int right_val = *right_ptr;
while (left_val < pivot) { // 找到第一个比pivot大的
left_ptr++;
left_val = *left_ptr;
}
while (right_val > pivot) { // 找到第一个比pivot小的
right_ptr--;
right_val = *right_ptr;
}
// do swap
if (left_ptr <= right_ptr) {
*left_ptr++ = right_val;
*right_ptr-- = left_val;
}
}
// recursive
int n_right = right_ptr - data;
int n_left = left_ptr - data;
// Launch a new block to sort the the left part.
if (left < (right_ptr - data)) {
cudaStream_t l_stream;
// 设置非阻塞流
cudaStreamCreateWithFlags(&l_stream, cudaStreamNonBlocking);
cdp_simple_quicksort<<<1, 1, 0, l_stream>>>(data, left, n_right,
depth + 1);
cudaStreamDestroy(l_stream);
}
// Launch a new block to sort the the right part.
if ((left_ptr - data) < right) {
cudaStream_t r_stream;
// 设置非阻塞流
cudaStreamCreateWithFlags(&r_stream, cudaStreamNonBlocking);
cdp_simple_quicksort<<<1, 1, 0, r_stream>>>(data, n_left, right,
depth + 1);
cudaStreamDestroy(r_stream);
}
}
动态并行准则和约束
尽管动态并行为我们提供了在 GPU 上移植快速排序等算法的机会，但仍有一些基本规则和准则需要遵循。
编程模型规则:基本上所有的 CUDA 编程模型规则都适用:
- 每个线程的内核启动是异步的。
- 仅允许每个块同步。
- 创建的流在一个块内共享。
- 事件可用于创建流间依赖关系。
内存一致性规则:
- 子内核在启动时看到父内核的状态。
- 父内核可以看到子内核所做的更改，但是只能在同步之后。
- 本地和共享内存通常是私有的，不能被父内核传递或访问。
指引:
- 理解每次内核启动都会增加延迟也很重要。随着时间的推移，从另一个内核内部启动一个内核的延迟随着新架构逐渐减少。
- 虽然发射吞吐量比来自主机的高一个数量级，但是可以对最大深度进行限制。最新一代卡允许的最大深度是 24。
- 从内核内部执行
cudaDeviceSynchronize()是一个非常昂贵的操作，应该尽可能避免。 - 在全局内存上预先分配了额外的内存，这样我们就可以在内核启动之前存储它们。
- 如果内核出现故障，错误只能从主机上看到。因此，建议您使用
-lineinfo标志和cuda-memcheck来定位错误的位置。