一、为什么要做GEMM优化？
在高性能领域，对于矩阵乘（GEMM）的优化是一个非常重要的课题。GEMM可以非常广泛地应用于航空航天、流体力学等科学计算领域，这也是之前HPC的主要应用场景。后来深度学习开展地如火如荼，由于对高算力的需要，也成为HPC的主要应用场景之一。这些年涌现了一系列的深度学习模型。模型里面最耗时的东西，包括卷积、全连接层、attention，都可以转换成GEMM操作。所以说，GEMM优化的重要性，怎么突出都不过分。
本篇文章主要是采纳了cutlass的行文思路，主要介绍GEMM中的数据分块和如何在多级存储进行数据搬运。这也是HPC优化的核心思想，怎么样让数据放在更近的存储上来掩盖计算的延时，从而减少存储墙的影响。文章分为四个方面进行叙述，首先介绍在global memory层面如何进行分块以及数据搬运，随后介绍在shared memory层面如何进行分块以及数据搬运，而后介绍在register层面如何进行分块以及避免bank冲突，最后介绍如何进行预读取prefetch以更好地掩盖访存时延。
二、从global memory到shared memory
假设有矩阵A,B，需要计算矩阵A和B的乘，即矩阵C。A、B、C三个矩阵的维度分别为，m∗k，k∗n，m∗n ，且三个矩阵中的数据都是单精度浮点数。对于C中每一个元素，C[i][j]，可以看作是A的一行和B的一列进行一次归约操作。采用最naive的GEMM算法，在GPU中，一共开启m∗n 个线程，每个线程需要读取矩阵A的一行与矩阵B的一列，而后将计算结果写回至矩阵C中。因而，完成计算一共需要从global memory中进行2mnk 次读操作和m*n次写操作。大量的访存操作使得GEMM效率难以提高，因而考虑global memory中进行分块，并将矩阵块放置到shared memory中。其示意图如下：
对global memory进行分块的GEMM算法示意图见上图右侧。首先将A、B、C三个矩阵划分为多个维度为，bm∗bk，bk∗bn，bm∗bn 的小矩阵块。三个矩阵形成，M∗K，K∗N，M∗N 的小矩阵网格。其中，M=m/bm，N=n/bn，K=k/bk 。随后在GPU中开启M∗N个block，每个block负责C中一个维度为bm∗bn 的小矩阵块的计算。计算中一共有K次迭代，每一次迭代都需要读取A中一个维度为bm∗bk 的小矩阵块和B中一个维度为bk∗bn 的小矩阵块，并将其放置在shared memory中。因而，完成C中所有元素的计算一共需要从global memory中读取M∗N∗K∗（bm∗bk+bk∗bn），即m∗n∗k（1/bn+1/bm） 个单精度浮点数。相比于naive的GEMM算法，访存量减少为原来的1/2∗(1/bm+1/bn) 。通过global memory中分块算法极大地减少了对global memory的访存量。并且，相比于naive算法，对global进行分块可以更充分地利用数据局部性。在naive算法中，每一个线程都需要直接从global memory中取数，其时延非常长，计算性能非常差。而进行分块后，将维度为，bm∗bk，bk∗bn 的小矩阵块先存储到shared memory之中。而后计算单元进行计算时可以直接从shared memory中取数，大大减少了访存所需要的时延。
三、从shared memory到register
随后，我们进一步考虑从shared memory到register的过程。在这里，只分析一个block中的计算。当进行K轮迭代中某一轮迭代时，GPU将维度为bm，bk和bk，bn的小矩阵块存储到shared memory中，而后各个线程将shared memory中的数据存入register中进行计算。
在不对shared memory分块时，一个block中含有bm∗bn个线程，每一个线程负责C中一个元素的计算。则一个block一共需要对shared memory进行2∗bm∗bn∗bk次读操作。而后考虑对shared memory进行分块，对bm∗bn 的小矩阵进行再一次划分，将其划分为多个维度为rm∗rn 的子矩阵。则一个block需要负责X∗Y 个子矩阵。其中，X=bm/rm ，Y=bn/rn。随后，在一个block中开启X∗Y 个线程，每个线程负责一个维度为rm∗rn 的子矩阵的计算。在计算中，一个block一共需要从shared memory读取X∗Y∗(rm+rn)∗bk ，即bm∗bn∗bk∗(1/rm+1/rn) 个单精度浮点数。相比于未分块的算法，对于shared memory中的访存量减少为原来的1/2∗(1/rm+1/rn) 。并且，由于将数据放入register中，可以直接对数据进行运算，减少了从shared memory中取数的时延。
四、register分块
在这里，我们考虑最后一层，即register中的计算，并且只分析一个thread。在完成以上的过程后，对于一个线程而言，它现在拥有：rm 个A矩阵的寄存器值，rn 个B矩阵的寄存器值，以及rm∗rn 个C矩阵的寄存器值。通过这些寄存器的值，需要计算rm∗rn 个数。这需要rm∗rn 条FFMA指令。
这个时候会涉及到寄存器的bank conflict。在NV的GPU中，每个SM不仅会产生shared memroy之间的bank 冲突，也会产生寄存器之间的bank冲突。这一点对于计算密集型的算子十分重要。像shared memory一样，寄存器的Register File也会被分为几个bank，如果一条指令的的源寄存器有2个以上来自同一bank，就会产生冲突。指令会重发射，浪费一个cycle。
我们假设对这个thread来说，rm=4,rm=4 。并且计算C的寄存器以一种非常naive的情况分配，如下图左侧所示。则需要产生16条FFMA指令，列举如下：
1
2
3
FFMA R0, R16, R20, R0
FFMA R1, R16, R21, R1
……
可以从中看出，这会产生大量的register bank冲突，所以需要对参与计算的寄存器重新进行分配和排布,如上图右侧所示。在有些地方，这种方式也可以叫做register分块。
五、数据的prefetch
最后，我们来讲讲如何通过对数据进行prefetch来减少访存的latency。我们再来回顾GEMM的过程，并且仔细地看看这个访存的latency到底是怎么导致的。对于一个block而言，需要计算一个bm∗bn 的矩阵块，这个时候需要进行K次迭代，每次迭代都需要先将来自A和B的两个小块送到shared memory中再进行计算。而从global中访存实际上是非常慢的，所以导致了latency。虽然GPU中可以通过block的切换来掩盖这种latency，但是由于分配的shared memory比较多，活跃的block并不太多，这种延时很难被掩盖。对于一个thread，需要计算一个rm∗rn 的小矩阵，但是必须先将数据从shared memory传到寄存器上，才能开始进行计算。所以导致了每进行一次迭代，计算单元就需要停下来等待，计算单元不能被喂饱。
为此，需要进行数据的Prefetch来尽可能地掩盖这种latency。思想也比较简单，需要多开一个buffer，进行读写分离。示意图如下。当block进行第2轮迭代时，需要对A2和B2进行计算，在计算单元进行计算的同时，我们将A3和B3提前放置到shared memory。而后，在进行第3轮迭代时，就可以直接对shared memory中的A3和B3进行计算，而不需要等待从global memory搬运到shared memory的时间。寄存器上的Prefetch也是同理。
