This problem, as well as the one in your other question, comes from the fact that automatic heap arrays are disabled when OpenMP is enabled. This means that without -fopenmp
, big arrays are automatically placed in the static storage (known as the .bss
segment) while small arrays are allocated on the stack. When you switch OpenMP support on, no automatic static allocation is used and your dense
arrays gets allocated on the stack of the routine. The default stack limits on OS X are very restrictive, hence the segmentation fault.
You have several options here. The first option is to make dense
have static allocation by giving it the SAVE
attribute. The other option is to explicitly allocate it on the heap by making it ALLOCATABLE
and then using the ALLOCATE
statement, e.g.:
REAL, DIMENSION(:,:,:), ALLOCATABLE :: dense
ALLOCATE(dense(256,256,256))
! Computations, computations, computations
DEALLOCATE(dense)
Newer Fortran versions support automatic deallocation of arrays without the SAVE
attribute when they go out of scope.
Note that your OpenMP directive is just fine and no additional data sharing clauses are necessary. You do not need to declare i
in a PRIVATE
clause since loop counters have predetermined private data-sharing class. You do not need to put the other variables in SHARED
clause as they are implicitly shared. Yet the operations that you do on dense
should either be synchronised with ATOMIC UPDATE
(or simply ATOMIC
on older OpenMP implementations) or you should use REDUCTION(+:dense)
. Atomic updates are translated to locked additions and should not incur much of a slowdown, compared to the huge slowdown from having conditionals inside the loop:
INTEGER :: xi, yi, zi
!$OMP PARALLEL DO PRIVATE(xi,yi,zi)
...
if (x1(i).gt.0.and.y1(i).gt.0.and.z1(i).gt.0) then
xi = int(x1(i))
yi = int(y1(i))
zi = int(z1(i))
!$OMP ATOMIC UPDATE
dense(xi,yi,zi) = dense(xi,yi,zi) &
+ dx1(i) * dy1(i) * dz1(i) * mpart
end if
...
Replicate the code with the proper changes for the other cases. If your compiler complains about the UPDATE
clause in the ATOMIC
construct, simply delete it.
REDUCTION(+:dense)
would create one copy of dense
in each thread, which would consume a lot of memory and the reduction applied in the end would grow slower and slower with the size of dense
. For small arrays it would work better than atomic updates.