Search code examples
fortranopenmp

Fortran + OpenMP code with a subroutine stops abruptly


I have a piece of experimental code that works perfectly with serial compilation and execution. When I compile it with openmp option on ifort (on ubuntu), the compilation goes on fine but the execution stops abruptly. The code is as follows:

!!!!!!!! module
module array
implicit none
  real(kind=8),allocatable :: y(:)
end module array

module nonarray
implicit none
 real(kind=8):: aa
end module nonarray

use nonarray; use array
implicit none
integer(kind=8):: iter,i
integer(kind=8),parameter:: id=1
real(kind=8),allocatable:: yt(:)

allocate(y(id)); allocate(yt(id)); y=0.d0; yt=0.d0

aa=4.d0 !!A SYSTEM PARAMETER

 !$OMP PARALLEL PRIVATE(y,yt,iter,i)
 !$OMP DO

   loop1: do iter=1,20 !! THE INITIAL CONDITION LOOP
    call random_number(y)!! RANDOM INITIALIZATION OF THE VARIABLE

      loop2: do i=1,10000  !! ITERATION OF THE SYSTEM
       call evolve(yt)
       y=yt
      enddo loop2     !! END OF SYSTEM ITERATION

     write(1,*)aa,yt

   enddo loop1 !!INITIAL CONDITION ITERATION DONE

 !$OMP ENDDO
 !$OMP END PARALLEL
stop
end

recursive subroutine evolve(yevl)
use nonarray; use array
implicit none
integer(kind=8),parameter:: id=1
real(kind=8):: xf
real(kind=8),intent(out):: yevl(id)

  xf=aa*y(1)*(1.d0-y(1))
  yevl(1)=xf

end subroutine evolve

For compilation I use the following command: ifort -openmp -fpp test.f90.

test.f90 being the name of the program. Any suggestions or help is highly appreciated.


Solution

  • I am not an OMP expert, but I think if the subroutine evolve should see a different (private) y in each thread, you should pass it directly from within the parallelized code block to the subroutine instead of importing it from an external module:

    module common
      use iso_fortran_env
      implicit none
    
      integer, parameter :: dp = real64
      real(dp) :: aa
    
    contains
    
      subroutine evolve(y, yevl)
        implicit none
        real(dp), intent(in) :: y(:)
        real(dp), intent(out):: yevl(:)
    
        yevl(1) = aa * y(1) * (1.0_dp - y(1))
    
      end subroutine evolve
    
    end module common
    
    
    program test
      use common
      implicit none
    
      integer :: iter, i
      real(dp), allocatable :: yt(:), y(:)
    
      allocate(yt(1), y(1))
      y(:) = 0.0_dp
      yt(:) = 0.0_dp
      aa = 4.0_dp
    
      !$OMP PARALLEL DO PRIVATE(y,yt,iter,i)
      loop1: do iter = 1, 20
        call random_number(y)
        loop2: do i = 1, 10000
          call evolve(y, yt)
          y = yt
        end do loop2
        write(*,*) aa, yt
      end do loop1
      !$OMP END PARALLEL DO
    
    end program test
    

    Just an additional warning: the code above worked with various compilers (nagfor 5.3.1, gfortran 4.6.3, ifort 13.0.1), but not with ifort 12.1.6. So, although I can't see any obvious problems with it, I may have messed up something.