第八课：GPU 如何完成矩阵乘法

导语：这一课用矩阵乘法把前面学过的线程、Block、Shared Memory、寄存器和计算强度串联起来。

矩阵乘法是理解 GPU 架构最典型的例子，因为它同时涉及：

大量线程
Block 与 SM
合并访存
Shared Memory
寄存器
数据复用
Tensor Core

深度学习中的全连接层、卷积和 Transformer，大量计算最终都会转化为矩阵乘法。

---

一、矩阵乘法在算什么

假设：

$$A\in \mathbb{R}^{M\times K}$$ $$B\in \mathbb{R}^{K\times N}$$

计算：

$$C=A\times B$$

结果矩阵：

$$C\in \mathbb{R}^{M\times N}$$

其中每个元素为：

$$C_{i,j}=\sum_{k=0}^{K-1}A_{i,k}B_{k,j}$$

例如：

A 的第 i 行
    ×
B 的第 j 列
    ↓
得到 C[i][j]

每个输出元素需要完成 $K$ 次乘法和约 $K$ 次加法。

---

二、最直接的线程映射方式

最简单的 GPU 实现是：

一个 Thread 计算 C 中的一个元素。

例如：

图示：二、最直接的线程映射方式

每个线程执行：

float sum = 0.0f;

for (int k = 0; k < K; ++k) {
    sum += A[row * K + k] * B[k * N + col];
}

C[row * N + col] = sum;

这里：

• row 决定输出矩阵的行；
• col 决定输出矩阵的列；
• sum 保存在当前线程的寄存器中。

---

三、使用二维 Grid 和 Block

矩阵本身是二维结构，所以通常使用二维线程组织。

例如：

dim3 block(16, 16);
dim3 grid(
    (N + block.x - 1) / block.x,
    (M + block.y - 1) / block.y
);

一个 Block 有：

$$16\times16=256$$

个线程。

每个 Block 负责计算 C 中一个 $16\times16$ 的区域。

线程计算自己的全局坐标：

int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;

对应关系为：

图示：三、使用二维 Grid 和 Block

---

四、朴素矩阵乘法 Kernel

一个简单实现如下：

__global__ void matmul_naive(
    const float* A,
    const float* B,
    float* C,
    int M,
    int N,
    int K)
{
    int row =
        blockIdx.y * blockDim.y + threadIdx.y;

    int col =
        blockIdx.x * blockDim.x + threadIdx.x;

    if (row < M && col < N) {
        float sum = 0.0f;

        for (int k = 0; k < K; ++k) {
            sum += A[row * K + k]
                 * B[k * N + col];
        }

        C[row * N + col] = sum;
    }
}

它在功能上是正确的，但性能通常不够好。

问题主要不是计算方式，而是：

同样的数据被不同线程反复从 Global Memory 读取。

---

五、朴素实现中的重复访存

假设一个 Block 计算 $16\times16$ 个输出元素。

Block 中同一行的 16 个线程：

Thread (row, 0)
Thread (row, 1)
...
Thread (row, 15)

它们计算 C 同一行、不同列的元素。

在某个固定的 $k$ 下，它们都需要：

$$A[row][k]$$

也就是说，一个 A 元素可能被多个线程重复使用。

类似地，同一列方向的线程会重复使用 B 中的元素。

但朴素实现可能让每个线程分别发出自己的加载请求：

Thread 0 读取 A[row][k]
Thread 1 读取 A[row][k]
Thread 2 读取 A[row][k]
...

缓存可能缓解一部分重复访问，但程序没有显式保证高效复用。

---

六、矩阵乘法的核心：数据复用

矩阵乘法适合 GPU，不只是因为输出元素多，还因为数据可以被大量复用。

一个 A 元素：

$$A[i][k]$$

会参与 C 第 $i$ 行中的多个输出元素计算。

一个 B 元素：

$$B[k][j]$$

会参与 C 第 $j$ 列中的多个输出元素计算。

所以理想方式是：

从 Global Memory 读取一次
           ↓
放进片上存储
           ↓
被多个线程重复使用

这就是矩阵分块，也就是 Tiling。

---

七、什么是矩阵分块 Tiling

把大矩阵划分为多个小矩阵块：

