矩阵乘法-分块计算

将整个矩阵分解为这样的小块，每次完成一对小块的计算，以提高Cache的命中率。
提示：



图中n=N/m
计算次序为A11*B11, A11*B12,…, A11*B1n,，由于反复使用A11，因此可以提高Cache的命中率。

/*
矩阵分开计算
C=A*B  --- C(i,j)等于A的第i行乘以第j列
*/
#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include <math.h>
#include <windows.h>

/*
生成n*n矩阵
*/
void GenerateMatrix(float *m, int n);
void PrintMatrix(float *p, int n);
void GeneralMul(float *A, float *B, float *C, int n);
void ClearMatrix(float *m, int n);
/*
矩阵分块计算
*/
void BlockCacul(float *A, float *B, float *C, int n, int thread_num, int m);
/*
两个矩阵的误差
*/
float diff(float *C1, float *C0, int n);

struct ARG {
float *A;
int ax, ay;
float *B;
int bx, by;
float *C;
int cx, cy;
int m;
int n;
};
int main(int argc, char **argv)
{
if (argc != 4)
{
printf("Usage: %s N thread_num M\n", argv[0]);
return 0;
}
int n=atoi(argv[1]);
int thread_num = atoi(argv[2]);
int m = atoi(argv[3]);
float *A = new float[n*n];
float *B = new float[n*n];
float *C = new float[n*n];
float *C0 = new float[n*n];

GenerateMatrix(A, n);
GenerateMatrix(B, n);

clock_t start;
float time_used;

ClearMatrix(C0, n);
start=clock();
GeneralMul(A, B, C0, n);
time_used = static_cast<float>(clock() - start)/CLOCKS_PER_SEC*1000;
printf("General:   time = %f\n", time_used);

ClearMatrix(C, n);
start=clock();
BlockCacul(A, B, C, n, thread_num, m);
time_used = static_cast<float>(clock() - start)/CLOCKS_PER_SEC*1000;
printf("Block:  time = %f\n", time_used);
printf("Difference of two result: %f\n", diff(C0, C, n));

delete [] A;
delete [] B;
delete [] C;
delete [] C0;
return 0;
}

void ClearMatrix(float *m, int n)
{
for (int i=0; i<n; i++)
{
for (int j=0; j<n; j++)
m[i*n+j]=0.0;
}
}
/*
普通矩阵相乘
*/
void GeneralMul(float *A, float *B, float *C, int n)
{
for (int i=0; i<n; i++)
{
for (int j=0; j<n; j++)
{
float *p = C+i*n+j;
for (int k=0; k<n; k++)
{
*p += A[i*n+k]*B[k*n+j];
}
}
}
}

DWORD WINAPI Mul_Fun(LPVOID arg)
{
struct ARG *p = (struct ARG *)arg;

float *A = p->A;
float *B = p->B;
float *C = p->C;
int m = p->m;
int n = p->n;
for (int i=0; i<m; i++)
{
for (int j=0; j<m; j++)
{
float *t = C+(i+p->cx)*n+p->cy+j;

for (int k=0; k<m; k++)
{
*t += A[(p->ax+i)*n+p->ay+k]*B[(p->bx+k)*n+p->by+j];
}
}
}

return 0;
}

void BlockCacul(float *A, float *B, float *C, int n, int thread_num, int m)
{
//m = static_cast<int>(sqrt(m));
struct ARG *args = new struct ARG[thread_num];
HANDLE *handle = new HANDLE[thread_num];
int t=0;
int i;
for (i=0; i<thread_num; i++)
{
args[i].A = A;
args[i].B = B;
args[i].C = C;
args[i].m = m;
args[i].n = n;
}
//分成n/m x n/m块
//A i行j列
for (i=0; i<n; i+=m)
{
for (int j=0; j<n; j+=m)
{
//B j行k列
for (int k=0; k<n; k+=m)
{
args[t].ax = i;
args[t].ay = j;
args[t].bx = j;
args[t].by = k;
args[t].cx = i;
args[t].cy = k;
if (t<thread_num)
{
handle[t] = CreateThread(NULL, 0, Mul_Fun, (LPVOID)(&args[t]), 0, 0 );
t++;
}
if (t==thread_num)
{
for (int ii=0; ii<t; ii++)
WaitForMultipleObjects(thread_num, &handle[ii],TRUE,INFINITE);
t = 0;
}
}
}
}

}
void GenerateMatrix(float *p, int n)
{
srand(time(NULL)+rand());
for (int i=0; i<n*n; i++)
{
*p = static_cast<float>(rand())/ (static_cast<float>(rand())+ static_cast<float>(0.55));
p++;
}
}

float diff(float *C1, float *C0, int n)
{
float rst=0.0;
float t;

for (int i=0; i<n; i++)
{
for (int j=0; j<n; j++)
{
t = C1[i*n+j]-C0[i*n+j];
if (t<0)
t = -t;
rst += t;
}
}
return rst;
}

void PrintMatrix(float *p, int n)
{
for (int i=0; i<n; i++)
{
for (int j=0; j<n; j++)
{
printf("%.2f\t", p[i*n+j]);
}
printf("\n");
}
printf("\n");
}