Subroutine inside OpenMP PARALLEL DO - Program Crash

Question

I have read both Calling an internal subroutine inside OpenMP region and Global Variables in Fortran OpenMP. My understanding (from here) is that:

Variables in the argument list inherit their data scope attribute from the calling routine.
COMMON blocks or module variables in Fortran are shared, unless declared THREADPRIVATE.
SAVE variables in Fortran are shared.
All other local variables are private.

The following is a simplified version of my code:

!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(j,dummy1,dummy2,dummy3,dummy4)
DO j=1,ntotal
  dummy1 = 0.0d0
  dummy2 = foo(j)
  CALL kernel(dummy1,dummy1,dummy2,dummy3,dummy4)  
  Variable(j) = dummy3 + dummy4
END DO 
!$OMP END PARALLEL DO

The subroutine kernel then takes IN dummy1, and dummy2 and puts OUT dummy3 and dummy4. I compile with:

 -fopenmp -fno-automatic -fcheck=all

I get:

Fortran runtime error: Recursive call to nonrecursive procedure 'kernel'

Which from here I understand is expected. When I compile without the -fcheck sometimes the code passes the subroutine call without incident but the majority of the time it will crash with no error. I am guessing this is because my subroutine is not thread safe. All the arguments passed to the subroutine should be private and individual for each thread. The trimmed down subroutine is as follows:

SUBROUTINE kernel(r,dx,hsml,w,dwdx)   

  USE Initial_Parameters
  IMPLICIT NONE 

  ! DATA DICTIONARY: DECLARE CALLING PARAMETER TYPES AND DEFINITIONS
  REAL(KIND=dp), INTENT(IN)                 ::  r           
  REAL(KIND=dp), DIMENSION(dim), INTENT(IN) ::  dx          
  REAL(KIND=dp), INTENT(IN)                 ::  hsml        
  REAL(KIND=dp), INTENT(OUT)                ::  w           
  REAL(KIND=dp), DIMENSION(dim), INTENT(OUT)::  dwdx        
  ! DATA DICTIONARY: DECLARE LOCAL VARIABLE TYPES AND DEFINITIONS
  INTEGER                                   ::  i, d  
  REAL(KIND=dp)                             ::  q, dw
  REAL(KIND=dp)                             ::  factor      


  ! Kernel functions are funcitons of q, the distance between particles
  ! divided by the smoothing length  
  q = r/hsml 
  ! Preset the kernel to zero
  w = 0.e0
  ! Preset the derivative of the kernel to zero
  DO d=1,dim         
    dwdx(d) = 0.e0
  END DO   

  IF (skf == 1) THEN     

    ! If the problem is one dimensional then,
    IF (dim == 1) THEN
      ! The coefficient, alpha = factor is given by:
      factor = 1.e0/hsml
    ! If the problem is two dimensional then,
    ELSE IF (dim == 2) THEN
      ! The coefficient, alpha = factor is given by:
      factor = 15.e0/(7.e0*pi*hsml*hsml)
    ! If the problem is two dimensional then,
    ELSE IF (dim == 3) THEN
      ! The coefficient, alpha = factor is given by:
      factor = 3.e0/(2.e0*pi*hsml*hsml*hsml)
    ! If the dimension value is not 1, 2 or 3 then there is a problem.
    ELSE
       WRITE(*,*)' >>> Error <<< : Wrong dimension: Dim =',dim
       STOP
    END IF

    ! Smoothing function for 1st range of q.                                         
    IF (q >= 0 .AND. q <= 1.e0) THEN
      ! The smoothing function is given by:
      w = factor * (2./3. - q*q + q*q*q / 2.)
      ! For each dimension work out the gradient of the smoothing function
      DO d = 1, dim
        dwdx(d) = factor * (-2.+3./2.*q)/hsml**2 * dx(d)       
      END DO   

    ! Smoothing function for 2nd range of q.  
    ELSE IF (q > 1.e0 .AND. q <= 2) THEN  
      ! Smoothing function is equal to:        
      w = factor * 1.e0/6.e0 * (2.-q)**3 
      ! Gadient of the smoothing function in each dimension.
      DO d = 1, dim
        dwdx(d) =-factor * 1.e0/6.e0 * 3.*(2.-q)**2/hsml * (dx(d)/r)        
      END DO   

    ! Smoothing function and gradient for all other values of q is zero.
    ELSE
      ! Smoothing function is equal to: 
      w=0.
      ! Gadient of the smoothing function in each dimension.
      DO d= 1, dim
        dwdx(d) = 0.
      END DO             
    END IF     

END SUBROUTINE kernel

The local variables should be private and all the argumets which have been passed are private. The module parameters are shared but this is fine. Could you please explain why this crashes?

d_1999 d_1999 · Accepted Answer · 2016-06-01T15:31:18

With -fno-automatic the local variables in kernel will be implicitly SAVEd. The description here notes

Local variables with the SAVE attribute declared in procedures called from a parallel region are implicitly SHARED.

As such, kernel is indeed not thread-safe (as I understand it).

Note also in your example you pass dummy1 as the first and second arguments to kernel, but your definition of this routine specifies the first argument (r) in a scalar whilst the second (dx) is an array of length dim. I'm not sure if this is just a product of your minimal example or of your real code, but this could cause issues. Are you declaring kernel inside a module and then using the module? This will generate interfaces which should help catch this sort of thing.

Subroutine inside OpenMP PARALLEL DO - Program Crash

1 Answers