深入理解计算机系统CacheLab-PartB实验报告

实验答案托管在我的GitHub上

Cache Lab的Part B是我卡了比较久的实验，在7月份做完Part A之后我卡在了Part B的第二个矩阵优化，之后进度一直缓慢。直到几天之前团队分享，我才把这个实验重新捡了回来，最终将第二个矩阵转置优化到了1500+的miss数（尽管仍然没有达到满分），后参考了网上的思路完成了Part B，总结如下。

实验简介

Cache Lab - Understanding Cache Memories主要是有关缓存的实验，对应于书本的第6章：存储器层次结构。主要阐明了缓存对于C语言程序的性能影响。

本实验的第二部分要求优化一个简单的矩阵变换函数，使其具有尽可能晓得缓存不命中数。

关于本实验的具体介绍详见实验讲义。

实验要求 - Part B

在Part B中你需要在trans.c中写一个矩阵转置函数，该函数要能有尽可能少的缓存不命中(miss)数。

令A代表一个矩阵，A(i, j)代表矩阵A的第i行第j列的元素。
那么，令B为矩阵A的转置，则对于A中的任意一个元素，满足A(i, j) = B(j, i)。

一个朴素的矩阵转置函数如下：

void trans(int M, int N, int A[N][M], int B[M][N]) {
    int i, j, tmp;

    for (i = 0; i < N; i++) {
        for (j = 0; j < M; j++) {
            tmp = A[i][j];
            B[j][i] = tmp;
        }
    }

}

该函数是正确的，但却并不是高效的，因为其访问模式导致了相当多的缓存不命中。

你的任务是完成一个相似的函数，并对于不同大小的矩阵块（32×32，64×64，61×67），最小化该函数的缓存不命中数。

对于Part B有以下的要求和限制：

你所写的程序在编译时不能有任何的警告
在每一个你所写的转置函数中，你最多只能定义12个int型的局部变量
你不能使用long型的变量或是使用任何的技巧来使得在一个变量中存入多余一个的值
你的转置函数不能使用递归
如果你使用辅助函数，那么在同一时刻，你的调用栈上也不能出现超过12个局部变量
你的转置函数不能修改矩阵A，但是，你能任意的修改矩阵B
你不能定义任何的数组或是任何使用malloc的变量

实验过程 - Part B

缓存简介

局部性

局部性分为时间局部性和空间局部性——在一个具有良好的时间局部性的程序中，被引用过一次的内存位置可能在不远的将来再被多次引用。在一个具有良好的空间局部性的程序中，如果一个内存位置被引用了，那么程序很可能在不远的将来引用附近的一个内存位置。

一般来说，有良好局部性的程序比局部性差的程序运行的更快。现代计算机系统的各个层次，设计都利用了局部性。

在硬件层，通过高速缓存存储器来保存最近被引用的代码和数据，提高对主存的访问速度。
在操作系统层面，将虚拟内存作为主存和磁盘间高速缓存。
在网络中，主机可以作为代理缓存经常被访问的网页，浏览器也可以缓存最近访问过的页面。

存储器层次结构

我们可以依次地组织不同的存储器，形成一个典型的存储器层次结构。如下图所示。

在该结构中，自高层向底层走，存储设备变得更慢、更便宜和更大。每一层的存储器都作为下一级存储器的缓存。

高速缓存

早期的计算机系统的存储器结构只有3层：CPU寄存器、DRAM主存储器和磁盘存储。随着CPU和主存之间不断加大的差距，系统设计者被迫在寄存器文件和主存之间加入了一个小的SRAM高速缓存处理器，称为L1高速缓存（一级缓存），L1高速缓存的访问速度几乎和寄存器一样快，典型地是约4个时钟周期。

随着CPU和主存之间的性能差距进一步增大，系统设计者在L1缓存和主存之间又加入了一个更大的高速缓存，称为L2高速缓存，可以在约10个时钟周期内被访问到。

有些现代系统还包括有一个更大的高速缓存，称为L3高速缓存，它位于L2高速缓存和主存之间，可以在大约50个时钟周期内被访问到。

组织结构

假定一个计算机系统，其每个存储器地址为m位，那么形成了M = 2 ^ m个不同的地址。

这样一个机器的高速缓存被组织成一个有S = 2 ^ s个高速缓存组（Cache Set）的数组，每个组包含E个高速缓存行（Cache Line）。每个行是由一个B = 2 ^ b字节的数据块（Block）组成的，一个有效位（valid bit）指明这个行是否包含有意义的信息。t = m - (b + s)个标记位（tag bit），它们唯一地标识存储在这个高速缓存行中的块。

