LLM Systems 04: GPU Acceleration

0. Recap

Basic GPU CUDA operations
- memory allocation
- data movement
- creating threads and running on SMs
  - specifying number of threads and number of blocks in a grid
- referring to data in GPU memory within a thread
  - using building index variables to refer to the data

1. Tiling

A100 80G PCIe GPU 规格

FP32	19.5 TFLOPS
Tensor Core Float 32	156 TFLOPS
GPU Memory Bandwidth	1,935 GB/s

每秒 FLOPS 数不仅仅受限于计算单元理论上限，也可能受限于内存访问操作
计算访存比（Compute-to-global-memory-access-ratio）
- 每字节 global memory 访问对应的 FLOPs 数量
对于朴素矩阵乘法实现（每个线程负责输出矩阵一个元素）
- man loop: for (int k = 0; k < N; k++) Pvalue += a[row * N + k] * b[k * N + col];
- 每次循环完成乘加两次浮点运算，访问两个 float 共 8 字节
  - compute-to-global-memory-access: 2 / (2 * 4) = 0.25 FLOP/B

不同存储访问效率

register：约 1 时钟周期
constant cache：5
shared memory/L1 cache：5
L2 cache：200 ~ 300
global memory：500

对于如下 kernel 代码：

1	`C[i] = A[i] + B[i];`

其对应的 GPU 指令：

ld.global.f32 %f1, [%rd1]; // Load A[i]    500 cycle
ld.global.f32 %f2, [%rd2]; // Load B[i]    500 cycle

add.f32 %f3, %f1, %f2; // Perform fp32 addition    1 cycle

st.global.f32 [%rd3], %f3; // Store result

从 global memory 读取数据的时间远长于实际计算的时间

GPU 内存模型

不同 memory type 变量声明：

变量声明 (Variable declaration)	内存 (Memory)	作用域 (Scope)	生命周期 (Lifetime)
`int var;`	Register	Thread	Grid
`int varArr[N];`	Local	Thread	Grid
`__device__ __shared__ int SharedVar;`	Shared	Block	Grid
`__device__ int GlobalVar;`	Global	Grid	Application
`__device__ __constant__ int constVar;`	Constant	Grid	Application

Tile Matrix Multiply kernel

主要步骤：

load the first tile of each input matrix to shared memory
- each thread in the thread block loads one element
- wait for loading using __syncthreads()
each thread in the thread block computes the partial sum from the tiles in shared memory, threads wait for each other to finish
- wait for completing compute using __syncthreads()
load the next tile of each input matrix to shared memory (iterate over the reduction dimension tile by tile)
write the final result to global memory

2. Memory Parallelism

Memory Restriction

GPU	NVIDIA H100	NVIDIA A100
FP32	67 teraFLOPS	19.5 teraFLOPS
Memory	80GB HBM3	80GB HBM2e
Memory Bandwidth	3.35TB/s	2TB/s
SMs	132	104
Shared memory per SM	256K	192K
Registers per SM	64K	64K

shared memory 和 register 资源尤其有限

global memory 的访问要注意做到合并访问（Coalesced）：

合并访问：高效

非合并访问：低效

矩阵转置

naive 实现：

/* macro to index a 1D memory array with 2D indices in row-major order */
/* ld is the leading dimension, i.e. the number of columns in the matrix */

#define INDX( row, col, ld ) ( ( (row) * (ld) ) + (col) )

__global__ void naive_cuda_transpose(int m, float *a, float *c) {
    int myCol = blockDim.x * blockIdx.x + threadIdx.x;
    int myRow = blockDim.y * blockIdx.y + threadIdx.y;
    if( myRow < m && myCol < m ) {
        c[INDX( myCol, myRow, m )] = a[INDX( myRow, myCol, m )];
    }
}

