多个GPU标准C++并行编程加速计算的优点

这是 标准并行编程 系列的第二篇文章，讲述在标准语言中使用并行性来加速计算的优点。

用标准语言并行性开发加速代码

多个 GPU 标准 C ++并行编程，第 2 部分

将应用程序移植到 GPU 的难度因情况而异。在最佳情况下，您可以通过调用现有的 GPU 优化库来加速关键代码部分。例如，当模拟软件的构建块由 BLAS 线性代数函数组成时，可以使用 cuBLAS 对其进行加速。

但在许多代码中，你不能四处做一些手工工作。在这些场景中，您可以考虑使用特定于域的语言，例如 CUDA 来针对特定的加速器。或者，您可以使用基于指令的方法，如 OpenMP 或 OpenACC ，以保持原始语言，并使用相同的代码针对主机和各种类型的设备。

标准并行性

随着本机形式的并行在C++、FORTRAN和 Python 编程语言的现代版本中的出现，现在可以利用类似的高级方法而不需要语言扩展。

C ++中的标准并行性

我们的重点是 C ++语言，它作为 C ++ 17 标准，在标准库中提供了许多算法的并行版本。底层编程模型是前面提到的两种方法的混合体。它的工作方式类似于基于库的方法，因为 C ++提供了用于排序、搜索和累积和的常见任务的并行算法，并且可以在即将到来的版本中添加对特定领域特定算法的支持。此外，以通用for_each和transform_reduce算法的形式提供了并行手写循环的形式。

C ++并行算法通过语言的本地语法来代替非标准扩展来表示并行性。通过这种方式，它们保证了所开发软件的长期兼容性和可移植性。这篇文章还表明，获得的性能通常与 CUDA C ++等传统方法获得的性能相当。

这种方法的新颖之处在于它无缝地集成到现有的代码库中。此方法允许您保留软件体系结构，并有选择地加快关键组件的性能。

将 C ++项目移植到 GPU 可以简单到通过调用for_each或transform_reduce来替换所有的循环，如果它包含了一个还原。

我们将通过典型的重构步骤来克服当前 C ++编译器不兼容问题。这篇文章列出了出于性能原因所需的修改，这些修改更具普遍性，原则上独立于编程形式。这包括以允许合并内存访问的方式重新构造数据的要求。

对于当前编译器， C ++并行算法只针对单个 GPU ，而需要明确的 MPI 并行性来针对多个 GPU 。为此，重用现有并行 CPU 代码的 MPI 后端非常简单，我们将介绍一些实现最先进性能的指导原则。

我们将讨论这些实施规则、指南和最佳实践，并以 Palabos 为例进行说明， Palabos 是一个基于晶格玻尔兹曼方法（ LBM ）的计算流体力学软件库。 PalabOS 在 2021 中被移植到多个 GPU 硬件中，只有几个月的工作，并且说明了对于原来的代码很难适应 GPU 的建议重构步骤的需要。由于广泛使用面向对象的数据结构和编码机制，原始代码对 GPU 的适应性较差。

您知道使用 C ++标准并行性允许混合算法，其中一些算法在 GPU 上执行，但有些算法保持在 CPU 上。这完全取决于它们是否适合 GPU 执行，或者仅仅取决于它们的 GPU 端口的进度状态。这一特性突出显示了 C ++标准并行性的主要优点之一，与保持原有体系结构和大部分软件代码的完整性保持一致。

使用 C ++标准并行编程的 GPU 程序设计

如果这是您第一次听到 C ++并行算法，您可能想读 用标准语言并行性开发加速代码 ，它介绍了 C ++中的标准语言并行性的主题， FORTRAN 和 Python 。

基本概念非常简单。您可以通过执行策略来获得许多标准的 C ++算法在主机上或设备上并行运行，作为一个额外的参数提供给算法。在本文中，我们使用par_unseq execution policy，它表示对不同元素的计算是完全独立的。

以下代码示例执行并行操作，将std：：vector《double》的所有元素乘以 2 ：

for_each（execution：：par_unseq， begin（v）， end（v）， ［］（double& x） { x *= 2.0;

}）;

该算法由nvc++ compiler和-stdpar option编译，在 GPU 上执行。根据编译器、编译器选项和并行算法的实现，还可以在多核 CPU 或其他类型的加速器上获得多线程执行。

此示例使用通用的for_each算法，该算法以函数对象的形式将任何元素操作应用于向量 v 。在本例中，它是一个内联 lambda 表达式。可以使用算法transform_reduce而不是for_each来指定额外的缩减操作。

在for_each算法调用中，调用 lambda 函数时会引用连续的容器元素。但有时，为了访问外部数据数组或实现非本地模板，还必须知道元素的索引。

这可以通过在 C ++ C ++ 17 中的推力库（包括NVIDIA HPC SDK ）和std：：ranges：：views：：iota中提供的counting_iterator迭代来完成，或者在 C ++ 20 中更新。在 C ++ 17 中，最简单的解决方案是从当前元素的地址推导索引。

使用 C ++标准并行性的雅可比示例

为了说明这些概念，下面是一个代码示例，它使用并行 STL 计算非局部模具操作和误差估计的缩减操作。它执行雅可比迭代，计算每个矩阵元素的四个最近邻的平均值：

void jacobi_iteration(vector const& v, vector& tmp) { double const* vptr = v.data(); double *tmp_ptr = tmp.data(); double l2_error = transform_reduce(execution::par_unseq, begin(v), end(v), 0., plus, [=](double& x) { int i = &x - vptr; // Compute index of x from its address. auto [iX, iY] = split(i); double avg = 0.25 * ( vptr[fuse(iX-1, iY)] + vptr[fuse(iX+1, iY)] + vptr[fuse(iX, iY-1)] + vptr[fuse(iX, iY+1)] ); tmp_ptr[i] = avg; return (avg – x) * (avg – x); } );
)

这里，split表示将线性索引i分解为 x 坐标和 y 坐标，fuse则相反。如果域是一个统一的nx-by-ny矩阵， Y 索引在内存中顺序运行，则定义如下：

fuse = ［ny］（int iX， int iY） { return iY + ny * iX; }

split = ［ny］（int i） { return make_tuple（i / ny， i % ny）; }

当算法同时执行时，使用临时向量存储计算出的平均值可以保证确定性结果。

此代码的完整和通用版本可从 gitlab.com/unigehpfs/paralg GitLab 存储库获得。该存储库还包括一个混合版本（ C ++标准并行和 MPI ），它是围绕本文提供的建议构建的，在多个 GPU 上高效运行。

您可能已经注意到，没有明确的语句将数据从主机传输到设备，然后再传输回来。 C ++标准实际上不提供任何这样的语句，并且任何并行算法在一个设备上的实现必须依赖于自动存储器传输。通过NVIDIA HPC SDK ，这是通过 CUDA 统一内存，从 CPU 和 GPU 访问的单个存储器地址空间来实现的。如果代码在 CPU 上访问此地址空间中的数据，然后在 GPU 上访问，则内存页会自动迁移到访问处理器。

对于使用 GPU 加速 CPU 应用程序， CUDA 统一内存特别有用，因为它使您能够专注于以增量方式逐个函数移植应用程序的算法，而无需担心内存管理。

另一方面，隐藏数据传输的性能开销很容易抵消 GPU 的性能优势。通常，在 GPU 上生成的数据应尽可能保存在 GPU 内存中，通过并行算法调用表示其所有操作。这包括数据后处理，如数据统计计算和可视化。如本文 Part 2 所示，它还包括 MPI 通信的数据打包和解包。

按照这篇文章的建议，将代码移植到 GPU 变得非常简单，只需通过调用并行算法来替换所有时间关键的循环和相关的数据访问。不过，最好记住， GPU 通常比 CPU 拥有更多的内核，并且应该暴露在更高级别的并行性中。例如，在下一节介绍的流体动力学问题中，流体域被均匀的、类似矩阵的网格部分覆盖。

在原始的 CPU 代码中，每个 CPU 核心按顺序处理一个或多个网格部分，如图 1 顶部所示。至于 GPU ，网格部分元素上的环路应该并行，以完全占据 GPU 核心。

示例： Lattice Boltzmann 软件和 Palabos

LBM 采用显式时间步格式求解流体流动方程，涵盖了广泛的应用。这包括经过复杂几何形状的流动，如多孔介质、多相流、可压缩超音速流等。

LBM 通常比其他解算器在数值网格的每个节点上分配更多变量。当经典的不可压缩 Navier-Stokes 解算器仅用三个变量表示速度分量，外加一个临时压力项时， LBM 方法通常需要 19 个变量，称为populations。因此， LBM 的内存占用空间要高出 5-6 倍。如果 Navier Stokes 解算器使用临时变量，或者如果系统中添加了进一步的物理量（如密度和温度），实际里程可能会有所不同。

因此，丰富的内存访问和较低的运算强度是 LBM 的特点。在集群级 GPU ，像 NVIDIA V100 和 NVIDIA A100 一样，性能完全受限于内存访问，甚至对于计算密集和复杂的 LBM 方案也是如此。

以 NVIDIA A100 40 GB GPU为例，它具有1555 Gb/s的内存带宽。在每一个明确的时间步长，每个节点访问19个变量或［EZX27 ］，每个都占用八个字节，每个都是双精度的。它们被计数两次：一次用于从 GPU 内存到 GPU 内核的数据传输，另一次用于在计算操作后写回 GPU 内存。

假设一个完美的内存子系统和最大的数据重用， LBM 的峰值吞吐量性能为每秒处理 1555 /（ 19 * 8 * 2 ）= 51.1 亿个网格节点。在 LBM 术语中，通常使用每秒千兆晶格节点更新（ GLUPS ），例如 5.11 GLUPS 。

然而，在现实生活中的应用程序中，每个节点都会额外读取一些信息来管理域异常情况。在 Palabos 中，这是节点标记的 32 位整数和额外数据数组的 64 位索引，有效地将峰值性能降低到 4.92 GLUPS 。

该模型提供了一种简单的方法来估计 LBM 代码可以达到的最佳峰值性能，因为缓存中不适合足够大的网格。我们在整个帖子中使用这个模型来证明用 C ++并行算法获得的性能是一样好的。在几个百分点的差距之外，无论是 CUDA 、 OpenMP ，还是任何其他 GPU 形式主义，都不能做得更好。

LBM 巧妙地区分了由局部collision step表示的计算和封装在streaming step中的内存传输操作。以下代码示例显示了具有矩阵式拓扑结构的结构化网格的典型时间迭代：

for (int i = 0; i < N; ++i) { // Fetch local populations "f" from memory. double f_local[19]; for (int k = 0; k < 19; ++k) { f_local[k] = f[i][k]; } collide(f_local); // Execute collision step. // Write data back to neighboring nodes in memory (streaming). auto [iX, iY, iZ] = split(i); for (int k = 0; k < 19; ++k) { int nb = fuse(iX+c[k][0], iY+c[k][1], iZ+c[k][2]); ftmp[nb][k] = f_local[k]; }
}

与前一节中的 Jacobi 迭代一样，该函数将计算出的数据写入temporary array ftmp，以避免多线程执行期间出现争用情况，这使其成为演示本文概念的理想候选。有 替代就地算法 可以避免内存复制。然而，它们更复杂，因此不太适合用于说明目的。

自然过程的模拟和建模 课程介绍了LBM。有关如何使用C++并行算法开发 GPU 的LBM代码的更多信息，请参见多核格子玻尔兹曼模拟的跨平台编程模型。

在本文中，我们使用 开源 LBM 库 Palabos 来展示如何将现有的 C ++库用并行算法移植到多 GPU 。乍一看， Palabos 似乎不适合 GPU 端口，因为它强烈依赖面向对象的机制。然而，在 Palabos 的案例中，我们将介绍几种变通方法，这些方法只需对代码体系结构进行表面的更改即可实现最先进的性能。

从面向对象设计转向面向数据设计

为了服务于大型社区，Palabos强调多态性和其他面向对象技术。包含数据（种群）和方法（局部碰撞模型）的对象代表每个网格节点。这为开发新模型提供了一个方便的API，并提供了一个灵活的机制来调整模型从一个单元到另一个单元的物理行为或数值方面。

然而，由于数据布局效率低下、执行路径复杂，以及对虚拟函数调用的依赖，这种面向对象的方法不太适合在 GPU 上执行。以下几节将教您如何通过采用开发模型，以 GPU 友好的方式重构代码，我们在总括术语下称之为data-oriented programming。

摆脱基于类的多态性

图 2 的左半部分展示了 Palabos 中网格节点上冲突模型的典型代码执行链。算法的不同组件被叠加起来，并通过虚拟函数调用，包括底层的数值 LBM 算法（ RR ）、附加物理（“ Smagorinsky ”）和附加数值方面（左边界）。

这个面向对象设计的教科书案例变成了代码 GPU 端口的责任。这个问题的出现是因为当前版本的 HPCSDK 不支持 C ++并行算法中的 GPU 的虚函数调用。一般来说，出于性能原因，在 GPU 上应该避免这种类型的设计，因为它限制了执行路径的可预测性。

简单的解决方法是将执行链收集到单个函数中，这些函数按顺序显式调用各个组件，并用唯一的标记标识它们。

在 Palabos 中，这个独特的标记是在序列化机制的帮助下生成的，该机制最初是为了支持动态自适应仿真的检查点和网络通信而开发的。这表明，如果重构软件项目的体系结构足够灵活，那么 GPU 端口的大部分重构工作都是自动完成的。

现在，您可以为每个网格节点提供一个标记，标识完整冲突步骤的代码，并用一个大的 switch 语句表示冲突步骤：

switch（tag） { case rr_les： fun_rr_les（f_local）; break; case rr_les_BCleft： fun_rr_les_BCleft（f_local）; break; …

}

随着 switch 语句变大，由于生成的内核在 GPU 内存中占用的空间，它可能会遇到性能问题。

另一个问题是软件项目的可维护性。目前，如果不修改这个 switch 语句（它是库核心的一部分），就不可能提供新的碰撞模型。这两个问题的解决方案在于，在编译时使用 C ++模板机制在最终用户应用程序中生成具有选定数量的实例的切换语句。 Palabos-GPU 资源页详细介绍了这种技术。

重新安排内存以鼓励合并的内存访问

面向对象的设计还会导致内存布局无法在 GPU 的多核架构上高效处理。当每个节点对 LBM 方案的 19 个局部种群进行分组时，数据以结构阵列（ AoS ）的形式结束（图 3 ）。为了简单起见，每个节点只显示四个总体。

AoS 数据布局导致性能不佳，因为它阻止了流式处理步骤中的合并内存访问，而流式处理步骤由于非本地模板，在内存访问方面是算法最关键的部分。

数据应该以阵列结构（ SoA ）的方式对齐，在流式处理步骤中进行通信的给定类型的所有群体都在连续的内存地址上对齐。在这种重新安排之后，即使是一个相对简单的 LBM 算法也能获得接近典型 GPU 内存带宽的 80% 。

面向数据的设计意味着您将中心重要性赋予数据的结构和布局，并围绕此结构构建数据处理算法。面向对象的方法通常采用反向路径。

GPU 端口的第一步应该是了解应用程序的理想数据布局。就 LBM 而言， GPU 上 SoA 布局的优越性是众所周知的事实。内存布局和内存遍历算法的细节在之前发布的案例研究 开源 STLBM 代码 中进行了测试。

结论

在这篇文章中，我们讨论了使用 C ++标准并行编程编写 GPU 应用程序的基本技术。我们还提供了晶格玻尔兹曼方法和 Palabos 应用程序的背景信息，我们在案例研究中使用了这些信息。最后，我们讨论了两种方法，可以重构源代码，使其更适合在 GPU 上运行。

在下一篇文章中，我们继续使用这个应用程序，并讨论如何在NVIDIA上运行时在C++应用程序中获得高性能。我们还演示了如何通过MPI扩展应用程序以使用多个 GPU 。

关于作者

Jonas Latt 是瑞士日内瓦大学计算机科学系的副教授。他从事高性能计算和计算流体力学的研究，并在包括地球物理、生物医学和航空航天领域在内的跨学科领域进行应用。他是 lattice Boltzmann 复杂流动模拟开源软件 Palabos 的最初开发者和当前共同维护者。他以前在日内瓦大学获得物理学和计算机科学博士学位，并通过塔夫斯大学（波士顿，美国）和综合理工学校 F.E.EdRaelde 洛桑 EPFL （瑞士）的研究，并作为 CFD 公司 FuluKIT 的联合创始人，对流体力学感兴趣。

Christophe Guy Coreixas 是一名航空工程师， 2014 年毕业于 ISAE-SUPAERO （法国图卢兹）。 2018 年，他在 CERFACS 从事面向行业应用的可压缩晶格玻尔兹曼方法研究时获得了博士学位（流体动力学）。作为日内瓦大学计算机科学系的博士后，克里斯多夫现在开发格子玻尔兹曼模型来模拟航空、多物理和生物医学流程。

Gonzalo Brito 是 NVIDIA 计算性能与 HPC 团队的高级开发技术工程师，工作于硬件和软件的交叉点。他热衷于让加速计算变得更容易实现。在加入NVIDIA 之前，冈萨洛在 RWTH 亚琛大学空气动力学研究所开发了多物理方法，用于颗粒流。

Jeff Larkin 是 NVIDIA HPC 软件团队的首席 HPC 应用程序架构师。他热衷于高性能计算并行编程模型的发展和采用。他曾是 NVIDIA 开发人员技术小组的成员，专门从事高性能计算应用程序的性能分析和优化。 Jeff 还是 OpenACC 技术委员会主席，曾在 OpenACC 和 OpenMP 标准机构工作。在加入NVIDIA 之前，杰夫在位于橡树岭国家实验室的克雷超级计算卓越中心工作。

审核编辑：郭婷