Красно-черный Гаусс Зейдель и OpenMP
Я пытался доказать свою точку зрения с OpenMP по сравнению с MPICH, и я подготовил следующий пример, чтобы продемонстрировать, насколько легко было добиться высокой производительности в OpenMP.
Итерация Гаусса-Зейделя разбита на два отдельных цикла, так что в каждом цикле каждая операция может выполняться в любом порядке, и между каждой задачей не должно быть никакой зависимости. Таким образом, теоретически каждый процессор никогда не должен ждать, пока другой процесс выполнит какую-либо синхронизацию.
Проблема, с которой я сталкиваюсь, заключается в том, что я, независимо от размера проблемы, нахожу, что есть только слабое ускорение двух процессоров, а с более чем двумя процессорами это может быть даже медленнее.
Во многих других линейных параллельных программах я могу получить очень хорошее масштабирование, но это сложно.
Я боюсь, что не смогу «объяснить» компилятору, что операция, которую я выполняю над массивом, является поточно-ориентированной, так что она не может быть действительно эффективной.
Смотрите пример ниже.
Кто-нибудь знает, как сделать это более эффективным с OpenMP?
void redBlackSmooth(std::vector<double> const & b,
std::vector<double> & x,
double h)
{
// Setup relevant constants.
double const invh2 = 1.0/(h*h);
double const h2 = (h*h);
int const N = static_cast<int>(x.size());
double sigma = 0;
// Setup some boundary conditions.
x[0] = 0.0;
x[N-1] = 0.0;
// Red sweep.
#pragma omp parallel for shared(b, x) private(sigma)
for (int i = 1; i < N-1; i+=2)
{
sigma = -invh2*(x[i-1] + x[i+1]);
x[i] = (h2/2.0)*(b[i] - sigma);
}
// Black sweep.
#pragma omp parallel for shared(b, x) private(sigma)
for (int i = 2; i < N-1; i+=2)
{
sigma = -invh2*(x[i-1] + x[i+1]);
x[i] = (h2/2.0)*(b[i] - sigma);
}
}
Дополнение:
Теперь я также попытался с использованием необработанного указателя, и он ведет себя так же, как и контейнер STL, поэтому можно исключить, что это какое-то псевдокритическое поведение, поступающее из STL.
Решение
Прежде всего, убедитесь, что x вектор выровнен по границам кэша. Я провел некоторый тест, и у меня получилось что-то вроде 100% улучшения с вашим кодом на моей машине (core duo), если я форсирую выравнивание памяти:
double * x;
const size_t CACHE_LINE_SIZE = 256;
posix_memalign( reinterpret_cast<void**>(&x), CACHE_LINE_SIZE, sizeof(double) * N);
Во-вторых, вы можете попытаться назначить больше вычислений каждому потоку (таким образом, вы можете держать разделенные строки кэша), но я подозреваю, что openmp уже делает что-то подобное под капотом, так что это может быть бесполезно с большим N.
В моем случае эта реализация намного быстрее, когда x не выровнен по кешу.
const int workGroupSize = CACHE_LINE_SIZE / sizeof(double);
assert(N % workGroupSize == 0); //Need to tweak the code a bit to let it work with any N
const int workgroups = N / workGroupSize;
int j, base , k, i;
#pragma omp parallel for shared(b, x) private(sigma, j, base, k, i)
for ( j = 0; j < workgroups; j++ ) {
base = j * workGroupSize;
for (int k = 0; k < workGroupSize; k+=2)
{
i = base + k + (redSweep ? 1 : 0);
if ( i == 0 || i+1 == N) continue;
sigma = -invh2* ( x[i-1] + x[i+1] );
x[i] = ( h2/2.0 ) * ( b[i] - sigma );
}
}
В заключение, у вас определенно есть проблема с кешированием, но, учитывая работу openmp (к сожалению, я не знаком с ней), этого должно быть достаточно для работы с правильно распределенными буферами.
Другие решения
Я думаю, что основная проблема связана с типом структуры массива, которую вы используете. Попробуем сравнить результаты с векторами и массивами. (Массивы = c-массивы с использованием оператора new).
Вектор и размеры массива N = 10000000. Я заставляю функцию сглаживания повторяться, чтобы поддерживать время выполнения> 0,1 сек.
Vector Time: 0.121007 Repeat: 1 MLUPS: 82.6399
Array Time: 0.164009 Repeat: 2 MLUPS: 121.945
MLUPS = ((N-2)*repeat/runtime)/1000000 (Million Lattice Points Update per second)
MFLOPS вводят в заблуждение, когда дело доходит до расчета сетки. Несколько изменений в базовом уравнении могут привести к высокой производительности при одном и том же времени выполнения.
Модифицированный код:
double my_redBlackSmooth(double *b, double* x, double h, int N)
{
// Setup relevant constants.
double const invh2 = 1.0/(h*h);
double const h2 = (h*h);
double sigma = 0;
// Setup some boundary conditions.
x[0] = 0.0;
x[N-1] = 0.0;
double runtime(0.0), wcs, wce;
int repeat = 1;
timing(&wcs);
for(; runtime < 0.1; repeat*=2)
{
for(int r = 0; r < repeat; ++r)
{
// Red sweep.
#pragma omp parallel for shared(b, x) private(sigma)
for (int i = 1; i < N-1; i+=2)
{
sigma = -invh2*(x[i-1] + x[i+1]);
x[i] = (h2*0.5)*(b[i] - sigma);
}
// Black sweep.
#pragma omp parallel for shared(b, x) private(sigma)
for (int i = 2; i < N-1; i+=2)
{
sigma = -invh2*(x[i-1] + x[i+1]);
x[i] = (h2*0.5)*(b[i] - sigma);
}
// cout << "In Array: " << r << endl;
}
if(x[0] != 0) dummy(x[0]);
timing(&wce);
runtime = (wce-wcs);
}
// cout << "Before division: " << repeat << endl;
repeat /= 2;
cout << "Array Time:\t" << runtime << "\t" << "Repeat:\t" << repeat
<< "\tMLUPS:\t" << ((N-2)*repeat/runtime)/1000000.0 << endl;
return runtime;
}
Я не изменил ничего в коде, кроме как тип массива. Для лучшего доступа к кешу и блокирования вы должны посмотреть на выравнивание данных (_mm_malloc).