1

I aim to compute a simple N-body program on C++ and I am using OpenMP to speed things up with the computations. At some point, I have nested loops that look like that:

int N;
double* S          = new double[N];
double* Weight     = new double[N];
double* Coordinate = new double[N];
 ...

#pragma omp parallel for
for (int i = 0; i < N; ++i)
{
   for (int j = 0; j < i; ++j)
   {
       double K = Coordinate[i] - Coordinate[j];
       S[i]    += K*Weight[j];
       S[j]    -= K*Weight[i];
   }
}

The issue here is that I do not obtain exactly the same result when removing the #pragma ... I am guessing it has to do with the fact that the second loop is dependent on the integer i, but I don't see how to get past that issue

Jeffrey
  • 11,063
  • 1
  • 21
  • 42
Johan Valentin
  • 13
  • 2
  • 1
    The dependent inner loop would only be a problem in terms of correct results if you were using the `collapse` clause. It may however affect performance, if OpenMP uses `static` scheduling per default, as different threads might have very different workloads. Therefore you should make sure to use `dynamic` scheduling, e.g. with a `schedule(dynamic)` clause. It's explained in depth under [this question](https://stackoverflow.com/q/10850155/10107454). – paleonix Jul 30 '21 at 07:19

1 Answers1

2

The problem is that there is a data race during updating S[i] and S[j]. Different threads may read from/write to the same element of the array at the same time, therefore it should be an atomic operation (you have to add #pragma omp atomic) to avoid data race and to ensure memory consistency:

for (int j = 0; j < i; ++j)
{
    double K = Coordinate[i] - Coordinate[j];
    #pragma omp atomic
    S[i]    += K*Weight[j];    
    #pragma omp atomic
    S[j]    -= K*Weight[i];
}
Laci
  • 2,738
  • 1
  • 13
  • 22
  • Thank you a lot it really helped decrease the error. Their is still a 1e-15 relative error between the unparalleled and the paralleled loops but I think it has more to do with floating addition error than the parallelisation. – Johan Valentin Jul 30 '21 at 13:39
  • You are welcome. This small error is not surprising because the order of additions are different, so the rounding errors are also different. – Laci Jul 30 '21 at 15:48
  • If you are interested in the details of rounding error, please read this: https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html – Laci Jul 30 '21 at 16:06