杂记1：缓存，局部性，循环

——“由剑剑同学的两段代码引发的讨论”
跟Jay同学讨论循环不变量的优化问题，在email里比较难说清楚，还是发到blog上来方便贴代码。

剑剑同学 写道

某型CPU的一级数据缓存大小为16K字节，cache块大小为64字节；二级缓存大小为256K字节，cache块大小为4K字节，采用二路组相联。经测试，下面两段代码运行时效率差别很大，请分析哪段代码更好，以及可能的原因。
为了进一步提高效率，你还可以采取什么办法？
A段代码：

int matrix[1023][15];
const char *str = "this is a str";
int i, j, tmp, sum = 0;

tmp = strlen(str);
for (i = 0; i < 1023; i++) {
    for (j = 0; j < 15; j++) {
        sum += matrix[i][j] + tmp;
    }
}

B段代码

int matrix[1025][17];
const char *str = "this is a str";
int i, j, sum = 0;

for (i = 0; i < 17; i++) {
    for (j = 0; j < 1025; j++) {
        sum += matrix[j][i] + strlen(str);
    }
}

这条问题和大家分享下，希望大家一起讨论。呵呵。

这个问题的主要关注点很明显是关于存储器层次（memory hierarchy）与缓存（caching）的。先看看相关背景。

存储技术在几个不同层次上发展，其中存储密度高、单价便宜的存储器速度比较慢，速度快的存储器的存储密度则相对较低且价格昂贵。为了在性能与价格间找到好的平衡点，现代计算机系统大量采用了缓存机制，使用小容量的高速存储器为大容量的低速存储器提供缓存。
最快的存储器是CPU里的各种寄存器，其次是在CPU芯片内的L1缓存，再次是在CPU芯片内或者离CPU很近的L2缓存，然后可能还有L3缓存，接着到主内存，后面就是各种外部存储设备如磁盘之类，最后还有诸如网络存储器之类的更慢的存储器。
L1缓存可能会成对出现，一个用于指令，另一个用于数据。L2缓存和后面的缓存则设计得更通用些。由于主内存比磁盘快很多但相对来说价格昂贵许多，而同时运行多个程序所需要的存储空间通常不能直接被主内存满足，所以现代操作系统一般还有虚拟内存，使用磁盘作为主内存的扩充。虚拟内存也可以反过来看作“所有虚拟内存都是在磁盘上的，将其中活跃的一些放在物理内存里是一种优化”（Eric Lippert如是说）。

如果要访问的数据位于存储器层次的较低层，则数据是逐层传递到CPU的。例如，程序要访问某个地址的数据，在L1缓存里没有发现（称为L1缓存不命中，L1 cache miss），则跑到L2缓存找；找到的话，会先把这一数据及邻近的一块数据复制到L1缓存里，然后再从L1缓存把需要的数据传给CPU。

为了充分利用缓存机制，程序应该有良好的局部性（locality）。局部性指的是程序行为的一种规律性：在程序运行中的短时间内，程序访问数据位置的集合限于局部范围。局部性有两种基本形式：时间局部性（temporal locality）与空间局部性（spatial locality）。由于指令也可以看作数据的一种特殊形式，因而局部性对指令来说也有效。
时间局部性指的是反复访问同一个位置的数据：如果程序在某时刻访问了存储器的某个地址，则程序很可能会在短时间内再次访问同一地址。例如，在执行一个循环，则循环的代码就有好的时间局部性；又例如在循环里访问同一个变量，则对该变量的访问也有好的时间局部性。
空间局部性指的是反复访问相邻的数据：如果程序在某时刻访问了存储器的某个地址，则程序很可能会在短时间内访问该地址附近的存储器空间。例如按顺序执行的指令就有良好的空间局部性；又例如按存储顺序挨个遍历数组，也有良好的空间局部性。

还有很多很重要的背景信息，这里就不详细写了。我主要是读《Computer Organization and Architechture: Designing for Performance, 5th Edition》和《Computer Systems: A Programmer's Perspective》学习的。大一的时候也好好上了计算机组成与结构的课，用的课本就是前一本书，还能记得一些。

