内存访问密集型高性能计算
在并行计算和 GPU 优化中，我们通常会遇到两类主要瓶颈：计算限制（Compute-dominated）与访存限制（Memory-bound）。在之前的实验中，我们主要面对的是计算密集型任务。本章我们将视角转向到访存密集型工作负载。以**波动方程模拟（Wave Simulation），**其实不用理解其具体的物理/数学含义，只需要知道我们将模拟波随时间的变化，每个时间点，每个点的值，都依赖于这个点和这个点的邻居在上个时刻的点。可见需要频繁访问非局部数据。
对应的 cpu 代码如下：
template <typename Scene> void wave_cpu_step(float t, float *u0, float const *u1) {
constexpr int32_t n_cells_x = Scene::n_cells_x;
constexpr int32_t n_cells_y = Scene::n_cells_y;
constexpr float c = Scene::c;
constexpr float dx = Scene::dx;
constexpr float dt = Scene::dt;
for (int32_t idx_y = 0; idx_y < n_cells_y; ++idx_y) {
for (int32_t idx_x = 0; idx_x < n_cells_x; ++idx_x) {
int32_t idx = idx_y * n_cells_x + idx_x;
bool is_border =
(idx_x == 0 || idx_x == n_cells_x - 1 || idx_y == 0 ||
idx_y == n_cells_y - 1);
float u_next_val;
if (is_border || Scene::is_wall(idx_x, idx_y)) {
u_next_val = 0.0f;
} else if (Scene::is_source(idx_x, idx_y)) {
u_next_val = Scene::source_value(idx_x, idx_y, t);
} else {
constexpr float coeff = c * c * dt * dt / (dx * dx);
float damping = Scene::damping(idx_x, idx_y);
u_next_val =
((2.0f - damping - 4.0f * coeff) * u1[idx] -
(1.0f - damping) * u0[idx] +
coeff *
(u1[idx - 1] + u1[idx + 1] + u1[idx - n_cells_x] +
u1[idx + n_cells_x]));
}
u0[idx] = u_next_val;
}
}
}
接下来看看上述这类**访存限制（Memory-bound）**问题应该如何在GPU上进行优化。
GPU 的内存层次结构 (Memory Hierarchy)
为了优化内存密集型应用，我们必须深入了解目标硬件的内存架构。以 NVIDIA RTX 4000 Ada 为例，它拥有 4 层内存层次结构：
DRAM（全局内存/显存）：容量 20 GB，带宽约 360 GB/sec。所有的数据通常都从这里出发，但它的速度相对最慢。
L2 缓存 (L2 Cache)：48 MB，带宽大幅提升至 ~2.5 TB/sec。它是所有 SM（流式多处理器）共享的，常规的全局内存读写都会经过这里。
L1 缓存与共享内存 (L1 Cache / Shared Memory)：每个 SM 拥有 128 KB 的 SRAM，带宽极高（全 GPU 聚合带宽可达 13.4 TB/sec）。
寄存器堆 (Register File)：每个 Warp 调度器有 64 KB 的 SRAM，速度最快。
L1 缓存的特殊之处
与常规 CPU 不同，GPU 的 L1 缓存并不保证跨 SM 的一致性（Not Coherent）。因此，绝大多数常规的内存读写实际上会绕过（Bypass） L1 缓存。L1 在 GPU 中主要用于以下四个场景：
缓存 CUDA 的线程局部（Thread-local）内存（如 C 语言堆栈）。
缓存**只读（Read-only）**的全局内存（可以通过给指针加
const restrict或使用__ldg内联指令来触发）。缓存编译器认为未来可能发生变化但仍值得缓存的数据。
作为软件管理的 Scratchpad（即 CUDA 中的 Shared Memory）：这是我们需要在代码中显式分配和管理的一块超高速存储。
各级存储访问延迟的直观感受
在更底层（PTX 虚拟汇编）的视角下，GPU 提供了不同的 Load 指令来控制缓存行为：
ld.global.ca：缓存在所有层级（L1 和 L2）。ld.global.cg：只缓存在 L2 缓存（绕过 L1）。ld.global.cv：不缓存（每次读取都会使得 L2 对应的 cache line 失效并重新拉取）。
可以利用这些命令结合实验感受各级存储访问的延迟差距。
__attribute__((optimize("O0"))) __global__ void l1_mem_latency(
unsigned long *time_start,
unsigned long *time_end,
data_type *array_1,
data_type *array_2) {
unsigned long start_time, end_time;
data_type value1 = 0.0f, result;
unsigned long temp_addr;
asm volatile(
// L1 cache setup - load a zero and create offset address
"ld.global.ca.u64 %0, [%5];\n\t"
"add.u64 %0, %0, %5;\n\t"
// warm-up
"ld.global.ca.f32 %2, [%0];\n\t"
"mov.u64 %1, %%clock64;\n\t"
"ld.global.ca.f32 %2, [%0];\n\t"
"mov.u64 %4, %%clock64;\n\t"
"st.global.f32 [%5], %2;\n\t"
: "=l"(temp_addr), "=l"(start_time), "=f"(result), "+f"(value1), "=l"(end_time)
: "l"(array_2)
: "memory");
*time_start = start_time;
*time_end = end_time;
array_1[0] = result;
}
////////////////////////////////////////////////////////////////////////////////
// L2 Cache Memory Latency
__attribute__((optimize("O0"))) __global__ void l2_mem_latency(
unsigned long *time_start,
unsigned long *time_end,
data_type *array_1,
data_type *array_2) {
unsigned long start_time, end_time;
data_type value1 = 0.0f, result;
unsigned long temp_addr;
asm volatile(
// L2 cache setup - load a zero and create offset address
"ld.global.cg.u64 %0, [%5];\n\t"
"add.u64 %0, %0, %5;\n\t"
"membar.gl;\n\t"
// warm-up（进入 L2）
"ld.global.cg.f32 %2, [%0];\n\t"
"mov.u64 %1, %%clock64;\n\t"
"ld.global.cg.f32 %2, [%0];\n\t"
"mov.u64 %4, %%clock64;\n\t"
"st.global.f32 [%5], %2;\n\t"
: "=l"(temp_addr), "=l"(start_time), "=f"(result), "+f"(value1), "=l"(end_time)
: "l"(array_2)
: "memory");
*time_start = start_time;
*time_end = end_time;
array_1[0] = result;
}
////////////////////////////////////////////////////////////////////////////////
// Global Memory Latency
__attribute__((optimize("O0"))) __global__ void global_mem_latency(
unsigned long *time_start,
unsigned long *time_end,
volatile data_type *array_1,
volatile data_type *array_2) {
unsigned long start_time, end_time;
data_type value1 = 0.0f, result;
asm volatile(
// Measure memory load latency directly - no warm-up access
"membar.gl;\n\t"
"mov.u64 %0, %%clock64;\n\t"
"ld.global.cv.f32 %1, [%3];\n\t"
"mov.u64 %4, %%clock64;\n\t"
"st.global.f32 [%3], %1;\n\t"
: "=l"(start_time), "=f"(result), "+f"(value1), "+l"(array_2), "=l"(end_time)
:
: "memory");
*time_start = start_time;
*time_end = end_time;
array_1[0] = result;
}
./mem-latency
global_mem_latency latency = 475 cycles
l2_mem_latency latency = 6 cycles
l1_mem_latency latency = 6 cycles
内存访问合并
同一 Warp（32 个cuda线程）的内存访问应当是连续的。
Coalesced Load（合并加载）：当 Warp 内的各个线程读取连续的内存地址时，GPU 可以通过单次（或极少数几次）内存事务（Transaction）取回所有数据。
Non-coalesced Load（非合并加载/步长加载）：如果各个线程访问的地址步长很大（Stride），内存请求会被打散成大量的离散事务，导致内存延迟飙升，带宽利用率断崖式下跌。
Bank conflict
在 GPU 的 Shared Memory（共享内存） 优化中，Bank Conflict（存储体冲突） 是最常见的性能杀手之一。
为了实现极高的带宽，Shared Memory 被平均分成 32 个等大小的内存模块，称为 Banks（存储体）。
在逻辑上，连续的 4 字节（
float或int32）被轮流映射到这 32 个 Bank 中。规则：第 n 个地址映射到第 n % 32 个 Bank。
在一个 Warp（32 个线程）执行内存指令时，如果所有线程访问的地址分别指向 32 个不同的 Banks，那么这些访问可以 100% 并行完成。
如果 Warp 中有两个或更多线程请求的地址落在 同一个 Bank 中，这些请求就无法同时处理。
硬件必须将这些冲突的请求串行化（Serializing）。例如，如果有 2 个线程冲突，耗时就会翻倍（2-way conflict）。
最典型的例子是跨步访问（Strided Access）：
如果线程 i 访问
shared_data[i * 2]，那么线程 0 访问 Bank 0，线程 1 访问 Bank 2… 看起来没问题。但如果步长是 32 的倍数（例如访问二维数组的列，且行宽是 32），所有 32 个线程都会请求同一个 Bank，导致严重的 32-way conflict，性能瞬间跌至 1/32。
优化技巧
最经典的技巧是 Padding（填充）：
在定义二维 Shared Memory 数组时，故意把列宽增加 1。
例如：原本是
shared float data[32][32]，改为data[32][33]。原理：这样每一行的起始地址在 Bank 中的偏移都会错开 1 位，原本纵向对齐到同一个 Bank 的元素，现在会分布在不同的 Bank 中，从而完美消除冲突。
内存访问密集型高性能实战
Naive 实现
先看 naive gpu实现
template <typename Scene>
__global__ void wave_gpu_naive_step(
float t,
float *u0, /* pointer to GPU memory */
float const *u1 /* pointer to GPU memory */
) {
constexpr int32_t n_cells_x = Scene::n_cells_x;
constexpr int32_t n_cells_y = Scene::n_cells_y;
constexpr float c = Scene::c;
constexpr float dx = Scene::dx;
constexpr float dt = Scene::dt;
const int idx_x = blockIdx.x * blockDim.x + threadIdx.x;
const int idx_y = blockIdx.y * blockDim.y + threadIdx.y;
const int idx = idx_y * n_cells_x + idx_x;
if (idx_x >= n_cells_x || idx_y >= n_cells_y) {
return;
}
bool is_border =
(idx_x == 0 || idx_x == n_cells_x - 1 || idx_y == 0 ||
idx_y == n_cells_y - 1);
float u_next_val;
if (is_border || Scene::is_wall(idx_x, idx_y)) {
u_next_val = 0.0f;
} else if (Scene::is_source(idx_x, idx_y)) {
u_next_val = Scene::source_value(idx_x, idx_y, t);
} else {
constexpr float coeff = c * c * dt * dt / (dx * dx);
float damping = Scene::damping(idx_x, idx_y);
u_next_val =
((2.0f - damping - 4.0f * coeff) * u1[idx] -
(1.0f - damping) * u0[idx] +
coeff *
(u1[idx - 1] + u1[idx + 1] + u1[idx - n_cells_x] +
u1[idx + n_cells_x]));
}
u0[idx] = u_next_val;
}
template <typename Scene>
std::pair<float *, float *> wave_gpu_naive(
float t0,
int32_t n_steps,
float *u0, /* pointer to GPU memory */
float *u1 /* pointer to GPU memory */
) {
const int n_cells_x = Scene::n_cells_x;
const int n_cells_y = Scene::n_cells_y;
dim3 blockDim(block_dim_x, block_dim_y);
dim3 gridDim(
(n_cells_x + block_dim_x - 1) / block_dim_x,
(n_cells_y + block_dim_y - 1) / block_dim_y
);
for (int32_t step = 0; step<n_steps; ++step) {
float t = t0 + step * Scene::dt;
wave_gpu_naive_step<Scene><<<gridDim, blockDim>>>(t, u0, u1);
std::swap(u0, u1);
}
return {u0, u1};
}
最直观的 GPU 实现方式是将上述算法直接翻译为 Kernel 函数：
为网格中的每个像素分配一个线程。
每个线程在 Kernel 内从 Global Memory 读取自己的当前值、历史值和四个邻居的值。
计算结果，写回 Global Memory。
痛点：对于计算每个像素点，我们需要从 Global Memory 读取至少 5 个值并写入 1 个值。在这个过程中，邻居像素被相邻线程重复读取了极多次。面对极高的全局内存带宽压力，计算单元被迫处于饥饿状态等待数据。
分析与对比
Small scale tests (on scene 'DoubleSlitSmallScale'):
CPU sequential implementation:
run time: 603.79 ms
GPU naive implementation:
run time: 4.30 ms
correctness: 3.13e-06 relative RMSE
GPU shared memory implementation:
run time: 1.85 ms
correctness: 4.04e-06 relative RMSE
CPU -> GPU naive speedup: 140.52x
CPU -> GPU shared memory speedup: 326.59x
GPU naive -> GPU shared memory speedup: 2.32x
Large scale tests (on scene 'DoubleSlit'):
GPU naive implementation:
run time: 2363.80 ms
GPU shared memory implementation:
run time: 1500.22 ms
correctness (w.r.t. GPU naive): 9.01e-05 relative RMSE
GPU naive -> GPU shared memory speedup: 1.58x