一种基于划分与归并的并行快速排序算法∗

贾思禹

（天津工业大学计算机科学与软件学院 天津 300387）

1 引言

排序操作对于大数据分析显得尤为重要，这种重要性尤其体现在数据挖掘，商业化应用程序，大规模科学计算等领域。在所有的排序算法中，快速排序（QuickSort）算法［1］作为一种技术人员广为熟知的优秀算法，相对于现存的其它各类排序算法而言，在实际生产过程中，具有速度快，效率高等优点［2～3］，且近年来被国内外学者从各种不同的角度进行了改进与优化。

2 研究背景

现如今，大多数的并行体系结构计算机都曾被用于执行快速排序算法，早期研究显示［4～5］排序算法可以执行在一种名为multi-ring 的可重构多处理器网络中。而近年来又出现了大量增强经典快速排序算法的并行快速排序算法。2003 年，P Tsigas和Y Zhang［6］提出了一种实现在具有cache 一致性机制的共享地址空间微处理器之上的并行快速排序算法，这一算法中数据可以适应L1 cache。Codish M［7］等从标准排序函数库的排序算法网络入手，对排序算法底层实现进行优化，性能提升明显。D Koch 和J Torresen［8］在2011 年提出了一种实现在现场可编程逻辑门阵列（FPGA）之上的利用运行时重构的高性能排序体系结构，在这一架构上运行的排序实验结果显示，使用FPGA 片上存储器时可以达到2GB/s 的排序数据吞吐量，在使用外部存储器的时候也可以达到1GB/s 的数据吞吐量。图形处理器（GPU）内部的成千上万个处理核心都可作为协处理器用于进行基于单指令多数据流（SIMD）并行机制的快速排序算法，例如运行在英特尔Xeon Phi 处理器上的Bitonic-Merge 排序算法［9］，又例如E Manca［10］等又将快速排序算法实现在了CUDA 并行计算架构之上，CUDATM是一种由NVIDIA 推出的通用并行计算架构，该架构使GPU 能够解决复杂的计算问题，他们实现的两种算法：基于GPU 的快速排序以及基于CUDA 的快速排序，他们的实验结果表明基于CUDA 的快速排序速度是基于GPU 的快速排序的四倍以上。

近年来，国内学者也从各种不同的角度对基本快速排序进行了改进与优化。

庹清等［11］从混合快速排序与其他排序的角度进行优化，其算法思想是将基本快速排序与计数排序混合，从而构造高效的排序方法。其算法效率远高于一般快速排序，缺点是在降低算法时间复杂度的同时提高了空间复杂度，其算法的本质是用更多的空间来换取更快的速度。

汤亚玲等［12］从概率学的角度，先对待排序数据进行统计并选取出最佳的枢轴元素，以确保快速排序的每一次划分位置都正好处于待排序序列的正中间。其算法的本质是先选取合适的枢轴，再进行基本快速排序，而选取合适的枢轴本身就是一个很浪费时间的操作，因此其方法只能在某些特定情况下提高排序效率，而在另一些情况下反而效率低于基本快速排序。

2.1 OpenMP并行函数库

OpenMP并行函数库是一款非常知名的并行函数库，使用它可以在多核处理器架构之上开发并行应用程序。它为共享存储多核处理器平台提供了一套基于线程并行机制的应用程序编程接口（API）。这些API 中包括了一系列编译器指令，库函数，以及环境变量等，用以支持在多种并行架构上使用FORTRAN 或C/C++语言开发并行应用程序。OpenMP使用fork-join 模型作为其多线程执行模型，所谓fork-join 模型是指在一个并行区域里，多个线程并发地执行指令。使用OpenMP 的最大优点是所有的处理器核心都可以共享地访问到相同的存储器数据，即所有的处理器核心访问存储单元有着相同的延迟和带宽，换言之，OpenMP相对于其它非共享存储并行计算解决方案（如集群计算、网格计算）来说用于通信的额外开销极小，而且几乎没有通信延迟。

