本文针对NVIDIA T4Turing图灵架构GPU完整拆解CUDA矩阵乘法优化工程。从朴素矩阵乘法存在的显存带宽瓶颈入手逐层介绍Tiling分块算法、Shared Memory共享内存复用、线程级并行、指令级并行ILP等一系列经典优化手段附带编译、性能测试、NCU性能剖面全套实操教程结合Roofline模型评估硬件理论上限与实测性能差距所有代码为工程原生CUDA代码。项目释义矩阵乘法是深度学习、科学计算、仿真模拟最基础核心算子。LLM的Transformer、CNN卷积、数值求解器底层全部依赖大规模矩阵乘。朴素CUDA矩阵乘法直接从全局显存读取矩阵A、B做计算全局显存访问延迟极高带宽成为性能天花板GPU大量计算单元闲置算力利用率很低。本项目基于Turing架构T4显卡一步步优化矩阵乘法Kernel。核心优化手段Tiling分块/切片把超大矩阵切成小块Tile适配GPU片上共享内存Shared Memory把Tile加载到低延迟片上共享内存重复复用数据减少全局显存访问线程级并行利用GPU大量线程每个线程负责输出矩阵C的单个元素指令级并行ILP调整指令排布让GPU多发射单元同时执行多条独立计算指令掩盖访存延迟同时支持和cuBLAS官方库做性能对标使用NCUNVIDIA Compute Profiler采集指标绘制Roofline模型评估算子是否逼近硬件理论极限。核心目标提升GPU算力利用率把矩阵乘法从“显存带宽受限”逐步推向“计算受限”。行业技术知识点Turing T4硬件基础T4是Turing架构SM流式多处理器内置Tensor Core但本项目聚焦基础FP32矩阵乘法优化。全局显存容量大、带宽高但访问延迟高是最大瓶颈Shared Memory共享内存在SM片上带宽远高于全局显存延迟极低但容量很小每个SM的共享内存资源有限WarpGPU调度基本单元32个线程组成1个Warp同一个Warp内线程指令同步执行线程块Block多个Block分配到不同SM执行Block内所有线程共享同一块Shared Memory。Tiling分块原理矩阵乘法CA×BCA\times BCA×B假设A维度M×KB维度K×NC维度M×N。朴素版本计算C[i][j]时每次循环都从全局显存读取A[i][k]、B[k][j]。同一个A/B元素会被多次读取反复访问全局显存。Tiling思路将A、B矩阵切分成尺寸相同的小块Tile。一个Block负责计算输出矩阵C中一块TileBlock先把A、B对应的Tile从全局显存一次性载入Shared Memory所有线程在片上共享内存中循环计算Tile计算完成再加载下一块Tile一次全局显存加载的数据会被多次复用大幅降低全局显存访问次数。经典公式全局显存读取次数减少约1/Tile_Size。Tile尺寸不能无限增大会超出Shared Memory容量限制。线程级并行 TLPGPU天然大规模多线程。每个线程负责输出C矩阵里一个元素。大量线程并发执行当部分线程等待显存数据时硬件调度器切换其他就绪Warp执行用计算掩盖访存延迟。指令级并行 ILP单个线程内编译器/手动调整指令顺序让多条无依赖的算术指令并行发射。普通代码访存指令阻塞后续计算流水线停顿。ILP优化拆分计算预先加载多组数据让GPU在等待一组数据时执行另一组浮点乘法累加指令充分利用SM内多条运算流水线。Roofline 屋顶线模型Roofline模型用来评估GPU算子性能上限两条边界算力上限水平直线GPU理论最大FLOPS带宽上限斜线由显存带宽和算子算术强度决定。算术强度 总浮点运算量 / 总字节访存量算术强度低算子受显存带宽限制落在斜线区域算术强度高算子受GPU计算能力限制水平段优化目标提升算术强度让工作点尽可能靠近Roofline顶点。cuBLASNVIDIA官方高度调优矩阵乘库作为性能基准对标。cuBLAS已经融合Tiling、寄存器分块、TensorCore等全套优化用来衡量自研Kernel距离工业最优实现还有多大差距。整体架构运行流水线生成测试矩阵 A(M×K), B(K×N) ↓ cudaMemcpy H2D ↓ 配置 grid / block 维度 ↓ 执行 CUDA 矩阵乘法内核 ↓ cudaMemcpy D2H ↓ CPU 参考结果对比误差 1e-3 ↓ 输出 GFLOPS / 耗时 / Roofline 数据If you need the complete source code, please add the WeChat number (c17865354792)第一阶段Naive 内核基线最朴素的实现就是把三重循环直接映射到 GPU 线程上每个线程负责输出矩阵 C 的一个元素// mmpy_kernel.cu - naive 版本__global__voidmatmul_naive(constfloat*A,constfloat*B,float*C,intM,intN,intK){introwblockIdx.y*blockDim.ythreadIdx.y;intcolblockIdx.x*blockDim.xthreadIdx.x;if(rowM||colN)return;floatsum0.0f;for(intk0;kK;k){sumA[row*Kk]*B[k*Ncol];}C[row*Ncol]sum;}这个实现有四个致命问题全局内存重复加载每个 C[i][j] 都要从全局内存读取 A 的第 i 行K 个元素和 B 的第 j 列K 个元素。对于 N×N 矩阵A 的每个元素会被读取 N 次B 同理。总全局内存访问量约为2×N³次 float 读取。B 矩阵列访问不合并B[k * N col]中相邻线程的col相差 1但k在内层循环中变化导致 warp 内线程访问的地址跨度为 N无法合并成高效的 128 字节内存事务。零数据复用加载到寄存器的数据只用一次就丢弃没有利用共享内存做中转。无指令级并行单累加器sum编译器难以重排指令掩盖延迟。第二阶段Tiling 共享内存核心思想把矩阵切成 TILE_SIZE×TILE_SIZE 的小块。每个线程块负责计算输出矩阵 C 的一个 TILE_SIZE×TILE_SIZE 子块。在计算前先把 A 和 B 对应的子块从全局内存加载到共享内存然后线程块内所有线程从共享内存读取数据进行计算。量化收益设 TILE_SIZE Tnaive 版本每个输出元素需要 2N 次全局内存读取tiling 版本每次迭代加载 2T² 个元素到共享内存这批数据被 T² 个输出元素复用摊到每个输出只需 2 次全局读。总共 N/T 次迭代每个输出总全局读 ≈2N/T节省因子为T 倍。// mmpy_kernel.cu - tiling 版本#defineTILE_SIZE32__global__voidmatmul_tiled(constfloat*A,constfloat*B,float*C,intM,intN,intK){__shared__floatAs[TILE_SIZE][TILE_SIZE];__shared__floatBs[TILE_SIZE][TILE_SIZE];intbxblockIdx.x,byblockIdx.y;inttxthreadIdx.x,tythreadIdx.y;introwby*TILE_SIZEty;intcolbx*TILE_SIZEtx;floatsum0.0f;// 沿 K 维度分块迭代for(intt0;t(KTILE_SIZE-1)/TILE_SIZE;t){// 协作加载 A 的子块到共享内存intaColt*TILE_SIZEtx;if(rowMaColK)As[ty][tx]A[row*KaCol];elseAs[ty][tx]0.0f;// 协作加载 B 的子块到共享内存intbRowt*TILE_SIZEty;if(bRowKcolN)Bs[ty][tx]B[bRow*Ncol];elseBs[ty][tx]0.0f;// 等待整个 tile 加载完成__syncthreads();// 从共享内存计算部分乘积#pragmaunrollfor(intk0;kTILE_SIZE;k){sumAs[ty][k]*Bs[k][tx];}// 等待所有线程读完再覆盖共享内存__syncthreads();}if(rowMcolN)C[row*Ncol]sum;}关键细节解析1. 两个__syncthreads()缺一不可第一个__syncthreads()保证所有线程把 tile 加载完成后才开始计算第二个保证所有线程读完当前 tile 后才能进入下一轮迭代覆盖共享内存。去掉任何一个都会导致数据竞争。2. 边界填充零矩阵维度不一定是 TILE_SIZE 的整数倍。越界的线程写入 0.0f保证累加结果不受影响。3.#pragma unroll提示编译器展开内层 k 循环减少分支开销给指令调度器更多空间做指令级并行。4. 线程块配置// setGrid.cu dim3 block(TILE_SIZE, TILE_SIZE); // 32x32 1024 线程/块 dim3 grid((N TILE_SIZE - 1) / TILE_SIZE, (M TILE_SIZE - 1) / TILE_SIZE);TILE_SIZE32 时每个块 1024 线程正好是 warp32 线程的整数倍SM 调度效率高。第三阶段消除 Bank Conflict问题本质共享内存被划分为 32 个 bank每个 bank 宽度 4 字节32 位。当 warp 内多个线程同时访问同一 bank 的不同地址时访问会被串行化。在上面的 tiling 实现中Bs[k][tx]的访问模式是warp 内 32 个线程的tx分别为 0~31访问Bs[k][0]~Bs[k][31]这 32 个地址恰好分布在 32 个不同 bank——无冲突很好。但As[ty][k]在内层循环中k变化而ty固定。当不同 warp 的线程访问时如果TILE_SIZE是 32 的倍数As[ty][k]中k相同的线程会命中同一 bank产生冲突。解决方案Padding最简单有效的方法是在共享内存第二维加 1__shared__ float As[TILE_SIZE][TILE_SIZE 1]; __shared__ float Bs[TILE_SIZE][TILE_SIZE 1];这样As[ty][k]的实际地址偏移为ty * (TILE_SIZE1) k打破了 32 的整数倍对齐bank 冲突从 32-way 降为 0。代价多占 TILE_SIZE 个 float对于 TILE_SIZE32 仅多 128 字节几乎可以忽略。第四阶段线程级并行1D Block Tiling / 寄存器分块前面每个线程只算 C 的一个元素。现在让每个线程算多个元素提高寄存器复用率减少共享内存访问次数。#define BM 128 // Block 在 M 维的大小 #define BN 128 // Block 在 N 维的大小 #define BK 8 // 共享内存 K 维分块 #define TM 8 // 每个线程在 M 维算 8 个元素 #define TN 8 // 每个线程在 N 维算 8 个元素 __global__ void matmul_1d_tiled(const float *A, const float *B, float *C, int M, int N, int K) { __shared__ float sA[BM][BK]; __shared__ float sB[BK][BN]; int bx blockIdx.x, by blockIdx.y; int tx threadIdx.x, ty threadIdx.y; int tid tx ty * blockDim.x; // 每个线程负责 TM x TN 个输出元素 float rC[TM][TN] {0.0f}; // 线程在 block 内的行列映射 int threadRow tid / (BN / TN); int threadCol tid % (BN / TN); for (int k0 0; k0 (K BK - 1) / BK; k0) { // 协作加载 A 的 BM x BK 子块 for (int i 0; i BM; i blockDim.x * blockDim.y / BK) { int r i tid; int c tx; // tx 遍历 BK if (by * BM r M k0 * BK c K) sA[r][c] A[(by * BM r) * K k0 * BK c]; else sA[r][c] 0.0f; } // 协作加载 B 的 BK x BN 子块 for (int i 0; i BN; i blockDim.x * blockDim.y / BK) { int r tx; // tx 遍历 BK int c i tid; if (k0 * BK r K bx * BN c N) sB[r][c] B[(k0 * BK r) * N bx * BN c]; else sB[r][c] 0.0f; } __syncthreads(); // 每个线程累加 TM x TN 个元素 for (int k 0; k BK; k) { #pragma unroll for (int i 0; i TM; i) { float aVal sA[threadRow * TM i][k]; #pragma unroll for (int j 0; j TN; j) { rC[i][j] aVal * sB[k][threadCol * TN j]; } } } __syncthreads(); } // 写回全局内存 #pragma unroll for (int i 0; i TM; i) { #pragma unroll for (int j 0; j TN; j) { int r by * BM threadRow * TM i; int c bx * BN threadCol * TN j; if (r M c N) C[r * N c] rC[i][j]; } } }收益每个线程用 64 个寄存器TM×TN64 个 float缓存部分结果共享内存中的数据被复用 TM×TN 次共享内存访问次数降低为原来的 1/64。第五阶段指令级并行ILP与双缓冲指令级并行通过#pragma unroll和多个独立累加器让编译器调度出更多独立指令流掩盖访存延迟float sum0 0.0f, sum1 0.0f, sum2 0.0f, sum3 0.0f; #pragma unroll for (int k 0; k TILE_SIZE; k 4) { sum0 As[ty][k0] * Bs[k0][tx]; sum1 As[ty][k1] * Bs[k1][tx]; sum2 As[ty][k2] * Bs[k2][tx]; sum3 As[ty][k3] * Bs[k3][tx]; } sum sum0 sum1 sum2 sum3;四个累加器之间无依赖GPU 可以并行发射这些 FMAD 指令提高计算单元利用率。双缓冲Double Buffering用两组共享内存缓冲区让加载下一个 tile和计算当前 tile重叠__shared__ float As[2][TILE_SIZE][TILE_SIZE 1]; __shared__ float Bs[2][TILE_SIZE][TILE_SIZE 1]; // 预加载第一个 tile load_tile(A, B, As[0], Bs[0], 0, ...); for (int t 0; t numTiles; t) { int cur t % 2; int nxt (t 1) % 2; // 计算当前 tile compute_tile(As[cur], Bs[cur], sum); // 异步加载下一个 tile简化示意实际需用 cp.async 或拆分同步 if (t 1 numTiles) load_tile(A, B, As[nxt], Bs[nxt], t 1, ...); __syncthreads(); }双缓冲可以隐藏约 40%~50% 的全局内存加载延迟。线程网格配置setGrid.cu// setGrid.cu void setGridDim(int M, int N, int K, dim3 grid, dim3 block) { const int TILE 32; block dim3(TILE, TILE); grid dim3((N TILE - 1) / TILE, (M TILE - 1) / TILE); }配置原则线程数设为 32 的倍数warp 大小常用 128、256、1024每个 SM 驻留的 block 数要足够多保证 occupancyTILE_SIZE 过大会降低 occupancy过小则数据复用不足通常 16~64 之间实验测试矩阵生成与正确性验证genMatrix.cpp utils.cu// genMatrix.cpp - 生成随机测试矩阵voidgenMatrix(float*A,float*B,intM,intN,intK){srand(42);for(inti0;iM*K;i)A[i](float)rand()/RAND_MAX*2.0f-1.0f;for(inti0;iK*N;i)B[i](float)rand()/RAND_MAX*2.0f-1.0f;}// CPU 参考实现voidmatmul_cpu(constfloat*A,constfloat*B,float*C,intM,intN,intK){for(inti0;iM;i)for(intj0;jN;j){floatsum0.0f;for(intk0;kK;k)sumA[i*Kk]*B[k*Nj];C[i*Nj]sum;}}// 误差校验boolverify(constfloat*C_gpu,constfloat*C_ref,intM,intN){floateps1e-3f;for(inti0;iM*N;i){floatdifffabsf(C_gpu[i]-C_ref[i]);floatreffmaxf(fabsf(C_ref[i]),1e-6f);if(diff/refeps){printf(Mismatch at [%d]: gpu%f ref%f\n,i,C_gpu[i],C_ref[i]);returnfalse;}}printf(All elements match (eps%e)\n,eps);returntrue;}环境配置与编译运行硬件与软件要求项目说明GPUNVIDIA T4Turingsm_75或其他 CUDA GPUCUDA Toolkit11.x 及以上编译器nvcc构建工具make性能分析Nsight Compute构建命令# 构建基础版本make-Cbuild_T4# 构建带 cuBLAS 对比的版本make-Cbuild_T4cublastest1# 清理make-Cbuild_T4 cleanMakefile 中关键编译选项NVCC_FLAGS -O3 -archsm_75 --use_fast_math -Xptxas-v # -O3 最高优化等级 # -archsm_75 T4 目标架构 # --use_fast_math 启用快速数学函数__sinf, __expf 等 # -Xptxas-v 输出寄存器/共享内存使用情况运行测试# 运行不同矩阵规模./mmpy-n256./mmpy-n512./mmpy-n1024./mmpy-n2048./mmpy-n4096# 指定 M x K x N./mmpy-m1024-k512-n2048权限修复chmodx tools/run_ncu.shchmodx mmpy查看 GPU 信息nvidia-smi lspci|grep-i--colorvga|3d|2dNsight Compute 分析./tools/run_ncu.sh# 或手动运行ncu--setfull ./mmpy-n1024重点关注指标sm__throughput_average— SM 计算吞吐dram__bytes.sum— 全局内存访问量l1tex__t_sectors_pipe_lsu_mem_global_op_ld.sum— L1 缓存事务smem__transactions— 共享内存事务数smem__bank_conflicts— Bank conflict 次数achieved_occupancy— 实际占用率总结CUDA 矩阵乘法优化是 GPU 异构计算最经典的案例。T4 图灵架构下最大的性能敌人不是计算单元速度而是全局显存延迟与带宽。Tiling 分块 共享内存复用本质就是把数据搬到离计算核心更近的片上存储再搭配线程级并行掩盖访存等待、指令级并行挖掘单线程流水线潜力一步步把矩阵乘算子从带宽瓶颈推向计算瓶颈。Roofline 模型作为评估工具可以量化每一步优化带来的收益避免盲目调参。这套工程适合底层 GPU 算子学习完整覆盖从朴素实现到多层优化、性能剖面、基准对标全链路是理解 GPU 内存层级、并行编程思想的绝佳实践。 本文针对NVIDIA T4Turing图灵架构GPU完整拆解CUDA矩阵乘法优化工程。从朴素矩阵乘法存在的显存带宽瓶颈入手逐层介绍Tiling分块算法、Shared Memory共享内存复用、线程级并行、指令级并行ILP等一系列经典优化手段附带编译、性能测试、NCU性能剖面全套实操教程结合Roofline模型评估硬件理论上限与实测性能差距所有代码为工程原生CUDA代码。项目释义矩阵乘法是深度学习、科学计算、仿真模拟最基础核心算子。LLM的Transformer、CNN卷积、数值求解器底层全部依赖大规模矩阵乘。朴素CUDA矩阵乘法直接从全局显存读取矩阵A、B做计算全局显存访问延迟极高带宽成为性能天花板GPU大量计算单元闲置算力利用率很低。本项目基于Turing架构T4显卡一步步优化矩阵乘法Kernel。核心优化手段Tiling分块/切片把超大矩阵切成小块Tile适配GPU片上共享内存Shared Memory把Tile加载到低延迟片上共享内存重复复用数据减少全局显存访问线程级并行利用GPU大量线程每个线程负责输出矩阵C的单个元素指令级并行ILP调整指令排布让GPU多发射单元同时执行多条独立计算指令掩盖访存延迟同时支持和cuBLAS官方库做性能对标使用NCUNVIDIA Compute Profiler采集指标绘制Roofline模型评估算子是否逼近硬件理论极限。核心目标提升GPU算力利用率把矩阵乘法从“显存带宽受限”逐步推向“计算受限”。行业技术知识点Turing T4硬件基础T4是Turing架构SM流式多处理器内置Tensor Core但本项目聚焦基础FP32矩阵乘法优化。全局显存容量大、带宽高但访问延迟高是最大瓶颈Shared Memory共享内存在SM片上带宽远高于全局显存延迟极低但容量很小每个SM的共享内存资源有限WarpGPU调度基本单元32个线程组成1个Warp同一个Warp内线程指令同步执行线程块Block多个Block分配到不同SM执行Block内所有线程共享同一块Shared Memory。Tiling分块原理矩阵乘法CA×BCA\times BCA×B假设A维度M×KB维度K×NC维度M×N。朴素版本计算C[i][j]时每次循环都从全局显存读取A[i][k]、B[k][j]。同一个A/B元素会被多次读取反复访问全局显存。Tiling思路将A、B矩阵切分成尺寸相同的小块Tile。一个Block负责计算输出矩阵C中一块TileBlock先把A、B对应的Tile从全局显存一次性载入Shared Memory所有线程在片上共享内存中循环计算Tile计算完成再加载下一块Tile一次全局显存加载的数据会被多次复用大幅降低全局显存访问次数。经典公式全局显存读取次数减少约1/Tile_Size。Tile尺寸不能无限增大会超出Shared Memory容量限制。线程级并行 TLPGPU天然大规模多线程。每个线程负责输出C矩阵里一个元素。大量线程并发执行当部分线程等待显存数据时硬件调度器切换其他就绪Warp执行用计算掩盖访存延迟。指令级并行 ILP单个线程内编译器/手动调整指令顺序让多条无依赖的算术指令并行发射。普通代码访存指令阻塞后续计算流水线停顿。ILP优化拆分计算预先加载多组数据让GPU在等待一组数据时执行另一组浮点乘法累加指令充分利用SM内多条运算流水线。Roofline 屋顶线模型Roofline模型用来评估GPU算子性能上限两条边界算力上限水平直线GPU理论最大FLOPS带宽上限斜线由显存带宽和算子算术强度决定。算术强度 总浮点运算量 / 总字节访存量算术强度低算子受显存带宽限制落在斜线区域算术强度高算子受GPU计算能力限制水平段优化目标提升算术强度让工作点尽可能靠近Roofline顶点。cuBLASNVIDIA官方高度调优矩阵乘库作为性能基准对标。cuBLAS已经融合Tiling、寄存器分块、TensorCore等全套优化用来衡量自研Kernel距离工业最优实现还有多大差距。整体架构运行流水线生成测试矩阵 A(M×K), B(K×N) ↓ cudaMemcpy H2D ↓ 配置 grid / block 维度 ↓ 执行 CUDA 矩阵乘法内核 ↓ cudaMemcpy D2H ↓ CPU 参考结果对比误差 1e-3 ↓ 输出 GFLOPS / 耗时 / Roofline 数据If you need the complete source code, please add the WeChat number (c17865354792)第一阶段Naive 内核基线最朴素的实现就是把三重循环直接映射到 GPU 线程上每个线程负责输出矩阵 C 的一个元素// mmpy_kernel.cu - naive 版本__global__voidmatmul_naive(constfloat*A,constfloat*B,float*C,intM,intN,intK){introwblockIdx.y*blockDim.ythreadIdx.y;intcolblockIdx.x*blockDim.xthreadIdx.x;if(rowM||colN)return;floatsum0.0f;for(intk0;kK;k){sumA[row*Kk]*B[k*Ncol];}C[row*Ncol]sum;}这个实现有四个致命问题全局内存重复加载每个 C[i][j] 都要从全局内存读取 A 的第 i 行K 个元素和 B 的第 j 列K 个元素。对于 N×N 矩阵A 的每个元素会被读取 N 次B 同理。总全局内存访问量约为2×N³次 float 读取。B 矩阵列访问不合并B[k * N col]中相邻线程的col相差 1但k在内层循环中变化导致 warp 内线程访问的地址跨度为 N无法合并成高效的 128 字节内存事务。零数据复用加载到寄存器的数据只用一次就丢弃没有利用共享内存做中转。无指令级并行单累加器sum编译器难以重排指令掩盖延迟。第二阶段Tiling 共享内存核心思想把矩阵切成 TILE_SIZE×TILE_SIZE 的小块。每个线程块负责计算输出矩阵 C 的一个 TILE_SIZE×TILE_SIZE 子块。在计算前先把 A 和 B 对应的子块从全局内存加载到共享内存然后线程块内所有线程从共享内存读取数据进行计算。量化收益设 TILE_SIZE Tnaive 版本每个输出元素需要 2N 次全局内存读取tiling 版本每次迭代加载 2T² 个元素到共享内存这批数据被 T² 个输出元素复用摊到每个输出只需 2 次全局读。总共 N/T 次迭代每个输出总全局读 ≈2N/T节省因子为T 倍。// mmpy_kernel.cu - tiling 版本#defineTILE_SIZE32__global__voidmatmul_tiled(constfloat*A,constfloat*B,float*C,intM,intN,intK){__shared__floatAs[TILE_SIZE][TILE_SIZE];__shared__floatBs[TILE_SIZE][TILE_SIZE];intbxblockIdx.x,byblockIdx.y;inttxthreadIdx.x,tythreadIdx.y;introwby*TILE_SIZEty;intcolbx*TILE_SIZEtx;floatsum0.0f;// 沿 K 维度分块迭代for(intt0;t(KTILE_SIZE-1)/TILE_SIZE;t){// 协作加载 A 的子块到共享内存intaColt*TILE_SIZEtx;if(rowMaColK)As[ty][tx]A[row*KaCol];elseAs[ty][tx]0.0f;// 协作加载 B 的子块到共享内存intbRowt*TILE_SIZEty;if(bRowKcolN)Bs[ty][tx]B[bRow*Ncol];elseBs[ty][tx]0.0f;// 等待整个 tile 加载完成__syncthreads();// 从共享内存计算部分乘积#pragmaunrollfor(intk0;kTILE_SIZE;k){sumAs[ty][k]*Bs[k][tx];}// 等待所有线程读完再覆盖共享内存__syncthreads();}if(rowMcolN)C[row*Ncol]sum;}关键细节解析1. 两个__syncthreads()缺一不可第一个__syncthreads()保证所有线程把 tile 加载完成后才开始计算第二个保证所有线程读完当前 tile 后才能进入下一轮迭代覆盖共享内存。去掉任何一个都会导致数据竞争。2. 边界填充零矩阵维度不一定是 TILE_SIZE 的整数倍。越界的线程写入 0.0f保证累加结果不受影响。3.#pragma unroll提示编译器展开内层 k 循环减少分支开销给指令调度器更多空间做指令级并行。4. 线程块配置// setGrid.cu dim3 block(TILE_SIZE, TILE_SIZE); // 32x32 1024 线程/块 dim3 grid((N TILE_SIZE - 1) / TILE_SIZE, (M TILE_SIZE - 1) / TILE_SIZE);TILE_SIZE32 时每个块 1024 线程正好是 warp32 线程的整数倍SM 调度效率高。第三阶段消除 Bank Conflict问题本质共享内存被划分为 32 个 bank每个 bank 宽度 4 字节32 位。当 warp 内多个线程同时访问同一 bank 的不同地址时访问会被串行化。在上面的 tiling 实现中Bs[k][tx]的访问模式是warp 内 32 个线程的tx分别为 0~31访问Bs[k][0]~Bs[k][31]这 32 个地址恰好分布在 32 个不同 bank——无冲突很好。但As[ty][k]在内层循环中k变化而ty固定。当不同 warp 的线程访问时如果TILE_SIZE是 32 的倍数As[ty][k]中k相同的线程会命中同一 bank产生冲突。解决方案Padding最简单有效的方法是在共享内存第二维加 1__shared__ float As[TILE_SIZE][TILE_SIZE 1]; __shared__ float Bs[TILE_SIZE][TILE_SIZE 1];这样As[ty][k]的实际地址偏移为ty * (TILE_SIZE1) k打破了 32 的整数倍对齐bank 冲突从 32-way 降为 0。代价多占 TILE_SIZE 个 float对于 TILE_SIZE32 仅多 128 字节几乎可以忽略。第四阶段线程级并行1D Block Tiling / 寄存器分块前面每个线程只算 C 的一个元素。现在让每个线程算多个元素提高寄存器复用率减少共享内存访问次数。#define BM 128 // Block 在 M 维的大小 #define BN 128 // Block 在 N 维的大小 #define BK 8 // 共享内存 K 维分块 #define TM 8 // 每个线程在 M 维算 8 个元素 #define TN 8 // 每个线程在 N 维算 8 个元素 __global__ void matmul_1d_tiled(const float *A, const float *B, float *C, int M, int N, int K) { __shared__ float sA[BM][BK]; __shared__ float sB[BK][BN]; int bx blockIdx.x, by blockIdx.y; int tx threadIdx.x, ty threadIdx.y; int tid tx ty * blockDim.x; // 每个线程负责 TM x TN 个输出元素 float rC[TM][TN] {0.0f}; // 线程在 block 内的行列映射 int threadRow tid / (BN / TN); int threadCol tid % (BN / TN); for (int k0 0; k0 (K BK - 1) / BK; k0) { // 协作加载 A 的 BM x BK 子块 for (int i 0; i BM; i blockDim.x * blockDim.y / BK) { int r i tid; int c tx; // tx 遍历 BK if (by * BM r M k0 * BK c K) sA[r][c] A[(by * BM r) * K k0 * BK c]; else sA[r][c] 0.0f; } // 协作加载 B 的 BK x BN 子块 for (int i 0; i BN; i blockDim.x * blockDim.y / BK) { int r tx; // tx 遍历 BK int c i tid; if (k0 * BK r K bx * BN c N) sB[r][c] B[(k0 * BK r) * N bx * BN c]; else sB[r][c] 0.0f; } __syncthreads(); // 每个线程累加 TM x TN 个元素 for (int k 0; k BK; k) { #pragma unroll for (int i 0; i TM; i) { float aVal sA[threadRow * TM i][k]; #pragma unroll for (int j 0; j TN; j) { rC[i][j] aVal * sB[k][threadCol * TN j]; } } } __syncthreads(); } // 写回全局内存 #pragma unroll for (int i 0; i TM; i) { #pragma unroll for (int j 0; j TN; j) { int r by * BM threadRow * TM i; int c bx * BN threadCol * TN j; if (r M c N) C[r * N c] rC[i][j]; } } }收益每个线程用 64 个寄存器TM×TN64 个 float缓存部分结果共享内存中的数据被复用 TM×TN 次共享内存访问次数降低为原来的 1/64。第五阶段指令级并行ILP与双缓冲指令级并行通过#pragma unroll和多个独立累加器让编译器调度出更多独立指令流掩盖访存延迟float sum0 0.0f, sum1 0.0f, sum2 0.0f, sum3 0.0f; #pragma unroll for (int k 0; k TILE_SIZE; k 4) { sum0 As[ty][k0] * Bs[k0][tx]; sum1 As[ty][k1] * Bs[k1][tx]; sum2 As[ty][k2] * Bs[k2][tx]; sum3 As[ty][k3] * Bs[k3][tx]; } sum sum0 sum1 sum2 sum3;四个累加器之间无依赖GPU 可以并行发射这些 FMAD 指令提高计算单元利用率。双缓冲Double Buffering用两组共享内存缓冲区让加载下一个 tile和计算当前 tile重叠__shared__ float As[2][TILE_SIZE][TILE_SIZE 1]; __shared__ float Bs[2][TILE_SIZE][TILE_SIZE 1]; // 预加载第一个 tile load_tile(A, B, As[0], Bs[0], 0, ...); for (int t 0; t numTiles; t) { int cur t % 2; int nxt (t 1) % 2; // 计算当前 tile compute_tile(As[cur], Bs[cur], sum); // 异步加载下一个 tile简化示意实际需用 cp.async 或拆分同步 if (t 1 numTiles) load_tile(A, B, As[nxt], Bs[nxt], t 1, ...); __syncthreads(); }双缓冲可以隐藏约 40%~50% 的全局内存加载延迟。线程网格配置setGrid.cu// setGrid.cu void setGridDim(int M, int N, int K, dim3 grid, dim3 block) { const int TILE 32; block dim3(TILE, TILE); grid dim3((N TILE - 1) / TILE, (M TILE - 1) / TILE); }配置原则线程数设为 32 的倍数warp 大小常用 128、256、1024每个 SM 驻留的 block 数要足够多保证 occupancyTILE_SIZE 过大会降低 occupancy过小则数据复用不足通常 16~64 之间实验测试矩阵生成与正确性验证genMatrix.cpp utils.cu// genMatrix.cpp - 生成随机测试矩阵voidgenMatrix(float*A,float*B,intM,intN,intK){srand(42);for(inti0;iM*K;i)A[i](float)rand()/RAND_MAX*2.0f-1.0f;for(inti0;iK*N;i)B[i](float)rand()/RAND_MAX*2.0f-1.0f;}// CPU 参考实现voidmatmul_cpu(constfloat*A,constfloat*B,float*C,intM,intN,intK){for(inti0;iM;i)for(intj0;jN;j){floatsum0.0f;for(intk0;kK;k)sumA[i*Kk]*B[k*Nj];C[i*Nj]sum;}}// 误差校验boolverify(constfloat*C_gpu,constfloat*C_ref,intM,intN){floateps1e-3f;for(inti0;iM*N;i){floatdifffabsf(C_gpu[i]-C_ref[i]);floatreffmaxf(fabsf(C_ref[i]),1e-6f);if(diff/refeps){printf(Mismatch at [%d]: gpu%f ref%f\n,i,C_gpu[i],C_ref[i]);returnfalse;}}printf(All elements match (eps%e)\n,eps);returntrue;}环境配置与编译运行硬件与软件要求项目说明GPUNVIDIA T4Turingsm_75或其他 CUDA GPUCUDA Toolkit11.x 及以上编译器nvcc构建工具make性能分析Nsight Compute构建命令# 构建基础版本make-Cbuild_T4# 构建带 cuBLAS 对比的版本make-Cbuild_T4cublastest1# 清理make-Cbuild_T4 cleanMakefile 中关键编译选项NVCC_FLAGS -O3 -archsm_75 --use_fast_math -Xptxas-v # -O3 最高优化等级 # -archsm_75 T4 目标架构 # --use_fast_math 启用快速数学函数__sinf, __expf 等 # -Xptxas-v 输出寄存器/共享内存使用情况运行测试# 运行不同矩阵规模./mmpy-n256./mmpy-n512./mmpy-n1024./mmpy-n2048./mmpy-n4096# 指定 M x K x N./mmpy-m1024-k512-n2048权限修复chmodx tools/run_ncu.shchmodx mmpy查看 GPU 信息nvidia-smi lspci|grep-i--colorvga|3d|2dNsight Compute 分析./tools/run_ncu.sh# 或手动运行ncu--setfull ./mmpy-n1024重点关注指标sm__throughput_average— SM 计算吞吐dram__bytes.sum— 全局内存访问量l1tex__t_sectors_pipe_lsu_mem_global_op_ld.sum— L1 缓存事务smem__transactions— 共享内存事务数smem__bank_conflicts— Bank conflict 次数achieved_occupancy— 实际占用率总结CUDA 矩阵乘法优化是 GPU 异构计算最经典的案例。T4 图灵架构下最大的性能敌人不是计算单元速度而是全局显存延迟与带宽。Tiling 分块 共享内存复用本质就是把数据搬到离计算核心更近的片上存储再搭配线程级并行掩盖访存等待、指令级并行挖掘单线程流水线潜力一步步把矩阵乘算子从带宽瓶颈推向计算瓶颈。Roofline 模型作为评估工具可以量化每一步优化带来的收益避免盲目调参。这套工程适合底层 GPU 算子学习完整覆盖从朴素实现到多层优化、性能剖面、基准对标全链路是理解 GPU 内存层级、并行编程思想的绝佳实践。