矩阵乘法的优化

源于微博上一个同学发了一条微博：“如何对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  ---------- 

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

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