我用c++做了一些科学计算，并尝试利用OpenMP来并行化一些循环。到目前为止，这在8线程的Intel i7-4770上运行得很好。

<标题> 设置

我们有一个小型工作站，由两个英特尔cpu (E5-2680v2)在一个主板上组成。代码工作，只要它运行在1个CPU与尽可能多的线程，我喜欢。但是一旦我使用了第二个CPU，我就会不时地观察到不正确的结果(大约每运行代码50 -100次)。即使我只使用2个线程并将它们分配给两个不同的cpu，也会发生这种情况。因为我们有5个这样的工作站(都是相同的)，我在每一个工作站上运行代码，都显示了这个问题。

工作站运行在OpenSuse 13.1，内核3.11.10-7。这个问题存在于g++ 4.8.1和4.9.0，以及Intel的icc 13.1.3.192(尽管icc的问题不经常发生，但它仍然存在)。

<标题> 症状

故障现象描述如下:

• 我有一个大数组std::complex: std::complex<double>* mFourierValues;
• 在循环中，我访问并设置每个元素。每次迭代访问一个不同的元素，所以我没有并发访问(我检查了这个):mFourierValues[idx] = newValue;
• 如果我将设置的数组值与之后的输入值进行比较，大致为mFourierValues[idx] == newValue，则此检查不时失败(尽管不是每次结果都不正确)。

