c - OpenMP：并行快速排序

Question

我尝试使用 OpenMP 在分区部分和 QuickSort 部分中并行化 QuickSort。我的C代码如下：

#include "stdlib.h"
#include "stdio.h"
#include "omp.h"

// parallel partition
int ParPartition(int *a, int p, int r) {
    int b[r-p];
    int key = *(a+r); // use the last element in the array as the pivot
    int lt[r-p]; // mark 1 at the position where its element is smaller than the key, else 0
    int gt[r-p]; // mark 1 at the position where its element is bigger than the key, else 0
    int cnt_lt = 0; // count 1 in the lt array
    int cnt_gt = 0; // count 1 in the gt array
    int j=p;
    int k = 0; // the position of the pivot
    // deal with gt and lt array
    #pragma omp parallel for
    for ( j=p; j<r; ++j) {
        b[j-p] = *(a+j);
        if (*(a+j) < key) {
            lt[j-p] = 1;
            gt[j-p] = 0;
        } else {
            lt[j-p] = 0;
            gt[j-p] = 1;
        }
    }
    // calculate the new position of the elements
    for ( j=0; j<(r-p); ++j) {
        if (lt[j]) {
            ++cnt_lt;
            lt[j] = cnt_lt;
        } else
            lt[j] = cnt_lt;
        if (gt[j]) {
            ++cnt_gt;
            gt[j] = cnt_gt;
        } else
            gt[j] = cnt_gt;
    }
    // move the pivot
    k = lt[r-p-1];
    *(a+p+k) = key;
    // move elements to their new positon
    #pragma omp parallel for 
    for ( j=p; j<r; ++j) {
        if (b[j-p] < key)
            *(a+p+lt[j-p]-1) = b[j-p];
        else if (b[j-p] > key)
            *(a+k+gt[j-p]) = b[j-p];
    }
    return (k+p);
}

void ParQuickSort(int *a, int p, int r) {
    int q;
    if (p<r) {
        q = ParPartition(a, p, r);
        #pragma omp parallel sections
        {
        #pragma omp section
        ParQuickSort(a, p, q-1);
        #pragma omp section
        ParQuickSort(a, q+1, r);
        }
    }
}

int main() {
    int a[10] = {5, 3, 8, 4, 0, 9, 2, 1, 7, 6};
    ParQuickSort(a, 0, 9);
    int i=0;
    for (; i!=10; ++i)
        printf("%d\t", a[i]);
    printf("\n");
    return 0;
}

对于main函数中的例子，排序结果为：

0   9   9   2   2   2   6   7   7   7

我使用gdb进行调试。在早期的递归中，一切顺利。但是在某些递归中，它突然搞砸了，开始重复元素。然后生成上面的结果。

有人可以帮我找出问题所在吗？

Answer 1

我对我的第一个评论感到抱歉。与您的问题无关。我没有找到您问题的真正问题（可能您的移动元素有问题）。根据您的意见，我写了一个类似的程序，它工作正常.（我也是 OpenMP 的新手）。

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


int partition(int * a, int p, int r)
{
    int lt[r-p];
    int gt[r-p];
    int i;
    int j;
    int key = a[r];
    int lt_n = 0;
    int gt_n = 0;

#pragma omp parallel for
    for(i = p; i < r; i++){
        if(a[i] < a[r]){
            lt[lt_n++] = a[i];
        }else{
            gt[gt_n++] = a[i];
        }   
    }   

    for(i = 0; i < lt_n; i++){
        a[p + i] = lt[i];
    }   

    a[p + lt_n] = key;

    for(j = 0; j < gt_n; j++){
        a[p + lt_n + j + 1] = gt[j];
    }   

    return p + lt_n;
}

void quicksort(int * a, int p, int r)
{
    int div;

    if(p < r){ 
        div = partition(a, p, r); 
#pragma omp parallel sections
        {   
#pragma omp section
            quicksort(a, p, div - 1); 
#pragma omp section
            quicksort(a, div + 1, r); 

        }
    }
}

int main(void)
{
    int a[10] = {5, 3, 8, 4, 0, 9, 2, 1, 7, 6};
    int i;

    quicksort(a, 0, 9);

    for(i = 0;i < 10; i++){
        printf("%d\t", a[i]);
    }
    printf("\n");
    return 0;
}

Answer 2

我决定发布这个答案是因为：

1. 接受的答案是错误的，并且这些天用户似乎不活跃。有一个竞争条件

    #pragma omp parallel for
    for(i = p; i < r; i++){
        if(a[i] < a[r]){
            lt[lt_n++] = a[i]; //<- race condition lt_n is shared
        }else{
            gt[gt_n++] = a[i];  //<- race condition gt_n is shared
        }   
    }   
   

2. 尽管如此，即使它是正确的，这个问题的现代答案是使用 OpenMPtasks而不是sections.

我正在为社区提供这种方法的完整可运行示例，包括测试和分析。

#include <assert.h>
#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#include <omp.h>

#define TASK_SIZE 100

unsigned int rand_interval(unsigned int min, unsigned int max)
{
    // https://stackoverflow.com/questions/2509679/
    int r;
    const unsigned int range = 1 + max - min;
    const unsigned int buckets = RAND_MAX / range;
    const unsigned int limit = buckets * range;

    do
    {
        r = rand();
    } 
    while (r >= limit);

    return min + (r / buckets);
}

