今天以矩阵乘法为例分享一点软件优化的技巧。

    矩阵乘法在很多领域都有应用，最著名的要数这几年流行的深度学习和卷积神经网络了。对于常见的CNN网络，比如MobileNet，计算量最大的莫过于卷积层，而卷积层中，最主要的运算就是矩阵乘法。

    为了便于理解后面的代码，先借一张图来复习一下矩阵点积（我们要讨论的乘法形式）的计算方式。

    简单说，对于矩阵A和矩阵B的点积C，计算C中0号元素的方法就是把A中0行与B中0列的各个元素相乘，再把各个乘机求和，即：

    1*1 + 4*2 + 6*9 + 10*3 = 93

    顺便说一下，就为了加速上面这样的乘法和加法交替进行的计算，今天的大多数处理器中都支持所谓的乘加指令，常常称为FMA（Fused Multuply-Add）。

    理解了基本规则后，只要学过一点C语言的人，都可以把上述过程实现出来，下面是最直观的一种实现（代码来自VTune的示例程序）。

void multiply1(int msize, int tidx, int numt, TYPE a[][NUM], TYPE b[][NUM], TYPE c[][NUM], TYPE t[][NUM])

{     int i,j,k;    // Naive implementation     for(i=tidx; i<msize; i=i+numt) {         for(j=0; j<msize; j++) {             for(k=0; k<msize; k++) {                     c[i][j] = c[i][j] + a[i][k] * b[k][j];            }        }    } }

    在上面代码中，有一行注释，Naive implementation，意思是这是最简单的实现，或者说，很单纯的实现，使用网络流行语，说的难听一点，就是比较白痴的实现。我们不妨就将其称为小白实现吧。

     编译运行上面的小白实现，会发现慢的不行。对于2048*2048的A、B矩阵，加上多线程支持，把Kaybe Lake i7的8个逻辑CPU全用上，也要60s左右。

    在强大的i7上面，1分钟都超过了，太慢了。

    为什么这么慢呢？上VTune一看就知道了。

    在VTune颜色体系中，红色常常代表着“出问题”的地方，看右面的微管线（uPipe）示意图，一大片红色，上面写的Memory Bound代表内存约束。形象一点说，强大的i7工厂因为内存方面的制约，只发挥了1/4左右的产能（绿色部分）。

    为什么这样呢？看左侧的详细数据，可以看到大量的LLC（Last Level Cache） Miss。

    不是说英特尔的Cache预测算法很牛么？怎么还有如此高的LLC Miss呢？

    简单说，问题在于访问内存的方式上，注意观察最内层的循环：

            for(k=0; k<msize; k++) {

                    c[i][j] = c[i][j] + a[i][k] * b[k][j];    

    这个内层循环的执行次数是最多的，在这个循环中，c[i][j]有读有写，a[i][k]和b[k][j]的访问方式都是读。

    思考相邻的两次循环，c[i][j]保持不变，a[i][k]与a[i][k+1]是相邻的，大多数时候应该不会cache miss，而b[k][j]与b[k+1][j] 则不然了，二者相隔一行，距离2047个元素，太远了，显然违反了cache的局部性（locality）原则。   

    长话短说，对于这样的经典问题，有个经典的解决方法，名叫“循环交换”，简单讲就是调换内外层循环的顺序，在保持代码逻辑正确的前提下，让内存访问尽可能满足局部性原则。

    对于上面的小白实现，可以将其调整为如下形式：

void multiply2(int msize, int tidx, int numt, TYPE a[][NUM], TYPE b[][NUM], TYPE c[][NUM], TYPE t[][NUM])

{     int i,j,k;// Step 2: Loop interchange    for(i=tidx; i<msize; i=i+numt) {         for(k=0; k<msize; k++) {             for(j=0; j<msize; j++) {                 c[i][j] = c[i][j] + a[i][k] * b[k][j];            }        }    } }

    仔细观察上面的multiply2函数，特别是最内层循环，再思考相邻两次循环：c[i][j] 和 a[i][k]都不变，第二次访问时一定会命中cache，而b[k][j]和b[k][j+1]是地址相邻的，命中的概率也非常高。

    把修改后的代码编译运行，结果喜人。

    睁大眼睛啊，在相同的软硬件环境下，只有1.02秒了，只是原来的60/1。或者说快了60倍。

    再用VTune看一下，管线图上，红色部分小多了，内存限制从原来的70%降低到12%了。同时，绿色部分大大增加，8端口的x86核心有6成以上时间都使用了3个以上（见左侧小字）。

    在软件调优中，局部性是个很重要的概念。又可以细分为时间局部性和空间局部性，借用一下老雷自己的一张讲义吧。

    简单理解，空间局部性就是让要用的代码和数据尽可能紧凑，尽可能小，这样就可以更容易fit到cache中，时间局部性，可以保持要访问的内容保持在cache中，形象一点说，就是保持cache的温度。二者都有利于提高cache的命中率，降低cache miss。

    在真是的软件环境中，实现局部性会面临各方面的调整。或者说，很多因素都会破坏局部性，比如，硬件终端，软件异常，C++中的多态特征，链表等复杂结构体，操作系统的线程切换，第三方库的影响，NUMA等等。

    考虑常见的干扰因素，把上面的实现再经过一番异化后，可以做到红色部分全无，一片绿色，强大的x86核心畅行无阻。

    要做到这个效果，需要绑定核心，并且增加prefetch指令，向cache单元预定将要访问的内存。能用上SIMD指令也会提高使用率，增加数学计算（相对于内存访问）的速度。

    时间关系，以后再详细介绍后面几个优化步骤，条件是能有足够多的读者。希望大家读了之后，如果喜欢的话， 那么请转发和分享，满5000阅读量后，立刻更新下一篇。

    或者，如果觉得读文章不够干脆，那么可以参加2019年1月份在魔都上海举行的，2天1晚，与老雷面对面，深度体验软件优化的无限挑战和无穷魅力。

***********************************************************

正心诚意，格物致知，以人文情怀审视软件，以软件技术改变人生。

欢迎关注格友公众号