有限差分方法 - 拉普拉斯算子第二部分

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

2023年1月4日 作者：Justin Chang, Rajat Arora, Thomas Gibson, Sean Miller, Ossian O’Reilly

在之前的拉普拉斯算子文章中，我们开发了一种基于HIP实现的有限差分模板代码，专门用于拉普拉斯算子。初步实现发现该代码受限于内存带宽，也就是说其运行时间受限于我们移动数据到全局内存和从全局内存中提取数据的速率。此外，目前的内存访问模式需要多次访问全局内存来加载所有数据，因此如果我们缓存更多的数据，执行时间可以减少。我们将性能指标（FOM）定义为_有效内存带宽_，即理论上的内存传输量除以实际执行时间。我们HIP实现的FOM目前达到了单个MI250X GCD峰值的50%[1]，但我们的分析表明，如果将实际内存传输量降低到理论值，我们的FOM可以达到峰值的至少71%[1]. 在这篇文章中，我们将介绍两种常见的优化技术，这些技术可以应用于内核，以帮助实现这一目标：

1. 循环分块以显著减少内存加载
2. 重新排序内存访问模式以改善缓存性能

回顾

在上一篇文章中，我们讨论了基于HIP实现拉普拉斯算子的中心有限差分模板。回顾一下，拉普拉斯算子的形式是标量场u(x,y,z)的梯度的散度：

∇⋅∇u=∇²u=∂²u/∂x² + ∂²u/∂y² + ∂²u/∂z²

最初的HIP实现如下所示：

template <typename T>
__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 = threadIdx.y + blockIdx.y * blockDim.y;
    int k = threadIdx.z + blockIdx.z * blockDim.z;

    // Exit if this thread is on the boundary
    if (i == 0 || i >= nx - 1 ||
        j == 0 || j >= ny - 1 ||
        k == 0 || k >= nz - 1)
        return;

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

    // Compute the result of the stencil operation
    f[pos] = u[pos] * invhxyz2
           + (u[pos - 1]     + u[pos + 1]) * invhx2
           + (u[pos - nx]    + u[pos + nx]) * invhy2
           + (u[pos - slice] + u[pos + slice]) * invhz2;
}

template <typename T>
void laplacian(T *d_f, T *d_u, int nx, int ny, int nz, int BLK_X, int BLK_Y, int BLK_Z, T hx, T hy, T hz) {

    dim3 block(BLK_X, BLK_Y, BLK_Z);
    dim3 grid((nx - 1) / block.x + 1, (ny - 1) / block.y + 1, (nz - 1) / block.z + 1);
    T invhx2 = (T)1./hx/hx;
    T invhy2 = (T)1./hy/hy;
    T invhz2 = (T)1./hz/hz;
    T invhxyz2 = -2. * (invhx2 + invhy2 + invhz2);

    laplacian_kernel<<<grid, block>>>(d_f, d_u, nx, ny, nz, invhx2, invhy2, invhz2, invhxyz2);
}

以及在单个MI250X GCD上的相应性能：

$ ./laplacian_dp_kernel1
Kernel: 1
Precision: double
nx,ny,nz = 512, 512, 512
block sizes = 256, 1, 1
Laplacian kernel took: 2.64172 ms, effective memory bandwidth: 808.148 GB/s

报告的808.148 GB/s 达到了我们在上一篇文章中设定目标性能的69.4%。我们将继续使用`rocprof`来帮助我们评估以下优化的有效性。

在上一篇文章中，我们可以推断出`f`设备数组是高效存储的，因为`WRITE_SIZE`指标与理论值一致。然而，报告的`FETCH_SIZE`几乎是理论值的两倍，因此如果该数据量可以减少，性能可能会有所提升。需要记住的是，从`u`数组加载的七个元素用于更新`f`数组中的每一个条目，但这些`u`元素最多可以被相邻的网格点重用六次。重要的是，从波阵面角度来看，每个`u`条目的加载是一个连续的64个条目块。我们当前的线程块配置（256 × 1 × 1）沿着`x`方向加载四个波的元素，为单个波中的线程提供了缓存和重用`x`方向元素用于`x - 1`和`x + 1`模板计算的机会。然而，沿着`y`和`z`方向加载波必须小心实施以最大化空间局部性并重用相邻波和线程块的值。因此，我们重点通过循环分块优化其中一个方向的加载。 

循环平铺

