CUDA动态并行实现快速排序

简介

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

一些排序算法属于分治算法的范畴。这些算法适用于并行性，并适合 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);
    }
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

在本文中，我们将利用动态并行（从 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 );
1
2
3
4

这是非常重要的一步，因为，否则，内核的启动可能会被序列化。

CUDA快速排序

对于我们的快速排序实现，我们将利用动态并行性递归启动 GPU 内核。实施快速排序的主要步骤如下:

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

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

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;
 }
1
2
3
4
5
6
7
8

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);
     }
 }
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

4. 执行代码:使用以下命令用nvcc编译器编译您的应用:

nvcc -o quick_sort --gpu-architecture=sm_70 -rdc=true quick_sort.cu 
1

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 
#include 

#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
*/
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211

耗时测试对比

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

• 测试500万数据

• 测试50万数据

• 测试5万数据

• 测试5千数据

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

动态并行准则和约束

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

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

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

内存一致性规则:

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

指引:

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