图示：七、什么是矩阵分块 Tiling

假设 Tile 大小为：

$$16\times16$$

一个 Block 计算 C 中一个 $16\times16$ 的输出 Tile。

为了计算这个输出 Tile，需要依次加载：

• A 的一个 $16\times16$ Tile；
• B 的一个 $16\times16$ Tile。

然后让 Block 中的线程重复使用这些数据。

---

八、为什么要使用 Shared Memory

Shared Memory 是 Block 内线程共享的片上存储器。

对于一个输出 Tile，可以执行：

第一步：
把 A 的 Tile 从 Global Memory 搬到 Shared Memory

第二步：
把 B 的 Tile从 Global Memory 搬到 Shared Memory

第三步：
Block 内所有线程使用这些数据进行乘加

第四步：
继续加载下一组 Tile

数据路径变成：

图示：八、为什么要使用 Shared Memory

相对于每次乘法都访问显存，这种方式显著减少了 Global Memory 流量。

---

九、一个 Tile 如何计算

假设 Tile 大小为 16。

一个 Block 负责：

$$C_{\text{tile}}$$

计算公式可以写成：

$$C_{\text{tile}} = \sum_t A_{\text{tile},t} \times B_{\text{tile},t}$$

例如 K 维度为 64，则需要分成 4 个阶段：

图示：九、一个 Tile 如何计算

每个阶段都会对最终结果做部分累加。

---

十、Tiled Matrix Multiplication 的代码

简化实现如下：

#define TILE 16

__global__ void matmul_tiled(
    const float* A,
    const float* B,
    float* C,
    int M,
    int N,
    int K)
{
    __shared__ float As[TILE][TILE];
    __shared__ float Bs[TILE][TILE];

    int row =
        blockIdx.y * TILE + threadIdx.y;

    int col =
        blockIdx.x * TILE + threadIdx.x;

    float sum = 0.0f;

    int num_tiles = (K + TILE - 1) / TILE;

    for (int t = 0; t < num_tiles; ++t) {
        int a_col = t * TILE + threadIdx.x;
        int b_row = t * TILE + threadIdx.y;

        if (row < M && a_col < K) {
            As[threadIdx.y][threadIdx.x] =
                A[row * K + a_col];
        } else {
            As[threadIdx.y][threadIdx.x] = 0.0f;
        }

        if (b_row < K && col < N) {
            Bs[threadIdx.y][threadIdx.x] =
                B[b_row * N + col];
        } else {
            Bs[threadIdx.y][threadIdx.x] = 0.0f;
        }

        __syncthreads();

        for (int k = 0; k < TILE; ++k) {
            sum += As[threadIdx.y][k]
                 * Bs[k][threadIdx.x];
        }

        __syncthreads();
    }

    if (row < M && col < N) {
        C[row * N + col] = sum;
    }
}

---

十一、第一次同步在保证什么

线程共同加载 Shared Memory：

As[threadIdx.y][threadIdx.x] = ...;
Bs[threadIdx.y][threadIdx.x] = ...;

每个线程只负责加载部分元素。

但后面的计算会读取其他线程加载的数据：

sum += As[threadIdx.y][k]
     * Bs[k][threadIdx.x];

因此必须保证：

整个 Tile 已经加载完成，才能开始计算。

所以需要：

__syncthreads();

执行过程：

图示：十一、第一次同步在保证什么

---

十二、第二次同步在保证什么

计算结束后，下一轮循环会覆盖：

As
Bs

但不同线程完成当前 Tile 计算的时间可能不完全一致。

如果某些线程提前覆盖 Shared Memory，其他线程可能仍在读取旧 Tile。

所以第二个：

__syncthreads();

保证所有线程都完成当前 Tile 的使用后，才开始加载下一 Tile。

完整阶段是：

加载 Tile
    ↓
同步
    ↓
使用 Tile 计算
    ↓
同步
    ↓
覆盖为下一 Tile

---

十三、Shared Memory 如何提高数据复用

一个 $16\times16$ 的 A Tile 有：

$$256$$

个元素。

加载到 Shared Memory 后，其中每个 A 元素可能被同一输出行的 16 个线程使用。

B Tile 也是类似。

因此：

每个元素从显存读取一次
        ↓
