Since OpenMP 4.5, arrays can be reduced using pragmas. But, I faced a problem when I tried to use it on arrays of derived type because !$OMP DECLARE REDUCTION is not found for this type. I understand the problem. After many tries, I finally found a program that seems to work. This is the code :
The module (module.f90):
module my_mod
implicit none
type :: my_type_t
integer :: x
integer :: y
end type my_type_t
interface operator(+)
module procedure add_my_type_t_vec
end interface
!$OMP DECLARE REDUCTION (+ : my_type_t : omp_out = omp_out + omp_in) INITIALIZER (omp_priv = omp_orig)
contains
function add_my_type_t_vec(a,b) result(res)
type(my_type_t), dimension(:),intent(in) :: a, b
type(my_type_t), dimension(size(a)) :: res
res(:)%x = a(:)%x + b(:)%x
res(:)%y = a(:)%y + b(:)%y
end function add_my_type_t_vec
end module my_mod
And the program (pgm.f90) that use it :
program parallel
use omp_lib
use my_mod
implicit none
type(my_type_t), dimension(:), allocatable :: private_var, shared_var
integer :: i
allocate(private_var(5), shared_var(5))
shared_var%x=0
shared_var%y=0
!$OMP PARALLEL NUM_THREADS(4) PRIVATE(private_var,i) REDUCTION (+ : shared_var)
!$OMP DO
do i=1,1000
private_var%x = 1
private_var%y = 2
shared_var = shared_var + private_var
end do
!$OMP END DO
!$OMP END PARALLEL
print *, "Total : ", shared_var
end program parallel
Compilation with : gfortran -fopenmp module.f90 pgm.f90
To find this solution, I used the following answer https://stackoverflow.com/a/61148205/7462275.
It works but I am not sure that I will not have undefined behaviour in certain circumstances. And I do not fully understand the syntax and the mechanism of the DECLARE REDUCTION and INITIALIZER pragmas. Indeed, I did a copy/paste of the answer mentioned above and there is not explanation in it. Thanks for advice and explanations.
The issue with INITIALIZER (omp_priv = omp_orig) is that it only works, if the reduction variable is 0 before entering the region with the reduction. The private instance of the reduction variable should be initialized to the neutral element of the reduction operation. For + the neutral element is 0. For min it would be "positive infinity" or whatever is the closest representation for the data type.
In this specific case, we need an array of 0 elements that has the right size which we can achieve by adding an initializer function that copies the size of the array:
module my_mod
implicit none
type :: my_type_t
integer :: x
integer :: y
end type my_type_t
interface operator(+)
module procedure add_my_type_t_vec
end interface
!$OMP DECLARE REDUCTION (+ : my_type_t : omp_out = omp_out + omp_in) INITIALIZER (omp_priv = init_my_type_vec(omp_orig))
contains
function add_my_type_t_vec(a,b) result(res)
type(my_type_t), dimension(:),intent(in) :: a, b
type(my_type_t), dimension(size(a)) :: res
res%x = a%x + b%x
res%y = a%y + b%y
end function add_my_type_t_vec
function init_my_type_vec(vec_in) result(vec_out)
type(my_type_t), dimension(:),intent(in) :: vec_in
type(my_type_t), dimension(size(vec_in)) :: vec_out
vec_out%x=0
vec_out%y=0
end function init_my_type_vec
end module my_mod
Now, modifying the application code to enter the reduction with a different value shows the difference:
program parallel
use omp_lib
use my_mod
implicit none
type(my_type_t), dimension(:), allocatable :: private_var, shared_var
integer :: i
allocate(private_var(5), shared_var(5))
shared_var = my_type_t(5,10)
!$OMP PARALLEL NUM_THREADS(4) PRIVATE(private_var,i) REDUCTION (+ : shared_var)
!$OMP DO
do i=1,1000
private_var%x = 1
private_var%y = 2
shared_var = shared_var + private_var
end do
!$OMP END DO
!$OMP END PARALLEL
print *, "Total : ", shared_var
end program parallel
Prints with the original initializer:
Total : 1025 2050 1025 2050 1025 2050 1025 2050 1025 2050
and with the modified initializer:
Total : 1005 2010 1005 2010 1005 2010 1005 2010 1005 2010
For further reading on how to use reductions I suggest to look at the OpenMP Examples document (openmp.org/wp-content/uploads/openmp-examples-6.0.pdf) to find more udr C and Fortran examples.