为什么转置512×512矩阵，会比513×513矩阵慢很多？

原博客：http://blog.jobbole.com/28219/

       很普通的一个求矩阵转置的程序。但是，当MATSIZE取512和513的时候，出现了非常有意思的结果：513竟然比512快。更进一步的研究发现，size=512的时候，运算速度会比同数量级的其它数字慢很多很多。这是怎么一回事呢？

#define SAMPLES 1000
#define MATSIZE 513
 
#include <time.h>
#include <stdio.h>
int mat[MATSIZE][MATSIZE];
 
void transpose()
{
    int aux ;     
    int i,j;       
    for (i = 0 ; i < MATSIZE ; i++ )
        for (j = 0 ; j < MATSIZE ; j++ )
        {
            aux = mat[i][j];
            mat[i][j] = mat[j][i];
            mat[j][i] = aux;
        }
    return;
}
 
int tst_main()
{
   //initialize matrix
    int i,j;  
    int t, elapsed;
    for (i = 0 ; i < MATSIZE ; i++ )
        for (j = 0 ; j < MATSIZE ; j++ )
        {
            mat[i][j] = i+j;
        }
 
    t = clock();
    for (i = 0 ; i < SAMPLES ; i++ )
    {
       transpose();
    }
    elapsed = clock() - t;
    printf("Average for a matrix of:%d:%d\n", MATSIZE, elapsed / SAMPLES);
    return 0;
}
int main()
{
    tst_main();
}

      很容易就联想到，造成这个问题的原因是CPU Cache.原作者没有给出他的CPU型号，但是如今的pc几乎都是采用的set associative的cache结构，下面用64bits CPU 的2-way set associate来做例子，讲解一下cache的工作原理。

      CPU 2-way set associate Cache结构如上图所示。Cache由Cache Set组成，2-way 表示每个Set里面有2条记录（Cache Line），每条记录包括16个字节（如上图）。CPU得遍历这两条记录，比较前面50位的Tag，如果Tag一样，并且Valid bit(V)=1，那么恭喜你，你的数据在Cache里面，如果遍历完这两条记录，还是没有找到Tag的话，那么很遗憾，你的数据不在Cache里，得从内存里读。从内存里获取相应的数据，然后把它存到对应的Cache Set里，如果Set里有空位的话最好，如果没有的话，用LRU来替换。因为一个Set里只有2条数据，所以实现LRU仅仅需要一个额外bit就可以了，非常高效。

      依据cache的这种结构，一个内存地址（Memory address），可以划分为tag，block，word，byte 4个部分，其中CPU Cache的大小决定了block的位数。如上图10bits的block，对应了1024个Cache Set(共计8KB)，内存地址的block固定了，就必须存储在相应的Set里面，这样可以把查询cache的事件从O(n)缩短为O(1)。

可能有人会觉得奇怪，为什么block不是取的最前面的10bits，这当然是有道理的，通常在内存里数据都是连续存放的，就是说，同一段程序用的数据，他们前10位几乎都是一样的，如果用前10位来定位block，那么collision的发生率非常高，cache效率非常低下，所以才选了后面的10位来定位block。

      举个例子，block是1023（1111111111），你的数据就放在第1023个set里面。CPU首先遍历第1023个Cache Set里的两条记录，比较前面50位的tag，如果tag一样，并且Valid bit(V)=1，那么说明数据在Cache里面，接着就可以通过内存地址的word和byte域来读取Cache Line（16Bytes）中的数据了。否则得从内存读取数据。

好了，背景知识介绍的差不多了，让我们回到这个问题上来。为什么512大小的矩阵，会比其它数字慢那么多？

      让我们来计算一下，512×512的int矩阵，在内存里是连续存放的。每个cache line是16bytes，对应4个int，所以一个n阶矩阵的row可以填充n/4个cache set。假设第一个数据a[0][0]正好对应cache set 0，那么其中每一个数据a[i][j]对应的cache set是(512*i+j)/4%1024=(128i+j/4)%1024。可以看到，前面的系数正好可以整除。很不巧的是，在进行矩阵转置的运算时，在第2个for循环中，我们需要依次访问每一个row中对应i的值。这样会造成下面的结果：假设i=0，

set(a[0][0])=0;

set(a[1][0])=128;

set(a[2][0])=256;

set(a[3][0])=384;

set(a[4][0])=512;

set(a[5][0])=640;

set(a[6][0])=768;

set(a[7][0])=896;

set(a[8][0])=0;

set(a[9][0])=128;

…

set(a[15][0])=896;

      到a[15][0]的时候刚好填完整个cache的所有128整数倍的set，当读取a[16][0]的时候，将会发生replace，把a[0][0]从cache里移除。这样，当源程序的i=1时，将完全重复i=0的计算过程，每次取数据都需要先从memory读到cache中来，cache的作用完全没有体现。

而当size=513的时候，事情就不一样了，mat[i][j]对应的cache set是(513*i+j)/4%1024，前面的系数除不尽了，每递增4次结果会比size=512时偏差1,例如：

set(a[0][0])=0;

set(a[1][0])=128;

set(a[2][0])=256;

set(a[3][0])=384; 

set(a[4][0])=513;

set(a[5][0])=641;

…

      这样就很微妙的把Cache Set给错开了。a[16][0]不在第0行而是第4行，不会覆盖之前的数据。即使将全部的a[0-15][i]都读入cache，也不会发生碰撞。之后，由于一个cache有4个word，a[0-15][i+1]，a[0-15][i+2]，a[0-15][i+3]也同时被读进cache里了，所以计算i+1,i+2,i+3时，仅仅需要读对应行的数据就可以了，同一行的数据都是连续的，所以碰撞率很低。这个计算过程很好的利用了cache，如果不考虑其他因素（实际上，这个已经是影响运行时间的最大因素了），理论上我们可以节省75%的运行时间，可以看到，这个理论预测是和提问者给的数据相符合的。

      总之，当你的data size是128的整数倍的时候，得特别小心，搞不好cache collision就把你的程序给拖慢了呢。