简介

排序是任何应用的基本构造块的关键算法之一。有许多排序算法已经被广泛研究，常见的排序算法时间和空间复杂度如下：

一些排序算法属于分治算法的范畴。这些算法适用于并行性，并适合 GPU 等架构，在这些架构中，要排序的数据可以被划分进行排序。快速排序就是这样一种算法。如前所述，Quicksort 属于“分而治之”的范畴。这是一个三步走的方法，如下所示:

1. 从需要排序的数组中选择一个元素。该元素充当枢轴元素。
2. 第二步是划分所有元素的位置。小于枢轴的所有元素都向左移动，大于或等于枢轴的所有元素都向右移动。这一步也称为分区。
3. 递归地执行步骤 1 和 2，直到所有的子数组都被排序。

3 种快排枢轴元素选择方法：

• 随机（rand函数）：几乎不会出现最坏情况的复杂性
• 固定（队首、队尾）
• 三数取中（队首、队中和队尾的中间数）

与其他排序算法(如需要额外存储的合并排序)相比，快速排序的内存负载和要求更低。快速排序的更实际的实现使用随机化版本。随机化版本的预期时间复杂度为 $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内核，可以让递归算法更加清晰易懂，也更加容易理解。下图简化了解释:

如果我们将这个概念转化为快速排序的执行方式，它看起来会是这样的:

扫描二维码关注公众号，回复： 14655686 查看本文章

深度 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 内核。实施快速排序的主要步骤如下:

1. CPU 启动第一个内核:内核一个块一个线程启动。左边的元素是数组的开始，而右边是数组的最后一个元素(基本上是整个数组):

int main(int argc, char **argv)
{
    
     ...
    cdp_simple_quicksort<<< 1, 1 >>>(data, left, right, 0);
}

2. 极限检查:在从内核内部启动内核之前，检查两个标准。首先，检查我们是否已经达到硬件允许的最大深度限制。其次，我们需要检查子数组中要排序的元素数量是否小于扭曲大小(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;
 }

3. 划分:如果满足前面的条件，那么将数组划分为两个子数组，并启动两个新的内核，一个用于左数组，另一个用于右数组。并从从内核内部启动一个内核:

__global__ void cdp_simple_quicksort( unsigned int *data, int left, int right, int depth ) {
    
    
...
while(left_ptr <= right_ptr)
 {
    
    
     // Launch a new block to sort the left part.
     if(left < (right_ptr-data))
     {
    
     // Create a new stream for the eft sub array
        cdp_simple_quicksort<<< 1, 1, 0, s >>>(data, left, n_right, depth+1);
     }
    // Launch a new block to sort the right part.
    if((left_ptr-data) < right)
     {
    
    //Create stream for the right sub array
         cdp_simple_quicksort<<< 1, 1, 0, s1 >>>(data, n_left, right, depth+1);
     }
 }

4. 执行代码:使用以下命令用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 上实现动态并行的一个关键参数。

完整代码如下：

#include <cstdio>
#include <iostream>

#include "cuda_runtime_api.h"


#define MAX_DEPTH 16
#define INSERTION_SORT 32

#define GPU_CHECK(ans)                                                         \
    {
      
       GPUAssert((ans), __FILE__, __LINE__); }
inline void GPUAssert(cudaError_t code, const char *file, int line,
                      bool abort = true) {
    
    
    if (code != cudaSuccess) {
    
    
        fprintf(stderr, "GPUassert: %s %s %d\n", cudaGetErrorString(code), file,
                line);
        if (abort)
            exit(code);
    }
};

__device__ void selection_sort(unsigned int *data, int left, int right) {
    
    
    for (int i = left; i <= right; ++i) {
    
    
        unsigned min_val = data[i];
        int min_idx = i;

        // Find the smallest value in the range [left, right].
        for (int j = i + 1; j <= right; ++j) {
    
    
            unsigned val_j = data[j];

            if (val_j < min_val) {
    
    
                min_idx = j;
                min_val = val_j;
            }
        }

        // Swap the values.
        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);
    }
}

// Call the quicksort kernel from the host.
void run_qsort(unsigned int *data, unsigned int nitems) {
    
    
    // Prepare CDP for the max depth 'MAX_DEPTH'.
    GPU_CHECK(cudaDeviceSetLimit(cudaLimitDevRuntimeSyncDepth, MAX_DEPTH));

    int left = 0;
    int right = nitems - 1;
    cdp_simple_quicksort<<<1, 1>>>(data, left, right, 0);
    GPU_CHECK(cudaDeviceSynchronize());
}

// Initialize data on the host.
void initialize_data(unsigned int *dst, unsigned int nitems) {
    
    
    // Fixed seed for illustration
    srand(2047);

    // Fill dst with random values
    for (unsigned i = 0; i < nitems; i++)
        dst[i] = rand() % nitems;
}

