利用矩阵分块优化算法
本文摘要（由AI生成）：
文章探讨了存储系统对程序执行时间的影响，特别强调了CPU cache在提高处理器性能中的关键作用。由于CPU与内存性能差距，CPU内部引入多级cache以减少访存暂停时间。通过软件优化技术，如改善时间局部性，可减少cache失效次数，提升访问效率。文章以数组处理为例，分析了按行或按列存储的局限，提出分块算法优化时间局部性。文章还讨论了cache容量对性能的影响，提出基于子矩阵计算的代码改写方法，通过引入块参数和函数内联进一步优化性能。总之，文章强调了CPU cache的重要性，并详细阐述了分块算法如何通过挖掘程序局部性来减少cache失效，提高程序执行效率，对数值计算库性能优化具有重要意义。
存储系统对于程序执行时间有显著影响。处理器由于访存导致的暂停时间受到失效率和失效代价的影响。众所周知，为了弥补CPU与内存两者之间的性能差异，就在CPU内部引入了CPU cache，也称高速缓存。CPU cache通常分为大小不等的三级缓存，分别是L1 cache、L2 cache和L3 cache。如图1所示。
▲ 图1 三级cache
许多软件优化技术通过重用cache中的数据来大幅度提高处理器性能，通过改善程序的时间局部性来提升访问效率。换言之，不要频繁替换cache中的数据。
处理数组时，如果能将数组元素按照访问顺序存放在存储器中，则能够获得性能上的好处。但是，假设同时处理多个数组，一些数组按行访问，一些数组按列访问。按行存储(称为行优先)或者按列存储(称为列优先)数组都不能解决问题，这是因为在程序的每个循环体中行访问和列访问同时会被使用到。因而，分块算法针对子矩阵(submatrice)或者数据块来进行操作，并不针对数组中完整的一行或一列进行操作。它的目标是，在替换之前对已在cache中的数据进行尽可能多的访问，这就是说，提高程序的时间局部性以减少cache失效。例如，在数值计算库Lapack的DGEMM中的内层循环中
for (int j = 0; j < n; ++j) {
double cij = C[i+j*n]; /* cij = C[i][j] */
for( int k = 0; k < n; k++ ) {
cij += A[i+k*n] * B[k+j*n]; /* cij += A[i][k]*B[k][j] */
C[i+j*n] = cij; /* C[i][j] = cij */
}
}
首先读入数组B的所有 个元素，重复地读入数组A中某一行中的 个元素，最后将结果写入C数组某一行中对应的 个元素。图2给出了对这三个数组的访问快照，深色阴影部分表示最近被访问过，浅色阴影部分表示较早被访问过，白色表示还未被访问。
▲ 图2 三个数组 C、A和 B 的快照(N=6，i=1)。采用不同的阴影对各个数组元素的访问时机进行表示:白色表示尚未被访问过，浅色阴影表示较早被访问过，深色阴影表示最近被访问过。与图4 相比，数组 A和B的素被重复读入以计算数组 C的新元素。变量 i、j和 k 用来进行数组访问，对应行或者列的变化
容量失效的次数明显与 和cache的容量有关。如果cache中可以存入三个 的矩阵，假设没有其他冲突，那么一切完美。例如将矩阵大小设为32x32，那么情况正好如此。每个矩阵有32x32=1024个元素，每个元素的大小为8字节，三个矩阵的大小就为24KiB，这可以很容易存入的容量为32KiB的cache中。如果cache能够存入一个 的矩阵和一行 个元素，至少数组A中的第 行数据和数组B可以一直保留在cache中。如果cache容量小了，则数组B和数组C都可能发生失效。最坏情况下， 个操作需要访问 个存储字。
为保证需要访问的数组元素都尽可能在cache中，原始代码需要改写为基于子矩阵的计算方式。这样，我们需要调用图3中的DGEMM版本，该版本就是在大小为的矩阵上重复计算。 也被称为块参数。该函数增加了3个新参数 si、sj和sk用来描述数组A、B和C的子矩阵的起始点。函数doblock的两个内部循环以 为步长进行计算，并不是按照数组B和C的全长进行计算。gcc编译器通过函数内联(function inlining)消除了函数调用的所有开销。也就是说在程序中直接插入函数代码，以避免通常的参数传递和现场保存操作。
▲ 图3 DGEMM的cache分块版本。假设数组C初始化为0。函数doblock 以普通的DGEMM为基础，增加了新参数来描述BLOCKSIZE大小的子矩阵的起始位置。
图4给出使用分块思想对三个数组进行访问的示例。仅对于容量失效来说，需要访问的内存数据字总数为 。(相比未分块前)这个数据得到了改善，原因在于参数 。由此可见，分块思想挖掘出程序的时间局部性和空间局部性，比如数组A的访问得益于空间局部性，而数组B的访问得益于时间局部性。
▲ 图4 数组C、A和B的访问(BLOCKSIZE =3)。注意，与图2 相比，访问的元素数量变少
图5中给出了采用cache分块对未优化DGEMM性能的影响。矩阵规模最大时，未优化程序的性能折半。采用cache分块的版本，即使矩阵规模达到960x960=3232阵规模性能也仅仅降低了不到10%。
▲ 图5 未优化DGEMM与 cache 分块 DGEMM的性能比较，矩阵维度从32x32增加至960x960
