1

I want to do some element-wise calculation on arrays in Fortran 90, while parallelize my code with openmp. I have now the following code :

program test
implicit none

integer,parameter :: n=50
integer :: i
integer(8) :: t1,t2,freq
real(8) :: seq(n),r(n,n,n,n)
real(8),dimension(n,n,n,n) :: x

call system_clock(COUNT_RATE=freq)

seq=[(i,i=1,n)]
x=spread(spread(spread(seq,2,n),3,n),4,n)

call system_clock(t1)

!$omp parallel workshare
! do some array calculation
r=atan(exp(-x))
!$omp end parallel workshare

call system_clock(t2)

print*, sum(r)
print '(f6.3)',(t2-t1)/real(freq)

end program test

I want now to replace the static arrays x and r with allocatable arrays, so I type :

real(8),dimension(:,:,:,:),allocatable :: x,r
allocate(x(n,n,n,n))
allocate(r(n,n,n,n))

but that the program run in serial without errors and the compiler doesn't take account of the line "!$omp parallel workshare".

What options should I use to parallelize in this case? I have tried with omp parallel do with loops but it is much slower.

I am compiling my code with gfortran 5.1.0 on windows :

gfortran -ffree-form test.f -o main.exe -O3 -fopenmp -fno-automatic
x1hgg1x
  • 115
  • 10
  • *"but it seems that !$omp parallel workshare doesn't work for allocatable array. "* Why? What does mean *"it seems"*? What does mean *"it doesn't work"*? Don't use these phrase which don't mean anything. Tell us what have you tried and what errors have you encountered. BTW, `real(8)` is ugly: http://stackoverflow.com/questions/838310/fortran-90-kind-parameter/856243#856243 http://stackoverflow.com/questions/3170239/fortran-integer4-vs-integer4-vs-integerkind-4/3170438#3170438 – Vladimir F Героям слава Mar 20 '16 at 15:30
  • It means that the program run in serial without errors and the compiler doesn't take account of the line "!$omp parallel workshare". – x1hgg1x Mar 20 '16 at 15:55
  • 1
    Related: https://stackoverflow.com/questions/17812003/parallelization-of-elementwise-matrix-multiplication/17832699#17832699 – Alexander Vogt Mar 20 '16 at 16:32

1 Answers1

4

I have come across this issue in gfortran before. The solution is to specify the array in the following form: 

!$omp parallel workshare
! do some array calculation
r(:,:,:,:) = atan(exp(-x))
!$omp end parallel workshare

Here is the reference.

Alexander Vogt
  • 17,879
  • 13
  • 52
  • 68
  • It works for me. In fact I have already tried it with static arrays but I saw no difference so I have removed it. So the difference is only seen for allocatable array. – x1hgg1x Mar 20 '16 at 16:31
  • @VladimirF On my machine with gfortran 5.3.1 the issue is still present. – Alexander Vogt Mar 20 '16 at 16:41
  • 3
    Perhaps the reason why the version in the answer works while the `r=atan(exp(-x))` version doesn't is that r is allocatable _and_ is being assigned to as a whole. This means that a reallocation could be required, which might trigger the deactivation of parallel workshare. On the other hand, using `r(:,:,:,:)` forces the expression result to be assigned to the _existing_ array r - possibly raising an error if bounds are incompatible depending on your build configuration. – Javier Martín Mar 20 '16 at 16:45
  • @JavierMartín That sounds plausible to me... Do you have a reference to support this? – Alexander Vogt Mar 20 '16 at 17:07
  • @AlexanderVogt no, that is why I stated it as a hypothesis with "perhaps". Sorry, but I'm not familiar enough with OpenMP to pinpoint things to the standard, or to undefined behaviour therein. However, it seems that "workshare" allows the compiler quite a bit of leeway in the division of work, so _perhaps_ the moment it sees something that might spell trouble (like a possible reallocation) it switches back to "single". – Javier Martín Mar 20 '16 at 17:15