OpenMP程序运行开始后只有主线程会立即执行，当执行到使用OpenMP 指令定义的并行区域时，多个线程才能并行执行并行区域中的代码，而当程序执行完并行区域的代码后，进入串行区域后，主线程继续单独执行代码，直至遇到下一个并行区域，即在两个并行区域之间，除了主线程执行代码外，其他任何线程将不执行。fork-join 模型示意图如图1所示。

图1 fork-join模型示意图

2008 年5 月，OpenMP3.0 版本问世，其中Task结构在此版本中首次被提出，此结构用于定义一个显式的任务。当一个线程在执行过程中遇到Task结构时，相关的结构代码块就会生成为一个显式的任务。此外Task 结构允许嵌套，也就是说可以在外围Task 结构中定义内嵌Task 结构。2013 年7月，OpenMP4.0版本问世，此时的OpenMP已经可以支持单指令多数据流（Single Instruction Multiple Data，SIMD），SIMD 是指多核处理器可以同时的对多组数据执行相同的操作。SIMD 指令集基本原理如图2所示。

图2 SIMD指令集基本原理

图2 中R1 和R2 分别为两个SIMD 专用寄存器，每个寄存器整体长度为128 位，用于存储需要处理的操作数。一条指令同时启动4 个整型数相加操作，两个寄存器在同一指令周期内分别提供四个32 位操作数经过硬件加法器相加后，将结果保存到R3 寄存器中，R3 寄存器同样也是SIMD 专用128 位寄存器，此时仅有一条加法指令，但却对四组数据执行了加法操作，这便是单指令多数据流的基本原理。

2.2 现有研究成果及比较

DUHU MAN［13］等开发了一种可以兼容标准库快速排序函数qsort 的psort 算法，他们的算法将待排序数据集划分为多个数据块，并对每个数据块进行qsort 排序，但在归并时使用了额外空间，其算法效率虽然远高于一般的快速排序，但是在减低算法时间复杂度的同时提高了空间复杂度。

Mahafzah B A［14］使用pthread 并行函数库实现了一种并行快速排序算法，pthread 并行函数库会在程序启动时创建一束线程，随后将工作分配到线程束上。然而这种方式需要相当多的线程指定代码，而且不能保证线程随着可用处理器的数量而合理地进行扩充。相对地，OpenMP 不会将工作事先锁定在设定的线程数量中，修改代码中的配置信息即可。

Rashid L［15］等致力于优化并行过程中产生的数据依赖对性能的影响，他们所提出的算法虽然将快速排序的性能提升了4.17倍，但在测试过程中并没有考虑处理器的超线程技术（Hyper-Threading，HT）对于算法执行效率的影响。

3 并行快速排序算法

下面将对本文中提出的基于划分与归并的并行快速排序算法进行详细介绍，详细介绍了该算法的实现过程。本文中提出的快速排序算法的基本思想是将一个未排序的原始数组经过划分操作划分为多个部分有序数组，这些部分有序数组随后会经过标准库快速排序函数qsort 以OpenMP 任务的形式分别进行独立排序，直至最后原数组整体有序。基于划分与归并的并行快速排序算法的框图如图3所示。

图3 基于划分与归并的并行快速排序算法框图

3.1 划分函数详细设计

初期，一个未经排序的数组a=a0，a1，…an-1，在基准点（pivot）位置被划分为两个独立的子数组，基准点位置由MedianOfThree（）函数所确定，而不是像一般快速排序函数一样简单地使用第一个元素作为基准点，也没有采用随机选取基准点的方式，因为这种方式在最坏的情况下，也可能导致快速排序的时间复杂度降低至O（n2）。三数取中法（MedianOfThree）的基本做法是使用最左端、最右端和中心位置上的三个元素的中值作为枢纽元，三数取中法实现代码如下：

unsigned int mid=low+（（high-low）＞＞1）；

if（arr［mid］＞arr［high］）｛

SWAP（arr［mid］，arr［high］）；｝

if（arr［low］＞arr［high］）｛

SWAP（arr［low］，arr［high］）；｝

if（arr［low］＞arr［mid］）｛