目前，每个线程仅计算单个网格点的模版。如果让每个线程计算多个网格点的模版会发生什么？这会使每个线程需要更多的加载指令，但由于线程分配的连续性，这些加载更有可能重用缓存的 u 值。这种优化称为循环平铺，将减少启动的线程块数量并增加每个线程计算的模版数量。

在深入优化的拉普拉斯核之前，让我们先通过一个小例子来进一步解释循环平铺的概念和好处：假设我们想执行步幅计算 f[pos] = u[pos - 1] + u[pos] + u[pos + 1]。在不进行平铺的情况下，每个线程需要执行三个加载和一个存储指令。如果我们能让每个线程只加载一个 u 元素并重用已经加载并存储在寄存器中的其他两个 u 元素，这将很理想。换句话说，我们应该尽量减少每个线程的每次存储指令的加载数量。如果我们按某个因子平铺，比如选两个，那么每个线程将执行两次存储，但考虑加载次数：

f[pos]     = u[pos - 1]  + u[pos]     + u[pos + 1];
f[pos + 1] = u[pos]      + u[pos + 1] + u[pos + 2];

注意到 u[pos] 和 u[pos + 1] 出现了两次，这意味着我们只需要加载它们一次。这个观察结果使我们可以重用之前加载的值。为了更清楚地说明这一点，我们可以引入两个变量：

double u0 = u[pos];
double u1 = u[pos + 1];
f[pos]     = u[pos - 1] + u0  + u1;
f[pos + 1] = u0         + u1  + u[pos + 2];

结果，我们显式地将加载指令从 6 次减少到了 4 次（即减少了 33%）。我们现在大约有每次存储 2 次加载。请参阅下表，了解循环平铺因子与加载存储比率的关系。

平铺因子	加载次数	存储次数	比率
1	3	1	3.00
2	4	2	2.00
4	6	4	1.50
8	10	8	1.25
16	18	16	1.13

请记住，增加平铺因子会增加内核中的寄存器使用，这会减少其占用率。如果编译器耗尽寄存器，则寄存器溢出到全局内存，这可能会对性能产生负面影响。

回到 3D 拉普拉斯核，有三个方向可以进行循环平铺。让我们演示 y 方向的平铺。如下面的图示所示，考虑重用模式：

图 1: 每个线程在 xy 平面上循环平铺的说明。加载和重用的网格点数量取决于平铺宽度。

让我们用 m 表示平铺因子，它是编译时确定的用户定义的宏变量。由于此内核的代码修改比较复杂，因此我们将其分为两个阶段：设置和计算。首先，我们应用以下更改：

内核 1 设置（之前）	内核 2 设置（之后）
int j = threadIdx.y + blockIdx.y * blockDim.y; ... // Exit if this thread is on the boundary if (i == 0 || i >= nx - 1 || j == 0 || j >= ny - 1 || k == 0 || k >= nz - 1) return; ... dim3 grid((nx - 1) / block.x + 1, (ny - 1) / block.y + 1, (nz - 1) / block.z + 1);	#define m 1 ... int j = m*(threadIdx.y + blockIdx.y * blockDim.y); ... // Exit if this thread is on the xz boundary if (i == 0 || i >= nx - 1 || k == 0 || k >= nz - 1) return; ... dim3 grid((nx - 1) / block.x + 1, (ny - 1) / (block.y * m) + 1, (nz - 1) / block.z + 1);

内核 2 引入了宏变量 m，用于在 y 方向上划分网格维度。当前设定 m=1 没有平铺效果 - 试验其他值需要重新编译代码。我们去掉了需要在线程退出时与边界重叠的情况。

在我们的下一组代码修改中，我们重点重写平铺因子 m 的 f[pos] = ... 计算。因为我们在 y 方向上平铺，每个线程按 nx 步幅前进。这些修改非常复杂，因此我们将其分为四个步骤：

1. 在主要计算内核中添加一个 for 循环
2. 引入大小为 m 的数组以积累模版点的运行和

3. 将 u 元素加载和 f 元素存储分成单独的循环。
4. 引入一个变量来保存 u 元素以便重用。

下面是应用四个步骤之前和之后的代码片段。我们首先通过在一个for循环中进行我们的模板评估：

Kernel 1计算 - 步骤0（之前）	Step 1（之后）
f[pos] = u[pos] * invhxyz2 + (u[pos - 1] + u[pos + 1]) * invhx2 + (u[pos - nx] + u[pos + nx]) * invhy2 + (u[pos - slice] + u[pos + slice]) * invhz2;	for (int n = 0; n < m; n++) if (j + n > 0 && j + n < ny - 1) f[pos + n*nx] = u[pos + n*nx] * invhxyz2 + (u[pos - 1 + n*nx] + u[pos + 1 + n*nx]) * invhx2 + (u[pos - nx + n*nx] + u[pos + nx + n*nx]) * invhy2 + (u[pos - slice + n*nx] + u[pos + slice + n*nx]) * invhz2;

