使用 openmp 有效计算向量元素之间的交互

Question

我需要计算对象向量中所有元素i,j 之间的交互。在大小为N 的向量中，这相当于(N*(N-1))/2 的计算，并且会天真地在嵌套的 for 循环中解决，如下所示：

for( unsigned int i = 0; i < vector.size()-1; i++ ) {
  for ( unsigned int j = i+1; j < vector.size(); j++ ) {
    //compute interaction between vector[i] and vector[j]
  }
}

困难在于尝试使用 OpenMP 并行化加速该过程。随着i 的增加，内部循环中的计算次数线性减少。据我了解，#pragma omp parallel for 将循环除以使用的线程数。尽管外部循环会被平均划分，但实际计算不会。例如，长度为 257 的向量将进行 (257*256)/2=32896 次计算。如果 OpenMP 平均分割外循环（线程 1 的 i=0...127，线程 2 的 i=128...255），线程 1 必须计算 24640 次交互，而线程 2 必须计算 8256 次交互，取大约 75% 的时间，总效率为 62%。在 4 个线程之间拆分外循环将花费约 44% 的时间，效率约为 57%。我可以验证这是 MCVE 的问题

#include <iostream>
#include <unistd.h>
#include <omp.h>
#include <vector>
#include <ctime>

int main()
{
  timespec sleepTime;
  sleepTime.tv_sec = 0;
  sleepTime.tv_nsec = 1e6; // 1 ms                             
  std::vector< int > dummyVector(257,0);
  #pragma omp parallel for
  for(unsigned int i = 0; i < dummyVector.size()-1; i++ ) {
    for(unsigned int j = i+1; j < dummyVector.size(); j++ ) {
      // calculate( dummyVector[i], dummyVector[j] );
      nanosleep(&sleepTime,NULL);
    }
  }
  return 0;
}

使用 nanosleep 模拟我的交互，2 线程和 4 线程版本分别耗时 75% 和 44%

[me@localhost build]$ export OMP_NUM_THREADS=1
[me@localhost build]$ time ./Temp

real    0m38.242s ...
[me@localhost build]$ export OMP_NUM_THREADS=2
[me@localhost build]$ time ./Temp

real    0m28.576s ...
[me@localhost build]$ export OMP_NUM_THREADS=4
[me@localhost build]$ time ./Temp

real    0m16.715s ... 

如何更好地平衡线程间的计算？有没有办法告诉 OpenMP 不连续地拆分外循环？

---

为了将嵌套的 for 循环移出 omp 并行块，我尝试预先计算所有可能的索引对，然后遍历这些对

  std::vector< std::pair < int, int > > allPairs;
  allPairs.reserve((dummyVector.size()*(dummyVector.size()-1))/2);
  for(unsigned int i = 0; i < dummyVector.size()-1; i++ ) {
    for(unsigned int j = i+1; j < dummyVector.size(); j++ ) {
      allPairs.push_back(std::make_pair<int,int>(i,j));
    }
  }

  #pragma omp parallel for
  for( unsigned int i = 0; i < allPairs.size(); i++ ) {
    // calculate( dummyVector[allPairs[i].first], 
    //   dummyVector[allPairs[i].second] ); 
    nanosleep(&sleepTime,NULL);
  }

这确实有效地平衡了跨线程的计算，但它引入了索引对的不可避免的串行构造，随着N 的增长，这将损害我的运行时间。我还能做得更好吗？

Answer 1

正如@HighPerformanceMark 所建议的，解决方案在于调度 OpenMP 并行 for 循环。 Lawrence Livermore OpenMP tutorial 对不同的选项有很好的描述，但一般语法是 #pragma parallel for schedule(type[,chunk])，其中块参数是可选的。如果您未指定时间表，则默认值是特定于实现的。对于libgomp，默认是STATIC，将循环迭代均匀连续地划分，导致这个问题的负载均衡很差。

另外两个调度选项以稍高的开销为代价解决了负载平衡问题。第一个是 DYNAMIC，它在线程完成之前的工作时动态地为每个线程分配一个块（默认块大小为 1 次循环迭代）。因此代码看起来像这样

#pragma omp parallel for schedule( dynamic )
  for(unsigned int i = 0; i < dummyVector.size()-1; i++ ) {
    for(unsigned int j = i+1; j < dummyVector.size(); j++ ) {
      // calculate( dummyVector[i], dummyVector[j]);
    }
  }

因为内部循环的计算成本是结构化的（随着i 的增加而线性减少），所以 GUIDED 计划也很有效。它还动态地为每个线程分配工作块，但它从更大的块开始，并随着计算的继续而减小块大小。分配给线程的第一个迭代块的大小为number_iterations/number_threads，每个后续块的大小为remaining_iterations/number_threads。但是，这确实需要颠倒外部循环的顺序，以便初始迭代包含最少的工作量。

#pragma omp parallel for schedule( guided )
  for(unsigned int i = dummyVector.size()-1; i > 0; i-- ) {
    for(unsigned int j = i; j < dummyVector.size(); j++ ) {
      // calculate( dummyVector[i], dummyVector[j] ); 
    }
  }

Answer 2

我建议查看其他答案（使用 omp parallel for 的调度指令）。

还有一种选择，有一种可能是计算一个线性索引，然后从中提取i和j，如下：

如果您必须同时计算 i->j 和 j->i 交互：

#pragma omp parallel for
for(size_t u=0; u<vector.size()*vector.size(); ++u) {
   size_t i = u/vector.size();
   size_t j = u%vector.size();
   if(i != j) {
      compute interaction between vector[i] and vector[j]
   }
}

如果交互是对称的（如您的情况）：

也有类似的公式，但难度更大。您需要生成以下 (i,j) 对的序列：

(1,0) (2,0) (2,1) (3,0) (3,1) (3,2) (4,0) (4,1) (4,2) (4 ,3) ...

对于索引 i，对 (i,j) 的关联序列的长度为 i，因此将一对 (i,j) 转换为线性索引 u 的公式为：

u = i(i-1)/2 + j

现在需要“反转”这个公式，并检索整数 i 和 j

当 j=0 时检索 i 的值：

i^2 - i - 2*u = 0

解二次方程得到：

i = (1 + (int)(sqrt(1 + 8*u))) / 2

推导出j的值：

j = u - i*(i-1)/2

是的，这样做的方式非常复杂，我绝对更喜欢其他不易出错的解决方案！

edit1： 将 u,i,j 声明为 size_t（而不是 int）以避免溢出（如果 vector.size() 大于 2 的 32 次方，仍然会发生，但这留下了合理的空间） .感谢 EOF 的发言。

edit2：如果交互是对称的，那么它比我想象的更微妙（我最初的建议相当于负载不平衡的嵌套循环，请参阅 cmets）。

edit3：反转对称交互情况的公式