Parallelizing pairwise gravitational force calculation with OpenMP in Fortran 90

Question

I am trying to parallelize the calculation of gravitational forces in my program using OpenMP. The calculation of the distances (R and R2) is no problem but the forces/accelerations (A) come out wrong. I know that it has something to do with race conditions in the summation. I have experimented a bit with atomic and critical constructs but could not find a solution. Also, I'm not sure which variables should be private and why.

Does someone with more experience in using OpenMP have a suggestion on how to correct this in the following code example?

A = 0.0
!$omp parallel do
do i = 1, Nobj
    do j = i + 1, Nobj
        R2(i,j) = (X(j,1) - X(i,1))**2 &
                + (X(j,2) - X(i,2))**2 &
                + (X(j,3) - X(i,3))**2
        R(i,j) = sqrt(R2(i,j))
        do k = 1, 3
            A(i,k) = A(i,k) + ((mass_2_acc(i,j) / R2(i,j)) * ((X(j,k) - X(i,k)) / R(i,j)))
            A(j,k) = A(j,k) + ((mass_2_acc(i,j) / R2(i,j)) * ((X(i,k) - X(j,k)) / R(i,j)))
        enddo
    enddo
    A(i,:) = A(i,:) * G / mass_acc(i)
enddo
!$omp end parallel do

Answer 1

You are modifying A(j,k) - neither j nor k are "local" to thread as the thread-parallel index is i. What I mean is neither of those index ranges are restricted to a particular thread, all threads will update all A(j,k) hence the race condition.

Things you can do - split up R and A calculations or not use symmetry to update A.

Also, Fortran is column major and you are traversing outer index first which is bad for performance.