【代码札记】基于Cache的矩阵乘积算法性能改进

标题挺唬人，其实没啥。可就这么一点小东西，让我三年来第一次体会到大学的价值，与正经程序员的专业所在。

背景展开目录

今天看了看体系结构课程的实验，之前这门课上的都是一些理论课，内容就跟计算机组成原理差不多，又是选修课，老师凑活着讲，我凑活着听。

其中第二个实验是「基于 Cache 的矩阵乘积算法性能改进」。看了看实验要求，就是写一个矩阵乘法的代码，一种是行列相乘累加，另一种则是要先计算转置矩阵，然后利用转置矩阵计算乘法，本质上还是行列相乘累加，而且还额外多了一个求转置。一个三层循环（行列遍历、相乘累加）优化成了一个两层循环（行列遍历求转置）加一个三层循环（行列遍历、相乘累加），不论怎么看都像是负优化，怎么可能快呢。估计是我理解错了，遂去谷歌。

结果展开目录

结果我真是万没想到，还真是这样的，不光代码是这样，结果也是有优化的。我在 WSL2 里面使用 GCC 8.3.0 编译如下代码：

#include <time.h>
#include <stdio.h>
#include <stdlib.h>

int matrix_size = 200;
int optimized = 1;

int main() {
    // 定义变量
    float *a,*b,*c;
    long int i,j,k,size,m;
    double duration;
    clock_t start, finish;
    size=matrix_size; 
    m=size*size;
    a=(float*)malloc(sizeof(float)*m);
    b=(float*)malloc(sizeof(float)*m);
    c=(float*)malloc(sizeof(float)*m);

    // 初始化
    for(i=0;i<size;i++)
        for(j=0;j<size;j++)
        {
            a[i*size+j]=(float)(rand()%1000/100.0);
            b[i*size+j]=(float)(rand()%1000/100.0);
            c[i*size+j]=(float)(rand()%1000/100.0);
        }

    // 开始计时
    start = clock();
    // 运算
    if (optimized) {
        // 计算转置
        for(i=0;i<size;i++)
            for(j=0;j<size;j++)
            {
                b[i * size + j] = c[ j * size + i];
            }
        // 计算
        for(i = 0;i < size;i++)
            for(j=0;j<size;j++)
            {
                c[i*size+j]=0;
                for(k=0;k<size;k++)
                    c[i*size+j]+=a[i*size+k]*b[j*size+k];
            }
    } else {
        for(i=0;i<size;i++)
            for(j=0;j<size;j++)
            {
                c[i*size+j] = 0;
                for(k=0;k<size;k++)
                    c[i*size+j]+=a[i*size+k]*b[k*size+j];
            }
    }

    // 停止计时
    finish = clock();

    // 计算用时
    duration = (double)(finish - start) / CLOCKS_PER_SEC;
    printf( "%fms\n", duration * 1000);
    return 0;
}

调整 matrix_size 和 optimized，得出了以下数据：

好家伙，我这波直接好家伙，WSL2 的单核性能太差了！（不对

好家伙，我这波直接好家伙，优化前后差距太大了！

后来我又用 Java 试了一下（Java11 AdoptOpenJDK HotSpot ）：

package demo.structure;

import java.util.Random;

public class Main {

    public static final int matrix_size = 3000;
    public static final boolean optimized = true;

    public static void main(String[] args) {
        double[] a = new double[matrix_size * matrix_size];
        double[] b = new double[matrix_size * matrix_size];
        double[] c = new double[matrix_size * matrix_size];

        long start, finish;
        Random random = new Random();


        // 初始化
        for (int i = 0; i < matrix_size; i++) {
            for (int j = 0; j < matrix_size; j++) {
                a[i * matrix_size + j] = random.nextInt(32767) % 1000 / 100.0;
                b[i * matrix_size + j] = random.nextInt(32767) % 1000 / 100.0;
                c[i * matrix_size + j] = random.nextInt(32767) % 1000 / 100.0;
            }
        }

        // 开始计时
        start = System.currentTimeMillis();
        // 运算
        if (optimized) {
            // 计算转置
            for (int i = 0; i < matrix_size; i++) {
                for (int j = 0; j < matrix_size; j++) {
                    b[i * matrix_size + j] = c[j * matrix_size + i];
                }
            }
            // 计算
            for (int i = 0; i < matrix_size; i++) {
                for (int j = 0; j < matrix_size; j++) {
                    c[i * matrix_size + j] = 0;
                    for (int k = 0; k < matrix_size; k++) {
                        c[i * matrix_size + j] += a[i * matrix_size + k] * b[j * matrix_size + k];
                    }
                }
            }
        } else {
            for (int i = 0; i < matrix_size; i++) {
                for (int j = 0; j < matrix_size; j++) {
                    c[i * matrix_size + j] = 0;
                    for (int k = 0; k < matrix_size; k++) {
                        c[i * matrix_size + j] += a[i * matrix_size + k] * b[k * matrix_size + j];
                    }
                }
            }
        }

        // 停止计时
        finish = System.currentTimeMillis();

        // 计算用时
        System.out.println((finish - start) + "ms");
    }
}

没想到离底层硬件那么遥远的 Java 还是受体系结构影响颇深，我挺惊讶的。

惊讶之余我还想讨论一下原理。

原理展开目录