一个通用的高速缓存组织结构可以用元组(S, E, B, m)来描述。
高速缓存的大小C指的是所有块（不包括标记位和有效位）的大小的和，C = S × E × B。

一个m位的地址L被组织为3个部分，分别是高t位的标记位，中间s位的组索引以及低b位的块偏移。

L中的s个组索引是一个到S个组的数组的索引，从0开始。其指示了该地址所指示的字应当存储到哪个组中。
L中的t个标记位指明了该组中的哪一行包含这个字。当且仅当有效位为1且行的标记位与地址的标记位相同时，组中的行才包含了这个字。
块偏移b指明了在B个字节的数据块中的字偏移。

高速缓存写

高速缓存写分为直写和写回。

直写，当缓存更新后立刻将相应的高速缓存块写回到下一层存储器中。直写的缺点在于每次写都会引起总线流量。
写回，尽可能地推迟更新，只有当替换算法需要驱逐更新过的块时，才执行写入操作。写回极大地降低了总线流量。但是高速缓存必须对于每一个高速缓存行维护一个额外的修改位（dirty bit）来表明这个高速缓存块是否会修改过。

实验思路

首先要明确，尽管矩阵的转置本身导致对于A矩阵（原始矩阵）的读和B矩阵（转置矩阵）的写不可能同时为连续的（即不可能同时存在连续读和连续写——对A矩阵行的连续读必然导致对B矩阵列的非连续写）。
但只要矩阵的大小小于缓存的总大小，那么在理想的情况下，在最初的强制不命中（即缓存为空导致的不命中）后，整个矩阵都会被加载进入缓存。在这之后的所有对于B矩阵的不连续写的引用都会命中。

在该实验中，缓存采用的是直接映射高速缓存，s = 5，b = 5，E = 1。对于该缓存，总共存在32个组，每个组共32个字节，可以装入8个int型变量，是非常有限的缓存，主要需要解决以下两个问题：

直接映射缓存所带来的冲突不命中。观察程序中矩阵存储的位置即可以发现，矩阵A和矩阵B的同一行实际上被映射到了同一个缓存组。当进行对角线的引用时，一定会发生缓存的冲突不命中。需要仔细地处理对角线上的元素。
所需优化的矩阵的总大小超出了缓存的总大小。必然导致程序的访存效率低下。

为了解决第一个问题，我们需要仔细地考虑对于矩阵元素访问的顺序。至于第二个问题，我们采用矩阵的分块技术来降低不命中数。

矩阵分块的目的在于将大的、不能完全加载进入缓存的大矩阵分块成小的、可以完全加载进入缓存的小矩阵块来处理。小矩阵块具有良好的局部性，性能显著增加。
但同时也要注意，分块使得程序的可阅读性大大降低，因此一般只在常用的库函数中采用分块优化。

第一部分 32 × 32矩阵优化

第一部分满分的要求是300个misses以内，misses超过600则0分。

首先对矩阵进行分块处理。为了完全利用每一个缓存快（32个字节）采用8 × 8分块。然后处理对角线的问题。这里我采用的方法是无论是哪一个矩阵分块，均从该矩阵分块的对角线开始处理。同时对于A矩阵（原始矩阵）按列优先（不连续读），对于B矩阵（转置矩阵）按行优先（连续写）。

通过优先处理对角线(a, a)的元素，保证了B矩阵的第a行被载入缓存中，接下来对于A矩阵的列优先处理保证了B矩阵的第a行缓存被充分利用。

对于32 × 32的矩阵，总共存在1024次读和1024次写。对于非对角线的分块（总共12个），其缓存不命中率是1/8（仅强制不命中），对于对角线的分块（总共4个），其写的缓存不命中率是1/8（强制不命中），其读的缓存不命中率为1/4（强制不命中和冲突不命中各一半）。

因此，理论上优化之后的总缓存不命中数为2048 × 0.75 × 0.125 + 1024 × 0.25 × 0.125 + 1024 × 0.25 × 0.25 = 288次。

第一部分优化之后的代码如下：

for (a = 0 ; a < N ; a += 8) {
    for (b = 0 ; b < M ; b += 8) {
        for (c = b ; c < b + 8 ; ++c) {
            for (d = a + c - b ; d >= a ; --d) {
                B[c][d] = A[d][c];
            }
            for (d = a + c - b + 1 ; d < a + 8 ; ++d) {
                B[c][d] = A[d][c];
            }
        }
    }
}

实际测试的缓存不命中数为287次，与理论值几乎一致。

第二部分 64 × 64矩阵优化

