Anatomy of High-Performance Many-Threaded Matrix Multiplication 笔记_anatomy of high-performance matrix multiplication

文章主要是对GEMM中几种潜在的多线程实现方法进行性能分析。看本篇文章需要先看前置论文：
Anatomy of High-Performance Matrix Multiplication
Anatomy of High-Performance Many-Threaded Matrix Multiplication

GEMM分块执行流程（单线程）

上图是单线程下GEMM中各分块所在内存层次的位置。它对应的是列主序下最优的情况。行主序则是将A的小分块放在L1，B的分块放在L2。

具体实现的伪代码如下图所示：

最外层的两层循环是用于选择A的一大列和B的一大行，也即jc和pc。这两个大循环内部的一次循环可以计算C的“一层皮”。ic则是在A的一大列和B的一大行确定下来之后，在A的一大列中选择一个分块，分块的高度为mc。

jr是在A中的分块确定下来之后，在B的一大行中选择一个小分块（上图红色）。

ir是在B中的分块（上图红色）确定下来之后，选择A分块中的小分块（上图蓝色），并在内部做register kernel。

单线程实现中潜在的并行机会

1. register kernel内的并行

就是在计算C中mr x nr分块时并行，但是作者不建议这么细粒度的并行，原因有：

1. 不同的线程需要对计算结果进行累加，才能得到一个C中的mr x nr分块。线程的reduce操作开销太大。
2. 对于更新一个C中mr x nr分块，每个线程执行的计算量太少了，无法掩盖将更新结果的开销（将结果写回内存中C存放的位置）。

所以应该将register kernel作为并行的最小计算单位。

2. ir循环上的并行

如上图所示。这里为了方便起见，将细长条的分块称为sliver。在这个循环层次中，实际上是对A的一个分块block和B的一个sliver进行计算（也就是遍历A分块block中的sliver，分别与B的sliver进行计算）。因此每个线程负责对一个A的sliver（蓝色）和B的sliver（红色）进行计算，那么这些线程共享B的sliver（根据上面的内存层次图，B的sliver是放在L1 cache中的）。这其实就是对多个register kernel进行并行计算。

在这个循环层次上，这个循环的迭代次数为ceil(mc / mr)（ceil为向上取整）。因为在这个循环层次上主要就是执行A的block和B的sliver的乘法运算，而主要的访存开销来自于将所有线程共享的B的sliver从L3 cache读入到L1 cache中，所以为了掩盖这个访存开销，这个循环层次上的计算时间必须足够长。换句话说，就是ceil(mc / mr)必须足够大，在这个循环层次上的并行才有收益。

3. jr循环上的并行

如上图所示。这个循环层次上的并行，比上一个层次上的并行又提升了一个维度。这次是在循环选择B中的sliver时进行并行。（也就是遍历B的sliver，分别与A分块block进行计算）。每个线程负责对一个A的block和B的sliver（红色）进行计算，那么这些线程共享A的block。

在这个循环层次上，这个循环的迭代次数为ceil(nc / nr)（ceil为向上取整）。在这个循环层次上主要的访存开销为对所有线程共享的A的block进行packing，并且在packing的过程中将block从memory加载到L2 cache中。所以同理，为了掩盖这个访存开销，这个循环层次上的计算时间必须足够长。换句话说，就是ceil(nc / nr)必须足够大，在这个循环层次上的并行才有收益。

还有一个问题，在某些CPU上，L2 cache是每个核心独有的（例如现在正在测试的鲲鹏920），那么如果在实现过程中对packing过程进行并行，那么不同核心的L2 cache中会保存着A block的一部分（因为每个线程在packing的过程中，将它负责packing的那部分放入到L2 cache中）。所以在进行计算的过程中，缓存一致性协议cache coherency protocol的实现方式会影响到不同L2 cache间传送数据的效率，从而影响计算的效率。如果不同L2 cache间传送数据的开销较小，register kernel（micro-kernel）计算量能够掩盖，那么L2 cache共享不会影响到计算的效率。

4. ic循环上的并行

如上图所示。这个循环层次上的并行，比上一个层次上的并行又提升了一个维度，这是在循环选择A中的block时进行并行。每个线程拥有一个独占的A的block，它存在于各自的L2 cache中，以及一个共享的B的panel（也就是B的一行），它存在于共享的L3 cache或者memory中。当 m < mc（mc是每个block的高度）* thread_num 时，并行的效果会收到影响，因为各线程间负载不均衡。所以这个并行方法需要m足够大。

因为本身L3 cache是多个核心共享的，所以不用担心各线程对B的panel进行重排（加载到L3 cache）时所采用的缓存一致性协议的开销。

5. pc循环上的并行

如上图所示。这是在计算C的“每一层皮”上进行并行。每个线程选取A的一列和B的一行，计算得到C的“一层皮”。多个线程的计算结果进行叠加，就是最终的计算结果。因为每一个线程都需要对C进行更新，这里会涉及到线程间的竞争（race），因此可能需要使用操作系统中的锁机制（lock）来保证结果的正确，或者是每个线程暂存自己的计算结果，最后再进行规约（reduce）。因为涉及到线程间的竞争，所以这个并行方法的性能较差，除了在一些极端情况下：