在 Shared Memory 中被使用多次

一个阶段中，Block 加载：

$$2\times16\times16=512$$

个 FP32 数据。

但会完成：

$$16\times16\times16=4096$$

次乘法，以及相应的加法。

这显著提高了计算强度。

---

十四、为什么累加结果保存在寄存器

每个线程计算自己的输出元素：

float sum = 0.0f;

sum 通常保存在寄存器中。

在整个 K 循环中：

sum += ...

都不需要把中间结果写回显存。

最终只执行一次：

C[row * N + col] = sum;

数据流是：

读取 A、B
    ↓
乘法
    ↓
结果在寄存器 sum 中累加
    ↓
全部计算完成
    ↓
一次写回 C

这比每次乘加都访问 Global Memory 高效得多。

---

十五、合并访存在这个 Kernel 中如何体现

加载 A 时：

A[row * K + a_col]

同一行中相邻线程的 threadIdx.x 连续，因此：

Thread 0 → A[row][0]
Thread 1 → A[row][1]
Thread 2 → A[row][2]
...

地址连续，容易形成合并访存。

加载 B 时：

B[b_row * N + col]

同一行中相邻线程的 col 连续：

Thread 0 → B[b_row][col]
Thread 1 → B[b_row][col+1]
...

也容易形成合并访存。

所以优秀的 Tiling 设计同时实现：

Global Memory 合并访存
+
Shared Memory 数据复用
+
寄存器累加

---

十六、为什么不能把整个矩阵放进 Shared Memory

Shared Memory 容量有限。

例如一个大矩阵可能有几百 MB，甚至几 GB，而一个 SM 的 Shared Memory 只有有限的片上容量。

因此只能分块：

大矩阵
    ↓
每次搬一小块
    ↓
计算完后换下一块

这也是计算机存储层次中常见的思想：

用小容量高速存储器反复缓存正在使用的数据。

---

十七、Tile 大小是不是越大越好

不是。

增大 Tile 可以增加数据复用，但也会消耗更多资源。

假设 Tile 为 $T\times T$。

Shared Memory 大约需要存放：

$$2T^2$$

个元素。

如果使用 FP32，所需空间约为：

$$2T^2\times4\text{ Bytes}$$

例如：

TILE = 16

$$2\times16^2\times4 =2048\text{ Bytes}$$

TILE = 32

$$2\times32^2\times4 =8192\text{ Bytes}$$

这还只是最简单的两个 Tile，不包括双缓冲、Padding 和其他数据。

Tile 增大还可能导致：

• Block 线程数量增加；
• 每线程寄存器使用增加；
• Shared Memory 使用增加；
• 同时驻留的 Block 数量下降；
• Occupancy 下降。

因此 Tile 大小需要在数据复用和并发度之间权衡。

---

十八、为什么每个线程不只计算一个输出元素

前面的实现中：

一个线程 → 一个 C 元素

更高性能的实现常让一个线程计算多个输出元素，例如：

一个线程 → 2×2 输出微块

线程内部维护多个累加器：

float c00 = 0.0f;
float c01 = 0.0f;
float c10 = 0.0f;
float c11 = 0.0f;

这样从 Shared Memory 读取一个 A 或 B 数据后，可以用于多个乘加。

例如：

读取一个 A 值
    ↓
分别与多个 B 值相乘
    ↓
更新多个输出累加器

这称为 Register Tiling 或 Thread-Level Tiling。

---

十九、现代 GEMM 的多级分块

高性能矩阵乘法通常不只是简单的 Block Tile。

它可能包含多个层次：

图示：十九、现代 GEMM 的多级分块

例如：

一个 Block 计算 128×128 输出 Tile

Block 内：
Warp 0 计算其中一个区域
Warp 1 计算另一个区域
Warp 2 ...

每个 Warp 再划分：

一个 Warp 计算 32×64 Tile

每个线程维护更小的输出微块：

一个 Thread 计算若干输出元素

这种分层结构能够充分利用：

• Global Memory 带宽；
• Shared Memory 带宽；
• 寄存器；
• Warp 调度；
• Tensor Core。

---

二十、普通 CUDA Core 如何做矩阵乘法

使用普通 CUDA Core 时，核心指令通常是 FMA：

