C语言、C++矩阵乘法优化

一个正经的开头

本来博主就想写这样一个博客，但是又怕写不好所以一直不敢，直到今天我有个同学说他需要在安卓上跑神经网络，而网上很多开源框架不仅体积大，开发、执行效率也不高，所以想让我给他个基于CPU的矩阵库。

于是乎

矩阵乘法的基本实现方法

首先，矩阵乘法在神经网络的计算中属于占用时间非常多的一个环节，所以把矩阵计算的速度优化好，是能够提升神经网络训练速度的。矩阵的点乘、加法之类的计算很简单，这个没什么好优化的，顶多加个多线程，但是矩阵乘法就不一样了，由于CPU和内存的硬件特点，使得矩阵乘法存在非常大的优化空间。

如果看到这你还感兴趣的话，就跟博主一起慢慢地优化基于CPU的矩阵乘法计算函数吧~

基于数学定义的基本实现方法如下

void MatMul_Lv1(double* C, double *A, double *B, int rowa, int cola, int colb){    int i, j, k;    for (i = 0; i < rowa; i++)     {        for (j = 0; j < colb; j++)        {            for (k = 0; k < cola; k++)             {                C[i*colb + j] += A[i*cola + k] * B[k*colb + j];            }        }    }}

现在为了观测算法的性能，这里统计一下基本操作的次数，即 C[?] += A[?] * B[?] 作为一次基本操作，同时为了简便我们的测试，这里统统使用方阵做测试，虽然有些特殊情况不适用，但是在这只想表达下核心思想。用横坐标表示矩阵维度大小，50代表为50×50的方阵，以此类推。基本实现的性能图如下：

图的纵坐标为基本操作的次数，见上面的定义，单位为10亿次每秒（1G）。可以看到基于数学定义的基本实现方法性能很不好，是很慢的，并且随着矩阵维度的变大，性能越来越低，最后收敛于0.2G，这就是说，计算一个维度为1000的方阵乘法，需要1.5秒左右。现在我们变更一下函数的一个小地方，性能图如下：

可以看到，矩阵乘法函数的计算性能得到了有效地提升，改动后的代码如下：

void MatMul_Lv2(double* C, double *A, double *B, int rowa, int cola, int colb){    int i, j, k;    for (i = 0; i < rowa; i++)    {        for (k = 0; k < cola; k++)        {            for (j = 0; j < colb; j++)            {                C[i*colb + j] += A[i*cola + k] * B[k*colb + j];            }        }    }}

别说你没看到改动在那里……

这里我们把 k 层循环和 j 层循环互换了位置，放心吧，这不会造成错误的结果，但是我们的矩阵乘法计算速度变快了，这是为什么呢！（原谅我被领导带坏了，他发邮件总用感叹号- -）

这里要说到一个叫做 Cache 的东东，也可以叫缓存，现在的CPU通常有三级缓存，即 L1 Cache 、L2、L3…等，L1缓存最小但也最快，通常为32Kb或者64Kb，后面两级以此类推。这个要说起来比较复杂，后面再扯，现在只需要记住一个规则，就是CPU在从内存读取数据时，会同时缓存他周围64字节的数据（常见结构都是64字节，即8个double）。那么现在再来看下方法一和方法二的区别，可以看到方法一的结果矩阵C是按顺序读的，左矩阵A是按顺序读的，但是！右矩阵B不是按顺序读的，是有间隔的。再看看方法二，你会发现所有的矩阵都是按顺序读的，所以方法二更符合缓存的特性，虽然方法一中可以使用寄存器保存中间变量来减少读写次数，但是还是没有方法二快。

可能有的童鞋已经注意到了，不管是方法一还是方法二，都在矩阵维度为1100左右时性能会存在下降，这就是缓存无法保存所有数据造成的性能下降。这里是由于超出L3缓存造成的下降，博主的CPU是AMD Ryzen，三缓是19Mb，如果你用的是Intel的CPU，大约会在维度为700时就下降，Intel的CPU通常在3~8Mb左右。从L1读数据，需要几个时钟周期，而从内存读取，需要几百个时钟周期，要知道一个时钟周期，通常已经可以完成一次我们定义的基本运算了（C+=A*B）

网上现有的大多数C语言矩阵优化也就到此为止了，但是因为这个同学有点牛逼，不卖他点人情我怕以后血亏，所以现在让我们来看看进一步的优化方法~

基于CPU的向量化指令

