程序性能优化探讨（5）——高速缓存、存储器山与矩阵乘法优化

        这一节内容将综合（3）和（4），讨论高速缓存相关的程序优化。

一、牛B完了的存储器山

        一个程序从存储系统中读数据的速率被称为读吞吐量或读带宽。如果一个程序在s秒的时间段内读n个字节，那么读吞吐量就是n/s，一般用MB/s作为单位。

        第（4）节中讨论的时间局部性和空间局部性。从缓存块大小来权衡，时间局部性和空间局部性似乎刚好成反比，但当我们全面的讨论整个存储层次性能时会发现，其他两种局部性可能各自有各自的规律。

        如果我们把一次性读取的数据量的大小称为工作集，把上一次和下一次读取数据的之间存储距离称为步长，如果我想获得工作集和步长这两种因素影响下程序的读吞吐量，会是什么结果呢？教材中给出实现代码，但我不打算在这里全部贴出来，因为详细的代码分析涉及到教材后面章节里的K次最优测量方法←_←内啥，我自己还没看懂，准确的说是暂时没耐着性子看下去，好像天书，等哪天有心情看懂了，再补上代码实现也不迟，那我们就先把存储器山的图贴出来吧：

        我们看到，存储器山由三个坐标确定点位，底面右边的坐标是工作集大小（wording set size），按字节计数；底面左边的坐标是步长（stride），按字计数。这里可以把字理解成4字节；垂直的坐标就是读吞吐量，根据吞吐量的大小判断程序性能的优劣。注意到这个存储器山针对的CPU型号已经写在上面了，这个CPU里面分别有大小2*32kb、块大小64字节的高速缓存L1，指令和数据缓存都是这么大。另外还有6M大小的高速缓存L2。

        首先把步长作为常数，来看工作集32k大小时的情况。一次性处理32k字节的数据，由于L1有32k的指令缓存和32k数据缓存，因此，32k工作集的数据是完全可以存入L1缓存的，因此可以看到此时它的吞吐量处于最优状态。奇怪的是，当工作集坐标往右减小时，山峰没有继续升高，而是急速降低，你可以明显的看到，从32k之后的下降陡坡，这是为什么？其实，这里需要了解测试程序的结构：

#define MINBYTES (1 << 12)  /* 最小工作集 ranges from 4 KB */