那么回到开头的题目。两段代码有一些特征是相同的，包括：
(1) 它们都使用了一个int矩阵，而且行的宽度比列的长度要短。
(2) 它们都含有一个char指针，指向的是一个字符串字面量。这意味着对该字符串调用strlen()总是会得到同一个值，而且该值在编译时可计算。
(3) 它们都遍历了整个矩阵，并且对矩阵中每个元素的值求和。

两段代码主要的差异是：
(1) 遍历顺序不同。A按行遍历，B按列遍历。
(2) 内外循环的分布不同。使用两层的嵌套循环来遍历这个数组，则：按行遍历的话，外层循环次数等于行数，内层循环次数等于列数；按列遍历则正好相反。A的外层循环比内层循环次数多很多；而B的内层循环次数比外层多。
(3) 循环中是否重复求值。A在遍历矩阵前预先对字符串调用了strlen()，将结果保存在一个临时变量里；遍历矩阵时访问临时变量来获取该值；B在遍历矩阵时每轮都调用strlen()。
(4) A与B的矩阵大小不同，A较小而B较大。

========================================================================

《Computer Systems: A Programmer's Perspective》的6.2.1小节，Locality of References to Program Data介绍了程序数据的局部性。其中提到一个概念：访问连续存储空间中每隔k个的元素，称为stride-k reference pattern。连续访问相邻的元素就是stride-1访问模式，是程序中空间局部性的重要来源。一般来说，随着stride的增大，空间局部性也随之降低。

C的二维数组在内存中是按行优先的顺序储存的。许多其它编程语言也是如此，但并非全部；FORTRAN就是一种典型的例外，采用列优先的存储顺序。在C中，按行遍历二维数组，遍历顺序就与存储顺序一致，因而是stride-1访问模式。如果按列访问一个int matrix[M][N]，则是stride-(N*sizeof(int))访问模式。
由此可知，差异(1)使得A段代码比B段代码有更好的空间局部性，因而应该能更好的利用缓存层次。

结合两段代码中矩阵的大小来看看缓存不命中的状况。题目没有提到“某型CPU”上int的长度是多少，也没有提到内存的寻址空间有多大。这里把两者都假设是32位的来分析。

32位的int意味着sizeof(int)等于4B。则A段代码的matrix大小为sizeof(int)*15*1023 == 61380B，小于128KB；matrix的每行大小为sizeof(int)*15 == 60B，小于64B。B段代码的matrix大小为sizeof(int)*17*1025 == 69700B，也小于128KB；每行大小为sizeof(int)*17 == 68B，大于64B。

题中L1缓存大小是16KB，每条cache line大小是64B，也就是说一共有256条cache line。没有说明L1与L2缓存的映射方式，假设是直接映射。
L2缓存大小是256KB，每条cache line大小4KB，也就是说一共有64条cache line。因为L2缓存是二路组相联，所以这些cache line被分为每两条cache line一组，也就是分为32组。同样因为是二路组相联，所以主内存中地址连续的数据在把L2缓存的一半填满之后，要继续填就要开始出现冲突了。幸好L2缓存有256KB，两段代码中都能顺利装下各自的matrix。假设两层缓存都采用LRU（least recently used）算法来替换缓存内容。

A段代码中，matrix的一行可以完整的放在一条L1 cache line里，一条L2 cache line可以装下68行多一些。观察其遍历的方式。假设两层缓存刚开始都是“冷的”，访问matrix[0][0]时它尚未被加载到L2缓存，并且假设matrix[0][0]被映射到一条L2 cache line的起始位置（意味着matrix在4KB对齐的地址上）。
这样，在第一轮内层循环时访问matrix[0][0]，会发生一次L1缓存不命中和一次L2缓存不命中，需要从主内存读4KB到L2缓存，再将其中64字节读到L1缓存。第二轮内层循环时，访问matrix[0][1]在L1缓存命中。第三轮也是L1缓存命中。直到读到matrix[1][1]的时候才会再发生一次L1缓存不命中，此时L2缓存命中，又从L2缓存读出64字节复制到L1缓存。重复这个过程，直到遍历了68行多一些的时候，又会发生一次L2缓存不命中，需要从主内存读数据。tmp与sum变量有良好的时间局部性，应该能一直在寄存器或者L1缓存中。以此类推，可以算出：每轮外层循环都执行15轮内层循环，遍历了matrix的一整行；每遍历16行会发生大约15次L1缓存不命中（如果matrix[0][0]不是被映射到cache line的开头的话，会发生16次）；每遍历1023行会发生大约15次L2缓存不命中。加起来，A段代码在循环中大概会遇到15次L2缓存不命中，960次L1缓存不命中。