$$c=a\times b+c$$

每个线程不断执行：

sum = fma(a, b, sum);

对一个 Warp 而言：

32 个线程
各自完成标量 FMA

每个线程有自己的：

• A 操作数；
• B 操作数；
• 累加器；
• 输出元素。

---

二十一、Tensor Core 如何做矩阵乘法

Tensor Core 不只是执行一个标量乘加，而是执行小矩阵块的乘加：

$$D=A\times B+C$$

这类指令常称为：

MMA：Matrix Multiply-Accumulate

与普通 CUDA Core 相比：

CUDA Core：
执行许多标量 FMA

Tensor Core：
一次处理一个小矩阵块

但 Tensor Core 的操作通常由整个 Warp 协作完成，而不是简单地让一个线程独立计算完整小矩阵。

一个 Warp 中的线程共同持有矩阵片段：

图示：二十一、Tensor Core 如何做矩阵乘法

---

二十二、Tensor Core 不是自动使用的

代码中出现矩阵乘法，不代表一定会使用 Tensor Core。

需要满足一些条件，例如：

• 使用支持的数据类型；
• 矩阵维度和布局合适；
• 编译器或库生成 MMA 指令；
• 数据对齐满足要求；
• 使用对应的 CUDA 库或接口。

常见支持格式包括：

FP16
BF16
TF32
FP8
INT8
FP4

具体支持情况取决于 GPU 架构。

实际开发中通常通过：

• cuBLAS；
• cuBLASLt；
• CUTLASS；
• 深度学习框架；

来使用 Tensor Core，而不是手工编写所有底层指令。

---

二十三、为什么 Tensor Core 仍然需要复杂的数据搬运

即使 Tensor Core 算得非常快，也必须不断获得数据。

如果数据供应不足：

Tensor Core 计算很快
        ↓
很快把当前数据算完
        ↓
等待下一批数据

因此高性能 Tensor Core Kernel 必须精心设计：

图示：二十三、为什么 Tensor Core 仍然需要复杂的数据搬运

Tensor Core 越快，数据供应问题反而越明显。

这也是现代 GPU 设计会重点强化：

• HBM 带宽；
• L2 Cache；
• Shared Memory；
• 异步数据搬运；
• Tensor Memory Accelerator；

等结构的原因。

---

二十四、什么是双缓冲

简单 Tiling 的过程是：

加载 Tile 0
等待加载完成
计算 Tile 0

加载 Tile 1
等待加载完成
计算 Tile 1

加载和计算大部分时间是串行的。

更高级的实现会使用两组缓冲区：

Buffer 0：当前正在计算
Buffer 1：提前加载下一 Tile

执行过程：

计算 Tile 0 的同时
加载 Tile 1

计算 Tile 1 的同时
加载 Tile 2

这叫 Double Buffering，双缓冲。

目的就是让：

数据搬运
+
矩阵计算

尽可能重叠。

---

二十五、异步拷贝的作用

现代 GPU 可以提供更高效的数据搬运方式，让 Global Memory 到 Shared Memory 的传输与计算重叠。

简化理解：

传统方式：
线程发出加载
等待数据返回
写入 Shared Memory

异步方式：
启动搬运任务
继续执行已有数据的计算
之后等待新数据准备好

这有助于减少：

计算单元等待数据

但异步拷贝仍需要：

• 正确的流水阶段；
• 缓冲区管理；
• 同步；
• 数据依赖控制。

---

二十六、矩阵边界如何处理

如果 M、N、K 不是 Tile 大小的整数倍，最后一块可能不完整。

例如：

$$M=1000$$

Tile 大小为 16，则最后一组线程可能对应：

行 992～1007

但只有 992～999 有效。

因此加载时需要边界判断：

if (row < M && a_col < K) {
    As[...] = A[...];
} else {
    As[...] = 0.0f;
}

写回时也要判断：

if (row < M && col < N) {
    C[...] = sum;
}

这些边界分支通常只影响边缘 Block，整体影响往往不大。

---

二十七、Shared Memory Bank Conflict 是否会出现

在内层计算中：

As[threadIdx.y][k]
Bs[k][threadIdx.x]

不同线程的访问模式需要仔细分析。

常见情况包括：