• C本身的size太小，前面的几种并行方式不能充分并行（线程创建、启动和回收的开销大于计算时间）；
• 对结果进行reduce的开销明显小于计算开销。

6. jc循环上的并行

如上图所示。每一个线程对整一个A和B一行上的分块进行计算。线程间共享A，A放在memory中；单个线程分别占有B一行上的一个分块，这个分块放在L3 cache中。这种并行方式比较特殊，适合multi-socket的架构，它一般支持non-uniform memory access (NUMA)。每一个NUMA node都有自己独立的L3 cache，刚好用于存放线程自己独立的B一行上的一个分块。

案例分析——Intel Xeon Phi

1. 硬件架构细节

这一块不多说，可以直接看原文。总之大致的意思是有60个核心，每个核心支持4个硬件线程。核心内部共享L1 cache。每个核心有两条流水线，其中只有一条流水线用于执行SIMD浮点运算操作，而每个线程每隔一个时钟周期只能向每个管道发出一条指令。所以每个核心至少拥有2个硬件线程，才能使SIMD浮点运算单元流水线满流。论文中直接在每个核心中使用了4个硬件线程，总共有240个线程。

2. BLIS的参数选择

Intel Xeon Phi一个核心的四个线程共享L1 cache，因此如果需要在m和n维度同时并行时，这意味着在L1 cache中至少拥有A的2个sliver和B的2个sliver。而为了将这么多的数据存放到L1 cache中，那只能减少kc的值。而如果减少了kc的值，那么register block的计算量又不足以掩盖更新C中mr x nr 小分块的访存开销。因此论文中提出了一个方法，在L2 cache中固定一个区域作为虚拟L1 cache。当然访存的开销还是相同的。

register block的大小为30 x 8，也就是mr为30，nr为8。mc为120，kc为240。因为没有L3 cache，B的panel只能放在memory中，所以nc只受到了memory的限制。nc直接选择为14400，它是下面所有测试矩阵中n的最大值。

3. 并行的循环选择

• ir循环（不选择）：mc为120，mr为30，ceil(mc / mr) = 4，循环次数不足，因此不考虑在这个循环上并行。
• jr循环（选择）：nc为14400，nr为8，ceil(nc / nr) = 1800，循环次数足够多，考虑在这个循环上并行。
• ic循环（选择）：mc为120，m很大时，循环次数足够多，可以考虑在这个循环上并行。
• pc循环（不选择）：不考虑在这个循环上并行，理由如前面所说。
• jc循环（不选择）：因为本身Intel Xeon Phi没有L3 cache，所以在这个循环上并行的效果没有在jr循环上并行的效果好。

4. 核心内并行

因为核心内部的线程共享L1 cache，所以可以考虑使用同步（Synchronization）来控制核心内部线程的读数据，使得每次只需要读入一个元素数据，四个线程都可以一起使用，减少所需要的内存带宽。例如在jr循环上的并行，线程间共享A的block，那么访问到A的某个元素时，可以先让一个线程将其读入到L1 cache中，然后让其他的线程用过这个元素之后，才将其换出L1 cache。这里就需要一个线程间的同步来控制。

5. 核心间并行

因为每个核心都有自己的L2 cache，所以按照前面前面的说法，可以考虑在ic循环上实现核心间的并行。不过当m比较小时，ic循环上的并行效果较差，这时就可以考虑在jr循环上实现核心间的并行。这是一个动态调整的过程。

6. 性能结果

结果如上图所示。纵坐标是每个核心的GFLOPS。每个核心的计算峰值为17.6GFLOPS。每个核心使用了4个线程，总共有60个核心。因此总共有240个线程。

jr: 240 way表示将240个线程全部用到jr循环的并行上。从图中看效果并不好，原因有：

1. jr循环上的线程共享A的block，通过缓存一致性协议来进行共享。线程过多导致共享的开销过大（extra memory traffic）。
2. 线程过多，每个线程只是负责对A的block与少量B的sliver进行计算，每个线程的计算量较少，无法掩盖packing A（loading A）的开销。
3. A分块本身比较小，如果所有线程一起对A分块进行packing时，每个线程所负责packing的部分过小，并行效果自然就不好。

因此考虑使用ic : 60 way, jr : 4 way。在核心间使用ic循环上的并行。但是如果m不能被mc * thread_num整除时，会出现负载均衡的问题。因此并行策略进一步调整为ic : 15 way, jr : 16 way。从上图中可以看出，它的性能曲线更平滑一些。

图9是BLIS与MKL的性能比较。其中上图使用的并行策略是ic : 15 way, jr : 16 way，因为在m，n变化时，这种策略会表现得更好。而下图使用的并行策略是ic : 60 way, jr : 4 way，因为在m比较大时，这种并行策略会表现得更好。可以看出两者的性能差距并不大，这说明BLIS本身有足够好的扩展性。

同时从图9的下图可以看出，当K略大于kc（240）的倍数时，性能会下降。通过分析可以发现，性能的下降主要是因为register kernel计算更新C中mr x nr小分块时，每次都会剩余一小部分k未计算。这一部分的计算量较小，对C中mr x nr小分块进行更新（register kernel计算）的内存开销无法被足够的计算量所掩盖（更新mr x nr小分块的访存开销是固定的，mr x nr / 8），也就是计算访存比较低，从而导致了性能的下降。