使用OpenMP计算直方图

Question

我希望能够将这段代码并行化以获得最佳性能。"histogram" 存储了特定颜色出现的次数（有 10 种不同的颜色，因此直方图的大小为 10）。"img" 是一个存储某些图像信息的数组。在 img 的每个索引中都存储着一种颜色（int 值，范围为 0..9）。以下是该代码：

for( i=0; i<N1; i++ ){
  for( j=0; j<N2; j++ ){
    histogram[ img[i][j] ]  = histogram[ img[i][j] ] + 1;
  }
}

我尝试过这个方法，但性能非常差（比串行执行还要差）：

#pragma omp parallel for schedule(static, N1/nthreads) private(i,j)
for(i=0; i<N1; i++){
  for(j=0; j<N2; j++)
  {
    #pragma omp atomic
    histogram[img[i][j]]++;
  }
}

有什么建议吗？谢谢。

- Ortzi

Answer 1

使用OpenMP 4.5，您可以在C语言中简单地使用数组归约：

int histogram[BINS] = {0};
#pragma omp parallel for reduction(+:hist)
for(i=0; i<N1; i++) {
   for(j=0; j<N2; j++) {    
       histogram[img[i][j]]++;
   }
}

- Zulan

---

这可能是目前最好的答案，因为它完全反映了您想要实现的内容，并让编译器/运行库来完成艰苦的工作。更少出错，可能更有效率。 - Jorge Bellon

@JorgeBellón 关于所提出的编辑。您知道是否实际上需要该部分规范吗？我只测试了数组和GCC 7可以正常工作。标准说：“如果列表项是数组或数组部分，...”，但也说：“C / C ++：在减少子句中出现的列表项可能包括数组部分。” 所以我不确定百分之百。此外，如果编译器知道长度，hist[:]应该与hist[:BINS]相同，对吗？ - Zulan

嗨@Zalan。你可能是完全正确的。我试图找到规范的那部分，但恐怕我无法做到。如果您能向我引用该章节，我将不胜感激。我记得不允许使用指针，但我在某个地方读到他们在4.5中更改了这个限制，以允许整个数组被缩小。 - Jorge Bellon

@JorgeBellón 我指的是第205页的第14行和2.4节。 - Zulan

Answer 2

我已经详细介绍了如何进行这项工作，您可以点击此处查看： 使用OpenMP在不使用关键段的情况下以并行方式填充直方图（数组约简）

这与数组的约简相同。在C/C++中，OpenMP没有内置支持此功能（但Fortran有），因此您需要自己完成它。

简单的解决方案是创建直方图的私有版本，在并行填充它们，然后在关键段中将它们合并为一个直方图。在您的情况下，您可以像这样操作：

int i, histogram[10];
for(i=0; i<10; i++) histogram[i] = 0;
#pragma omp parallel
{
    int i, j, histogram_private[10];
    for(i=0; i<10; i++) histogram_private[i] = 0;
    #pragma omp for nowait
    for(i=0; i<N1; i++) {
       for(j=0; j<N2; j++) {    
           histogram_private[img[i][j]]++;
       }
    }      
    #pragma omp critical 
    {
        for(i=0; i<10; i++) histogram[i] += histogram_private[i];
    }
}

同时进行合并也是可能的，但这更加复杂。有关详细信息，请参阅我提到的第一个链接。

- Z boson

---

我尝试了这个，但是出现了错误信息：“段错误（'core' generated)”。它发生在“#pragma omp for nowait”部分。 - Ortzi

@Ortzi。这对我来说很好用。您可以在http://coliru.stacked-crooked.com/a/c63d92bf5389fff8上查看结果，甚至可以编辑/编译代码。 - Z boson

尝试使用长度为5000的N1和N2。这是问题所在：http://coliru.stacked-crooked.com/a/56bdf5e3bf756cea - Ortzi

@Ortzi，我不能为你做所有的事情！你现在正试图在栈上分配一个5000x5000的int数组！那是95MB！你必须使用malloc来处理这样大的数组。 - Z boson

@Ortzi，这里，我使用malloc修复了大数组的问题。http://coliru.stacked-crooked.com/a/7e862724c4a66278 - Z boson

Answer 3

你想创建一种“缩减”，使每个线程都有自己的直方图数组，并且必须在第二个循环中合并所有的直方图.... 请参见下面的伪代码：

histogram = new int[256];
histogram_thread = new int[nbthread * 256];

#pragma omp parallel 
for(i=0; i<N1; i++){
  current_thread_id = omp_get_thread_num();
  for(j=0; j<N2; j++)
  {
    histogram_thread[current_thread_id*256 + img[i][j]]++;
  }
}

//merge
for(unsigned int ui = 0 ; ui < 256 ; ++ui)
{ 
   for(int t=0; t<nbthread ; ++t) 
   {
      histogram [i] += histogram_thread[t * 256 + i];
   }
}

delete [] histogram_thread;

- alexbuisson