SWAP（arr［mid］，arr［low］）；｝

return mid；

首先，为取到中心位置下标采用了右移位运算实现，这是因为在大多数处理器上，执行除法指令速度较慢，需要30 个时钟周期左右，而移位操作只需要1 个或几个时钟周期。而且在对大样本排序的过程中，需要对大量的子数组进行三数取中操作，故采用移位操作将大大提升程序执行效率。其次，本文所提出算法的实现代码中所有的交换操作均采用了宏定义函数的方式进行实现，即#define SWAP（a，b）｛a=a+b；b=a-b；a=a-b｝，这是由于交换操作在算法实现的过程中大量使用，如果采用普通函数的方式来实现的话，每当程序流程执行到交换函数时，为了切换并保护现场，交换函数就会被压入系统栈，而执行结束之后又会被弹出系统栈，处理器堆栈操作将会被反复执行，而宏函数只是在源代码文件上的文本替换，没有一个切换并保护现场的过程，因此可显著提高代码的执行效率。

在确定了基准点位置p 之后，就划分出了基准点左侧的左子数组和基准点右侧的右子数组。随后左右子数组将在双线程下同时进行一般的划分操作，使用相同的基准点ap。当左右子数组划分操作结束以后将返回划分结束时的索引位置，最终的原数组将被基准点和左右结束索引划分为四个子数组，具体代码实现如下：

p=MedianOfThree（arr，iL，jR）；

unsigned int d=p；

iR=p-1；

jL=p+1；

#pragma omp parallel num_threads（2）private（d）

｛

#pragma omp task

id=partition（arr，d，iL，iR）；

#pragma omp task

jd=partition（arr，d，jL，jR）；

｝

iL=id+1；

jR=jd-1；

其中iL，iR jL，jR 分别代表了左右子数组的最左端、最右端下标，局部变量d 则保存了基准点位置p变量的副本，使用pragma 编译器指示字指示编译器执行OpenMP 并行编程，总线程数为2，因为d在两个线程之间需要进行互斥访问，所以将其设置为线程私有化变量。在两个OpenMP 任务中将左右两个子数组划分为四个子数组，id，jd 保存划分结束的位置下标，用于左右子数组的再划分。

3.2 归并函数详细设计

首先，划分操作结束后，根据快速排序函数划分操作的特点可知第二部分子数组（大于部分）和第三部分子数组（小于等于部分）的长度是未必相等的，可能出现长度大于、小于、等于三类情况，因此在进行交换操作之前，应先对两个字数组的长度进行比较，根据比较结果进行不同的处理。由于长度计算操作频繁出现，因此本文所提出的算法中，计算数组长度操作同样采用宏函数LEN实现，用于提高代码执行效率，代码实现如下。

#define LEN（x，y）（（y）-（x）+1）

其中x，y为子数组的两端下标，y ≥x。以其中一种情况为例，当第二部分子数组较短时，应首先将需要交换的长度length 设置为较短子数组的长度，将较长子数组中超出长度length 部分元素中下标最大的元素同基准点元素进行数据交换，并更新基准点位置变量，再使用一个临时局部变量temp存储更新后的第三部分子数组的最左端位置。随后将循环地进行交换操作，从第二部分子数组最左端位置和temp位置同时开始，循环进行length次数据交换即可完成两子数组的交换操作。其它两种情况以此类推，实现代码如下：

unsigned int temp=0；

unsigned int length=0；

if（LEN（iL，p-1）＜LEN（p+1，jR））｛

length=LEN（iL，p-1）；

SWAP（arr［p］，arr［jR-length］）；

temp=p+1；

｝else if（LEN（iL，p-1）＞LEN（p+1，jR））｛

temp=p+1；

length=LEN（p+1，jR）；

SWAP（arr［p］，arr［iL+length］）；

p=iL+length；

｝else｛

temp=p+1；

length=LEN（p+1，jR）；

｝

unsigned int i=0；

#pragma omp parallel for simd num_threads（2）

for（i=0；i ＜length；i++）｛

SWAP（arr［iL+i］，arr［temp+i］）；