读取 a（合并访问）：
- a[INDX(myRow, myCol, m)] 中，连续的 threadIdx.x 对应连续的 myCol。
- 同一 warp 内的线程访问 a 的同一行连续元素 → 读取是合并的（coalesced）。
写入 c（非合并访问）：
- c[INDX(myCol, myRow, m)] 中，连续的 myCol 对应 c 的连续行。
- 同一 warp 内的线程写入 c 的同一列不同行 → 地址跨度为 m（矩阵宽度）。
- 写入是非合并的（uncoalesced），会导致严重的带宽浪费。
优化方向：
- 使用 shared memory tiling（分块转置）。
- 以合并方式从 a 读取一个 tile 到 shared memory，再以合并方式从 shared memory 写入 c，从而将 global memory 的读写都变成合并访问。

利用 shared memory 优化版本实现：

核心思想：

解决方案：将矩阵分块（tile），每个 block 负责一个 tile。
步骤：
1. 以合并方式从全局内存读取一个 tile 到 shared memory（连续的 threadIdx.x 访问连续列）。
2. 在 shared memory 中完成转置（或任意重排）。shared memory 是片上 SRAM，不需要合并访问，延迟低。
3. 以合并方式将转置后的 tile 从 shared memory 写回全局内存。
效果：全局内存的读和写都变成合并访问，只有 shared memory 内部存在非连续访问，但代价小得多。

Coalesced 矩阵转置代码：

#define INDX( row, col, ld ) ( ( (row) * (ld) ) + (col) )

__global__ void smem_cuda_transpose(int m, float *a, float *c) {
	// shared memory tile
    __shared__ float smemArray[THREADS_PER_BLOCK_X][THREADS_PER_BLOCK_Y];

    const int tileCol = blockDim.x * blockIdx.x;
    const int tileRow = blockDim.y * blockIdx.y;

    // 合并读取：连续 threadIdx.x 访问 a 的连续列
    // 同时完成转置操作：smemArray[threadIdx.x][threadIdx.y]
    smemArray[threadIdx.x][threadIdx.y] = 
        a[INDX(tileRow + threadIdx.y, tileCol + threadIdx.x, m)];

    __syncthreads();

    // 合并写入：连续 threadIdx.x 访问 c 的连续列
    c[INDX(tileCol + threadIdx.y, tileRow + threadIdx.x, m)] = 
        smemArray[threadIdx.y][threadIdx.x];
}

关键点分析：

读取 a：threadIdx.x 连续 → tileCol + threadIdx.x 连续 → 同一行连续列 → 合并。
写入 c：threadIdx.x 连续 → tileRow + threadIdx.x 连续 → 同一行连续列 → 合并。
转置在 shared memory 内完成：
- 写入 shared memory：smemArray[threadIdx.x][threadIdx.y]
- 读取 shared memory：smemArray[threadIdx.y][threadIdx.x]
- 两者互为转置，且 shared memory 无需合并，非连续访问代价低。
为什么 smemArray 声明为 [X][Y] 而不是 [Y][X]？
- 这个 shared memory 数组保存的是 转置后的 tile。
- 写入时：smem[tx][ty] = A[ty][tx]，即 A 的行变成 smem 的列。
- 读取时：smem[ty][tx] = C[ty][tx]，再写回 C 的同一行。
- 这种布局使得：
  - 从 A 读：tx 连续 → A 同一行连续列 → 合并。
  - 写回 C：tx 连续 → C 同一行连续列 → 合并。
- 转置存储让全局内存两端的访问方向都被“转正”，只有 shared memory 内部发生非连续访问，代价小得多。
- 若按常规 [Y][X] 声明，读取 smem 时 tx 会映射到行，导致写回 C 时地址跨步，丧失合并。
__syncthreads()：确保 tile 完全加载到 shared memory 后，再进行写回。

以一个 tile 为例，第一阶段从 A 读到 shared memory：

A tile (global memory, row-major)              Shared memory: smem[tx][ty]
         tx=0  tx=1  tx=2  tx=3                ty=0   ty=1   ty=2   ty=3
ty=0     A00   A01   A02   A03          tx=0   A00    A10    A20    A30
ty=1     A10   A11   A12   A13          tx=1   A01    A11    A21    A31
ty=2     A20   A21   A22   A23          tx=2   A02    A12    A22    A32
ty=3     A30   A31   A32   A33          tx=3   A03    A13    A23    A33

读取: 固定 ty，tx 连续 → A 的同一行连续列 → 合并访问
写入 shared: smem[tx][ty] = A[ty][tx] → 转置存储