所以症状看起来像我并发地访问元素而没有任何同步。但是，当我将索引存储在std::vector(具有适当的#pragma omp critical)中时，

所有索引必须唯一且在正确的范围内。<标题>

经过几天的调试，我越来越怀疑有别的事情发生了，我的代码是正确的。在我看来，当cpu与主存同步缓存时，似乎发生了一些奇怪的事情。

因此，我的问题是:

• OpenMP甚至可以用于这样的系统吗?(我还没有找到一个来源说没有)
• 是否有这种情况下已知的bug(我没有发现任何在bug跟踪器)?
• 你认为问题可能在哪里?
  • 我的代码(似乎运行良好的1个CPU多核!)，
  • 编译器(gcc, icc两个!)，
  • 操作系统，
  • 硬件(所有5个工作站的缺陷?)

<标题> 代码

[编辑:旧代码删除，见下文]

使用最小示例编辑

好了，我终于能够生成一个更短(并且自一致)的代码示例。

关于代码

• 预留一些内存空间。对于堆栈上的数组，可以这样访问:complex<double> mAllElements[tensorIdx][kappa1][kappa2][kappa3]。也就是说，我有3个秩3张量(tensorIdx)每个张量代表一个三维数组，由kappa1、kappa2和kappa3索引。
• 我有4个嵌套循环(超过所有4个索引)，而kappa1循环是一个得到并行化(并且是最外面的一个)。它们位于DoComputation()。
• 在main()中，我调用DoComputation()一次以获得一些参考值，然后我调用它几次并比较结果。它们应该完全匹配，但有时却不匹配。

不幸的是，代码仍然大约190行长。我试着进一步简化它(只有一个秩为1的张量，等等)，但我永远无法重现这个问题。我猜这似乎是因为内存访问是非对齐的(tensorIdx上的循环是最内层的)(我知道，这远非最佳。)

此外，在适当的地方需要一些延迟，以重现错误。这就是nops()调用的原因。没有它们，代码运行得快得多，但到目前为止还没有显示出问题。

请注意，我再次检查了关键部分CalcElementIdx()，并认为它是正确的(每个元素被访问一次)。我还运行了valgrind的memcheck, helgrind和drd(使用正确重新编译的libgome)，这三个都没有错误。

输出

每隔一秒到第三次启动程序，我就会得到一个或两个不匹配。示例输出:

41      Is exactly 0
42      Is exactly 0
43      Is exactly 0
44      Is exactly 0
45      348496
46      Is exactly 0
47      Is exactly 0
48      Is exactly 0
49      Is exactly 0

对于gcc和icc都是如此。

我的问题

我的问题是:下面的代码看起来正确吗?(除了明显的设计缺陷。)(如果它太长，我会尝试进一步减少它，但如上所述，到目前为止我失败了。)

使用

编译代码

g++ main.cc -O3 -Wall -Wextra -fopenmp

或

icc main.cc -O3 -Wall -Wextra -openmp

两个版本都显示了在2个cpu上运行总共40个线程时所描述的问题。我不能在1个CPU(和尽可能多的线程)上观察到这个bug。

// File: main.cc
#include <cmath>
#include <iostream>
#include <fstream>
#include <complex>
#include <cassert>
#include <iomanip>
#include <omp.h>
using namespace std;

// If defined: We add some nops in certain places, to get the timing right.
// Without them, I haven't observed the bug.
#define ENABLE_NOPS
// The size of each of the 3 tensors is: GRID_SIZE x GRID_SIZE x GRID_SIZE
static const int GRID_SIZE = 60;
//=============================================
// Produces several nops. Used to get correct "timings".
//----
template<int N> __attribute__((always_inline)) inline void nop()
{
    nop<N-1>();
    asm("nop;");
}
//----
template<> inline void nop<0>() { }
//----
__attribute__((always_inline)) inline void nops()
{
    nop<500>(); nop<500>(); nop<500>(); nop<500>(); nop<500>(); nop<500>(); nop<500>(); nop<500>(); nop<500>();
}


//=============================================
/*
Memory layout: We have 3 rank-3-tensors, i.e. 3 arrays of dimension 3.
The layout looks like this: complex<double> allElements[tensorIdx][kappa1][kappa2][kappa3];
The kappas represent the indices into a certain tensor, and are all in the interval [0; GRID_SIZE-1].
*/
class MemoryManagerFFTW
{
public:
  //---------- Constructor ----------
  MemoryManagerFFTW()
  {
    mAllElements = new complex<double>[GetTotalNumElements()];
  }
  //---------- Destructor ----------
  ~MemoryManagerFFTW() 
  { 
    delete[] mAllElements; 
  }
  //---------- SetElement ----------
  void SetElement(int tensorIdx, int kappa1, int kappa2, int kappa3, const complex<double>& newVal)
  {
    // Out-of-bounds error checks are done in this function.
    const size_t idx = CalcElementIdx(tensorIdx, kappa1, kappa2, kappa3);
    // These nops here are important to reproduce the bug.
#if defined(ENABLE_NOPS)
    nops();
    nops();
#endif
    // A flush makes the bug appear more often.
    // #pragma omp flush
    mAllElements[idx] = newVal;
    // This was never false, although the same check is false in DoComputation() from time to time.
    assert(newVal == mAllElements[idx]);
  }
  //---------- GetElement ----------
  const complex<double>& GetElement(int tensorIdx, int kappa1, int kappa2, int kappa3)const
  {  
    const size_t idx = CalcElementIdx(tensorIdx, kappa1, kappa2, kappa3);
    return mAllElements[idx];
  }

  //---------- CalcElementIdx ----------
  size_t CalcElementIdx(int tensorIdx, int kappa1, int kappa2, int kappa3)const
  {
    // We have 3 tensors (index by "tensorIdx"). Each tensor is of rank 3. In memory, they are placed behind each other.
    // tensorStartIdx is the index of the first element in the tensor.
    const size_t tensorStartIdx = GetNumElementsPerTensor() * tensorIdx;
    // Index of the element relative to the beginning of the tensor. A tensor is a 3dim. array of size GRID_SIZE x GRID_SIZE x GRID_SIZE
    const size_t idxInTensor = kappa3 + GRID_SIZE * (kappa2 + GRID_SIZE * kappa1);
    const size_t finalIdx = tensorStartIdx + idxInTensor;
    assert(finalIdx < GetTotalNumElements());
    return finalIdx;
  }

  //---------- GetNumElementsPerTensor & GetTotalNumElements ----------
  size_t GetNumElementsPerTensor()const { return GRID_SIZE * GRID_SIZE * GRID_SIZE; }
  size_t GetTotalNumElements()const { return NUM_TENSORS * GetNumElementsPerTensor(); }

public:
  static const int NUM_TENSORS = 3; // The number of tensors.
  complex<double>* mAllElements; // All tensors. An array [tensorIdx][kappa1][kappa2][kappa3]
};


//=============================================
void DoComputation(MemoryManagerFFTW& mSingleLayerManager)
{
  // Parallize outer loop.
  #pragma omp parallel for
  for (int kappa1 = 0; kappa1 < GRID_SIZE; ++kappa1)
  {
    for (int kappa2 = 0; kappa2 < GRID_SIZE; ++kappa2)
    {
      for (int kappa3 = 0; kappa3 < GRID_SIZE; ++kappa3)
      {    
#ifdef ENABLE_NOPS
        nop<50>();
#endif
        const double k2 = kappa1*kappa1 + kappa2*kappa2 + kappa3*kappa3;
        for (int j = 0; j < 3; ++j)
        {
          // Compute and set new result.
          const complex<double> curElement = mSingleLayerManager.GetElement(j, kappa1, kappa2, kappa3);
          const complex<double> newElement = exp(-k2) * k2 * curElement;
          mSingleLayerManager.SetElement(j, kappa1, kappa2, kappa3, newElement);
          // Check if the results has been set correctly. This is sometimes false, but _not_ always when the result is incorrect.
          const complex<double> test = mSingleLayerManager.GetElement(j, kappa1, kappa2, kappa3);
          if (test != newElement)
            printf("Failure: (%g, %g) != (%g, %g)n", test.real(), test.imag(), newElement.real(), newElement.imag());
        }
      }
    }
  }
}

//=============================================
int main()
{
  cout << "Max num. threads: " << omp_get_max_threads() << endl;
  // Call DoComputation() once to get a reference-array.
  MemoryManagerFFTW reference;
  for (size_t i = 0; i < reference.GetTotalNumElements(); ++i)
    reference.mAllElements[i] = complex<double>((double)i, (double)i+0.5);
  DoComputation(reference);
  // Call DoComputation() several times, and each time compare the result to the reference.
  const size_t NUM = 1000;
  for (size_t curTry = 0; curTry < NUM; ++curTry)
  {
    MemoryManagerFFTW mSingleLayerManager;
    for (size_t i = 0; i < mSingleLayerManager.GetTotalNumElements(); ++i)
      mSingleLayerManager.mAllElements[i] = complex<double>((double)i, (double)i+0.5);
    DoComputation(mSingleLayerManager);
    // Get the max. difference. This *should* be 0, but isn't from time to time.
    double maxDiff = -1;
    for (size_t i = 0; i < mSingleLayerManager.GetTotalNumElements(); ++i)
    {
      const complex<double> curDiff = mSingleLayerManager.mAllElements[i] - reference.mAllElements[i];
      maxDiff = max(maxDiff, max(curDiff.real(), curDiff.imag()));
    }
    if (maxDiff != 0)
      cout << curTry << "t" << maxDiff << endl;
    else
      cout << curTry << "t" << "Is exactly 0" << endl;
  }
  return 0;
}

<标题>编辑

从下面的评论和Zboson的回答可以看出，内核3.11.10-7有一个bug。在更新到3.15.0-1之后，我的问题消失了，代码正常工作。