在编译时已知的值`m`的for循环引入具有与循环展开相似的效果，其中编译器将最小化循环的开销。请记住，`m`不能太大，否则编译器会将寄存器溢出到全局内存。

此时，内核在技术上已被平铺，然而，为了最小化负载-存储比，还需要进行一些代码修改。接下来我们创建一个累加数组：

Step 1（之前）	Step 2（之后）
for (int n = 0; n < m; n++) if (j + n > 0 && j + n < ny - 1) f[pos + n*nx] = u[pos + n*nx] * invhxyz2 + (u[pos - 1 + n*nx] + u[pos + 1 + n*nx]) * invhx2 + (u[pos - nx + n*nx] + u[pos + nx + n*nx]) * invhy2 + (u[pos - slice + n*nx] + u[pos + slice + n*nx]) * invhz2;	T Lu[m] = {0}; for (int n = 0; n < m; n++) if (j + n > 0 && j + n < ny - 1) { Lu[n] = u[pos + n*nx] * invhxyz2; Lu[n] += (u[pos - 1 + n*nx] + u[pos + 1 + n*nx]) * invhx2; Lu[n] += (u[pos - nx + n*nx] + u[pos + nx + n*nx]) * invhy2; Lu[n] += (u[pos - slice + n*nx] + u[pos + slice + n*nx]) * invhz2; f[pos + n*nx] = Lu[n]; }

累加数组`Lu`暂时保存计算的总和。需要注意的是，我们保持了模板计算的原始顺序——首先加载x方向的模板，然后是y方向的模板，最后是z方向的模板。我们最终将重新审视这个顺序，但接下来的步骤是将加载和存储步骤分开：

Step 2（之前）	Step 3（之后）
T Lu[m] = {0}; for (int n = 0; n < m; n++) if (j + n > 0 && j + n < ny - 1) { Lu[n] = u[pos + n*nx] * invhxyz2; Lu[n] += (u[pos - 1 + n*nx] + u[pos + 1 + n*nx]) * invhx2; Lu[n] += (u[pos - nx + n*nx] + u[pos + nx + n*nx]) * invhy2; Lu[n] += (u[pos - slice + n*nx] + u[pos + slice + n*nx]) * invhz2; f[pos + n*nx] = Lu[n]; }	T Lu[m] = {0}; for (int n = 0; n < m; n++) if (j + n > 0 && j + n < ny - 1) { Lu[n] = u[pos + n*nx] * invhxyz2; Lu[n] += (u[pos - 1 + n*nx] + u[pos + 1 + n*nx]) * invhx2; Lu[n] += (u[pos - nx + n*nx] + u[pos + nx + n*nx]) * invhy2; Lu[n] += (u[pos - slice + n*nx] + u[pos + slice + n*nx]) * invhz2; } for (int n = 0; n < m; n++) if (j + n > 0 && j + n < ny - 1) f[pos + n*nx] = Lu[n];

将加载和存储分成独立的for循环使所有`Lu`元素能够在写入`f`之前同时累积模板计算。第四个也是最后一个改变则是通过重用加载的`u`元素跨不同模板来显式地删除加载指令。每次迭代的`n`仍然会加载`x`方向和`z`方向的模板以计算`Lu[n]`，但现在可以潜在地重用`u`元素来计算属于`Lu[n-1]`和/或`Lu[n+1]`的`y`方向模板：

Step 3（之前）

Kernel 2计算 - 步骤4（之后）

T Lu[m] = {0}; for (int n = 0; n < m; n++) if (j + n > 0 && j + n < ny - 1) { Lu[n] = u[pos + n*nx] * invhxyz2; Lu[n] += (u[pos - 1 + n*nx] + u[pos + 1 + n*nx]) * invhx2; Lu[n] += (u[pos - nx + n*nx] + u[pos + nx + n*nx]) * invhy2; Lu[n] += (u[pos - slice + n*nx] + u[pos + slice + n*nx]) * invhz2; } for (int n = 0; n < m; n++) if (j + n > 0 && j + n < ny - 1) f[pos + n*nx] = Lu[n];

