使用 OpenMP 原子捕获操作获取粒子 3D 直方图并创建索引的竞争条件

Question

使用 OpenMP 原子捕获操作获取粒子 3D 直方图并创建索引的竞争条件

And*_*And 3 c++ multithreading atomic openmp histogram

我的完整代码中有一段代码：

\n

const unsigned int GL=8000000;\nconst int cuba=8;\nconst int cubn=cuba+cuba;\nconst int cub3=cubn*cubn*cubn;\nint Length[cub3];\nint Begin[cub3];\nint Counter[cub3];\nint MIndex[GL];\nstruct Particle{\n  int ix,jy,kz;\n  int ip;\n};\nParticle particles[GL];\nint GetIndex(const Particle & p){return (p.ix+cuba+cubn*(p.jy+cuba+cubn*(p.kz+cuba)));}    \n...\n#pragma omp parallel for\nfor(int i=0; i<cub3; ++i) Length[i]=Counter[i]=0;\n#pragma omp parallel for\nfor(int i=0; i<N; ++i)\n{\n  int ic=GetIndex(particles[i]);\n  #pragma omp atomic update\n  Length[ic]++;\n}\nBegin[0]=0;\n#pragma omp single\nfor(int i=1; i<cub3; ++i) Begin[i]=Begin[i-1]+Length[i-1];\n#pragma omp parallel for\nfor(int i=0; i<N; ++i)\n{\n  if(particles[i].ip==3)\n  {\n    int ic=GetIndex(particles[i]);\n    if(ic>cub3 || ic<0) printf("ic=%d out of range!\\n",ic);\n    int cnt=0;\n  #pragma omp atomic capture\n    cnt=Counter[ic]++;\n    MIndex[Begin[ic]+cnt]=i;\n  }\n}\n

Run Code Online (Sandbox Code Playgroud)\n

如果要删除

\n

#pragma omp parallel for\n

Run Code Online (Sandbox Code Playgroud)\n

代码工作正常，输出结果始终相同。\n但是使用此编译指示，代码中存在一些未定义的行为/竞争条件，因为每次它都会给出不同的输出结果。\n如何解决此问题？

\n

更新：任务如下。有很多带有一些随机坐标的粒子。需要将粒子的数组粒子在坐标系的每个单元（笛卡尔立方体，例如1\xc3\x971\xc3\x971 cm）中的索引输出到数组MIndex中。因此，在MIndex的开头，应该在坐标系的第一个单元格中的粒子的数组粒子中存在索引，然后 - 在第二个单元格中，然后 - 在第三个单元格中，依此类推。区域MIndex中给定单元内的索引顺序并不重要，可以是任意的。如果可能的话，需要并行进行，可以使用原子操作。

\n
有一种直接的方法：并行遍历所有坐标单元，并在每个单元中检查所有粒子的坐标。但对于大量细胞和粒子来说，这似乎很慢。有更快的方法吗？是否可以并行遍历粒子数组一次并使用原子操作填充MIndex数组，就像上面的代码片段中所写的那样？
\n

Answer 1

Pet*_*des 5

如果您想要一个可以高效工作的算法（不需要共享计数器上的原子 RMW，这将是一场灾难，请参见下文），您可能无法让编译器为您自动并行化标量代码。但您也许可以使用 OpenMP 作为启动线程和获取线程 ID 的方式。

保留初始直方图中的每线程计数数组，在第二遍中使用

（更新：这可能行不通：我if(particles[i].ip==3)之前没有注意到源代码中的。我假设它Count[ic]会像Length[ic]串行版本中一样高。如果情况不是这样，这个策略可能会留下空白或其他东西。但正如 Laci 所说指出，也许您希望在计算长度时首先进行检查，那么就可以了。）

手动对第一个直方图进行多线程处理（进入Length[]），每个线程处理已知的i值范围。保留这些每个线程的长度，即使您对它们求和并通过前缀求和来构建 Begin[]。

Length[thread][ic]该立方体中的粒子数量也是如此，超出了i该线程处理的值范围。（并且将在第二个循环中再次循环：关键是我们以相同的方式将线程之间的粒子划分两次。理想情况下，同一线程在同一范围内工作，因此 L1d 缓存中的事物可能仍然很热。）

将其预处理到每个线程Begin[][]数组中，以便每个线程知道 MIndex 中的位置将每个存储桶中的数据放入其中。

// pseudo-code, fairly close to actual C
for(ic < cub3) {
   // perhaps do this "vertical" sum into a temporary array
   // or prefix-sum within Length before combining across threads?
   int pos = sum(Length[0..nthreads-1][ic-1]) + Begin[0][ic-1];

   Begin[0][ic] = pos;
   for (int t = 1 ; t<nthreads ; t++) {
       pos += Length[t][ic];   // prefix-sum across threads for this cube bucket
       Begin[t][ic] = pos;
   }
}

Run Code Online (Sandbox Code Playgroud)

这有一个非常糟糕的缓存访问模式，特别是在cuba=8makeLength[t][0]和Length[t+1][0]4096 字节彼此分开的情况下。（因此 4k 别名是一个可能的问题，缓存冲突未命中也是如此）。

也许每个线程都可以将自己的 Length 切片前缀求和到 Begin 切片中，1. 用于缓存访问模式（以及局部性，因为它刚刚写入了这些长度），2. 为该工作获得一些并行性。