// Verify the results.
void print_results(int n, unsigned int *results_d) {
    
    
    unsigned int *results_h = new unsigned[n];
    GPU_CHECK(cudaMemcpy(results_h, results_d, n * sizeof(unsigned),cudaMemcpyDeviceToHost));
    std::cout << "Sort data : ";
    for (int i = 1; i < n; ++i){
    
    
        std::cout << results_h[i] << " ";
        // if (results_h[i - 1] > results_h[i]) {
    
    
        //     std::cout << "Invalid item[" << i - 1 << "]: " << results_h[i - 1]
        //               << " greater than " << results_h[i] << std::endl;
        //     exit(EXIT_FAILURE);
        // }
    }
    std::cout << std::endl;
    delete[] results_h;
}

void check_results(int n, unsigned int *results_d)
{
    
    
    unsigned int *results_h = new unsigned[n];
    GPU_CHECK(cudaMemcpy(results_h, results_d, n*sizeof(unsigned), cudaMemcpyDeviceToHost));

    for (int i = 1 ; i < n ; ++i)
        if (results_h[i-1] > results_h[i])
        {
    
    
            std::cout << "Invalid item[" << i-1 << "]: " << results_h[i-1] << " greater than " << results_h[i] << std::endl;
            exit(EXIT_FAILURE);
        }

    std::cout << "OK" << std::endl;
    delete[] results_h;
}

// Main entry point.

int main(int argc, char **argv) {
    
    
    int num_items = 100;
    bool verbose = true;

    // Find/set device and get device properties
    int device = 1;
    cudaDeviceProp deviceProp;
    GPU_CHECK(cudaGetDeviceProperties(&deviceProp, device));

    if (!(deviceProp.major > 3 ||
          (deviceProp.major == 3 && deviceProp.minor >= 5))) {
    
    
        printf("GPU %d - %s  does not support CUDA Dynamic Parallelism\n Exiting.",
            device, deviceProp.name);
        return 0;
    }

    // Create input data
    unsigned int *h_data = 0;
    unsigned int *d_data = 0;

    // Allocate CPU memory and initialize data.
    std::cout << "Initializing data:" << std::endl;
    h_data = (unsigned int *)malloc(num_items * sizeof(unsigned int));
    initialize_data(h_data, num_items);

    if (verbose) {
    
    
       std::cout << "Raw  data : ";
        for (int i = 0; i < num_items; i++)
            std::cout << h_data[i] << " ";
    }
    std::cout << std::endl;

    // Allocate GPU memory.
    GPU_CHECK(cudaMalloc((void **)&d_data, num_items * sizeof(unsigned int)));
    GPU_CHECK(cudaMemcpy(d_data, h_data, num_items * sizeof(unsigned int),cudaMemcpyHostToDevice));

    // Execute
    std::cout << "Running quicksort on " << num_items << " elements" << std::endl;
    run_qsort(d_data, num_items);

    // print result
    print_results(num_items, d_data);
    // check result
    std::cout << "Validating results: ";
    check_results(num_items, d_data);
    free(h_data);
    GPU_CHECK(cudaFree(d_data));

    return 0;
}

/* 
编译：nvcc -o quicksort_cuda --gpu-architecture=sm_70 -rdc=truequicksort_cuda.cu
*/

耗时测试对比

测试自己实现c++版、cuda版、Thrust库的sort函数三种排序算法的时间耗时，计算耗时时间未考虑主机和设备之间数据复制耗时，纯统计排序耗时时间。

• 测试500万数据

• 测试50万数据

• 测试5万数据

• 测试5千数据

由上面4组数据可以看出，本文实现的quicksort 比thrust库sort快，且在大数据排序，cuda实现排序比c++快很多！

动态并行准则和约束

尽管动态并行为我们提供了在 GPU 上移植快速排序等算法的机会，但仍有一些基本规则和准则需要遵循。

编程模型规则:基本上所有的 CUDA 编程模型规则都适用:

• 每个线程的内核启动是异步的。
• 仅允许每个块同步。
• 创建的流在一个块内共享。
• 事件可用于创建流间依赖关系。

内存一致性规则:

• 子内核在启动时看到父内核的状态。
• 父内核可以看到子内核所做的更改，但是只能在同步之后。
• 本地和共享内存通常是私有的，不能被父内核传递或访问。

指引:

• 理解每次内核启动都会增加延迟也很重要。随着时间的推移，从另一个内核内部启动一个内核的延迟随着新架构逐渐减少。
• 虽然发射吞吐量比来自主机的高一个数量级，但是可以对最大深度进行限制。最新一代卡允许的最大深度是 24。
• 从内核内部执行cudaDeviceSynchronize()是一个非常昂贵的操作，应该尽可能避免。
• 在全局内存上预先分配了额外的内存，这样我们就可以在内核启动之前存储它们。
• 如果内核出现故障，错误只能从主机上看到。因此，建议您使用-lineinfo标志和cuda-memcheck来定位错误的位置。