使用共享内存和分块矩阵优化矩阵乘法

一、朴素矩阵乘法实现

线程分工

CUDA 矩阵乘法中，一个线程只负责计算并输出结果矩阵 P 的一个元素。这条规则在朴素版和分块版中完全一致，两者的差异只在于数据的读取方式。

对于 Width × Width 的方阵，线程坐标到矩阵坐标的映射为：

Row = blockIdx.y × TILE_WIDTH + threadIdx.y // 对应 P 的行
Col = blockIdx.x × TILE_WIDTH + threadIdx.x // 对应 P 的列

每个线程计算 P[Row][Col]，需要对 M 的第 Row 行与 N 的第 Col 列做完整点积：

完整代码

#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>

#define TILE_WIDTH 16

// 朴素矩阵乘法核函数
__global__ void MatrixMulNaive(float *M, float *N, float *P, int Width) {

    // 当前线程负责的 P 矩阵坐标
    int Row = blockIdx.y * TILE_WIDTH + threadIdx.y;
    int Col = blockIdx.x * TILE_WIDTH + threadIdx.x;

    if (Row < Width && Col < Width) {
        float Pvalue = 0.0f;

        // 点积：每次循环都直接读取全局内存
        for (int k = 0; k < Width; k++) {
            Pvalue += M[Row * Width + k] * N[k * Width + Col];
        }

        P[Row * Width + Col] = Pvalue;
    }
}

与分块版的区别仅在核函数：主函数的内存管理、启动配置完全相同，差异全部集中在核函数内部的数据读取方式。

二、朴素实现为什么慢？

2.1 全局内存的访问代价

GPU 的内存层次结构从快到慢依次为：

寄存器 → 共享内存（片上）→ L1/L2 缓存 → 全局内存（DRAM，片外）

朴素版的点积循环每轮读取 M[Row][k] 和 N[k][Col] 各一次，两次都是全局内存访问。对于 Width=1024 的矩阵，每个线程的内层循环执行 1024 次，每次都要忍受几百个时钟周期的延迟等待。

2.2 数据被重复加载

考虑计算 P 矩阵同一行的所有元素：

• 线程 (Row=0, Col=0) 需要读取 M[0][0], M[0][1], ..., M[0][Width-1]
• 线程 (Row=0, Col=1) 需要读取完全相同的一行 M[0][0], M[0][1], ..., M[0][Width-1]
• ……
• 线程 (Row=0, Col=Width-1) 再次读取同一行

一行 M 数据被同一线程块内的 TILE_WIDTH 个线程各自独立地从全局内存读取了一遍，没有任何复用。对于 N 矩阵也存在同样的问题：同一列被多个线程重复加载。

全局内存带宽被大量浪费在重复搬运同一份数据上。

2.3 访存比（Arithmetic Intensity）极低

访存比 = 浮点运算量 ÷ 全局内存流量，单位 FLOP/B，衡量每字节内存带宽对应多少有效计算。

朴素版中：

• 每次循环读取 2 个 float（M 和 N 各一个）= 8 字节
• 完成 1 次乘法 + 1 次加法 = 2 FLOP
• 访存比 = 2 ÷ 8 = 0.25 FLOP/B

A100 GPU 的峰值 FP32 算力约 19500 GFLOPS，全局内存带宽约 2000 GB/s，硬件本身能支撑的访存比上限约为 9.75 FLOP/B。

朴素版的 0.25 FLOP/B 远低于此阈值，意味着 GPU 大部分时间在等待内存数据，计算单元严重空转——这类程序称为"内存瓶颈（Memory-Bound）"程序。

三、分块优化核心思路

分块算法不改变线程分工，只改变数据访问方式：

1. 将数据从全局内存提前搬进共享内存；
2. 计算阶段全程只访问共享内存，不再触碰全局内存；
3. 块内所有线程协作加载，实现「一次加载，全块共用」。