B段代码中，matrix的一行无法完整的放在一条L1 cache line里，一条L2 cache line可以装下60行多一些。遍历的元素相隔一行。同样假设两层缓存刚开始都是“冷的”，则遍历过程中刚开始每访问matrix的一个元素都会发生1次L1缓存不命中，每访问60行多一些就会发生1次L2缓存不命中。等遍历完了matrix的第一列之后，经过了1轮外层循环（1025轮内层循环）；此时两层缓存都已经热起来，整个matrix都被缓存到L2中；根据LRU算法，matrix[0][0]已经不在L1缓存中。照此观察，后续的遍历过程中都不会再出现L2缓存不命中，但每访问一个元素仍然会发生一次L1缓存不命中。加起来，B段代码在循环中大概会遇到18次L2缓存不命中，17425次L1缓存不命中。

遍历顺序与矩阵大小结合起来，使A段代码发生L1缓存不命中的次数远小于B段代码的，而两者的L2缓存不命中次数差不多。因此，从缓存的角度看，A段代码会比B段代码执行得更有效率。

========================================================================

然后再看看在循环中调用strlen()的问题。从源码表面上看，B段代码的每轮内层循环中都要调用一次strlen()，其中要遍历一次str字符串。strlen()本身的时间开销是O(L)的（L为字符串长度），放在M×N的嵌套循环里调用，会带来O（L×M×N）的时间开销，相当可观。
但前面也分析过，题中两段代码都是对字符串字面量调用strlen()，是编译时可以计算的量，所以会被编译器优化为常量。事实上VC和GCC都会将这种情况下的strlen()的调用优化为常量。所以这题里在循环中调用strlen()并不会带来额外的开销——因为编译出来的代码里就不会在循环里调用strlen()了。
即使不是对字符串字面量调用strlen()，如果str在循环中没有改变，那么strlen(str)的结果也应该是循环不变量，理论上B段代码可以被编译器自动优化为A段代码的形式，将strlen()的调用外提。不过在许多例子里，VC与GCC似乎都没能成功的进行这种优化。以后会找个实际例子来看看。

========================================================================

这个题目里的代码还不仅涉及缓存层次的问题，还涉及到指令执行的问题。现代CPU一般都支持指令流水线（instruction pipelining）和预测性执行(speculative execution)。通过将一条指令拆分为多个可以并行执行的阶段，CPU的一个执行核心可以在一个时钟周期内处理多条指令；通过预先将后面的指令读进CPU执行，CPU可以预测将来的执行结果。为了能尽可能多的预测执行结果，CPU会对分支指令也做预测，猜测其会进行跳转（branch taken）还是不跳转（branch not-taken）。实际执行到跳转指令的时候，并不是“发现需要跳转到某地址”，而是“印证先前就发现的跳转的猜测”。如果猜中了，则执行结果就会从一个缓存写到寄存器中；如果猜错了，就只能刷掉之前猜测的执行结果，重新读取指令，重新开始流水线的执行，从而带来相当的开销。对分支的预测称为branch prediction，猜错的情况称为branch misprediction。分支预测有许多算法，多数都会考虑某条分支指令上一次或多次的跳转情况。

为了让循环能正常结束，循环一般都有循环条件。这样就至少有一个条件跳转。可以想象，重复多次的循环，控制其结束的条件分支，除了最后一次都应该是向同一个目标跳转的。这样，每个循环至少会导致一次分支预测错误。计算循环条件本身也有一定开销，与分支预测错误一起，都是循环的固有开销。
在嵌套循环中，无论是内层循环还是外层循环，都是循环，固有开销是避免不了的。把重复次数多的循环放在内层与外层会导致总的循环次数的不同。开头的题目中，如果A与B的matrix都统一为1024*16的大小，则A总共要执行1024+1024*16 = 17408次循环，而B总共要执行16 + 16*1024 = 16400次循环。显然，把重复次数多的循环放在内层比放在外层需要执行的循环次数少，相应需要付出的循环固有开销也小。

