CUDA 性能优化 Q&A


Q1:优化一个 CUDA Kernel 可以从哪些角度入手?

性能优化没有万能药,但有清晰的分析框架。一个 kernel 的性能瓶颈通常来自以下三个方向之一,找到瓶颈再针对性优化,效果远好于盲目调参。

优化层次全景

1
2
3
4
5
6
7
8
9
10
11
┌─────────────────────────────────────────────────┐
│ 算法层 选择计算复杂度更低的算法(最高回报) │
├─────────────────────────────────────────────────┤
│ 内存层 减少全局内存访问、提升合并、利用共享内存 │
├─────────────────────────────────────────────────┤
│ 并行层 提高 Occupancy、消除 Warp Divergence │
├─────────────────────────────────────────────────┤
│ 指令层 Loop Unroll、向量化加载、减少低效指令 │
├─────────────────────────────────────────────────┤
│ 调度层 多 Stream、CUDA Graph、算子融合 │
└─────────────────────────────────────────────────┘

确定瓶颈的第一步:Roofline 模型

用"算术强度"(Arithmetic Intensity = FLOP / Byte)判断 kernel 是计算受限还是内存带宽受限

  • 内存带宽受限(Memory-bound):算术强度低,每访问一字节只做少量运算。优化方向:减少全局内存访问次数,提升合并访问,用共享内存缓存热点数据。
  • 计算受限(Compute-bound):算术强度高,运算量大。优化方向:用 Tensor Core、FP16、ILP,提高指令级并行。

常用工具

  • Nsight Compute:深度分析单个 kernel 的寄存器使用、内存吞吐、Occupancy、Warp 效率等指标
  • Nsight Systems:分析整个应用的时间线,找到 CPU-GPU 交互的瓶颈
  • nvcc --resource-usage:编译期查看寄存器和共享内存用量

主要优化方向速查