T center; T Lu[m] = {0}; for (int n = 0; n < m; n++) { center = u[pos + n*nx]; Lu[n] = center * invhxyz2 + (u[pos - 1 + n*nx] + u[pos + 1 + n*nx]) * invhx2; if (n == 0) Lu[n] += u[pos - nx + n*nx] * invhy2; if (n > 0) Lu[n-1] += center * invhy2; if (n < m - 1) Lu[n+1] += center * invhy2; if (n == m - 1) Lu[n] += u[pos + nx + n*nx] * invhy2; Lu[n] += (u[pos - slice + n*nx] + u[pos + slice + n*nx]) * invhz2; } for (int n = 0; n < m; n++) if (j + n > 0 && j + n < ny - 1) f[pos + n*nx] = Lu[n];

下面是捕获所有上述代码修改的完整kernel 2实施：

// Tiling factor
#define m 1
template <typename T>
__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;

    // Exit if this thread is on the xz boundary
    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;

    // Each thread accumulates m stencils in the y direction
    T Lu[m] = {0};

    // Scalar for reusable data
    T center;

    // Loop tiling
    for (int n = 0; n < m; n++) {
        center = u[pos + n*nx]; // store for reuse

        // x direction
        Lu[n] += center *invhxyz2
              + (u[pos - 1 + n*nx] + u[pos + 1 + n*nx]) * invhx2;

        // y - 1, first n
        if (n == 0) Lu[n] += u[pos - nx + n*nx] * invhy2;

        // 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, last n
        if (n == m - 1) Lu[n] += u[pos + nx + n*nx] * invhy2;

        // z - 1 and z + 1
        Lu[n] += (u[pos - slice + n*nx] + u[pos + slice + n*nx]) * invhz2;
    }

    // Store only if thread is inside y boundary
    for (int n = 0; n < m; n++)
      if (j + n > 0 && j + n < ny - 1)
        f[pos + n*nx] = Lu[n];
 }

template <typename T>
void laplacian(T *d_f, T *d_u, int nx, int ny, int nz, int BLK_X, int BLK_Y, int BLK_Z, T hx, T hy, T hz) {

    dim3 block(BLK_X, BLK_Y, BLK_Z);
    dim3 grid((nx - 1) / block.x + 1, (ny - 1) / (block.y * m) + 1, (nz - 1) / block.z + 1);
    T invhx2 = (T)1./hx/hx;
    T invhy2 = (T)1./hy/hy;
    T invhz2 = (T)1./hz/hz;
    T invhxyz2 = -2. * (invhx2 + invhy2 + invhz2);

    laplacian_kernel<<<grid, block>>>(d_f, d_u, nx, ny, nz, invhx2, invhy2, invhz2, invhxyz2);
}

需要注意的是，这个内核目前写成了`ny`必须为`block.y * m`的整数倍。让我们实验一些与我们选择的问题大小兼容的`m`值，看看循环平铺是否有任何好处。

	Speedup	% of target
Kernel 1 - Baseline	1.00	69.4%
Kernel 2 - Loop tiling m=1	1.00	69.4%
Kernel 2 - Loop tiling m=2	0.98	68.3%
Kernel 2 - Loop tiling m=4	0.94	65.5%
Kernel 2 - Loop tiling m=8	0.92	64.0%
Kernel 2 - Loop tiling m=16	0.29	20.1%

令人奇怪的是，所测试的 m 都没有带来显著的加速。事实上，增加平铺因子反而加剧了性能问题。让我们检查一下 FETCH_SIZE 和 L2CacheHit 指标以获得进一步的见解：

	FETCH_SIZE (GB)	Fetch efficiency (%)	L2CacheHit (%)
Theoretical	1.074	-	-
Kernel 1 - Baseline	2.014	53.3	65.0
Kernel 2 - Loop tiling m=1	2.014	53.3	65.0
Kernel 2 - Loop tiling m=2	1.848	58.1	60.5
Kernel 2 - Loop tiling m=4	1.880	57.1	57.0
Kernel 2 - Loop tiling m=8	1.820	59.0	56.0
Kernel 2 - Loop tiling m=16	5.637	19.1	40.9

从表中可以看出，取回效率只是略有提高，而 L2 缓存命中率显著下降。我们怀疑原因是累加步骤遵循了与初始内核相同的访问模式。也就是说，我们首先计算 x 方向的模板，然后计算可能可重用的 y - 1 和 y + 1 模板，最后是 z - 1 和 z + 1 模板。从内存地址的角度来看，读访问模式前后跳跃，这可能产生了一些问题，如下所述并解决。

重新排序读取访问模式

