通过矩阵乘法看内存访问对CPU运算速度的影响

关于Intel C++编译器和Visual C++编译器的差异块可见“Intel和Microsoft C++编译器在矩阵乘法测试例子中运行时间的差异”，从速度上考量这里仅测试Intel C++编译器的情形。矩阵乘法有普通的按定义的方法和块方法，测试结果表明后者可达到前者的两倍速度。速度和“通过加法运算看内存访问对CPU运算速度的影响”中的结果差不多。程序代码如下：

#include<iostream>
#include<ctime>
#include<cmath>
using namespace std;
template<class T,int TILE_WIDTH>
void MatrixMulTile(T *A,T *B, T *C, int Width)
{
T As[TILE_WIDTH][TILE_WIDTH];
T Bs[TILE_WIDTH][TILE_WIDTH];
T Cvalue(0);
for(int by=0;by<Width/TILE_WIDTH;++by)
for(int bx=0;bx<Width/TILE_WIDTH;++bx)
{
for(int ty=0;ty<TILE_WIDTH;++ty)
for(int tx=0;tx<TILE_WIDTH;++tx)
{
C[(by*TILE_WIDTH+ty)*Width+bx*TILE_WIDTH+tx] = 0;
}
for(int m=0;m<Width/TILE_WIDTH;++m)
{
for(int ty=0;ty<TILE_WIDTH;++ty)
for(int tx=0;tx<TILE_WIDTH;++tx)
{
As[ty][tx]=A[by*TILE_WIDTH*Width+m*TILE_WIDTH+ty*Width+tx];
Bs[tx][ty]=B[m*TILE_WIDTH*Width+bx*TILE_WIDTH+ty*Width+tx];
}
for(int ty=0;ty<TILE_WIDTH;++ty)
for(int tx=0;tx<TILE_WIDTH;++tx)
{
Cvalue=C[(by*TILE_WIDTH+ty)*Width+bx*TILE_WIDTH+tx];
for(int k=0;k<TILE_WIDTH;++k)
Cvalue+=As[ty][k]*Bs[tx][k];
C[(by*TILE_WIDTH+ty)*Width+bx*TILE_WIDTH+tx]=Cvalue;
}
}
}
}
template<class T>
void MatrixMul(const T *A, const T *B, T *C, int Width) {
int i, j, k;
for(i=0; i<Width; i++)
for(j=0; j<Width; j++){
T s=0;
for(k=0; k<Width; k++)
s+=A[i*Width+k]*B[k*Width+j];
C[i*Width+j]=s;
}
}

template<class T,int TILE_WIDTH>
void test(int Width)
{
T *h_A=(T*)malloc(Width*Width*sizeof(T));
T *h_B=(T*)malloc(Width*Width*sizeof(T));
T *h_C=(T*)malloc(Width*Width*sizeof(T));
T *h_C_ref=(T*)malloc(Width*Width*sizeof(T));
T error_norm=0, ref_norm=0;

for(int i=0; i<Width*Width; i++) {
h_A[i]=rand()/(T)RAND_MAX;
h_B[i]=rand()/(T)RAND_MAX;
}
clock_t t0=clock();
MatrixMul<T>(h_A, h_B, h_C_ref,Width);
float cpu_t1=static_cast<float>(clock()-t0)/1000.0;

t0=clock();
MatrixMulTile<T,TILE_WIDTH>(h_A,h_B,h_C,Width);
float cpu_t2=static_cast<float>(clock()-t0)/1000.0;

error_norm=0;
ref_norm=0;
for(int i=0; i<Width*Width; i++) {
T diff=h_C_ref[i]-h_C[i];
error_norm+=diff*diff;
ref_norm+=h_C_ref[i]*h_C_ref[i];
}
printf("Test %s/n", (sqrt((error_norm+0.0)/ref_norm)<1E-6) ? "PASSED" : "FAILED");

if(sizeof(T)==4)
printf("单精度:/n");
else
printf("双精度:/n");
printf("矩阵阶数为%4d,", Width);
printf("普通方法: %.6fs(%.3fGflops),", cpu_t1, 1e-9*Width*Width*Width*2/cpu_t1);
printf("%dx%d块方法: %.6fs(%.3fGflops)/n", TILE_WIDTH,TILE_WIDTH,cpu_t2, 1e-9*Width*Width*Width*2/cpu_t2);
free(h_A);
free(h_B);
free(h_C);
free(h_C_ref);
}
int main(int argc,char*argv[])
{
int Width=atoi(argv[1]);
const unsigned int TILE_WIDTH=4;
test<float,TILE_WIDTH<<1>(Width);
test<float,TILE_WIDTH<<2>(Width);
test<float,TILE_WIDTH<<3>(Width);
test<float,TILE_WIDTH<<4>(Width);
test<float,TILE_WIDTH<<5>(Width);
test<double,TILE_WIDTH<<0>(Width);
test<double,TILE_WIDTH<<1>(Width);
test<double,TILE_WIDTH<<2>(Width);
test<double,TILE_WIDTH<<3>(Width);
test<double,TILE_WIDTH<<4>(Width);
}

