I am currently learning to use OpenMP with Fortran, and as practice I tried to generate two large arrays of random reals and multiply them together. The first thing I tried to do was to generate the arrays in parallel. Here I encountered a perplexing issue. At first my code looked like this:
program omp_test
use omp_lib
implicit none
!Variable definitions.
integer,parameter :: n=1000 !Matrix size.
real*8,dimension(n,n) :: A,B !Matrices.
integer :: row,col !Loop variables.
write(*,*) "Before OMP assignment."
!Parallel loop to assign a random number to every element in A and B.
!$OMP PARALLEL DO
do row = 1,n
do col=1,n
A(row,col)=rand()
B(row,col)=rand()
end do
end do
!$OMP END PARALLEL DO
write(*,*) "After OMP assignment."
end program omp_test
This code (compiled with gfortran and -fopenmp -O3 flags) results in the output
Segmentation fault (core dumped)
and nothing else, not even a backtrace. If I change the matrix size to 512x512 it still crashes the same way, but 511x511 works fine (changing OMP_STACKSIZE does not seem to help). It also works fine if I omit the -fopenmp flag during compilation. Note that not even the write statement executes. If I remove the assignment loop, or remove the OMP flags there are no problems, so the existance of the loop causes the program to crash before it reaches the loop. However, if I change the A and B assignment in the beginning to
real*8,dimension(:,:),allocatable :: A,B
and allocate them directly after the first write statement with
allocate(A(n,n))
allocate(B(n,n))
it runs fine for all matrix sizes. I suspect rand()s use of SAVE might be the culprit somehow, but I do not understand why it would work for allocatable arrays if that was the issue.
Does anyone know what the issue is and how to solve it?
