OpenMP 的黄金法则是，在外部作用域中定义的所有变量（有一些例外）默认在并行区域中共享。由于 2008 年之前的 Fortran 中没有局部作用域（即没有BLOCK ... END BLOCK在早期版本中），所有变量（除了threadprivate）是共享的，这对我来说很自然（与伊恩布什不同，我不太喜欢使用default(none)然后在各种复杂的科学代码中重新声明所有 100 多个局部变量的可见性）。
以下是确定每个变量的共享类的方法：
-
N- 共享，因为它在所有线程中应该是相同的，并且它们只读取它的值。
-
ii- 它是循环的计数器，受工作共享指令的约束，因此其共享类预先确定为private。在 a 中明确声明它并没有什么坏处PRIVATE条款，但这并不是真正必要的。
-
jj- 循环的循环计数器，不受工作共享指令的约束，因此jj应该private.
-
X- 共享，因为所有线程都引用并且只能从中读取。
-
distance_vector- 显然应该是private因为每个线程都作用于不同的粒子对。
-
distance, distance2, and coff- 同上。
-
M- 应该出于与以下相同的原因共享X.
-
PE- 充当累加器变量（我猜这是系统的势能）并且应该是归约操作的主题，即应该放入REDUCTION(+:....) clause.
-
A- 这个很棘手。它可以被共享和更新A(jj,:)使用同步构造进行保护，或者您可以使用归约（与 C/C++ 不同，OpenMP 允许在 Fortran 中对数组变量进行归约）。A(ii,:)永远不会被多个线程修改，因此不需要特殊处理。
With reduction over A in place, each thread would get its private copy of A and this could be a memory hog, although I doubt you would use this direct O(N2) simulation code to compute systems with very large number of particles. There is also a certain overhead associated with the reduction implementation. In this case you simply need to add A to the list of the REDUCTION(+:...) clause.
通过同步构造，您有两种选择。您可以使用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.
With CRITICAL无需取消循环向量化：
!$OMP CRITICAL (forceloop)
A(jj,:)=A(jj,:)+(M(ii)/coff)*(distance_vector)
!$OMP END CRITICAL (forceloop)
命名关键区域是可选的，在这种特殊情况下有点不必要，但通常它允许分离不相关的关键区域。
哪个更快？展开与ATOMIC or CRITICAL？这取决于很多事情。通常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
我假设你的系统是三维的。
话虽如此，我同意伊恩·布什的观点，你需要重新思考位置和加速度矩阵在内存中的布局方式。正确的缓存使用可以增强您的代码，并且还允许某些操作，例如X(:,ii)-X(:,jj)进行矢量化，即使用矢量 SIMD 指令实现。