优化后，没有显著的加速效果。虽然增加块因子减少了加载与存储的比率，但这不一定会直接减少全局数据的传输量。要减少`FETCH_SIZE`，需减少L2缓存和全局内存之间的数据传输。随着加载与存储比率的增加，发送到L1缓存的读取请求数量应减少，同样适用于L2（假设L1缓存命中率不变）。由于我们观察到`FETCH_SIZE`保持不变并且`L2CacheHit`降低，说明优化减少了对L2缓存的压力（发送给它的请求减少），但未能改善从全局内存加载到L2缓存的数据重用。为了理解之前的内核为何未能实现L2数据重用的最优效果，让我们在`m=2`的情况下可视化3D模板及其读取访问模式：

图2：块因子为`m=2`时的三维有限差分模板。黑色数字代表Kernel 2中每个线程访问`u`元素的顺序。

图3：Kernel 2中单个线程对数组`u`的内存访问模式，块因子为`m=2`。数字和黑色箭头对应线程访问`u`元素的顺序。`n=0`和`n=1`行表示计算模板时所需的`u`元素。第一个访问的元素（`u[pos]`）在`n=0`迭代中加载，并在`n=1`迭代中重用作为`y-1`元素。同样，第七次访问的元素（`u[pos+nx]`）在`n=1`迭代中加载，并在`n=0`迭代中重用作为`y+1`元素。

我们马上看到一个问题。线程经常需要在`u`数组的内存空间内“向后”跳转。在访问每个网格点的`z+1`元素后，线程需要“向后”跳转以访问下一个`n`迭代的`x`方向元素。在内存地址中频繁地来回访问`u`元素可能会提前将可重用数据从缓存中逐出。我们更希望重新排列内核中的指令，只使用一个方向，即按升序访问内存地址：

图4：块因子为`m=2`的三维有限差分模板。黑色数字代表提议的内核中每个线程访问`u`元素的顺序。

图5：提议中的单个线程对数组`u`的内存访问模式，块因子为`m=2`。数字和黑色箭头对应线程访问`u`元素的顺序。`n=0`和`n=1`行表示计算模板时所需的`u`元素。第5次访问的元素（`u[pos]`）在`n=0`迭代中加载，并在`n=1`迭代中重用作为`y-1`元素。同样，第8次访问的元素（`u[pos + nx]`）在`n=1`迭代中加载，并在`n=0`迭代中重用作为`y+1`元素。

在这种新方法下，我们首先访问所有`z-1`元素，然后是`n=0`迭代的一个`y-1`元素，接着是所有`x`方向的元素，`n=m-1`迭代的一个`y+1`元素，最后访问所有`z+1`元素。现在，每个线程按升序内存地址访问所有所需的`u`元素。这需要内核的大幅重写，因此我们首先呈现完整实现：

// 块因子
#define m 1
template <typename T>
__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; // 边界检查

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

        // x
        center = u[pos + n*nx]; // 存储以供重用
        Lu[n] += center * invhxyz2;

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

        // 重用: 前一个n的y+1
        if (n > 0) Lu[n-1] += center * invhy2;

        // 重用: 下一个n的y-1
        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)
        f[pos + n*nx] = Lu[n];
}

template <typename T>
void laplacian(T *d_f, T *d_u, int nx, int ny, int nz, int BLK_X, int BLK_Y, int BLK_Z, T hx, T hy, T hz) {

    dim3 block(BLK_X, BLK_Y, BLK_Z);
    dim3 grid((nx - 1) / block.x + 1, (ny - 1) / (block.y * m) + 1, (nz - 1) / block.z + 1);
    T invhx2 = (T)1./hx/hx;
    T invhy2 = (T)1./hy/hy;
    T invhz2 = (T)1./hz/hz;
    T invhxyz2 = -2. * (invhx2 + invhy2 + invhz2);

    laplacian_kernel<<<grid, block>>>(d_f, d_u, nx, ny, nz, invhx2, invhy2, invhz2, invhxyz2);
} 

接下来让我们深入探讨内核中的计算步骤。首先，我们访问所有`z - 1`网格点，然后是一个`y - 1`：

    // 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

请注意，引入了条件运算符，以确保只有在`n = 0`网格点不限于`y`边界时才计算`y - 1`模板。`z - 1`和`y - 1`元素都不会在线程级别重用。

