一个矩阵乘法的问题

问题：1024阶双精度浮点矩阵相乘，矩阵满秩

经典代码：

for (i = 0; i < N; i++)
{
    for (j = 0; j < N; j++)
    {
        for (k = 0; k < N; k++)
        {

            c[i*N+j] = c[i*N+j] + a[i*N+k]*b[k*N+j];
        }
    }
}

这是比较经典的方式，发现一个问题，1024阶的时间居然比1025阶的时间多不少，很令人费解，于是加了一些类似gettimeofday这种函数统计计算每一个行，每一个结果的时间，发现在1024的时候没隔一定个数的元素，时间都会有一个比较大的增加（一般是加倍），而1025则不会，个人的猜想是cache造成的影响，1024的情况，每行正好是cache line size的整数倍，而1025则不是，可能会有一些预取，而在这种乘法访问顺序的情况下，会出现“跳着”访问，这样这些预取就会使得性能有所提升。

对原有乘法做一些改变：

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

师兄说这是一个比较经典的矩阵相乘方法，自习想了想才想通，字面上看是把二三层的循环调换了一下顺序，其实把加法均摊，实现了对a数组、b数组都顺序访问，这样就是cache性能提升了很多，整体性能就提升了，而且这样修改后，1025和1024的时间也正常了，阶数多的时间长。

完整代码：

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

//#define N 1023
int main(int argc, char* argv[])
{
    double *a, *b, *c;
    int i, j, k;
    int N;
    struct timeval start, end;
    N = atoi(argv[1]);
    a = (double*)malloc(N*N*sizeof(double));
    b = (double*)malloc(N*N*sizeof(double));
    c = (double*)malloc(N*N*sizeof(double));
    for (i = 0; i < N; i++)
    {
        for (j = 0; j < N; j++)
        {
            a[i*N+j] = 2.0;
            b[i*N+j] = 1.0;
            c[i*N+j] = 0;
        }
    }

    gettimeofday(&start, NULL);

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

    gettimeofday(&end, NULL);
    printf("line %d\t%d\n", i, 1000000*(end.tv_sec - start.tv_sec) + end.tv_usec - start.tv_usec);
    free(a);
    free(b);
    free(c);
    return 0;
}