第二部分的满分要求是misses小于1300，当misses大于2000则零分。

第二部分对于misses的要求限制的非常严格，同时如果采用第一部分的8 × 8分块方式会出人意料地带来大量的misses。下面具体分析8 × 8分块导致misses增多的原因。

实验采用的缓存为直接映射高速缓存，s = 5， b = 5， E = 1。对于任意一个地址，其从低地址开始的第0-4位为块偏移b，第5-9为组索引s。

对于32 × 32的矩阵M来说，M[0][1]和M[1][1]之间总共间隔32个int型元素，也就是0x80个字节，也就是说，同一列相邻行的元素之间的地址间隔为0x80 = 0x100|00000。对于8 × 8的矩阵分块而言，其8行可以全部被加载进入缓存中而不发生任何冲突不命中。

然而，对于64 × 64的矩阵，其同一列相邻行的元素之间的地址间隔为0x100 = 0x1000|00000。对于8 × 8的矩阵分块而言，其第1、2、3、4行的元素会和第5、6、7、8行的元素占用相同的高速缓存组，进而出现严重的冲突不命中现象。

使用4 × 4的矩阵分块又无法充分利用每一个高速缓存行（32个字节=8个int数据），仍然无法达到所要求的misses数。

经过大量的尝试，我使用了以下的方法进行矩阵转置的优化：

仍然按照8 × 8对矩阵进行分块，只是在8 × 8的分块内部再按照4 × 4进一步分块，得到左上、右上、左下、右下4个子块。

紧接着依次按照左上、右上、右下、左下的方式处理4个子块（A矩阵）。
对于左上、右下这两个可能出现对角线元素的块，按照第一个矩阵优化的方式进行处理。
右上、左下子块不能简单地按照A矩阵不连续读，B矩阵连续写的方式处理。原因是对于对角线上的8 × 8分块来说，A、B矩阵的左上子块和右下子块占用了相同的高速缓存组，存在着严重的冲突不命中风险。
因此对于右上、左下的子块，我们按照下图的方式处理。

图中，上方的矩阵为矩阵A，下方的矩阵为矩阵B，小方块代表矩阵的元素，黑色方块的表明加载进入临时变量/已写入的元素。红色的线表明该行被缓存。

利用8个临时变量，将左上子块的前两行加载进入临时变量中，考虑到之前的缓存条件，该次加载的缓存命中。
将一个小的2×2的矩阵转置写入矩阵B右下子块的前两行，无论是否为对角线上的分块，该次写入一定会发生缓存不命中，同时将B矩阵的前两行载入高速缓存行。
将矩阵A左上子块的后两行的2×2的矩阵加载进入空闲的4个临时变量中，同之前加载相似，该次加载的缓存命中。
将刚刚加载的2×2矩阵转置写入B矩阵右下子块前两行的剩余位置，由于之前这两行已经加载进入了高速缓存行，故该次写入的缓存全部命中。
将矩阵A后两行的剩余元素加载进入空闲的4个临时变量中，缓存命中。
将8个临时变量中的元素依次转置写入矩阵B右下子块的最后两行中，同2相似，写入一定会发生缓存不命中，同时将矩阵B的后两行载入高速缓存行。

A矩阵左下子块的转置操作类似，这里不再赘述。

在A矩阵右上子块转置完成后，紧接着执行的是右下子块的转置，此时，对于非对角线上的分块而言，写入时的缓存必定命中。对于对角线上的分块，则会发生缓存不命中。

因此，对于64 × 64的矩阵，总共进行4096次读和4096次写，对于非对角线的分块（总共56个），对于A矩阵（原始矩阵而言），其左上、右下分块的不命中率为1/4，左下、右上分块的不命中率为0；对于B矩阵（转置矩阵而言），其左上、右上、左下分块的不命中率均为1/4，右下分块的不命中率为0。对于对角线上的矩阵，其B矩阵不命中率上升至1/4，对于A矩阵，其左上、右下的不命中率上升至1/2。

因此，理论上优化之后的总缓存不命中数为4096 × (8/64) × (1/2 × 1/4 + 0 × 1/4 + 0 × 1/4 + 1/2 × 1/4) + 4096 × (8/64) × 1/4 + 4096 × (56/64) × (1/4 × 1/4 + 0 × 1/4 + 0 × 1/4 + 1/4 × 1/4) + 4096 × (56/64) × (1/4 × 1/4 + 1/4 × 1/4 + 0 × 1/4 + 1/4 × 1/4) = 1376次。

该优化方法对应的代码如下：

