是否可以对派生的 Fortran 类型的元素进行 OpenMP 缩减？

Question

我正在尝试修改 Fortran 代码 (Gfortran) 以使用 OpenMP。它是一种基于粒子的代码，其中数组的索引可以对应于粒子或粒子对。该代码使用派生类型为每个粒子存储多个矩阵。遇到需要使用存储在此派生类型中的矩阵的循环是很常见的。这个矩阵可以被多个线程访问。该循环还需要对该派生类型中的元素进行缩减。我目前必须编写一个临时数组才能进行此缩减，然后将派生类型的元素设置为等于此临时缩减数组。如果不使用 OpenMP，则不需要临时数组。

问题：是否可以对派生类型的元素进行归约？我不认为我可以对整个派生类型进行缩减，因为我需要访问派生类型中的某些元素才能工作，这意味着它需要共享。（通过阅读 specification 我了解到，当使用 REDUCTION 时，会创建每个列表项的私有副本。）

完成下面的最小工作示例。它可能会更小，但我担心删除更多组件可能会过度简化问题。

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

我能想到的唯一其他选择是删除所有派生类型，并将其替换为类似于示例中 reduc_ten1 的大型数组。

Answer 1

很遗憾，您想要的是不可能的。至少如果我正确理解你的（对我来说非常复杂！）代码。

问题是您有一个派生类型数组，每个派生类型都有一个数组。你不能引用它。

考虑这个玩具示例：

  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 给出了更具体的错误：

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)
--------------------------------^

请记住，完全没有 OpenMP 甚至不允许这样做：

parts%mat = 0

英特尔：

 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)

你必须这样做：

  do i = 1, n
    parts(i)%mat = 0
  end do

上面Intel报错的原因很相似

实际上reduction子句中不允许使用派生类型组件，只能使用变量名。这就是 gfortran 报告语法错误的原因。它不希望那里有任何 % 。 Intel再次给出更清晰的错误信息。

但是可以解决这个问题，比如将其传递给子例程并在那里进行归约。

是否可以对派生的 Fortran 类型的元素进行 OpenMP 缩减？

Is it possible to do OpenMP reduction over an element of a derived Fortran type?

fortran

openmp

gfortran

reduction

derived-types