有限差分法 - 拉普拉斯方程（第4部分）

Finite difference method - Laplacian part 4 — ROCm Blogs (amd.com)

2023年7月18日 作者： Justin Chang, Thomas Gibson, Sean Miller.

在前三篇文章中，我们展示了使用HIP实现的拉普拉斯算子的有限差分法，并进行了四种不同的优化。在这些代码修改过程中，由于减少了对全局内存的总读取次数，我们观察到了渐进的性能提升。然后我们应用了进一步的优化，在MI250X GPU的单个GCD上，在`512`×`512`×`512`点的网格上达到预期的性能目标。下面展示的是最终的HIP内核代码，我们称之为Kernel 5：

// 瓦片因子
#define m 8
// 启动界限
#define LB 256
template <typename T>
__launch_bounds__(LB)
__global__ void laplacian_kernel(T * f, const T * u, int nx, int ny, int nz, T invhx2, T invhy2, T invhz2, T invhxyz2) {

    int i = threadIdx.x + blockIdx.x * blockDim.x;
    int j = m*(threadIdx.y + blockIdx.y * blockDim.y);
    int k = threadIdx.z + blockIdx.z * blockDim.z;

    // 如果线程在xz边界，则退出
    if (i == 0 || i >= nx - 1 ||
        k == 0 || k >= nz - 1)
        return;

    const int slice = nx * ny;
    size_t pos = i + nx * j + slice * k;

    // 每个线程在y方向上累加m个模板
    T Lu[m] = {0};

    // 可重用数据的标量
    T center;

    // z - 1，循环瓦片
    for (int n = 0; n < m; n++)
        Lu[n] += u[pos - slice + n*nx] * invhz2;

    // y - 1
    Lu[0]   += j > 0 ? u[pos - 1*nx] * invhy2 : 0; // bound check

    // x 方向，循环瓦片
    for (int n = 0; n < m; n++) {
        // x - 1
        Lu[n] += u[pos - 1 + n*nx] * invhx2;

        // x
        center = u[pos + n*nx]; // store for reuse
        Lu[n] += center * invhxyz2;

        // x + 1
        Lu[n] += u[pos + 1 + n*nx] * invhx2;

        // reuse: y + 1 for prev n
        if (n > 0) Lu[n-1] += center * invhy2;

        // reuse: y - 1 for next n
        if (n < m - 1) Lu[n+1] += center * invhy2;
    }

    // y + 1
    Lu[m-1]  += j < ny - m ? u[pos + m*nx] * invhy2 : 0; // bound check

    // z + 1，循环瓦片
    for (int n = 0; n < m; n++)
      Lu[n] += u[pos + slice + n*nx] * invhz2;

    // 仅当线程在y边界内时存储
    for (int n = 0; n < m; n++)
      if (n + j > 0 && n + j < ny - 1)
        __builtin_nontemporal_store(Lu[n],&f[pos + n*nx]);
}

第4部分的目标是探讨这一代码在各种AMD GPU和问题规模上的性能表现。作为此次研究的一部分，我们将通过总结经验教训来结束拉普拉斯系列，并提供用户可以应用于类似问题的代码优化建议。本文其余部分的结构如下：

1. 在多个AMD GPU和各种问题规模上运行Kernel 5
2. 检查L2缓存命中性能，并识别其退化的根本原因
3. 提出规避L2缓存大小限制的技术
4. 回顾在本拉普拉斯系列中进行的所有代码优化。

在不同硬件和尺寸上的性能研究

除了 MI250X GPU，我们还考虑了以下 AMD GPU：

1. RX 6900 XT

2. RX 7900 XTX

3. MI50

4. MI100

5. MI210

6. MI250

我们将对 Kernel 5 进行一个缩放研究，从问题规模 nx,ny,nz = 256, 256, 256 开始，跨上述所有 GPU。每个维度会以 256 的倍数增加，直到达到 nx,ny,nz = 1024, 1024, 1024 的全局问题规模。请记住，评判标准 (FOM) 定义为有效内存带宽：

