使用 Openmp 进行前缀扫描的流压缩（或数组打包）答案

【问题标题】：Stream compaction (or Array Packing) with prefix scan using Openmp使用 Openmp 进行前缀扫描的流压缩（或数组打包）
【发布时间】：2015-01-01 15:54:33
【问题描述】：

我正在使用 openmp 来并行化我的代码。我有一个原始数组：

A=[3,5,2,5,7,9,-4,6,7,-3,1,7,6,8,-1,2]

还有一个标记数组：

M=[1,0,1,0,0,0,1,0,0,1,1,0,0,0,1,1]

使用数组 M i 可以在这个打包数组中压缩我的原始数组：

A=[3,2,-4,-3,1,-1,2]

我想使用多线程方法来解决这个问题。 C++ 的库“推力”解决了这个问题，但我无法为 Fortran 找到类似的工具。是否有一个库，比如 C++ 的“推力”，我可以用来执行流压缩？或者，有没有一种我可以使用 fortran 和 openmp 自己编写的算法来解决这个问题？

【问题讨论】：

我认为你很难编写一个 OpenMP 程序来超越 A = pack(A,M==1)。我认为让多个线程写入A 的开销将扼杀分配packing 工作的任何加速。但我期待被证明是错误的。 Thrust 如何解决问题？
我可以而且也许应该在我之前的评论中添加我知道没有库可以在 Fortran 中实现 pack 内在函数的并行版本。我想你可能会发现从 Thrust 调用 C++ 例程很容易。
如果您的向量非常长，您可以尝试在OMP do 循环中将其分成几块，并在每个子集上使用pack。您需要独立存储生成的子集并在最后合并它们。
首先，感谢您的回复。我不知道 Thrust 如何解决这个问题，但我读过这个库中有很多 API 在 GPU 应用程序的多线程中执行这些操作（减少、前缀和、重新排序等）。这里 (http.developer.nvidia.com/GPUGems3/gpugems3_ch39.html) 似乎有必要进行并行前缀扫描，然后进行分散，但我不明白如何使用 openmp 在 fortran 中编写它。包内在函数是串行的，所以我不知道我是否可以通过大数组和大量线程（MIC 或 GPU）获得更好的性能。
我相信，没有证据证明，您的计算机在主机 RAM 和协处理器 RAM 之间移动非常大的阵列所花费的时间将超过并行化此 pack 可能获得的任何好处荷兰国际集团运营。如果你走这条路，我很想看看你的结果。我认为您指向我们的 URL 中的大多数材料都显示了在 GPU 的 RAM 中获得数据后使用并行算法的优势，但是对于将数据转移到那里和转移回来所花费的时间相当腼腆。

标签： fortran openmp thrust prefix-sum stream-compaction

【解决方案1】：

是否有一个库，比如 C++ 的“thrust”，可以用来执行流压缩？

从 Fortran 调用推力例程应该不难（如果您愿意编写一点 C++ 代码）。此外，推力可以针对 OMP 后端而不是 GPU 后端。

或者，是否有一种我可以使用 fortran 和 openmp 编写自己的算法来解决这个问题？

基本的并行流压缩算法如下。我们将假设最初为您的数据数组中的每个元素分配一个线程。

在M 数组上执行parallel prefix sum (inclusive scan)：

 M=[1,0,1,0,0,0,1,0,0,1,1,0,0,0,1,1]
sM=[1,1,2,2,2,2,3,3,3,4,5,5,5,5,6,7]

然后每个线程将检查其在M 数组中的元素，如果该元素不为零，则将其在A 数组中的对应元素复制到输出数组（我们称之为O）：
```
 M=[1,0,1,0,0,0, 1,0,0, 1,1,0,0,0, 1,1]
sM=[1,1,2,2,2,2, 3,3,3, 4,5,5,5,5, 6,7]
 A=[3,5,2,5,7,9,-4,6,7,-3,1,7,6,8,-1,2]
 O=[3,  2,      -4,    -3,1,      -1,2]
```

如果您在 OMP 中执行此操作，则需要在步骤 1 和步骤 2 之间设置 OMP 屏障。步骤 2 中的工作相对简单且完全独立，因此您可以使用 OMP 并行执行循环，并将工作分解以您希望的任何方式。第 1 步会很复杂，我建议按照你和我链接的章节中提供的大纲进行操作。那里的 OMP 代码将需要各种障碍，但可以并行化。

正如 cmets 中已经提到的，如果这是您想要并行化的唯一部分工作，我不推荐使用 GPU，因为将数据传输到 GPU 或从 GPU 传输数据的成本可能会超过您可能获得的任何并行执行时间收益。但正如我已经提到的，推力可以针对 OMP 实现而不是 GPU 实现。可能值得一试。

关于 fortran 的推力，你需要的大部分是here。这无疑是 CUDA fortran，但唯一的区别应该是不使用设备属性，并且使用推力::主机向量而不是推力::设备向量（至少，开始使用）。

【讨论】：