类别 具体手段
内存访问 合并访问、共享内存 Tiling、向量化加载(float4)、只读缓存(__ldg
并行度 调整 blockDim 提高 Occupancy、消除 Warp Divergence、使用 grid-stride 循环
指令效率 #pragma unroll、使用 intrinsic 函数、避免低效的整数除法取模
数据精度 FP16/BF16 替换 FP32(在精度允许时)、Tensor Core 加速矩阵运算
计算/传输重叠 多 Stream 流水线、双缓冲(Double Buffering)、异步拷贝(LDGSTS/TMA)
启动开销 CUDA Graph 消除重复启动开销、算子融合减少 kernel 数量

Q2:什么是 Instruction-Level Parallelism(ILP)?如何利用?

什么是 ILP

指令级并行(ILP,Instruction-Level Parallelism)是指在单个线程内,多条独立的指令可以被流水线化或同时发射执行,而不必等待前一条指令的结果。

GPU 的 Warp 调度器通过在 Warp 间切换来隐藏内存延迟(线程级并行,TLP)。但当活跃 Warp 数有限(低 Occupancy),或者所有 Warp 都在等待同一资源时,TLP 不足以填满流水线。此时,ILP 可以作为补充手段——让单个线程内的独立操作互相隐藏延迟。

如何利用 ILP

方法一:让每个线程处理多个数据元素

最直接的 ILP 来源是让单个线程负责多个独立的运算,编译器可以将这些运算流水线化:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
// 低 ILP:每个线程处理 1 个元素
__global__ void kernel_low_ilp(float* a, float* b, float* c, int N) {
int i = threadIdx.x + blockDim.x * blockIdx.x;
c[i] = a[i] + b[i];
}

// 高 ILP:每个线程处理 4 个元素,4 个加法相互独立,可以并行发射
__global__ void kernel_high_ilp(float* a, float* b, float* c, int N) {
int i = (threadIdx.x + blockDim.x * blockIdx.x) * 4;
float r0 = a[i+0] + b[i+0]; // 这四条加法互不依赖
float r1 = a[i+1] + b[i+1];
float r2 = a[i+2] + b[i+2];
float r3 = a[i+3] + b[i+3];
c[i+0] = r0; c[i+1] = r1; c[i+2] = r2; c[i+3] = r3;
}

方法二:向量化加载(float4、int4)

float4 一次加载 128 位(4 个 float),既触发向量化的内存访问(更高效的内存事务),也在单条指令内提供了 4 路并行计算机会:

1
2
3
4
5
__global__ void vectorized(float4* a, float4* b, float4* c, int N) {
int i = threadIdx.x + blockDim.x * blockIdx.x;
float4 va = a[i], vb = b[i];
c[i] = {va.x+vb.x, va.y+vb.y, va.z+vb.z, va.w+vb.w};
}

方法三:Loop Unrolling(见 Q4)

展开循环使编译器能看到更多独立指令,从而重排调度,提高 ILP。

ILP 与 Occupancy 的权衡

每个线程处理更多元素意味着每个线程可能需要更多寄存器,这会降低 Occupancy(每 SM 能容纳的线程数变少)。ILP 和 Occupancy 需要协同调优——有时候低 Occupancy + 高 ILP 反而比高 Occupancy + 低 ILP 更快,需要通过实际 profiling 确认。


Q3:如何减少全局内存访问延迟(Latency Hiding)?

全局内存访问延迟高达数百个时钟周期。GPU 隐藏这一延迟的核心机制是多线程并发:当一个 Warp 等待内存时,调度器立刻切换到另一个就绪 Warp。但除了依赖 Occupancy,还有更多主动手段。

手段一:提高 Occupancy,给调度器更多选择

SM 上活跃的 Warp 越多,在等待内存的 Warp 越多时,调度器越有可能找到就绪的 Warp 来填充空档。优化方法见基础概念 Q8。

手段二:合并全局内存访问(最重要)

非合并访问不仅带宽利用率低(最差情况 12.5%),还会产生更多内存事务,每个事务都要承受完整的延迟。合并访问把 32 个线程的请求合并为少量 32 字节事务,延迟次数大幅减少。

1
2
3
4
5
// 合并访问:连续线程访问连续地址
data[threadIdx.x + blockDim.x * blockIdx.x]

// 用 __ldg() 利用只读纹理缓存,对只读数据有额外缓存加速
float val = __ldg(&data[i]);

手段三:共享内存缓存(减少对全局内存的重复访问)

对于会被多次读取的数据,先一次性加载到共享内存,之后所有访问从共享内存(低延迟)读取,彻底避免重复访问全局内存:

1
2
3
4
5
6
7
__global__ void kernel(float* global_data, int N) {
__shared__ float smem[BLOCK_SIZE];
// 一次性全局内存加载(合并)
smem[threadIdx.x] = global_data[blockDim.x * blockIdx.x + threadIdx.x];
__syncthreads();
// 后续多次使用 smem[...],无需再访问全局内存
}

手段四:异步拷贝(LDGSTS / TMA,CC ≥ 8.0)

传统的 shared[i] = global[i] 会占用寄存器作为中转(全局内存 → 寄存器 → 共享内存),且必须等待完成才能继续。

CC 8.0 引入的 LDGSTS(Load-Global-Store-Shared)可以直接从全局内存异步写入共享内存,绕过寄存器,线程在此期间可以继续执行计算。这是双缓冲(Double Buffering)技术的硬件基础。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#include <cuda/pipeline>
__global__ void async_copy_kernel(int* global_out, int const* global_in, int N) {
__shared__ int smem[BLOCK];
auto pipeline = cuda::make_pipeline();

for (int i = 0; i < N; i += BLOCK) {
// 异步发起拷贝,线程不阻塞
cuda::memcpy_async(smem, global_in + i, sizeof(int)*BLOCK, pipeline);
pipeline.producer_commit();
// 等待异步拷贝完成
pipeline.consumer_wait();
__syncthreads();
compute(global_out, smem);
pipeline.consumer_release();
}
}

手段五:L2 缓存持久化(CC ≥ 8.0)

对于反复访问的热点数据(如大型权重矩阵),可以设置 L2 持久化窗口,让这部分数据在 L2 中高优先级驻留,减少对 DRAM 的访问次数。


Q4:什么是 Loop Unrolling?在 CUDA 中如何使用 #pragma unroll

什么是 Loop Unrolling

循环展开(Loop Unrolling)是一种编译器优化技术,将循环体复制若干次,减少循环控制开销(分支判断、计数器更新),并暴露出更多独立指令,给编译器提供更大的指令调度空间,从而提高 ILP 和流水线利用率。

展开前:

1
2
3
for (int i = 0; i < 4; i++) {
c[i] = a[i] + b[i]; // 每次迭代:加载a、加载b、做加法、存储c、更新i、判断i<4
}

展开后(4 路展开):

1
2
3
4
5
c[0] = a[0] + b[0];   // 4 对加法互相独立
c[1] = a[1] + b[1]; // 编译器可以重排,隐藏内存延迟
c[2] = a[2] + b[2];
c[3] = a[3] + b[3];
// 消除了循环控制指令(分支、计数器)

CUDA 中的 #pragma unroll

CUDA 编译器(NVCC)对行程数(trip count)固定的小循环会自动展开,但可以用 #pragma unroll 显式控制展开行为:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
__device__ void example(int* p1, int* p2) {

// 完全展开:不带参数,循环必须有编译期确定的行程数
#pragma unroll
for (int i = 0; i < 12; i++)
p1[i] += p2[i] * 2;

// 部分展开:展开 4 次(参数必须是编译期常量表达式)
constexpr int UNROLL = 4;
#pragma unroll UNROLL
for (int i = 0; i < 12; i++)
p1[i] += p2[i] * 4;

// 禁用展开:参数为 1
#pragma unroll 1
for (int i = 0; i < 12; i++)
p1[i] += p2[i] * 8;
}

在 GEMM Tiling 中的典型使用

矩阵乘法的内循环(K 方向的累加)是 #pragma unroll 最常见的使用场景:

1
2
3
4
5
6
7
8
__global__ void gemm(float* A, float* B, float* C, int K) {
float acc = 0.0f;
#pragma unroll 8 // 展开 8 次,让编译器流水线化 8 个 FMA 操作
for (int k = 0; k < K; k++) {
acc += A[row * K + k] * B[k * N + col];
}
C[row * N + col] = acc;
}

使用注意事项

  • #pragma unroll 只作用于紧跟其后的那一个循环
  • 展开过度会增加寄存器压力,可能降低 Occupancy——当 K 很大时,完全展开不现实。
  • 展开因子通常取 4、8、16,具体值需 profiling 确定。
  • 若行程数不是展开因子的整数倍,编译器会自动处理余数循环。

Q5:向量化加载(如 float4int4)的意义是什么?

问题背景

每次全局内存访问产生一个内存事务(至少 32 字节)。如果每个线程只访问 4 字节(float),则 32 个线程恰好用完 128 字节(4 × 32 = 128),需要 4 次 32 字节事务。这已经是合并访问的最佳情况。

但这里还有进一步提升空间:用 float4 类型,每个线程一次读取 16 字节(4 × 4B),32 个线程共读取 512 字节,需要 16 次 32 字节事务。乍看事务数增加了,但好处是:

  1. 减少了指令数量:4 个独立的 float 加载需要 4 条内存指令,float4 一次加载只需 1 条。
  2. 内存系统效率更高:内存控制器更喜欢连续的大块读取,减少地址解码开销。
  3. 线程可以处理更多元素:每个线程一次处理 4 个数据,天然提升 ILP。

使用方法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
// 普通加载:每次 4 字节,32 个线程产生 128 字节,4 次事务
__global__ void scalar_add(float* a, float* b, float* c, int N) {
int i = threadIdx.x + blockDim.x * blockIdx.x;
c[i] = a[i] + b[i];
}

// 向量化加载:每次 16 字节,32 个线程处理 4×N 的数据
__global__ void vector_add(float4* a, float4* b, float4* c, int N) {
int i = threadIdx.x + blockDim.x * blockIdx.x;
if (i < N / 4) {
float4 va = a[i];
float4 vb = b[i];
c[i] = make_float4(va.x+vb.x, va.y+vb.y, va.z+vb.z, va.w+vb.w);
}
}
// 调用时:线程数/4,因为每线程处理 4 个元素
scalar_add<<<(N+255)/256, 256>>>(a, b, c, N);
vector_add<<<(N/4+255)/256, 256>>>((float4*)a, (float4*)b, (float4*)c, N);

对齐要求

float4 要求数据的起始地址对齐到 16 字节cudaMalloc 分配的内存默认满足此要求,但若使用指针偏移,需确保偏移量是 16 字节的整数倍。

1
2
3
4
5
// 检查对齐
assert((uintptr_t)ptr % 16 == 0);

// 安全的类型转换
float4* ptr4 = reinterpret_cast<float4*>(floatPtr);

常用向量类型

类型 大小 元素数 内存对齐
float2 8B 2 × float 8B
float4 16B 4 × float 16B
int4 16B 4 × int 16B
double2 16B 2 × double 16B
half2 4B 2 × half 4B

Q6:如何优化矩阵乘法(GEMM)?Tiling 策略如何实现?

朴素的矩阵乘法(C = A × B)对全局内存的访问量是 O(N³),但每个元素只参与一次 FMA,算术强度极低,是典型的内存带宽受限问题。Tiling(分块)是解决这一问题的核心方法。

朴素实现的问题

对于 N×N 矩阵,每个输出元素 C[i][j] 需要读取 A 的第 i 行(N 个元素)和 B 的第 j 列(N 个元素)。全局内存总读取量 = 2N³ 个元素,而实际执行 N³ 次 FMA,算术强度 = 0.5 FLOP/Byte——远低于现代 GPU 的算力/带宽比,导致 90% 以上的时间在等待内存。

Tiling 的核心思路

把输出矩阵 C 分成若干 TILE_SIZE × TILE_SIZE 的小块,每个线程块负责计算一个输出 Tile。计算该输出 Tile 所需的 A 和 B 的数据,分批次(沿 K 方向)加载到共享内存,然后在共享内存内完成这一批次的乘法累加,大幅减少全局内存访问次数。

数据复用分析:
输出 Tile 大小为 T×T,每个元素需要 K 次乘加。块内所有线程共享从 A 加载的一列 Tile(T×t 个元素)和从 B 加载的一行 Tile(t×T 个元素)。每个元素被复用 T 次,将全局内存访问量从 O(N³) 降低到 O(N³/T)。

标准 Tiling GEMM 实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
#define TILE 32

__global__ void tiled_gemm(float* A, float* B, float* C, int M, int N, int K) {
__shared__ float As[TILE][TILE]; // A 的 Tile 缓存
__shared__ float Bs[TILE][TILE]; // B 的 Tile 缓存

int row = blockIdx.y * TILE + threadIdx.y;
int col = blockIdx.x * TILE + threadIdx.x;
float acc = 0.0f;

// 沿 K 方向分批次处理
for (int t = 0; t < (K + TILE - 1) / TILE; t++) {

// 协作加载 A 的当前 Tile(合并访问)
As[threadIdx.y][threadIdx.x] =
(row < M && t*TILE + threadIdx.x < K)
? A[row * K + t*TILE + threadIdx.x] : 0.0f;

// 协作加载 B 的当前 Tile(合并访问)
Bs[threadIdx.y][threadIdx.x] =
(t*TILE + threadIdx.y < K && col < N)
? B[(t*TILE + threadIdx.y) * N + col] : 0.0f;

__syncthreads(); // 确保加载完成再计算

// 在共享内存内计算当前 Tile 的贡献(无全局内存访问)
#pragma unroll
for (int k = 0; k < TILE; k++) {
acc += As[threadIdx.y][k] * Bs[k][threadIdx.x];
}

__syncthreads(); // 确保计算完成再加载下一批
}

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

启动配置:

1
2
3
dim3 block(TILE, TILE);                         // 32×32 = 1024 线程/块
dim3 grid((N+TILE-1)/TILE, (M+TILE-1)/TILE);
tiled_gemm<<<grid, block>>>(A, B, C, M, N, K);

进一步优化方向

优化 说明
更大 Tile 增大 TILE(如 64、128),提高数据复用率,但增加共享内存用量
每线程多输出元素 每个线程计算 2×2 或 4×4 的输出 Tile,提高 ILP,减少同步次数
双缓冲 加载下一批 Tile 与计算当前批次重叠(见 Q7)
向量化加载 float4 加载共享内存,减少加载指令数
Tensor Core 用 WMMA API 或 cuBLAS,利用硬件矩阵乘法单元(见 Q9)
Padding 消除 Bank 冲突 __shared__ float As[TILE][TILE+1]

Q7:什么是双缓冲(Double Buffering / Ping-Pong Buffer)技术?

问题根源

标准 Tiling GEMM 的执行流程是:加载 → 同步 → 计算 → 同步 → 加载 → …,加载和计算完全串行。当加载全局内存时,计算单元空闲等待;当计算时,内存系统空闲等待。

双缓冲的思路

准备两份共享内存缓冲区(Buffer 0 和 Buffer 1),交替使用:

  • 计算 Buffer 0 中的数据的同时,异步加载下一批数据到 Buffer 1
  • 计算 Buffer 1 中的数据的同时,异步加载下一批数据到 Buffer 0

这样计算和内存传输完全重叠,消除了等待开销。

1
2
3
4
5
时间轴:
传统: [Load0] [Sync] [Compute0] [Load1] [Sync] [Compute1] ...
双缓冲:[Load0] [Sync] [Compute0+Load1] [Sync] [Compute1+Load2] ...
^^^^^^^^^^^^^^^^^
计算与加载并行

使用 CUDA Pipeline API 实现双缓冲

CC 8.0 引入的 cuda::pipeline + LDGSTS 是实现双缓冲的现代方式,无需占用寄存器作中转:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
#include <cuda/pipeline>

__global__ void double_buffer_gemm(float* A, float* B, float* C,
int M, int N, int K) {
// 两份缓冲
__shared__ float As[2][TILE][TILE];
__shared__ float Bs[2][TILE][TILE];

auto pipeline = cuda::make_pipeline();
int cur = 0; // 当前计算用的 buffer
float acc = 0.0f;

// 预加载第 0 个 Tile
pipeline.producer_acquire();
cuda::memcpy_async(As[cur], &A[...], TILE*TILE*sizeof(float), pipeline);
cuda::memcpy_async(Bs[cur], &B[...], TILE*TILE*sizeof(float), pipeline);
pipeline.producer_commit();

for (int t = 1; t < numTiles; t++) {
int next = 1 - cur;

// 异步加载下一个 Tile 到 next buffer
pipeline.producer_acquire();
cuda::memcpy_async(As[next], &A[...next tile...],
TILE*TILE*sizeof(float), pipeline);
cuda::memcpy_async(Bs[next], &B[...next tile...],
TILE*TILE*sizeof(float), pipeline);
pipeline.producer_commit();

// 等待 cur buffer 加载完成
pipeline.consumer_wait();
__syncthreads();

// 在 cur buffer 上计算(与 next buffer 的加载重叠!)
#pragma unroll
for (int k = 0; k < TILE; k++)
acc += As[cur][threadIdx.y][k] * Bs[cur][k][threadIdx.x];

pipeline.consumer_release();
cur = next; // 切换缓冲
}

// 处理最后一个 Tile
pipeline.consumer_wait();
__syncthreads();
#pragma unroll
for (int k = 0; k < TILE; k++)
acc += As[cur][threadIdx.y][k] * Bs[cur][k][threadIdx.x];
pipeline.consumer_release();

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

双缓冲需要两倍的共享内存,在共享内存紧张时是个约束。可扩展到 N 级缓冲(multi-stage pipeline),进一步隐藏延迟,但需要更多共享内存。


Q8:如何利用 Shared Memory 做 Tile 分块以提升 GEMM 性能?

这是 Q6 的深化,重点讲 Tile 设计的量化分析和关键细节。

量化分析:Tile 大小与带宽节省

设矩阵为 N×N,Tile 大小为 T×T。每个输出 Tile 的计算需要从 A 加载 T×K 元素,从 B 加载 K×T 元素,共 2KT 个全局内存访问,但产生 T² 个输出元素,每个输出元素涉及 K 次 FMA。

  • 朴素实现:每个输出元素加载 2K 个数据 → 全局内存访问量 = 2KN²
  • Tiling(Tile=T):每个 Tile 的 T² 个输出共享 2KT 次加载 → 全局内存访问量 = 2KN²/T

Tile 越大,节省越多。T=32 时节省 32 倍,T=64 时节省 64 倍。

设计 Tile 的关键约束

1
共享内存用量 = 2 × T × T × sizeof(float)  (两个 Tile 的缓存)

以 CC 9.0(最大 228 KB 共享内存/SM)为例:

  • T=32:2×32×32×4 = 8 KB,非常宽松
  • T=64:2×64×64×4 = 32 KB,合理
  • T=128:2×128×128×4 = 128 KB,接近上限,可行但只能单块/SM

Occupancy 约束: 共享内存越多,每个 SM 能并发的线程块越少,需要在 Tile 大小和 Occupancy 之间权衡。

矩阵转置与 Bank 冲突规避

标准 As[TILE][TILE] 的写操作(按列访问)会引发严重 Bank 冲突。加一列 Padding 解决:

1
2
__shared__ float As[TILE][TILE + 1];  // +1 列 padding,错开 Bank 对齐
__shared__ float Bs[TILE][TILE + 1];

这是所有生产级 GEMM 实现的标配。

每线程处理多个输出元素(Register Tiling)

进一步提升:不是每线程只算 1 个输出,而是每线程算 RY × RX 个输出(寄存器 Tile),大幅减少 __syncthreads() 次数,提高计算密度:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
// 每个线程计算 4×4 的输出小块
float acc[4][4] = {0};
// 一次 smem 加载后,在 4×4 的寄存器 Tile 内做乘加
for (int k = 0; k < TILE; k++) {
float a[4], b[4];
#pragma unroll
for (int r = 0; r < 4; r++) a[r] = As[threadIdx.y*4 + r][k];
#pragma unroll
for (int c = 0; c < 4; c++) b[c] = Bs[k][threadIdx.x*4 + c];
#pragma unroll
for (int r = 0; r < 4; r++)
for (int c = 0; c < 4; c++)
acc[r][c] += a[r] * b[c];
}

这正是 CUTLASS 等高性能 GEMM 库的核心优化思路。


Q9:Tensor Core 是什么?如何通过 WMMA API 或 cuBLAS 使用?

Tensor Core 是什么

Tensor Core 是 NVIDIA 从 Volta 架构(CC 7.0)开始引入的专用矩阵乘法硬件单元。它能在一个时钟周期内完成一次小矩阵的乘加操作(MMA,Matrix Multiply-Accumulate),运算形式为:

D=A×B+CD = A \times B + C

对于 FP16 精度,单个 Tensor Core 一个周期可完成 16×16×16 的矩阵乘加(4096 次 FMA 运算)。相比传统 CUDA Core 的标量 FMA,Tensor Core 的矩阵吞吐量提升了约 8~16 倍(取决于架构)。

各架构 Tensor Core 支持的数据类型

CC 架构 支持类型
7.5 Turing FP16, INT8, INT4
8.0 Ampere (A100) FP64, TF32, BF16, FP16, INT8
9.0 Hopper + FP8
10.0 Blackwell + FP6, FP4

方式一:WMMA API(Warp Matrix Multiply-Accumulate)

WMMA 是 CUDA 暴露 Tensor Core 的 C++ API,以 Warp 为单位操作,需要整个 Warp 协作:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#include <mma.h>
using namespace nvcuda;

__global__ void wmma_gemm(half* A, half* B, float* C, int M, int N, int K) {
// 声明 Fragment(矩阵分片,内部布局对程序员不透明)
wmma::fragment<wmma::matrix_a, 16, 16, 16, half, wmma::col_major> a_frag;
wmma::fragment<wmma::matrix_b, 16, 16, 16, half, wmma::row_major> b_frag;
wmma::fragment<wmma::accumulator, 16, 16, 16, float> c_frag;

// 初始化累加器为零
wmma::fill_fragment(c_frag, 0.0f);

// 沿 K 方向循环
for (int k = 0; k < K; k += 16) {
// 从全局/共享内存加载 16×16 的矩阵分片
wmma::load_matrix_sync(a_frag, A + row * K + k, K);
wmma::load_matrix_sync(b_frag, B + k * N + col, N);
// 执行矩阵乘加:c_frag = a_frag * b_frag + c_frag
wmma::mma_sync(c_frag, a_frag, b_frag, c_frag);
}

// 将结果写回
wmma::store_matrix_sync(C + row * N + col, c_frag, N, wmma::mem_row_major);
}

WMMA 支持的矩阵大小(m-n-k):

矩阵 A 类型 矩阵 B 类型 累加器类型 矩阵尺寸
half half float 16×16×16, 32×8×16, 8×32×16
half half half 16×16×16
unsigned char unsigned char int 16×16×16
double double double 8×8×4(CC≥8.0)

使用限制:

  • 必须在条件分支外调用,或确保条件对整个 Warp 完全相同(否则死锁)。
  • Fragment 的内部布局对程序员不透明,不能在不同架构间跨链接传递(会导致未定义行为)。

方式二:cuBLAS(推荐用于生产)

直接使用 WMMA API 正确高效地实现 GEMM 非常复杂,生产代码强烈建议使用 cuBLAS:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#include <cublas_v2.h>

cublasHandle_t handle;
cublasCreate(&handle);

// 设置 Tensor Core 模式
cublasSetMathMode(handle, CUBLAS_TENSOR_OP_MATH);

// FP16 输入,FP32 累加的矩阵乘法
// C = alpha * A * B + beta * C
cublasGemmEx(handle,
CUBLAS_OP_N, CUBLAS_OP_N,
M, N, K,
&alpha,
A, CUDA_R_16F, M, // A: FP16
B, CUDA_R_16F, K, // B: FP16
&beta,
C, CUDA_R_32F, M, // C: FP32(累加器)
CUBLAS_COMPUTE_32F,
CUBLAS_GEMM_DEFAULT_TENSOR_OP);

cublasDestroy(handle);

cuBLAS 会根据矩阵大小和 GPU 型号自动选择最优算法,包括 Tensor Core 的使用策略,是大多数场景下最优的选择。

方式三:CUTLASS

CUTLASS 是 NVIDIA 开源的高性能矩阵运算模板库,提供比 cuBLAS 更细粒度的控制(自定义 Tile 大小、Epilogue、Fused 操作),适合需要定制化 GEMM 的场景。


Q10:什么是 Half-precision(FP16)计算?与 FP32 相比有何优劣?

FP16 的格式

FP16(IEEE 754 半精度浮点数)使用 16 位存储:

字段 位数
符号位 1
指数位 5(偏置值 15)
尾数位 10

数值表示范围约 6×10⁻⁸ 到 6.5×10⁴,精度(十进制有效位)约 3 位。

FP16 vs FP32 对比

维度 FP32 FP16
位宽 32 位 16 位
内存占用 0.5×
内存带宽需求 0.5×
CUDA Core 吞吐量(非 TC)
Tensor Core 吞吐量 8~16×(架构相关)
数值范围 ≈±3.4×10³⁸ ≈±65504
精度 ~7 位十进制 ~3 位十进制

CUDA 中使用 FP16

1
2
3
4
5
6
7
8
9
10
11
12
13
#include <cuda_fp16.h>

__global__ void fp16_kernel(__half* a, __half* b, __half* c, int N) {
int i = threadIdx.x + blockDim.x * blockIdx.x;
// __half 类型的加法
c[i] = __hadd(a[i], b[i]);

// 或使用 __half2 同时处理两个 half(SIMD 风格)
__half2* a2 = (__half2*)a;
__half2* b2 = (__half2*)b;
__half2* c2 = (__half2*)c;
c2[i] = __hadd2(a2[i], b2[i]); // 一条指令处理 2 个 half
}

BF16(Brain Float 16)

BF16 是另一种 16 位格式(CC 8.0 及以上支持),与 FP32 共享相同的指数位宽(8 位),因此数值范围与 FP32 完全相同,只是精度较低(7 位尾数)。深度学习训练中,BF16 比 FP16 更受欢迎,因为它不需要损失缩放(Loss Scaling)来防止溢出。

TF32(TensorFloat-32)

TF32 是 Ampere(CC 8.0)引入的专为 Tensor Core 设计的格式:

  • 指数位:8 位(与 FP32 相同,范围不缩减)
  • 尾数位:10 位(与 FP16 相同,精度低于 FP32)

TF32 在 Tensor Core 中自动用于 FP32 精度的矩阵乘法,无需任何代码改动,只需启用 Tensor Core 模式,是 AI 训练中精度和性能的良好平衡点。

使用建议

  • 推理任务:FP16 或 INT8 是首选,显存减半、速度翻倍,精度损失在可接受范围内(需量化校准)。
  • 训练任务:混合精度训练(AMP)——前向传播用 FP16,梯度累加用 FP32,Tensor Core 自动加速,配合损失缩放(FP16)或不需损失缩放(BF16)。
  • 科学计算:通常需要 FP64(双精度),不可替换。

Q11:什么是 Flash Attention?其核心优化思路是什么?

Flash Attention 不在 CUDA Programming Guide 中直接介绍,但它是基于 CUDA 内存模型的经典算法优化案例,其思路完全建立在本篇已介绍的概念上。

传统 Attention 的问题

标准自注意力(Self-Attention)计算流程:

1
2
3
S = Q × Kᵀ          // [N, N] 注意力得分矩阵
P = softmax(S) // [N, N]
O = P × V // [N, d] 输出

其中 N 是序列长度,d 是头维度。问题在于:N×N 的注意力矩阵(当 N=8192 时,FP16 需要 128 MB)必须写到全局显存,再读回来做 softmax,然后再次写到显存,再读回来做与 V 的乘积。

内存访问量 = O(N²),当 N 变大时,内存访问成为绝对瓶颈,而不是矩阵乘法本身。

Flash Attention 的核心优化:IO-Aware Tiling

Flash Attention(Dao et al., 2022)的关键洞察:注意力计算是内存带宽受限的,N×N 矩阵从不需要完整地存在显存中

核心思路——在 SRAM(共享内存)中分块计算 softmax,从不将完整的 N×N 矩阵写到显存:

  1. 分块(Tiling):将 Q、K、V 分成小块,每次只把一小块加载进共享内存。
  2. Online Softmax:利用数学技巧(log-sum-exp 递推),在遍历 K、V 块的同时,在线计算 softmax,不需要先看完所有 K 才能计算分母。
  3. 融合(Fused Kernel):把"Q×Kᵀ → softmax → × V"的三步操作合并到一个 kernel 中,中间结果全在共享内存,对全局显存的读写量从 O(N²) 降低到 O(N)。

Online Softmax 的数学基础

标准 softmax 计算 P[i] = exp(S[i]) / sum_j(exp(S[j])) 需要先遍历所有 S[j] 求分母。Online softmax 通过维护运行中的最大值 m 和归一化因子 l,允许分块处理:

1
2
3
4
对每个新的块 B:
m_new = max(m_old, max(B))
l_new = exp(m_old - m_new) * l_old + sum(exp(B - m_new))
O_new = (exp(m_old - m_new) * l_old * O_old + exp(B - m_new) * V_B) / l_new

这样每次只需加载一小块 K 和 V,中间结果(O 的累积值)保存在寄存器或共享内存中,最终写一次输出 O 到全局内存。

性能提升

指标 标准 Attention Flash Attention
显存占用 O(N²) O(N)
全局内存读写量 O(N²) O(N)
速度(A100, FP16) 2~4×
最大支持序列长度 ~4K(显存受限) ~100K(显存不再是瓶颈)

Flash Attention 是"IO 感知算法"和"算子融合"思想的完美体现——通过重新设计算法的内存访问模式,而不是更换硬件,获得了 2~4 倍的性能提升。


Q12:如何减少 Kernel Launch 的开销?

Kernel Launch 开销的来源

每次从 CPU 调用 kernel<<<grid, block>>>(args) 都有约 5~20 微秒的 CPU 端开销(驱动程序验证参数、提交命令到 GPU 命令队列、通知 GPU 调度器等)。对于需要执行数千次短小 kernel 的场景(如深度学习推理中的每一层),这些开销会累积为显著的性能损失。

方案一:CUDA Graph(最推荐)

CUDA Graph 将一系列 kernel 启动和内存操作预先"录制"成一个图,然后以极低开销反复执行这个图:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
cudaGraph_t graph;
cudaGraphExec_t instance;
bool graphCreated = false;

for (int step = 0; step < NSTEP; step++) {
if (!graphCreated) {
// 第一次:录制 Graph
cudaStreamBeginCapture(stream, cudaStreamCaptureModeGlobal);
for (int k = 0; k < NKERNEL; k++) {
shortKernel<<<blocks, threads, 0, stream>>>(out, in);
}
cudaStreamEndCapture(stream, &graph);
cudaGraphInstantiate(&instance, graph, NULL, NULL, 0);
graphCreated = true;
}
// 后续每次:以极低开销执行预编译的 Graph
cudaGraphLaunch(instance, stream);
cudaStreamSynchronize(stream);
}

Graph 一旦实例化,每次 cudaGraphLaunch 的 CPU 开销只有 ~2 微秒(相比直接启动多个 kernel 的 n × 10 微秒)。

CUDA Graph 还支持在不重新录制的情况下更新参数(cudaGraphExecUpdate),适合每次迭代只有数值变化(指针不变)的场景。

方案二:算子融合(见 Q13)

将原本 N 个 kernel 的工作合并到 1 个 kernel 中,从根本上减少 kernel 启动次数。

方案三:Persistent Kernel(常驻 Kernel)

让 kernel 永久驻留在 GPU 上,通过原子操作或共享内存的工作队列持续拉取新任务,彻底消除启动延迟:

1
2
3
4
5
6
7
8
9
10
__global__ void persistent_kernel(WorkQueue* queue) {
while (true) {
Task task = queue->dequeue_atomic(); // 原子取任务
if (task.type == TERMINATE) break;
process(task);
}
}
// 只启动一次,持续运行
persistent_kernel<<<numSMs, blockSize>>>(queue);
// CPU 端只需向 queue 添加任务,不需要再次调用 kernel 启动

方案四:grid-stride 循环

已在 Q2 介绍,固定线程块数量,每个线程块处理多批数据,减少总线程块数从而减少调度开销。


Q13:如何做 Kernel Fusion(算子融合)?有何意义?

为什么需要算子融合

深度学习框架(PyTorch、TensorFlow)通常将每个操作(ReLU、LayerNorm、矩阵乘法)实现为独立的 kernel。每个 kernel 都要:

  1. 从全局显存读取输入
  2. 执行计算
  3. 将结果写回全局显存
  4. 下一个 kernel 再读取这些中间结果

对于逐元素操作(activation、normalization、bias add),计算量极小,但每次都要经历完整的"读显存 → 计算 → 写显存"循环,内存带宽成为严重瓶颈,且 kernel 启动开销不可忽视。

融合的收益

将多个顺序执行的 kernel 合并为一个融合 kernel:

  • 消除中间结果的显存读写:中间值直接在寄存器或共享内存中传递
  • 减少 kernel 启动次数:N 次启动 → 1 次启动
  • 提高算术强度:计算和内存访问比例改善

典型例子: GEMM + Bias Add + ReLU

1
2
3
4
5
6
7
8
9
10
11
12
13
14
// 未融合:3 个 kernel,中间结果两次写显存两次读显存
matmul(A, B, C); // 写 C
bias_add(C, bias, C); // 读 C,写 C
relu(C, C); // 读 C,写 C

// 融合:1 个 kernel,中间结果在寄存器中
__global__ void fused_gemm_bias_relu(float* A, float* B, float* bias,
float* C, int M, int N, int K) {
float acc = 0.0f;
// ... GEMM 计算(Tiling)...
acc += bias[col]; // Bias Add(寄存器内完成)
C[row*N + col] = max(0.0f, acc); // ReLU(寄存器内完成)
// 只写一次全局内存
}

Flash Attention 就是终极融合案例

Flash Attention 把"GEMM(Q×Kᵀ)→ Scale → Mask → Softmax → GEMM(×V)"五个操作融合为一个 kernel,中间结果全在 SRAM,全局内存读写从 O(N²) 降到 O(N),这是算子融合的最佳实践。

实现融合的工具

工具 说明
手写融合 kernel 最高性能,最大灵活性,开发成本高
CUTLASS Epilogue 在 GEMM 的写回阶段直接插入 Bias/ReLU/Gelu 等操作
torch.compile / triton 自动/半自动融合,Triton 可用 Python 描述 kernel 由编译器生成高效代码
cuDNN Graph API 声明式地描述算子图,cuDNN 自动选择融合策略
CUDA Graph 不改变单个 kernel 实现,但减少多 kernel 的调度开销(不是真正的融合,但效果类似)

融合的适用场景

场景 是否适合融合
逐元素操作链(activation + normalization) 非常适合,收益大
GEMM + Epilogue 适合,CUTLASS 标准支持
多个独立 GEMM 不适合(每个 GEMM 已是计算密集型),用 Stream 并行更好
内存带宽受限的规约 适合(减少中间写回次数)