BLIS
本文在《Anatomy of High-Performance Matrix Multiplication》基础上实现多线程下的高性能矩阵乘法。对于尺寸为 $m\times k \times n$ 的 $C=AB+C$， 列主序算法具体如下：
$n_c$, $k_c$, $m_c$ 分别是 $n$, $k$, $m$ 三个维度上的分块大小。最外面两层循环选择矩阵A的某一列和矩阵B的某一行，第三层循环由 $i_c$ 对A的某一列的分块的索引，分块的大小为 $m_c \times k_c$。第四层循环选择 $\widetilde{B}$ 中的某一个小条（sliver），用 $j_r$ 索引。第五层循环选择 $\widetilde{A}$ 中的某一个小条，用 $i_r$ 索引。可以发现 $j_r$ 和 $i_r$ 这两层循环是在GotoBLAS中实现的inner-kernel。
其中，$\widetilde{A}$ 被放在L2 Cache，$\widetilde{B}$ 被放在L3 Cache，$\widetilde{B}$ 的一个小条（sliver）放在L1 Cache，$C_{aux}$ 的一个子块放在寄存器中。
OPPORTUNITIES FOR PARALLELISM
Parallelism within the micro-kernel
在micro-kernel内并行，即对计算 $m_r \times n_r$ 大小的分块并行化，文章提出这种并行方法太低效，原因有两点：
- 该过程需要多线程对寄存器内的结果进行累加，线程计算结果归并的开销太大。
- 计算量太小，不足以隐藏将寄存器内结果写回内存的开销。
Parallelizing the first loop around the micro-kernel
在 $i_r$ 循环上并行。每一个线程获取一个 $\widetilde{A}$ 的小条（蓝色），所有线程共享同一个 $\widetilde{B}$ 上的小条（红色），每个线程单独计算一个寄存器内的 $C_{aux}$ 的分块（绿色）。该并行方法需要考虑两个问题：
- 该层循环迭代次数固定为 $\lceil m_c/m_r \rceil$，在实际案例中，这种并行度很小。
- 该循环主要的访存开销为将 $\widetilde{B}$ 的一个小条从L3 Cache拷贝到L1 Cache，需要有足够计算量来隐藏这个开销，因此要确保 $\lceil m_c/m_r \rceil$ 足够大，这与上一点矛盾。
Parallelizing the second loop around the micro-kernel
在 $j_r$ 循环上并行。每一个线程获取 $\widetilde{B}$ 的一个小条，所有线程共享 $\widetilde{A}$。该层循环迭代次数为 $\lceil n_c/n_r \rceil$，在实际案例中，这个值通常会很大，因此该层循环并行度很高，可以为隐藏访存开销提供足够的计算量。由于 $\widetilde{A}$ 存在L2 Cache，需要考虑线程是否共享L2 Cache两种情况：
- 当线程共享L2 Cache时，一个 $\widetilde{A}$ 块和多个 $\widetilde{B}$ 的小条都存在L2 Cache种，此时可能需要调整 $\widetilde{A}$ 和 $\widetilde{B}$ 的尺寸。但实际情况中，$\widetilde{B}$ 的一个小条大小都比较小。
- 当线程不共享L2 Cache时，$\widetilde{A}$ 将被pack到多个L2 Cache中，如果pack过程时多线程并行的，那么每个L2 Cache都会存有 $\widetilde{A}$ 的一部分。而这些不同的 $\widetilde{A}$ 部分跨L2 Cache传输时将会受到缓存一致性协议的影响，从而产生额外的开销。这个过程可能与micro kernel的执行重叠。
Parallelizing the third loop around the inner-kernel
在 $i_c$ 循环上并行。每个线程获取一个不同的 $\widetilde{A}$ 块，所有线程共享一个 $\widetilde{B}$。每个线程执行GEBP_OPT1。需要注意当 $m < m_c \times number_ of_ threads $，线程将会负载不均衡。同时仍然需要关注L2 Cache是否共享的情况：
- 当线程共享L2 Cache时，所有线程获取的 $\widetilde{A}$ 存放在L2 Cache中，此时可能需要调整 $\widetilde{A}$ 的尺寸。
- 当线程不共享L2 Cache时，每个 $\widetilde{A}$ 块将会驻留在自己的缓存中。
- 当L3 Cache不共享时，每个L3 Cache将会保存 $\widetilde{B}$ 的一部分，此时需要依靠缓存一致性协议实现跨L3 Cache的数据拷贝。
Parallelizing the fourth loop around the inner-kernel
在 $p_c$ 循环上并行。每个线程获取矩阵A的一列和矩阵B的一行，计算得到矩阵C分块的结果，此时将涉及到多个线程对C的同一子块进行累加，导致线程之间的竞争，或者线程之间将计算结果进行归并。因此一般不在这层循环进行并行，除非以下两种情况：
- 当C很小时，只有通过并行化这个循环才能达到较高的并行性能
- 线程计算结果的归并开销相对于其它计算开销更小
Parallelizing the outer-most loop
在 $j_c$ 循环上并行。每个线程获取矩阵B的一个子块，所有线程共享整个矩阵A。当L3共享时，多个矩阵B的子块将驻留在L3 Cache中，受到空间的限制，矩阵B的尺寸不能太大。这种并行方法适合multi-socket架构，它一般支持NUMA机制，即每一个节点都独立的L3 Cache用于存放矩阵B的块。
INTEL XEON PHI
Architectural Details
Xeon Phi处理器有60个core，每个core有512 KB L2 Cache和32 KB L1 Cache。每个core有四个硬件线程，共享L1 Cache。一个core能够利用两个流水线，在每个时钟周期内调度两条指令。其中一个可用于执行向量浮点指令或向量内存指令，另一个可能只用于执行标量指令或预取指令。一个线程每隔一个时钟周期只能向每个管道发出一条指令。
The BLIS implementation on the Intel Xeon Phi
在Intel Xeon架构中，四个硬件线程共享一个L1 Cache，如果在m和n两个维度上并行，则必须要有2个 $\widetilde{A}$ 的小条和2个 $\widetilde{B}$ 的小条驻留在L1 Cache。为了满足L1 Cache的空间限制，需要减小 $k_c$。这又会导致micro kernel的计算量变小，不够分摊将 $m_r \times n_r$ 的 $C_{aux}$ 子块写回内存的开销。为了解决这一问题，文章给出的方法是将一部分L2 Cache视为虚拟的L1 Cache，其访存开销与L2 Cache的相同。
文章参数设置为 $m_r=30$, $n_r=8$, $m_c=120$, $k_c=240$, $n_c=14400$。
The BLIS implementation on the Intel Xeon Phi
文章利用fork-join模型对多个循环进行并行。
- $i_r$ 循环：$m_c/m_r=120/30=4$，不适合并行化。
- $j_r$ 循环：$n_c/n_r$ 较大，并行度很高，适合并行。四个硬件线程共享一个L2 Cache中的 $\widetilde{A}$ 块。
- $i_c$ 循环：该循环步长为 $m_c=120$，当m很大时，并行度将很高。同时，如果每个core的四个线程不在该循环上并行，那么一个L2 Cache将保留1个 $\widetilde{A}$ 块，因此无需修改 $\widetilde{A}$ 的尺寸。
- $p_c$ 循环：上文已经谈到该循环不适合并行。
- $j_c$ 循环：由于Intel Xeon没有L3 Cache，所有矩阵B块的拷贝都将保存在内存中，因此该层循环不适合并行处理。
Parallelism within cores And Parallelism between cores
从上文的分析得出，在单个core内四个线程并行执行 $j_r$ 循环，四个线程共享L2 Cache中的 $\widetilde{A}$ 块。如果这四个线程是同步执行的，那么它们将并发的访问 $\widetilde{A}$ 的一个小条。这是将由最早调度的线程将 $\widetilde{A}$ 的一个小条读到L1 Cache中，当其他所有线程都使用完后，再替换出去，这样就降低了带宽。如果四个线程是异步执行的，那么也可以设计同步机制（如barrier）来达到相同的效果。
由于每个core都有独立的L2 Cache，core之间适合并行 $i_c$ 循环。但是当m较小时，并行性能较低，因此core之间将同时并行 $j_r$ 循环。
Performance results
从结果看出，当240（$60 \times 4$）个线程全部用来并行 $j_r$ 循环时（$j_r: \text{240 way}$），性能表现较差。主要原因有以下几点：
- 每个线程执行 $\widetilde{A}$ 与若干 $\widetilde{b}$ 的小条的计算，计算量较小，不足以平摊访存开销。
- 如果对 $\widetilde{A}$ 的pack操作并行化，每一个线程工作量很小。
- 所有线程都共享一个 $\widetilde{A}$ 块，多个core的独立L2 Cache将会保存 $\widetilde{A}$ 的一部分，缓存一致性将影响性能。
$i_c: \text{60 way}; j_r: \text{4 way}$ 表示60个core并行 $i_c$ 循环，每个core内的4个线程并行 $j_r$ 循环，曲线呈现锯齿状是出现了负载不均衡的问题。$i_c: \text{15 way}; j_r: \text{16 way}$ 是将4个core视为一体，16个线程并行 $j_r$ 循环，可以看到负载不均衡明显减弱。
将BLIS与MKL库相比，两者性能接近。可以发现当k略大于 $k_c=240$ 的倍数时，会有一个明显的性能下降。这是因为分块不均导致末尾有一个很小的C块会在一次新的迭代中被计算，计算量太小无法隐藏访存的开销。