第5章 内存与数据局部性 - 局部性的艺术 (Memory and Data Locality - The Art of Locality)
一、导读
5.1 本章学习目标
通过本章的学习，读者将能够：
深入理解数据局部性的两个维度：时间局部性和空间局部性，以及它们在GPU架构中的体现
掌握GPU内存访问模式的分类和特征，能够识别合并访问、跨步访问、随机访问等模式
学会分析kernel的内存访问行为，使用Nsight Compute工具定位内存瓶颈
掌握数据重排和结构转换技巧（AoS到SoA），优化内存访问效率
理解缓存行和内存事务的概念，能够设计缓存友好的数据结构
掌握共享内存作为软件管理缓存的使用方法，实现高效的数据复用
学会使用预取技术隐藏内存延迟，提高流水线效率
能够综合运用局部性优化策略，显著提升kernel的整体性能
5.2 为什么本章如此重要
数据局部性是计算机体系结构中的核心概念，在GPU编程中尤为重要。GPU拥有极高的计算吞吐量，但如果数据供应跟不上，强大的计算单元就会闲置等待。内存墙（Memory Wall）问题在GPU上比CPU更加严重，因为GPU的线程数量远超CPU，对内存带宽的需求呈指数级增长。
理解并优化数据局部性意味着能够让数据在需要的时候出现在正确的位置。良好的局部性可以提高缓存命中率，减少昂贵的全局内存访问；可以实现合并访问，最大化内存带宽利用率；可以减少冗余数据传输，提高整体能效。
在实际项目中，许多性能低下的kernel并非因为算法复杂度高，而是因为内存访问模式不佳。一个简单的矩阵遍历，如果按列访问而不是按行访问，性能可能下降数倍。一个粒子模拟程序，如果粒子数据在内存中分散存储，缓存命中率会极低。这些问题的根源都在于数据局部性。
本章将从理论和实践两个层面深入探讨数据局部性，帮助读者建立起系统的优化思维。我们不仅会讲解局部性的基本原理，还会提供具体的优化技巧和代码示例，让读者能够将所学知识直接应用到实际项目中。
5.3 核心问题引导
什么是时间局部性和空间局部性？它们在GPU编程中如何体现？
GPU的内存事务是如何工作的？一个warp的内存请求如何被处理？
什么是缓存行（cache line）？缓存行大小对性能有何影响？
AoS（Array of Structures）和SoA（Structure of Arrays）有什么区别？何时应该进行转换？
如何分析一个kernel的内存访问模式？有哪些工具可以帮助诊断？
什么是预取（prefetching）？软件预取在GPU上是否有效？
数据重排（data tiling/blocking）如何提高局部性？如何选择最优的tile大小？
在多GPU系统中，数据局部性有哪些新的挑战和优化机会？
二、核心概念详解
5.1 数据局部性的理论基础
数据局部性（Data Locality）是指程序在运行过程中访问数据的模式具有某种规律性，使得最近访问的数据或附近的数据很可能在不久的将来再次被访问。局部性原理是缓存机制有效的根本原因，也是高性能程序设计的基础。
时间局部性（Temporal Locality）：
时间局部性指的是如果一个数据项被访问，那么在不久的将来它很可能再次被访问。典型的时间局部性场景包括：
- 循环变量：循环计数器在每次迭代中都被访问
- 累加器：归约操作中的累加变量在多次迭代中被反复更新
- 查找表：频繁使用的常量表或映射表
利用时间局部性的方法是将被频繁访问的数据保存在快速的存储层次中（如寄存器、共享内存、缓存），避免重复从慢速存储器中加载。
// Good temporal locality: reusing loaded data
__global__ void goodTemporal(const float* input, float* output, int N) {
__shared__ float shared[256];
int tid = threadIdx.x;
int idx = blockIdx.x * blockDim.x + threadIdx.x;
// Load once into shared memory
if (idx < N) {
shared[tid] = input[idx];
}
__syncthreads();
// Reuse the same data multiple times
float result = 0.0f;
for (int i = 0; i < 10; i++) {
result += shared[(tid + i) % 256] * i;
}
if (idx < N) {
output[idx] = result;
}
}在这个例子中，数据从全局内存加载到共享内存后，被同一Block内的多个线程重复使用10次，充分利用了时间局部性。
空间局部性（Spatial Locality）：
空间局部性指的是如果一个数据项被访问，那么它附近的数据项也很可能在不久的将来被访问。典型的空间局部性场景包括：
- 数组遍历：顺序访问数组元素
- 矩阵运算：访问相邻的行或列元素
- 图像处理：处理相邻像素
利用空间局部性的方法是确保相邻的数据在内存中也相邻存储，并且以连续的顺序访问它们。这样可以充分利用缓存行和内存事务的批量传输特性。
// Good spatial locality: sequential access
__global__ void goodSpatial(const float* input, float* output, int N) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < N) {
// Adjacent threads access adjacent memory locations
output[idx] = input[idx] * 2.0f;
}
}
// Poor spatial locality: strided access
__global__ void poorSpatial(const float* input, float* output, int N, int stride) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx * stride < N) {
// Threads access elements far apart
output[idx] = input[idx * stride] * 2.0f;
}
}在goodSpatial中，相邻线程访问相邻的内存地址，实现了完美的空间局部性。在poorSpatial中，线程访问的元素相隔stride个位置，空间局部性很差。
局部性与内存层次的关系：
GPU的内存层次结构正是为了利用局部性而设计的：
- 寄存器：利用线程内部的时间局部性，保存频繁使用的变量
- 共享内存：利用Block内的时间和空间局部性，缓存复用数据
- L1缓存：利用线程内的空间局部性，自动缓存最近的访问
- L2缓存：利用全局的空间局部性，缓存跨SM的共享数据
- 全局显存：容量最大但速度最慢，是所有数据的最终存储位置
当程序表现出良好的局部性时，数据更可能停留在快速的存储层次中，从而获得更高的性能。反之，如果局部性差，数据会频繁地在不同层次之间迁移，造成性能损失。
5.2 GPU内存访问模式分析
理解GPU如何处理内存请求是优化访问模式的前提。我们需要从硬件的角度审视内存访问的机制。
内存事务（Memory Transaction）：
GPU的全局内存访问以内存事务为单位进行。当一个warp发起内存加载请求时，硬件会将32个线程的请求合并为尽可能少的内存事务。
对于32位数据类型（如float、int），理想的内存事务是：
- warp内的32个线程访问连续的128字节内存区域（32 × 4字节）
- 硬件发出一个128字节的内存事务
- 所有线程的请求在一次事务中得到满足
这种情况下，带宽利用率为100%。
如果线程访问的地址不连续，硬件可能需要发出多个内存事务来满足同一个warp的请求。例如：
- 如果thread i访问地址
base + i * 8（步长为8），则32个线程访问的范围是256字节 - 硬件可能需要发出2个128字节的事务
- 带宽利用率为50%（128字节有用数据 / 256字节传输）
在最坏情况下，每个线程都需要独立的内存事务，带宽利用率仅为3.125%（4字节 / 128字节）。
合并访问的条件：
对于NVIDIA GPU，合并访问的条件可以总结为：
对齐：访问的起始地址应该是128字节对齐的（或至少是32字节对齐）
连续：warp内线程访问的地址应该是连续的
同段：所有地址应该落在同一个或相邻的128字节段内
具体来说，对于32位数据类型，如果thread k访问地址A + k * 4（k = 0, 1, ..., 31），且A是128字节对齐的，则这是一个完美的合并访问。
访问模式分类：
| 模式 | 描述 | 带宽利用率 | 示例 |
|---|---|---|---|
| 合并访问 | 连续地址 | 100% | data[threadIdx.x] |
| 跨步访问 | 固定步长 | 1/stride | data[threadIdx.x * stride] |
| 广播访问 | 相同地址 | 100%* | data[0]（所有线程） |
| 随机访问 | 无规律 | 很低 | data[randomIdx[threadIdx.x]] |
*广播访问虽然所有线程访问相同地址，但硬件可以优化为单次读取并广播给所有线程。
二维数据的访问模式：
对于存储在行优先（row-major）顺序的二维数组，按行访问具有良好的空间局部性，按列访问则较差：
// Row-major storage: M rows x N columns
// Element (i, j) is at index i * N + j
// Good: accessing by row (consecutive columns)
__global__ void rowAccess(float* matrix, int N) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
int col = idx % N;
// Adjacent threads access adjacent columns in the same row
}
// Poor: accessing by column (stride = N)
__global__ void colAccess(float* matrix, int N) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
int row = idx / N;
// Adjacent threads access elements N apart
}5.3 数据结构优化：AoS vs SoA
数据结构的选择对内存访问模式有深远影响。AoS（Array of Structures）和SoA（Structure of Arrays）是两种常见的数据组织方式，它们在不同的访问模式下表现各异。
AoS（Array of Structures）：
AoS将相关的数据字段组织在一起，形成一个结构体，然后将结构体放入数组中。这是面向对象编程中常见的模式。
struct Particle {
float x, y, z; // Position
float vx, vy, vz; // Velocity
float mass;
};
Particle particles[N]; // AoS layout在内存中，AoS的布局为：
[x0, y0, z0, vx0, vy0, vz0, mass0, x1, y1, z1, vx1, vy1, vz1, mass1, ...]SoA（Structure of Arrays）：
SoA将每个字段分别存储在独立的数组中，然后通过索引关联。
struct Particles {
float x[N], y[N], z[N];
float vx[N], vy[N], vz[N];
float mass[N];
};在内存中，SoA的布局为：
[x0, x1, x2, ..., xN-1, y0, y1, y2, ..., yN-1, z0, z1, z2, ..., zN-1, ...]性能对比：
当kernel只需要访问结构体的部分字段时，SoA通常表现更好：
// AoS: loading all fields even if only position is needed
__global__ void updatePositionAoS(Particle* particles, int N, float dt) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < N) {
// Loads entire struct (28 bytes) but only uses x, y, z (12 bytes)
particles[idx].x += particles[idx].vx * dt;
particles[idx].y += particles[idx].vy * dt;
particles[idx].z += particles[idx].vz * dt;
}
}
// SoA: loading only needed fields
__global__ void updatePositionSoA(Particles* p, int N, float dt) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < N) {
// Only loads x, y, z arrays, better bandwidth utilization
p->x[idx] += p->vx[idx] * dt;
p->y[idx] += p->vy[idx] * dt;
p->z[idx] += p->vz[idx] * dt;
}
}在AoS版本中，即使只需要更新位置，硬件也需要加载整个结构体（包括速度和质量字段），浪费了带宽。在SoA版本中，只加载需要的数组，带宽利用率更高。
选择准则：
- 使用SoA：当kernel通常只访问部分字段，且需要高带宽时
- 使用AoS：当kernel通常访问所有字段，或需要保持数据封装时
- 混合布局：有时可以采用折中方案，如将相关字段分组为小的SoA
转换方法：
如果现有代码使用AoS，可以通过以下方式转换为SoA：
重新定义数据结构，将字段分离为独立数组
修改数据初始化代码，将数据从AoS格式复制到SoA格式
更新所有访问数据的kernel代码
测试验证正确性和性能提升
这种转换可能需要较大的代码改动，但对于性能敏感的应用，收益通常是显著的。
5.4 缓存行与内存块
缓存行（Cache Line）：
缓存是GPU内存层次中的重要组成部分，而缓存行是缓存的基本单位。当GPU从全局内存加载数据到缓存时，不是加载单个字，而是加载整个缓存行。
NVIDIA GPU的缓存行大小通常为128字节。这意味着：
- 当访问某个地址时，包含该地址的128字节块会被加载到缓存中
- 后续访问同一缓存行内的其他地址可以直接从缓存中获取，无需再次访问全局内存
- 如果访问跨越多个缓存行，则需要多次内存事务
利用缓存行：
为了充分利用缓存行，应该确保：
顺序访问：按内存顺序访问数据，使得预取的缓存行中的数据能被充分利用
数据紧凑：将经常一起访问的数据放在一起，使其落在同一个缓存行内
避免填充浪费：不要在数据结构中添加不必要的填充字节，这会浪费缓存行空间
// Good: compact structure fits well in cache lines
struct CompactParticle {
float x, y, z; // 12 bytes, tightly packed
};
// Wasteful: padding wastes cache space
struct PaddedParticle {
float x;
double padding; // Unnecessary padding
float y;
float z;
};缓存替换策略：
当缓存满时，需要决定哪些数据被替换出去。GPU通常使用LRU（Least Recently Used）或类似的策略，替换最近最少使用的数据。
理解缓存替换策略有助于预测哪些数据会保留在缓存中：
- 最近访问的数据更可能留在缓存中
- 大量一次性访问的数据会挤出有用的缓存内容
- 循环访问的工作集如果超过缓存容量，会导致频繁的缓存失效
工作集大小：
工作集（Working Set）是指程序在一段时间内访问的数据总量。如果工作集可以放入L1缓存，性能会非常好；如果能放入L2缓存，性能也不错；如果超出L2缓存，性能会显著下降。
估算工作集大小的方法：
工作集 = 活跃线程数 × 每线程访问的唯一数据量例如，如果一个kernel有10000个活跃线程，每个线程访问4个唯一的float值，则工作集约为10000 × 4 × 4 = 160KB。如果L2缓存为8MB，工作集完全可以容纳，缓存命中率会很高。
5.5 数据分块与Tiling
数据分块（Tiling/Blocking）是一种经典的局部性优化技术，通过将大数据集划分为小块，使得每个块可以放入快速存储器中，从而提高数据复用率。
矩阵乘法中的Tiling：
前面章节已经展示了使用共享内存的矩阵乘法优化，其核心思想就是Tiling：
#define TILE_SIZE 16
__global__ void matMulTiled(const float* A, const float* B, float* C,
int M, int N, int K) {
__shared__ float As[TILE_SIZE][TILE_SIZE];
__shared__ float Bs[TILE_SIZE][TILE_SIZE];
int tx = threadIdx.x;
int ty = threadIdx.y;
int row = blockIdx.y * TILE_SIZE + ty;
int col = blockIdx.x * TILE_SIZE + tx;
float sum = 0.0f;
for (int t = 0; t < (K + TILE_SIZE - 1) / TILE_SIZE; t++) {
// Load tile of A
int aCol = t * TILE_SIZE + tx;
if (row < M && aCol < K) {
As[ty][tx] = A[row * K + aCol];
} else {
As[ty][tx] = 0.0f;
}
// Load tile of B
int bRow = t * TILE_SIZE + ty;
if (bRow < K && col < N) {
Bs[ty][tx] = B[bRow * N + col];
} else {
Bs[ty][tx] = 0.0f;
}
__syncthreads();
// Compute using tile data
for (int k = 0; k < TILE_SIZE; k++) {
sum += As[ty][k] * Bs[k][tx];
}
__syncthreads();
}
if (row < M && col < N) {
C[row * N + col] = sum;
}
}在这个实现中：
- 每个TILE×TILE的子矩阵从全局内存加载一次到共享内存
- 然后被TILE个线程复用TILE次
- 全局内存访问减少了TILE倍（对于16×16的tile，减少16倍）
选择Tile大小：
Tile大小的选择需要考虑以下因素：
共享内存容量：Tile不能太大，否则无法放入共享内存
Occupancy：较大的Tile使用更多共享内存，可能降低Occupancy
数据复用率：较大的Tile提供更高的复用率，但边际收益递减
Bank Conflict：Tile大小应避免导致bank conflict
经验法则：从16×16或32×8开始测试，根据性能测量调整。
多维Tiling：
对于三维或更高维的数据，可以进行多维Tiling：
// 3D convolution with 3D tiling
#define TILE_X 8
#define TILE_Y 8
#define TILE_Z 8
__global__ void conv3DTiled(...) {
__shared__ float shared[TILE_Z + 2][TILE_Y + 2][TILE_X + 2];
// +2 for halo regions needed for convolution
}额外的边界（halo）用于存储卷积所需的邻域数据。
5.6 预取与延迟隐藏
预取（Prefetching）：
预取是指在数据被实际需要之前，提前将其加载到更快的存储层次中。预取可以是硬件自动进行的，也可以是软件显式控制的。
硬件预取：
现代GPU包含硬件预取器，可以检测顺序访问模式并提前加载数据。硬件预取的优势是透明性，程序员无需修改代码即可受益。但硬件预取只能识别简单的模式（如顺序访问、固定步长访问），对于复杂的访问模式无效。
软件预取：
在CPU编程中，软件预取通常使用特殊的预取指令（如x86的PREFETCH）。在GPU上，软件预取的效果有限，因为：
GPU已经有大量的活跃线程来隐藏延迟
显式的预取加载会占用寄存器和内存带宽
GPU的硬件预取器已经相当智能
然而，在某些情况下，软件预取仍然有用：
// Software prefetching example: load next tile while processing current
__global__ void prefetchExample(const float* input, float* output, int N) {
extern __shared__ float shared[];
int tid = threadIdx.x;
int blockSize = blockDim.x;
// Process data in chunks
for (int chunk = 0; chunk < N / blockSize; chunk++) {
int currentIdx = chunk * blockSize + tid;
int nextIdx = (chunk + 1) * blockSize + tid;
// Load current data
float current = input[currentIdx];
// Prefetch next chunk (if not last)
if (nextIdx < N && tid < blockSize) {
// This could be optimized with async copy in newer architectures
shared[tid] = input[nextIdx];
}
__syncthreads();
// Process current
output[currentIdx] = process(current);
__syncthreads();
}
}异步内存复制：
在较新的GPU架构（如Ampere及以后）中，引入了异步内存复制指令（cp.async），允许在后台进行内存加载，同时执行其他计算。这是一种更高效的预取机制：
// Using async copy (requires compute capability 8.0+)
__global__ void asyncCopyExample(...) {
// Initiate async load
cp.async.cg.shared.global[sharedAddr] = &globalAddr[size];
// Do other work while load is in progress
// Wait for async load to complete
cp.async.commit_group();
cp.async.wait_group(0);
}延迟隐藏的其他方法：
除了预取，还有其他方法可以隐藏内存延迟：
增加Occupancy：更多的活跃warp可以在某个warp等待内存时切换到其他warp
指令级并行：编译器重新排序指令，使得长延迟指令的结果在使用前有时间完成
多Kernel重叠：使用CUDA流并发执行多个小kernel
三、程序员的视角
3.1 在实际项目中的应用
粒子系统模拟：在物理引擎中，粒子的位置和速度需要频繁更新。使用SoA布局可以显著提高带宽利用率，因为位置更新kernel只需加载位置数组，无需加载质量等其他字段。邻近搜索使用空间哈希或网格分块，提高缓存命中率。
深度学习推理：卷积神经网络的权重通常以特定的内存布局存储（如NHWC或NCHW），以匹配卷积操作的访问模式。Winograd算法通过数学变换减少乘法次数，但增加了数据重排的开销，需要在计算节省和数据移动成本之间权衡。
数据库查询加速：列式存储（Columnar Storage）本质上是SoA布局，适合OLAP查询，因为查询通常只涉及少数几列。GPU加速的连接操作使用哈希表构建和探测，哈希表的大小需要适应L2缓存以获得最佳性能。
实时光线追踪：BVH（Bounding Volume Hierarchy）树的遍历涉及大量的指针追逐，缓存友好性至关重要。使用紧凑的节点表示和广度优先遍历可以提高缓存命中率。
3.2 常见陷阱与误区
忽视数据结构的影响：继续使用AoS布局，即使kernel只访问部分字段，导致带宽浪费。应该在性能关键路径上考虑SoA转换。
过度优化局部性：对于计算密集型kernel，内存可能不是瓶颈，过度优化局部性可能带来代码复杂度增加而收益甚微。应该先进行性能剖析，确定真正的瓶颈。
错误的Tile大小：选择的Tile大小导致共享内存不足或bank conflict。应该使用Nsight Compute检查共享内存使用和bank conflict情况。
假设硬件预取万能：依赖硬件预取处理复杂的访问模式，实际上预取器只能处理简单模式。对于不规则访问，需要显式优化。
忽视写操作的成本：只关注读操作的优化，忽略写操作同样需要内存事务。写回结果时也应确保合并访问。
3.3 与其他层级的联系
与算法设计的联系：局部性优化不仅仅是底层实现细节，也影响算法选择。例如，外部排序算法（如归并排序）比内部排序算法（如快速排序）更适合大规模数据，因为它们更好地利用了I/O局部性。
与编译器的联系：编译器可以进行一些自动的局部性优化，如循环交换（loop interchange）、循环分块（loop tiling）。但编译器的优化能力有限，程序员仍需要进行高层次的数据结构优化。
与硬件演进的联系：随着缓存容量的增加和预取器的智能化，一些传统的局部性优化可能变得不那么重要。但基本原则（如合并访问、数据复用）仍然适用。
四、代码实战
4.1 AoS到SoA转换性能对比
#include <cuda_runtime.h>
#include <stdio.h>
#define N 10000000
// AoS structure
struct ParticleAoS {
float x, y, z;
float vx, vy, vz;
float mass;
float padding; // Align to 32 bytes
};
// SoA structure
struct ParticleSoA {
float* x;
float* y;
float* z;
float* vx;
float* vy;
float* vz;
float* mass;
};
// AoS kernel
__global__ void updateVelocityAoS(ParticleAoS* particles, int N, float ax, float ay, float az, float dt) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < N) {
particles[idx].vx += ax * dt;
particles[idx].vy += ay * dt;
particles[idx].vz += az * dt;
}
}
// SoA kernel
__global__ void updateVelocitySoA(ParticleSoA particles, int N, float ax, float ay, float az, float dt) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < N) {
particles.vx[idx] += ax * dt;
particles.vy[idx] += ay * dt;
particles.vz[idx] += az * dt;
}
}
// Benchmark
void benchmarkAoSvsSoA() {
int aosSize = N * sizeof(ParticleAoS);
ParticleAoS* h_aos = (ParticleAoS*)malloc(aosSize);
// Initialize
for (int i = 0; i < N; i++) {
h_aos[i].x = h_aos[i].y = h_aos[i].z = 0.0f;
h_aos[i].vx = h_aos[i].vy = h_aos[i].vz = 1.0f;
h_aos[i].mass = 1.0f;
}
// Allocate AoS on device
ParticleAoS* d_aos;
cudaMalloc((void**)&d_aos, aosSize);
cudaMemcpy(d_aos, h_aos, aosSize, cudaMemcpyHostToDevice);
// Allocate SoA on device
ParticleSoA d_soa;
cudaMalloc((void**)&d_soa.x, N * sizeof(float));
cudaMalloc((void**)&d_soa.y, N * sizeof(float));
cudaMalloc((void**)&d_soa.z, N * sizeof(float));
cudaMalloc((void**)&d_soa.vx, N * sizeof(float));
cudaMalloc((void**)&d_soa.vy, N * sizeof(float));
cudaMalloc((void**)&d_soa.vz, N * sizeof(float));
cudaMalloc((void**)&d_soa.mass, N * sizeof(float));
// Copy to SoA format
float *tempVx = (float*)malloc(N * sizeof(float));
float *tempVy = (float*)malloc(N * sizeof(float));
float *tempVz = (float*)malloc(N * sizeof(float));
for (int i = 0; i < N; i++) {
tempVx[i] = h_aos[i].vx;
tempVy[i] = h_aos[i].vy;
tempVz[i] = h_aos[i].vz;
}
cudaMemcpy(d_soa.vx, tempVx, N * sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(d_soa.vy, tempVy, N * sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(d_soa.vz, tempVz, N * sizeof(float), cudaMemcpyHostToDevice);
int blockSize = 256;
int gridSize = (N + blockSize - 1) / blockSize;
cudaEvent_t start, stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
float ax = 0.1f, ay = 0.2f, az = 0.3f, dt = 0.01f;
// Benchmark AoS
cudaEventRecord(start);
for (int i = 0; i < 100; i++) {
updateVelocityAoS<<<gridSize, blockSize>>>(d_aos, N, ax, ay, az, dt);
}
cudaEventRecord(stop);
cudaEventSynchronize(stop);
float msAoS = 0;
cudaEventElapsedTime(&msAoS, start, stop);
// Benchmark SoA
cudaEventRecord(start);
for (int i = 0; i < 100; i++) {
updateVelocitySoA<<<gridSize, blockSize>>>(d_soa, N, ax, ay, az, dt);
}
cudaEventRecord(stop);
cudaEventSynchronize(stop);
float msSoA = 0;
cudaEventElapsedTime(&msSoA, start, stop);
printf("=== AoS vs SoA Performance ===\n");
printf("AoS: %.3f ms (100 iterations)\n", msAoS);
printf("SoA: %.3f ms (100 iterations)\n", msSoA);
printf("Speedup: %.2fx\n", msAoS / msSoA);
// Cleanup
cudaFree(d_aos);
cudaFree(d_soa.x);
cudaFree(d_soa.y);
cudaFree(d_soa.z);
cudaFree(d_soa.vx);
cudaFree(d_soa.vy);
cudaFree(d_soa.vz);
cudaFree(d_soa.mass);
free(h_aos);
free(tempVx);
free(tempVy);
free(tempVz);
cudaEventDestroy(start);
cudaEventDestroy(stop);
}
int main() {
benchmarkAoSvsSoA();
return 0;
}4.2 缓存友好的矩阵遍历
#include <cuda_runtime.h>
#include <stdio.h>
// Row-major traversal (good spatial locality)
__global__ void rowMajorTraversal(const float* matrix, float* output,
int rows, int cols) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
int totalElements = rows * cols;
if (idx < totalElements) {
// Sequential access in row-major order
output[idx] = matrix[idx] * 2.0f;
}
}
// Column-major traversal (poor spatial locality for row-major stored matrix)
__global__ void colMajorTraversal(const float* matrix, float* output,
int rows, int cols) {
int col = blockIdx.x * blockDim.x + threadIdx.x;
int row = blockIdx.y * blockDim.y + threadIdx.y;
if (col < cols && row < rows) {
// Access by column: stride = cols
int idx = row * cols + col;
output[idx] = matrix[idx] * 2.0f;
}
}
// Tiled traversal (optimized with shared memory)
#define TILE_SIZE 16
__global__ void tiledTraversal(const float* matrix, float* output,
int rows, int cols) {
__shared__ float tile[TILE_SIZE][TILE_SIZE];
int tx = threadIdx.x;
int ty = threadIdx.y;
int baseRow = blockIdx.y * TILE_SIZE;
int baseCol = blockIdx.x * TILE_SIZE;
int row = baseRow + ty;
int col = baseCol + tx;
// Load tile with coalesced access
if (row < rows && col < cols) {
tile[ty][tx] = matrix[row * cols + col];
}
__syncthreads();
// Process from shared memory
if (row < rows && col < cols) {
output[row * cols + col] = tile[ty][tx] * 2.0f;
}
}
void benchmarkMatrixTraversal(int rows, int cols) {
size_t size = rows * cols * sizeof(float);
float *h_matrix = (float*)malloc(size);
float *h_output = (float*)malloc(size);
for (int i = 0; i < rows * cols; i++) {
h_matrix[i] = (float)i;
}
float *d_matrix, *d_output;
cudaMalloc((void**)&d_matrix, size);
cudaMalloc((void**)&d_output, size);
cudaMemcpy(d_matrix, h_matrix, size, cudaMemcpyHostToDevice);
cudaEvent_t start, stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
// Row-major
int totalElements = rows * cols;
int blockSize = 256;
int gridSize = (totalElements + blockSize - 1) / blockSize;
cudaEventRecord(start);
for (int i = 0; i < 50; i++) {
rowMajorTraversal<<<gridSize, blockSize>>>(d_matrix, d_output, rows, cols);
}
cudaEventRecord(stop);
cudaEventSynchronize(stop);
float msRow = 0;
cudaEventElapsedTime(&msRow, start, stop);
// Column-major
dim3 blockSize2D(16, 16);
dim3 gridSize2D((cols + 15) / 16, (rows + 15) / 16);
cudaEventRecord(start);
for (int i = 0; i < 50; i++) {
colMajorTraversal<<<gridSize2D, blockSize2D>>>(d_matrix, d_output, rows, cols);
}
cudaEventRecord(stop);
cudaEventSynchronize(stop);
float msCol = 0;
cudaEventElapsedTime(&msCol, start, stop);
// Tiled
dim3 tileSize(TILE_SIZE, TILE_SIZE);
dim3 gridTiles((cols + TILE_SIZE - 1) / TILE_SIZE,
(rows + TILE_SIZE - 1) / TILE_SIZE);
cudaEventRecord(start);
for (int i = 0; i < 50; i++) {
tiledTraversal<<<gridTiles, tileSize>>>(d_matrix, d_output, rows, cols);
}
cudaEventRecord(stop);
cudaEventSynchronize(stop);
float msTiled = 0;
cudaEventElapsedTime(&msTiled, start, stop);
printf("=== Matrix Traversal Benchmark (%dx%d) ===\n", rows, cols);
printf("Row-major: %.3f ms\n", msRow / 50);
printf("Column-major: %.3f ms\n", msCol / 50);
printf("Tiled: %.3f ms\n", msTiled / 50);
cudaFree(d_matrix);
cudaFree(d_output);
free(h_matrix);
free(h_output);
cudaEventDestroy(start);
cudaEventDestroy(stop);
}
int main() {
benchmarkMatrixTraversal(2048, 2048);
return 0;
}4.3 数据分块优化卷积
#include <cuda_runtime.h>
#include <stdio.h>
#define FILTER_RADIUS 2
#define FILTER_SIZE (2 * FILTER_RADIUS + 1)
#define TILE_SIZE 16
// Naive 2D convolution
__global__ void convolutionNaive(const float* input, float* output,
const float* filter,
int width, int height) {
int x = blockIdx.x * blockDim.x + threadIdx.x;
int y = blockIdx.y * blockDim.y + threadIdx.y;
if (x < width && y < height) {
float sum = 0.0f;
for (int fy = -FILTER_RADIUS; fy <= FILTER_RADIUS; fy++) {
for (int fx = -FILTER_RADIUS; fx <= FILTER_RADIUS; fx++) {
int px = min(width - 1, max(0, x + fx));
int py = min(height - 1, max(0, y + fy));
int filterIdx = (fy + FILTER_RADIUS) * FILTER_SIZE +
(fx + FILTER_RADIUS);
sum += input[py * width + px] * filter[filterIdx];
}
}
output[y * width + x] = sum;
}
}
// Optimized convolution with shared memory tiling
__global__ void convolutionTiled(const float* input, float* output,
const float* filter,
int width, int height) {
// Shared memory with halo region
__shared__ float shared[TILE_SIZE + 2 * FILTER_RADIUS]
[TILE_SIZE + 2 * FILTER_RADIUS];
int tx = threadIdx.x;
int ty = threadIdx.y;
int baseX = blockIdx.x * TILE_SIZE;
int baseY = blockIdx.y * TILE_SIZE;
int x = baseX + tx;
int y = baseY + ty;
// Load input into shared memory with halo
// Each thread may load multiple elements
for (int ly = ty; ly < TILE_SIZE + 2 * FILTER_RADIUS; ly += blockDim.y) {
for (int lx = tx; lx < TILE_SIZE + 2 * FILTER_RADIUS; lx += blockDim.x) {
int px = min(width - 1, max(0, baseX + lx - FILTER_RADIUS));
int py = min(height - 1, max(0, baseY + ly - FILTER_RADIUS));
shared[ly][lx] = input[py * width + px];
}
}
__syncthreads();
// Compute convolution from shared memory
if (tx < TILE_SIZE && ty < TILE_SIZE && x < width && y < height) {
float sum = 0.0f;
for (int fy = 0; fy < FILTER_SIZE; fy++) {
for (int fx = 0; fx < FILTER_SIZE; fx++) {
sum += shared[ty + fy][tx + fx] * filter[fy * FILTER_SIZE + fx];
}
}
output[y * width + x] = sum;
}
}
void benchmarkConvolution(int width, int height) {
size_t imageSize = width * height * sizeof(float);
size_t filterSize = FILTER_SIZE * FILTER_SIZE * sizeof(float);
float *h_input = (float*)malloc(imageSize);
float *h_output = (float*)malloc(imageSize);
float *h_filter = (float*)malloc(filterSize);
// Initialize
for (int i = 0; i < width * height; i++) {
h_input[i] = (float)(i % 256);
}
// Gaussian filter
float sigma = 1.0f;
float sum = 0.0f;
for (int fy = -FILTER_RADIUS; fy <= FILTER_RADIUS; fy++) {
for (int fx = -FILTER_RADIUS; fx <= FILTER_RADIUS; fx++) {
int idx = (fy + FILTER_RADIUS) * FILTER_SIZE + (fx + FILTER_RADIUS);
h_filter[idx] = expf(-(fx * fx + fy * fy) / (2 * sigma * sigma));
sum += h_filter[idx];
}
}
for (int i = 0; i < FILTER_SIZE * FILTER_SIZE; i++) {
h_filter[i] /= sum;
}
float *d_input, *d_output, *d_filter;
cudaMalloc((void**)&d_input, imageSize);
cudaMalloc((void**)&d_output, imageSize);
cudaMalloc((void**)&d_filter, filterSize);
cudaMemcpy(d_input, h_input, imageSize, cudaMemcpyHostToDevice);
cudaMemcpy(d_filter, h_filter, filterSize, cudaMemcpyHostToDevice);
cudaEvent_t start, stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
// Naive
dim3 blockSize1(16, 16);
dim3 gridSize1((width + 15) / 16, (height + 15) / 16);
cudaEventRecord(start);
for (int i = 0; i < 30; i++) {
convolutionNaive<<<gridSize1, blockSize1>>>(d_input, d_output,
d_filter, width, height);
}
cudaEventRecord(stop);
cudaEventSynchronize(stop);
float msNaive = 0;
cudaEventElapsedTime(&msNaive, start, stop);
// Tiled
dim3 blockSize2(TILE_SIZE, TILE_SIZE);
dim3 gridSize2((width + TILE_SIZE - 1) / TILE_SIZE,
(height + TILE_SIZE - 1) / TILE_SIZE);
cudaEventRecord(start);
for (int i = 0; i < 30; i++) {
convolutionTiled<<<gridSize2, blockSize2>>>(d_input, d_output,
d_filter, width, height);
}
cudaEventRecord(stop);
cudaEventSynchronize(stop);
float msTiled = 0;
cudaEventElapsedTime(&msTiled, start, stop);
printf("=== Convolution Benchmark (%dx%d, filter %dx%d) ===\n",
width, height, FILTER_SIZE, FILTER_SIZE);
printf("Naive: %.3f ms\n", msNaive / 30);
printf("Tiled: %.3f ms\n", msTiled / 30);
printf("Speedup: %.2fx\n", msNaive / msTiled);
cudaFree(d_input);
cudaFree(d_output);
cudaFree(d_filter);
free(h_input);
free(h_output);
free(h_filter);
cudaEventDestroy(start);
cudaEventDestroy(stop);
}
int main() {
benchmarkConvolution(1024, 1024);
return 0;
}五、深入思考
5.1 历史背景
数据局部性的概念最早由Peter Denning在1968年提出，他发现了程序访问内存的"局部性原理"。这一发现为缓存机制奠定了理论基础，成为现代计算机体系结构的核心原则之一。
在GPU发展的早期阶段，程序员主要关注计算吞吐量，内存优化相对次要。但随着GPU应用于更广泛的领域，内存瓶颈日益凸显。NVIDIA在Fermi架构中引入L1/L2缓存，标志着GPU开始更加重视自动化的局部性优化。
近年来，随着AI工作负载的兴起，数据局部性的研究重点转向了张量操作的优化。Tensor Core的使用需要特定的数据布局（如半精度、特定的矩阵分块），这对传统的局部性优化提出了新的挑战。
5.2 前沿方向
非易失性内存：Intel Optane等非易失性内存技术模糊了内存和存储的界限。未来的GPU可能直接访问持久化内存，数据局部性的概念将扩展到存储层次。
近存计算：HBM-PIM等技术将计算单元集成到内存芯片中，从根本上改变数据移动的范式。在这种架构下，数据局部性可能意味着将计算推向数据所在的位置。
机器学习辅助预取：使用轻量级的机器学习模型预测未来的内存访问模式，指导硬件预取器。这种方法可以处理比传统预取器更复杂的模式。
量子内存：虽然仍处于早期研究阶段，量子内存的访问模式与经典内存完全不同。局部性的概念可能需要重新定义。
5.3 延伸阅读
Peter J. Denning, "The Locality Principle", Communications of the ACM, 2005
NVIDIA CUDA Best Practices Guide - Memory Optimizations
"Cache-Oblivious Algorithms", Frigo et al., 1999
NVIDIA Developer Blog: "Using Shared Memory for CUDA C/C++"
六、实验笔记
实验1：AoS vs SoA性能对比
目的：量化数据结构对内存带宽利用率的影响。
步骤：
实现AoS和SoA两种粒子数据结构
编写仅更新速度的kernel
使用Nsight Compute测量带宽利用率
比较两种布局的执行时间
预期：SoA应该有更高的带宽利用率，因为只加载需要的字段。加速比取决于结构体大小和被访问字段的比例。
实验2：Tile大小对卷积性能的影响
目的：找到卷积操作的最优tile大小。
步骤：
实现不同tile大小（8x8, 16x16, 32x32）的卷积kernel
测量每种配置的執行时间和共享内存使用
分析occupancy与性能的关系
预期：中等tile大小（16x16）通常在性能和资源使用之间取得平衡。过小的tile共享内存开销占比大，过大的tile降低occupancy。
七、本章小结
本章深入探讨了数据局部性在GPU编程中的重要性。我们从时间局部性和空间局部性两个维度出发，分析了局部性原理如何指导内存访问优化。GPU的内存层次结构正是为了利用局部性而设计的：寄存器捕获线程内的时间局部性，共享内存利用Block内的数据复用，缓存自动捕获最近的访问模式。
合并访问是GPU内存优化的核心要求，它确保warp内的线程以连续的方式访问内存，最大化带宽利用率。我们详细分析了合并访问的条件和非合并访问的性能惩罚，提供了识别和优化不良访问模式的方法。
数据结构的选择对局部性有深远影响。AoS布局适合面向对象的编程风格，但在GPU上可能导致带宽浪费；SoA布局虽然编程复杂度稍高，但能显著提高带宽利用率。在性能关键的应用中，从AoS到SoA的转换往往是值得的。
数据分块（Tiling）是一种经典的局部性优化技术，通过将大数据集划分为适合快速存储器的小块，实现数据的高效复用。矩阵乘法和卷积是Tiling的典型应用场景，合理的tile大小选择需要在共享内存使用和occupancy之间权衡。
预取技术可以进一步隐藏内存延迟，但在GPU上的效果因架构而异。较新的GPU提供了异步内存复制指令，为软件控制的预取提供了更强大的工具。
通过三个代码示例，我们实践了AoS/SoA对比、缓存友好的矩阵遍历和分块卷积优化，将理论知识转化为实际的编程技能。
八、习题
局部性分类：以下访问模式分别体现了时间局部性还是空间局部性？
a) 循环中反复使用同一个变量
b) 顺序遍历数组
c) 哈希表查找
d) 矩阵的行优先遍历
合并访问计算：对于double类型（8字节）数组，warp内的thread i访问data[i]。这是否是合并访问？需要多少个128字节的内存事务？如果改为float类型呢？
AoS/SoA转换：给定结构体struct { int id; float x, y, z; double timestamp; }，大小为32字节。如果kernel只访问x、y、z字段，从AoS转换为SoA可以减少多少带宽消耗？
Cache行分析：假设缓存行为128字节，一个kernel顺序访问float数组。每次缓存失效会加载多少个float到缓存中？如果后续访问是随机的，缓存命中率大约是多少？
Tile大小选择：对于共享内存为112KB的SM，如果使用float shared[TILE_SIZE + 4][TILE_SIZE + 4]（含边界），最大的TILE_SIZE是多少？考虑4字节对齐。
预取有效性：在一个流式处理应用中，数据只访问一次然后丢弃。硬件预取是否有效？软件预取是否有意义？说明理由。
工作集估算：一个kernel有2048个活跃线程，每个线程访问8个唯一的float值。工作集大小是多少？如果L2缓存为6MB，工作集是否可以完全容纳？
矩阵转置优化：编写一个使用共享内存优化的矩阵转置kernel，确保读写都是合并访问。与非优化版本比较性能。
稀疏矩阵局部性：CSR格式的稀疏矩阵向量乘法（SpMV）的内存访问模式有什么特点？为什么SpMV通常难以优化局部性？提出一种可能的改进方法。
综合优化：设计一个高效的并行直方图均衡化算法。说明如何使用共享内存减少全局内存访问，如何优化数据布局以提高局部性，以及如何选择合适的block和tile大小。