for (a = 0 ; a < M ; a += 8) {
    for (b = 0 ; b < N ; b += 8) {
        for (c = a ; c < a + 4 ; ++c) {
            for (d = b + c - a ; d >= b ; -- d) {
                B[c][d] = A[d][c];
            }
            for (d = b + c - a + 1; d < b + 4 ; ++d) {
                B[c][d] = A[d][c];
            }
        }

        tmp0 = A[b][a + 4];
        tmp1 = A[b][a + 5];
        tmp2 = A[b][a + 6];
        tmp3 = A[b][a + 7];
        tmp4 = A[b + 1][a + 4];
        tmp5 = A[b + 1][a + 5];
        tmp6 = A[b + 1][a + 6];
        tmp7 = A[b + 1][a + 7];
        B[a + 4][b] = tmp0;
        B[a + 4][b + 1] = tmp4;
        B[a + 5][b] = tmp1;
        B[a + 5][b + 1] = tmp5;
        tmp0 = A[b + 2][a + 4];
        tmp4 = A[b + 2][a + 5];
        tmp1 = A[b + 3][a + 4];
        tmp5 = A[b + 3][a + 5];
        B[a + 4][b + 2] = tmp0;
        B[a + 4][b + 3] = tmp1;
        B[a + 5][b + 2] = tmp4;
        B[a + 5][b + 3] = tmp5;
        tmp0 = A[b + 2][a + 6];
        tmp4 = A[b + 2][a + 7];
        tmp1 = A[b + 3][a + 6];
        tmp5 = A[b + 3][a + 7];
        B[a + 6][b] = tmp2;
        B[a + 6][b + 1] = tmp6;
        B[a + 6][b + 2] = tmp0;
        B[a + 6][b + 3] = tmp1;
        B[a + 7][b] = tmp3;
        B[a + 7][b + 1] = tmp7;
        B[a + 7][b + 2] = tmp4;
        B[a + 7][b + 3] = tmp5;

        for (c = a + 4; c < a + 8 ; ++c) {
            for (d = b + c - a ; d >= b + 4 ; --d) {
                B[c][d] = A[d][c];
            }
            for (d = b + c - a + 1 ; d < b + 8 ; ++d) {
                B[c][d] = A[d][c];
            }
        }

        tmp0 = A[b + 6][a];
        tmp1 = A[b + 6][a + 1];
        tmp2 = A[b + 6][a + 2];
        tmp3 = A[b + 6][a + 3];
        tmp4 = A[b + 7][a];
        tmp5 = A[b + 7][a + 1];
        tmp6 = A[b + 7][a + 2];
        tmp7 = A[b + 7][a + 3];
        B[a + 2][b + 6] = tmp2;
        B[a + 2][b + 7] = tmp6;
        B[a + 3][b + 6] = tmp3;
        B[a + 3][b + 7] = tmp7;
        tmp2 = A[b + 4][a + 2];
        tmp3 = A[b + 4][a + 3];
        tmp6 = A[b + 5][a + 2];
        tmp7 = A[b + 5][a + 3];
        B[a + 2][b + 4] = tmp2;
        B[a + 2][b + 5] = tmp6;
        B[a + 3][b + 4] = tmp3;
        B[a + 3][b + 5] = tmp7;
        tmp2 = A[b + 4][a];
        tmp3 = A[b + 4][a + 1];
        tmp6 = A[b + 5][a];
        tmp7 = A[b + 5][a + 1];
        B[a][b + 4] = tmp2;
        B[a][b + 5] = tmp6;
        B[a][b + 6] = tmp0;qizhongqizhong
        B[a][b + 7] = tmp4;
        B[a + 1][b + 4] = tmp3;
        B[a + 1][b + 5] = tmp7;
        B[a + 1][b + 6] = tmp1;
        B[a + 1][b + 7] = tmp5;
    }
}

实际测试的缓存不命中数是1379次。

到了这里我已经优化到了极限，但是依然没有达到满分，最后参考了网上的满分答案。
这里要注意到，实验讲义中说的很清楚，不允许修改矩阵A，但是矩阵B可以任意修改。因此，我们可以通过在矩阵B中暂存转置的结果来充分利用缓存，进一步降低缓存不命中数。思路如下图。

按行加载矩阵A，并且将其存入矩阵B。依次执行4次，直到整个分块的上半部分处理完毕。其中，每行的前4个元素被正确转置，后四个元素被暂存至矩阵B的右上分块。
对于分块的下半部分的第一行，先将矩阵B的右上分块的4个元素载入至临时变量，然后从矩阵A中的左下分块读取第一列并转置进入矩阵右上分块的第一行，然后将读出的4个元素存入矩阵B右下分块的第一行，最后再将矩阵A右下分块第一列转置送入矩阵B右下分块的第一行。
按照2的方式依次处理完下半部分的所有行。