现在很多神经网络框架都是基于GPU做计算的，无非就是利用了算法可以轻松向量化的优势，虽然CPU在这方面远不如GPU，但是还是有办法利用这个特性的。现在我们将对基本操作（C+=A*B）做向量化。

向量化需要掌握 intrinsic 或者内嵌汇编的知识，但是呢，对于大多数搞科研的孩子们来说，其实不是太关心编程的技巧，那么，有没有办法在不会这些东西的情况下使用向量化指令呢？有！！！

接下来请坐稳扶好

void MatMul_Lv3_C(double* C, double *A, double *B, int len){    double R = A[0];    for (int i = 0; i < len; i++) C[i] += R * B[i];}void MatMul_Lv3(double* C, double *A, double *B, int rowa, int cola, int colb){    int i, k;    for (i = 0; i < rowa; i++)    {        for (k = 0; k < cola; k++)        {            MatMul_Lv3_C(C + i*colb, A + i*cola + k, B + k*colb, colb);        }    }}

如果你使用的是visual studio，那么请打开高性能指令集开关，路径为[项目属性->C/C++->代码生成->启动增强指令集]。接着上述代码便会成功“暗示”编译器，于是乎：

没错，就只是把 j 层单独写成函数这么无脑的操作就有这么大的提升。原因很简单，CPU内部有多媒体寄存器，这类寄存器的特点是很长，特别长那种，嗯！！！通用寄存器只有64位，常见多媒体寄存器有128位和256位的，也就是可以同时计算2个或4个double，所以能够提速。现在高端型号已经有512位的了。

可以看到，向量化指令稍微提升了下矩阵乘法的部分性能，但是超出缓存的部分还是会有很严重的性能损失，那么接下来我们继续优化维度小的部分~

CPU的缓存一致性

CPU为了保证任务不出错，当有一个缓存的数据被改写时，总线就会广播，其他位置对应到该数据的缓存位置就会失效（数据进缓存是有特定位置的，不是随便放的，这次进缓存在这里，下次来缓存，座位不变哦）。也就是说我一个六核心CPU其中一个核心进行一次基础计算（C+=A*B）时，自身+其他五个核心的所有L1、L2、L3缓存都- - - -简直浪费生命！

为了减少向内存中写回数据，我们对方法三进行改进，代码如下：

void MatMul_Lv4_C4(double* C, double *A, double *B, int len){    double a0 = A[0], a1 = A[1], a2 = A[2], a3 = A[3];    double *b0 = B, *b1 = B + len, *b2 = B + len * 2, *b3 = B + len * 3;    for (int i = 0; i < len; i++) {        C[i] += a0*b0[i] + a1*b1[i] + a2*b2[i] + a3*b3[i];    }}void MatMul_Lv4_C1(double* C, double *A, double *B, int len){    double R = A[0];    for (int i = 0; i < len; i++) C[i] += R * B[i];}void MatMul_Lv4(double* C, double *A, double *B, int rowa, int cola, int colb){    int i, k;    for (i = 0; i < rowa; i++)    {        for (k = 0; k < cola - cola % 4; k += 4) {            MatMul_Lv4_C4(C + i*colb, A + i*cola + k, B + k*colb, colb);        }        for (; k < cola; k++) {            MatMul_Lv4_C1(C + i*colb, A + i*cola + k, B + k*colb, colb);        }    }}

这里我们的点乘一次性计算四个元素，即四次写回操作合成了1次，其他三次写回是在寄存器上操作的，性能图如下：

呐，现在相对于基于数学定义的基本实现，是不是快了很多。现在，虽然小维度还有油水可挖，但是不多了，可以愉快地优化大维度的部分了~

缓存的合理复用

看了前面的分析，我们已经知道了，当矩阵维度过大时，缓存容量不足会造成矩阵乘法的性能下降，那么，我们只需要让缓存的容量一直够用不就行了嘛~代码如下：

