Search code examples
fortranopenmpreduction

OpenMP reduction on user defined Fortran type containing allocatable array


I want to do an OpenMP reduction on a user defined Fortran type. I know OpenMP does not support Fortran types in reduction clauses but it is possible to define own reductions. This is done in the following example. This also works and does what it is expected to

 module types 
  !!! your type this can contain scalars and arrays
  type my_type
    Real*8,allocatable,dimension( : )  ::x
  end type

  !!! makes it possible to use the plus symbol for the reduction staement
  !!! replaces my_add by plus symbol
  interface operator(+)
     module procedure :: my_add
  end interface

 !$omp declare reduction (+ : my_type : omp_out = omp_out + omp_in) initializer (omp_priv = my_type ([0,0,0,0,0,0,0,0,0,0]))

 contains


  function my_add( a1 , a2 )
    type( my_type ),intent( in ) :: a1, a2
    type( my_type )              :: my_add
    my_add%x          =   a1%x + a2%x
    return
  end function my_add
 end module types 






program main
  use types
  use omp_lib
  type(my_type) :: my_var

  ! Initialize the reduction variable before entering the OpenMP region
  Allocate( my_var%x( 1:10 ) )  
  my_var%x = 0d0

  !$omp parallel reduction (+ : my_var) num_threads(4)
    my_var%x = omp_get_thread_num() + 6
    print*,omp_get_thread_num()
  !$omp end parallel

  print *, "sum of x is ", my_var%x
end program

My problem is now the allocatable array.

Because I hard coded the array initializer for the OpenMP reduction statement as initializer (omp_priv = my_type ([0,0,0,0,0,0,0,0,0,0])) I have to put 10 zeros there since the array is allocated with a length of 10. Is it possible to this with a variable name N (length of array)??


Solution

  • Inside the reduction initializer clause we have limited access to variables, making an array constructor of variable length difficult. However, we have available to us the Fortran version of the C++ approach.

    We can use the variable omp_orig to refer to the "storage of the original variable to be reduced":

    !$omp declare reduction (+ : my_type : omp_out = omp_out + omp_in) &
    !$omp&   initializer (omp_priv=omp_orig)
    

    The assignment statement here successfully allocates the array component of each private copy.