CUDA系列学习(六) 从并行排序方法理解并行化思维——冒泡、归并、双调排序的GPU实现

本文主要是介绍CUDA系列学习(六) 从并行排序方法理解并行化思维——冒泡、归并、双调排序的GPU实现,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

转载自女神 —— Rachel-Zhang 的专栏

http://blog.csdn.net/abcjennifer/article/details/47110991


在第五讲中我们学习了GPU三个重要的基础并行算法: Reduce, Scan 和 Histogram,分析了 其作用与串并行实现方法。 在第六讲中,本文以冒泡排序 Bubble Sort、归并排序 Merge Sort 和排序网络中的双调排序 Bitonic Sort 为例, 讲解如何从数据结构课上学的串行并行排序方法转换到并行排序,并附GPU实现代码。

在并行方法中,我们将考虑到并行方法需要注意的特点进行设计,希望大家在读完本文后对GPU上并行算法的设计有一些粗浅的认识。需注意的特点如下: 
1. 充分发挥硬件能力(尽量不要有空闲且一直处于等待状态的SM) 
2. 限制branch divergence(见CUDA系列学习(二)) 
3. 尽量保证内存集中访问(即防止不命中)

( 而我们在数据结构课上学习的sort算法往往不注意这几点。)




CUDA系列学习目录:

CUDA系列学习(一)An Introduction to GPU and CUDA

CUDA系列学习(二)CUDA memory & variables - different memory and variable types

CUDA系列学习(三)GPU设计与结构QA & coding练习

CUDA系列学习(四)Parallel Task类型 与 Memory Allocation

CUDA系列学习(五)GPU基础算法: Reduce, Scan, Histogram





I. Bubble Sort

冒泡排序,相信大家再熟悉不过了。经典冒泡排序通过n轮有序冒泡(n为待排序的数组长度)得到有序序列, 其空间复杂度O(1), 时间复杂度O(n^2)。

那么如何将冒泡排序算法改成并行算法呢? 这里就需要解除一些依赖关系, 比如是否能解除n轮冒泡间的串行依赖 & 是否能解除每一轮冒泡内部的串行依赖, 使得同样的n^2次冒泡操作可以通过并行, 降低step complexity。 
1996年, J Kornerup针对这些问题提出了odd-even sort算法,并在论文中证明了其排序正确性。 


I.1 从Bubble Sort到Odd-even Sort

先来看一下odd-even sort的排序方法: 


这里写图片描述 
图1.1

上图为odd-even sort的基本方法。 
奇数步中, array中奇数项array[i]与右边的item(array[i + 1])比较; 
偶数步中, array中奇数项array[i]与左边的item(array[i - 1]) 比较;

这样,同一个step中的各个相邻比较就可以并行化了

PS: 对于array中有偶数个元素的数组也是一样: 


这里写图片描述 
图1.2



I.2 Odd-even Sort复杂度

在odd-even sort的算法下, 原本O(n^2)的总比较次数不变,但是由于并行,时间复杂度降到O(n), 即

step complexity = O(n) 
work complexity = O(n^2)

code详见 < Bubble sort and its variants >






II. Merge Sort

看过odd-even sort后,我们来看如何将归并排序并行化。数据结构课上我们都学过经典归并排序: 基于divide & conquer 思想, 每次将一个array分拆成两部分, 分别排序后合并两个有序序列。 可以通过 T(n)=2T(n/2)+n 得到, 其complexity = O(nlogn)。 和I.1节类似, 我们看看merge sort中的哪些步是可以并行的。

这里可以将基于merge sort的大规模数据排序分为三部分。 经过divide步之后, 数据分布如图所示: 


这里写图片描述 
图2.1

最下端的为 大量-短序列 的合并; 
中间一块为中等数量-中等长度序列 的合并; 
最上端的为少量-长序列 的合并;

我们分这三部分进行并行化。 之后大家会明白为啥要这么分~



II.1 Step 1: Huge number of small tasks

这一部分,每一部分序列合并的代价很小, 而有众多这样的任务。 所以我们采取给每个merge一个thread去执行, thread 内部就用串行merge的方法。 


II.2 Step 2: Mid number of mid task

在这一阶段, 有中等数量的task, 每个task的工作量也有一定增长。 所以我们采取给每个merge一个SM去执行, 每个block上运行的多个thread并行merge的方法。 和step1中的主要区别就是block内部merge从串行改成了并行。 那么怎样去做呢?

如下图所示,假如现在有两个长为4个元素的array, 想对其进行排序, 将merge排序结果index 0 - 7 写入数据下方方块。


这里写图片描述 
图2.2

