【MPI、OpenMP】并行快速排序(C语言)

本文记录了使用MPI与OpenMP两种并行计算方法实现快速排序算法，题目是专业实验课上老师给的，主要分享一下自己的做法，希望大家不吝赐教（使用的语言是C语言，有例子+图阐述原理，代码注释很全）。

快速排序算法并行实现

1. 快速排序算法

本篇文章主要针对快排算法的并行实现，所以大家如果有忘记快排算法的，可以参考一下这篇文章（快速排序算法原理及实现（单轴快速排序、三向切分快速排序、双轴快速排序）），这篇文章里面有动图，感觉还是挺不错的，有需要的朋友可以看看，其次再贴一个C语言实现串行快速排序算法的程序（快速排序算法（QSort，快排）及C语言实现），并行算法的实现核心还是围绕这个方法进行的。

• 核心功能代码
  实现方式有挺多种的，我使用的是把数据data的第一个元素作为基准进行快排的，代码如下：

int Partition(int* data, int start, int end)   //划分数据
{
    int temp = data[start];   //以第一个元素为基准
    while (start < end) {
        while (start < end && data[end] >= temp)end--;   //找到第一个比基准小的数
        data[start] = data[end];
        while (start < end && data[start] <= temp)start++;    //找到第一个比基准大的数
        data[end] = data[start];
    }
    data[start] = temp;   //以基准作为分界线
    return start;
}

2. 并行快速排序原理

此处主要讲MPI并行方式的原理，OpenMP的就想办法凑一些可以并行的代码就行。

2.1 MPI快排并行原理

除了快排本身的核心代码，该算法的核心问题是如何处理各个进程之间的通信。

2.1.1 进程数为2的整数次方

下面以四个进程为例,数据（下图中的星号 * 表示未排序的数据，不同颜色表示不同进程）可以被分发m = log₂4 = 2次，分发原则为：

• 0进程首先分发数据给 0+2^(m−1) = 2 进程，令m = m-1 = 1；
• 进而0进程剩余的数据分发给 0+2^(m−1) = 1 进程；同理2进程将从进程0接收的数据分发给 2+2^(m−1) = 3 进程；
• 无进程可以分配了，则四个进程对各自数据进行快排（排完序的数据用▲表示）；
• 最后每个进程将快排的结果传给分配给自己数据的进程（此处：1给0，3给2，2给0，0输出 ）。

下图为分发数据与发送结果的过程

2.1.2 进程数为2的非整数次方

当总进程数不是2的整数次方时，处理逻辑与上面稍有不同。
再以总进程数为6为例，数据可以被分发m = log₂6 = 3次(向上取整)，分发原则为：

• 0进程首先分发数据给 0+2^(m−1) = 4 进程，令m=m-1=2;
• 0进程剩余的数据分发给 0+2^(m−1) = 2 进程； 4进程将从进程0接收的数据分发给 4+2^(m−1) = 6 ，但是6进程不存在，因此在4进程内，使用循环一直令m=m-1，直到 4+2^(m−1) 进程是存在的，此处需要将m减到1，才满足 4+2^(m−1) = 5，5进程存在，所以4进程将数据分给5进程;
• 而0到3这四个进程的数据分配，和上一个例子仅有四个进程的原则一致;
• 无进程可以分配了，则六个进程对各自数据进行快排;
• 最后每个进程将快排的结果传给分配给自己数据的进程 （此处：1给0，3给2，2给0；5给4，4给0，0输出）。

下图为分发数据的过程

下图为发送结果的过程

2.2 OpenMP快排并行原理

omp的原理也不难理解，凑出可并行的区域即可，当然也要分析题目，先从理论上判断是否适合多线程并行。
下面代码的思路就是使用omp的sections设置并行域，需要注意的是，设置sections的时候，如果没有添加parallel，表示sections中的几个section是串行执行的；而加上parallel后表示每个section内部是串行的，而section之间是并行的（可以参考文章：OpenMp之sections用法）。

void quickSort(int* data, int start, int end)  //并行快排
{
    if (start < end) {
        int pos = Partition(data, start, end);
        #pragma omp parallel sections    //设置并行区域
        {
            #pragma omp section          //该区域对前部分数据进行排序
            quickSort(data, start, pos - 1);
            #pragma omp section          //该区域对后部分数据进行排序
            quickSort(data, pos + 1, end);
        }
    }
}

3. 完整代码与结果（含注释）

3.1 MPI

• 代码

#include"mpi.h"
#include<time.h>
#include <stdlib.h>
#include <stdio.h>

int Partition(int* data, int start, int end)   //划分数据
{
    int temp = data[start];   //以第一个元素为基准
    while (start < end) {
        while (start < end && data[end] >= temp)end--;   //找到第一个比基准小的数
        data[start] = data[end];
        while (start < end && data[start] <= temp)start++;    //找到第一个比基准大的数
        data[end] = data[start];
    }
    data[start] = temp;   //以基准作为分界线
    return start;
}

void QuickSort(int* data, int start, int end)  //串行快排
{
    if (start < end) {    //未划分完
        int r = Partition(data, start, end);   //继续划分，进行递归排序
        QuickSort(data, start, r - 1);
        QuickSort(data, r + 1, end);
    }
}

//求2的n次方
int exp2(int n)
{
    int i = 1;
    while (n-- > 0) i *= 2;
    return i;
}

//求以2为底n的对数，向下取整
int log2(int n)
{
    int i = 1, j = 2;
    while (j < n) {
        j *= 2;
        i++;
    }
    return i;
}