题目问到要进一步提高效率应该采用什么办法。从前面的分析看，A在缓存方面有优势但在指令执行方面有劣势。如果要改进，可以把A中的matrix转置为int matrix[15][1023]，使行的宽度比列的长度长。这样在按行遍历时重复次数较多的就从外层循环移到了内层，扭转了A段代码在指令执行方面的劣势。

========================================================================

前面的分析都属于“理想分析”，现实中我们写的程序在实际机器上到底是怎么执行的，那简直就是magic。虽然Eric说别把东西想象成magic，但这里我没办法……

例如说，我们不知道题目中的程序一共开了多少个线程。既然题目没说“某型CPU”是多核的，假设它是单核的，那么多个线程都要共享同一个L1和L2缓存，留给A段代码用的缓存到底有多少呢？就算不考虑线程的多少，操作系统也有些核心数据会尽量一直待在高速缓存里，留给应用程序的缓存有多少呢？

既然我们知道要遍历连续的数据，那与其让它逐渐进入缓存，还不如先一口气都放进缓存，后面实际访问数据的时候就不会遇到缓存不命中。这叫做预取（prefetch）。在x86上有专门的指令prefetch-*来满足预取的需求，如非时间性的prefetchnta与时间性的prefetcht0、prefetcht1等等。编译器有没有为代码生成预取指令？使用预取之后缓存不命中的状况能减少多少？不针对具体情况都没办法回答。毕竟有些CPU实现的时候干脆就忽略指令中的预取，又或者编译器生成了很糟糕的预取指令反而降低了程序性能，这些极端的可能性都存在。

另外一个要考虑的因素是，应用程序构建在操作系统之上，而操作系统一般有采用分页的虚拟内存。像32位Windows的页大小就是4KB。matrix有60KB左右，无法完整放在一页里。页在映射到物理内存的时候，并不保证在matrix跨越不同页仍然保持在物理内存中地址的连续性。所以matrix是否能理想的缓存到L2缓存而不发生冲突，其实不好说。

CPU支持的指令集与其实际执行的方式也不完全一致。像x86这样的指令集早就成为“遗留接口”了，实际硬件用类似RISC的方式去实现了CISC的x86指令集，通过指令级并行执行（instruction-level parallelism，ILP）来提高CPU的吞吐量。
x86一个很讨厌的地方就是它可用的通用寄存器（general purpose register）的数量太少了，32位GPR只有8个。那么少的寄存器，指令是怎么并行起来的呢？其实那8个GPR也是假象，CPU可以通过寄存器重命名（register renaming）的方式让一些指令可以直接把计算结果传给下一条指令而不需要实际经过寄存器。预测性执行的结果也不是直接写到寄存器，而是等分支预测被确认正确后才写进去。这样就能够预测性执行多条指令而不破坏“当前”的CPU状态。

It's magic...应用程序员一般也不会需要关心这种magic般的细节。在合适的抽象层次上选用合适的算法，用清晰的方式把代码组织起来，远比关心这种细节要重要得多。不过如果要写编译器的话，这些细节就是恶魔了。Devil is in the details……

========================================================================

Jay同学对编译器处理循环和strlen()的方式感兴趣。下一篇简单分析一下strlen()的特性。

评论

14 楼 RednaxelaFX 2009-06-22  

幸存者 写道

但是具体的循环次数我不知道是怎么算出来的。

顺带：循环轮数就是把外层循环执行的次数，加上它与每轮外层循环中内层循环执行的次数的乘积。循环的固有overhead包括每轮都有的条件检查，和每个循环的分支预测错误。或许还有些其它的？
一个循环如果要终止，总得通过某种方法；要么是条件分支，要么是中断或异常之类。中断之类的不在预测的考虑范围内，因为无从预测……