做法: 
对于array中的每个数字, 看两个位置: 
1. 自己所在序列的位置: 看它前面有几个元素 
2. 另一个序列的位置: 另一个序列中有多少个元素比它小(采用二分搜索)

这样做来, 第一步O(1)可得到, 第二步二分查找O(logn)可得到; 整个merge过程用shared memory存结果。



II.3 Step 3: Small number of huge task

第三个环节中,也就是merge的顶端(最后一部分),每个merge任务的元素很多, 但是merge任务数很少。 这种情况下, 如果采用Step2的方法, 最坏的情况就是只有一个很大的task在跑, 此时只有1个SM在忙, 其他空闲SM却没法利用, 所以这里我们尝试将一个任务分给多个SM执行。

方法: 如图2.3所示, 将每个序列以256个元素为单位分段, 得到两个待merge序列In1和In2。然后,对这些端节点排序,如EABFCGDH。 如step2中的方法, 我们计算每个段节点在另一个短序列(长256)中的位置, 然后只对中间那些不确定的部分进行merge排序, 每个merge分配一个SM。

如E~A的部分, 
1. 计算出E在In1中的位置posE1, A在In2中的位置posA2 
2. merge In1中 posE1~A 和 In2中 E~posA2的元素


这里写图片描述 
图2.3



II.4 Merge sort in GPU

以上面step 1的merge sort为例, 其gpu代码中kernel函数如下:

其中temp为排好序的序列, 每次排序两个大小为sortedsize的block,为temp赋值2 * sortedsize个元素。

所以实际上,sortedsize就是一个sorted block的大小。

__global__ void mergeBlocks(int *a, int *temp, int sortedsize)
{int id = blockIdx.x;int index1 = id * 2 * sortedsize;int endIndex1 = index1 + sortedsize;int index2 = endIndex1;int endIndex2 = index2 + sortedsize;int targetIndex = id * 2 * sortedsize;int done = 0;while (!done){if ((index1 == endIndex1) && (index2 < endIndex2))temp[targetIndex++] = a[index2++];else if ((index2 == endIndex2) && (index1 < endIndex1))temp[targetIndex++] = a[index1++];else if (a[index1] < a[index2])temp[targetIndex++] = a[index1++];elsetemp[targetIndex++] = a[index2++];if ((index1 == endIndex1) && (index2 == endIndex2))done = 1;}
}
  • 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
  • 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

主函数中,定义block大小并调用kernel function:

        int blocks = BLOCKS/2;int sortedsize = THREADS;while (blocks > 0){mergeBlocks<<<blocks,1>>>(dev_a, dev_temp, sortedsize);cudaMemcpy(dev_a, dev_temp, N*sizeof(int), cudaMemcpyDeviceToDevice);blocks /= 2;sortedsize *= 2;}cudaMemcpy(a, dev_a, N*sizeof(int), cudaMemcpyDeviceToHost);
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11






III. Bitonic Sort



III.1 Bitonic Sequence 双调序列

不同于以上两种排序方法, 现在我们要接触的 双调排序 是排序网络方法中的一种。 想起当年在浙大面试某导师的实验室时就是让实现的双调排序, 并不断优化, 不过当时土得一坨, 就没听说过这个算法。。。 最后写出个多线程就结束了, 后来也没再整理。 现在我们来看看bitonic sort是个什么鬼。

双调排序是排序网络中最快的方法之一。所谓的排序网络是data-independent的排序, 即网络比较顺序与数据无关的排序方法, 所以特别适合硬件做并行化。

在了解双调排序算法之前,我们先来看看什么是双调序列。 双调序列是一个先单调递增后单调递减 或者 先单调递减后单调递增的序列。



III.2 双调排序算法

假如我们现在拿到了双调序列, 怎样对它按照从小到大进行排序呢?形象一点来看, 我们将一个双调序列切成两半, 每一段的单调性统一, 然后如下图图3.1所示, 将两段叠放起来, 进行两两比较, 这样一定能够在左右两段分别得到一个双调序列(想想为什么得到的是两个双调序列), 且左边的双调序列中元素全部小于右侧得到的双调序列的所有元素。 迭代这个过程, 每次都能将序列二分成两个子双调序列, 直到这个子双调序列的长度为2, 也就变成了一个单调子序列, 这个过程排序后原先的长双调序列就变为有序了 。 整个过程如下图图3.2所示。


这里写图片描述 
图3.1


这里写图片描述 
图3.2



III.3 任意序列生成双调序列