然后在最后一个循环中MIndex，每个线程都可以int pos = --Length[t][ic]从长度中派生出唯一的 ID。（就像您对所做的那样Count[]，但没有将另一个每线程数组引入为零。）

Length 的每个元素都将返回零，因为同一个线程正在查看它刚刚计算的相同点。Begin[t][ic]如果位置计算正确，MIndex[...] = i商店就不会发生冲突。错误共享仍然是可能的，但它是一个足够大的数组，点往往会分散在周围。

不要过度使用线程数，特别是当cuba线程数大于 8 时。长度/开始预处理工作量随线程数而变化，因此最好为不相关的线程或任务留出一些空闲 CPU完成一些吞吐量。OTOH，这cuba=8意味着每个线程数组只有 4096 字节（太小而无法并行化归零，顺便说一句），实际上并没有那么多。

（编辑之前的先前回答更清楚地说明了发生了什么。）

这基本上是一个直方图吗？如果每个线程都有自己的计数数组，您可以在最后将它们相加（您可能需要手动执行此操作，而不是让 OpenMP 为您执行此操作）。但似乎您还需要这个计数在每个体素中都是唯一的，以便正确更新 MIndex？这可能是一个令人震惊的事情，就像要求调整每个 MIndex 条目（如果可能的话）。

更新后，您将在 Length[] 中绘制直方图，以便可以加快该部分的速度。

原子 RMW 对于您的代码来说是必要的，性能灾难

共享计数器的原子增量会更慢，并且在 x86 上可能会严重破坏内存级并行性。在 x86 上，每个原子 RMW 也是一个完整的内存屏障，在发生之前耗尽存储缓冲区，并阻止稍后的加载从开始到发生之后。

与使用非原子访问的单个线程相反，单个线程可能会缓存未命中多个元素Counter，Begin并且元素未完成。MIndex（由于无序 exec，下一次迭代的 load / inc / store forCounter[ic]++可以执行加载，同时 forBegin[ic]和/或Mindex[]store 存在未完成的缓存未命中。）

允许宽松原子增量的 ISA 可能能够有效地做到这一点，例如 AArch64。（同样，OpenMP 可能无法为您做到这一点。）

即使在 x86 上，拥有足够的（逻辑）核心，您仍可能获得一定的加速，特别是如果Counter访问足够分散，核心不会不断争夺相同的缓存线。不过，您仍然会得到大量缓存行在内核之间跳跃，而不是在 L1d 或 L2 中保持热状态。（虚假分享是一个问题，

也许软件预取可以提供帮助，例如prefetchw（写入预取）计数器 5 或 10i次迭代后。

不过，即使有memory_order_seq_cst增量，哪个点按哪个顺序进行也不是确定性的。无论哪一个线程首先递增，Counter[ic]都是将 thatcnt与 that关联起来的线程i。

替代访问模式

也许让每个线程扫描所有点，但只处理其中的一个子集，且子集不相交。Counter[]因此，任何给定线程接触的元素集仅由该线程接触，因此增量可以是非原子的。

按范围过滤p.kz可能最有意义，因为这是索引中最大的乘数，因此每个线程“拥有”一个连续的Counter[].

但如果你的观点分布不均匀，你就需要知道如何分解事情以大致平均地分配工作。而且您不能将其划分得更细（如 OMP 动态调度），因为每个线程都会扫描所有点：这会增加过滤工作量。

也许几个固定分区是一个很好的权衡，可以在不引入大量额外工作的情况下获得一些并行性。

回复：您的编辑

你已经循环遍历整个点数组了吗Length[ic]++;？使用再次执行相同的直方图工作似乎是多余的Counter[ic]++;，但如何避免它并不明显。

计数数组很小，但如果完成后不需要两者，您可能只需递减 Length 即可为体素中的每个点分配唯一索引。至少第一个直方图可以受益于每个线程的不同计数数组的并行化，以及在最后垂直相加。由于计数数组对于 L1d 缓存来说足够小，因此应该可以与线程完美地扩展。

顺便说一句，for() Length[i]=Counter[i]=0;太小了，不值得并行化。对于cuba=8，它8*8*16 * sizeof(int)= 4096 字节，只是一页，所以它只是两个小的内存集。

（当然，如果每个线程都有自己单独的 Length 数组，则每个线程都需要将其归零）。如果一长串点都在同一个存储桶中，那么它足够小，甚至可以考虑展开每个线程可能有 2 个计数数组，以隐藏存储/重新加载串行依赖关系。最后组合计数数组是一项工作，#pragma omp simd或者只是正常的自动向量化，gcc -O3 -march=native因为它是整数工作。

对于最后一个循环，您可以将点数组分成两半（将一半分配给每个线程），并让一个线程通过从向下计数来获取唯一 ID --Length[i]，另一个线程从中的 0 开始向上计数Counter[i]++。通过不同的线程查看不同的点，这可以使您的速度提高 2 倍。（MIndex 存储的模争用。）

要做的不仅仅是向上和向下计数，您还需要从整个 Length 数组中获得的信息……但您暂时拥有这些信息。请参阅顶部部分

归档时间：	3 年，8 月前
查看次数：	214 次
最近记录：	3 年，8 月前