Sort排序的CUDA计算优化加速-冒泡\并归\双调
CUDA sort
- 基于reduce_sum的思想，利用cuda对排序进行优化是一个常见的进阶编程问题
冒泡排序
- Bubble Sort冒泡排序，空间复杂度O(1)，时间复杂度O(n^2)
- 1996年，J Kornerup 提出了odd-even sort算法
- odd-even sort的基本方法：
- 奇数步中， array中奇数项array[i]与右边的item(array[i + 1])比较；
- 偶数步中， array中奇数项array[i]与左边的item(array[i - 1]) 比较；
- 这样，同一个step中的各个相邻比较就可以并行化了。
- 在odd-even sort的算法下， 原本O(n^2)的总比较次数不变，但是由于并行，时间复杂度降到O(n)
#include<stdio.h>
#include<stdlib.h>
#include<cuda_runtime.h>
#include<iostream>
#include"cuda.h"
#define N 1024
// 冒泡排序 CPU实现
void cpuBubbleSort(int* nums) {
for (int i = 0; i < N-1; i++) {
// 每次冒泡将最大值推到末尾
for (int j = 0; j < N-i-1; j++) {
if (nums[j] > nums[j+1]) {
std::swap(nums[j], nums[j+1]);
}
}
}
}
// 冒泡排序 odd-even sort算法
// 原本O(n^2)的总比较次数不变 但是由于并行 时间复杂度降到O(n)
__global__ void gpuBubbleSort(int* nums) {
int i = 0;
int tmp = 0;
int tid = blockIdx.x * blockDim.x + threadIdx.x; // 此处仅一个block
// 交换相邻元素-奇数步骤
// array中奇数项array[i]与右边的item(array[i + 1])比较
if((2*tid + 1) < N)
{
if(nums[2*tid]>nums[2*tid + 1])
{
tmp = nums[2*tid];
nums[2*tid] = nums[2*tid + 1];
nums[2*tid + 1] = tmp;
}
}
__syncthreads();
// 交换相邻元素-偶数步骤
// array中奇数项array[i]与左边的item(array[i - 1]) 比较
if((2*tid +2) < N)
{
if(nums[2*tid+1]>nums[2*tid + 2])
{
tmp = nums[2*tid + 1];
nums[2*tid + 1] = nums[2*tid + 2];
nums[2*tid + 2] = tmp;
}
}
__syncthreads();
}
并归排序
- CUDA sort使用GPU加速并归排序的计算
- 归并操作是将两个有序序列合并成一个有序序列的过程， nlog(n) 的复杂度，这通常涉及以下步骤：
- 申请一个足够大的临时数组，用于存放合并后的序列。
- 设置两个指针，分别指向两个待合并序列的起始位置。
- 比较两个指针所指向的元素，将较小的元素复制到临时数组中，并移动相应的指针。
- 重复上述步骤，直到某个序列的元素全部被复制到临时数组中。
- 将另一个序列中剩余的元素（如果有）直接复制到临时数组的末尾。
- 并归排序是如何实现：
// 归并排序 CPU实现 nlog(n)复杂度
// 合并两个有序子数组
void merge(int *nums, int left, int mid, int right) {
int temp[N] = {0};
int i = left, j = mid+1, k = 0;
while (i <= mid && j <= right) {
temp[k++] = (nums[i] <= nums[j]) ? nums[i++] : nums[j++];
}
// 拷贝较长的剩余元素
while (i <= mid) temp[k++] = nums[i++];
while (j <= right) temp[k++] = nums[j++];
for (int p = 0; p < k; p++) {
nums[left + p] = temp[p];
}
}
void cpuMergeSort(int *nums, int left, int right) {
if (left >= right) return;
int mid = left + (right - left) / 2;
// 对左右子数组进行排序
cpuMergeSort(nums, left, mid);
cpuMergeSort(nums, mid+1, right);
// 合并已排序的子数组
merge(nums, left, mid, right);
}
// GPU实现合并两个有序块
__global__ void gpuMergeSort(int* nums, int n, int blockSize) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
int start = idx * 2 * blockSize;
int mid = start + blockSize;
int end = min(start + 2 * blockSize, N);
if (mid >= end) return;
int i = start, j = mid, k = 0;
int* temp = new int[end - start];
while (i < mid && j < end) {
temp[k++] = (nums[i] <= nums[j]) ? nums[i++] : nums[j++];
}
while (i < mid) temp[k++] = nums[i++];
while (j < end) temp[k++] = nums[j++];
for (int p = 0; p < k; p++) {
nums[start + p] = temp[p];
}
delete[] temp;
}
- 并归排序可以参考这篇文章
双调排序
- 为什么要使用双调排序而不是并归排序？
- 双调排序 Bitonic Sort，整体复杂度为 n*log^2(n)
- 归并排序的缺点： 当归并排序到后期时，待归并的段主要有少数的“大段”构成，而这几个大段的归并比较难运用到多核多线程的优势。
- 双调序列是一个先单调递增后单调递减（或者先单调递减后单调递增）的序列。
- 即存在一个索引i（0≤i≤n-1），使得序列A[0…n-1]满足：
- A[0] < A[1] < A[2] … < A[i-1] < A[i] > A[i+1] > A[i+2] > … > A[n-1]
- 或者反过来，先递减后递增。
- 双调排序的基本原理
- 构建双调序列：
- 对于任意给定的序列，首先通过一系列的比较和交换操作，将其转换为一个双调序列。
- 这个过程通常是通过将序列分成较小的子序列，对每个子序列进行排序（升序或降序），然后按照特定的方式合并这些子序列来实现的。
- 双调归并：
- 一旦构建了双调序列，就可以使用双调归并过程来对其进行排序。
- 双调归并是基于Batcher定理：
- 该定理指出，将任意一个长为2n的双调序列A分为等长的两半X和Y，将X中的元素与Y中的元素一一按原序比较（即a[i]与a[i+n]（i<n）比较），将较大者放入MAX序列，较小者放入MIN序列。则得到的MAX和MIN序列仍然是双调序列，并且MAX序列中的任意一个元素不小于MIN序列中的任意一个元素。
- 通过递归地应用这个过程，最终可以得到一个有序序列。
- 构建双调序列：
- 一个双调序列如何排序？
- 假设我们有一个双调序列，则我们根据Batcher定理，将该序列划分成2个双调序列，然后继续对每个双调序列递归划分，得到更短的双调序列，直到得到的子序列长度为1为止。这时的输出序列按单调递增顺序排列。
- 见下图：升序排序，具体方法是，把一个序列(1…n)对半分，假设n=2^k，然后1和n/2+1比较，小的放上，接下来2和n/2+2比较，小的放上，以此类推；然后看成两个(n/2)长度的序列，因为他们都是双调序列，所以可以重复上面的过程；总共重复k轮，即最后一轮已经是长度是2的序列比较了，就可得到最终的排序结果。
- 任意序列如何变成一个双调序列？
- Bitonic Sort（双调排序）
- 对于两个元素x，y，如果x<=y，则x，y都位于双调序列的递增部分，而递减部分没有元素，如果x>=y，则x，y都位于双调序列的递减部分，而递增部分没有元素，于是x和y构成一个双调序列。因此，任何无序的序列都是由若干个只有2个元素的双调序列连接而成。
- 于是，对于一个无序序列，我们按照递增和递减顺序合并相邻的双调序列，按照双调序列的定义，通过连接递增和递减序列得到的序列是双调的。最终，我们可以将若干个只有2个元素的双调序列合并成1个有n个元素的双调序列。
以16个元素的array为例：
相邻两个元素合并形成8个单调性相反的单调序列
两两序列合并，形成4个双调序列，分别按相反单调性排序
4个长度为4的相反单调性单调序列，相邻两个合并，生成两个长度为8的双调序列，分别排序
2个长度为8的相反单调性单调序列，相邻两个合并，生成1个长度为16的双调序列，排序
- 以下的例子中，先比较相邻的2个元素，按照 递增->递减->递增->递减 排序
- 再比较4个元素，按照 递增->递减 排序，4个元素内部进行对比排序
- 非2的幂次长度序列排序？
- 这样的双调排序算法只能应付长度为2的幂的数组。那如何转化为能针对任意长度的数组呢？
- 一个直观的方法就是使用padding。即使用一个定义的最大或者最小者来填充数组，让数组的大小填充到2的幂长度，再进行排序。最后过滤掉那些最大（最小）值即可。
- 这种方式会使用到额外的空间，而且有时候padding的空间比较大（如数组长度为1025个元素，则需要填充到2048个，浪费了大量空间）。但是这种方法比较容易转化为针对GPU的并行算法。所以一般来说，并行计算中常使用双调排序来对一些较小的数组进行排序
// CPU实现双调排序 dir调整升序降序 表示sorting direction
// 需要限制输入数组为2次幂 可以padding数组以实现填充长度
void bitonicMerge(int* nums, int low, int cnt, bool dir)
{
if (cnt > 1) {
int mid = cnt / 2;
for (int i = low; i < low + mid; i++) {
if (dir == (nums[i] > nums[i + mid])) {
std::swap(nums[i], nums[i + mid]);
}
}
bitonicMerge(nums, low, mid, dir);
bitonicMerge(nums, low + mid, mid, dir);
}
}
void cpuBitonicSort(int* nums, int low, int cnt, bool dir) {
if (cnt > 1) {
int mid = cnt / 2; // mid表示每次处理的长度
cpuBitonicSort(nums, low, mid, true); // 前半升序
cpuBitonicSort(nums, low + mid, mid, false); // 后半降序
bitonicMerge(nums, low, cnt, dir);
}
}
// GPU实现双调排序核函数
__global__ void gpuBitonicSort(int* nums, int stride, int inner_stride) {
int idx = threadIdx.x + blockIdx.x * blockDim.x;
int flipper = inner_stride/2;
int idx_other = idx ^ flipper; // 按位异或 计算对称索引
if (idx < idx_other) { // 操纵左侧的半部分
if ((idx & stride) == 0) { // 升序比较
if (nums[idx] > nums[idx_other]) {
int temp = nums[idx];
nums[idx] = nums[idx_other];
nums[idx_other] = temp;
}
}
else { // 降序比较
if (nums[idx] < nums[idx_other]) {
int temp = nums[idx];
nums[idx] = nums[idx_other];
nums[idx_other] = temp;
}
}
}
}
Sort算子实现完整测试代码：
#include<stdio.h>
#include<stdlib.h>
#include<cuda_runtime.h>
#include<iostream>
#include"cuda.h"
#define N 1024
// 冒泡排序 CPU实现
void cpuBubbleSort(int* nums) {
for (int i = 0; i < N-1; i++) {
// 每次冒泡将最大值推到末尾
for (int j = 0; j < N-i-1; j++) {
if (nums[j] > nums[j+1]) {
std::swap(nums[j], nums[j+1]);
}
}
}
}
// 冒泡排序 odd-even sort算法
// 原本O(n^2)的总比较次数不变 但是由于并行 时间复杂度降到O(n)
__global__ void gpuBubbleSort(int* nums) {
int i = 0;
int tmp = 0;
int tid = blockIdx.x * blockDim.x + threadIdx.x; // 此处仅一个block
// 交换相邻元素-奇数步骤
// array中奇数项array[i]与右边的item(array[i + 1])比较
if((2*tid + 1) < N)
{
if(nums[2*tid]>nums[2*tid + 1])
{
tmp = nums[2*tid];
nums[2*tid] = nums[2*tid + 1];
nums[2*tid + 1] = tmp;
}
}
__syncthreads();
// 交换相邻元素-偶数步骤
// array中奇数项array[i]与左边的item(array[i - 1]) 比较
if((2*tid +2) < N)
{
if(nums[2*tid+1]>nums[2*tid + 2])
{
tmp = nums[2*tid + 1];
nums[2*tid + 1] = nums[2*tid + 2];
nums[2*tid + 2] = tmp;
}
}
__syncthreads();
}
// 归并排序 CPU实现 nlog(n)复杂度
// 合并两个有序子数组
void merge(int *nums, int left, int mid, int right) {
int temp[N] = {0};
int i = left, j = mid+1, k = 0;
while (i <= mid && j <= right) {
temp[k++] = (nums[i] <= nums[j]) ? nums[i++] : nums[j++];
}
// 拷贝较长的剩余元素
while (i <= mid) temp[k++] = nums[i++];
while (j <= right) temp[k++] = nums[j++];
for (int p = 0; p < k; p++) {
nums[left + p] = temp[p];
}
}
void cpuMergeSort(int *nums, int left, int right) {
if (left >= right) return;
int mid = left + (right - left) / 2;
// 对左右子数组进行排序
cpuMergeSort(nums, left, mid);
cpuMergeSort(nums, mid+1, right);
// 合并已排序的子数组
merge(nums, left, mid, right);
}
// GPU实现合并两个有序块
__global__ void gpuMergeSort(int* nums, int n, int blockSize) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
int start = idx * 2 * blockSize;
int mid = start + blockSize;
int end = min(start + 2 * blockSize, N);
if (mid >= end) return;
int i = start, j = mid, k = 0;
int* temp = new int[end - start];
while (i < mid && j < end) {
temp[k++] = (nums[i] <= nums[j]) ? nums[i++] : nums[j++];
}
while (i < mid) temp[k++] = nums[i++];
while (j < end) temp[k++] = nums[j++];
for (int p = 0; p < k; p++) {
nums[start + p] = temp[p];
}
delete[] temp;
}
// CPU实现双调排序 dir调整升序降序 表示sorting direction
// 需要限制输入数组为2次幂 可以padding数组以实现填充长度
void bitonicMerge(int* nums, int low, int cnt, bool dir)
{
if (cnt > 1) {
int mid = cnt / 2;
for (int i = low; i < low + mid; i++) {
if (dir == (nums[i] > nums[i + mid])) {
std::swap(nums[i], nums[i + mid]);
}
}
bitonicMerge(nums, low, mid, dir);
bitonicMerge(nums, low + mid, mid, dir);
}
}
void cpuBitonicSort(int* nums, int low, int cnt, bool dir) {
if (cnt > 1) {
int mid = cnt / 2; // mid表示每次处理的长度
cpuBitonicSort(nums, low, mid, true); // 前半升序
cpuBitonicSort(nums, low + mid, mid, false); // 后半降序
bitonicMerge(nums, low, cnt, dir);
}
}
// GPU实现双调排序核函数
__global__ void gpuBitonicSort(int* nums, int stride, int inner_stride) {
int idx = threadIdx.x + blockIdx.x * blockDim.x;
int flipper = inner_stride/2;
int idx_other = idx ^ flipper; // 按位异或 计算对称索引
if (idx < idx_other) { // 操纵左侧的半部分
if ((idx & stride) == 0) { // 升序比较
if (nums[idx] > nums[idx_other]) {
int temp = nums[idx];
nums[idx] = nums[idx_other];
nums[idx_other] = temp;
}
}
else { // 降序比较
if (nums[idx] < nums[idx_other]) {
int temp = nums[idx];
nums[idx] = nums[idx_other];
nums[idx_other] = temp;
}
}
}
}
int main() {
int i=0;
// CPU内存
std::cout << " nums: " << std::endl;
int nums[N] = {0};
for(i=0; i<N; i++){
nums[i] = rand()%N;
std::cout << nums[i] << " "; // 输出顺序可能不正确（竞争条件）
}
// CPU实现
// cpuBubbleSort(nums);
// std::cout << " cpu bubble sort nums: " << std::endl;
// for(i=0; i<N; i++){
// std::cout << nums[i] << " "; // 输出顺序可能不正确（竞争条件）
// }
// cpuMergeSort(nums, 0, N-1);
// std::cout << " cpu merge sort nums: " << std::endl;
// for(i=0; i<N; i++){
// std::cout << nums[i] << " "; // 输出顺序可能不正确（竞争条件）
// }
// cpuBitonicSort(nums, 0, N, true);
// std::cout << " cpu bitonic sort nums: " << std::endl;
// for(i=0; i<N; i++){
// std::cout << nums[i] << " "; // 输出顺序可能不正确（竞争条件）
// }
// GPU内存
int *dev_nums = NULL;
cudaMalloc((void**)&dev_nums, N*sizeof(int));
cudaMemcpy(dev_nums, nums, N*sizeof(int), cudaMemcpyHostToDevice);
// GPU实现
// CUDA版冒泡排序由于并行竞争问题无法正确运行 对数组的写需要进行原子操作 效率极低
// gpuBubbleSort<<<1,N/2>>>(dev_nums);
// 并归 需要分块逐步合并
// 当归并排序到后期 待归并的段主要有少数的“大段”构成 较难运用到多核多线程的优势
// for (int blockSize = 1; blockSize < N; blockSize *= 2) {
// int threads = 32;
// int blocks = (N + (2*blockSize * threads) - 1) / (2*blockSize * threads);
// gpuMergeSort<<<blocks, threads>>>(dev_nums, N, blockSize);
// }
// 双调
for (int stride = 2; stride <= N; stride *= 2) { // 序列长度翻倍
for (int inner_stride = stride; inner_stride >= 2; inner_stride /= 2) { // 步长逐步减半
dim3 grid_dim((N / 256 == 0)? 1 : N / 256);
dim3 block_dim((N / 256 == 0)? N : 256);
gpuBitonicSort<<<grid_dim, block_dim>>>(dev_nums, stride, inner_stride);
}
}
cudaMemcpy(nums, dev_nums, N*sizeof(int), cudaMemcpyDeviceToHost);
cudaFree(dev_nums);
std::cout << " gpu sort nums: " << std::endl;
for(i=0; i<N; i++){
std::cout << nums[i] << " "; // 输出顺序可能不正确（竞争条件）
}
system("pause");
return 0;
}
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
