I am trying to adapt a Fortran code (Gfortran) to make use of OpenMP. It is a particle based code where the index of arrays can correspond to particles or pairs. The code uses a derived type to store a number of matrices for each particle. It is very common to come across loops which require the use of a matrix stored in this derived type. This matrix may be accessed by multiple threads. The loop also requires a reduction over an element in this derived type. I currently have to write a temporary array in order to do this reduction and then I set the element of the derived type equal to this temporary reduction array. If not using OpenMP no temporary array is needed.
Question: Is it possible to do a reduction over an element of a derived type? I don't think I can do a reduction over the entire derived type as I need to access some of the elements in the derived type to do work, which means it needs to be SHARED. (From reading the specification I understand that when using REDUCTION a private copy of each list item is created.)
Complete minimal working example below. It could be more minimal but I feared that removing more components might over simplify the problem.
PROGRAM TEST_OPEN_MP
USE, INTRINSIC :: iso_fortran_env
USE omp_lib
IMPLICIT NONE
INTEGER, PARAMETER :: dp = REAL64
INTEGER, PARAMETER :: ndim=3
INTEGER, PARAMETER :: no_partic=100000
INTEGER, PARAMETER :: len_array=1000000
INTEGER :: k, i, ii, j, jj
INTEGER, DIMENSION(1:len_array) :: pair_i, pair_j
REAL(KIND=dp), DIMENSION(1:len_array) :: pair_i_r, pair_j_r
REAL(KIND=dp), DIMENSION(1:no_partic) :: V_0
REAL(KIND=dp), DIMENSION(1:ndim,1:no_partic) :: disp, foovec
REAL(KIND=dp), DIMENSION(1:ndim,1:len_array) :: dvx
REAL(KIND=dp), DIMENSION(1:2*ndim,1:len_array):: vec
REAL(KIND=dp), DIMENSION(1:ndim) :: disp_ij,temp_vec1,temp_vec2
REAL(KIND=dp), DIMENSION(1:ndim,1:ndim) :: temp_ten1,temp_ten2
REAL(KIND=dp), DIMENSION(1:no_partic,1:ndim,1:ndim):: reduc_ten1
REAL(KIND=dp) :: sum_check1,sum_check2,cstart,cend
TYPE :: matrix_holder !<-- The derived type
REAL(KIND=dp), DIMENSION(1:ndim,1:ndim) :: mat1 !<-- The first element
REAL(KIND=dp), DIMENSION(1:ndim,1:ndim) :: mat2 !<-- The second element, etc.
END TYPE matrix_holder
TYPE(matrix_holder), DIMENSION(1:no_partic) :: matrix
! Setting "random" values to the arrays
DO k = 1, no_partic
CALL random_number(matrix(k)%mat1(1:ndim,1:ndim))
CALL random_number(matrix(k)%mat2(1:ndim,1:ndim))
END DO
CALL random_number(pair_i_r)
CALL random_number(pair_j_r)
CALL random_number(disp)
CALL random_number(vec)
CALL random_number(dvx)
CALL random_number(V_0)
disp = disp*10.d0
vec = vec*100.d0
dvx = dvx*200.d0
V_0 = V_0*10d0
pair_i = FLOOR(no_partic*pair_i_r)+1
pair_j = FLOOR(no_partic*pair_j_r)+1
! Doing the work
cstart = omp_get_wtime()
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP& PRIVATE(i,j,k,disp_ij,temp_ten1,temp_ten2,temp_vec1,temp_vec2,ii,jj), &
!$OMP& REDUCTION(+:foovec,reduc_ten1), SCHEDULE(static)
DO k= 1, len_array
i = pair_i(k)
j = pair_j(k)
disp_ij(1:ndim) = disp(1:ndim,i)-disp(1:ndim,j)
temp_vec1 = MATMUL(matrix(i)%mat2(1:ndim,1:ndim),&
vec(1:ndim,k))
temp_vec2 = MATMUL(matrix(j)%mat2(1:ndim,1:ndim),&
vec(1:ndim,k))
DO jj=1,ndim
DO ii = 1,ndim
temp_ten1(ii,jj) = -disp_ij(ii) * vec(jj,k)
temp_ten2(ii,jj) = disp_ij(ii) * vec(ndim+jj,k)
END DO
END DO
reduc_ten1(i,1:ndim,1:ndim)=reduc_ten1(i,1:ndim,1:ndim)+temp_ten1*V_0(j) !<--The temporary reduction array
reduc_ten1(j,1:ndim,1:ndim)=reduc_ten1(j,1:ndim,1:ndim)+temp_ten2*V_0(i)
foovec(1:ndim,i) = foovec(1:ndim,i) - temp_vec1(1:ndim)*V_0(j) !<--A generic reduction vector
foovec(1:ndim,j) = foovec(1:ndim,j) + temp_vec1(1:ndim)*V_0(i)
END DO
!$OMP END PARALLEL DO
cend = omp_get_wtime()
! Checking the results
sum_check1 = 0.d0
sum_check2 = 0.d0
DO i = 1,no_partic
matrix(i)%mat2(1:ndim,1:ndim)=reduc_ten1(i,1:ndim,1:ndim) !<--Writing the reduction back to the derived type element
sum_check1 = sum_check1+SUM(foovec(1:ndim,i))
sum_check2 = sum_check2+SUM(matrix(i)%mat2(1:ndim,1:ndim))
END DO
WRITE(*,*) sum_check1, sum_check2, cend-cstart
END PROGRAM TEST_OPEN_MP
The only other alternative I can think of would be to remove all the derived types and replace these with large arrays similar to reduc_ten1 in the example.
Unfortunately, what you want is not possible. At least if I understood your (very complicated for me!) code correctly.
The problem is that you have an array of derived types each have an array. You cannot reference that.
Consider this toy example:
type t
real :: mat(3)
end type
integer, parameter :: n = 100, nk = 1000
type(t) :: parts(n)
integer :: i
real :: array(3,n,nk)
do k = 1, nk
array(:,:,nk) = k
end do
do i = 1, n
parts(i)%mat = 0
end do
!$omp parallel do reduction(+:parts%mat)
do k = 1, nk
do i = 1, n
parts(i)%mat = parts(i)%mat + array(:,i,nk)
end do
end do
!$omp end parallel do
end
Intel Fortran gives a more concrete error:
reduction6.f90(23): error #6159: A component cannot be an array if the encompassing structure is an array. [MAT]
!$omp parallel do reduction(+:parts%mat)
--------------------------------------^
reduction6.f90(23): error #7656: Subobjects are not allowed in this OpenMP* clause; a named variable must be specified. [PARTS]
!$omp parallel do reduction(+:parts%mat)
--------------------------------^
Remember that it is not even allowed to do this, completely without OpenMP:
parts%mat = 0
Intel:
reduction6.f90(21): error #6159: A component cannot be an array if the encompassing structure is an array. [MAT]
gfortran:
Error: Two or more part references with nonzero rank must not be specified at (1)
You must do this:
do i = 1, n
parts(i)%mat = 0
end do
The reason for the error reported by Intel above is very similar.
Actually no derived type components are allowed in the reduction clause, only variable names can be used. That is the reason for the syntax error reported by gfortran. It does not expect any % there. Intel again gives a clearer error message.
But one could make a workaround around that, like passing it to a subroutine and do the reduction there.