void fillupRandomly (int *m, int size, unsigned int min, unsigned int max){
    for (int i = 0; i < size; i++)
    m[i] = rand_interval(min, max);
} 


void init(int *a, int size){
   for(int i = 0; i < size; i++)
       a[i] = 0;
}

void printArray(int *a, int size){
   for(int i = 0; i < size; i++)
       printf("%d ", a[i]);
   printf("\n");
}

int isSorted(int *a, int size){
   for(int i = 0; i < size - 1; i++)
      if(a[i] > a[i + 1])
        return 0;
   return 1;
}


int partition(int * a, int p, int r)
{
    int lt[r-p];
    int gt[r-p];
    int i;
    int j;
    int key = a[r];
    int lt_n = 0;
    int gt_n = 0;

    for(i = p; i < r; i++){
        if(a[i] < a[r]){
            lt[lt_n++] = a[i];
        }else{
            gt[gt_n++] = a[i];
        }   
    }   

    for(i = 0; i < lt_n; i++){
        a[p + i] = lt[i];
    }   

    a[p + lt_n] = key;

    for(j = 0; j < gt_n; j++){
        a[p + lt_n + j + 1] = gt[j];
    }   

    return p + lt_n;
}

void quicksort(int * a, int p, int r)
{
    int div;

    if(p < r){ 
        div = partition(a, p, r); 
        #pragma omp task shared(a) if(r - p > TASK_SIZE) 
        quicksort(a, p, div - 1); 
        #pragma omp task shared(a) if(r - p > TASK_SIZE)
        quicksort(a, div + 1, r); 
    }
}

int main(int argc, char *argv[])
{
        srand(123456);
        int N  = (argc > 1) ? atoi(argv[1]) : 10;
        int print = (argc > 2) ? atoi(argv[2]) : 0;
        int numThreads = (argc > 3) ? atoi(argv[3]) : 2;
        int *X = malloc(N * sizeof(int));
        int *tmp = malloc(N * sizeof(int));

        omp_set_dynamic(0);              /** Explicitly disable dynamic teams **/
        omp_set_num_threads(numThreads); /** Use N threads for all parallel regions **/

         // Dealing with fail memory allocation
        if(!X || !tmp)
        { 
           if(X) free(X);
           if(tmp) free(tmp);
           return (EXIT_FAILURE);
        }

        fillupRandomly (X, N, 0, 5);

        double begin = omp_get_wtime();
        #pragma omp parallel
        {
            #pragma omp single
             quicksort(X, 0, N);
        }   
        double end = omp_get_wtime();
        printf("Time: %f (s) \n",end-begin);
    
        assert(1 == isSorted(X, N));

        if(print){
           printArray(X, N);
        }

        free(X);
        free(tmp);
        return (EXIT_SUCCESS);

    return 0;
}

如何运行：

该程序接受三个参数：

1. 数组的大小；
2. 是否打印数组，0为否，否则为是；
3. 并行运行的线程数。

迷你基准

在 4 核机器中：输入 100000

1 Thread  -> Time: 0.784504 (s)
2 Threads -> Time: 0.424008 (s) ~ speedup 1.85x
4 Threads -> Time: 0.282944 (s) ~ speedup 2.77x

Answer 3

我已经在生产环境中实现了并行快速排序，尽管使用并发进程（即 fork() 和 join()）而不是 OpenMP。我还找到了一个相当不错的 pthread 解决方案，但就最坏情况的运行时间而言，并发进程解决方案是最好的。首先让我说您似乎没有为每个线程制作输入数组的副本，因此您肯定会遇到可能破坏数据的竞争条件。

本质上，发生的事情是您在共享内存中创建了一个数组N#pragma omp parallel sections ，当您执行 a 时，您将产生与#pragma omp section's 一样多的工作线程。每次工作线程尝试访问和修改 a 的元素时，都会执行一系列指令：“从给定地址读取N的第 n 个值”、“修改N的第 n 个值”、“写入N的第 n 个值返回给定地址”。由于您有多个没有锁定或同步的线程，读取、修改和写入指令可以由多个处理器以任意顺序执行，因此线程可能会覆盖彼此的修改或读取未更新的值。

我发现的最佳解决方案（经过数周的测试和基准测试，我想出了许多解决方案）是将列表log(n)细分，其中n是处理器的数量。例如，如果您有一台四核机器 (n = 4)，则将列表细分 2 次 (log(4) = 2)，选择作为数据集中位数的枢轴。重要的是枢轴是中位数，否则您可能会遇到选择不当的枢轴导致列表在进程之间分布不均的情况。然后每个进程对其本地子数组进行快速排序，然后将其结果与其他进程的结果合并。这被称为“hyperquicksort”，从最初的 github 搜索中，我发现了这个. 我不能保证那里的代码，也不能发布我写的任何代码，因为它受 NDA 保护。

顺便说一句，最好的并行排序算法之一是 PSRS（常规采样并行排序），它使列表大小在进程之间更加平衡，不会在进程之间不必要地通信密钥，并且可以在任意数量的并发进程上工作（它们不一定是 2) 的幂。