矩阵相乘并行算法.docx
- 文档编号:6513832
- 上传时间:2023-01-07
- 格式:DOCX
- 页数:29
- 大小:278.06KB
矩阵相乘并行算法.docx
《矩阵相乘并行算法.docx》由会员分享,可在线阅读,更多相关《矩阵相乘并行算法.docx(29页珍藏版)》请在冰豆网上搜索。
矩阵相乘并行算法
并行处理技术
课程设计分析报告
课程设计题目
矩阵相乘并行算法设计
姓名
廖杰
学号
M201372880
专业
计算机技术
任课教师
金海石宣化
所在学院
计算机科学与技术学院
报告提交日期
2014-01-13
一、实验目的
1、学习使用集群;
2、掌握并行处理或分布计算的编程方法;
3、学会以并行处理的思想分析问题。
二、实验要求
1、自行生成矩阵作为算法的输入;
2、使用并行处理技术编程,例如:
MPI、OpenMP、MR;
3、矩阵大小至少为1000*1000;
4、加速比越大成绩越高。
三、实验内容
3.1、矩阵的划分:
对于矩阵相乘的并行算法,可以有三种:
对矩阵按行划分、按列划分和棋盘式分块划分。
和按行或列划分相比,棋盘式划分可以开发出更高的并行度。
对于一个n×n的方阵,棋盘划分最多可以使用n^2个处理器进行并行计算,但使用按行或列分解最多可以使用n个。
对矩阵相乘采用棋盘式划分的算法通常称作Cannon算法。
A)行列划分
又叫带状划分(StripedPartitioning),就是将矩阵整行或者整列分成若干个组,每个组指派给一个处理器。
下图所例为4个CPU,8×8矩阵的带状划分。
在带状划分情况下,每个CPU将会均匀分配到2行(列)数据。
8×8矩阵变成了一个1×4或4×1的分块矩阵,每个CPU所属的分块矩阵大小为8×2或2×8。
B)棋盘划分
就是将矩阵分成若干个子矩阵,每个子矩阵指派给一个处理器,此时任一处理器均不包含整行或者整列。
下图所示即为4个处理器情况下8×8矩阵的棋盘划分,其中处理器阵列为2×2,每个处理器分配到的子矩阵大小为4×4。
矩阵划分成棋盘状可以和处理器连成二维网孔相对应。
对于一个n×n维矩阵和p×p的二维处理器阵列,每个处理器均匀分配有(n/p)×(n/p)=n^2/p^2个元素。
使用棋盘式划分的矩阵相乘算法一般有两种,Cannon算法和Summa算法。
SUMMA算法能够计算m*l的A矩阵和l*n的B矩阵相乘(m、l、n可不相等),而cannon算法只能实现n*n的A矩阵和n*n的B矩阵相乘,具有很大的局限性。
3.2、算法原理
A)行划分法
假设是M*N,计算前,将矩阵N发送给所有从进程,然后将矩阵M分块,将M中数据按行分给各从进程,在从进程中计算M中部分行数据和N的乘积,最后将结果发送给主进程。
这里为了方便,有多少进程,就将M分了多少块,除最后一块外的其他数据块大小都相等,最后一块是剩下的数据,大小大于等于其他数据块大小,因为矩阵行数不一定整除进程数。
最后一块数据在主进程中计算,其他的在从进程中计算。
定义两个矩阵M和N,N所有进程都需要,M可以只在主进程中定义。
其他的变量视主进程和从进程需要按要求定义在合适的位置。
代码参见附录部分。
B)Cannon算法
Cannon算法的基本思想可以如下表示:
假设两个矩阵A和B相乘,把A和B矩阵划分成p个方块,进程的编号从
到
,并在最初把子矩阵
和
分配给
。
虽然第i行的每个进程需要全部
的个子矩阵
,但我们还是能调度第i行
个进程的计算,使得每个进程在任何时刻都是用不同的
。
每完成一次矩阵乘法,这些块在各进程之间被轮流使用,似的每次轮流之后每个进程都可以得到新的
。
对列使用同样的调度,则在任何时刻,任何进程至多拥有每个矩阵的一个块,在所有进程中,改算法需要的总内存量为
。
下图为此算法中不同进程上子矩阵乘法的调度过程。
假如矩阵C=A*B,则C的的计算公式如下:
进程P存储分块矩阵这一部分。
块矩阵乘法要计算所有匹配的
和
,然而只有在主对角线的才是匹配的。
因此需要采用循环移动分块矩阵的方法来使每个进程
都有一对可以直接相乘的匹配的块,具体方法如下:
(1)将排第i行的
块循环左移i个位置,将第
列.块循环上移j个位置;
(2)进程
执行乘一加运算,然后将移动得到的
块循环左移1个位置,将移动得到的
块循环上移1个位置;
(3)重复第2步(
一1)次,每次移动后
进程执行乘一加运算。
经过以上操作后就可以得到矩阵C的解。
代码请参见附录部分
C)Summa算法
SUMMA算法首先将A,B和C划分为相同大小的矩阵,对应放在mesh_r×mesh_c的二维mesh上。
但SUMMA算法将矩阵乘法分解为一系列的秩nb修正,即各处理器中的A和B分别被分解为nb大小的列块和行块进行相乘,前面所说的分块尺寸就是指nb的大小。
算法中,广播实现为逻辑处理器行环或列环上的流水线传送,达到了计算与通信的重叠.具体描述如算法1所示。
C=0
fori=0tok-1stepnbdo
curcol=i×c/n
currow=i×r/m
ifmycol=currol向本行广播A从imod(k/c)列开始的nb列,存于A′
ifmyrow=currow向本列广播B从imod(k/r)行开始的nb行,存于B′
C=A′×B′
endfor
SUMMA算法的核心思想是:
各处理器收集来自同一行处理器中A矩阵子块的所有列和同一列处理器中B矩阵子块的所有行,然后将行列相乘后累加,形成一个C矩阵的分块矩阵。
最后由rank=0的处理器将其他处理器的数据收集起来,形成最终的矩阵C。
SUMMA算法相较于cannon算法的优势只要体现在SUMMA算法能够计算m*l的A矩阵和l*n的B矩阵相乘(m、l、n可不相等),而cannon算法只能实现n*n的A矩阵和n*n的B矩阵相乘,具有很大的局限性。
代码参见附录部分。
3.3、程序运行结果对比分析
A)统一的实验条件
矩阵大小:
1000*1000;
矩阵数字范围:
0~10;
矩阵数字分布是否随机:
是;
分配的进程数:
9;
B)实验进程数解释
由于Cannon算法本身局限性,要使用Cannon算法,必须满足进程数为整数的平方,比如1、4、9、16等。
在本次的实验环境之下,经过多次对比分析,发现对于分行还是分块算法,进程数安排在8~15可以得到最理想的运行速度:
进程数目过小则每个进程单独运算的时间过多,进程数过大则选路时间(进程与进程之间的通信时间)过长。
而对比要求每个算法的进程相同,故此处选择进程数目为9.
C)算法运行时间对比
Cannon算法运行时间如下:
分行法运行时间如下:
串行算法运行时间如下:
由于Summa算法与Cannon算法思路几乎相同,而且在算法预处理阶段要比Cannon算法更耗时,故没有做多余的实验。
算法
分行
Cannon
串行
时间
1.218810s
0.76s
10.420s
显而易见,单纯的运用分行算法所花费的时间是最短的。
D)结果分析
Cannon算法相对于简单的行划分并行处理算法,其优势仅仅在于并行度可以更高(可达到N*N个进程,N为矩阵宽),但在并行度相同的情况下,其多出的预处理过程、矩阵发送与结果回收机制会占用更多的时间。
3.4、程序调优
A)行划分算法优化
1.循环优化
在预估计矩阵大小为10的倍数的基础上,对每一个步长为1的循环做处理,改为步长为10的循环,将十次循环体全部压缩在一次循环中,从而大量减少了循环的判别时间,提升循环运算速度。
例如在单个线程在计算部分结果时,采用的循环为:
for(i=0;i for(j=0;j DATAtemp=0; for(k=0;k temp+=buffer[i*width+k]*n[j*width+k]; temp+=buffer[i*width+k+1]*n[j*width+k+1]; temp+=buffer[i*width+k+2]*n[j*width+k+2]; temp+=buffer[i*width+k+3]*n[j*width+k+3]; temp+=buffer[i*width+k+4]*n[j*width+k+4]; temp+=buffer[i*width+k+5]*n[j*width+k+5]; temp+=buffer[i*width+k+6]*n[j*width+k+6]; temp+=buffer[i*width+k+7]*n[j*width+k+7]; temp+=buffer[i*width+k+8]*n[j*width+k+8]; temp+=buffer[i*width+k+9]*n[j*width+k+9]; } ans[i*width+j]=temp; } } 在将循环次数压缩的同时,为了进一步减少循环的运算量,在每一个步长为10的循环之前做预处理,避免循环体中的重复运算。 例如在主进程在接受其他进程时,将结果矩阵整合的过程: for(k=1;k { MPI_Recv(ans,line*width,MPI_INT,k,2,MPI_COMM_WORLD,&status); for(i=0;i { count=i*k*width;//将i*k*width提前算好,减少了下一步循环的重复运算 count1=i*width; for(j=0;j p[count+j]=ans[count1+j]; p[count+j+1]=ans[count1+j+1]; p[count+j+2]=ans[count1+j+2]; p[count+j+3]=ans[count1+j+3]; p[count+j+4]=ans[count1+j+4]; p[count+j+5]=ans[count1+j+5]; p[count+j+6]=ans[count1+j+6]; p[count+j+7]=ans[count1+j+7]; p[count+j+8]=ans[count1+j+8]; p[count+j+9]=ans[count1+j+9]; } } } 2.节省空间 在进行矩阵工作量划分并传送的时候,为每一个进程开辟仅仅是自己所需要大小的空间,例如在9进程的环境下,每个进程所需要接受的缓存空间为B矩阵大小以及大约1/9大小A矩阵。 内存开辟: buffer=(DATA*)malloc(sizeof(DATA)*width*line); 矩阵A分块传输: for(k=1;k { for(i=k;i { count=i/numprocs*width; count1=i*width; for(j=0;j { buffer[count+j]=m[count1+j]; buffer[count+j+1]=m[count1+j+1]; buffer[count+j+2]=m[count1+j+2]; buffer[count+j+3]=m[count1+j+3]; buffer[count+j+4]=m[count1+j+4]; buffer[count+j+5]=m[count1+j+5]; buffer[count+j+6]=m[count1+j+6]; buffer[count+j+7]=m[count1+j+7]; buffer[count+j+8]=m[count1+j+8]; buffer[count+j+9]=m[count1+j+9]; } } MPI_Send(buffer,line*width,MPI_INT,k,1,MPI_COMM_WORLD); 同样的方式也运用在运行空间的开辟上。 这样做不仅仅是内存空间的节约,同时也减少了进程之间的数据传输量,大大节省了进程之间的协作时间! B)稀疏矩阵处理 虽然程序并未对稀疏矩阵进行优化,但是还是试着对程序的输入数据模式进行更改,体验一下稀疏矩阵运算的速度提升有多快。 voidmagten(DATA*a,intwidth) { inti,j; srand((unsigned)time(NULL)); for(i=0;i for(j=0;j a[i*width+j]=(rand()%MUL)*(rand()%2)*(rand()%3); a[i*width+j+1]=(rand()%MUL)*(rand()%2)*(rand()%3); a[i*width+j+2]=(rand()%MUL)*(rand()%2)*(rand()%3); a[i*width+j+3]=(rand()%MUL)*(rand()%2)*(rand()%3); a[i*width+j+4]=(rand()%MUL)*(rand()%2)*(rand()%3); a[i*width+j+5]=(rand()%MUL)*(rand()%2)*(rand()%3); a[i*width+j+6]=(rand()%MUL)*(rand()%2)*(rand()%3); a[i*width+j+7]=(rand()%MUL)*(rand()%2)*(rand()%3); a[i*width+j+8]=(rand()%MUL)*(rand()%2)*(rand()%3); a[i*width+j+9]=(rand()%MUL)*(rand()%2)*(rand()%3); } } } 如上面所示,对于每一个生成的数据都再一次进行乘法运算,其中乘数是0的概率为2/3. 运行结果如下: 可以看出,稀疏矩阵的乘法由0.76s变为0.75s,仅仅是短暂的提升。 提升计时显示的精度,可以看到,对于稀疏矩阵的处理要比普通矩阵快0.015s,提速约为2% 3.5、算法最优速度 对算法运用不同的进程数目运算进行了大量重复试验,最终得出在进程数大概为12的时候,本算法的运行速度最快。 最终结果如下: 发出工作量时间为0.138184s 运算时间为0.495569s 接收答案时间为0.025461s 总运算时间0.659240s 四、总结分析 串行\并行 行划分 串行 运行时间 0.659240s 10.420s 进程数目 12 1 加速比 15.80 1 效率 15.80/12=1.317 1 从效率大于1上可以看出,本次课程设计做出的算法为超线性加速,这主要得益于对循环体的优化。 附录 Cannon算法 #include #include"mpi.h" #include #include #include #include #defineMUL10 MPI_Statusstatus; double**A,**B,**C;//C=A*B double*a,*b,*c;//各个进程的缓冲区 intn;//矩阵的行列数 intnp;//每个进程控制的小矩阵的行列数 intp,rank;//进程个个数、当前进程的编号,笛卡尔进程编号 double*tempa,*tempb; voidProduceABC();//在根处理器中生成矩阵AB,初始化矩阵C voidPrintABC();//输出结果 voidScatterAB();//分发矩阵AB中的元素到各个进程中 voidMainProcess();//cannon算法的主过程 voidcollectC();//收集结果矩阵C voidMutiply();//矩阵相乘 voidPrintab(); voidPrintc(); intmain(intargc,char*argv[]) { inti; doublestarttime,endtime; MPI_Init(&argc,&argv); MPI_Comm_size(MPI_COMM_WORLD,&p); MPI_Comm_rank(MPI_COMM_WORLD,&rank); if(rank==0) { printf("pleaseinputtherawnumber: n="); fflush(stdout); scanf("%d",&n); printf("\n"); } MPI_Bcast(&n,1,MPI_DOUBLE,0,MPI_COMM_WORLD); //n=atoi(argv[1]); np=n/(int)sqrt(p); a=(double*)malloc(np*np*sizeof(double)); b=(double*)malloc(np*np*sizeof(double)); c=(double*)malloc(np*np*sizeof(double)); memset(c,0,np*np*sizeof(double)); tempa=(double*)malloc(np*np*sizeof(double)); tempb=(double*)malloc(np*np*sizeof(double)); if(rank==0) { //在根处理器中为矩阵ABC分配空间 A=(double**)malloc(n*sizeof(double*)); B=(double**)malloc(n*sizeof(double*)); C=(double**)malloc(n*sizeof(double*)); for(i=0;i { A[i]=(double*)malloc(n*sizeof(double)); B[i]=(double*)malloc(n*sizeof(double)); C[i]=(double*)malloc(n*sizeof(double)); } ProduceABC();//在根处理器中随机生成矩阵AB,初始化矩阵C ScatterAB();//分发矩阵AB中的元素到各个进程中 } else { MPI_Recv(a,np*np,MPI_DOUBLE,0,1,MPI_COMM_WORLD,&status); MPI_Recv(b,np*np,MPI_DOUBLE,0,2,MPI_COMM_WORLD,&status); } starttime=MPI_Wtime();//开始时间 MainProcess();//cannon算法的主过程 if(rank==0) { collectC();//收集结果矩阵C PrintABC();//输出结果 endtime=MPI_Wtime(); printf("timeused: %lf\n",endtime-starttime); for(i=0;i { free(A[i]); free(B[i]); free(C[i]); } free(A); free(B); free(C); } else { MPI_Send(c,np*np,MPI_DOUBLE,0,1,MPI_COMM_WORLD); } free(a); free(b); free(c); free(tempa); free(tempb); MPI_Finalize(); return0; } voidProduceABC()//在根处理器中生成矩阵AB { inti,j; srand((unsigned)time(NULL)); for(i=0;i for(j=0;j { A[i][j]=rand()%MUL; B[i][j]=rand()%MUL; C[i][j]=0.0; } } voidPrintABC()//输出结果 { printf("A[0][0]=%f\nB[0][0]=%f\nC[0][0]=%f\n",A[0][0],B[0][0],C[0][0]); } voidScatterAB()//分发矩阵AB中的元素到各个进程中 { intimin,imax,jmin,jmax; intsp; inti,j,k,m; for(k=0;k { /*计算相应处理器所分得的矩阵块在总矩阵中的坐标范围*/ sp=(int)sqrt(p); imin=(k/sp)*np; imax=imin+np-1; jmin=(k%sp)*np; jmax=jmin+np-1; /*rank=0的处理器将A,B中的相应块拷至tempa,tempb,准备向其他处理器发送*/ m=0; for(i=imin;i<=imax;i++) { for(j=jmin;j<=jmax;j++) { tempa[m]=A[i][j]; tempb[m]=B[j][i];//矩阵B按列优先存储 m++; } } /*根处理器将自己对应的矩阵块从tempa,tempb拷至a,b*/ if(k==0) { memcpy(a,tempa,np*np*sizeof(double)); memcpy(b,tempb,np*np*sizeof(double)); } else/*根处理器向其他处理器发送tempa,tempb中相关的矩阵块*/ { MPI_Send(tempa,np*np,MPI_DOUBLE,k,1,MPI_COMM_WORLD); MPI_Send(tempb
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 矩阵 相乘 并行 算法