前言
第七章学了卷积，这章学模板（Stencil）。你可能会问：这俩不是差不多吗？确实很像，但有微妙区别。卷积 强调的是信号处理中的滑动加权和，权重来自滤波器；模板 强调的是科学计算中的邻域更新，比如 PDE 求解、物理仿真。从 CUDA 优化角度看，技术是相通的，但模板计算有自己的特点：通常是多次迭代 的，需要处理时间步进 和数据交换 。
📦 配套资源 ：本系列文章配有完整的 GitHub 仓库 ，包含每章的练习题解答、CUDA 代码实现和详细注释。所有代码都经过测试，可以直接运行。
模板计算基础
什么是模板计算
模板计算（Stencil Computation）：用固定的邻域模式 更新网格中的每个点。
一维热传导 （最简单的例子）：
1 2 时间 t: [..., A, B, C, ...] 时间 t+1: B' = (A + B + C) / 3
每个点的新值是自己和邻居的加权平均。这就是模板：一个描述"如何从邻居计算自己"的模式。
模板的形状
常见的模板形状：
1D 三点模板 ：
1 2 3 [ -1 0 +1 ] ↓ ↓ ↓ A B C
2D 五点模板（冯·诺依曼邻域） ：
1 2 3 4 5 [0,-1] ↑ [-1,0] ← [0,0] → [+1,0] ↓ [0,+1]
2D 九点模板（摩尔邻域） ：
1 2 3 [-1,-1] [0,-1] [+1,-1] [-1, 0] [0, 0] [+1, 0] [-1,+1] [0,+1] [+1,+1]
3D 七点模板 ：
模板 vs 卷积
特性
卷积
模板
权重
来自滤波器，任意值
通常是固定系数
迭代
通常单次
多次时间步进
边界
零填充常见
周期/固定边界更常见
应用
信号/图像处理
PDE/物理仿真
读写模式
只读输入，写输出
读旧值，写新值
从 GPU 优化角度，两者的核心技术是相同的：常量内存 存系数，共享内存 Tiling 减少全局内存访问。
2D 模板：热传导方程
物理背景
二维热传导方程：
∂ T ∂ t = α ( ∂ 2 T ∂ x 2 + ∂ 2 T ∂ y 2 ) \frac{\partial T}{\partial t} = \alpha \left( \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} \right)
∂ t ∂ T = α ( ∂ x 2 ∂ 2 T + ∂ y 2 ∂ 2 T )
用有限差分离散化后：
T i , j n + 1 = T i , j n + α Δ t ( T i + 1 , j n + T i − 1 , j n + T i , j + 1 n + T i , j − 1 n − 4 T i , j n Δ x 2 ) T_{i,j}^{n+1} = T_{i,j}^n + \alpha \Delta t \left( \frac{T_{i+1,j}^n + T_{i-1,j}^n + T_{i,j+1}^n + T_{i,j-1}^n - 4T_{i,j}^n}{\Delta x^2} \right)
T i , j n + 1 = T i , j n + α Δ t ( Δ x 2 T i + 1 , j n + T i − 1 , j n + T i , j + 1 n + T i , j − 1 n − 4 T i , j n )
简化为五点模板：
1 T_new[i][j] = c0 * T[i][j] + c1 * (T[i-1][j] + T[i+1][j] + T[i][j-1] + T[i][j+1])
其中 c0 = 1 - 4*alpha*dt/dx²，c1 = alpha*dt/dx²。
朴素实现
1 2 3 4 5 6 7 8 9 10 __global__ void stencil_2d_basic(float *in, float *out, int N) { int i = blockIdx.x * blockDim.x + threadIdx.x; int j = blockIdx.y * blockDim.y + threadIdx.y; if (i > 0 && i < N-1 && j > 0 && j < N-1) { out[i*N + j] = C0 * in[i*N + j] + C1 * (in[(i-1)*N + j] + in[(i+1)*N + j] + in[i*N + (j-1)] + in[i*N + (j+1)]); } }
问题 ：每个点读取 5 次全局内存，相邻点的访问有大量重叠。
数据复用分析
考虑一行线程计算：
1 2 3 线程 j: 读取 T[i][j-1], T[i][j], T[i][j+1], T[i-1][j], T[i+1][j] 线程 j+1: 读取 T[i][j], T[i][j+1], T[i][j+2], T[i-1][j+1], T[i+1][j+1] ↑ ↑ 共享
相邻线程共享了 T[i][j] 和 T[i][j+1]。如果把一个 Tile 的数据加载到共享内存，这些共享访问都变成快速的片上访问。
Tiled 模板实现
输入输出 Tile 关系
与卷积类似，计算 OUT_TILE_SIZE × OUT_TILE_SIZE 的输出，需要 (OUT_TILE_SIZE + 2*r) × (OUT_TILE_SIZE + 2*r) 的输入。
对于五点模板（r=1）：
1 2 输入 Tile：(T + 2) × (T + 2) 输出 Tile：T × T
线程与数据映射
有两种策略：
策略 1：线程数 = 输入 Tile 大小
1 2 3 Block: (T+2) × (T+2) 线程 每线程加载 1 个输入元素 只有内部 T × T 线程计算输出
策略 2：线程数 = 输出 Tile 大小
1 2 3 Block: T × T 线程 部分线程负责加载 Halo 元素 所有线程都计算输出
书中推荐策略 1 ，因为加载更简单，且边界线程虽不计算但仍能并行执行。
完整实现
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 #define IN_TILE_SIZE 18 // 16 + 2 #define OUT_TILE_SIZE 16 __global__ void stencil_2d_tiled(float *in, float *out, int N) { __shared__ float tile[IN_TILE_SIZE][IN_TILE_SIZE]; // 线程在输入 Tile 中的位置 int tx = threadIdx.x, ty = threadIdx.y; // 对应的全局输入位置 int in_row = blockIdx.y * OUT_TILE_SIZE + ty - 1; // -1 是因为 Halo int in_col = blockIdx.x * OUT_TILE_SIZE + tx - 1; // 加载（含边界检查） if (in_row >= 0 && in_row < N && in_col >= 0 && in_col < N) { tile[ty][tx] = in[in_row * N + in_col]; } else { tile[ty][tx] = 0.0f; // 边界外填 0 } __syncthreads(); // 只有内部线程计算输出 if (tx > 0 && tx < IN_TILE_SIZE - 1 && ty > 0 && ty < IN_TILE_SIZE - 1) { int out_row = blockIdx.y * OUT_TILE_SIZE + ty - 1; int out_col = blockIdx.x * OUT_TILE_SIZE + tx - 1; if (out_row < N && out_col < N) { out[out_row * N + out_col] = C0 * tile[ty][tx] + C1 * (tile[ty-1][tx] + tile[ty+1][tx] + tile[ty][tx-1] + tile[ty][tx+1]); } } }
启动配置
1 2 3 4 5 dim3 block(IN_TILE_SIZE, IN_TILE_SIZE); // 18×18 = 324 线程 dim3 grid((N + OUT_TILE_SIZE - 1) / OUT_TILE_SIZE, (N + OUT_TILE_SIZE - 1) / OUT_TILE_SIZE); stencil_2d_tiled<<<grid, block>>>(d_in, d_out, N);
性能分析
指标
朴素实现
Tiled 实现
每输出元素
5 次全局读取
18²/16² ≈ 1.27 次
加速比
1×
~4×
共享内存
0
18×18×4 = 1.3 KB
线程粗化（Thread Coarsening）
为什么粗化
Tiled 模板有个问题：边界线程只加载不计算，效率不高。
1 2 3 Block 线程数：18×18 = 324 计算线程数：16×16 = 256 利用率：256/324 = 79%
如果模板半径更大，利用率更低。
线程粗化 的思路：让每个线程计算多个输出元素，分摊边界开销。
Z 方向粗化
对于 3D 模板，常用的技巧是在 Z 方向展开：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 #define IN_TILE_SIZE 18 #define OUT_TILE_SIZE 16 #define Z_COARSEN 8 // 每线程处理 8 层 __global__ void stencil_3d_coarsened(float *in, float *out, int Nx, int Ny, int Nz) { __shared__ float tile[3][IN_TILE_SIZE][IN_TILE_SIZE]; // 只需 3 层 int tx = threadIdx.x, ty = threadIdx.y; int ox = blockIdx.x * OUT_TILE_SIZE + tx - 1; int oy = blockIdx.y * OUT_TILE_SIZE + ty - 1; // 预加载前两层 load_layer(tile[0], in, ox, oy, 0, Nx, Ny, Nz); load_layer(tile[1], in, ox, oy, 1, Nx, Ny, Nz); // 滑动窗口遍历 Z 方向 for (int z = 1; z < Nz - 1; z++) { // 加载下一层 load_layer(tile[(z+1) % 3], in, ox, oy, z+1, Nx, Ny, Nz); __syncthreads(); // 计算当前层 if (tx > 0 && tx < IN_TILE_SIZE - 1 && ty > 0 && ty < IN_TILE_SIZE - 1) { int idx = z * Nx * Ny + oy * Nx + ox; out[idx] = compute_stencil(tile, tx, ty, z); } __syncthreads(); } }
寄存器粗化
更激进的做法：把滑动窗口存在寄存器里：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 __global__ void stencil_3d_register(float *in, float *out, int Nx, int Ny, int Nz) { __shared__ float tile[IN_TILE_SIZE][IN_TILE_SIZE]; int tx = threadIdx.x, ty = threadIdx.y; int ox = blockIdx.x * OUT_TILE_SIZE + tx - 1; int oy = blockIdx.y * OUT_TILE_SIZE + ty - 1; // 寄存器存 Z 方向的三个值 float z_prev, z_curr, z_next; z_prev = load_element(in, ox, oy, 0, Nx, Ny, Nz); z_curr = load_element(in, ox, oy, 1, Nx, Ny, Nz); for (int z = 1; z < Nz - 1; z++) { z_next = load_element(in, ox, oy, z+1, Nx, Ny, Nz); // 加载 XY 平面到共享内存 tile[ty][tx] = z_curr; __syncthreads(); // 计算 if (valid_output_thread) { out[idx] = C0 * tile[ty][tx] + C1 * (tile[ty-1][tx] + tile[ty+1][tx] + tile[ty][tx-1] + tile[ty][tx+1]) + C2 * (z_prev + z_next); // Z 方向用寄存器 } // 滑动窗口 z_prev = z_curr; z_curr = z_next; __syncthreads(); } }
优势 ：Z 方向的访问完全在寄存器，无需共享内存。
边界条件处理
常见边界条件
1. 固定边界（Dirichlet）
1 2 3 if (is_boundary) { out[idx] = BOUNDARY_VALUE; // 固定值 }
2. 周期边界（Periodic）
1 2 int wrapped_x = (x + N) % N; // 环绕 int wrapped_y = (y + N) % N;
3. 零梯度边界（Neumann）
1 2 3 // 边界处用相邻内部值 if (x == 0) x = 1; if (x == N-1) x = N-2;
边界处理策略
方法 1：条件分支
1 2 3 4 5 if (i > 0 && i < N-1 && j > 0 && j < N-1) { // 内部点计算 } else { // 边界处理 }
问题：分支发散。
方法 2：Padding
预先在数据周围加一圈边界值，kernel 内部无需判断。
1 2 原始数据 N×N → 填充后 (N+2)×(N+2) 所有计算都在 [1, N] 范围内
方法 3：分离边界 Kernel
1 2 3 4 5 // 先处理内部 stencil_interior<<<blocks, threads>>>(in, out, N); // 再处理边界 stencil_boundary<<<1, boundary_threads>>>(in, out, N);
内部 kernel 无分支，边界 kernel 线程少但逻辑复杂。
时间步进与双缓冲
迭代模式
模板计算通常需要多次迭代：
1 2 3 4 5 for (int t = 0; t < num_steps; t++) { stencil_kernel<<<grid, block>>>(d_in, d_out, N); cudaDeviceSynchronize(); swap(d_in, d_out); // 交换输入输出 }
关键 ：不能原地更新！读取邻居时可能读到已更新的值，导致错误。所以需要双缓冲 ：读旧写新，然后交换。
时间阻塞（Temporal Blocking）
每次 kernel 启动有开销。如果能在一次 kernel 中计算多个时间步，就能减少启动开销。
原理 ：
1 2 3 4 5 输入 Tile 计算 T 个时间步： 时间 0：需要 (OUT + 2r*T) 的输入 时间 1：需要 (OUT + 2r*(T-1)) 的输入 ... 时间 T：输出 OUT 的结果
随着时间推进，有效数据区域向内收缩。
实现复杂度高 ，通常在高性能计算库中使用，手写较难。
3D 模板优化
挑战
3D 模板的共享内存需求更大：
1 2 2D: (T+2)² × 4 = 324 × 4 = 1.3 KB 3D: (T+2)³ × 4 = 8000 × 4 = 32 KB // 可能超限！
解决方案
1. 减小 Tile 大小
1 2 #define TILE_3D 8 // (8+2)³ = 1000 元素，4 KB
代价：降低数据复用率。
2. 2.5D 分解
只加载 XY 平面的 Tile，Z 方向逐层处理：
1 2 3 4 5 6 __shared__ float tile[TILE_Y][TILE_X]; // 只存一层 XY for (int z = 0; z < Nz; z++) { // 加载当前层 // 计算（Z 方向邻居从全局内存读） }
Z 方向的访问仍走全局内存，但 XY 方向利用共享内存。
3. 寄存器滑动窗口
如前所述，Z 方向 3 个值存寄存器，XY 平面存共享内存。
指令级优化
循环展开
对于小模板，手动展开消除循环开销：
1 2 3 4 5 6 7 8 9 10 11 // 展开前 for (int di = -1; di <= 1; di++) { for (int dj = -1; dj <= 1; dj++) { sum += tile[ty+di][tx+dj] * weight[di+1][dj+1]; } } // 展开后 sum = tile[ty-1][tx-1] * W00 + tile[ty-1][tx] * W01 + tile[ty-1][tx+1] * W02 + tile[ty ][tx-1] * W10 + tile[ty ][tx] * W11 + tile[ty ][tx+1] * W12 + tile[ty+1][tx-1] * W20 + tile[ty+1][tx] * W21 + tile[ty+1][tx+1] * W22;
编译器通常会自动展开，但显式展开保证效果。
常量预计算
1 2 3 4 5 6 // 编译时常量 #define C0 0.25f #define C1 0.125f // 或使用 constexpr（C++11） __device__ constexpr float c0 = 0.25f;
避免运行时计算系数。
向量化加载
对于对齐的访问模式，使用向量类型：
1 2 float4 data = *reinterpret_cast<float4*>(&in[idx]); // 一次加载 4 个 float
需要数据对齐（16 字节边界）。
实战：Jacobi 迭代求解器
问题设置
求解泊松方程 ∇²u = f：
1 // Jacobi 迭代：u_new = (u_left + u_right + u_up + u_down - h²f) / 4
CUDA 实现
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 #define BLOCK_SIZE 16 #define TILE_SIZE 18 // BLOCK_SIZE + 2 __global__ void jacobi_iteration(float *u, float *u_new, float *f, int N, float h2) { __shared__ float tile[TILE_SIZE][TILE_SIZE]; int tx = threadIdx.x + 1; // 偏移 1，留出 Halo int ty = threadIdx.y + 1; int gx = blockIdx.x * BLOCK_SIZE + threadIdx.x; int gy = blockIdx.y * BLOCK_SIZE + threadIdx.y; // 加载中心 if (gx < N && gy < N) { tile[ty][tx] = u[gy * N + gx]; } // 加载 Halo if (threadIdx.x == 0 && gx > 0) { tile[ty][0] = u[gy * N + gx - 1]; } if (threadIdx.x == BLOCK_SIZE - 1 && gx < N - 1) { tile[ty][TILE_SIZE - 1] = u[gy * N + gx + 1]; } if (threadIdx.y == 0 && gy > 0) { tile[0][tx] = u[(gy - 1) * N + gx]; } if (threadIdx.y == BLOCK_SIZE - 1 && gy < N - 1) { tile[TILE_SIZE - 1][tx] = u[(gy + 1) * N + gx]; } __syncthreads(); // 计算（跳过边界） if (gx > 0 && gx < N - 1 && gy > 0 && gy < N - 1) { float f_val = f[gy * N + gx]; u_new[gy * N + gx] = 0.25f * (tile[ty-1][tx] + tile[ty+1][tx] + tile[ty][tx-1] + tile[ty][tx+1] - h2 * f_val); } } // Host 端迭代 void jacobi_solve(float *u, float *f, int N, int max_iter, float tol) { float *d_u, *d_u_new, *d_f; // ... 分配内存 ... dim3 block(BLOCK_SIZE, BLOCK_SIZE); dim3 grid((N + BLOCK_SIZE - 1) / BLOCK_SIZE, (N + BLOCK_SIZE - 1) / BLOCK_SIZE); for (int iter = 0; iter < max_iter; iter++) { jacobi_iteration<<<grid, block>>>(d_u, d_u_new, d_f, N, h*h); cudaDeviceSynchronize(); // 交换指针 float *temp = d_u; d_u = d_u_new; d_u_new = temp; // 可选：检查收敛 } }
收敛检查
每隔几步检查残差：
1 2 3 4 5 __global__ void compute_residual(float *u, float *f, float *residual, int N, float h2) { // 计算 ||∇²u - f|| // 使用规约求和 }
不需要每步都检查，开销太大。
性能调优建议
Tile 大小选择
1 2 3 4 5 共享内存限制：96 KB / SM（Ampere） 每 Block 可用：~48 KB 2D Tile: sqrt(48K/4) ≈ 110 → 实际用 32、48、64 3D Tile: cbrt(48K/4) ≈ 23 → 实际用 8、12、16
需要考虑占用率：Tile 太大导致每 SM 的 Block 数减少。
Nsight 指标
关注：
指标
意义
目标
Memory Throughput
带宽利用
接近峰值
Compute Throughput
计算利用
越高越好
Shared Memory Usage
共享内存
不超限
Warp Occupancy
占用率
>50%
常见瓶颈
带宽受限 ：Tiling 优化、寄存器缓存计算受限 ：向量化、指令级并行延迟受限 ：增加占用率、预取
模板优化总结
核心技术
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ┌─────────────────────────────────────────────────────────────┐ │ Level 3: 高级优化 │ │ - 时间阻塞（多步合并） │ │ - 自动调优（搜索最优参数） │ ├─────────────────────────────────────────────────────────────┤ │ Level 2: 数据局部性 │ │ - 共享内存 Tiling │ │ - 寄存器滑动窗口 │ │ - 线程粗化 │ ├─────────────────────────────────────────────────────────────┤ │ Level 1: 基础优化 │ │ - 双缓冲避免竞争 │ │ - 边界条件处理 │ │ - 合并内存访问 │ └─────────────────────────────────────────────────────────────┘
模板 vs 卷积决策
场景
推荐技术
单次滤波
卷积 + 常量内存
迭代求解
模板 + 双缓冲
大半径模板
时间阻塞
3D + 大规模
2.5D + 寄存器滑动
小结
第八章围绕模板计算，在第七章卷积的基础上增加了：
迭代结构 ：模板计算通常是多步迭代的。双缓冲是基本技巧——读旧写新，然后交换。
线程粗化 ：让每个线程计算多个输出，分摊边界线程的开销。Z 方向粗化是 3D 模板的常用技巧。
寄存器优化 ：对于滑动窗口模式，把时间维度或空间维度的小邻域存寄存器，比共享内存更快。
边界处理 ：Dirichlet、Neumann、周期边界各有策略。Padding 预处理能简化 kernel 逻辑，但需要额外内存。
3D 挑战 ：共享内存需求爆炸。2.5D 分解、寄存器滑动是解决方案。
模板计算是科学计算的核心，热传导、流体力学、电磁仿真都用到。掌握这套优化技术，就能高效实现各种 PDE 求解器。下一章进入另一个重要模式——并行直方图。
参考资料：
Hwu, W., Kirk, D., & El Hajj, I. (2022). Programming Massively Parallel Processors: A Hands-on Approach (4th Edition). Morgan Kaufmann.
CUDA C++ Programming Guide - Shared Memory Datta, K., et al. (2008). Stencil computation optimization and auto-tuning on state-of-the-art multicore architectures . SC08.
本文 GitHub 仓库 : https://github.com/psmarter/PMPP-Learning
