gemm优化

最终测试结果：固定条件时比cuBlas小优

测试尺寸： 2048 * 2048 * 2048

测试环境 cpu: i7-13700K

gpu: RTX4070ts

cuda： 13.1

vs: Microsoft Visual Studio Community 2022

最终结果

固定频率2100MHz时我们的实现比cuBlas更好(nvidia-smi -lgc 2100,2100)：

MyKernel GPUexecution time:6783.493164ms，TFLOPS:25.326

cuBLAS GPU execution time:6831.067383 ms，TFLOPS:25.150

不固定频率(通过nvidia-smi dmon 观察频率会掉到2100MHz)：

MyKernel GPU execution time: 6571.359375 ms， TFLOPS: 26.144

cuBLAS GPU execution time: 5911.912109 ms， TFLOPS: 29.060

性能分析

理论计算：计算强度和机器平衡点

AI = 内存量 / FMA 计算量比如我的kernel中为64x64的共享内存时，从内存读到共享内存的数据量为： As + Bs = 2 × 64² × 4B = 32 KB（每线程 8 条 float4 = 128 B）计算量：64³ = 262,144 次 FMA = 524,288 FLOP（每线程 4×4×64 = 1024 FMA）

那么 AI = 524288 FLOP ÷ 32768 B = 16 FLOP/Byte

而我的硬件为4070ts，机器平衡点（峰值算力 ÷ 显存带宽） = 44.10 TFLOPS ÷ 672.3 GB/s = 65.4

只有 AI > 平衡点时计算才是瓶颈，所以理论上我这里瓶颈应该是带宽

实际测试：ncu

ncu -k myGemm –launch-skip 8 –launch-count 1 –set full .\gemm.exe

或只用关键字： ncu -k myGemm –metrics sm__throughput.avg.pct_of_peak_sustained_elapsed,gpu__compute_memory_throughput.avg.pct_of_peak_sustained_elapsed,smsp__warp_issue_stalled_long_scoreboard_per_warp_active.pct,smsp__warp_issue_stalled_barrier_per_warp_active.pct,dram__bytes.sum .\gemm.exe

得到结果

  myGemm(float *, float *, float *, float, float, int, int, int) (320, 320, 1)x(16, 16, 1), Context 1, Stream 7, Device 0, CC 8.9
    Section: GPU Speed Of Light Throughput
    ----------------------- ----------- ----------------
    Metric Name             Metric Unit     Metric Value
    ----------------------- ----------- ----------------
    DRAM Frequency                  Ghz            10.19
    SM Frequency                    Ghz             2.17
    Elapsed Cycles                cycle    2,271,298,602
    Memory Throughput                 %            81.41
    DRAM Throughput                   %            81.41
    Duration                          s             1.05
    L1/TEX Cache Throughput           %            79.48
    L2 Cache Throughput               %            35.22
    SM Active Cycles              cycle 2,268,501,054.44
    Compute (SM) Throughput           %            52.35
    ----------------------- ----------- ----------------

可以看到确实加载内存的利用率Memory Throughput很大，计算利用率Compute (SM) Throughput不够就是说主要还是卡内存，需要优化内存访问，符合理论分析。

优化内容

分片缓存1：显存到共享内存

最简单的乘法实现无疑是每个线程计算一个结果数组的元素。但是这样的话采取上面的理论分析，其实可以看到，单次计算的计算量无法得到足够的内存带宽。我们需要让这个比例更协调。

最重要的思想就是分片。对于矩阵A:m x k 和 B：k x n来说，

A的每一行都会用来和B的所有n列都计算一次乘积。如果每次都重新读取A的那一行，无疑会导致n-1次多余的读取，浪费大量内存带宽

因此我们可以通过缓存这个A行供给B的列来减少这种情况. 同理也要缓存B的列来供给A的行。而显然我们无法把完整的矩阵存入高速存储空间。因此只能通过切块方法达到优化效果。

我们的数据最开始放在显存，因此首先往上一级缓存就是到共享内存（不讨论L1等无法控制的缓存）。

我们设定BM, BN, BK对应缓存矩阵A，B的大小为 BM x BK, BK x BN。一般设kernel为比较大的16 x 16. 为了方便表示，最初让取数据操作每个线程取1个数据存入，因此可以让BM=BN=BK=16. 每次取完，然后每个线程使用缓存矩阵计算临时结果。滚动k对应维度累计求和即可。

此时即可达到cuBlas的60%效果了

分片缓存2：共享内存到寄存器

此时继续分析仍然是内存受限，而共享内存之上还有一级寄存器可用。因此我们可以继续做缓存，用寄存器缓存计算结果C. 让每个线程的计算可以计算对应C的一个TM x TN的矩阵块。而根据上述的理论计算时，我的BM=BN=BK=64, 而TM = TN = 4。理论计算显示需要让计算乘以4倍才比较均衡，因此，我们需要调整BM=BN=128，而此时共享内存就不足了，因此必须让BK=8， TM=TN=8 从而减轻负载，通过多次迭代移动来遍历K来达到类似效果。

向量化内存读写

此时会有一个优化点是，我们按float操作的话，会形成多次for带来的多次指令。我们可以通过float4变量来一次性实现类似效果。而这样的话，因为每轮K迭代都是BM * BK = 1024 = 4 * 16 * 16 = 4 x block线程数。因此刚好每个线程取一个float4即可。非常完美。

同时，从共享内存读取到寄存器也可以类似操作，减少读共享内存次数75%.

此番修改后即可得到80%的cuBlas性能。

杂项：预取缓存1， Swizzle 增强L2缓存命中

此时读显存、计算两个步骤为完全串行，可以优化成两个步骤，读下一批次的float4，计算当前批次。从而形成流水线。

此时继续观察每个block（也即sm分片对应的计算结果区域）的计算结果区域是按顺序依次排列的，也即从0-x按列优先，然后逐行遍历。这样有个问题是B数组的读取每次都是从DRAM读取，无法缓存。通过重新映射方法来再在显存和共享内存之间做一层缓存（L2）。 block顺序从依次，改成按group分组，满group则先跳到下一“行”，而不是继续这一列。比如按照BM, BN, BK分块后，原本结果数组是每行分别有16个块的。那么我们按group = 8设置，第9个块不再处理第一行的第9个块，而是处理第2行的第一个块，此时计算这个块就会用到处理第一行第一列block时的B的列数据（在L2）。从而减少DRAM读次数，达到加速效果。

此时可以得到 90% 的cuBlas性能。

关于GPU频率影响

此时到了最奇妙的地方，因为我们无论如何都差10%无法达到cuBlas的效果，而通过nvidia-smi dmon观察频率可以发现，我们的计算过程中降频非常厉害，而cublas稍微少一些，这是我最不能理解的地方。如果我们锁定频率，此时，就会惊奇的发现，我们的性能已经超过cuBlas了~。