第二阶段从 shared memory 写到 C：

Shared memory: smem[tx][ty]                    C tile (global memory, row-major)
         ty=0   ty=1   ty=2   ty=3              tx=0   tx=1   tx=2   tx=3
tx=0     A00    A10    A20    A30        ty=0   A00    A10    A20    A30
tx=1     A01    A11    A21    A31        ty=1   A01    A11    A21    A31
tx=2     A02    A12    A22    A32        ty=2   A02    A12    A22    A32
tx=3     A03    A13    A23    A33        ty=3   A03    A13    A23    A33

读取: 固定 ty，tx 连续 → smem[ty][tx] 连续 → 高效
写入: C[tileCol+ty][tileRow+tx] → 固定 ty，tx 连续 → C 同一行连续列 → 合并

从线程视角看不同阶段线程如何 load/store/与 shared memory 交互：

以 4×4 tile、块 (0,0) 为例。每个线程 (tx, ty) 只做四件事：

线程 (tx, ty)	从 A 读取	写入 shared memory	从 shared memory 读取	写入 C
(0,0)	A[0][0]	smem[0][0] = A[0][0]	smem[0][0] = A[0][0]	C[0][0] = A[0][0]
(1,0)	A[0][1]	smem[1][0] = A[0][1]	smem[0][1] = A[1][0]	C[0][1] = A[1][0]
(2,0)	A[0][2]	smem[2][0] = A[0][2]	smem[0][2] = A[2][0]	C[0][2] = A[2][0]
(3,0)	A[0][3]	smem[3][0] = A[0][3]	smem[0][3] = A[3][0]	C[0][3] = A[3][0]
(0,1)	A[1][0]	smem[0][1] = A[1][0]	smem[1][0] = A[0][1]	C[1][0] = A[0][1]
(1,1)	A[1][1]	smem[1][1] = A[1][1]	smem[1][1] = A[1][1]	C[1][1] = A[1][1]

只看前两行，就能发现规律：

读取 A：线程 (tx, ty) 读 A[ty][tx]。
- 同一行 ty 的线程（tx=0,1,2,3）读 A 的同一行连续列 → 合并。
写入 smem：线程 (tx, ty) 写 smem[tx][ty]。
- 这是转置存储：A 的行变成 smem 的列。
读取 smem：线程 (tx, ty) 读 smem[ty][tx]。
- 对于 ty 相同的一组线程，协作读取 shared memory 的一行元素，对应 A 的一列的 tile。
写入 C：线程 (tx, ty) 写 C[ty][tx]。
- ty 相同的同一行线程写 C 的同一行连续列 → 合并。
每个线程只做四件事：
1. 读 A：线程 (tx, ty) 读 A[ty][tx] → tx 连续 → 合并。
2. 写 smem：线程 (tx, ty) 写 smem[tx][ty] → 转置存储（A 的行变成 smem 的列）。
3. 读 smem：线程 (tx, ty) 读 smem[ty][tx] → 等价于读 A[tx][ty]，正好是写入 C 所需的值。
4. 写 C：线程 (tx, ty) 写 C[ty][tx] → tx 连续 → 合并。

Bank Conflict

shared memory 有 32 个 banks，每个 bank 宽 4 字节。
一个 warp 有 32 个线程。如果这 32 个线程访问的地址落在同一个 bank 的不同位置，硬件无法并行处理，只能串行访问，导致性能下降。这就是 bank conflict。
如果 32 个线程访问的地址恰好落在32 个不同的 bank，则可以完全并行，无冲突。
冲突的访问模式：smemArray[threadIdx.x][threadIdx.y]
- 地址 = tx * 32 + ty → bank = ty。
- 同一 warp 内 ty 固定，所有线程 bank 相同 → 32-way conflict，慢。
无冲突的访问模式：smemArray[threadIdx.y][threadIdx.x]
- 地址 = ty * 32 + tx → bank = tx。
- 同一 warp 内 tx 从 0 到 31 变化 → 命中 32 个不同 bank → 无冲突，快。
矩阵转置中的应用：
- 转置时 shared memory 的读写会自然产生 smem[tx][ty] 和 smem[ty][tx] 两种模式。
- 一种会产生 bank conflict，另一种不会。
- 常见优化：将 smemArray[TILE][TILE] 声明为 smemArray[TILE][TILE + 1]（padding），打破 bank 对齐，消除冲突。

