11

我有一个更新矩阵 A 的循环，我想让它成为 openmp 但我不确定哪些变量应该共享和私有。我原以为只有 ii 和 jj 会起作用，但事实并非如此。我想我也需要一个 !$OMP ATOMIC UPDATE 某处......

该循环仅计算 N 和 N-1 个粒子之间的距离并更新矩阵 A。

            !$OMP PARALLEL DO PRIVATE(ii,jj)
            do ii=1,N-1
                    do jj=ii+1,N
                            distance_vector=X(ii,:)-X(jj,:)
                            distance2=sum(distance_vector*distance_vector)
                            distance=DSQRT(distance2)
                            coff=distance*distance*distance
                            PE=PE-M(II)*M(JJ)/distance
                            A(jj,:)=A(jj,:)+(M(ii)/coff)*(distance_vector)
                            A(ii,:)=A(ii,:)-(M(jj)/coff)*(distance_vector)
                    end do
            end do
            !$OMP END PARALLEL DO
4

2 回答 2

25

OpenMP 的黄金法则是，在外部范围中定义的所有变量（有一些排除项）默认在并行区域中共享。由于在 2008 年之前的 Fortran 中没有本地作用域（即BLOCK ... END BLOCK在早期版本中没有），所有变量（除了变量threadprivate）都是共享的，这对我来说很自然（不像 Ian Bush，我不喜欢使用default(none)然后重新声明各种复杂科学代码中所有 100 多个局部变量的可见性）。

以下是如何确定每个变量的共享类别：

  • N- 共享，因为它在所有线程中应该是相同的，并且它们只读取它的值。
  • ii- 它是循环的计数器，受制于工作共享指令，因此其共享类预定为private。在子句中明确声明它并没有什么坏处PRIVATE，但这并不是真正必要的。
  • jj- 循环的循环计数器，不受工作共享指令的约束，因此jj应该是private.
  • X- 共享，因为所有线程都引用并仅从中读取。
  • distance_vector- 显然应该是private因为每个线程都在不同的粒子对上工作。
  • distance, distance2, 和coff- 同上。
  • M- 出于与X.
  • PE- 充当累加器变量（我猜这是系统的势能）并且应该是归约操作的主题，即应该放在一个REDUCTION(+:....)子句中。
  • A- 这个很棘手。它可以通过A(jj,:)同步构造共享和更新以保护，或者您可以使用归约（OpenMP 允许在 Fortran 中对数组变量进行归约，这与 C/C++ 不同）。A(ii,:)不会被多个线程修改，因此不需要特殊处理。

随着减少A到位，每个线程都将获得它的私有副本，A这可能会占用内存，尽管我怀疑您是否会使用这种直接 O(N 2 ) 模拟代码来计算具有大量粒子的系统。还存在与缩减实现相关的一定开销。在这种情况下，您只需要添加A到REDUCTION(+:...)子句列表中。

使用同步构造，您有两种选择。您可以使用ATOMIC构造或CRITICAL构造。由于ATOMIC仅适用于标量上下文，您必须“取消向量化”赋值循环并ATOMIC分别应用于每个语句，例如：

!$OMP ATOMIC UPDATE
A(jj,1)=A(jj,1)+(M(ii)/coff)*(distance_vector(1))
!$OMP ATOMIC UPDATE
A(jj,2)=A(jj,2)+(M(ii)/coff)*(distance_vector(2))
!$OMP ATOMIC UPDATE
A(jj,3)=A(jj,3)+(M(ii)/coff)*(distance_vector(3))

您也可以将其重写为循环——不要忘记声明循环计数器private。

CRITICAL无需对循环进行非矢量化：

!$OMP CRITICAL (forceloop)
A(jj,:)=A(jj,:)+(M(ii)/coff)*(distance_vector)
!$OMP END CRITICAL (forceloop)

命名关键区域是可选的，在这种特殊情况下有点不必要，但通常它允许分离不相关的关键区域。

哪个更快？用ATOMICorCRITICAL展开 这取决于很多事情。通常CRITICAL要慢得多，因为它通常涉及对 OpenMP 运行时的函数调用，而原子增量（至少在 x86 上）是使用锁定加法指令实现的。正如他们常说的，YMMV。

概括地说，你的循环的工作版本应该是这样的：

!$OMP PARALLEL DO PRIVATE(jj,kk,distance_vector,distance2,distance,coff) &
!$OMP& REDUCTION(+:PE)
do ii=1,N-1
   do jj=ii+1,N
      distance_vector=X(ii,:)-X(jj,:)
      distance2=sum(distance_vector*distance_vector)
      distance=DSQRT(distance2)
      coff=distance*distance*distance
      PE=PE-M(II)*M(JJ)/distance
      do kk=1,3
         !$OMP ATOMIC UPDATE
         A(jj,kk)=A(jj,kk)+(M(ii)/coff)*(distance_vector(kk))
      end do
      A(ii,:)=A(ii,:)-(M(jj)/coff)*(distance_vector)
   end do
end do
!$OMP END PARALLEL DO

我假设您的系统是 3 维的。

说了这么多，我赞同 Ian Bush，你需要重新思考位置和加速度矩阵是如何在内存中布局的。适当的缓存使用可以提升您的代码，还可以允许某些操作，例如X(:,ii)-X(:,jj)向量化，即使用向量 SIMD 指令实现。

于 2013-01-17T14:30:13.130 回答
3

如所写，您将需要一些同步以避免竞争条件。考虑 2 线程的情况。假设线程 0 以 ii=1 开头，因此认为 jj=2,3,4,...。线程 1 以 ii=2 开头，因此认为 jj=3,4,5,6。因此，正如所写的那样，线程 0 可能正在考虑 ii=1,jj=3，而线程 1 同时正在查看 ii=2,jj=3。这显然会导致生产线出现问题

                        A(jj,:)=A(jj,:)+(M(ii)/coff)*(distance_vector)

因为两个线程具有相同的 jj 值。所以是的，您确实需要将更新同步到 A 以避免比赛，尽管我必须承认我的好方法对我来说并不是立即显而易见的。如果有什么事情发生在我身上，我会考虑并编辑。

但是我还有其他 3 条评论：

1）您的内存访问模式很糟糕，我希望纠正这个问题至少可以提供与任何 openmp 一样多的速度，而且麻烦少得多。在 Fortran 中，您希望以最快的速度降低第一个索引 - 这确保内存访问在空间上是本地的，从而确保充分利用内存层次结构。鉴于这是在现代机器上获得良好性能的最重要的事情，您应该真正努力做到这一点。因此，如果您可以排列数组以便将上述内容写为

        do ii=1,N-1
                do jj=ii+1,N
                        distance_vector=X(:,ii)-X(:jj)
                        distance2=sum(distance_vector*distance_vector)
                        distance=DSQRT(distance2)
                        coff=distance*distance*distance
                        PE=PE-M(II)*M(JJ)/distance
                        A(:,jj)=A(:,jj)+(M(ii)/coff)*(distance_vector)
                        A(:,ii)=A(:,ii)-(M(jj)/coff)*(distance_vector)
                end do
        end do

请注意这是如何降低第一个索引的，而不是您所拥有的第二个索引。

2) 如果你确实使用 openmp，我强烈建议你使用 default(None)，它有助于避免讨厌的错误。如果你是我的学生之一，你会因为不这样做而失去大量分数！

3) Dsqrt 是过时的——在现代 Fortran（即 1967 年之后的任何东西）中，除了一些晦涩的情况外，sqrt 已经足够好，而且更灵活

于 2013-01-16T15:15:50.183 回答