我尝试使用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进行调试。在早期的递归中,一切进展顺利。但是在一些递归中,它突然在开始重复元素时出错了。然后生成上述结果。
谁能帮我找出问题在哪里吗?我决定张贴这个答案是因为:
-
接受的答案是错误的,用户最近似乎不活跃。
上有一个竞态条件#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 } }
-
尽管如此,即使它是正确的,这个问题的现代答案是使用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;
}
如何运行:
程序接受三个参数:
- 数组的大小;
- 打印或不打印数组,0表示否,否则是;
- 并行运行的线程数。
迷你基准
在4核机器中:用
输入1000001 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的幂)。