｝

其中，在循环操作上使用了simd 指令，此时循环执行过程中，多次迭代过程将会使用处理器的simd 指令进行并行执行，如果处理器不支持SIMD指令集，则编译器会将此指令作为注释处理，不会影响程序的正常执行。

3.3 Qsort函数

该函数递归地实现了对排序过程中的每个子数组进行划分与归并操作，递归结束的条件为每个分块的大小已经小于切割长度变量u，之所以定义切割长度变量u，目的是为了防止过度划分而影响排序操作的整体性能，根据快速排序原理可知，对于一个长度不大且进行过划分操作的数组再使用标准库快速排序函数进行排序，速度提升明显。Qsort函数实现代码如下。

if（iL+u ＜jR）｛

p=ParallelPartition（arr，iL，jR）；

#pragma omp task

Qsort（arr，iL，p-1，u）；

#pragma omp task

Qsort（arr，p+1，jR，u）；

｝else｛

qsort（arr，LEN（iL，jR），sizeof（DATA_TYPE），compare）；｝

其中ParallelPartition 仅对划分操作与归并操作依次调用，因此不再过多赘述，由以上代码可知当一个划分块的长度小于u 时，程序会递归地对划分块进行划分与归并操作，直至划分块长度大于等于u 时，退出递归，对划分块使用标准库快速排序函数进行排序。

4 实验设计

4.1 硬件选择

本文中所设计的实验是将算法在三款处理器上进行测试，这三款处理器分别为Intel Core i7-4790、Intel Core i7-2630QM 以 及Intel Core 2 Quad Q9400，表1 对这些多核处理器的详细参数进行了说明。

表1 多核处理器详细参数表

4.2 实验平台详细设计

4.2.1 样本生成详细设计

在本文所设计的实验中，主要对三种数据类型的大样本数据进行排序，这三种数据类型分别是unsigned int，unsigned long long 和double，即无符号32 位整型数据，无符号64 位整型数据以及64 位浮点型数据，数据均由C++语言编写的生成函数生成得来。

值得一提的是，数据生成函数可以生成三种类型的数据，仅需要经过简单的宏常量配置即可实现切换，且其中任意一种类型的数据在本文提出的算法中均有三种表示形式，下表2 分别列出了三种数据类型的三种表示形式。

表2 基本数据类型表示法

不仅出现在随机数生成函数中，在本文提及的所有函数中几乎都实现了使用宏常量切换测试数据类型的设计方式，这种设计方式大大提升了算法的通用性。实验测试数据的数据类型之所以采用生成类型和生成码同时描述是因为，在某些情况下，仅使用单一形式宏常量描述，在编译预处理语句中使用可能会产生问题。

4.2.2 执行线程数与切割块长度控制

本文所提出的算法中，执行线程数以及大样本切割分块长度都可以使用宏常量进行设置。

1）执行线程数控制说明

OpenMP 规范定义了用于控制OpenMP 程序执行线程数的环境变量OMP_NUM_THREADS，以及控制启用或禁用线程数动态调整的环境变量OMP_DYNAMIC，这两个环境变量均可通过显式的库函数omp_set_num_threads 以及omp_set_dynamic对其值进行动态调整。在默认条件下，线程数动态调整被定义为禁用，但由于本文所设计的实验需要在不同的线程数条件下进行，需要手动启动线程数动态调整机制。

2）切割块长度控制说明

本文所设计的实验对大样本切割分块长度进行设置的目的，是为了使算法在不同的多核处理器架构下都可以对性能产生积极影响，当切割分块长度小于等于高速缓存（cache）容量时，算法运行过程中的递归和迭代过程的性能就会得到大大的提升。所谓cache是低容量、高速度的存储器，通常集成在CPU内部，由于主存的数据传输速度远远低于CPU 运算速度，所以高速缓存就显得必不可少，因此在大多数情况下，高速缓存能够缓解主存带来的影响。现代处理器通常设有至少两级cache，分别被称作L1 cache 和L2 cache。L1 cache 在一般情况下被分为两个部分，即指令cache（L1I）和数据cache（L1D）。而L2 cache 并未做此类区分，存储数据和存储指令通常是相同的。当CPU 发出读请求（加载）时，传输一个数据项到寄存器中，L1 cache逻辑检查该数据项是否已经存放在其中。如果在，则称为cache命中，且请求能够被立即响应，然而在cache未命中的情况下，数据需要从外层cache中取出，最糟糕的情况下需要从主存中获得。如果所有的cache 都被占用，则需要有硬件实现的算法来进行数据替换，用新的数据替换cache 中原有的数据。实验中所使用的三种多核处理器的架构示意图如下图4所示。

