Я делаю некоторые научные вычисления на С++ и пытаюсь использовать OpenMP для параллелизации некоторых из циклов. До сих пор это хорошо работало, например. на Intel i7-4770 с 8 потоками.
У нас есть небольшая рабочая станция, которая состоит из двух процессоров Intel (E5-2680v2) на одной материнской плате. Код работает до тех пор, пока он работает на одном процессоре с таким количеством потоков, сколько мне нравится. Но как только я использую второй процессор, время от времени я наблюдаю неправильные результаты (примерно каждые 50-100-й раз я запускаю код). Это происходит даже тогда, когда я использую только 2 потока и назначаю их двум различным процессорам. Поскольку у нас есть 5 из этих рабочих станций (все одинаковы), я запускал код на каждом из них, и все это показывает эту проблему.
Рабочая станция работает на OpenSuse 13.1, ядро 3.11.10-7. Проблема существует с g++ 4.8.1 и 4.9.0 и с Intel icc 13.1.3.192 (хотя проблема не возникает часто с icc, но она по-прежнему существует).
Симптом может быть описан следующим образом:
std::complex<double>* mFourierValues;mFourierValues[idx] = newValue;mFourierValues[idx] == newValue, эта проверка будет время от времени отказываться (хотя не каждый раз, когда результаты заканчиваются некорректно).Итак, симптом выглядит так, что я обращаюсь к элементам одновременно без каких-либо синхронизаций. Однако, когда я сохраняю индексы в std::vector (с правильным #pragma omp critical),
все знаки уникальны и находятся в правильном диапазоне.
После нескольких дней отладки мое подозрение растет, что происходит еще что-то, и что мой код верен. Для меня это выглядит как-то странно, когда CPU синхронизируют кеши с основной памятью.
Поэтому мои вопросы:
[Изменить: старый код удален, см. ниже]
ОК, наконец, я смог создать более короткий (и самосогласованный) пример кода.
complex<double> mAllElements[tensorIdx][kappa1][kappa2][kappa3]. То есть У меня 3 ранга-3-тензора (tensorIdx). Каждый тензор представляет собой трехмерный массив, индексированный kappa1, kappa2 и kappa3.kappa1 - это тот, который получает parallized (и является самым внешним). Они расположены в DoComputation().main() я вызываю DoComputation() один раз, чтобы получить некоторые ссылочные значения, а затем я вызываю его несколько раз и сравниваю результаты. Они должны соответствовать точно, но иногда они этого не делают.К сожалению, код по-прежнему составляет около 190 строк. Я попытался упростить его (только 1 тензор ранга 1 и т.д.), Но потом я никогда не смог воспроизвести проблему. Я думаю, это происходит потому, что обращения к памяти не выровнены (цикл выше tensorIdx является самым внутренним) (я знаю, это далеко не оптимально).
Кроме того, некоторые задержки были необходимы в соответствующих местах, чтобы воспроизвести ошибку. В этом причина вызовов nops(). Без них код работает намного быстрее, но пока не показал проблему.
Обратите внимание, что я снова проверил критическую часть, CalcElementIdx(), и считаю ее правильной (каждый элемент обращается один раз). Я также запускал valgrind memcheck, helgrind и drd (с надлежащим перекомпилированным libgomp), и все трое не выдавали ошибок.
Каждые секунды до третьего запуска программы я получаю один или два несоответствия. Пример вывода:
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
Обе версии показывают описанную проблему при запуске на двух процессорах с общим количеством 40 потоков. Я не мог заметить ошибку на 1 CPU (и столько потоков, сколько мне нравится).
// 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 произошла ошибка. После обновления до 3.15.0-1 моя проблема исчезла, и код работает так, как должен.
Проблема возникла из-за ошибки в ядре ядра ядра ядра 3.11.10-7. Ошибка может быть связана с тем, как ядро обрабатывает недействительность кэша TLB, как отметил Христо Ильев. Я догадался, что ядро может быть проблемой, потому что я читал, что в Linux Kernel 3.15 для систем NUMA будут внесены некоторые улучшения, поэтому я решил, что ядро версия важна для систем NUMA.
Когда OP обновил ядро Linux своей системы NUMA до 3.15.0-1, проблема исчезла.