好,III.2中讲了怎样对双调序列进行排序, 那问题来了,怎样从任意序列生成双调序列呢? 这里可以看看本文最后的参考文献3, 写得很详细。 这个过程叫Bitonic merge, 实际上也是divide and conquer的思路。 和III.2中的思路正相反, 我们可以将两个相邻的,单调性相反的单调序列看作一个双调序列, 每次将这两个相邻的,单调性相反的单调序列merge生成一个新的双调序列, 然后排序(同III.2)。 这样只要每次两个相邻长度为n的序列的单调性相反, 就可以通过连接得到一个长度为2n的双调序列。 n开始为1, 每次翻倍,直到等于数组长度, 就只需要一遍单方向(单调性)排序了。


这里写图片描述 
图3.3

以16个元素的array为例, 
1. 相邻两个元素合并形成8个单调性相反的单调序列, 
2. 两两序列合并,形成4个双调序列,分别按相反单调性排序 
3. 4个长度为4的相反单调性单调序列,相邻两个合并,生成两个长度为8的双调序列,分别排序 
4. 2个长度为8的相反单调性单调序列,相邻两个合并,生成1个长度为16的双调序列,排序


总算讲完了,那么肿么实现呢? 我们看这个过程需要控制哪些地方? 如上图所示, 我们可以将len=16的array的双调排序分成4部分,每一部分结束都会形成若干长度为 i 的单调序列。 在每一部分中,用 j 表示比较的间隔,如下图所示每一时刻的i和j。


这里写图片描述 
图3.4 



III.4 双调排序的并行实现

本着“talk is cheap, show me the code”的优良作风, 拿粗来双调排序的GPU实现代码share如下:

/** Author: Rachel* <zhangruiqing01@baidu.com>** File: bitonic_sort.cu* Create Date: 2015-08-05 17:10:44**/#include<iostream>
#include<stdio.h>
#include<stdlib.h>
#include"gputimer.h"
#include<time.h>
#define NThreads 8
#define NBlocks 4#define Num NThreads*NBlocksusing namespace Gadgetron;__device__ void swap(int &a, int &b){int t = a;a = b;b = t;
}__global__ void bitonic_sort(int* arr){extern __shared__ int shared_arr[];const unsigned int tid = blockIdx.x * blockDim.x + threadIdx.x;//const unsigned int tid = threadIdx.x;shared_arr[tid] = arr[tid];__syncthreads();//for(int i=2; i<=blociDim.x; i<<=1){for(unsigned int i=2; i<=Num; i<<=1){for(unsigned int j=i>>1; j>0; j>>=1){unsigned int tid_comp = tid ^ j;if(tid_comp > tid){if((tid & i)==0){ //ascendingif(shared_arr[tid]>shared_arr[tid_comp]){swap(shared_arr[tid],shared_arr[tid_comp]);}}else{ //desendingif(shared_arr[tid]<shared_arr[tid_comp]){swap(shared_arr[tid],shared_arr[tid_comp]);}}}__syncthreads();}}arr[tid] = shared_arr[tid];
}int main(int argc, char* argv[])
{GPUTimer timer;int* arr= (int*) malloc(Num*sizeof(int));//init array valuetime_t t;srand((unsigned)time(&t));for(int i=0;i<Num;i++){arr[i] = rand() % 1000; }//init device variableint* ptr;cudaMalloc((void**)&ptr,Num*sizeof(int));cudaMemcpy(ptr,arr,Num*sizeof(int),cudaMemcpyHostToDevice);for(int i=0;i<Num;i++){printf("%d\t",arr[i]);}printf("\n");dim3 blocks(NBlocks,1);dim3 threads(NThreads,1);timer.start();bitonic_sort<<<blocks,threads,Num*sizeof(int)>>>(ptr);//bitonic_sort<<<1,Num,Num*sizeof(int)>>>(ptr);timer.stop();cudaMemcpy(arr,ptr,Num*sizeof(int),cudaMemcpyDeviceToHost);for(int i=0;i<Num;i++){printf("%d\t",arr[i]);}printf("\n");cudaFree(ptr);return 0;
}
  • 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
  • 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



code中, 
tid^j用于单方向判断, 防止同一元素比较两次; 
tid & i == 0 用于判断这个部分应该是单增还是单减, 因为方向在每个长为i的单调序列中是一致的, 所以选用i判断单调方向。


参考文献: 
1. Bubble sort and its variants 
2. nvidia的mergesort实现 
3. 我用过的浅显易懂的Bitonic sort文档


这篇关于CUDA系列学习(六) 从并行排序方法理解并行化思维——冒泡、归并、双调排序的GPU实现的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/791342

相关文章

hdu1043(八数码问题,广搜 + hash(实现状态压缩) )