对于一个8×8的分块而言，过程1处理了分块的上半部分，共执行了32次读和32次写。对于对角线上的分块，其读不命中率为1/8，写不命中率为1/4；对于非对角线上的分块，其读不命中率为1/8，写不命中率为1/8。过程2和3处理了分块的下半部分包括将矩阵B的右上子块移动到正确位置，将矩阵A的左下子块转置到B的右上子块以及矩阵A右下子块的转置，共执行了48次读和48次写。对于对角线上的分块，其读不命中率为1/3，写不命中率为1/4；对于非对角线上的分块，其读不命中率为1/12，写不命中率为1/12。

因此对于56个非对角线分块以及8个对角线分块，理论上优化后的总缓存不命中数为(32 × 1/8 + 32 × 1/8) × 56 + (32 × 1/8 + 32 × 1/4) × 8 + (48 × 1/12 + 48 × 1/12) × 56 + (48 × 1/3 + 48 × 1/4) × 8 = 1216次。

该优化方法的代码如下：

for (a = 0 ; a < N ; a += 8) {
    for (b = 0 ; b < M ; b += 8) {
        for (c = b ; c < b + 4 ; ++c) {
            tmp0 = A[c][a];
            tmp1 = A[c][a + 1];
            tmp2 = A[c][a + 2];
            tmp3 = A[c][a + 3];
            tmp4 = A[c][a + 4];
            tmp5 = A[c][a + 5];
            tmp6 = A[c][a + 6];
            tmp7 = A[c][a + 7];
            B[a][c] = tmp0;
            B[a + 1][c] = tmp1;
            B[a + 2][c] = tmp2;
            B[a + 3][c] = tmp3;
            B[a][c + 4] = tmp4;
            B[a + 1][c + 4] = tmp5;
            B[a + 2][c + 4] = tmp6;
            B[a + 3][c + 4] = tmp7;
        }
        for (c = a + 4 ; c < a + 8 ; ++c) {
            tmp0 = B[c - 4][b + 4];
            tmp1 = B[c - 4][b + 5];
            tmp2 = B[c - 4][b + 6];
            tmp3 = B[c - 4][b + 7];

            B[c - 4][b + 4] = A[b + 4][c - 4];
            B[c - 4][b + 5] = A[b + 5][c - 4];
            B[c - 4][b + 6] = A[b + 6][c - 4];
            B[c - 4][b + 7] = A[b + 7][c - 4];

            B[c][b] = tmp0;
            B[c][b + 1] = tmp1;
            B[c][b + 2] = tmp2;
            B[c][b + 3] = tmp3;

            B[c][b + 4] = A[b + 4][c];
            B[c][b + 5] = A[b + 5][c];
            B[c][b + 6] = A[b + 6][c];
            B[c][b + 7] = A[b + 7][c];
        }
    }
}

实际测试的缓存不命中数是1219次。

第三部分 61 × 67矩阵优化

由于61 × 67的矩阵不是方阵，不方便定量分析。同时限制放的比较宽松，满分misses小于2000，misses大于3000零分。因此无需考虑处理对角线，仅尝试换用不同的边长分块即可。16 × 16的分块已可以保证满分。

第三部分优化之后的代码如下：

for (a = 0 ; a < N ; a += 16) {
    for (b = 0 ; b < M ; b += 16) {
        for (c = b ; (c < b + 16) && (c < M) ; ++c) {
            for (d = a ; (d < a + 16) && (d < N) ; ++d) {
                    B[c][d] = A[d][c];
            }
        }
    }
}

实际测试的缓存不命中数是1847次。

实验结果

自动评分脚本给出的输出如下：

Cache Lab summary:
                        Points   Max pts      Misses
Csim correctness          27.0        27
Trans perf 32x32           8.0         8         287
Trans perf 64x64           8.0         8        1219
Trans perf 61x67          10.0        10        1847
          Total points    53.0        53

实验总结

这次实验花费了我很久时间，最后还是参考了网上的解法并且花了很大精力去理解，实验给出的缓存条件非常苛刻，但同时也方便了定量分析。完成后确实大大加深了我对于缓存的理解。

这次实验中我认为比较重要/难的地方：一是对于缓存的理解以及矩阵元素在缓存中的排布问题；二是位于对角线上的分块不仅是内部而且在分块间存在的冲突不命中问题；三是虽然是矩阵，但是也要明确统一控制行列的变量。