其实原理没什么难的，对于计算机专业的人来说，想必诸如 CPU 缓存（Cache）之类的词都已经听烂了。为了照顾一下非计算机专业的读者，简单的说缓存的作用就是弥补 CPU 与内存的速度差距。通过参数我们可以看得出来，CPU 的频率一般都是 3GHz 左右，低压 U 可能是 2.8GHz，有一些不错的 CPU 甚至能上 5GHz，而内存，就我的笔记本来说还是 2400Hz，将他们换算成时间的话，3.5GHz 的 CPU 工作周期大约是 0.3 纳秒（10 的 - 9 次方），而 2400Hz 的内存工作周期是 0.4 毫秒，是 400000 纳秒。考虑到 DDR 技术，可以认为是 200000 纳秒响应一次 CPU，那内存给 CPU 一次数据，CPU 处理完成后就要空等下一个 200000 纳秒，大约 66 万多 CPU 周期。在这期间 CPU 什么也没法做（因为要等内存给它指令和数据），所以为了提高 CPU 利用率，设计了缓存这个东西。CPU 一次向内存请求一大片连续的数据，除了它要用的，附近的数据也被传送过来，存储在缓存中，之后 CPU 需要使用的时候先去缓存里查找，若有就直接用快得多的缓存，没有再向内存索要数据。

再来看我们的程序，未经优化的算法如下：

for(i=0;i<size;i++)
    for(j=0;j<size;j++)
    {
        c[i*size+j] = 0;
        for(k=0;k<size;k++)
            c[i*size+j]+=a[i*size+k]*b[k*size+j];
    }

可以看到实际的计算最终发生在最内层循环，分别要访问 a、b、c 三个矩阵。对 a 矩阵的访问是 i*size+k，对应来说 i 作为最外层循环的变量，变化是最小的，因此这个访问的偏移量也比较稳定，又由于 k 是连续变化，因此对 a 的访问可以认为是顺序访问。而对 c 的访问，在内层循环就是不变的。而 b 就是另一个故事了：k*size+j，由于 j 是中间层的循环变量，不如 i 稳定，但也比 k 稳定，因此每一次循环对 b 的访问总是能在索引上差出一个 size 的大小。如果 size 足够小，或许能够使得相邻的两次访问命中 Cache，但是下一次就很难说了。如果 size 不够小，每一次访问都没能命中缓存，因而每一次访问都必须向内存请求数据。整个程序的运行速度就被拖慢了。

这一点在现代程序中更为明显：某一时刻下 CPU 并非全力执行你的程序，正相反，CPU 会在 Typora、IDEA、Clion、Word、Notepad++ 等好多个程序之间来回切换，因此某一单个程序的缓存很容易就被调换出去，为新来的程序腾位置。因此有必要让程序访问的数据尽可能呆在一起，使 CPU 更多的命中缓存，减少等待数据的停顿，为程序运行提速。

下面再来看优化后的算法就一目了然了：

// 计算转置
for(i=0;i<size;i++)
    for(j=0;j<size;j++)
    {
        b[i * size + j] = c[ j * size + i];
    }
// 计算
for(i = 0;i < size;i++)
    for(j=0;j<size;j++)
    {
        c[i*size+j]=0;

        for(k=0;k<size;k++)
            c[i*size+j]+=a[i*size+k]*b[j*size+k];
    }

首先计算了一个转置矩阵，这个没问题，然后计算中对 b 的访问因为转置操作而变成了 j*size+k，这样一来对 b 的访问就成为了和 a 一样的顺序访问，因此 CPU 命中缓存的次数大大提高。于是虽然代码写了一个双层循环个一个三层循环，可实际跑起来却比单独的三层循环要快得多。

总结展开目录

现在网络资源的触手可得，使得大家学习编程变得十分容易，基本上可以说会用电脑的，再不济都能学个 Python，并且还能写的不差。这种分界线的模糊一度使我觉得编程这个东西甚至都没什么必要在学校学，自己在家买本书敲敲电脑就学了，哪还用得上老师？实际上我的 Java 就是初中看书自学的。直到今天下午，我还习以为常地认为现代编译器已经能够很好的优化人们写出来的像屎一样的代码了，再之后，这个实验一下子让我耳目一新。

或许是近些年来硬件的发展加上软件基建的优化，使得电脑性能大大提升，就连早年间肉的要死的 Java 现在都能和 C++ 试比高了，在这样的表象之下，认为「代码怎么写都一样，编译器能优化好」的人应当不在少数。今天这个实验虽然不是什么破天荒的东西，但确实让我体会到了受过系统教育的程序员的专业所在：大家都是写代码的，他写出来的代码怎么运行怎么产生效应，这个过程他心里有数，出了问题知道去哪里改去哪里优化，这就是专业。

也许大部分沉迷写 Python 的人一辈子都不会去考虑自己写的代码与 CPU 的底层实现有什么关系，他们只会觉得 Cpython 处理了一切脏活累活，并让他们的代码熠熠生辉。让他们优化性能也无非只是堆叠数学上的最优算法，真到了没办法再提升的地步时，就说性能已经到了极限，电脑只能跑这么快了。

人生苦短，我选 Java。每次写代码都能发现自己的无知与傲慢（那当然还是得感谢无情的编译器和 JVM，又是报错又是抛异常的），每次又都能在学习新知识的过程中感受到充实与快乐。人生如此，夫复何求？

- 全文完 -

---

【代码札记】基于 Cache 的矩阵乘积算法性能改进 由 天空 Blond 采用 知识共享 署名 - 非商业性使用 - 相同方式共享 4.0 国际 许可协议进行许可。
本许可协议授权之外的使用权限可以从 https://skyblond.info/about.html 处获得。