利用康拓展开将一个排列映射成一个自然数,然后就变成了普通的广搜题。 #include<iostream>#include<algorithm>#include<string>#include<stack>#include<queue>#include<map>#include<stdio.h>#include<stdlib.h>#include<ctype.h>#inclu

认识、理解、分类——acm之搜索

普通搜索方法有两种:1、广度优先搜索;2、深度优先搜索; 更多搜索方法: 3、双向广度优先搜索; 4、启发式搜索(包括A*算法等); 搜索通常会用到的知识点:状态压缩(位压缩,利用hash思想压缩)。

【数据结构】——原来排序算法搞懂这些就行,轻松拿捏

前言:快速排序的实现最重要的是找基准值,下面让我们来了解如何实现找基准值 基准值的注释:在快排的过程中,每一次我们要取一个元素作为枢纽值,以这个数字来将序列划分为两部分。 在此我们采用三数取中法,也就是取左端、中间、右端三个数,然后进行排序,将中间数作为枢纽值。 快速排序实现主框架: //快速排序 void QuickSort(int* arr, int left, int rig

【C++】_list常用方法解析及模拟实现

相信自己的力量,只要对自己始终保持信心,尽自己最大努力去完成任何事,就算事情最终结果是失败了,努力了也不留遗憾。💓💓💓 目录   ✨说在前面 🍋知识点一:什么是list? •🌰1.list的定义 •🌰2.list的基本特性 •🌰3.常用接口介绍 🍋知识点二:list常用接口 •🌰1.默认成员函数 🔥构造函数(⭐) 🔥析构函数 •🌰2.list对象

【Prometheus】PromQL向量匹配实现不同标签的向量数据进行运算

✨✨ 欢迎大家来到景天科技苑✨✨ 🎈🎈 养成好习惯,先赞后看哦~🎈🎈 🏆 作者简介:景天科技苑 🏆《头衔》:大厂架构师,华为云开发者社区专家博主,阿里云开发者社区专家博主,CSDN全栈领域优质创作者,掘金优秀博主,51CTO博客专家等。 🏆《博客》:Python全栈,前后端开发,小程序开发,人工智能,js逆向,App逆向,网络系统安全,数据分析,Django,fastapi

让树莓派智能语音助手实现定时提醒功能

最初的时候是想直接在rasa 的chatbot上实现,因为rasa本身是带有remindschedule模块的。不过经过一番折腾后,忽然发现,chatbot上实现的定时,语音助手不一定会有响应。因为,我目前语音助手的代码设置了长时间无应答会结束对话,这样一来,chatbot定时提醒的触发就不会被语音助手获悉。那怎么让语音助手也具有定时提醒功能呢? 我最后选择的方法是用threading.Time

Android实现任意版本设置默认的锁屏壁纸和桌面壁纸(两张壁纸可不一致)

客户有些需求需要设置默认壁纸和锁屏壁纸  在默认情况下 这两个壁纸是相同的  如果需要默认的锁屏壁纸和桌面壁纸不一样 需要额外修改 Android13实现 替换默认桌面壁纸: 将图片文件替换frameworks/base/core/res/res/drawable-nodpi/default_wallpaper.*  (注意不能是bmp格式) 替换默认锁屏壁纸: 将图片资源放入vendo

usaco 1.3 Mixing Milk (结构体排序 qsort) and hdu 2020(sort)

到了这题学会了结构体排序 于是回去修改了 1.2 milking cows 的算法~ 结构体排序核心: 1.结构体定义 struct Milk{int price;int milks;}milk[5000]; 2.自定义的比较函数,若返回值为正,qsort 函数判定a>b ;为负,a<b;为0,a==b; int milkcmp(const void *va,c

C#实战|大乐透选号器[6]:实现实时显示已选择的红蓝球数量

哈喽,你好啊,我是雷工。 关于大乐透选号器在前面已经记录了5篇笔记,这是第6篇; 接下来实现实时显示当前选中红球数量,蓝球数量; 以下为练习笔记。 01 效果演示 当选择和取消选择红球或蓝球时,在对应的位置显示实时已选择的红球、蓝球的数量; 02 标签名称 分别设置Label标签名称为:lblRedCount、lblBlueCount

浅谈主机加固,六种有效的主机加固方法

在数字化时代,数据的价值不言而喻,但随之而来的安全威胁也日益严峻。从勒索病毒到内部泄露,企业的数据安全面临着前所未有的挑战。为了应对这些挑战,一种全新的主机加固解决方案应运而生。 MCK主机加固解决方案,采用先进的安全容器中间件技术,构建起一套内核级的纵深立体防护体系。这一体系突破了传统安全防护的局限,即使在管理员权限被恶意利用的情况下,也能确保服务器的安全稳定运行。 普适主机加固措施: