这是标准并行编程系列文章的第二篇，讲述了在标准语言中使用并行来加速计算的优势。第一篇文章：《英伟达是如何做 GPU 编程的（一）》
将应用程序移植到 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 （如果它包含归约的话）来替换所有的 for 循环。
我们将带你通过完成典型的重构步骤来克服当前 C++ 编译器不兼容的问题。这篇文章列出了出于性能原因所进行的修改，这些修改更具普遍性，原则上与编程形式无关。它还包含了以支持合并内存访问方式重构数据的需求。
对于当前的编译器，C++ 并行算法只针对单个 GPU，而显式的 MPI 并行则需要针对多个 GPU。为此，重用现有并行 CPU 代码的 MPI 后端是很简单的，我们将提供一些准则，以实现最先进的性能。
我们将讨论这些实施规则、指导方针和最佳实践，并以 Palabos 为例进行说明，Palabos 是一个基于格子玻尔兹曼方法（Lattice Boltzmann Method，LBM）的计算流体动力学软件库。Palabos 在 2021 年只用了几个月的时间就被移植到了多 GPU 硬件上，这很好地说明了建议的重构步骤的必要性，因为原始代码对 GPU 的适应性很差。原始代码之所以对 GPU 的适应性较差，是由于广泛使用面向对象的数据结构和编码机制。
通过本文，你将会了解到使用 C++ 标准并行能允许采用混合算法，其中一些算法在 GPU 上执行，而另一些则保留在 CPU 上。这完全取决于它们对 GPU 执行的适应性，或者仅仅取决于它们 GPU 移植的进度状态。这一特性是 C++ 标准并行的主要突出优势之一，与保持全局体系结构和大部分软件代码基础完整的能力相当。
使用 C++ 标准并行进行 GPU 编程
如果这是你第一次听说 C++ 并行算法，你可能需要阅读一下 《英伟达是如何做 GPU 编程的（一）》，它介绍了 C++、Fortran 和 Python 中标准语言并行的相关主题。
基本概念非常简单。在执行策略的帮助下，可以让许多标准的 C++ 算法在主机或设备上并行运行，执行策略是作为算法的额外参数提供的。在本文中，我们将使用 par_unseq 执行策略，它表示不同元素的计算是完全独立的。
下面的代码示例执行了一个并行操作，将 std::vector<double> 中的所有元素乘以 2：
for_each(execution::par_unseq, begin(v), end(v), [](double& x) {
x *= 2.0;
});该算法使用 nvc++ 编译器和 -stdpar选项进行编译，在 GPU 上执行。根据编译器、编译器选项和并行算法的实现，你还可以在多核 CPU 或其他类型的加速器上实现多线程执行。
这个例子使用了通用的 for_each 算法，它以函数对象的形式应用于向量 v 中的任何元素操作。在本例中，它是一个内联 lambda 表达式。还可以使用 transform_reduce 算法而不是 for_each 来指定额外的归约（reduction）操作。
在 for_each 算法调用中，lambda 函数被调用时引用了连续的容器元素。但有时，为了访问外部数据数组或实现非局部模板（non-local stencil），还必须要知道元素的索引。
这可以通过对 counting_iterator 进行迭代来实现，counting_iterator 由 C++ 17（包含在 NVIDIA HPC SDK 中）中的 Thrust 库及 C++20 或更新版本中的 std::ranges::views::iota提供。在 C++17 中，最简单的解决方案是通过当前元素的地址推导出索引。
使用 C++ 标准并行的雅克比示例
为了说明这些概念，下面是一个代码示例，它使用并行 STL 计算非局部模板操作（non-local stencil operation）和误差估计的归约操作。它执行雅克比（Jacobi）迭代，计算每个矩阵元素中四个最近邻值的平均值：
void jacobi_iteration(vector<double> const& v, vector<double>& 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>, [=](double& x)
{
int i = &x - vptr; // 通过地址计算 x 的索引
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++ 标准实际上并未提供任何这样的语句，并且任何并行算法在一个设备上的实现必须依赖于自动内存传输（automatic memory transfers）。在 NVIDIA HPC SDK 中，这是通过 CUDA 统一内存实现的，CPU 和 GPU 都可以访问一个单一的内存地址空间。如果你的代码先在 CPU 上访问这个地址空间中的数据，然后再在 GPU 上访问该数据，那么内存页会被自动迁移到该访问处理器中。
对于使用 GPU 来加速 CPU 的应用程序，CUDA 统一内存特别有用，因为它使你能够专注于以增量的方式来移植应用程序算法，即逐个逐个函数的移植，而无需担心内存管理。
另一方面，隐藏数据传输的性能开销很容易抵消 GPU 的性能优势。通常，在 GPU 上生成的数据应尽可能地保存在 GPU 内存中，方法是通过并行算法调用来表达其所有操作。这包括数据后置处理（post-processing），如数据统计计算和可视化。如本文的第 2 部分所示，它还包括用于 MPI 通信的数据打包和解包操作。
按照本文的建议，将代码移植到 GPU 就变得非常简单了，只需通过调用并行算法来替换所有时序要求严格（time-critical）的循环和相关的数据访问即可。不过，请记住，GPU 通常比 CPU 拥有更多的内核，并且应该具有更高级别的并行性。例如，在下一节介绍的流体动力学问题中，流体域被均匀的矩阵状的网格分块覆盖了（图 1）。
图 1 在 Palabos 中，将流体域细分为矩阵状网格分块，以及在 CPU 和 GPU 上各自的并行化策略
在原始的 CPU 代码中，每个 CPU 内核按顺序处理一个或多个网格分块，如图 1 的顶部所示。至于 GPU，网格分块上元素的循环应该并行化以充分占用 GPU 内核。
示例：格子玻尔兹曼软件和 Palabos
格子玻尔兹曼方法（LBM）采用显式的时间步骤方案来求解流体流动方程，涵盖了广泛的应用。它包含了流经复杂几何形状的流动，如多孔介质、多相流动、可压缩超音速流动等。
LBM 通常在数值网格的每个节点上分配比其他求解器更多的变量。经典的不可压缩的 Navier-Stokes 求解器表示状态仅用 4 个变量，3 个标识速度分量的变量，再加上 1 个标识临时压力项的变量，LBM 方式通常需要 19 个变量，称为 populations 。因此，LBM 的内存占用要高出 5 到 6 倍。如果 Navier-Stokes 求解器使用临时变量，或者如果在系统中添加了其他物理参数（例如密度和温度），那么实际里程可能会有所不同。
因此，丰富的内存访问和较低的运算强度是 LBM 的特征。在集群级 GPU（如 NVIDIA V100 和 NVIDIA A100）上，性能完全受限于内存访问，即使是对于计算密集型和复杂的 LBM 方案也是如此。
以 NVIDIA A100 40 GB GPU 为例，它的内存带宽为 1555GB/s。在每个显式的时间步骤中，每个节点访问 19 个变量或 populations ，每个变量或 populations 以双精度形式占用 8 个字节。它们被计数两次：一次用于从 GPU 内存到 GPU 内核的数据传输，另一次用于计算操作完成后回写 GPU 内存。
假设一个完美的内存子系统和最大的数据重用，LBM 的峰值吞吐量性能为每秒处理 1555/（1982）=51.1 亿个网格节点。在 LBM 术语中，通常使用每秒千兆晶格节点更新（Giga Lattice-node Updates Per Second，GLUPS），如 5.11 GLUPS。
不过，在实际应用程序中，每个节点还需要额外读取一些信息来管理函数域的异常情况。在 Palabos 中，节点标记是一个 32 位整数，另一个数据数组是一个 64 位索引，它们能有效地将峰值性能降低到 4.92 GLUPS。
该模型提供了一种简单的方法来估计 LBM 代码可以达到的最佳峰值性能，因为缓存中不能容纳足够大的网格。在这篇文章中，我们使用这个模型来证明用 C++ 并行算法能所获得的性能是最好的。除了几个百分点的差距，无论是 CUDA、OpenMP，还是任何其他的 GPU 形式，都不能做得更好。
LBM 巧妙地区分了由局部碰撞步骤（local collision step）表示的计算和封装在流式处理步骤（local collision step）中的内存传输操作。以下代码示例展示了具有矩阵式拓扑结构的结构化网格的典型时间迭代：
for (int i = 0; i < N; ++i) {
// 从内存中获取本地 populations“f”。
double f_local[19];
for (int k = 0; k < 19; ++k) {
f_local[k] = f[i][k];
}
collide(f_local); // 执行碰撞步骤
// 将数据回写到内存中的相邻节点 (流)。
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）迭代一样，该函数将计算出的数据写入临时数组 ftmp ，以避免多线程执行期间出现竞争条件，这使它成为演示本文概念的理想候选项。还存在其他避免内存重复的原地算法（In-place Algorithm）。然而，这些更为复杂，因此不太适用于演示的目的。
在 《自然过程的模拟与建模》的课程 中有对 LBM 的介绍。有关如何使用 C++ 并行算法为 GPU 开发 LBM 代码的更多信息，请参见 《多核格子玻尔兹曼仿真的跨平台编程模型》。
在本文中，我们使用 开源的 LBM 库 Palabos 来展示如何将现有的 C++ 库移植到具有并行算法的多 GPU 上。乍一看，Palabos 似乎不适合 GPU 移植，因为它严重依赖于面向对象的机制。然而，在 Palabos 的示例中，我们将介绍几种变通的方案，这些变通的方案只需对代码体系架构进行一些表面的更改即可获得最先进的性能。
从面向对象设计转向面向数据设计
为了服务于大型社区，Palabos 强调多态性和其他面向对象的技术。用一个同时包含数据（populations）和方法（本地碰撞模型）的对象来表示每个网格节点。这为新模型的开发提供了一个便利的 API，并提供了一种灵活的机制来调整模型从一个单元到另一个单元的物理行为或数值相位。
然而，这种面向对象的方法不适合在 GPU 上执行，因为其数据布局效率低，执行路径复杂，并且依赖于虚拟函数调用。下面几节内容将会教你如何通过采用开发模型，以对 GPU 友好的方式重构代码，我们将其称为面向数据的编程。
摆脱基于类的多态性
图 2 的左侧展示了 Palabos 中网格节点上碰撞模型的典型代码执行链。该算法的不同组件被叠加起来，并通过虚拟函数调用，包括底层的数值 LBM 算法（RR）、附加物理场（“Smagorinsky”）和附加数值相位（Left Boundary）。
图 2 提取单个整数标记以识别虚函数调用链
这个面向对象设计的教科书级案例变成了为 GPU 移植代码的一个不利因素。这个问题的出现是因为当前版本的 HPC SDK 不支持在 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++ 模板机制在最终用户应用程序中生成具有选定数量 case 的 switch 语句。这种技术的更多详细信息请参见 Palabos-GPU 资源页面。
重新排列内存以鼓励合并内存访问
面向对象的设计还会导致内存布局无法在 GPU 的多核架构上高效处理。由于每个节点对 LBM 方案的 19 个局部 populations 进行了分组，数据最终以结构数组（array-of-structure，AoS）的形式排列（图 3）。为了简单起见，每个节点只显示四个 populations。
图 3 虽然在面向对象的设计中，结构数组（array-of-structure，AOS）的对齐更为自然，但数组结构（structure-of-array，SOA）改进了 LBM 流处理步骤中的合并内存访问
AoS 数据布局会导致性能低下，因为它阻止了流式处理步骤中的合并内存访问，由于是非局部模板（non-local stencil），所以流式处理步骤是算法中内存访问方面最关键的部分。
数据应该以数组结构（SoA）的方式对齐，在流式处理步骤期间通信的所有给定类型的 populations 都在连续的内存地址上对齐了。在重新排列之后，即使是一个相对简单的 LBM 算法也能获得接近典型 GPU 内存带宽的 80%。
面向数据的设计意味着你将重点放在数据的结构和布局上，并围绕该结构构建数据处理算法。面向对象的方法通常采用相反的路径。
GPU 移植的第一步应该是了解应用程序的理想数据布局。在 LBM 的案例中，SoA 布局在 GPU 上的优势是众所周知的。内存布局和内存遍历算法的细节在之前的一个案例研究中已经进行了测试，该案例研究是以开源 STLBM 代码的形式发布的。
结 论
在本文中，我们讨论了使用 C++ 标准并行编程编写 GPU 应用程序所涉及的基本技术。我们还提供了格点玻尔兹曼方法和 Palabos 应用程序的背景信息，我们在案例研究中使用了这些信息。最后，我们讨论了两种重构源代码以使其更适合在 GPU 上运行的方法。
在下一篇文章中，我们将继续讨论这个应用程序，并讨论当 C++ 应用程序运行在 NVIDIA GPU 上时，如何获得高性能。我们还将演示如何扩展应用程序以使用多个 GPU 和 MPI。
作者介绍
Jonas Latt
Jonas Latt 是瑞士日内瓦大学计算机科学系的副教授。他从事高性能计算和计算流体力学的研究，并在包括地球物理、生物医学和航空航天领域在内的跨学科领域进行应用。他是格子玻尔兹曼复杂流动模拟开源软件 Palabos 的最初开发者和当前共同维护者。此前，他曾在日内瓦大学获得了物理学和计算机科学博士学位，并通过在塔夫茨大学（美国波士顿）和洛桑联邦理工学院 EPFL（瑞士）的研究而对流体动力学很感兴趣，他同时也是 CFD 公司 FlowKit 的联合创始人。
Christophe Guy Coreixas
Christophe Guy Coreixas 是一名航空工程师，2014 年毕业于法国图卢兹的 ISAE-SUPAERO。2018 年，他在 CERFACS 从事面向工业应用的可压缩格子玻尔兹曼方法研究时获得了博士学位 (流体动力学)。作为日内瓦大学计算机科学系的博士后，Christophe 现在开发了格子玻尔兹曼模型来模拟航空、多物理场和生物医学流动。
Gonzalo Brito
Gonzalo Brito 是英伟达计算性能与高性能计算（Compute Performance & HPC）团队的高级开发技术工程师，负责硬件和软件的交叉部分。他对让加速计算变得更容易实现充满热情。在加入英伟达之前，Gonzalo 在德国亚琛工业大学（RWTH Aachen University）空气动力学研究所开发了用于颗粒流动的多物理场方法。
Jeff Larkin
Jeff 是英伟达高性能计算（HPC）软件团队的首席 HPC 应用架构师。他对高性能计算并行编程模型的发展和采用充满热情。他曾是英伟达开发技术小组的成员，专门从事高性能计算应用的性能分析和优化。Jeff 还是 OpenACC 技术委员会的主席，他曾在 OpenACC 和 OpenMP 标准机构工作过。在加入英伟达之前，Jeff 曾在位于橡树岭国家实验室（Oak Ridge National Laboratory）的 Cray 超级计算卓越中心工作过。
本文转自 公众号：AI前线 ，作者Jonas Latt、Christophe Guy Coreixas、Gonzalo Brito、Jeff Larkin，点击阅读原文
