OpenMP issue can someone help?

[Kay Diederichs]

Am 14.12.2010 12:25, schrieb Archambault Fabien:
> Hi all,
>
> I really hope my question is not silly but I am completely lost with
> that issue. Concerning the background: this issue is part of a much
> bigger software but I was able to isolate the part which is in fault and
> I made you an example of the code for test purposes.
>
> I am using gfortran --version: GNU Fortran (GCC) 4.5.1 20100924 (Red Hat
> 4.5.1-4)
> Which is the version available in Fedora 14. I did not compile on my own
> gcc of gfortran.
>
> When you try to run the program with that compilation (no OpenMP flag):
> gfortran -c test.F -Wall && gfortran -o test test.o -Wall && ./test
>
> All works (no segmentation fault). I can assure you that also the
> results given are correct!
>
> Then if you try with the OMP part included:
> $ gfortran -c test.F -fopenmp -DOMP -Wall && gfortran -o test test.o
> -fopenmp -Wall && ./test
> Segmentation fault (core dumped)
>
> If you run with the check option:
> $ gfortran -c test.F -fopenmp -DOMP -Wall -fcheck=all && gfortran -o
> test test.o -fopenmp -Wall -fcheck=all && ./test
> At line 30 of file test.F
> Fortran runtime error: Index '0' of dimension 1 of array 'y' below lower
> bound of 1
>
> I am assuming there is an issue with the OpenMP stuff but, even if the
> code is not really optimized, I think the code is correct. Am I right?
> I have also the code working with ifort without issue (I know it is not
> a good comparison seeing how lax is that compiler).
>
> If someone as an idea about that, please help me,
> Sincerely,
> Fabien
>

Fabien,

you have misunderstood the meaning of the reduction clause. The OpenMP 
3.0 spec says (2.9.3.6): "The reduction clause specifies an operator and 
one or more list items. For each list item, a private copy is created in 
each implicit task, and is initialized appropriately for the operator. 
After the end of the region, the original list item is updated with the 
values of the private copies using the specified operator." (more 
detailed description follows further down)

This means that each thread will have its own copy of j, initialized to 0.
Since you use j as an array index, this leads to the segfault.
Furthermore, I bet the program does not what you want it to do - you 
want indices to be added. This works for the serial version, but not for 
the parallel one.
You could use something like the following (untested!):

!$omp parallel if ( n > nmax_omp ) private (i,j)
!$omp& shared ( y, a, inc, n )
!$omp do
       do i=1,n
         j = 1 + (i-1)*inc
         y(j) = a
       enddo
!$omp end do
!$omp end parallel

The intended use of the reduction clause is shown in A.35 of the spec.

HTH,

Kay