Summation error in openmp fortran

Question

I am trying to sum up of a variable with openmp with code given below.

    normr=0.0
 !$omp parallel default(private) shared(nelem,normr,cell_data,alphar,betar,k)
   !$omp do REDUCTION(+:normr)
     do ii=1,nelem
        nnodese=cell_data(ii)%num_vertex
        pe=cell_data(ii)%porder
        ndofe=cell_data(ii)%ndof             
        num_neighboure=cell_data(ii)%num_neighbour

        be=>cell_data(ii)%Force
        Ke=>cell_data(ii)%K
        Me=>cell_data(ii)%M
        pressuree=>cell_data(ii)%p
        Rese=>cell_data(ii)%Res
        neighbour_indexe=>cell_data(ii)%neighbour_index(:)

        Rese(:)=be(:)
        Rese(:)=Rese(:)-cmplx(-1.0,1.0*alphar/k)*matmul(Me(:,:),pressuree(:))
        Rese(:)=Rese(:)-cmplx(1.0,1.0*k*betar)*matmul(Ke(:,:),pressuree(:))

        do jj=1,num_neighboure  
           nbeindex=neighbour_indexe(jj)
           Knbe=>cell_data(ii)%neighbour(jj)%Knb
           pressurenb=>cell_data(nbeindex)%p
           ndofnb=cell_data(nbeindex)%ndof



           Rese(:)=Rese(:)-cmplx(1.0,1.0*k*betar)*matmul(Knbe(:,:),pressurenb(:))


           nullify(pressurenb)
           nullify(Knbe)
        end do   


        normr=normr+dot_product(Rese(:),Rese(:))

        nullify(pressuree)
        nullify(Ke)
        nullify(Me)
        nullify(Rese)
        nullify(neighbour_indexe)
        nullify(be)

     end do
  !$omp end do
!$omp end parallel

The result for summed variable, normr, is different for parallel and sequantial code. In one of the posts I have seen that inner loop variable should be defined inside the parallel construct(Why I don't know). I also changed the pointers to locall allocated variables but result did not changed. normr is a saved real variable.

Any suggestions and helps will be appreciated.

Best Regards,

Gokmen

Try to reduce the code to a [short, self contained, compilable example](http://sscce.org/). This will force you to analize the problem better and permit others to reproduce it. — Massimiliano, Jan 16 '14 at 06:20
If my answer to http://stackoverflow.com/questions/20993327/using-openmp-critical-and-ordered/20996615 doesn't explain your error then tell us how the parallel and sequential results differ. — High Performance Mark, Jan 16 '14 at 06:40
Given that the norm is an accumulation of non-negative numbers and the way reduction works, there shouldn't be a big difference in the final result with and without OpenMP. How much exactly do both results differ? — Hristo Iliev, Jan 16 '14 at 11:28

BHF · Answer 1 · 2014-01-16T06:52:30.557

normr can be different for the parallel and the sequential code, because the summation does not take place in the same order. Hence, the difference does not need to be an error and can be expected from the reduction operation.

Not being an error does not necessary mean not being a problem. One way around this would be to move the summation out of the parallel loop:

!$omp parallel default(private) shared(... keep_dot_product)
!$OMP do
do ii=1,nelem
   ! ...
   keep_dot_product(ii) = dot_product(Rese(:),Rese(:))
   ! ...
end do
!$omp end do
!$omp end parallel
normr = sum(keep_dot_product)

Summation error in openmp fortran

1 Answers1