高速缓存与矩阵乘法(一)

来源：数据结构、算法与应用 C++语言描述（原书第2版）

简单计算机模型

我们来看一个简单的计算机模型，它的存储由一个一级缓存 L1（level 1）、一个二级缓存 L2 和主存构成。算术和逻辑操作由算术和逻辑单元（ALU）对存储在寄存器（R）中的数据进行处理来完成。下图是这个计算机模型的一部分。

通常，主存的大小是几十或几百 MB；二级缓存的大小不足 1MB；一级缓存的大小是几十 KB；寄存器的数量在 8 和 32 之间。程序开始运行时，所有数据都在主存。

要执行一个算术运算，例如加法，首先把相加的数据从主存移到寄存器，然后把寄存器的数据相加，最后把结果写入主存。

我们把寄存器的数据相加所需要的时间作为一个周期。把一级缓存的数据送到一个寄存器所需要的时间是两个周期。如果需要的数据没有在一级缓存，而是在二级缓存，即一级缓存未命中，那么把需要的数据从二级缓存送到一级缓存和寄存器需要 10 个周期。当需要的数据没有在二级缓存，即二级缓存未命中时，把需要的数据从主存复制到二级缓存、一级缓存和寄存器需要 100 个周期。我们把写操作，甚至向主存的写操作，算作一个周期，因为不需要等到写操作完成之后再进行下一个操作。

缓存未命中对运行时间的影响

在我们的简化计算机模型中，语句 a=b+c 编译后的机器指令是

load a;load b;add; store c;

其中，load 操作把数据送到寄存器，store 操作把相加后的结果送到主存。add 和 store 操作共需要两个周期。两个 load 操作可能需要 4 个周期至 200 个周期不等，这取决于数据是否在缓存中，即缓存是否命中。因此，语句 a=b+c 所需要的总时间从 6 个周期到 202 个周期不等。在实际操作中，时间差别没有这么极端，因为可以把连续的缓存未命中所花费的时间交叉处理。

假定有两个类型相同的算术运算。第一个算术运算是 2000 次加法，它需要 4000 次 load 操作、2000 次 add 操作和 2000 次 store 操作。第二个算术运算是 1000 次加法。第一个算术运算的数据访问有 25% 的 load 操作出现一级缓存未命中，另有 25% 的 load 操作出现二级缓存未命中。在我们这个简化的计算机模型中，第一个算术运算所需要的时间是 2000*2（有 50% 的 load 操作是缓存命中）+1000*10（有 25% 的 load 操作出现一级缓存未命中）+1000*100（有 25% 的 load 操作出现二级缓存未命中）+2000*1（用于 adds）+2000*1（用于 stores）=118000 个周期。如果第二个算术运算有 100% 的二级缓存未命中，它的用时 =2000*100（有 100% 的二级缓存未命中）+1000*1（adds）+1000*1（stores）=202000 个周期。第二个运算的量是第一个的一半，但实际用时却比第一个多 76%。

为了减少缓存未命中的数量，从而减少程序的运行时间，计算机采用了一些策略，比如，把最近需要处理的数据预载到缓存中，当出现一个缓存未命中时，把需要的数据和相邻字节中的数据装入缓存中。当连续的计算机操作使用的是相邻字节的数据时，这个策略很有效。

虽然我们的讨论集中在如何用缓存来减少访问数据的时间问题上，但是我们也用缓存来减少访问指令的时间。

矩阵乘法

也许有人不相信，在一台商用计算机上，一个操作多的程序可能比一个操作少的程序实际用时要少。本节就是要让这些人相信，确有此事。

我们从一个实际的程序入手，程序 2-22 如下

// 程序 2-22
template<class T>
void squareMatrixMultiply(T** a, T** b, T** c, int n)
{ // 将 n × n 矩阵 a 和 b 相乘得到矩阵 c
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
        {
            T sum = 0;
            for (int k = 0; k < n; k++)
                sum += a[i][k] * b[k][j];
            c[i][j] = sum;
        }
}

它是把两个用二维数组描述的方阵相乘。计算如下所示：
$$c[i][j]=\sum _{k=1}^{n}a[i][k]\ast b[k][j],1⩽i⩽n,1⩽j⩽n$$
（你不理解矩阵乘法没关系，这不影响你对我们要说明的问题的理解。）程序 2-22 是一段标准代码，你在很多教科书上都可以找到。程序 4-4 是另一段代码，它和程序 2-22 一样，产生一个二维数组 c。我们来观察程序 4-4。它有两层嵌套的 for 循环，这是程序 2-22 所没有的，这使它对数组 c 的索引处理得更多一些。其余的操作都一样。

// 程序 4-4
void fastSquareMatrixMultiply(int** a, int** b, int** c, int n)
{
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            c[i][j] = 0;
    
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            for (int k = 0; k < n; k++)
                c[i][j] += a[i][k] * b[k][j];
}

你会发现，把程序 4-4 的三层嵌套 for 循环重新排列一下顺序，结果是不变的。我们把程序 4-4 的嵌套循环顺序称为 ijk。当我们把第二层和第三层的 for 循环交换次序，我们得到的嵌套循环顺序是 ikj。一共有 3!=6 种嵌套循环顺序。由 6 种嵌套循环顺序分别生成的函数都以同样的数量执行每一种类型的操作。因此你也许认为这些函数所需的运行时间也是相同的。但是错了。改变了循环的次序，也就改变了数据访问模式，进而改变了缓冲区未命中的数量，最终影响了运行时间。

在 ijk 顺序中，数组 a 和 c 的元素是按行访问的，数组 b 的元素是按列访问的。因为同行的元素在存储中是相邻的，而同列的元素在存储中是分开的，所以当数组很大，以至三个数组不能同时存储在二级缓存 L2 中的时候，访问数组 b 可能导致很多二级缓存未命中的事件。在 ikj 的顺序中，数组 a、b 和 c 的元素是按行访问的，因此二级缓存未命中的事件就比较少，因此所需时间也比较少。

上图给出了程序 2-22 和程序 4-4 分别使用 ijk 顺序和 ikj 顺序时的运行时间（以秒为单位）。下图显示的是标准运行时间，即一个函数的运行时间除以在 ikj 顺序下执行的时间。

多么神啊！ikj 顺序要比 ijk 顺序和程序 2-22 运行快。实际上，当 n=500 时，ikj 顺序所需要的时间仅是 ijk 顺序的 1/3，是程序 2-22 的 1/2；当 n=1000 时，比率近似是 7/16 和 1/4；当 n=2000 时，比率近似是 1/13 和 1/16。记住，按操作步数计算，ikj 顺序比程序 2-22 和 ijk 顺序所执行的操作要多。只有 ikj 顺序的运行时间是按照渐进分析的比率 $\Theta (n^3)$增长的。ijk 顺序和程序 2-22 的运行时间是受缓冲未命中事件所控制的，而不受执行步数所控制。

存储等级制对代码性能的影响随着程序语言、编译器、编译器选项和计算机配置的变化而变化。例如，2.4GHz Intel Pentium IV PC 的二级缓存比 1.7GHz PC 的二级缓存大一倍，上述图中所给出的矩阵乘法的时间是用后者实验得来的，如果用前者实验，那么当 n=500 时，比率大约是 9/16 和 2/5；当 n=1000 时，比率大约是 1/2 和 1/3；当 n=2000 时，比率大约是 1/4 和 1/5。