OpenMP:并行快速排序



我尝试使用OpenMP在分区部分和快速排序部分并行化快速排序。我的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("%dt", a[i]);
    printf("n");
    return 0;
}

对于main函数中的示例,排序结果为:

0   9   9   2   2   2   6   7   7   7

我使用gdb进行调试。在早期的递归中,一切进展顺利。但是在一些递归中,它突然在开始重复元素时出错了。然后生成上述结果。

谁能帮我找出问题在哪里吗?

我决定张贴这个答案是因为:

  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. 尽管如此,即使它是正确的,这个问题的现代答案是使用OpenMP tasks而不是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

我为我的第一个评论感到抱歉。这与你的问题无关。我还没有找到你问题的真正问题(也许你的移动元素有问题)。根据你的意见,我写了一个类似的程序,它有效很好。(我也是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("%dt", a[i]);
    }
    printf("n");
    return 0;
}

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

从本质上讲,正在发生的事情是您已经创建了一个数组 N 在共享内存中,当你做一个#pragma omp parallel sections,产生尽可能多的工作线程有#pragma omp section’s。每次工作线程试图访问和修改的元素时,它将执行一系列指令:"读的N值 N 从给定的地址","修改的N值 N ","写的N值 N 回到给定地址"。因为你有多个没有锁或同步的线程,读、修改和写指令可以由多个处理器以任意顺序执行,所以线程可以覆盖彼此的修改或读取未更新的值。

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

顺便说一下,最好的并行排序算法之一是PSRS(定期抽样并行排序),它使进程之间的列表大小更加平衡,进程之间没有不必要的密钥通信,并且可以在任意数量的并发进程上工作(它们不一定是2的幂)。

最新更新