error in OpenMP program execution

Question

I am trying to sum all the elements of an array which is initialized in the same code. As every element is independent from each other, I tried to perform the sum in parallel. My code is shown below:

int main(int argc, char** argv) 
{
  cout.precision(20);
  double sumre=0.,Mre[11];
  int n=11;
  for(int i=0; i<n; i++)
   Mre[i]=2.*exp(-10*M_PI*i/(1.*n));

  #pragma omp parallel for reduction(+:sumre)
  for(int i=0; i<n; i++)
  {
   sumre+=Mre[i];
  }
  cout<<sumre<<"\n";
}

which I compile and run with:

g++ -O3 -o sum sumparallel.cpp -fopenmp
./sum

respectively. My problem is that the output differs every time I run it. Sometimes it gives

2.1220129388411006488

or

2.1220129388411002047 Does anyone have an idea what is happening here?

Answer 1

Some of these comments hint at the problem here, but there could be two distinct issues

Double precision numbers do not have 20 decimal precision

If you want to print the maximum precision of sumre, use something like this

#include <float.h>
int maint(int argc, char* argv[])
{
  ...

  printf("%.*g", DBL_DECIMAL_DIG, number);

  return 0;
}

Floating point arithmetic is non-commutative

The effect of this property is roundoff error. In fact, the function that you have defined, a gaussian, is especially prone to roundoff for summation. Considering that workload distribution of OpenMP parallel for is undefined, you may receive different answers when you run it. To get around this you may use the kahan summation algorithm. An implementation with OpenMP would look something like this:

  ...
  double sum = 0.0, c = 0.0; 

  #pragma omp parallel for reduction(+:sum, +:c) 
  for(i = 0; i < n; i++) 
  { 
    double y = Mre[i] - c; 
    double t = sum + y; 
    c = (t - sum) - y; 
    sum = t; 
  } 

  sum = sum - c;  
  ...