运行结果如下:

规模为640

Test PASSED
单精度:
矩阵阶数为 640,普通方法: 0.188000s(2.789Gflops),8x8块方法: 0.250000s(2.097Gflops)
Test PASSED
单精度:
矩阵阶数为 640,普通方法: 0.172000s(3.048Gflops),16x16块方法: 0.125000s(4.194Gflops)
Test PASSED
单精度:
矩阵阶数为 640,普通方法: 0.203000s(2.583Gflops),32x32块方法: 0.078000s(6.722Gflops)
Test PASSED
单精度:
矩阵阶数为 640,普通方法: 0.187000s(2.804Gflops),64x64块方法: 0.078000s(6.722Gflops)
Test PASSED
单精度:
矩阵阶数为 640,普通方法: 0.187000s(2.804Gflops),128x128块方法: 0.079000s(6.637Gflops)
Test PASSED
双精度:
矩阵阶数为 640,普通方法: 0.297000s(1.765Gflops),4x4块方法: 0.343000s(1.529Gflops)
Test PASSED
双精度:
矩阵阶数为 640,普通方法: 0.297000s(1.765Gflops),8x8块方法: 0.281000s(1.866Gflops)
Test PASSED
双精度:
矩阵阶数为 640,普通方法: 0.313000s(1.675Gflops),16x16块方法: 0.156000s(3.361Gflops)
Test PASSED
双精度:
矩阵阶数为 640,普通方法: 0.312000s(1.680Gflops),32x32块方法: 0.141000s(3.718Gflops)
Test PASSED
双精度:
矩阵阶数为 640,普通方法: 0.297000s(1.765Gflops),64x64块方法: 0.156000s(3.361Gflops)

规模为1280

Test PASSED
单精度:
矩阵阶数为1280,普通方法: 1.328000s(3.158Gflops),8x8块方法: 3.187000s(1.316Gflops)
Test PASSED
单精度:
矩阵阶数为1280,普通方法: 1.313000s(3.194Gflops),16x16块方法: 1.437000s(2.919Gflops)
Test PASSED
单精度:
矩阵阶数为1280,普通方法: 1.328000s(3.158Gflops),32x32块方法: 0.765000s(5.483Gflops)
Test PASSED
单精度:
矩阵阶数为1280,普通方法: 1.360000s(3.084Gflops),64x64块方法: 0.703000s(5.966Gflops)
Test PASSED
单精度:
矩阵阶数为1280,普通方法: 1.391000s(3.015Gflops),128x128块方法: 0.718000s(5.842Gflops)
Test PASSED
双精度:
矩阵阶数为1280,普通方法: 2.750000s(1.525Gflops),4x4块方法: 7.360000s(0.570Gflops)
Test PASSED
双精度:
矩阵阶数为1280,普通方法: 2.688000s(1.560Gflops),8x8块方法: 3.859000s(1.087Gflops)
Test PASSED
双精度:
矩阵阶数为1280,普通方法: 2.641000s(1.588Gflops),16x16块方法: 1.796000s(2.335Gflops)
Test PASSED
双精度:
矩阵阶数为1280,普通方法: 2.625000s(1.598Gflops),32x32块方法: 1.375000s(3.050Gflops)
Test PASSED
双精度:
矩阵阶数为1280,普通方法: 2.687000s(1.561Gflops),64x64块方法: 1.234000s(3.399Gflops)