关键在于第 3 点。以 TILE_WIDTH=16 为例，块内 256 个线程共同加载一个 16×16 的 M 瓦片和一个 16×16 的 N 瓦片进共享内存，之后每个线程各自用这份共享数据完成 16 次乘加——原本每人各自加载 16 次全局内存，现在 256 个线程共享 1 次加载，全局内存访问量减少 16 倍（理论值）。

这种技术称为 Strip-Mining：将长向量点积拆分为多个等长的小段，每段对应一个瓦片，分阶段完成。

四、分块乘法完整代码

前提：Width 必须是 TILE_WIDTH 的整数倍。当前代码不含边界检查，生产代码请自行添加。

#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>

#define TILE_WIDTH 16

// 分块矩阵乘法核函数
__global__ void MatrixMulTiled(float *M, float *N, float *P, int Width) {

    // 共享内存：存放当前瓦片的 M、N 分块
    __shared__ float Mds[TILE_WIDTH][TILE_WIDTH];
    __shared__ float Nds[TILE_WIDTH][TILE_WIDTH];

    int tx = threadIdx.x;
    int ty = threadIdx.y;

    // 当前线程负责的 P 矩阵坐标（整个核函数运行期间不变）
    int Row = blockIdx.y * TILE_WIDTH + ty;
    int Col = blockIdx.x * TILE_WIDTH + tx;

    // 私有浮点累加器，跨所有瓦片轮次持续累加
    float Pvalue = 0.0f;

    // Strip-Mining：将长点积拆分为 Width/TILE_WIDTH 个瓦片轮次
    for (int ph = 0; ph < Width / TILE_WIDTH; ph++) {

        // 每线程从全局内存各取 1 个 M、1 个 N 元素写入共享内存
        Mds[ty][tx] = M[Row * Width + ph * TILE_WIDTH + tx];
        Nds[ty][tx] = N[(ph * TILE_WIDTH + ty) * Width + Col];

        // ① 同步：等待块内所有线程完成加载（防止 RAW）
        __syncthreads();

        // 核心计算：全程访问共享内存，不碰全局内存
        for (int k = 0; k < TILE_WIDTH; k++) {
            Pvalue += Mds[ty][k] * Nds[k][tx];
        }

        // ② 同步：等待所有线程用完本轮瓦片（防止 WAR）
        __syncthreads();
    }

    // 所有轮次累加完成，写回最终结果
    P[Row * Width + Col] = Pvalue;
}

五、逐模块详解

5.1 关键参数含义

5.2 线程到矩阵元素的映射

int Row = blockIdx.y * TILE_WIDTH + ty;
int Col = blockIdx.x * TILE_WIDTH + tx;

每个线程通过「块索引 × 块大小 + 块内偏移」唯一定位其负责的 P[Row][Col]，整个核函数运行期间该映射关系不变。

5.3 共享内存加载：每线程各取一个元素

Mds[ty][tx] = M[Row * Width + ph * TILE_WIDTH + tx];
Nds[ty][tx] = N[(ph * TILE_WIDTH + ty) * Width + Col];

每个线程写入 1 个 M 元素 + 1 个 N 元素。块内 TILE_WIDTH² 个线程协作，刚好填满两个 TILE_WIDTH × TILE_WIDTH 的瓦片。每个线程的 (tx, ty) 唯一，写入地址互不重叠，不存在冲突。

5.3.1 具体示例：从线程到全局地址的映射

下面通过一个具体例子来展示上述索引公式的实际计算过程。
假设 TILE_WIDTH = 16，线程块索引 (blockIdx.x, blockIdx.y) = (1, 2)，线程在块内的索引 (tx, ty) = (3, 5)。

• 第一步：确定该线程负责的输出元素**

Row = blockIdx.y * TILE_WIDTH + ty = 2 * 16 + 5 = 37
Col = blockIdx.x * TILE_WIDTH + tx = 1 * 16 + 3 = 19