void paraQuickSort(int* data, int start, int end, int m, int id, int nowID, int N)
{
    int i, j, r = end, length = -1;  //r表示划分后数据前部分的末元素下标，length表示后部分数据的长度
    int* t;
    MPI_Status status;
    if (m == 0) {   //无进程可以调用
        if (nowID == id) QuickSort(data, start, end);
        return;
    }
    if (nowID == id) {    //当前进程是负责分发的
        while (id + exp2(m - 1) > N && m > 0) m--;   //寻找未分配数据的可用进程
        if (id + exp2(m - 1) < N) {  //还有未接收数据的进程，则划分数据
            r = Partition(data, start, end);
            length = end - r;
            MPI_Send(&length, 1, MPI_INT, id + exp2(m - 1), nowID, MPI_COMM_WORLD);
            if (length > 0)   //id进程将后部分数据发送给id+2^(m-1)进程
                MPI_Send(data + r + 1, length, MPI_INT, id + exp2(m - 1), nowID, MPI_COMM_WORLD);
        }
    }
    if (nowID == id + exp2(m - 1)) {    //当前进程是负责接收的
        MPI_Recv(&length, 1, MPI_INT, id, id, MPI_COMM_WORLD, &status);
        if (length > 0) {   //id+2^(m-1)进程从id进程接收后部分数据
            t = (int*)malloc(length * sizeof(int));
            if (t == 0) printf("Malloc memory error!");
            MPI_Recv(t, length, MPI_INT, id, id, MPI_COMM_WORLD, &status);
        }
    }
    j = r - 1 - start;
    MPI_Bcast(&j, 1, MPI_INT, id, MPI_COMM_WORLD);
    if (j > 0)     //负责分发的进程的数据不为空
        paraQuickSort(data, start, r - 1, m - 1, id, nowID, N);   //递归调用快排函数，对前部分数据进行排序
    j = length;
    MPI_Bcast(&j, 1, MPI_INT, id, MPI_COMM_WORLD);
    if (j > 0)     //负责接收的进程的数据不为空
        paraQuickSort(t, 0, length - 1, m - 1, id + exp2(m - 1), nowID, N);   //递归调用快排函数，对前部分数据进行排序
    if ((nowID == id + exp2(m - 1)) && (length > 0))     //id+2^(m-1)进程发送结果给id进程
        MPI_Send(t, length, MPI_INT, id, id + exp2(m - 1), MPI_COMM_WORLD);
    if ((nowID == id) && id + exp2(m - 1) < N && (length > 0))     //id进程接收id+2^(m-1)进程发送的结果
        MPI_Recv(data + r + 1, length, MPI_INT, id + exp2(m - 1), id + exp2(m - 1), MPI_COMM_WORLD, &status);
}

int main(int argc, char* argv[])
{
    int* data;
    int rank, size;
    int i, j, m, r, n = atoi(argv[1]);   //随机数组的长度
    double start_time, end_time;
    MPI_Status status;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);  //当前进程的进程号
    MPI_Comm_size(MPI_COMM_WORLD, &size);  //总进程数
    if (rank == 0) {   //根进程生成随机数组
        data = (int*)malloc(n * sizeof(int));
        srand(time(NULL) + rand());   //随机数种子
        for (i = 0; i < n; i++)
            data[i] = (int)rand();   //获取n个随机整数
    }
    m = log2(size);  //第一次分发需要给第2^(m-1)个进程
    start_time = MPI_Wtime();
    MPI_Bcast(&n, 1, MPI_INT, 0, MPI_COMM_WORLD);  //广播n
    paraQuickSort(data, 0, n - 1, m, 0, rank, size);  //执行快排
    end_time = MPI_Wtime();

    if (rank == 0) {   //根进程输出并行时间
        for (i = 0; i < n && i < 10; i++)   //n太大时只输出前10个
            printf("%d ", data[i]);
        printf("\n并行时间：%lfs\n", end_time - start_time);
    }
    MPI_Finalize();
}

• 结果
  

3.2 OpenMP

• 代码

#include<stdio.h>
#include<omp.h>
#include<time.h>
#include<stdlib.h>

int Partition(int* data, int start, int end)   //划分数据
{
    int temp = data[start];   //以第一个元素为基准
    while (start < end) {
        while (start < end && data[end] >= temp)end--;   //找到第一个比基准小的数
        data[start] = data[end];
        while (start < end && data[start] <= temp)start++;    //找到第一个比基准大的数
        data[end] = data[start];
    }
    data[start] = temp;   //以基准作为分界线
    return start;
}

void quickSort(int* data, int start, int end)  //并行快排
{
    if (start < end) {
        int pos = Partition(data, start, end);
        #pragma omp parallel sections    //设置并行区域
        {
            #pragma omp section          //该区域对前部分数据进行排序
            quickSort(data, start, pos - 1);
            #pragma omp section          //该区域对后部分数据进行排序
            quickSort(data, pos + 1, end);
        }
    }
}

int main(int argc, char* argv[])
{
    int n = atoi(argv[2]), i;   //线程数
    int size = atoi(argv[1]);   //数据大小
    int* num = (int*)malloc(sizeof(int) * size);

    double starttime = omp_get_wtime();
    srand(time(NULL) + rand());   //生成随机数组
    for (i = 0; i < size; i++)
        num[i] = rand();
    omp_set_num_threads(n);   //设置线程数
    quickSort(num, 0, size - 1);   //并行快排
    double endtime = omp_get_wtime();

    for (i = 0; i < 10 && i<size; i++)//输出前十个元素
        printf("%d ", num[i]);
    printf("\n并行时间：%lfs\n", endtime - starttime);
    return 0;
}

• 结果