• 多个线程读取同一个 A 值，可以通过广播高效处理；
• 不同线程读取 B 的连续元素，通常映射到不同 Bank；
• 矩阵转置或特殊步长访问可能产生 Bank Conflict。

有时会给 Shared Memory 数组增加 Padding：

__shared__ float tile[32][33];

而不是：

__shared__ float tile[32][32];

多出的 1 列可以改变地址到 Bank 的映射，避免某些转置访问中的严重冲突。

---

二十八、矩阵布局为什么重要

矩阵通常有两种主要存储方式。

行优先 Row-Major

A[0][0], A[0][1], A[0][2], ...
A[1][0], A[1][1], A[1][2], ...

同一行元素连续。

列优先 Column-Major

A[0][0], A[1][0], A[2][0], ...
A[0][1], A[1][1], A[2][1], ...

同一列元素连续。

线程访问方向如果和内存布局匹配，更容易形成合并访存。

矩阵库通常允许指定：

A 是否转置
B 是否转置
矩阵是行优先还是列优先

但内部实现可能采用更复杂的打包或重排布局。

---

二十九、朴素实现和分块实现的比较

对比项	朴素实现	分块实现
Thread 工作	一个输出元素	一个或多个输出元素
A/B 数据来源	主要直接访问 Global Memory	先加载到 Shared Memory
数据复用	较少	较多
Global Memory 流量	较大	较小
同步	通常不需要	每个 Tile 阶段需要
实现复杂度	低	较高
性能	通常较低	通常较高

分块实现增加了：

• Shared Memory；
• 同步；
• 地址控制；
• 边界处理。

但换来了更高的数据复用率和计算强度。

---

三十、为什么矩阵乘法特别适合 GPU

矩阵乘法具备几个重要特征：

第一，输出元素数量很多，可以产生大量并行线程。

第二，各个输出元素基本独立，容易分配给不同 Thread、Warp 和 Block。

第三，A 和 B 中的数据可以重复使用，适合 Shared Memory 和寄存器分块。

第四，核心计算是规则的乘加操作，分支很少。

第五，内存访问模式规则，容易形成合并访存。

第六，可以使用低精度和 Tensor Core 获得极高吞吐率。

可以概括为：

并行度高
+
控制流规则
+
数据复用高
+
乘加密集
=
非常适合 GPU

---

三十一、从架构角度看完整数据流

一个高性能矩阵乘法大致经历：

图示：三十一、从架构角度看完整数据流

每一层都在解决不同问题：

层次	主要作用
Global Memory	保存完整矩阵
L2 Cache	缓解重复片外访问
Shared Memory	Block 级数据复用
Register	Thread/Warp 级数据复用与累加
CUDA/Tensor Core	完成实际乘加计算

---

本课核心结论

第一，最简单的矩阵乘法映射是一个线程计算一个输出元素。

第二，一个二维 Block 通常负责 C 中的一个输出 Tile。

第三，朴素实现会让相同的 A、B 数据被反复从 Global Memory 读取。

第四，Tiling 将 A 和 B 的小块加载到 Shared Memory，供 Block 内线程重复使用。

第五，__syncthreads() 保证 Tile 加载和使用阶段的正确性。

第六，每个线程的部分和通常保存在寄存器中，最后只写回一次。

第七，合并访存、Shared Memory 复用和寄存器累加共同提高矩阵乘法性能。

第八，Tile 不是越大越好，因为会消耗线程、Shared Memory 和寄存器资源。

第九，高性能实现还会使用 Warp Tiling、Thread Tiling 和双缓冲。

第十，普通 CUDA Core 执行标量 FMA，Tensor Core 执行小矩阵 MMA。

第十一，Tensor Core 算得很快，但依然依赖高效的数据搬运和分块。

第十二，现代高性能 GEMM 的本质是多级存储层次上的数据复用。

下一课适合学习：

第九课：GPU 如何执行神经网络卷积

主要包括：

卷积为什么可以转化为矩阵乘法
输入特征图、权重和输出特征图如何映射
Channel、Height、Width 分别是什么
直接卷积与 im2col
Weight Stationary、Output Stationary
Tensor Core 如何加速卷积
为什么深度可分离卷积不一定容易跑满 GPU