接下来，线程计算`x`方向上的模板：

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

        // x
        center = u[pos + n*nx]; // 存储以供重用
        Lu[n] += center * invhxyz2;

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

        // 重用: 前一个n的y + 1
        if (n > 0) Lu[n-1] += center * invhy2;

        // 重用: 下一个n的y - 1
        if (n < m - 1) Lu[n+1] += center * invhy2;
    }

同样，`x - 1`和`x + 1`点不会在线程级别重用，但中心元素`u[pos + n*nx]`最多可以重用两次，就像之前的内核一样。

然后，我们加载最终的`y + 1`点和所有`z + 1`点：

    // 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;

再次使用了条件运算符，以确保仅在`n = m - 1`网格点不位于`y`边界时才计算`y + 1`模板。

最后，仅在线程位于`y`边界内时才将结果写回内存：

    // 只有线程在y边界内时才存储结果
    for (int n = 0; n < m; n++)
      if (n + j > 0 && n + j < ny - 1)
        f[pos + n*nx] = Lu[n];

让我们现在使用相同的块因子进行实验，看看这种重新排序是否产生影响：

	速度提升	目标百分比
Kernel 1 - 基线	1.00	69.4%
Kernel 2 - 循环分块 m=1	1.00	69.4%
Kernel 2 - 循环分块 m=2	0.98	68.3%
Kernel 2 - 循环分块 m=4	0.94	65.5%
Kernel 2 - 循环分块 m=8	0.92	64.0%
Kernel 2 - 循环分块 m=16	0.29	20.1%
Kernel 3 - 重新排序加载 m=1	1.20	82.9%
Kernel 3 - 重新排序加载 m=2	1.28	88.9%
Kernel 3 - 重新排序加载 m=4	1.34	93.1%
Kernel 3 - 重新排序加载 m=8	1.37	94.8%
Kernel 3 - 重新排序加载 m=16	0.42	29.4%

即使在`m=1`的情况下，重新排序`u`元素的访问模式已经显著提升了性能。每个`m`的增量提升符合预期。让我们查看新内核的`rocprof`指标：

	FETCH_SIZE (GB)	提取效率 (%)	L2CacheHit (%)
理论值	1.074	-	-
Kernel 1 - 基线	2.014	53.3	65.0
Kernel 2 - 循环分块 m=1	2.014	53.3	65.0
Kernel 2 - 循环分块 m=2	1.848	58.1	60.5
Kernel 2 -循环分块 m=4	1.880	57.1	57.0
Kernel 2 - 循环分块 m=8	1.820	59.0	56.0
Kernel 2 - 循环分块 m=16	5.637	19.1	40.9
Kernel 3 - 重新排序加载 m=1	1.347	79.7	72.0
Kernel 3 - 重新排序加载 m=2	1.166	92.1	70.6
Kernel 3 - 重新排序加载 m=4	1.107	97.0	68.8
Kernel 3 - 重新排序加载 m=8	1.080	99.4	67.7
Kernel 3 - 重新排序加载 m=16	3.915	27.4	44.5

FETCH_SIZE这一指标显著减少，使我们接近理论极限。`L2CacheHit`命中率不仅提高了，而且超出了我们从基线内核原本得到的值。不过，当`m=16`时，我们观察到了缓存命中率的大幅下降以及取数大小的显著增加。对于所选问题，`m=8`的内核3是我们目前最好的内核，达到了目标有效内存带宽的近95%和超过99%的取数效率。

总结

结合两种优化后，`FETCH_SIZE`减少了近2倍。这表明我们的HIP内核可以为特定网格大小有效地加载数据。要实现这一点，我们首先通过循环分块显式计算多个模板，减少了每个存储指令的加载次数。然而，最初的实现并未提高性能。为了解决这一问题，我们重新排序了内存访问模式，以提高L2缓存命中率。现在的问题是我们是否“完成”了针对Laplacian的有限差分方法的初始HIP实现的优化。我们必须先解决以下几个悬而未决的问题：

1. 是否还有进一步提高性能的空间？我们已经优化了L2缓存与全局内存之间的数据移动，因此我们必须在其他领域寻找性能提升的机会，例如隐藏延迟。
2. 为什么`m=16`的性能显著下降？无论是否重新排序内存访问都会发生这种情况。或许解决潜在问题可以帮助我们更接近目标？
3. 其他架构和问题规模如何影响块因子的选择？到目前为止，我们的所有优化都针对单个MI250X GCD和问题规模为`nx,ny,nz = 512, 512, 512`。

本系列的下一篇文章将回答这些悬而未决的问题。

附带代码示例

如果您有任何问题或意见，请在GitHub上联系我们 讨论区

---

[1](1,2)

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