有限差分法是计算物理学中一种典型的网格离散化方法，广泛应用于从地球物理学（天气和石油天然气）和电磁学（半导体和天体物理学）到气体动力学（气流和等离子体）的各种应用中。网格计算代码的特点是需要访问网格点（模板）的局部邻域，以便计算单个网格点的值，这意味着算法的性能与其内存访问模式密切相关。GPU 由于其算法的规则结构、高度暴露的并行性和对 GPU 高内存带宽的有效利用，已被证明非常适合加速网格计算代码。在本篇博文中，我们将使用 AMD 的 Heterogeneous Interface for Portability (HIP) API 开发一个 GPU 加速的网格计算代码，这将让我们能够显式控制内存访问模式。在本次开发过程中，我们将讨论以下概念：
-
使用有限差分法离散化拉普拉斯算子
-
展示离散化的直接 HIP 实现
-
指导读者选择优化策略来提高网格计算代码的性能
-
演示如何使用
rocprof获取内核的性能指标
拉普拉斯算子
拉普拉斯算子 是一种微分算子，可以使用偏微分方程描述各种物理现象。例如，它可用于描述声波传播（波动方程）、热流（热方程）以及静电和引力势（泊松方程）。它在计算机视觉和网格生成领域也有应用。有关更多信息，请参阅 wikipedia。
作为可由拉普拉斯算子描述的物理现象的一个例子，下面的动画是热方程的二维数值模拟，显示了热量从热区域（凸起的字母）向冷区域（平面）的传递。
| 初始条件 | 模拟 |
|---|---|
|  |
在笛卡尔坐标系中，拉普拉斯算子表现为标量场 的梯度的散度
等距的网格点，并在每个空间方向上应用二阶有限差分近似。例如， 方向的二阶导数近似为：其中 是 x 方向上的网格间距（相邻网格点在 x 方向上的距离）。将其扩展到三维需要对所有空间方向应用此模板：
下面的代码示例演示了如何在主机代码中实现拉普拉斯算子模具
template <typename T>void laplacian_host(T *h_f, const T *h_u, T hx, T hy, T, hz, int nx, int ny,int nz) {#define u(i, j, k) h_u[(i) + (j) * nx + (k) * nx * ny]#define f(i, j, k) h_f[(i) + (j) * nx + (k) * nx * ny] // Skip boundary points for (k = 1; k < nz - 1; ++k) for (j = 1; j < ny - 1; ++j) for (i = 1; i < nx - 1; ++i) { // Apply stencil approximation of Laplacian to all interior points u_xx = ( u(i + 1, j, k) - 2 * u(i, j, k) + u(i - 1, j, k) ) / (hx * hx); u_yy = ( u(i, j + 1, k) - 2 * u(i, j, k) + u(i, j - 1, k) ) / (hy * hy); u_zz = ( u(i, j, k + 1) - 2 * u(i, j, k) + u(i, j, k - 1) ) / (hz * hz); f(i, j, k) = u_xx + u_yy + u_zz; }#undef u#undef f}主机代码使用 C 预处理器宏来更轻松地将数学表示式 转换为代码：u(i, j, k) – 您将在稍后详细了解此转换。函数 laplacian_host 接收输入 h_u[] 和输出 h_f[]，并期望它们每个至少包含 nx * ny * nz 个元素（每个方向上的网格点数）。为防止因越界错误导致内存访问冲突，循环边界排除了每个网格方向上的第一个和最后一个元素。最内层的循环将拉普拉斯算子模具应用于 h_u[] 并将结果存储在 h_f[] 中。为区分主机数组和设备数组，我们分配前缀 h_（主机）和 d_（设备）。当我们介绍设备代码时，这种区分将变得重要。
数据布局
三维数据布局方式是，i 方向上的网格点在内存中是连续的，而 k 方向上的网格点则以 nx * ny 的步长进行跨越。这种映射可以通过宏来体现
#define u(i, j, k) h_u[(i) + (j) * nx + (k) * nx * ny]i、j 和 k 周围的括号可确保正确展开诸如 u(i, j + 1, k) 等表达式，依此类推。
下图显示了应用于网格点 pos 处的拉普拉斯算子模具
图 1：三维空间中的有限差分模具。
由于跨步数据布局，相邻 j 和 k 方向上的网格点在内存中实际上相距甚远。例如，考虑 nx = ny = nz = 5 的情况
图 2：中心差分模具的跨步内存空间。在此我们看到 h_u[] 数组（红色）中的哪些条目需要更新输出数组 h_f[]。请注意，这些条目不是连续的。
此博客文章的目标之一是展示如何在 AMD GPU 上高效处理这种内存访问模式。
HIP 实现
首先，考虑一个大型立方体，其中 nx = ny = nz = 512。我们首先使用 double 精度为输入数组 d_u 和输出数组 d_f 分配内存
// Parametersusing precision = double;size_t nx = 512, ny = 512, nz = 512; // Cube dimensions
// Grid point spacingsprecision hx = 1.0 / (nx - 1), hy = 1.0 / (ny - 1), hz = 1.0 / (nz - 1);
// Input and output arraysprecision *d_u, *d_f;
// Allocate on devicesize_t numbytes = nx * ny * nz * sizeof(precision);hipMalloc((void**)&d_u, numbytes);hipMalloc((void**)&d_f, numbytes);d_u 输入数组使用二次测试函数初始化，该函数简化了验证正确性的任务（为简洁起见省略）。
以下代码片段展示了我们最初的拉普拉斯算子 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);}laplacian 主机函数负责以 BLK_X = 256, BLK_Y = 1, BLK_Z = 1 的线程块大小启动设备内核 laplacian_kernel。当前线程块大小 BLK_X * BLK_Y * BLK_Z = 256 是推荐的默认选择，因为它与硬件工作调度程序很好地映射，但其他选择可能会表现更好。总之，选择此参数以获得高性能可能高度依赖于问题规模、内核实现以及设备的硬件特性。我们将在未来的文章中更详细地讨论此主题。尽管此内核实现支持三维块大小，但 GPU 没有二维和三维的概念。定义三维组的可能性是作为一种方便的功能提供给程序员的。
HIP 实现适用于适合 GPU 内存的任何问题规模。对于不能被线程块大小整除的问题规模，我们会引入额外的组来处理余数。如果余数很小，这些额外的线程块的线程利用率可能会很低，即许多线程映射到计算域之外的区域，而那里没有任何工作要做。这种低线程利用率会负面影响性能。在本篇博客文章中，我们将重点关注问题规模能够被线程块大小整除的情况。要调度的块数由向上取整操作 (nx - 1) / block.x + 1, … 确定。如前所述，它会向上取整以确保覆盖整个计算域。
注意内核中的退出条件。当线程在边界上或边界外时，我们返回（退出），因为我们只计算内部点的有限差分。内核退出策略的实际工作方式是，由于波前以锁步方式执行指令，任何在边界上或边界外的线程在内核的其余部分执行期间都将被屏蔽。被屏蔽的线程在指令执行期间将处于空闲状态。让我们回到块大小不能整除问题规模的性能问题。如果我们选择，例如 nx=258，我们将在 x 方向上得到两个线程块，但第二个线程块中只有一个线程有工作要做。
主要的计算，
// 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;多次访问数组 u。对于每次访问，编译器将生成一个全局加载指令，该指令将访问全局内存中的数据，除非该数据已存在于缓存中。正如我们将看到的，这些全局加载指令可能对性能产生重大影响。
计算密集型 vs. 内存密集型
为了分析我们 HIP 内核的性能，了解它是计算密集型还是内存密集型很有帮助。计算密集型内核受计算单元吞吐量的限制，而内存密集型内核受内存带宽或内存访问延迟的限制。要对我们的内核进行分类，我们使用以下公式来研究其算术强度 (AI)
其中 指浮点运算的总数， 指主内存流量。如果 ，则称该算法是内存密集型的，其中 表示峰值浮点运算率， 表示峰值内存带宽。反之，如果 ，则称该算法是计算密集型的。请注意，我们忽略整数运算，只关注浮点运算。这种对浮点运算和内存流量的分析通常用于创建“天花板模型”。对于有兴趣了解更多关于天花板模型的读者，我们推荐参考 维基百科链接 及其附带的参考文献。
在上面的 laplacian_kernel 中，每个线程执行 10 次浮点运算，跨越 7 个 d_u 元素，并将计算结果存储到单个 d_f 元素中。但是，我们将假设理想条件，即无限缓存，这意味着每个线程获取一个 d_u 元素，并且该元素被其他线程重用于模具计算。每个线程的 和 分别等于 10 和 2 * sizeof(double) = 16，因此我们的 AI 为 0.63。像这样的模具代码通常具有较小的 与 比率，并且当我们取消无限缓存的假设时，这个比率会更小。
在本篇博客文章的剩余部分，让我们考虑来自 MI250X GPU 的单个图形计算卡 (GCD)，其中 = 24 TFLOPs/sec，高带宽内存 (HBM) 带宽 = 1.6 TBytes/sec。 比率是近似 AI 的 24 倍，这清楚地将我们的 HIP 内核归类为内存密集型，而不是计算密集型。因此，我们的优化将侧重于解决内存瓶颈。
性能
为了进行性能评估，我们可以使用有效内存带宽作为我们的衡量指标 (FOM)。粗略地说，有效内存带宽是数据传输的平均速率，考虑需要通过整个内存子系统传输的最少数据量：从高延迟的片外 HBM 到低延迟的寄存器。有效内存带宽定义为
effective_memory_bandwidth = (theoretical_fetch_size + theoretical_write_size) / average_kernel_execution_time;为了测量平均内核执行时间，我们使用函数 hipEventElapsedTime()（有关完整的实现细节，请参阅配套的 laplacian.cpp）。用于计算 nx * ny * nz 立方体的理论获取和写入大小（以字节为单位）的公式为
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);由于模具的形状，理论获取大小包括读取 d_u 中的所有网格点，除了立方体的 8 个角点和 12 个边点。在理论写入大小中，只有 d_f 中的内部网格点才会被写回内存。这些计算忽略了内存访问粒度，并且再次假设无限缓存。因此，它们为实际数据移动设置了下限。
以下是单个 MI250X GCD 上的性能数据
$ ./laplacian_dp_kernel1Kernel: 1Precision: doublenx,ny,nz = 512, 512, 512block sizes = 256, 1, 1Laplacian kernel took: 2.64172 ms, effective memory bandwidth: 808.148 GB/s对于经过优化的内核和能够使设备饱和的足够大的工作量，有效内存带宽应尽可能接近峰值理论 HBM 带宽。完成此任务需要满足两个目标：
-
减少到 HBM 的过度内存流量，使其尽可能接近理论极限，
-
使设备的 HBM 带宽饱和。
鉴于 MI250X 每个 GCD 的理论带宽为 1638.4 GB/s，基线内核实现了峰值理论极限的约 50 %[1]。是什么瓶颈阻碍了进一步提高性能？是由于过多的数据移动还是内存带宽饱和度不佳？例如，如果计算模具时数据重用性很差，可能会出现过多的数据移动。也许一些相邻的值被加载到缓存中，但被其他值替换，导致缓存未命中和额外的访问主内存。当内核生命周期内没有足够的内存请求在执行中时，就会发生内存带宽饱和度不足。一般来说，这可能有很多原因。也许每个波前请求的数据量很少，并且由于资源限制，没有足够的波前可以在每个 CU 上并发运行。另一种可能性是一些波前因等待数据依赖而停滞，并且不eligible发出内存指令。为了理解瓶颈是什么并解决它，我们需要收集性能计数器数据。
为此，我们使用 AMD ROCm™ 分析器 (rocprof) 来捕获三个重要指标：
FETCH_SIZE : The total kilobytes fetched from the video memory. This is measured with all extra fetches and any cache or memory effects taken into account.WRITE_SIZE : The total kilobytes written to the video memory. This is measured with all extra fetches and any cache or memory effects taken into account.L2CacheHit : The percentage of fetch, write, atomic, and other instructions that hit the data in L2 cache. Value range: 0% (no hit) to 100% (optimal).请注意，FETCH_SIZE 和 WRITE_SUM 的总和可能不完全代表所有内存流量。尽管如此，这些指标肯定可以为 HIP 内核的设计提供更大的见解。可以从 rocprof 查询其他指标（有关选项，请参阅 rocprof —list-basic 和 rocprof —list-derived）。rocprof 允许我们将请求组织在 ASCII 文本文件中。例如，在包含的工作中，我们使用 rocprof_input.txt
pmc: FETCH_SIZEpmc: WRITE_SIZEpmc: L2CacheHitkernel: laplacian_kernel这将捕获列出内核的性能计数器 (pmc) FETCH_SIZE、WRITE_SIZE 和 L2CacheHit。我们可以使用命令对我们的可执行文件运行此命令
$ rocprof -i rocprof_input.txt -o rocprof_output.csv laplacian_dp_kernel这将生成一个 rocprof_output.csv 文件，其中包含每个内核启动的请求指标。有关 rocprof 的更多信息，我们建议所有感兴趣的读者参考 ROCm™ Profiler 文档。
FETCH_SIZE 和 WRITE_SIZE 越接近 theoretical_fetch_size 和 theoretical_write_size，我们的 HIP 内核可能表现越好。L2CacheHit 接近 100% 是内核最优的另一个指标。下表比较了上述三个 rocprof 指标
| FETCH_SIZE (GB) | WRITE_SIZE (GB) | 获取效率 (%) | 写入效率 (%) | L2CacheHit (%) | |
|---|---|---|---|---|---|
| 理论值 | 1.074 | 1.061 | – | – | – |
| 内核 1 | 2.014 | 1.064 | 53.3 | 99.7 | 65.0 |
65% 的 L2CacheHit 率本身并不能提供太多见解。虽然 WRITE_SIZE 指标与其理论估计值相符，但 FETCH_SIZE 却接近翻倍。这一观察告诉我们，约 50% 的模具是重用的。如果我们能将模具重用提高到 100%，FETCH_SIZE 将会减少 0.9 GB 以上，从而使总内存流量减少 31%。假设 HBM 带宽饱和度保持不变，这可能会将内核执行时间缩短 1.44 倍。因此，有效内存带宽的更现实目标将是大约 808.148 GB/s * 1.44 = 1165 GB/s，这大约是每个 GCD 的峰值理论 HBM 带宽 1638.4 GB/s 的 71%[1]。
一般而言，有多种不同的方法可以优化有限差分内核。在下一篇文章中，我们将解释和讨论一些可以应用的优化，以减少过多的内存流量并达到新的有效内存带宽目标。
结论
至此，我们完成了使用 HIP 开发和优化拉普拉斯算子内核的第一部分。在开始任何优化工作之前，收集性能测量以建立基线并识别性能瓶颈非常重要。在这项工作中，我们确定基线内核由于其算术强度而受内存带宽限制。我们还发现，与我们的估计相比，内核从全局内存空间加载了显著更多的数据。这些条件告诉我们，减少过多的内存负载将对性能产生影响。在本系列的第二部分中，我们将介绍旨在减少全局内存数据移动的优化技术，这些技术不仅对模具内核有价值，而且对于使用规则的、非连续内存访问模式的其他带宽受限内核也很有价值。
如果您有任何问题或评论，请在 GitHub 讨论区 与我们联系
测试使用 ROCm 5.3.0-63 版本进行。基准测试结果不是经过验证的性能数据，仅用于展示代码修改的相对性能改进。实际性能结果取决于多种因素，包括系统配置和环境设置，不保证结果的可重现性。
