OpenMP over Summation

Question

I have been trying to apply OpenMP on a simple summation operation inside two nested loop, but it produced incorrect result so far. I have been looking around in here and here, also in here. All suggest to use reduction clause, but it does not work for my case by producing very large number which leads to segmentation fault.

I also tried this way posted in here and my own question here which has been solved. Both do not use reduction and simply just set summation variable as shared, but it also produces incorrect result. Is there something that I am missing? When to use reduction and not using that while facing summation operation?

Codes using reduction clause

index = 0
!$OMP PARALLEL DO PRIVATE(iy,ix) REDUCTION(:+index)
do iy = 1, number(2)
    do ix = 1, number(1)
        index = index + 1
        xoutput(index)=xinput(ix)
        youtput(index)=yinput(iy)
    end do
end do
!$OMP END PARALLEL DO

Code without using reduction clause

index = 0
!$OMP PARALLEL DO PRIVATE(iy,ix) SHARED(index)
do iy = 1, number(2)
    do ix = 1, number(1)
        index = index + 1
        xoutput(index)=xinput(ix)
        youtput(index)=yinput(iy)
    end do
end do
!$OMP END PARALLEL DO

Answer 1

I think you have a mis-conception of what the reduction clause does...

REDUCTION(+:index)

means, that you will have the correct sum index in the end. In each step of the iteration, each tread will have a different version with different values! So the reduction is not suitable to manage array indices during the parallel section.

Let me try to illustrate this...

The following loop

!$OMP PARALLEL DO PRIVATE(iy) REDUCTION(+:index)
do iy = 1, number(2)
  index = index + 1
end do
!$OMP END PARALLEL DO

is (more or less) equivalent to

!$OMP PARALLEL PRIVATE(iy, privIndex) SHARED(index)
!$OMP DO
do iy = 1, number(2)
  privIndex = privIndex + 1
end do
!$OMP END DO

!$OMP CRITICAL
index = index + privIndex
!$OMP END CRITICAL
!$OMP END PARALLEL

You can see that during the loop all threads work on different variables privIndex which are private to that thread, and calculate local (partial) sums. In the end, the total sum is taken, using a critical section to avoid race conditions.

This might not be what the compiler does, but it gives you an idea how a reduction works: at no point within the first loop does privIndex correspond to the correct index you would expect in the serial version.

---

As Vladimir suggests in his comment, you can calculate the index directly as you are only incrementing it in the inner loop:

!$OMP PARALLEL DO PRIVATE(iy,ix, index)
do iy = 1, number(2)
    do ix = 1, number(1)
        index = (iy-1)*number(1) + ix
        xoutput(index)=xinput(ix)
        youtput(index)=yinput(iy)
    end do
end do
!$OMP END PARALLEL DO