#define MAXBYTES (1 << 25)  /*最大工作集 ... up to 32 MB */
#define MAXSTRIDE 30           /*最大步长*/

    for (size = MAXBYTES; size >= MINBYTES; size >>= 1) {
for (stride = 1; stride <= MAXSTRIDE; stride++) {
   printf("%.1f\t", run(size, stride, Mhz));
}

        上面这个是生成存储器山的核心循环，第一层循环是遍历工作集大小，从32MB开始两倍变小，直到4kB为止；第二层循环是步进，从1到30字（4字节）。有了这两个测试关键参数传入存储器山生成函数run里，就能将这座山遍历出来。

        好了，当第一层size取值小于32k，而步长又大于20时，由于工作集size太小，因此每一个工作集处理时间会非常短，但正因为工作集太小，而步长又太大，数据集越不够连续，命中的情况就越遭，因此惩罚也就越大，主要的时间就消耗在不命中处罚和循环本身的开销上，因此山高峰背后那个明显下降，其实并不能反映出L1高速缓存的真实性能。

        我们再来看看这座山，很明显，工作集越小或者步长越小，对应的都是山峰越高，也就是吞吐量越好，从循环主题上看，无论是多大的步长，工作集越小，越可能在L1或者L2缓存中多停留，那么在第二层循环遍历时，现成缓存数据被重复调用的可能性就越大， 所以时间局部性就越好；（比如工作集大小为1，第一次读1，缓存123，第二次第三次就有现成的2、3用，但如果工作集大小为4，缓存装不下4个数，就算装下123下一次也用不上，那么每次都没现成的，并且还要往更慢的下层缓存要数据，当然时间局部性就很差了）

        另一方面，工作集不变时，你步长越小，在第二层遍历内部，上一次循环的数据被下一次循环共享的几率就越大，而步长越大，共享机会当然就越小，所以空间局部性就越好。（比如步长为2，那么第一次读1、2，缓存了1234，第二次读3、4时就有现成的；如果步长为4，那么每次几乎都没现成的了） 

        正因为有上面的分析，我们看到，空间局部性的斜坡相对连续，而时间局部性的斜坡就泾渭分明，体现了L1、L2和内存的读吞吐量之间巨大的性能差别。

        另外我们还注意到，L1数据缓存是32k，所以它在工作集超过32k时急速下降，符合预期。但L2大小是6M，那么我们预测的吞吐量也应该在工作集为6M时下降才对，但从我图上标注的两条线来看，下降处居然是工作集为4M的地方，why？按书上的解释是，L2不像L1那样把指令缓存和数据缓存独立区分，L2的数据和指令统一缓存在一起，因此虽然大小是6M，但你数据可不能独享这6M的空间，分2M给指令缓存也是可以理解的……

        还有个有趣的现象，我们单独来看L2区域的空间局部性。我们假定工作集大小固定在512kB，也就是L2区域的中间部位，我们看到，步长从0~16的部分，山体沿着步长坐标的方向有明显的下降趋势，这也符合之前对空间局部性的预期——步长越长，空间局部性越差！步长越长，被L1命中的可能性越低。但是从步长16开始往后走，坡度消失，几乎变成飞机场，这是为啥呢？原来，步长为16字，对应的就是64字节，这刚好是L1的块大小，也就是说，之前L1的一个块里就可能缓存了几个步进大小的数据，一次访问可能导致后面的数据在L1缓存块上命中，现在步进达到16，也就是64字节的步进，那么L1的每一个块都不可能再有命中的存在，都必须从L2得到服务，因此吞吐量完全取决于L2传送到L1的性能，因此保持不变。

        存储器山是反映特定系统的时间和空间局部性的山，对于高级别的程序员，一旦有了这样的信息，他就会尽可能的使得自己的程序中频繁使用的字是从L1中存取，同时还要让尽可能多的字从L1高速缓存中访问到。这就分别利用时间局部性和空间局部性。

二、简单应用

        作为关心性能的程序员，知道对存储器层次结构各个部分访问时间的粗略估计值是很重要的。根据上面这匹山，估计下列这些位置读出一个4字节所需的时间，以CPU周期为单位（T9300的频率是2.5GHz）：

        1、在芯片上的L1 d-cash；

        2、在芯片上的L2 d-cash；

        3、在主存上（工作集大小16M，步长=16），读吞吐率为80MB/s

        为什么要讲这个看似简单的题呢？也许从中能发现概念上的一个模糊点，至少我模糊得很，如果你也模糊，那恰好可以跟着理一理！

        第一个问， L1上读取4字节的时间，以CPU周期为单位。我们看L1的峰值吞吐率是10000MB/s，也就是10GB/s，而CPU的频率是2.5GHz，因此4个字节的读取时间就是（2.5/10）*4 = 1周期，也就是接近一周期。

        首先明确下，关于G和M的参数描述中是（1GHZ=10^3MHZ=10^6KHZ= 10^9HZ）的关系，而不是1024。OK，关于上面这个式子，你完全理解透了没？反正我刚开始是一头雾水，这里详细纠结下周期和频率的概念。

        下面是我面对纠结时的思考步骤：

                ①：CPU频率是2.5GHz，说明一秒钟运行2.5G次，那么一个周期就是1/2.5G秒，∵L1的峰值吞吐率是10GB每秒，现在CPU一个周期又是1/2.5G秒，那么一个周期的吞吐量就是两者相乘：10/2.5=4B，也就是说，CPU一个周期内的吞吐量是4字节——我去，这典型的是撞上答案！（当然，如果算出来是8B，我会用8/4算出2周期答案）

                ②：上面是先算出单位周期的吞吐量，然后与4字节相除，得出周期数，现在从频率出发来考虑：既然L1吞吐率是10GB每秒，那么每一个字节的耗时就是（1/10G）秒，4个字节耗时自然就是（4/10G）秒，既然算出了总耗时，接着就是把它转化成CPU周期数就OK了。怎么转换呢？考虑到CPU频率是2.5GHz，说明一秒钟运行2.5G次，现在我有（4/10G）秒，两者相乘，就能算出总的运行次数：(2.5G/10G)*4 = 1次，这意味着什么？意味着（4/10G）秒的这个时间量，刚好是CPU运行一次所需的时间，也就是CPU的周期！因此得出答案1周期！（当然，如果算出来是2次，那就说明答案是2周期！）

                ③：以上两种方法都得出正确答案，但有存在两个问题，其一，不可能每次算个时间周期数都要整这么麻烦；其二，感觉没有彻彻底底的理解清楚什么吞吐率频率和周期之间的关系，有云里雾里的错觉。那么接下来我就好好的整理了一下概念。

        周期T：每次循环的消耗多少时间、每次多少秒、Ts/次

        频率f： 每个单位时间多少次循环、每秒多少次 、f次/s

        上面是我根据以往常识总结的不太严谨但够用的概念，根据这两个概念可以得知，所谓的周期，它的单位是（秒/次），描述的是每次循环的耗时；所谓频率，单位是（次/s），描述每秒包含多少次循环，他们的单位和丈量标准刚好是相反的倒数关系，因此有f=1/T，T=1/f的换算公式，具体怎么理解这两个换算公式呢？我们先看第一个公式：

        f = 1/T，完整的写就是f(次/s) = 1(s)/T(s/次)，我们发现，T表示的是每次多少秒，它虽然是个时间概念，但它被限制在仅仅一次循环范围内的耗时秒数，（s/次）可以叫做"单位次时间"。好了，有了这个"单位次时间"，现在我关心对于一个确定的周期T，在总共1s的时间内会循环多少次呢？既然T衡量单位次（也就是一次）消耗的时间，现在用1s除以T，就相当于这个1s总时间被“单位次时间”砍成一节一节的段（一个总时间被一个特殊的时间量除），每个段都代表1次循环的时间量，有多少段就代表有多少次循环，好，这1s被砍了多少刀，也就是1s内有多少次循环，这恰恰就是每秒多少次——频率的概念！

        再来看1(s)/T(s/次)这个运算，“次”是分母T的单位中的分母，当两个s被约掉，“次”这个单位就会跑到分子上，最终答案就是(1/T)次，而我们的总时间是1s，因此也就是f(次/s)。注意，(1/T)计算结果本身的单位是“次”，只是因为分子刚好是1s，在定义频率概念时，就有了f(次/s) = 1/T。

        接下来是T = 1/f，完整的写就是T(s/次) = 1次/f(次/s)，f表示的是一秒钟有多少次，虽然它有"次"的概念，但它是被限制在仅仅1秒中内的次数，（次/s）可以叫做“单位时间次数”，好了，有个了这个“单位时间次数”，对于一个确定的频率f，在总共一次循环内会消耗多少秒时间呢？既然f衡量单位时间（也就是一秒）循环的次数，现在用1次除以f，就相当于这个1次的总循环次数被“单位时间次数”砍成一节一节的段（一个总次数被一个特殊的次数量除），每个段都代表1秒时间的次数，有多少段就代表有多少秒，好，这1次循环被砍了多少刀，也就是1次循环内有多少秒，这恰恰就是每次多少秒——周期的概念！

        同样的，1(次)/f(次/s)这运算，结果是(1/f)s，只是由于分子刚好是1次，所以周期T的单位才是(s/次)。这里注意的是，如果不止一次，那单位就应该是s了！

        有个问题，1s对应多少个周期呢？简单，用总时间除以单位次时间，就能得出周期数，也就是1/T，哈哈，恰恰就是我们f!，也就是说，1s对应的周期数是f，那2s对应的周期数是2f，ns对应的周期数就是nf！所以频率还可以理解成单位秒周期数。同理，1次对应多少时间呢？用总次数除以单位时间次，就能得出时间数，也就是1/f，恰恰就是周期T！也就是说，1次循环对应时间是T，那2次对应的时间就是2T，n对应的时间就是nT，所以周期还可以理解成单位周期时间数……这不废话么……

        好了，有了上面冗长的推理，终于对所谓的单位秒，单位次，单位时间等概念彻底理清了，物理中的很多除法公式也是类似，下面来看题。很明显，（2.5/10）*4 = 1周期这个式子是先计算出来一个字节的传输周期数，然后再计算4字节。怎么理解（2.5/10）计算的是传输1字节的周期数呢？有两种思维模式。其一，既然吞吐率是10GB/s，说明传输每个字节的时间就是1/10G，根据上面的结论，有了时间n，周期数就是nf，于是就有2.5G*(1/10G)；其二，既然频率是2.5GBHz，根据上面的结论，说明每秒钟的周期数就是1*f也就是2.5G，既然1秒的周期数有2.5G，现在要处理1B，但每秒钟的吞吐率是10GB，根本用不了1秒这么长的时间，因此用总周期数除以多余的10GB，2.5G/10G，得出的就是1B所需的周期数了……

         剩下斜字的部分属于本人任性的部分，闲着没事干可以细读。

        说到这，我再用自己的理解，来阐述下除法和分数的内涵。除法的本质到底是什么？比如8÷2 = 4，有两种理解方式：

        1、8本来是个整体，现在用长度为2的模具，从头切到尾，发现刚好能切成4份长度为模具长度的块，因此除法的结果是2——除法的结果，就是描述以模具长度标准来切割被除数后，所切下的模具长度的块的个数！那么3÷2呢？整体3被模具2切下一块后，剩余的1就不够切了，于是整数位还是模具长度的块的个数1，而剩下的1仍然要用来描述模具长度的个数，不够一块，那就有了1/2的概念，于是1.5的答案仍然是描述切割后模具长度的块的个数。

        2、8本来是个整体，被均分成两份，每份就是4这么多，因此除法的结果是4——除法的结果，就是描述以特定份数对目标数进行均匀切割后，每份的数量。那么1÷2呢？均分无法用整数实现，于是1/2就是答案本身。

        观察两种对除法的定义可以发现，定义1中的除法结果，是份，切割后剩余的份数；而定义2中的除法结果，是数量本身，和被除数一样。

        用定义1再来看1/2，我们发现，除法还可以理解成，将被除数这个整体，均分成除数份，比如1/2就是将1等分成两份，有了这两种理解方式？就可以解释3/(1/2)了：以（1/2）为模具，切割3，我们发现这个模具比我们长度为1的模具还要小一半，因此切割出来的个数肯定更大，相当于6个1/2的数量，因此答案就是6，这个例子相当于再描述：先把模具切小，然后再用来切除数！

        我们熟悉的路程S单位是m，时间T单位是s，速度v单位是m/s。这个速度单位明显就是个除法，表示1m/1s。如果S是6m，v是2(m/s)，时间是多少呢？答案简单是3s，这里面是有玄机的，为什么除法的结果是s这个单位呢？有两种解释，一种解释就是纯代数运算，6m/2m/s，可以换算成(6m/2m)*s，s作为分母的分母，自然可以转换到分子上来，而剩下的m可以通过约分消掉，得出答案……关键是我们该怎么去理解这个运算的本质呢？为啥一个m单位的量除以一个m/s单位的量，结果就是s？m/s到底是个啥东西？

        可以观察到，m/s是1m/1s这样的一个除法式子，1m被1s这个模子切割，切得动么？当然切不动，就像1/3也是切不动那样，所以m/s就是最终结果。好关键在于这种复合单位如何做除法，它背后的意义是什么？要理解这个类似除法的式子m/s，我们不妨换一个角度：

        一般的自然数1、2、3……如果代表数量，我们默认其单位为（1数量），比如7就代表7个（1数量）。现在来看这个复合除法：8/(6/3)，有两种方式去理解：

        1、把分母的6看成是单位为（1/3）的量，把分子8看成是单位为（1）的量，于是就有了8(1)/6(1/3)，由于除法结果的单位是在分子，因此有必要把分母的单位转换成（1），那么依靠通分转化：8*3(1)/6(1)——8(3)/6(1) = 4/3（3），我们得出4/3的单位为（3）！

        2、先把分母约分成2，此时2的单位已经是（1/3），我们再把分子单位进行转换：8（1）——8/3（3），于是就有(8/3)（3）/2（1）= 4/3（3），结果相等。

        为什么我一定要特别把答案构造成（3）单位呢？这是为了体现（1/3）和（3）这两个单位的在除法中的转换关系。我们按照上面1的理解，分母6的单位是（1/3），经过转换，将分母单位转换成（1），分子的单位变换成（3），出来的结果也是（3），这是不是很像路程和速度的计算公式(6m)/(2m)/s？原来分母的单位是（1/s），经过计算，单位1/s转换到分子变成单位s，最后计算出来的结果也是就是s了！

       

        我们再回到上面的题目，周期的概念是s/次，频率的概念是次/s，我们要计算耗时是多少周期，就是要按照Ts/次，把实际时间分割分割成一段一段的Ts，多少份就是多少个周期，然而我们发现，Ts是每次循环的时间，如果事先已经知道是多少次循环，那根本不需要关心实际时间。那到底是多少次呢？这道题的条件是频率500M次/s，吞吐率1000MB/s，既然都是按“每秒”来衡量，每秒500M次，每秒吞吐量1000MB，那么要计算1B（字节）所执行的次数，就自然有500M次/s/1000MB/s = (1/2)次/B，每个字节消耗1/2次，也就是1/2周期，那4个字节当然就是2周期。事实上，如果按照另一个算法，每个字节所耗时间(1/1000)（s/B），每个周期（每次）耗时(1/500)（s/次），你把两者一除会发现，答案的单位仍然是次/B。

三、矩阵乘法的优化

        这里我们讨论n×n矩阵的乘法：E=AB，比如n=3时，那么：

    e00e01e02e10e11e12e20e21e22

=  

a00a01a02a10a11a12a20a21a22
×

b00b01b02b10b11b12b20b21b22

e00 = a00×b00 + a01×b10 + a02×b20

e11 = a00×b01 + a01×b11 + a02×b21

……

        一般来说，第一个下标表示行，递增方向是从上到下；第二个下标表示列，递增方向是从左到右。对应C语言的二维数组，则访问顺序是先遍历每行的各列，然后再跳到下一行。

        总之，矩阵乘法的本质就是，目标矩阵的Exy的值，等于Ax(0~n)与B(0~n)y各自相乘后的加法结果。可以理解成A的x这行一排数，和B的y这一列数，拼成个十字架扔到Exy这个空空里去，而十字架的运算方式就是先乘再加，我们暂且称它为“交叉乘加”。

        英文里行是row，也就是常说的排，column是列，也就是常说的纵队。因此要实现矩阵乘法，分别用r、c来表示行列，用k来表示递增变量，出现六个版本，下面参看第一个版本：

typedef double array[MAXN][MAXN]

void rck(array A, array B, array E, int n) 
{
    int r, c, k;
    double sum;

for (r = 0; r < n; r++) 
    for (c = 0; c < n; c++) {
sum = 0.0;
for (k = 0; k < n; k++)
   sum += A[r][k]*B[k][c];
E[r][c] += sum;
    }

}

        这个版本应该是最直观明了的版本了。ABE都定义成数组，n是矩阵的位数。遇上多重循环，我本人习惯从最里面开始看起，sum += A[r][k]*B[k][c]，很明显，是在计算上面的矩阵乘法，某个E的结果是多个乘数相加的结果，而sum += A[r][k]*B[k][c]就是在做其中一个乘法。k是遍历n位，也就是说，把A确定的r行里的每一列数遍历完，同理B[k][c]是在确定的列纵队c中遍历每一排。两个十字交叉完成，于是对应的E[r][c]位置也就有了它该有的值。

        接下来是第二层循环，对于列c的遍历，重复里层循环，那么A完全相同的一行数据与B不同的列的数据进行交叉乘加，得出不同的E[r][c]，好了，当这层训循环结束时，B的所有列都和A的某一个行交叉乘加完成，此时会跳出到第一层循环，也就是遍历r的步骤。当A的每一行都重复之前同一行A数据交叉乘加B的每一列时，E[r][c]的每一个元素都被赋值，整个矩阵乘法结束。

        接下来看第二个版本：

vord crk(array A, array B, array E, int n) 
{
    int r, c, k;
    double sum;

for (c = 0; c < n; c++) 
    for (r = 0; r < n; r++) {
sum = 0.0;
for (k = 0; k < n; k++)
   sum += A[r][k]*B[k][c];
E[r][c] += sum;
    }
}

        很明显里层循环没有动，外层循环颠倒了顺序。具体分析方法和rck版本类似，先用A的每一行去交叉乘加B的某一个列，完成后调到第一层循环，变换B的列，然后重复里层循环，当B的列也被遍历完时，E[r][c]的每一个元素都被赋值，整个矩阵乘法结束。接着看第三个版本：

vord rkc(array A, array B, array E, int n) 
{
    int r, c, k;
    double m;
    
    for (r = 0; r < n; r++)
        for (k = 0; k < n; k++) {
            m = A[r][k];
            for (c = 0; c < n; c++)
                E[r][c] += m*B[k][c];
        }
}

        这个版本的算法有了明显的不同，我们先分析3×3矩阵，下图描绘的是r=0时的，第二三层循环遍历访问到的ABE各元素区域：

a00a01a02 b00b01b02    b10b11b12    b20b21b22
e00e01e02      

        e00 = a00×b00 + a01×b10 + a02×b20

        e01 = a00×b01 + a01×b11 + a02×b21

        e02 = a00×b02 + a01×b12 + a02×b22

        e10 = a10×b00 + a11×b10 + a12×b20

        r=0时，目标矩阵计算的就是第一排数据E[0][c]，由于e00~e02三个元素都是由B的第一排a00~a02去交叉乘加B的三个列，因此a00~a02三个元素各自都会被调用三次，该算法就索性让他们只调用一次。基本思想是让e00~e02分成三次计算，每次只计算每个e0y式子里的其中一个乘法，由里层循环实现e00~e02的单次乘法，再由第二层循环变更a0y，重复里层运算，直到每个e0y交叉乘加的三个乘法加拼装计算完成时，e00~e02就计算完成。根据上面的式子，相当于在计算e同行元素结果时，把交叉乘加式子整体从左到右竖着遍历。

e的其余行以此类推，无非就是遍历a1y和a2y。接着看krc版本：

void krc(array A, array B, array E, int n)
{
    int r, c, k;
    double m;
    for (k = 0; k < n; k++)
       for (r = 0; r < n; r++) {
    m = A[r][k];
    for (c = 0; c < n; c++)
        E[r][c] += m*B[k][c];
       }
}

        咋一看和rkc版本没区别，认真看发现第二层和第一层的循环顺序变了。那么这我们假设3×3，k=0的情形：

a00   b00b01b02a10      a20      
e00e01e02e10e11e12e20e21e22
        e00 = a00×b00 + a01×b10 + a02×b20

        e01 = a00×b01 + a01×b11 + a02×b21

        e02 = a00×b02 + a01×b12 + a02×b22

        e10 = a10×b00 + a11×b10 + a12×b20

        e11 = a10×b01 + a11×b11 + a12×b21

        e12 = a10×b02 + a11×b12 + a12×b22

        ……

        e22 = a20×b02 + a21×b12 + a22×b22

        k=0时，第一层循环内已经把e00~e22所有元素过了一遍，由于受k=0的限制，在第二层循环内A遍历的是第一列，第三层循环内B遍历的是第一行。这里很明显，k=0时，里层循环的每一步乘法运算，计算的都是exy是交叉乘加的第一个乘法，如上面的式子，相当于把所有exy交叉乘加式子按E的行顺序遍历完，r每递增一次就遍历三个exy。当第一层循环k递增时，A的列和B的行随即递增，计算上面式子的第二竖。

        这个算法实质上是把交叉乘加的这个十字架，彻底拆分，但基本思想和上面类似，都是要实现axy的多次利用，只是遍历的方向不同罢了。



void kcr(array A, array B, array E, int n)
{
    int r, c, k;
    double m;

    for (k = 0; k < n; k++)
        for (c = 0; c < n; c++) {
    m = B[k][c];
    for (r = 0; r < n; r++)
        E[r][c] += A[r][k]*m;
    }
}

        咋一看就是把AB调换位置，但实际上把cr调换的位置，假设3×3，k=0的情形：

a00   b00b01b02a10      a20      
e00e01e02e10e11e12e20e21e22

        我们发现，kcr版本和krc版本在元素遍历区域上几乎一样，只是里层循环遍历的是a00~a20，第二层循环遍历的是b00~b02，光从对交叉乘加的拆分思想上来看，几乎和kcr完全一样。只是对E而言你，遍历顺序变成了e00 、e10 、e20 、e01、e11……也就是按E的列顺序遍历。

void ckr(array A, array B, array E, int n)
{
    int r, c, k;
    double m;

    for (c = 0; c < n; c++)
        for (k = 0; k < n; k++) {
   m = B[k][c];
    for (r = 0; r < n; r++)
        E[r][c] += A[r][k]*m;
        }
}

        ck再次交换顺序，假设3×3，c=0的情形：

a00a01a02 b00  a10a11a12 b10  a20a21a22 b20  
e00  e10  e20  

        e00 = a00×b00 + a01×b10 + a02×b20

        e10 = a10×b00 + a11×b10 + a12×b20

        e20 = a20×b00 + a21×b10 + a22×b20

        e01 = a00×b01 + a01×b11 + a02×b21

        e11 = a10×b01 + a11×b11 + a12×b21

        e21 = a20×b01 + a21×b11 + a22×b21

        终于到了一轮外层循环就遍历完所有A元素的算法了。我们看到，c每递增一次，只能计算出E的一个列，而里层循环遍历的是A的一个列，第二层循环遍历的是b的一个列，因此从拆分交叉乘加的角度来看，里层循环拆分了十字架的横杠，第二层循环拆分的是竖杠。在看上面的式子，按E的列顺序，里层循环顺序是竖着的，而第二层循环顺序是横着的。

        好了，六种算法全部展示，接着就该分析性能了。很明显，最平凡调用的是里层循环，因此我们就以里层循环作为分析对象，统计不命中率：

版本

每次循环加载

每次循环的存储

每次循环A不命中率

每次循环B不命中率

每次循环E不命中率

每次循环的总不命中率

rck&crk

2

0

0.25

1.00

0.00

1.25

ckr&kcr

2

1

1.00

0.00

1.00

2.00

krc&rkc

2

1

0.00

0.25

0.25

0.50

        注意到六个版本可以划分成三个等价类，划分原则是里层循环的内容完全一致。我们这里有个前提：

        1、数组元素是double类型，并且sizeof(double)=8，

        2、只有一个高速缓存，块大小为32字节（B=32）

        3、n的实际值很大，使得矩阵的一行不能完全装进L1高速缓存中“中的一块或一行”

        4、所有局部变量都被编译器安排到寄存器中存储，不存在为局部变量本身进行加载和存储指令。

    “5”、

        （引号部分是教材上漏写的部分，很关键，很容易引起歧义）

    好了，先来分析rck&crk,他们的里层循环部分都是

    for (k = 0; k < n; k++)
       sum += A[r][k]*B[k][c];

    这里可以明显看出，对A的访问是线性的，利用了空间局部性，由于元素大小是8字节，L1缓存块大小又是32字节，因此对A的访问不命中率是0.25（具体原因在本章第(4)节有详细讲解）。而对于B的访问，由于是跳着行访问，是非线性的，而且n足够大L1缓存不可能装下B的整个一行，因此其不命中率就是1.0，所以rck&crk里层循环的总不命中率是1.25。

    接下来分析ckr&kcr，他们的里层循环部分都是

        for (r = 0; r < n; r++)
            E[r][c] += A[r][k]*r;

        显然，这里每次对A的访问仍然是跨行的。因此A的不命中率是1.0。同理，里层循环还用对E进行应用，由于变化的是r，因此E也是跨行访问的，命中率仍然是1.0，所以ckr&kcr的里层循环总不命中率是2.0

 

        最后分析krc&rkc，他们的里层循环部分都是

            for (c = 0; c < n; c++)
                E[r][c] += m*B[k][c];

        用同样分分析方法得出，B和E的访问不命中率都是0.25，因此里层循环总不命中率是0.5。

        注意到“每次循环加载“这一项，其实就是类似A[r][k]或B[r][k]的读取操作，当然都是两次。而”每次循环的存储“则是对E[r][c]的写入。由于rck&crk版本的里层循环都由临时变量sum存储加值，不会产生对E[r][c]的写入操作（存储），因此是0。

        好了，光从不命中率上来看，胜负已分，性能关系是krc&rkc > rck&crk > ckr&kcr，实际情况是不是这样呢？

        下面我要进行实际测试数据，我的CPU是酷睿2 p8400，查询参数：

        http://www.cpu-world.com/CPUs/Core_2/Intel-Core%202%20Duo%20Mobile%20P8400%20AV80577SH0513M.html

查询到L1高速缓存的资料，如果只考虑单核，数据缓存总大小C=32KB，块大小B=64字节，相联度E=4，组数就是S = 32KB/(64B*4) = 2^15/2^8 = 2^7 =128，组数有128组，因此这款CPU的单核L1高速数据缓存参数（S,E,B,m）= （128,4,64,32）。

Level 1 cache size  ? 2 x 32 KB 8-way set associative instruction caches
2 x 32 KB 8-way set associative write-back data cachesCache:L1 dataL1 instructionL2Size:2 x 32 KB2 x 32 KB3 MBAssociativity:8-way set
associative8-way set
associative12-way set
associativeLine size:64 bytes64 bytes64 bytesComments:Direct-mappedDirect-mappedNon-inclusive
Direct-mapped
Shared between all cores
        如果你不能根据这些资料熟练的推到并理解（S,E,B,m）= （128,4,64,32）这个结果，那还是回去看下本章的第（3）节内容，那里有详细分析。

        总之我么得到B = 64，和我们之前假设的32有出入，好吧，那不命中率等于0.25的地方事实上就得换成0.125（理由懒得解释了，有本章前面的知识铺垫），那么三个版本的不命中率分别应该是1.125、2和0.25。同时，为了增加n使得A、B的一行数据足够大到不能完全装进L1高速缓存"里的一块或者说一行"好了，现在块大小B = 64B，而数组元素sizeof(double) = 8B，那么n必须大于64B/8B = 8。

    好吧，按照教材的实现方式，n从25到400，以25递增，对6个版本进行性能测试，计算每次循环所消耗的周期数。运行结果如下：

[root@localhost matmult]# ./mm

matmult cycles/loop iteration

  n   rck   crk   ckr   kcr   krc   rkc

 25  0.13  0.46  0.23  0.02  0.01  0.00 

 50  0.08  0.06  0.05  0.01  0.04  0.01 

 75  0.04  0.01  0.03  0.02  0.04  0.03 

100  2.22  0.07  1.44  2.89  0.09  4.46 

125  0.72  2.30  3.04  5.77  0.71  1.91 

150  1.88  1.34  4.57  6.54  1.78  1.97 

175  1.56  1.66  9.69  8.43  1.95  2.36 

200  2.81  1.86 13.45 11.29  2.50  2.53 

225  2.70  2.70 18.04 14.86  2.70  2.90 

250  2.31  3.17 19.28 16.67  2.84  2.26 

275  3.73  3.55 19.64 17.28  3.15  2.91 

300  4.76  5.26 20.42 19.08  2.65  2.84 

325  5.00  5.25 20.40 19.26  3.16  3.16 

350  5.65  5.44 20.48 19.44  3.32  2.79 

375  7.28  7.47 20.61 17.47  3.26  3.16 

400  7.97  8.06 20.92 17.76  3.03  2.74 

     看着是不是有点晕？那我们就把这组数据转换成性能统计图来分析更加直观：

    这个图直观的反映了六个版本的执行效率，纵坐标表示每次里层循环所需的CPU周期数，越高说明耗的周期数越多，性能也就越差。当n超过250时，我惊奇的发现，性能图非常鲜明的把六个版本分成三个梯队，krc&rkc性能最佳，rck&crk其次，ckr&kcr最糟糕，回过头去看我们对6个版本的总不命中率分析，他们的结果竟然相同！

    为啥我要说惊人和竟然呢？理论和实际结论相同很奇怪么？当然奇怪了。因为教科书上得出的结论和我稍有不同，在教材里性能最优的是rck&crk！这再次体现了尽信书不如无书。实际测试结果很有意思。当教材作者得出rck&crk性能优于krc&rkc的结论时给出了解释，认为不命中率并不说明一切，krc&rkc虽然有最少的不命中率，却有额外的存储器访问：krc&rkc的每次里层循环需要引用E[r][c]\A[r][k]\B[k][c]这三个存储器，而rck&crk只有引用A[r][k]\B[k][c]两个存储器，因此存储器引用影响了最终性能……而经过我的实际测试结果，是不是可以得出结论：这个存储器影响因素消失了？是不是由于CPU更高级，使得存储器引用的效率得以大幅提升？

    如果学过CPU内部原理应该有所了解，CPU在进行运算时，可能会经历取指、译码、执行、访存、写回、更新PC等步骤，而其中的访存阶段可以将数据写入存储器，或者从存储器读取数据。而这里面的所谓存储器，有可能是L1缓存，也有可能是L2缓存甚至有可能是L3、内存、硬盘，因此访存被认为是CPU执行指令的过程中最可能耗费更多时间的步骤。因此L1~L3甚至内存到硬盘的大小与读取速度可能直接影响访存的效率，我只能推断，我运行测试时的电脑，从L1~L3到内存到硬盘可能都全面由于作者的测试电脑，那我能不能尝试下获得整个我的电脑存储器层次的性能对测试结果的影响呢？这至少能涉及到从L1到内存的部分，好的，接下来再理一下我的虚拟机linux配置的配置：

L1：32kB

L2:3072kB

L3：p8400没有L3

内存：1024MB

    也就是说，当我对存储器的引用增加到E\A\B三个时，n的大小将决定我的数据能缓存到那一级。比如，我要限制在L1中，那么n最大就不能超过对(32k/(8*3))开方那么多次，也就是n<=37；如果我要限制在L2中，那就是限制n不超过对（3072k/(8*3)）开方那么多次，也就是n<=362;如果同理，要限制在内存中，n<=11585，内啥，n大于2048时，我的小本本已经快跑糊了，要不是果断shutdown估计可怜p8400就要报销了，因此我只能分别测试L1和L2，然后呢？n<=37时结论已经有了，我们发现krc&rkc的性能仍然完爆其他对手，而在75~275之间，rck&crk昂首挺胸；也就是说，在数据借助到L2缓存时，L2的本身的性能劣势使得访存耗时增加。而当大于275时，krc&rkc几乎已不可战胜，即便当n>362，缓存铁定到了内存级别时，也不能改变这个趋势。

    好了，为什么分界点不在362而是在275呢？我猜呢，注意哈，是我猜，没有充足依据，应该是数据缓存还要处理其他数据，不只是三个矩阵独享……