矩阵乘法的优化

源于微博上一个同学发了一条微博：“如何对C语言下的矩阵乘法进行优化？我试了几种小技巧，但是发现怎么试都远远比不上MATLAB。”

然后他贴了一下他的基本实现方法如下图：



基于他的方法，我也做了一下程序实现，代码如下：

/*
* =====================================================================================
*
*       Filename:  matrix_multi.cpp
*
*    Description:
*
*        Version:  1.0
*        Created:  2013年01月15日 13时43分36秒
*       Revision:  none
*       Compiler:  gcc
*
*         Author:  Weidong Yang (cn), ywdong@mail2.sysu.edu.cn
*        Company:
*
* =====================================================================================
*/

#include <cstdio>
#include <cstdlib>
#include <ctime>

#define MAXNUM 1024
#define BASE 10
#define BLOCK (MAXNUM / 8)

double matrix[MAXNUM][MAXNUM];
double res[MAXNUM][MAXNUM];

double matrix_oneD[MAXNUM * MAXNUM];
double res_oneD[MAXNUM * MAXNUM];

char fileName[] = "Matrix_Data.txt";
clock_t start,  finish;
double duration;
bool isPrint = false;

void GenData()
{
srand(time(NULL));
FILE* fp = fopen(fileName, "w+");
if (fp)
{
for (int i = 0; i < MAXNUM; i++)
{
for (int j = 0; j < MAXNUM; j++)
{
double value = (double)rand() / 100000;
fprintf(fp, "%lf\t", value);
}
fprintf(fp, "\n");
}
}
fclose(fp);
}

void ReadData()
{
FILE* fp;
int index = 0;
if ((fp = fopen(fileName, "r")) != NULL)
{
for (int i = 0; i < MAXNUM; i++)
{
for (int j = 0; j < MAXNUM; j++)
{
double tmp;
fscanf(fp, "%lf", &tmp);
matrix[i][j] = tmp;
matrix_oneD[index++] = tmp;
}
}
printf("Data num is %d\n", index);
}
else
{
printf("Read file %s failed.\n", fileName);
}
fclose(fp);
}

void MM_Normal()
{
int i,j,k;
for (i = 0; i < MAXNUM; i++)
for (j = 0; j < MAXNUM; j++)
res[i][j] = 0.0;

start = clock();

for (i = 0; i < MAXNUM; i++)
for (j = 0; j < MAXNUM; j++)
for (k = 0; k < MAXNUM; k++)
res[i][j] += matrix[i][k] * matrix[k][j];

finish = clock();
duration = (double)(finish - start) / CLOCKS_PER_SEC;
printf("Normal method costs %lf seconds.\n", duration);
}

void MM_Block()
{
int i,j,k,i0,j0,k0;
for (i = 0; i < MAXNUM; i++)
for (j = 0; j < MAXNUM; j++)
res[i][j] = 0.0;

start = clock();

for (i = 0; i < MAXNUM; i += BLOCK)
for (j = 0; j < MAXNUM; j += BLOCK)
for (k = 0; k < MAXNUM; k += BLOCK)
for (i0 = i; i0 < (i + BLOCK); i0++)
for (j0 = j; j0 < (j + BLOCK); j0++)
for (k0 = k; k0 < (k + BLOCK); k0++)
res[i0][j0] += matrix[i0][k0] * matrix[k0][j0];

finish = clock();
duration = (double)(finish - start) / CLOCKS_PER_SEC;
printf("Block method costs %lf seconds.\n", duration);

}

void MM_Transpose()
{
int i,j,k;
for (i = 0; i < MAXNUM; i++)
for (j = 0; j < MAXNUM; j++)
res[i][j] = 0.0;

start = clock();

for (i = 0; i < MAXNUM; i++)
for (k = 0; k < MAXNUM; k++)
for (j = 0; j < MAXNUM; j++)
res[i][j] += matrix[i][k] * matrix[k][j];

finish = clock();
duration = (double)(finish - start) / CLOCKS_PER_SEC;
printf("Transpose method costs %lf seconds.\n", duration);
}

void MM_T_OneDimension()
{
int i,j,k;
for (i = 0; i < MAXNUM * MAXNUM; i++)
res_oneD[i] = 0.0;

start = clock();

for (i = 0; i < MAXNUM; i++)
for (k = 0; k < MAXNUM; k++)
for (j = 0; j < MAXNUM; j++)
res_oneD[i * MAXNUM + j] += matrix_oneD[i * MAXNUM + k] * matrix_oneD[k * MAXNUM + j];

finish = clock();
duration = (double)(finish - start) / CLOCKS_PER_SEC;
printf("Transpose OneDimention method costs %lf seconds.\n", duration);

}

void MM_Transpose_Block()
{
int i,j,k,i0,j0,k0;
for (i = 0; i < MAXNUM; i++)
for (j = 0; j < MAXNUM; j++)
res[i][j] = 0.0;

start = clock();

for (i = 0; i < MAXNUM; i += BLOCK)
for (k = 0; k < MAXNUM; k += BLOCK)
for (j = 0; j < MAXNUM; j += BLOCK)
for (i0 = i; i0 < (i + BLOCK); i0++)
for (k0 = k; k0 < (k + BLOCK); k0++)
for (j0 = j; j0 < (j + BLOCK); j0++)
res[i0][j0] += matrix[i0][k0] * matrix[k0][j0];

finish = clock();
duration = (double)(finish - start) / CLOCKS_PER_SEC;
printf("Transpose Block method costs %lf seconds.\n", duration);

}

void PrintResOneD()
{
if (!isPrint)
return;
int i,j;
for (i = 0; i < MAXNUM; i++)
{
for (j = 0; j < MAXNUM; j++)
printf("%lf ", res_oneD[i * MAXNUM + j]);
printf("\n");
}
printf("\n");
}

void PrintRes()
{
if (!isPrint)
return;
int i,j;
for (i = 0; i < MAXNUM; i++)
{
for (j = 0; j < MAXNUM; j++)
printf("%lf ", res[i][j]);
printf("\n");
}
printf("\n");
}

int main ( int argc, char *argv[] )
{
ReadData();

MM_Normal();
PrintRes();

MM_Transpose();
PrintRes();

MM_T_OneDimension();
PrintResOneD();

MM_Block();
PrintRes();

MM_Transpose_Block();
PrintRes();

return 0;
}		// ----------  end of function main  ----------

在我机子上实验结果如下：



（未完待续......）