MENU

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

December 11, 2020 • 瞎折腾

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

背景

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

其中第二个实验是「基于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_sizeoptimized,得出了以下数据:

矩阵大小未优化耗时(ms)优化耗时(ms)加速比
20044.02335.9861.223
500698.498545.8911.280
10007554.3174364.0171.731
150032430.60214646.7102.214
200077378.99835082.5802.206
2500204543.49767557.6533.028
3000343537.171116164.1932.957

好家伙,我这波直接好家伙,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");
    }
}
矩阵大小未优化耗时(ms)优化耗时(ms)加速比
20019200.95
5001971751.126
1000286214172.020
15001098148082.284
200059386142914.155
2500132121283924.653
3000258684477865.413

没想到离底层硬件那么遥远的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 处获得。

Last Modified: December 17, 2020
Archives QR Code
QR Code for this page
Tipping QR Code
Leave a Comment

已有 1 条评论
  1. 勘误:原文中的内存频率应该是2400MHz,后续的周期也要跟着变动。但大体思路是正确的:内存跟不上CPU的速度,因此需要Cache来暂存从内存获取的数据,以备后续需要时能够快速响应。(毕竟Cache是在CPU内部,而内存是放在主板上的)