消除 bank conflict 矩阵转置 kernel 代码：

#define INDX( row, col, ld ) ( ( (row) * (ld) ) + (col) )

__global__ void smem_cuda_transpose(int m, float *a, float *c) {
	// 关键：padding 1 个元素
    __shared__ float smemArray[THREADS_PER_BLOCK_X][THREADS_PER_BLOCK_Y+1];

    const int tileCol = blockDim.x * blockIdx.x;
    const int tileRow = blockDim.y * blockIdx.y;

    smemArray[threadIdx.x][threadIdx.y] = a[INDX( tileRow + threadIdx.y, tileCol + threadIdx.x, m )];
    __syncthreads();

    c[INDX( tileCol + threadIdx.y, tileRow + threadIdx.x, m )] = smemArray[threadIdx.y][threadIdx.x];
    return;
}

padding 后线程对 shared memory 访问模式示意图：

数组同一列由于 padding 的存在，不同行散落在不同 bank 上，消除了 shmem[tx][ty] 模式访问的 bank conflict。

3. Sparse Matrix Multiplication

CSR (Compressed Sparse Row) 格式

背景：稀疏矩阵中绝大多数元素为 0，直接使用稠密矩阵存储和计算会浪费大量内存和带宽。
CSR 使用三个一维数组（设矩阵有 m 行，非零元素个数为 nnz）：
- data：长度为 nnz，按行主序存储所有非零元素的值。
- col_index：长度为 nnz，存储每个非零元素对应的列索引。
- row_ptr：长度为 m + 1，存储每行第一个非零元素在 data / col_index 中的起始位置。row_ptr[i] 到 row_ptr[i+1]-1 即为第 i 行的非零元素。
示例（4×4 矩阵）：
- row_ptr = [0, 3, 4, 6, 7]
- col_index = [0, 2, 3, 1, 2, 3, 3]
- data = [a, b, c, d, e, f, g]
- 第 0 行的非零元素在 data[0..2]，对应列索引 [0, 2, 3]。

Sparse Matrix-Vector Multiplication 算法（串行视角）

计算 \(Y = A \times X\)。
对每一行 row：
- 获取该行的非零元素范围 [row_start, row_end)。
- 遍历该范围内的每个非零元素 ele，计算 dot += x[col_index[ele]] * data[ele]。
- 将结果写入 y[row]。

for(int row = 0; row < n; row++) {
	float dot = 0;
	int row_start = row_ptr[row];
	int row_end = row_ptr[row + 1];
	for(int el = row_start; el < row_end; el++)
	{
		dot += x[col_index[el]] * data[el];
	}
	y[row] += dot;
}

CUDA 实现（SpMVCSRKernel）

__global__ void SpMVCSRKernel(float *data, int *col_index, int *row_ptr, 
                              float *x, float *y, int num_rows) {
    int row = blockIdx.x * blockDim.x + threadIdx.x;
    if(row < num_rows) {
        float dot = 0;
        int row_start = row_ptr[row];
        int row_end = row_ptr[row + 1];
        for(int ele = row_start; ele < row_end; ele++) {
            dot += x[col_index[ele]] * data[ele];
        }
        y[row] += dot;
    }
}

线程映射：每个线程负责一行 row 的计算。
边界检查：if(row < num_rows) 确保不越界。

性能分析

优点：
- 内存占用与计算量正比于非零元素数 nnz，而非 O(m n)。
- 避免了零元素的无效计算。
缺点：
- 非合并访问：x[col_index[ele]] 的访问完全依赖 col_index，列索引不连续，导致对 x 的读取是随机的、非合并的。
- 负载不均衡：不同行的非零元素数量差异可能很大，有的线程计算量远大于其他线程，导致 warp 内线程负载不均衡，影响吞吐。
- 行间并行度受限：如果矩阵行数较少，难以将 GPU 的 SM 喂满。
优化方向：
- 调整线程映射（如一个 warp 处理一行，利用合并访问）。
- 使用 ELL / JDS 等针对 GPU 优化的稀疏格式

