GEMM/GEMV 优化
GEMM/GEMV 优化核心原则
- GEMV(矩阵向量乘法)是典型的内存瓶颈负载,优化重点是带宽利用率
- GEMM(矩阵矩阵乘法)是计算瓶颈负载,优化重点是计算单元利用率
- 分块(Tiling)是核心策略,平衡容量与计算密度
- 充分利用张量核心(Tensor Core)等专用硬件单元
阅读本文时,建议先判断算子属于 GEMV 还是 GEMM,再选择对应优化路径。GEMV 优先检查访存是否合并;GEMM 优先检查分块复用、寄存器累加和张量核心使用。
| 场景 | 主要瓶颈 | 优先检查项 |
|---|---|---|
| GEMV | 内存带宽 | 合并访存、向量化加载、向量缓存 |
| GEMM | 计算吞吐 | 分块复用、寄存器累加、张量核心 |
基础示例:从基础实现(Naive)到共享内存优化
下面的示例先给出直接实现,再给出共享内存分块实现。两者计算同一个矩阵乘法,差异在于是否把 A/B 的局部数据缓存到共享内存中复用。
- 基础实现
- 共享内存分块
基础实现(Naive,无共享内存)
每个线程直接从全局内存读取数据,计算 C 矩阵的一个元素:
// 矩阵存储结构
typedef struct {
int width;
int height;
float* elements;
} Matrix;
#define BLOCK_SIZE 16
__global__ void MatMulNaive(Matrix A, Matrix B, Matrix C)
{
// 每个线程计算 C 的一个元素
float Cvalue = 0;
int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;
if (row < C.height && col < C.width) {
for (int e = 0; e < A.width; ++e) {
Cvalue += A.elements[row * A.width + e] * B.elements[e * B.width + col];
}
C.elements[row * C.width + col] = Cvalue;
}
}
// 启动配置
dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);
dim3 dimGrid((B.width + dimBlock.x - 1) / dimBlock.x,
(A.height + dimBlock.y - 1) / dimBlock.y);
MatMulNaive<<<dimGrid, dimBlock>>>(A, B, C);
查看问题分析
问题 分析:
| 矩阵 | 全局内存读取次数 | 原因 |
|---|---|---|
| A | B.width 次 | 每个线程独立读取 A 的行元素 |
| B | A.height 次 | 每个线程独立读取 B 的列元素 |
示例:对于 1024×1024 矩阵,BLOCK_SIZE=16
- 全局内存读取:约 2 × 1024 × 1024 × 16 = 33M 次
- 带宽利用率极低
共享内存优化实现
使用分块(Tiling)策略,将数据加载到共享内存后重复使用:
// 增加 stride 字段以支持子矩阵
typedef struct {
int width;
int height;
int stride;
float* elements;
} Matrix;
#define BLOCK_SIZE 16
// 获取子矩阵
__device__ Matrix GetSubMatrix(Matrix M, int row, int col)
{
Matrix Msub;
Msub.width = BLOCK_SIZE;
Msub.height = BLOCK_SIZE;
Msub.stride = M.stride;
Msub.elements = &M.elements[M.stride * BLOCK_SIZE * row + BLOCK_SIZE * col];
return Msub;
}
__global__ void MatMulTiled(Matrix A, Matrix B, Matrix C)
{
int blockRow = blockIdx.y;
int blockCol = blockIdx.x;
// 每个线程块计算 C 的一个子矩阵
Matrix Csub = GetSubMatrix(C, blockRow, blockCol);
float Cvalue = 0;
int row = threadIdx.y;
int col = threadIdx.x;
__shared__ float As[BLOCK_SIZE][BLOCK_SIZE];
__shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE];
// 分块迭代
for (int m = 0; m < (A.width + BLOCK_SIZE - 1) / BLOCK_SIZE; ++m) {
// 获取 A、B 的子矩阵
Matrix Asub = GetSubMatrix(A, blockRow, m);
Matrix Bsub = GetSubMatrix(B, m, blockCol);
// 协作加载:每个线程加载一个元素
int a_global_row = blockRow * BLOCK_SIZE + row;
int a_global_col = m * BLOCK_SIZE + col;
int b_global_row = m * BLOCK_SIZE + row;
int b_global_col = blockCol * BLOCK_SIZE + col;
As[row][col] = (a_global_row < A.height && a_global_col < A.width)
? Asub.elements[row * Asub.stride + col]
: 0.0f;
Bs[row][col] = (b_global_row < B.height && b_global_col < B.width)
? Bsub.elements[row * Bsub.stride + col]
: 0.0f;
// 同步:确保所有数据加载完成
__syncthreads();
// 在共享内存中计算
for (int e = 0; e < BLOCK_SIZE; ++e) {
Cvalue += As[row][e] * Bs[e][col];
}
__syncthreads(); // 确保计算完成后再加载下一块数据
}
int c_global_row = blockRow * BLOCK_SIZE + row;
int c_global_col = blockCol * BLOCK_SIZE + col;
if (c_global_row < C.height && c_global_col < C.width) {
Csub.elements[row * Csub.stride + col] = Cvalue;
}
}
查看优化效果和关键设计
优化效果:
| 矩阵 | 全局内存读取次数 | 降低倍数 |
|---|---|---|
| A | B.width / BLOCK_SIZE 次 | 16 倍 |
| B | A.height / BLOCK_SIZE 次 | 16 倍 |
关键设计:
- 分块策略:将大矩阵分解为 BLOCK_SIZE×BLOCK_SIZE 的小块
- 数据复用:每个块内元素被 BLOCK_SIZE 个线程重复使用
- 协作加载:每个线程只加载一个元素,避免重复
- 同步原语:
__syncthreads()确保数据一致性