CUDA 矩阵乘法：分块 vs 朴素实现中的内存合并与共享内存优化

在 CUDA 编程中，矩阵乘法是 GPU 加速的经典应用场景。朴素实现虽然简单，但往往受限于全局内存的高延迟和非合并访问模式，导致吞吐量低下。相比之下，分块（tiled）方法通过引入共享内存和内存合并机制，能显著提升性能。本文将从内存访问效率入手，剖析两种实现的差异，并给出可落地的优化参数和监控要点，帮助开发者在实际项目中实现高效矩阵运算。

朴素实现的内存瓶颈

朴素矩阵乘法内核中，每个线程独立计算 C 矩阵的一个元素，通过循环遍历 A 和 B 的对应行与列。假设矩阵维度为 N × N，每个线程执行的伪代码如下：

__global__ void naive_matmul(float *A, float *B, float *C, int N) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    if (row < N && col < N) {
        float sum = 0.0f;
        for (int k = 0; k < N; k++) {
            sum += A[row * N + k] * B[k * N + col];
        }
        C[row * N + col] = sum;
    }
}

这种实现的核心问题是所有数据访问均来自全局内存（Global Memory），其延迟高达数百个时钟周期，且带宽利用率低。特别是，当 warp（32 个线程组）内的线程访问 A 矩阵的行元素时，由于 row 固定而 k 递增，访问地址是非连续的。例如，一个 warp 的线程可能同时读取 A[row*N + k] 中的不同 k，导致 32 个独立的内存事务，而不是合并成一个宽事务。这违背了内存合并（Coalescing）的原则：warp 内线程应访问连续、对齐的地址，以最小化事务数。

证据显示，在 NVIDIA GPU 上，非合并访问可将有效带宽降低至峰值的 10% 以下。对于 1024×1024 矩阵，朴素内核的 GFLOPS 通常仅为设备峰值的 1-5%。此外，B 矩阵的访问虽可能合并（col 递增），但整体重复读取同一数据（如多个线程读取 B 的同一列）进一步浪费带宽。结果是计算单元（SM）频繁空闲等待内存，整体吞吐量受限。

分块实现的共享内存优势

分块方法将大矩阵分解为小 tile（通常 16×16 或 32×32），每个线程块（block）仅处理一个 C tile。通过共享内存（Shared Memory）缓存 A 和 B 的对应 tile，块内线程可多次复用数据，减少全局内存访问。

典型分块内核引入双重循环：外层循环遍历 K 维度上的 tile，内层在共享内存中计算。核心代码片段如下：

#define TILE_SIZE 16
__global__ void tiled_matmul(float *A, float *B, float *C, int N) {
    __shared__ float sA[TILE_SIZE][TILE_SIZE];
    __shared__ float sB[TILE_SIZE][TILE_SIZE];
    int row = blockIdx.y * TILE_SIZE + threadIdx.y;
    int col = blockIdx.x * TILE_SIZE + threadIdx.x;
    float sum = 0.0f;
    for (int t = 0; t < (N + TILE_SIZE - 1) / TILE_SIZE; t++) {
        if (row < N && t * TILE_SIZE + threadIdx.x < N)
            sA[threadIdx.y][threadIdx.x] = A[row * N + t * TILE_SIZE + threadIdx.x];
        else
            sA[threadIdx.y][threadIdx.x] = 0.0f;
        if (col < N && t * TILE_SIZE + threadIdx.y < N)
            sB[threadIdx.y][threadIdx.x] = B[(t * TILE_SIZE + threadIdx.y) * N + col];
        else
            sB[threadIdx.y][threadIdx.x] = 0.0f;
        __syncthreads();  // 同步，确保加载完成
        for (int k = 0; k < TILE_SIZE; k++) {
            sum += sA[threadIdx.y][k] * sB[k][threadIdx.x];
        }
        __syncthreads();  // 同步前清空共享内存
    }
    if (row < N && col < N)
        C[row * N + col] = sum;
}

共享内存的优势在于其访问延迟仅为全局内存的 1/100，且位于 SM 内，支持块级复用。在加载阶段，warp 内线程加载 A tile 时，通过调整索引（如 threadIdx.x 对应列偏移），可实现合并访问：连续线程读取连续地址，减少事务数至 1-2 个。同样，B tile 的转置加载（threadIdx.y 对应行）也优化了合并。

性能证据：在 A100 GPU 上，对于 2048×2048 矩阵，朴素实现 GFLOPS 约 500，而分块可达 5000+，提升 10 倍。共享内存复用率高：每个 tile 元素被 TILE_SIZE 次使用，全球访问量降至原 1/TILE_SIZE。NVIDIA 文档指出，这种优化将内存流量从 O(N^3) 降至 O(N^3 / TILE_SIZE)，直接提升 Roofline 模型中的计算强度。

然而，共享内存有限（每 SM 约 64KB），需注意银行冲突（Bank Conflict）：若多个线程同时访问同一银行，访问串行化。解决方案：垫奇偶索引或使用向量类型（如 float4）交错访问。

可落地参数与优化清单

要实现最优吞吐量，需根据 GPU 架构（如 Ampere 或 Hopper）调整参数。以下是工程化建议：

Tile 大小选择：
- 起始：TILE_SIZE = 16（适用于大多数 SM，平衡占用率与复用）。
- 高级：32 或 64，对于矩阵 > 1024 时测试；过大可能导致共享内存溢出（检查 cudaFuncGetAttributes）。
- 参数：若 N % TILE_SIZE != 0，使用边界检查填充 0，避免越界。
块维度配置：
- blockDim: (TILE_SIZE, TILE_SIZE, 1)，确保每个块处理一个 C tile。
- gridDim: ( (N + TILE_SIZE - 1)/TILE_SIZE, 同上 )，覆盖整个矩阵。
- 占用率优化：目标 >50%；使用 nsight-compute 分析寄存器与共享内存使用，避免超过 SM 限制（A100: 64K 共享/SM）。
内存合并增强：
- 加载 A 时：使用 threadIdx.x 偏移，确保 warp 连续访问。
- 加载 B 时：转置布局（B 的行变为列），如上代码所示。
- 向量加载：用 float4 替换 float，减少事务 4 倍（需对齐地址）。
同步与监控要点：
- 必用 __syncthreads() 两次：加载后计算前，清空后下一迭代。
- 监控：用 nsys 捕获时间线，检查全局事务数（目标 < 朴素的 1/10）；ncu 分析 L1 命中率 (>80%) 和银行冲突 (0%)。
- 阈值：若吞吐 < 设备峰值 20%，检查 coalescing；< 50%，优化 tile。
回滚策略：
- 若分块失效（小矩阵），fallback 到朴素或 cuBLAS。
- 测试规模：从 512×512 起步，逐步至 8192×8192，记录 GFLOPS = 2*N^3 / 时间 (ms) / 1e9。

通过这些参数，在标准 CUDA（如 11.8+）中，分块实现可接近 cuBLAS 的 90% 性能。实际部署时，结合异步流（streams）进一步重叠计算与传输，提升端到端吞吐。

总之，分块 vs 朴素的对比凸显了 GPU 内存层次的重要性。开发者应优先从 coalescing 入手，逐步引入共享内存，实现从低效到高吞吐的跃升。（字数：1028）