即该线程负责计算 P[37][19]。

• 第二步：在不同轮次 ph 下，加载 M 和 N 的哪些全局元素？**

为什么这样设计能保证正确性？

每个线程计算 P[Row][Col] 的公式为：

分块后等价于：

各轮次使用的矩阵区域对应关系：

通过这种方式，长向量点积被准确拆分为多个独立瓦片的点积之和，而每个瓦片的数据由整个线程块协作加载至共享内存，大幅减少了全局内存访问。

一个瓦片轮次只完成了点积的一段计算（TILE_WIDTH 项），并非最终结果。全部 Width/TILE_WIDTH 轮的分段结果累加在 Pvalue 中，最终写入 P。矩阵元素是 float，乘加结果自然也需要用 float 存储。

第 0 轮：Pvalue += Mds[ty][0..15] · Nds[0..15][tx]  （前 16 项）
第 1 轮：Pvalue += Mds[ty][0..15] · Nds[0..15][tx]  （中 16 项）
 ……
第 ph 轮结束后 Pvalue = P[Row][Col] 完整点积结果

5.5 单瓦片内共享内存访问次数

每个线程内存访问次数为：

5.6 边界检查的缺失（重要）

⚠️ 当前分块代码假设 Width 为 TILE_WIDTH 的整数倍，没有边界检查。若 Width 不满足此条件，全局内存访问会越界，行为未定义。生产代码需在加载和写回时加入如下判断：

if (Row < Width && Col < Width)
    Mds[ty][tx] = M[Row * Width + ph * TILE_WIDTH + tx];

六、两处同步屏障

__syncthreads() 是块内线程同步原语：确保块内所有线程都到达该语句后，才允许任意线程继续执行。分块乘法中用了两次，各有其使命。

第一处：加载之后，计算之前

Mds[ty][tx] = ...;
Nds[ty][tx] = ...;
__syncthreads();   // ← 此处
for (int k = 0; k < TILE_WIDTH; k++) { ... }

依赖类型： 写后读（Read After Write，RAW），又称真依赖。

若不加： 执行快的线程可能在其他线程还未写完共享内存时就开始读取，读到的是未初始化的脏数据，乘加结果错误。

作用： 保证所有线程都完成加载（写完）之后，再让任何线程开始计算（读取）。

第二处：计算之后，下一轮加载之前

for (int k = 0; k < TILE_WIDTH; k++) { ... }
__syncthreads();   // ← 此处
// 下一轮 ph 循环开始，覆盖 Mds、Nds

依赖类型： 读后写（Write After Read，WAR），又称假依赖。

若不加： 下一轮快速线程已经开始往共享内存写入新瓦片的数据，而慢速线程还在读当前轮次的旧数据，旧数据被覆盖，计算结果被污染。

作用： 保证所有线程都用完当前瓦片（读完）之后，再允许任何线程加载下一个瓦片（写入）。

小结

两处缺一不可，共同构成分块算法的数据安全边界。

七、性能本质分析

7.1 访存比对比

分块版访存比提升 TILE_WIDTH 倍（理论值）。TILE_WIDTH 越大，提升越多，但受限于每个 SM 的共享内存总容量（A100 为 164 KB/SM）。

⚠️ 手写分块通常只能达到硬件峰值的 20%–35%。生产环境建议直接使用 cuBLAS 或 CUTLASS，它们通过 warp-level 指令、双缓冲、向量化等高级技巧逼近极限性能。

八、CPU 缓存 vs GPU 共享内存

分块优化并非 GPU 独有，CPU 也有类似的 cache blocking 技术，但实现机制有本质区别：

GPU 共享内存存在的核心意义在于：当片上缓存面对海量并发线程时，硬件自动管理不可靠，必须由程序员显式控制数据复用，才能保证计算效率。

参考资料：Programming Massively Parallel Processors (PMPP), 4th Edition — David Kirk & Wen-mei Hwu