4. cuBLAS

概述

cuBLAS：CUDA Basic Linear Algebra Subroutine library（CUDA 基础线性代数子程序库）。
一个轻量级库，专用于 **GEMM 操作，同时也提供基础的向量和矩阵运算。

生命周期管理

必须在使用前调用：

1	`cublasStatus_t cublasCreate(cublasHandle_t *handle)`

必须在结束后调用：

1	`cublasStatus_t cublasDestroy(cublasHandle_t handle)`

所有 cuBLAS API 调用都需传入 cublasHandle_t handle。

核心 API

1. 向量点积（Vector Dot Product）

cublasStatus_t cublasSdot(
    cublasHandle_t handle, int n,
    const float *x, int incx,
    const float *y, int incy,
    float *result)

计算 \(x \cdot y\)。
incx / incy：向量元素的步长（通常为 1）。

2. 矩阵-向量乘法（Matrix-Vector Product） 公式： \[ y = \alpha A \cdot x + \beta y \]

cublasStatus_t cublasSgemv(
    cublasHandle_t handle, cublasOperation_t trans,
    int m, int n,
    const float *alpha,
    const float *A, int lda,
    const float *x, int incx,
    const float *beta,
    float *y, int incy)

trans：是否转置 A（CUBLAS_OP_N 或 CUBLAS_OP_T）。
m / n：矩阵 A 的维度。
lda：leading dimension of A（A 的列数）。

3. 矩阵-矩阵乘法（Matrix-Matrix Multiplication） 公式： \[ C = \alpha A \cdot B + \beta C \]

cublasStatus_t cublasSgemm(
    cublasHandle_t handle, cublasOperation_t transa,
    cublasOperation_t transb,
    int m, int n, int k, const float *alpha,
    const float *A, int lda,
    const float *B, int ldb,
    const float *beta,
    float *C, int ldc)

transa / transb：是否转置 A / B。
m / n / k：矩阵维度。
lda / ldb / ldc：各矩阵的 leading dimension。

关键注意点

列主序（Column-Major）：cuBLAS 默认使用列主序存储。如果数据是行主序（C/C++ 默认），通常需要利用转置参数或手动调整索引，这是常见的坑。
句柄（handle）：所有 cuBLAS 函数都需要传入 handle，用于管理上下文和资源。
S 前缀：cublasSdot、cublasSgemv、cublasSgemm 中的 S 表示单精度浮点数（float）。对于双精度用 D（如 cublasDgemm），半精度用 H（如 cublasHgemm）。
性能：cuBLAS 内部高度优化，通常使用 Tensor Core、Tiling、多级流水线等高级技术，性能远超手写 kernel。

cuBLASLt

定位：cuBLAS 的灵活扩展，专为 DL/ML 中的 GEMM 设计。
与 cuBLAS 的区别：
- 描述符 + 启发式 API，而非固定参数。
- 支持任意数据布局、FP8/FP4 等低精度格式。
- 支持 epilogue 融合（bias + 激活）。
核心 API：
- cublasLtCreate / cublasLtDestroy
- cublasLtMatrixLayoutCreate
- cublasLtMatmulDescCreate
- cublasLtMatmulAlgoGetHeuristic（查询最优算法）
- cublasLtMatmul（执行）
调优要点：
- 启发式结果建议缓存重用。
- workspace 大小可调（可能影响性能）。
- 不同问题尺寸的最优算法可能不同，建议 autotune。

5. Summary

Tiling for efficient matrix computation
Coalesced memory access
Sparse matrix representation and multiplication
cuBLAS
- readily available vector, matrix-vector, matrix-matrix operations

ML System

#CMU 11868 #LLM System

LLM Systems 04: GPU Acceleration

https://arcsin2.cloud/posts/2026/09/1906032670/

作者

arcsin2

发布于

2026年9月19日

许可协议

PMPP Chapter 01: Introduction 上一篇

LLM Systems 03: GPU Programming-2 下一篇