利用高速缓存(Cache)的局部性优化矩阵乘法

-1-一、实验目的与要求：实验目的：1.增进对cache工作原理以及计算机存储体系的理解；2.体验程序中访存模式变化是如何影响cahce效率进而影响程序性能的过程；实验要求：1.自学课本高速缓存章节的相关知识；2.分析普通矩阵乘法访存时高速缓存缺失与矩阵大小n和cache容量的关系；3.将矩阵乘法分块实现来充分利用高速缓存的局部性，优化程序性能；4.测试分块版本矩阵计算程序在矩阵不同大小（32，64，128，160，480，960，….MAX），时与未分块版本的性能对比。注：最大值MAX使得矩阵尺寸增大到不能在cache中完全容纳这三个矩阵时值（计算机cache大小可自行查找相关处理器的资料或用相关测试软件）5.调整BLOCKSIZE大小，分析BLOCKSIZE大小对程序性能的影响；6.制作相关数据的对比图表，分析原因并在blackboard上提交实验报告。二、实验环境：硬件：桌面PC软件：linux（GCC采用“inline”功能去除调用开销）三、实验内容：按照下面的实验步骤及说明，完成相关实验并提交实验报告：（1）如下，普通矩阵乘法一般采用3重循环完成。voiddgemm(intn,double*A,double*B,double*C){for(inti=0;in;++i)for(intj=0;jn;++j){doublecij=C[i+j*n];/*cij=C[i][j]*/for(intk=0;kn;k++)cij+=A[i+k*n]*B[k+j*n];/*cij+=A[i][k]*B[k][j]*/C[i+j*n]=cij;/*C[i][j]=cij*/}}（2）分析普通矩阵乘法访存时高速缓存缺失与矩阵大小n和cache容量的关系。-2-（3）将矩阵乘法分块实现来充分利用高速缓存的局部性，优化程序性能。#defineBLOCKSIZE32voiddo_block(intn,intsi,intsj,intsk,double*A,double*B,double*C){}voiddgemm_block(intn,double*A,double*B,double*C){for(intsj=0;sjn;sj+=BLOCKSIZE)for(intsi=0;sin;si+=BLOCKSIZE)for(intsk=0;skn;sk+=BLOCKSIZE)do_block(n,si,sj,sk,A,B,C);}完成do_block子程序的设计和相关性能测试程序的代码。（4）测试分块版本矩阵计算程序在矩阵不同大小（32，64，128，160，480，960，….MAX），时与未分块版本的性能对比。注：最大值MAX使得矩阵尺寸增大到不能在cache中完全容纳这三个矩阵时值（计算机cache大小可自行查找相关处理器的资料或用相关测试软件）。（5）调整BLOCKSIZE大小，分析BLOCKSIZE大小对程序性能的影响。（6）制作相关数据的对比图表，分析原因并提交实验报告。四、实验步骤与过程：(给出程序分析和算法描述（流程图或文字）、程序核心代码。)实验步骤：1.分析普通矩阵乘法访存时高速缓存缺失与矩阵大小n和cache容量的关系：普通矩阵乘法：voiddgemm(intn,double*A,double*B,double*C){for(inti=0;in;++i)for(intj=0;jn;++j){doublecij=C[i+j*n];/*cij=C[i][j]*/for(intk=0;kn;k++)cij+=A[i+k*n]*B[k+j*n];/*cij+=A[i][k]*B[k][j]*/C[i+j*n]=cij;/*C[i][j]=cij*/}}-3-分析：计算机在实际计算上述普通矩阵乘法时，所计算矩阵C的每一个数据时，都要用到矩阵A的某行和矩阵B中的某列，而矩阵A、B和C都是存储在内存中的，又由于CPU的速度远远大于访问内存的速度，如果是直接从内存读取和写回计算数据，那么计算效率是非常低下的，由于访问内存会导致时延，CPU的计算资源被浪费，即计算效率低。为了提高计算速度，引入了cache机制，即先把存放在内存中的矩阵A、B的元素调入cache，这样寄存器可以先寻访cache，访问cache的速度要比访问内存的速度快，如果在cache中没有所需要的数据时，才需要访问内存。但是，矩阵A、B在实际应用中都包含大量的元素，数据量非常分庞大，也即，上述程序中n很大，而处理器中的cache往往很小，因此不能将整个矩阵全部放入cache中。因此需要将这些大的矩阵按照某种方法进行分块，使得分块后的小矩阵可以放入到cache中，但是分块又不能随意分，需要有一定的原则去分块，如果分块子矩阵太大，那么子矩阵还是不能全部放入cache中，如果分块子矩阵太小，那么为了计算一个大矩阵的数据，需要调入cache的子矩阵的次数会增加，因此需要选择合适的分块方法。2.分块实现矩阵乘法，利用cache的局部性，优化程序性能：a)安装Linux系统：b)查看Linux系统cache的大小：c)编写普通矩阵乘法运算的代码、利用cache局部性分块实现矩阵相乘的代码，如下：#includestdio.h#includememory.h#includetime.h#includestdlib.h#includesys/time.h-4-inlinevoiddgemm(intn,intblocksize,double*A,double*B,double*C){inti,j,k;for(i=0;iblocksize;i++){for(j=0;jblocksize;j++){doublecij=C[i*n+j];/*cij=C[i][j]*/for(k=0;kblocksize;k++)cij+=A[i*n+k]*B[k*n+j];/*cij+=A[i][k]*B[k][j]*/C[i*n+j]=cij;/*C[i][j]=cij*/}}}inlinevoiddo_block(intn,intblocksize,intsi,intsj,intsk,double*A,double*B,double*C){dgemm(n,blocksize,A+si*n+sk,B+sk*n+sj,C+si*n+sj);//printf(\n);//printf(%d%d%d\n,si,sj,sk);//for(inti=0;in;i++)//{//for(intj=0;jn;j++)//{//printf(%lf,C[si*n+sj]);//}//printf(\n);//}//printf(\n);}inlinevoiddgemm_block(intn,intblocksize,double*A,double*B,double*C){intsj,si,sk;for(sj=0;sjn;sj+=blocksize)for(si=0;sin;si+=blocksize)for(sk=0;skn;sk+=blocksize)-5-do_block(n,blocksize,si,sj,sk,A,B,C);}intmain(){intsizes[9]={64,128,256,480，512，960，1024，1536,1920};inti,j,th;for(th=0;th9;th++){intBLOCKSIZE=16;intN=sizes[th];if(N%BLOCKSIZE!=0)continue;double*A=NULL,*B=NULL,*C=NULL;A=(double*)malloc(sizeof(double)*N*N);B=(double*)malloc(sizeof(double)*N*N);C=(double*)malloc(sizeof(double)*N*N);if(C==NULL||A==NULL||B==NULL){printf(mallocerror!);return1;}memset(C,0,sizeof(double)*N*N);for(i=0;iN;i++){for(j=0;jN;j++){A[i*N+j]=1;B[i*N+j]=2;}}/********************计时*****************/structtimevalstartTime,endTime;floatTimeuse;gettimeofday(&startTime,NULL);/********************计时*****************///dgemm(N,N,A,B,C);dgemm_block(N,BLOCKSIZE,A,B,C);/********************计时*****************/gettimeofday(&endTime,NULL);Timeuse=1000000*(endTime.tv_sec-startTime.tv_sec)+(endTime.tv_usec--6-startTime.tv_usec);printf(meantimeuse=%fus\n,Timeuse/N/N/N);/********************计时*****************/printf(N:%dblocksize:%d\n,N,BLOCKSIZE);printf(%d%lf%lf\n\n,N,C[0],C[N*(N-1)+N-1]);free(A);free(B);free(C);}printf(successful!\n);getchar();return0;}intmain2(){intN=2048;intsizes[9]={64,128,256,480，512，960，1024，1536,1920};inti,j,th;for(th=0;th9;th++){intBLOCKSIZE=N/sizes[th];double*A=NULL,*B=NULL,*C=NULL;A=(double*)malloc(sizeof(double)*N*N);B=(double*)malloc(sizeof(double)*N*N);C=(double*)malloc(sizeof(double)*N*N);if(C==NULL||A==NULL||B==NULL){printf(mallocerror!);return1;}memset(C,0,sizeof(double)*N*N);for(i=0;iN;i++){for(j=0;jN;j++){A[i*N+j]=1;B[i*N+j]=2;}}-7-/********************计时*****************/structtimevalstartTime,endTime;floatTimeuse;gettimeofday(&startTime,NULL);/********************计时*****************///doublestart=clock_t();dgemm(N,N,A,B,C);//dgemm_block(N,BLOCKSIZE,A,B,C);/*doubleend=clock_t();printf(%lfsec\n,(end-start)/CLOCKS_PER_SEC*1000/N/N);*//********************计时*****************/gettimeofday(&endTime,NULL);gettimeofday(&endTime,NULL);gettimeofday(&endTime,NULL);gettimeofday(&endTime,NULL);gettimeofday(&endTime,NULL);gettimeofday(&endTime,NULL);Timeuse=1000000*(endTime.tv_