В общем, вы должны быть осторожны и подумать о том, чтобы сделать переменную приватной, когда есть присваивание переменной. В противном случае существует риск гонок между потоками, что приведет к случайному некорректному поведению.
Это происходит в нескольких ситуациях в вашем коде.
compute_it = ... (уже приватная переменная)
агрегат + = ... (особый случай, требующий сокращения)
u = ... (должно быть приватно)
гистограмма [u] + = ... (опять проблема редукции)
Назначение элементам массива также может быть проблемой, но это зависит от индексов. Если индексы зависят от потока и различны в каждом потоке, то это, как правило, правильно, за исключением случаев ложного совместного использования. Это то, что происходит в большинстве ваших назначений для массивов. Например, для new[i*DIM*DIM+j*DIM+k]=... все потоки будут иметь разные значения i (благодаря параллели для), будут затронуты разные части массива, и нет конкретной проблемы с параллелизмом.
Для присвоения histogrammy[u] ситуация другая, так как u зависит от данных и может быть одинаковым в разных потоках. Им можно управлять с уменьшением количества новых версий omp, но в противном случае необходимо выполнить локальное накопление, и в конце потока глобальный массив обновится в должным образом защищенном регионе.
Вот переработанная версия вашего кода (непроверенная, поскольку вы не предоставили рабочий пример). Я также добавил некоторые комментарии и модификации, не связанные с распараллеливанием. Проверьте комментарий с тройкой ///
void work_it_par(long *old, long *new) {
int i, j, k;
int u, v, w;
long compute_it;
long aggregate=1.0; //// really?
//// I am really surprised that you use a long.
//// A double seems more appropriate
long weNeedTheFunc = we_need_the_func();
long gimmieTheFunc = gimmie_the_func();
int marker = DIM-1;
/// #pragma omp parallel for private(i, j, k, compute_it)
# pragma omp parallel for private(i, j, k, compute_it) reduction(+:aggregate)
/// introduced a reduction on aggregate
for (i=1; i<marker; i++) {
for (j=1; j<marker; j++) {
for (k=1; k<marker; k++) {
compute_it = old[i*DIM*DIM+j*DIM+k] * weNeedTheFunc;
/// aggregate+= compute_it / gimmieTheFunc; /// race on shared var aggregate
/// solved by the reduction
aggregate += compute_it ; /// Unrelated to parallel processing,
/// but do not do a division in your loop
/// divisions are *expensive* and
/// denominator is always the same
}
}
}
aggregate /= gimmieTheFunc ; /// now we do the division, but just once
printf("AGGR:%ld\n",aggregate);
//#pragma omp parallel for private(i, j, u, v)
for (i=1; i<marker; i++) {
#pragma omp parallel for private(k)
for (j=1; j<marker; j++) {
for (k=1; k<marker; k++){
new[i*DIM*DIM+j*DIM+k]=0;
for (u=-1; u<=1; u++) {
for (v=-1; v<=1; v++) {
for (w=-1; w<=1; w++) {
new[i*DIM*DIM+j*DIM+k]+=old[(i+u)*DIM*DIM+(j+v)*DIM+(k+w)];
}
}
}
new[i*DIM*DIM+j*DIM+k]/=27;
}
}
}
///#pragma omp parallel for private(i, j)
#pragma omp parallel private(i, j, u) /// parallel region
{
int private_histogrammy[10]; /// used to accumulate in the threads
for (int ii=0; ii<10; ii++) private_histogrammy[ii]=0;
# pragma omp for /// a parallel for loop in the parallel region
for (i=1; i<marker; i++) {
for (j=1; j<marker; j++) {
for (k=1; k<marker; k++) {
/// u=(new[i*DIM*DIM+j*DIM+k]/100);
u=(new[i*DIM*DIM+j*DIM+k]); /// to reduce number of divisions
/// if (u<=0) u=0;
/// if (u>=9) u=9;
/// histogrammy[u]++;
if (u<=0) private_histogrammy[0]++;
else if (u>=900) private_histogrammy[9]++;
else private_histogrammy[u/100]++;
}
}
}
/// all is done update the global histogrammy
# pragma omp critical
/// update the shared array
/// probably faster with a critical section that updates globally
/// the (small) array than using atomic on array elements
/// but alternatives may be tested
for (int uu=0; uu<10; uu++) histogrammy[uu] += private_histogrammy[uu];
} /// end of parallel region
}