I have observed an unexpected (for me!) behavior of an openmp code which I am writing. The code structure is the following:
#pragma omp parallel for
for(int i=0;i<N;i++){
// lots of calculations that produce 3 integers i1,i2,i3 and 3 doubles d1,d2,d3
#pragma omp atomic
J1[i1] += d1;
#pragma omp atomic
J2[i2] += d2;
#pragma omp atomic
J3[i3] += d3;
}
I have compiled three different versions of this code:
1) with openmp (-fopenmp)
2) without openmp
3) with openmp, but without the 3 atomic operations (just as a test, since atomic operations are necessary)
When I run version 1) with environment variable OMP_NUM_THREADS=1, I observe a significant slowdown with respect to version 2); while version 3) runs as fast as version 2).
I would like to know the reason for this behavior (why do atomic operations slow the code down even if single-threded?!) and if it is possible to compile/rewrite the code in such a way that version 1) runs as fast as version 2).
I attach at the end of the question a working example which shows the aforementioned behavior. I compiled 1) with:
g++ -fopenmp -o toy_code toy_code.cpp -std=c++11 -O3
2) with:
g++ -o toy_code_NO_OMP toy_code.cpp -std=c++11 -O3
and 3) with:
g++ -fopenmp -o toy_code_NO_ATOMIC toy_code_NO_ATOMIC.cpp -std=c++11 -O3
The version of the compiler is gcc version 5.3.1 20160519 (Debian 5.3.1-20). The execution time of the 3 versions is:
1) 1 min 24 sec
2) 51 sec
3) 51 sec
Thanks in advance for any advice!
// toy_code.cpp
#include <stdio.h>
#include <iostream>
#include <stdlib.h>
#include <cmath>
#include <omp.h>
#define Np 1000000
#define N 1000
int main (){
double* Xp, *Yp, *J,*Jb;
Xp = new double[Np];
Yp = new double[Np];
J = new double [N*N];
Jb = new double [N*N];
for(int i=0;i<N*N;i++){
J[i]=0.0;
Jb[i]=0.0;
}
for(int i=0;i<Np;i++){
Xp[i] = rand()*1.0/RAND_MAX - 0.5;
Yp[i] = rand()*1.0/RAND_MAX - 0.5;
}
for(int n=0; n<2000; n++){
#pragma omp parallel for
for(int p=0;p<Np;p++){
double rx = (Xp[p]+0.5)*(N-1);
double ry = (Yp[p]+0.5)*(N-1);
int xindex = (int)floor(rx+0.5);
int yindex = (int)floor(ry+0.5);
int k;
k=xindex*N+yindex;
#pragma omp atomic
J[k]+=1;
#pragma omp atomic
Jb[k]+=1;
}
}
delete[] Xp;
delete[] Yp;
delete[] J;
delete[] Jb;
return 0;
}