theoretical_fetch_size = ((nx * ny * nz - 8 - 4 * (nx - 2) - 4 * (ny - 2) - 4 * (nz - 2) ) * sizeof(double);
theoretical_write_size = ((nx - 2) * (ny - 2) * (nz - 2)) * sizeof(double);
effective_memory_bandwidth = (theoretical_fetch_size + theoretical_write_size) / average_kernel_execution_time;

这是考虑整个内存子系统中需要移动的最少数据量的平均数据传输速率。目标是获得尽可能接近设备 HBM 带宽的 FOM。通过在 512 × 512 × 512 的点网格上对 Kernel 5 进行各种线程块尺寸、启动边界和瓦片因子的多次实验，我们发现以下参数为每个单独 GPU 提供了最高的 FOM：

GPU	线程块	启动边界	瓦片因子
RX 6900 XT	256 x 1 x 1	256	16
RX 7900 XTX	256 x 1 x 1	256	8
MI50	256 x 1 x 1	256	4
MI100	128 x 1 x 1	128	4
MI210	256 x 1 x 1	256	8
MI250	256 x 1 x 1	256	8
MI250X	256 x 1 x 1	256	8

此后，MI250 和 MI250X GPU 上的实验将仅在一个 GCD 上进行。下图展示了各种 AMD GPU 和问题规模下的 FOM[1]。

图1：Kernel 5 在各种 AMD GPU 和问题规模下的性能

请注意，RX 6900 XT GPU 没有足够的内存来运行最大的问题规模 nx,ny,nz = 1024, 1924, 1024 ，因此图中未显示。RX 6900 XT 和 RX 7900 XTX GPU 在所有尺寸上的性能相对一致，而所有 Instinct™ GPU 在问题规模跨越某个阈值时表现出性能下降。在接下来的章节中，我们将调查 Instinct™ GPU 性能下降的原因。

性能下降的根本原因

让我们重新关注MI250X GPU。为了理解性能下降的原因，我们需要收集`rocprof`统计数据。首先，我们将收集`FETCH_SIZE`并将其与每种`nx, ny, nz`组合的`theoretical_fetch_size`进行比较。理想情况下，`FETCH_SIZE`应该与理论值相等，因为每个线程加载的元素可以被相邻线程重复使用最多六次。需要注意的是，MI100、MI210、MI250和MI250X GPU的每个GCD都共享8MB的L2缓存，因此我们收集到的`L2CacheHit`将有助于确定缓存数据的重用情况。下图显示了我们在各种问题规模上的发现：

图2：在单个MI250X GCD上的Kernel 5在各种问题规模下的取数效率和L2缓存命中率

FOM（有效内存带宽）、取数效率和L2缓存命中率之间存在明显的相关性。回顾从Laplacian Part 1帖子得到的基线HIP内核在保持相对较高的L2缓存命中率（65%）的同时表现出了50%的取数效率。这让我们相信内核设计存在问题。在最后的实验中，取数效率和L2缓存命中率都稳定在33%左右，这表明性能的限制因素是L2缓存的大小。

从上图还可以看到一个有趣的现象：性能的下降仅发生在`nx`或`ny`增加时——nz的增加对性能几乎没有影响。从Kernel 5中的线程块配置和网格索引可以看出，所有来自单个xy平面的元素将在获取来自下一个xy平面的元素之前，从全局内存中获取。每个线程的模板计算在z方向上有`-nx*ny`、`0`和`nx*ny`的跳跃，因此需要三个xy平面的数据。如果三个xy平面过大而不能同时适应L2缓存，那么每个线程将循环缓存以从全局内存中获取额外的xy平面数据——例如，当只有一个xy平面适合（或部分适合）L2缓存时，线程需要循环通过两个额外的xy平面从全局内存中获取数据，因此取数效率下降到33%。

让我们进行第二次实验来进一步说明这一点。我们将检查不同的`nx`和`ny`值，以确定在性能开始下降之前的最佳z方向跨度。以约1 MB的起始z跨度（相当于`nx,ny = 512,256`）并固定`nx`为三个不同的大小，我们逐步增加`ny`，直到xy平面超过8MB（`nx*ny*sizeof(double)`）。见下图：

图3：在各种固定的`nx`值上，增加`nx*ny`跨度对单个MI250X GCD性能的影响

无论`nx`的大小，性能都会在单个xy平面超过约2.5 MB内存时下降。事实证明，2.5 MB足够小，可以使三个连续的xy平面适应8 MB的L2缓存，这解释了FOM值较高的原因。然而，跨度越大，更多的元素开始从缓存中掉出，从而降低L2缓存命中率并随后降低FOM。当`nx*ny`超过8 MB后，性能趋于平稳，因为此时只有一个平面适合（或部分适合）缓存——在这种情况下，`FETCH_SIZE`比理论估计值大了三倍。

关于Radeon™ GPU性能的说明

关于Radeon™ GPU性能的说明

让我们快速进行一项新的缩放研究，以验证两款Radeon™ GPU的LLC确实是影响较大规模问题性能的限制因素。为了确保不耗尽内存，我们从`nx,ny,nz = 1024,1024,32`开始，逐渐增加`ny`的值：

图4：随着`nx*ny`跨度增加，对RX 6900 XT和RX 7900 XTX GPU性能的影响

在上图中，当xy平面跨越某些阈值时，有效内存带宽仍会恶化。RX 6900 XT的性能在xy平面超过32 MB时恶化，而RX 7900 XTX的性能在xy平面超过25 MB时恶化。RX 7900 XTX的LLC比RX 6900 XT稍小，因此其性能开始下降的xy平面尺寸较小。无论使用哪种AMD GPU，当`nx * ny`太大时，性能都会下降。

解决缓存大小限制

那么我们如何解决这些较大问题时的缓存大小问题呢？这在像MI50这样只有4 MB L2缓存的硬件上尤其重要。即使`nz`很小，较大的`nx`和`ny`值，L2缓存也可能成为内存流量瓶颈。让我们考虑当`nx,ny,nz = 1024,1024,1024`时，在MI250X GCD上进行的几次实验。以下是解决L2缓存瓶颈的两种可能策略：

选项1：网格块配置

如前所述，当前的线程块和网格索引配置会先获取x方向的元素，然后是y方向和z方向。对于`nx,ny,nz = 1024,1024,1024`这样的大问题，这种方法会在任何z方向的模板被获取之前完全填满L2缓存。让我们首先通过将基线线程块配置从256 × 1 × 1修改为128 × 1 × 8来进行尝试。这需要我们将启动边界重置回默认值1024，但让我们看看这个简单的变化是否对性能有任何影响：

	FOM (GB/s)	获取效率 (%)	L2 缓存命中率
Kernel 5 - baseline	555.389	33.1	34.4
Kernel 5 - 128x1x8	934.111	79.8	60.6

有了显著的改进 - FOM提高了68%，获取效率接近80%，L2缓存命中率提高到超过60%。将线程块配置修改为包括z方向的元素使得HIP内核在这些较大问题规模上显著提高了性能。然而，获取效率仍然有大约20%的提升空间。

另一种修改方法是仍使用原始的256 × 1 × 1，但是配置HIP内核内的网格启动和索引。这需要在主机侧和设备侧对代码进行修改。对于主机侧，我们修改网格启动配置：

Kernel 5 (之前的配置)	Kernel 6 (之后的配置)
dim3 grid( (nx - 1) / block.x + 1, (ny - 1) / (block.y * m) + 1, (nz - 1) / block.z + 1);	dim3 grid( (ny - 1) / (block.y * m) + 1, (nz - 1) / block.z + 1, (nx - 1) / block.x + 1);

第二个变化涉及配置HIP内核内部的索引：

Kernel 5 (之前的配置)	Kernel 6 (之后的配置)
int i = threadIdx.x + blockIdx.x * blockDim.x; int j = m*(threadIdx.y + blockIdx.y * blockDim.y); int k = threadIdx.z + blockIdx.z * blockDim.z;	int i = threadIdx.x + blockIdx.z * blockDim.x; int j = m*(threadIdx.y + blockIdx.x * blockDim.y); int k = threadIdx.z + blockIdx.y * blockDim.z;

之前，`blockIdx.x`和`blockIdx.y`索引按x和y方向分别跨越`pos`，这将大小为`nx * ny`的xy平面加载到缓存中。现在，它们按y和z方向分别跨越`pos`，这将大小为`blockDim.x * ny`的xy平面加载到缓存中。索引顺序的重新排列使多个小块xy平面可以填充L2缓存，从而增加了某些模板计算有效重用其他线程块已经获取的数据的可能性。以下是性能数据：

	FOM (GB/s)	获取效率 (%)	L2 缓存命中率
Kernel 5 - 基线	555.389	33.1	34.4
Kernel 5 - 128x1x8	934.111	79.8	60.6
Kernel 6 - 网格索引	1070.04	95.4	66.2

这种修改后的网格索引实际上比之前修改线程块大小的性能更好。获取效率非常接近理论限制，L2缓存命中率略有改善。然而，这个解决方案也不是完美的 - 最后一个例子是用`blockDim.x = 256`运行的，因此每个子xy平面大约占用2 MB的内存，从而仍然允许三个不同的子xy平面适合L2缓存。较大的`nx * ny`平面将不可避免地面临与之前相同的问题。

选项 2：将整体问题拆分成子域

另一种更可靠的方法是将整体问题拆分成较小的子域，使得三个 xy 平面能够适合 L2 缓存，并串行执行内核启动。当前 xy 平面的大小 nx,ny = 1024,1024 是之前考虑的适合 L2 缓存的三个 nx,ny = 512,512 xy 平面的 4 倍，因此我们考虑将这个问题分解成 4 个子域。

让我们回到 Kernel 5。首先我们修改主机端代码：

Kernel 5 (之前)	Kernel 7 (之后)
dim3 grid( (nx - 1) / block.x + 1, (ny - 1) / (block.y * m) + 1 , (nz - 1) / block.z + 1); ... laplacian_kernel<<<grid, block>>>(d_f, d_u, nx, ny, nz, invhx2, invhy2, invhz2, invhxyz2);	int bny = (ny - 1) / 4 + 1; dim3 grid( (nx - 1) / block.x + 1, (bny - 1) / (block.y * m) + 1 , (nz - 1) / block.z + 1); ... for (int i = 0; i < 4; i++) laplacian_kernel<<<grid, block>>>(d_f, d_u, nx, ny, nz, invhx2, invhy2, invhz2, invhxyz2, bny*i);

代码中做了三处修改。首先，我们通过使用 bny = (ny - 1) / 4 + 1 而不是 ny 修改 grid.y 值，以表示我们只在 y 方向上遍历域的四分之一。接下来，我们添加一个 for 循环启动四个 HIP 内核，以在每个子域上计算模板。最后，我们通过添加 y 方向的偏移量来修改内核参数。设备内核需要如下修改：

Kernel 5 (之前)	Kernel 7 (之后)
__global__ void laplacian_kernel(T * f, const T * u, int nx, int ny, int nz, T invhx2, T invhy2, T invhz2, T invhxyz2) { ... int j = m*(threadIdx.y + blockIdx.y * blockDim.y);	__global__ void laplacian_kernel(T * f, const T * u, int nx, int ny, int nz, T invhx2, T invhy2, T invhz2, T invhxyz2, int sy) { ... int j = sy + m*(threadIdx.y + blockIdx.y * blockDim.y);

内核中的其他内容可以保持不变。我们来看看以下性能数据：

性能指标	FOM (GB/s)	取回效率 (%)	L2 缓存命中率
Kernel 5 - 基线	555.389	33.1	34.4
Kernel 5 - 128x1x8	934.111	79.8	60.6
Kernel 6 - 网格索引	1070.04	95.4	66.2
Kernel 7 - 子域拆分	1175.54	99.6	67.7

FOM 比之前更高。取回效率几乎达到 100%，这表明我们优化了内存传输。与之前的方法不同，这种方法能处理任何大小的 xy 平面——xy 平面越大，用户可以将问题拆分成更多的子域。如果您的内核要在 Radeon™ 和 Instinct™ GPU 上运行并且不同的 LLC 尺寸会显著影响性能的“临界点”，则这种方法特别有用。我们选择在 y 方向拆分问题大小，但是也可以轻松应用这种优化到 x 方向。

总结与回顾

在这一部分，我们首先检查了 Kernel 5 在各种 AMD GPU 和问题大小上的性能，并注意到当网格跨越某个大小阈值时，性能会急剧下降。通过实验确定，LLC（最后一级缓存）的大小是大尺寸 xy 平面的性能瓶颈。为了解决缓存大小问题，我们提出了两种不同的解决方法，这两种方法都只需对代码进行少量修改。

在整个有限差分方法拉普拉斯算子系列中, 我们从一个简单的 HIP 内核实现开始，逐步应用了几种不同的优化，极大地提高了内核的性能。我们现在将快速回顾在此工作中进行的所有不同优化:

1. 循环分块(Loop tiling) (Part 2): 如果已知内核会加载多个可重用值，可以考虑增加循环分块（即循环展开）。这需要重构代码的几个步骤，但最终减少了启动的线程块数量并增加了每个线程计算的模板数量。通过这种方法，可以大大减少 FETCH_SIZE。

2. 重新排序的访问模式(Reordered access patterns) (Part 2): 经常前后访问设备内存中的元素可能会过早地从缓存中驱逐可重用数据。内核应重新排序内核中的所有加载指令，以单向访问数据，即按升序地址访问内存。为了正确利用循环分块优化，需要进行此优化。

3. 启动界限(Launch bounds) (Part 3): 如果内核的寄存器使用量非常高，应用启动界限可能会使编译器在寄存器分配方面做出更合适的决策，从而改善寄存器压力。默认情况下启动线程数为 1024，但我们强烈建议将其设置为内核将使用的线程数量。

4. 非时间性加载(Nontemporal loads) (Part 3): 因为输出数组 f 的元素不会被其他线程或内核共享，所以可以使用内建的内在函数（intrinsics）将有限差分模板计算结果写回设备内存，并在从缓存中驱逐时优先级更高，从而为输入数组 u 腾出更多的缓存空间。

5. 网格块配置(Grid block configuration) (Part 4): 一个快速解决溢出 L2 缓存问题的方法是配置线程块或网格索引配置。这可以使 xy 平面的较小段首先填满缓存，从而让一些线程能够执行计算而不必多次从全局内存中提取数据。

6. 子域划分(Splitting into subdomains) (Part 4): 如果整个问题受到 LLC 大小的瓶颈限制，可以考虑将问题划分为多个子域，并依次启动 HIP 内核来处理每个子域。这种方法非常稳健，可以根据需要调整以适应任何 LLC 或问题大小。

在结构化网格上实现拉普拉斯算子的有限差分模板在 HIP 中相对简单。在未来的实验笔记中，我们将考虑其他偏微分方程，比如声波方程，这可能需要比这里展示的更高阶的模板。尽管如此，我们相信，在这个拉普拉斯算子系列中列出的所有优化策略都应该被利用，以获得模板内核的最佳性能。

伴随的代码示例

如果您有任何问题或意见，请在 GitHub 上留言讨论

---

[1]

测试使用 ROCm 版本 5.3.0-63 进行。基准测试结果不是经过验证的性能数据，仅用于演示代码修改的相对性能改进。实际性能结果取决于多个因素，包括系统配置和环境设置，结果的可重复性不保证。