3

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
Franky
  • 73
  • 4

1 Answers1

4

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
Alexander Vogt
  • 17,879
  • 13
  • 52
  • 68
  • So `reduction` works as some kind of synchronization to ensure the value of summation at end of the loop is correct? I just want to make sure, it is also incorrect to put `index` in private, isn't it? because the value of summation should be in sequence not distributed over threads? – Franky May 24 '15 at 15:45
  • Reduction works exactly as @AlexanderVogt describes it, private copies of the variable are reduced at the end of the construct. Private also does not help, you would have no synchronization and you would get some values of `index` duplicated. – Vladimir F Героям слава May 24 '15 at 15:46
  • 1
    @FrankyDjutanta Why don't you just compute `index` from `ix` and `iy`? – Vladimir F Героям слава May 24 '15 at 15:48
  • Oh ya I just took the simplest one as the code was initially developed for serial, didn't think of it before. Thank you very much both of you. Much appreciated. – Franky May 24 '15 at 16:12
  • @AlexanderVogt, you forgot to set number(1) in `private` since it is reused inside the loop. – Franky May 24 '15 at 16:24
  • @FrankyDjutanta making `number` private is not required. – Alexander Vogt May 24 '15 at 16:48
  • I agree @AlexanderVogt, I was wrong in checking the code before where I got incorrect result without putting `number(1)` in private and got the correct one with putting that in private, but it turned out that I commented `!$OMP` on the second try. I have tried the code correctly now, but it seems to delete the value of `xoutput` and `youtput` as I tried to do `write(*,*)` it showed nothing while it is okay when compiling in serial. – Franky May 24 '15 at 16:51