13 楼 RednaxelaFX 2009-06-22  

幸存者 写道

另外关于预测失败的那部分，当内层循环较大时预测失败较少应该是没什么问题了，但是具体的循环次数我不知道是怎么算出来的。按我的理解，每次循环首尾应该都是会预测失败一次的，即使只有循环末预测失败一次，但是外层循环也应该有预测失败，所以我无论怎么算似乎都和文中的结果不同。而且，即使预测失败其开销也未必就是多执行一次循环吧？

杂记里我没写具体会有多少次预测错误，只是很含糊的写了每个循环至少会在结束的时候预测错误。注意到我在杂记里没只写了循环的轮数，但没有把循环轮数与预测失败次数联系在一起，而是与循环个数联系在一起。
像这样：

for (i = 0; i < 10; i++) {
    // ...
}

我是说这种结构会至少导致一次条件分支的预测失败。

for (i = 0; i < 5; i++) {
    for (j = 0; j < 7; j++) {
        // ...
    }
}

这种结构则至少导致5+1=6次条件分支失败，因为执行了5个内层循环和一个外层循环（而不是按5×7=35轮内层循环和5轮外层循环来算）。

这个“至少”是怎么来的呢？其实不应该说“至少”的，因为没有提供预测的算法，没有足够根据。
如果是最简单的、不考虑先前的跳转记录的预测算法，可以选用“只向回跳”的预测算法，也就是如果遇到分支，则总是预测会跳往低地址的目标，因为那样符合循环的特征，而循环比较容易成为“热点”。这样的话循环结束的条件分支就肯定会预测错误。
但如果用其它算法呢？例如说简单的累加taken与not-taken，选取较多的一侧，如果外加初始优先向回跳的话，每个循环最后的条件跳转也会有一次预测错误。
但如果预测算法能发现循环模式的话呢？那说不定循环末尾就不会出现预测错误了。“至少”在这里不是严谨的说法，不过在不提供实际算法的时候，我也找不到什么好词来描述这个装抗……

你是如何计算预测错误的次数的呢？

12 楼 幸存者 2009-06-22  

发现之前写错一个字predication（其实应该是prediction），刚才查wikipedia发现还真有一个branch predication，和branch prediction是不同的两个概念。

另外关于预测失败的那部分，当内层循环较大时预测失败较少应该是没什么问题了，但是具体的循环次数我不知道是怎么算出来的。按我的理解，每次循环首尾应该都是会预测失败一次的，即使只有循环末预测失败一次，但是外层循环也应该有预测失败，所以我无论怎么算似乎都和文中的结果不同。而且，即使预测失败其开销也未必就是多执行一次循环吧？

11 楼 RednaxelaFX 2009-06-22  

幸存者 写道

话说CSAPP真是本好书，关于locality和branch predication书中都有详尽的解释。

想到一个事，如果有谁能用汇编写出性能比现代编译器编译出来的C程序的性能更高的代码，那就不止是magic了，那是god！

CSAPP里分支预测的介绍是不错，读起来感觉比较愉快。不过毕竟是比较原理方面的介绍，而且也已经出版好几年了。之前读到的一篇很有趣的分支预测算法CSAPP就没提到，是关于模拟神经分支预测的。诶，硬件我还是不熟悉，现代硬件怎么看都是magic……

有一种说法是，人总是能写出比编译器更高效的代码。要做到这点很简单：先用高级语言写出程序逻辑，让编译器编译，然后手工去调整汇编。编译器只能做到optimizing，做不到完全optimized，所以总能留下机会让人来做修改来提高性能。不过在指令特别复杂、寄存器特别多的时候怎么做才比较高效，光靠人脑确实是越来越难跟上了。

10 楼 幸存者 2009-06-22  

话说CSAPP真是本好书，关于locality和branch predication书中都有详尽的解释。

想到一个事，如果有谁能用汇编写出性能比现代编译器编译出来的C程序的性能更高的代码，那就不止是magic了，那是god！

9 楼 RednaxelaFX 2009-06-22  

night_stalker 写道

但是这里把 sum 返回，最后让 main() return 了 …… 应该不会被咔嚓吧？

