我有一个程序，可以将粒子分布到单元中的云网格中。简单地循环粒子总数(Ntot)并填充一个256^3网格(即每个粒子分布在8个单元上)。
% gfortran -fopenmp cic.f90 -o ./cic
哪个编译很好。但是当我运行它(./cic)时，我遇到了分段错误。我我的循环是一个经典的 omp do 问题。当我不在 openmp 中编译它时，该程序可以工作。
!$omp parallel do
do i = 1,Ntot
if (x1(i).gt.0.and.y1(i).gt.0.and.z1(i).gt.0) then
dense(int(x1(i)),int(y1(i)),int(z1(i))) = dense(int(x1(i)),int(y1(i)),int(z1(i))) &
+ dx1(i) * dy1(i) * dz1(i) * mpart
end if
if (x2(i).le.Ng.and.y1(i).gt.0.and.z1(i).gt.0) then
dense(int(x2(i)),int(y1(i)),int(z1(i))) = dense(int(x2(i)),int(y1(i)),int(z1(i))) &
+ dx2(i) * dy1(i) * dz1(i) * mpart
end if
if (x1(i).gt.0.and.y2(i).le.Ng.and.z1(i).gt.0) then
dense(int(x1(i)),int(y2(i)),int(z1(i))) = dense(int(x1(i)),int(y2(i)),int(z1(i))) &
+ dx1(i) * dy2(i) * dz1(i) * mpart
end if
if (x2(i).le.Ng.and.y2(i).le.Ng.and.z1(i).gt.0) then
dense(int(x2(i)),int(y2(i)),int(z1(i))) = dense(int(x2(i)),int(y2(i)),int(z1(i))) &
+ dx2(i) * dy2(i) * dz1(i) * mpart
end if
if (x1(i).gt.0.and.y1(i).gt.0.and.z2(i).le.Ng) then
dense(int(x1(i)),int(y1(i)),int(z2(i))) = dense(int(x1(i)),int(y1(i)),int(z2(i))) &
+ dx1(i) * dy1(i) * dz2(i) * mpart
end if
if (x2(i).le.Ng.and.y1(i).gt.0.and.z2(i).le.Ng) then
dense(int(x2(i)),int(y1(i)),int(z2(i))) = dense(int(x2(i)),int(y1(i)),int(z2(i))) &
+ dx2(i) * dy1(i) * dz2(i) * mpart
end if
if (x1(i).gt.0.and.y2(i).le.Ng.and.z2(i).le.Ng) then
dense(int(x1(i)),int(y2(i)),int(z2(i))) = dense(int(x1(i)),int(y2(i)),int(z2(i))) &
+ dx1(i) * dy2(i) * dz2(i) * mpart
end if
if (x2(i).le.Ng.and.y2(i).le.Ng.and.z2(i).le.Ng) then
dense(int(x2(i)),int(y2(i)),int(z2(i))) = dense(int(x2(i)),int(y2(i)),int(z2(i))) &
+ dx2(i) * dy2(i) * dz2(i) * mpart
end if
end do
!$omp end parallel do
迭代之间没有依赖关系。想法？
这个问题以及您另一个问题中的问题来自这样一个事实，即启用 OpenMP 时自动堆数组被禁用。这意味着在没有-fopenmp的情况下，大数组会自动放置在静态存储中(称为.bss段)，而小数组则被分配到堆栈上。当您打开 OpenMP 支持时，不会使用自动静态分配，并且您的dense数组将在例程堆栈上分配。OS X 上的默认堆栈限制非常严格，因此存在分段错误。
您在这里有几个选择。第一个选项是通过赋予denseSAVE属性来使其具有静态分配。另一种选择是通过使其ALLOCATABLE然后使用ALLOCATE语句在堆上显式分配它，例如：
REAL, DIMENSION(:,:,:), ALLOCATABLE :: dense
ALLOCATE(dense(256,256,256))
! Computations, computations, computations
DEALLOCATE(dense)
较新的 Fortran 版本支持在数组超出范围时自动释放没有SAVE属性的数组。
请注意，您的 OpenMP 指令很好，不需要额外的数据共享子句。不需要在PRIVATE子句中声明i，因为循环计数器具有预先确定的私有数据共享类。您不需要将其他变量放在子句SHARED因为它们是隐式共享的。然而，您在dense上执行的操作应该与ATOMIC UPDATE同步(或者只是在较旧的 OpenMP 实现上ATOMIC)，或者您应该使用REDUCTION(+:dense)。原子更新被转换为锁定的添加，与循环内有条件的巨大减速相比，不应该引起太大的减速：
INTEGER :: xi, yi, zi
!$OMP PARALLEL DO PRIVATE(xi,yi,zi)
...
if (x1(i).gt.0.and.y1(i).gt.0.and.z1(i).gt.0) then
xi = int(x1(i))
yi = int(y1(i))
zi = int(z1(i))
!$OMP ATOMIC UPDATE
dense(xi,yi,zi) = dense(xi,yi,zi) &
+ dx1(i) * dy1(i) * dz1(i) * mpart
end if
...
复制代码，并对其他情况进行适当的更改。如果您的编译器抱怨ATOMIC构造中的UPDATE子句，只需将其删除即可。
REDUCTION(+:dense)将在每个线程中创建dense的一个副本，这将消耗大量内存，并且最终应用的减少会随着dense的大小而变得越来越慢。对于小数组，它会比原子更新工作得更好。
有关如何使变量共享和私有的说明，请参阅 https://computing.llnl.gov/tutorials/openMP/#Clauses。
看起来所有变量都应该共享，除了循环变量i必须是私有的。 建议使用以下行：
!$omp parallel do default(shared) private(i)
这应该可以解决您的分段错误(假设我所有变量都正确
但是，存在不同的线程尝试同时覆盖dense的相同部分的风险，从而导致总数不正确。为了防止这种情况，您需要将每个作业包装在dense!$omp atomic或!$omp critical部分中。
但是，您可能会发现这些关键部分将导致线程花费大部分时间等待，因此您可能不会看到与纯串行代码相比有任何改进。
原则上，您可以通过使用reduction关键字声明dense来解决此问题，但不幸的是它不能用于数组。