从矩阵乘法的不同计算方式来看局部性原理

今天碰到的关于矩阵乘法不同情况下运算速度的问题，隐约记得是因为缓存的问题，后来突然想起来CSAPP那本书上讲过这个东西的，就是通过矩阵乘法三重循环的不同顺序来讲的局部性原理的，所以翻过来又看了一下。

两个矩阵A,B相乘得到C【为了方便起见，把它们都看成n*n的方阵】经典的做法就是用三重循环来实现，但是具体这三重循环如何进行排列，就非常的有讲究。
假设n是一个非常大的数,也就意味着如果跨行的话,必然不会缓存命中,假设只有一个告诉缓存,其块大小为32字节,也就说一个块只能存4个double

经典做法ijk和jik

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

这两种方法的复杂度分析是一样的.AB之中必有一个是每次都不命中的，剩下一个每四个不命中一次。

每次迭代A不命中的次数 每次迭代B不命中的次数 每次迭代C不命中的次数 总数 0.25 1 0 1.25 1 0.25 0 1.25

jki和kji版本

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

这两种方法是最差的方法，因为每次c和a两个都完全不命中

每次迭代A不命中的次数 每次迭代B不命中的次数 每次迭代C不命中的次数 总数 1 0 1 2

最好的方法kij和ikj版，即让最频繁变动的层不动

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

这样最里面的循环那层，每次访问都保证了空间局部性

每次迭代A不命中的次数 每次迭代B不命中的次数 每次迭代C不命中的次数 总数 0 0.25 0.25 0.5

所以可以看到，虽然矩阵乘法非常的简单，但是具体实现起来，其优化所要做到的细节还是非常多的，涉及到访存的问题还是比较复杂的，尤其是对于局部性的理解。