图4 Haswell和Sandy Bridge处理器架构示意图

图5 Yorkfield处理器架构示意图

在算法执行过程中，为了提高程序的执行效率，使大样本切割分块长度小于等于cache 容量便可以减少仿存次数，这对于性能的提高有着至关重要的作用。对于图4 所示的架构，应使得分块长度小于L2 cache，这是因为L1 cache 是CPU 的第一层高速缓存，它距离CPU寄存器最近，带宽最大，延迟最低，管理开销最小，不过高速缓存均由静态RAM组成，结构较复杂，在CPU 管芯面积不能太大的情况下，L1 级高速缓存的容量不可能做得太大，因此不能随L1 cache 将大样本切割为极小的分块。而L3 cache 虽然容量极大，但被四个处理器核心所共享，在算法执行的过程中自然会产生大量的额外开销用于四个处理器核之间的通信，因此也不应将分块大小以L3 cache作为标尺。而对于图5所示的架构，是由两个双核处理器组成的四核处理器芯片。每个双核处理器具有共享的L2 cache 和独立的L1 cache，共有两组双核L2高速缓存，因此应使得分块长度小于L2 cache，因为虽然双核处理器共享的L2 cache 也是存在通信开销的，但是要比四核共享cache 的开销小的多。至于不能随L1 cache 将大样本切割为极小的分块，理由同上。

4.2.3 实验平台其他设计

实验进行于Microsoft Windows 10 操作系统之上，使用集成开发环境Dev-C++内嵌的TDM-GCC 4.9.2 64-bit Release 编译器对算法代码进行了编译和连接，为了使算法可以支持SIMD指令集，代码中使用了OpenMP4.0 并行函数库。算法性能的评估方式是将其与串行标准库快速排序函数运行时间进行对比，从而得到加速比。实验数据集包括三种类型的数据，它们分别是无符号32 位整型（Uint32），无符号64 位整型（Uint64）以及双精度浮点型（Double）数据，这些数据集在开始排序前都已经由生成程序生成为了数据集文件，因此可以保证每种算法都可以对相同的数据集进行排序。所有的实验参数详见下表3。

表3 实验参数集

表3 中将所有平台下实验的切割块大小都定义为200K，这样做的目的，一方面是为了在硬件方面可以兼容所有的处理器缓存，另一方面也是为结果分析提前控制了变量。此外，线程数为1 的实验指的是标准库快速排序函数实验，它与切割块的大小是无关的，仅仅为了对本文实现的快速排序进行性能评估。

5 实验验证及结果分析

这一部分详细地描述了在不同的实验参数下，排序算法的并行执行时间、串行执行时间、加速比等，并对实验结果进行了相关的统计分析，各种处理器平台下的实验结果如下表4所示。

三组实验结果包含了本文所提出的算法运行在三种处理器平台上的并行执行时间，串行执行时间以及加速比。其中加速比是由串行执行时间除以并行执行时间所得的结果，它代表了本文所提出的算法对于标准库串行快速排序性能的提升。

表4 Intel Core i7-4790处理器平台实验结果

表5 Intel Core i7-2630QM处理器平台实验结果

表6 Intel Core 2 Quad Q9400处理器平台实验结果

5.1 加速比