今晚太困了……明天起来check一下 =v=

8 楼 night_stalker 2009-06-22  

但是这里把 sum 返回，最后让 main() return 了 …… 应该不会被咔嚓吧？

7 楼 RednaxelaFX 2009-06-22  

night_stalker 写道

原来是犯了傻x的错误…… 测了同一个…… 修改后的代码：
//...

擅自增加了个 k，不过不妨碍感性认识 ……

呵呵，这样不行的哦，关键的matrix数组会被优化掉的……要挂住它必须要在循环给matrix赋值之后还要使用它才行，例如说把matrix返回出去，把sum用printf()输出之类。杂记3里就有讲到这个，刚才貌似手滑把那篇的草稿发出来了 orz

6 楼 night_stalker 2009-06-22  

原来是犯了傻x的错误…… 测了同一个…… 修改后的代码：

#include <cstring>
#include <cstdio>

int A(char* str) {
	int matrix[1023][15];
	int i, j, tmp, sum = 0;

	tmp = strlen(str);
	for (int k = 0; k < 1000; k++) {
		for (i = 0; i < 1023; i++) {
			for (j = 0; j < 15; j++) {
				sum += matrix[i][j] + tmp;
			}
		}
	}
	return sum;
}

int B(char* str) {
	int matrix[1025][17];
	int i, j, sum = 0;

	for (int k = 0; k < 1000; k++) {
		for (i = 0; i < 17; i++) {
			for (j = 0; j < 1025; j++) {
				sum += matrix[j][i] + strlen(str);
			}
		}
	}
	return sum;
}

int main(int argc, char** argv) {
	if (argc == 2) {
		printf("A: ");
		return A(argv[0]);
	} else {
		printf("B: ");
		return B(argv[0]);
	}
}

require 'benchmark'
puts Benchmark::CAPTION
s = './a thisisastr'
b = s + ' b'
[s, b, s, b].each do |s_b|
  puts Benchmark.measure { system s_b }
end

      user     system      total        real
A:   0.000000   0.016000   0.077000 (  0.060000)
B:   0.000000   0.000000   0.748000 (  0.729000)
A:   0.000000   0.000000   0.061000 (  0.045000)
B:   0.000000   0.000000   0.749000 (  0.731000)

擅自增加了个 k，不过不妨碍感性认识 ……

5 楼 night_stalker 2009-06-21  

一开始把它们放到两个 cpp 去做的…… 差别巨大。

但是后来写成两个函数放到同一个文件后，动作就很一致了 -_-。不过 A 还是快一小丁点。

我把 argc 传给这两个函数做 sum 的初始值，汇编看不懂…… 写了个去循环版对比，似乎没优化掉循环。

4 楼 RednaxelaFX 2009-06-21  

night_stalker 写道

唔， 纯粹的 gcc a.cpp …… 可能会被优化掉，等我加个 volatile 试试

gcc a.cpp呃（摸胡子……）。至少开-O2嘛
别加volatile，那个可能会干扰代码顺序相关的优化。最简单的办法就是在这块代码前先time一下，之后再time一下，然后循环printf()把sum和matrix都挂住，顺带输出个time来看看……|||

本来这篇还有好几部分的。越写越长，眼看就要坑了。还是决定拆分，先发第一段……

3 楼 night_stalker 2009-06-21  

唔， 纯粹的 gcc a.cpp …… 可能会被优化掉，等我加个 volatile 试试

2 楼 RednaxelaFX 2009-06-21  

嘿嘿，NS老兄怎么测的？贴一下代码来看看～用什么编译器和什么参数？
之前同学刚问我这题的时候，我也是直觉觉得A快，然后就把A和B都包成函数，用-O2编译，准备运行之前觉得有点诡异，就先-S看了下生成的代码，发现A和B都变成了空函数……OTL 所以想看看你是用什么方法把代码留住的。用printf()了么？

1 楼 night_stalker 2009-06-21  

直觉觉得是 A 快 ……

本机测了下， A B A B 100 次：
   total        real
0.046000 (  0.031000)
0.155000 (  0.131000)
0.061000 (  0.030000)
0.155000 (  0.130000)