inline void MatMul_Lv5_C8(double* C, double* A, double* B, int colb, int len){    double a0 = A[0], a1 = A[1], a2 = A[2], a3 = A[3];    double a4 = A[4], a5 = A[5], a6 = A[6], a7 = A[7];    double *b0 = B, *b1 = B + colb, *b2 = B + colb * 2, *b3 = B + colb * 3;    double *b4 = B + colb * 4, *b5 = B + colb * 5, *b6 = B + colb * 6, *b7 = B + colb * 7;    register double R = 0;    for (int j = 0; j < len; j++) {        R = 0;        R += a0*b0[j] + a1*b1[j] + a2*b2[j] + a3*b3[j];        R += a4*b4[j] + a5*b5[j] + a6*b6[j] + a7*b7[j];        C[j] += R;    }}inline void MatMul_Lv5_C1(double* C, double* A, double* B, int len){    double R = A[0];    for (int j = 0; j < len; j++) C[j] += R * B[j];}inline void MatMul_Lv5(double* C, double* A, double* B, int rowa, int cola, int colb){    int i, j, k;    int r = 96;    for (j = 0; j < colb - colb%r; j += r) {        for (k = 0; k < cola - cola % 8; k += 8) {            for (i = 0; i < rowa; i++) {                MatMul_Lv5_C8(C + i*colb + j, A + i*cola + k, B + k*colb + j, colb, r);            }        }        if (cola % 8) {            for (k = cola - cola % 8; k < cola; k++) {                for (i = 0; i < rowa; i++) {                    MatMul_Lv5_C1(C + i*colb + j, A + i*cola + k, B + k*colb + j, r);                }            }        }    }    if (colb%r) {        for (k = 0; k < cola - cola % 8; k += 8) {            for (i = 0; i < rowa; i++) {                MatMul_Lv5_C8(C + i*colb + j, A + i*cola + k, B + k*colb + j, colb, colb%r);            }        }        if (cola % 8) {            for (k = cola - cola % 8; k < cola; k++) {                for (i = 0; i < rowa; i++) {                    MatMul_Lv5_C1(C + i*colb + j, A + i*cola + k, B + k*colb + j, colb%r);                }            }        }    }}

这里我们将方法四中连续使用A中四个元素改成了八个元素，因为64字节是八个double，然后我们调整了循环层之间的关系，我们固定矩阵B中的一个小块（8×96）的一个块，然后从矩阵A中不断地抽取需要跟这个块做计算的元素（一次八个double），当B中的小块不再需要计算时我们就换一个小块，这就是矩阵分块的思想。性能图如下：

看图~

简单地使用多线程进行加速

要单机使用多线程无非两种方法，要么自己用API开线程，要么用OpenMp这种无脑指令集，基于现有的状态，我们将方法5的函数用多线程执行。如果你想用OpenMp，那么，不好改！！！所以这里讲讲怎么自己开简单的API，如果你不是特别有兴趣，建议使用C++里的多线程库，如果真的很有兴趣，那么去研究CreateThread或者Linux下的API吧，这里为了方便，用C++的多线程库。代码如下：

inline void MatMul_Lv6(double* C, double* A, double* B, int rowa, int cola, int colb){    int n = 12, m = rowa / n, d = rowa%n, i;    std::thread *T = (std::thread*)calloc(n, sizeof(std::thread));    // 我不推荐使用new这种关键字，当然仅仅是我的建议    if (T == NULL) {        printf("\n=> in fun<MatMul_Lv5>, out of memory\n");        getchar(); exit(0);    }    for (i = 0; i < n - 1; i++) {        T[i] = std::thread(MatMul_Lv5, C + i*colb*m, A + i*cola*m, B, m, cola, colb);    }    T[i] = std::thread(MatMul_Lv5, C + i*colb*m, A + i*cola*m, B, m + d, cola, colb);    for (i = 0; i < n; i++) {        T[i].join();    }    free(T);}

接下来，面对疾风吧：

可以看到，除了小矩阵因为线程的开销而导致性能低下之外，随着维度越来越大，性能越来越好！这就是多线程的好处，多线程之间的切换可以一定程度抵消CPU从内存读取数据所产生的延迟。博主CPU是6C12T，所以开了12个线程，线程的数量要根据具体情况确定，不是越大越好。

有兴趣的童鞋可以再想想有哪些地方可以优化，是可以超越matlab的矩阵计算速度的哦~当然超不超越不重要，编程嘛，自己开心就好。用C\C++优化矩阵乘法的相关内容就到这了，感谢~

稍微总结一下

第一，也是最重要的一点，以上所有程序博主并没有检验，不保证正确性，仅仅提供思路作为参考。

第二，如果你是做科学计算的，在现在这个环境下，能用GPU就用GPU吧- -我已经头破血流了

第三，上面的优化方法并不是最优的，仅提供思路

第四，博主是业余爱好并非科班出身，如果有什么概念说错了，欢迎批评指正。

最后，感谢~~~