对比三组数据，不难发现Uint32数据类型的加速比一般优于Uint64和Double数据类型，而在三种处理器平台中，除Intel Core 2 Quad Q9400 不支持超线程技术（即Q9400虽有四核心但仅有四个物理线程同时执行）外，其余两种处理器平台均支持超线程技术（即四核心八线程），因此它们的加速比一般优于非超线程技术处理器。当然也有例外情况产生，这是由于OpenMP 并行函数库OMP_DYNAMIC 环境变量的定义所导致的，定义该环境变量之后就允许了对执行线程数的动态调整，在任务调度的过程中，会动态地将任务划分为多个部分，每个任务将不可中断地自始至终执行下去，同一个任务的不同部分会根据当前系统资源使用情况受到调度而同步地执行下去，而同时执行的线程个数以及任务不同部分的执行顺序是不得而知的。使用void omp_set_dynamic（int dynamic_threads）函数可以覆盖OMP_DYNAMIC 环境变量的定义，一旦参数dynamic_threads被定义为非零值，当前任务的动态调整即被启用了，此时线程组线程数量根据系统的资源状态动态调整，因此出现了加速比提升的例外情况。

对每组数据中的加速比按数据类型进行横向对比，不难发现加速比随着线程数的提高而不断提高，即线程数越高，算法的性能提升越明显。但在某些特殊情况下出现了加速比小于1 的情况，即算法执行效率低于同条件下标准库快速排序函数执行效率，经观察发现，这种情况大多数出现在定义的执行线程数低于处理器逻辑线程数的情况下，究其原因可总结为两点：1）由于Intel Core 2 Quad Q9400 处理器自身的特点（即两个处理器核心共享同一个L2 cache），故在双线程情况下，如果其中一个线程清空了cache 中的内容，那么另一个线程如果需要访问原cache 中的内容，只能等待到处理器将数据项再次加载到cache 中后，才能继续执行操作。2）本文中所提出的算法使用了SIMD指令处理循环部分操作，但是使用SIMD 指令并不总能带来性能的提升，使用SIMD指令，只会大大加快寄存器到寄存器操作的性能，但是会大大延长寄存器从内存子系统中获取新数据的时间，这一情况在数据长度过高（Uint64、Double数据类型长度均为64位）且执行线程数过低的情形下显得尤为突出。

5.2 并行执行时间

在三种处理器平台下，并行执行时间均与线程数的递增成反比，即并行执行线程数越高，并行执行时间反而越低，执行效率越高。随着线程数的增高，程序并行执行时间近似线性，算法可扩展性良好。以Uint32数据类型为例，三种处理器平台下算法的并行执行时间如下图6所示。

对每组数据中的加速比按并行执行时间进行纵向对比可以发现，并行执行时间与处理器主频有关，但在Intel Core i7-2630QM 与Intel Core 2 Quad Q9400 两处理器平台上在数据类型为Uint32 时又出现了相差不大的情况，综合考虑各方面因素，出现这一情况可能由于英特尔睿频加速技术在算法执行过程中的影响，英特尔睿频加速技术是英特尔酷睿i7 处理器和英特尔酷睿i5 处理器的独有特性。该技术可以智能地加快处理器速度，从而为高负载任务提供最佳性能——即最大限度地有效提升性能以匹配工作负载。CPU 会确定当前工作功率、电流和温度是否已达到最高极限，如仍有多余空间，CPU 会逐渐提高活动内核的频率，以进一步提高当前任务的处理速度，当程序只用到其中的某些核心时，CPU 会自动关闭其它未使用的核心，睿频加速技术无需用户干预，自动实现。

图6 三种处理器平台下算法的并行执行时间

6 结语

为了提升排序操作的性能，使用灵活的OpenMP 并行函数库以及C/C++语言标准库中提供的快速排序函数qsort 实现了一种可以运行于任意共享存储多核计算机上的并行快速排序算法。在循环操作上使用了simd 指令，循环执行过程中，多次迭代过程使用处理器的simd 指令进行并行执行，缩短了算法执行时间，提高了算法执行效率。并且在实验设计中，使大样本切割分块长度小于等于cache 容量减少了仿存次数，并行执行线程数量越多，速度提升越明显，提高了程序的执行效率。