并行计算矩阵分块乘法

1、目录一、题目及要求11、题目12、要求1二、设计算法、算法原理1三、算法描述、设计流程23.1 算法描述23.2 设计流程4四、源程序代码及运行结果61、超立方61.1 超立方的源程序代码61.2 运行结果112、网孔连接112.1 源程序代码112.2 运行结果183、在数学软件中的计算结果19五、算法分析、优缺点191、简单的并行分块乘法的过程为192、使用Cannon算法时的算法分析203、算法的优缺点21六、总结22参考文献23、题目及要求1、题目简单并行分块乘法：(1)情形1：超立方连接；(2)情形2:二维环绕网孔连接12334、,9-23-74)已知A=20121-143B=121

2、543,求C=AmB-271591013-5918-3-567>-11517172、要求(1)题目分析、查阅与题目相关的资料;(2)设计算法；(3)算法实现、代码编写；(4)结果分析和性能分析与改进；(5)论文撰写、答辩；二、设计算法、算法原理要考虑的计算问题是C=AB其中A与B分别是口父n矩阵A、B和C分成p=Jp父Jp的方块阵Aj,Bo和Cj,大小均为处理器编号为p。,。,，.£5Pj存放Aj，Bij和Cij通讯：每行处理器进行A矩阵块的多到多播送(得到Ak,k=。7p-1)每列处理器进行B矩阵块的多到多播送(得至UBkj,k=。1)p1乘-加运算：Pj做Cj=£

3、AkBkjk=e三、算法描述、设计流程3.1 算法描述超立方情形下矩阵的简单并行分块算法输入：待选路的信包在源处理器中输出：将原处理器中的信包送至其目的地Begin(1) fori=1tondoh=Sj-diendfor(2) i=1,V=S(3) whilei<ndo(3.(1) ifri=1then从当前节点V选路到节点为V1(3.(2) i=i+1endwhileEnd二维网孔情形下矩阵的简单并行分块算法输入：待选路的信包处于源处理器中输出：将各信包送至各自的目的地中Begin(1)沿x维将信包向左或向右选路至目的地的处理器所在的列(2)沿y维将信包向上或向下选路至目的地的处理器所

4、在的行分块乘法算法/输入：入於BnM;子快大小均为输出：CnnnBegin(1)fori=0top-1doforallpar-dopjifi>kthenAjA,ijmocendififj>kthenBijB(i+1)mod,jendifendforendforfori=0to.p-1doforallpjpar-doCij=Aj+BjendforEndforEnd3.2设计流程如图3-1以下是二维网孔与超立方连接设计流程二维网孔步骤：(1)先进行行播送；(2)再同时进行列播送；144图3-1二维网孔示意图超立方步骤：依次从低维到高维播送,d-立方，d=0,1,2,3,4算法流程如图所

5、示：图3-2算法流程四、源程序代码及运行结果1、超立方1.1超立方的源程序代码#include"stdio.h"#include"stdlib.h"#include"mpi.h"#defineintsizesizeof(int)#definefloatsizesizeof(float)#definecharsizesizeof(char)#defineA(x,y)Ax*K+y#defineB(x,y)Bx*N+y#defineC(x,y)Cx*N+y#definea(x,y)ax*K+y#defineb(x,y)bx*n+y#defi

6、nebuffer(x,y)bufferx*n+y#definec(l,x,y)cx*N+y+l*nfloat*a,*b,*c,*buffer;ints;float*A,*B,*C;intM,N,K,P;intm,n;intmyid;intp;FILE*dataFile;MPI_Statusstatus;doubletime1;doublestarttime,endtime;voidreadData()inti,j;starttime=MPI_Wtime();dataFile=fopen("yin.txt","r");fscanf(dataFile,&qu

7、ot;%d%d",&M,&K);A=(float*)malloc(floatsize*M*K);for(i=0;i<M;i+)for(j=0;j<K;j+)fscanf(dataFile,"%f",A+i*K+j);6)fscanf(dataFile,"%d%d",&P,&N);if(K!=P)(printf("theinputiswrongn");exit；)B=(float*)malloc(floatsize*K*N);for(i=0;i<K;i+)(for(j=0;j&

8、lt;N;j+)(fscanf(dataFile,"%f",B+i*N+j);)fclose(dataFile);printf("Inputoffile"yin.txt"n");printf("%dt%dn",M,K);for(i=0;i<M;i+)(for(j=0;j<K;j+)printf("%ft",A(i,j);printf("n");)printf("%dt%dn",K,N);for(i=0;i<K;i+)(for(j=0;j&

9、lt;N;j+)printf("%ft",B(i,j);printf("n");)C=(float*)malloc(floatsize*M*N);)intgcd(intM,intN,intgroup_size)(.inti;for(i=M;i>0;i-)(if(M%i=0)&&(N%i=0)&&(i<=group_size)returni;)return1;)voidprintResult()(inti,j;printf("nOutputofMatrixC=ABn");for(i=0;i&l

10、t;M;i+)(for(j=0;j<N;j+)printf("%ft",C(i,j);printf("n");endtime=MPI_Wtime();printf("n");printf("Wholerunningtime=%fsecondsn",endtime-starttime);printf("Distributedatatime=%fsecondsn",time1-starttime);printf("Parallelcomputetime=%fsecondsn"

11、;,endtime-time1);intmain(intargc,char*argv)(inti,j,k,l,group_size,mp1,mm1;MPI_Init(&argc,&argv);MPI_Comm_size(MPI_COMM_WORLD,&group_size);MPI_Comm_rank(MPI_COMM_WORLD,&myid);p=group_size;if(myid=0)(readData();if(myid=0)for(i=1;i<p;i+)(MPI_Send(&M,1,MPI_INT,i,i,MPI_COMM_WORLD);

12、MPI_Send(&K,1,MPI_INT,i,i,MPI_COMM_WORLD);MPI_Send(&N,1,MPI_INT,i,i,MPI_COMM_WORLD);else(MPI_Recv(&M,1,MPI_INT,0,myid,MPI_COMM_WORLD,&status);MPI_Recv(&K,1,MPI_INT,0,myid,MPI_COMM_WORLD,&status);MPI_Recv(&N,1,MPI_INT,0,myid,MPI_COMM_WORLD,&status);p=gcd(M,N,group_size

13、);m=M/p;n=N/p;if(myid<p)(a=(float*)malloc(floatsize*m*K);b=(float*)malloc(floatsize*K*n);c=(float*)malloc(floatsize*m*N);8if(myid%2!=0)buffer=(float*)malloc(K*n*floatsize);if(a=NULL|b=NULL|c=NULL)printf("Allocatespacefora,borcfail!");if(myid=0)for(i=0;i<m;i+)for(j=0;j<K;j+)a(i,j)=

14、A(i,j);for(i=0;i<K;i+)for(j=0;j<n;j+)b(i,j)=B(i,j);if(myid=0)for(i=1;i<p;i+)MPI_Send(&A(m*i,0),K*m,MPI_FLOAT,i,i,MPI_COMM_WORLD);for(j=0;j<K;j+)MPI_Send(&B(j,n*i),n,MPI_FLOAT,i,i,MPI_COMM_WORLD);free(A);free(B);elseMPI_Recv(a,K*m,MPI_FLOAT,0,myid,MPI_COMM_WORLD,&status);for(j

15、=0;j<K;j+)MPI_Recv(&b(j,0),n,MPI_FLOAT,0,myid,MPI_COMM_WORLD,&status);if(myid=0)time1=MPI_Wtime();for(i=0;i<p;i+)l=(i+myid)%p;for(k=0;k<m;k+)for(j=0;j<n;j+)for(c(l,k,j)=0,s=0;s<K;s+)c(l,k,j)+=a(k,s)*b(s,j);mm1=(p+myid-1)%p;mp1=(myid+1)%p;if(i!=p-1)if(myid%2=0)(MPI_Send(b,K*n,M

16、PI_FLOAT,mm1,mm1,MPI_COMM_WORLD);MPI_Recv(b,K*n,MPI_FLOAT,mp1,myid,MPI_COMM_WORLD,&status);else(for(k=0;k<K;k+)for(j=0;j<n;j+)buffer(k,j)=b(k,j);MPI_Recv(b,K*n,MPI_FLOAT,mp1,myid,MPI_COMM_WORLD,&status);MPI_Send(buffer,K*n,MPI_FLOAT,mm1,mm1,MPI_COMM_WORLD);if(myid=0)for(i=0;i<m;i+)f

17、or(j=0;j<N;j+)C(i,j)=*(c+i*N+j);if(myid!=0)MPI_Send(c,m*N,MPI_FLOAT,0,myid,MPI_COMM_WORLD);else(for(k=1;k<p;k+)(MPI_Recv(c,m*N,MPI_FLOAT,k,k,MPI_COMM_WORLD,&status);for(i=0;i<m;i+)for(j=0;j<N;j+)C(k*m+i),j)=*(c+i*N+j);if(myid=0)printResult();MPI_Finalize();if(myid<p)(free(a);free(

18、b);free(c);10if(myid=0)free(C);if(myid%2!=0)free(buffer);)return(0);)1.2运行结果Inputoffile''yin.tjct"441.0000002.000000201.000000-2.000000e.ooooao449.00000012,000000101.000000'ii.ooaooo21.00003071.000000-3.000000-2.0000001,0000003,00000051,00000033,K箕如-14.0000005,ODOaQO-SG.oaaooo3,oooa

19、oc54.000000-S.M州如7.0000004.0000003,0000009.ooooao7.000000-74.0000003.0000009.0000001".OO0OOODuuputofMatrixC=Afi3322.000000303-00000614.000000_£4二，”二二-5697,000000-270.000000549-OQOaOC1.70.000000-26.OOOCOO1828.aOCOOO191.000000297.000000-145£.OCOOOO55次WMOO-966,000000Wholerunningtiine=0r0

20、01000j&ccndaDistributedatatiire三d.DQloaosecondsparalleleciiDetime-口seccnds图4.14个处理器的运行结果2、网孔连接2.1源程序代码#include<stdlib.h>#include<string.h>#include<mpi.h>#include<time.h>#include<stdio.h>#include<math.h>/*全局变量声明*/float*A,*B,*C;/*float*a,*b,*c,*tmp_a,*tmp_b;/*a冲

21、区*/总矩阵,C=A*B*/、b、c表分块，tmp_a、tmp_b表缓intdg,dl,dl2,p,sp;/*dg:总矩阵维数;dl:矩阵块维数;dl2=dl*dl;p:处理器个数;sp=sqrt(p)*/intmy_rank,my_row,my_col;/*my_rank:处理器ID;(my_row,my_col):处理器应辑阵列坐温*/一一一一11MPI_Statusstatus;/* 函数名：get_index* 功能：处理外逻辑阵列坐标至rank号的转换* 输入：坐标、逻辑阵列维数* 输出：rank号* /intget_index(introw,intcol,intsp)(return

22、(row+sp)%sp)*sp+(col+sp)%sp;/*函数名：random_A_B*功能：随机生成矩阵A和B*/voidrandom_A_B()(inti,j;floatm;/srand(unsignedint)time(NULL);/*设随机数种子*/*随机生成A,B,并初始化C*/for(i=0;i<dg;i+)for(j=0;j<dg;j+)(scanf("%f",&m);Aij=m;Cij=0.0;m=0;for(i=0;i<dg;i+)for(j=0;j<dg;j+)(scanf("%f",&m);

23、Bij=m;m=0;/*函数名：scatter_A_B12功能：rank为0的处理器向其他处理器发送A、B矩阵的相关块*/voidscatter_A_B()(inti,j,k,l;intp_imin,p_imax,p_jmin,p_jmax;for(k=0;k<p;k+)(/*计算相应处理器所分得的矩阵块在总矩阵中的坐标范围*/p_jmin=(k%sp)*dl;p_jmax=(k%sp+1)*dl-1;p_imin=(k-(k%sp)/sp*dl;p_imax=(k-(k%sp)/sp+1)*dl-1;l=0;/*rank=0的处理器将A,B中的相应块拷至tmp_a,tmp_b,准备向其

24、他处理器发送*/一一for(i=p_imin;i<=p_imax;i+)(一一for(j=p_jmin;j<=p_jmax;j+)(tmp_al=Aij;tmp_bl=Bij;l+;)/*rank=0的处理器直接将自己对应的矩阵块从tmp_a,tmp_b拷至a,b*/if(k=0)(memcpy(a,tmp_a,dl2*sizeof(float);memcpy(b,tmp_b,dl2*sizeof(float);else/*rank=0的处理器向其他处理器发送tmp_a,tmp_b中相关的矩阵块*/(MPI_Send(tmp_a,dl2,MPI_FLOAT,k,1,MPI_COMM

25、_WORLD);MPI_Send(tmp_b,dl2,MPI_FLOAT,k,2,MPI_COMM_WORLD);/*13* 函数名:init_alignment* 功能：矩阵仄和B初始对准* /voidinit_alignment()(.MPI_Sendrecv(a,dl2,MPI_FLOAT,get_index(my_row,my_col-my_row,sp),1,tmp_a,dl2,MPI_FLOAT,get_index(my_row,my_col+my_row,sp),1,MPI_COMM_WORLD,&status);memcpy(a,tmp_a,dl2*sizeof(flo

26、at);/*将B中坐标为(i,j)的分块B(i,j)向上循环移动j步*/MPI_Sendrecv(b,dl2,MPI_FLOAT,get_index(my_row-my_col,my_col,sp),1,tmp_b,dl2,MPI_FLOAT,get_index(my_row+my_col,my_col,sp),1,MPI_COMM_WORLD,&status);memcpy(b,tmp_b,dl2*sizeof(float);./* 函数名：main_shift* 功能：分块矩M左移和上移，并计算分块c* /voidmain_shift()(.inti,j,k,l;for(l=0;l

27、<sp;l+)(/*矩阵块相乘，c+=a*b*/for(i=0;i<dl;i+)for(j=0;j<dl;j+)for(k=0;k<dl;k+)ci*dl+j+=ai*dl+k*bk*dl+j;/*将分块a左移1位*/MPI_Send(a,dl2,MPI_FLOAT,get_index(my_row,my_col-1,sp),1,MPI_COMM_WORLD);MPI_Recv(a,dl2,MPI_FLOAT,get_index(my_row,my_col+1,sp),1,MPI_COMM_WORLD,&status);/*将分块b上移1位*/MPI_Send(

28、b,dl2,MPI_FLOAT,get_index(my_row-1,my_col,sp),1,MPI_COMM_WORLD);14MPI_Recv(b,dl2,MPI_FLOAT,get_index(my_row+1,my_col,sp),1,MPI_COMM_WORLD,&status);/*函数名：collect_c*功能：rank为0的处理器从其余处理器收集分块矩阵c*/voidcollect_C()(.inti,j,i2,j2,k;intp_imin,p_imax,p_jmin,p_jmax;/*分块矩阵在总矩阵中顶点边界值*/*将rank为0的处理器中分块矩阵c结果赋给总矩

29、阵C对应位置*/for(i=0;i<dl;i+)for(j=0;j<dl;j+)Cij=ci*dl+j;for(k=1;k<p;k+)(/*将rank为0的处理器从其他处理器接收相应的分块c*/MPI_Recv(c,dl2,MPI_FLOAT,k,1,MPI_COMM_WORLD,&status);p_jmin=(k%sp)*dl;p_jmax=(k%sp+1)*dl-1;p_imin=(k-(k%sp)/sp*dl;p_imax=(k-(k%sp)/sp+1)*dl-1;i2=0;/*将接收到的c拷至C中的相应位置，从而构造出C*/for(i=p_imin;i<

30、;=p_imax;i+)(一一j2=0;for(j=p_jmin;j<=p_jmax;j+)(Cij=ci2*dl+j2;j2+;i2+;15/*函数名：print*功能：打印矩阵*输入：指向矩阵指针的指针，字符串*/voidprint(float*m,char*str)(inti,j;printf("%s",str);/*打印矩阵m*/for(i=0;i<dg;i+)(for(j=0;j<dg;j+)printf("%15.0f",mij);printf("n");printf("n");/*函

31、数名：main*功能：主过程，Cannon算法，矩阵相乘*输入：argc为命令行参数个数，argv为每个命令行参数组成的字符串数组*/intmain(intargc,char*argv)(inti;MPI_Init(&argc,&argv);/*启动MPI计算*/MPI_Comm_size(MPI_COMM_WORLD,&p);/*确定处理器个数*/MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);/*确定各自的处理器标识符*/sp=sqrt(p);/*确保处理器个数是完全平方数，否则打印错误信息，程序退出*/if(sp*sp!=p)

32、(if(my_rank=0)printf("Numberofprocessorsisnotaquadraticnumber!n");MPI_Finalize();exit(1);if(argc!=2)(16mpirun-npProcNumcannonMatrixDimensionn");if(my_rank=0)printf("usage:MPI_Finalize();exit(1);dg=atoi(argv1);/*总矩阵维数*/dl=dg/sp;/*计算分块矩阵维数*/dl2=dl*dl;/*计算处理器在逻辑阵列中的坐标*/my_col=my_ran

33、k%sp;my_row=(my_rank-my_col)/sp;/*为a、b、c分配空间*/a=(float*)malloc(dl2*sizeof(float);b=(float*)malloc(dl2*sizeof(float);c=(float*)malloc(dl2*sizeof(float);/*初始化c*/for(i=0;i<dl2;i+)ci=0.0;/*为tmp_a、tmp_b分配空间*/tmp_a=(float*)malloc(dl2*sizeof(float);tmp_b=(float*)malloc(dl2*sizeof(float);if(my_rank=0)(./

34、*rank为0的处理器为A、RC分配空间*/A=(float*)malloc(dg*sizeof(float*);B=(float*)malloc(dg*sizeof(float*);C=(float*)malloc(dg*sizeof(float*);for(i=0;i<dg;i+)(Ai=(float*)malloc(dg*sizeof(float);Bi=(float*)malloc(dg*sizeof(float);Ci=(float*)malloc(dg*sizeof(float);random_A_B();/*rank为0的处理器随机化生成AB矩阵*/scatter_A_B(

35、);/*rank为0的处理器向其他处理器发送A、B矩阵的相关块*/else/*rank不为0的处理器接收来自rank为0的处理器的相应矩阵分块*/(MPI_Recv(a,dl2,MPI_FLOAT,0,1,MPI_COMM_WORLD,&status);MPI_Recv(b,dl2,MPI_FLOAT,0,2,MPI_COMM_WORLD,&status);17init_alignment();/*A、B矩阵的初始对准*/main_shift();/*分块矩阵左移、上移，cannon算法的主过程*/if(my_rank=0)为0的处理器从其余处理器收集分块矩阵(.collect

36、_C();/*rank*/打印矩阵A*/打印矩阵B*/打印矩阵C*/print(A,"randommatrixA:n");/*print(B,"randommatrixB:n");/*print(C,"MatrixC=A*B:n");/*else(MPI_Send(c,dl2,MPI_FLOAT,0,1,MPI_COMM_WORLD);MPI_Barrier(MPI_COMM_WORLD);/*同步所有处理器*/MPI_Finalize();/*结束MPI计算*/return0;2.2运行结果1233420121-143-2S9S-3

37、-5675-23-%121SJ3LOL3-59-1151717randomjnaurixA:1233420121-3.43-27159e-3-567randcmmatrixB:g-231215431013Sq-1151717MatrixC=A1B;3322303-2297614-2701823-148S61240S49焚品559-5S97170151-956图4.24个处理器的运行结果183、在数学软件中的计算结果»尾口2,现£20121,-49;比75优门A二I233420121-143-2715S8-3-567»0=-233,-F4;12a1,54,3;101

38、,3,-5,3,-11,51,7,17B=9-23-7412154习1013-5e-1151717>>C=A*Bc二3322303-26297614-2701828-14886124口5493866559-5697170101-936图4.3在MATLAB中的运行结果五、算法分析、优缺点1、简单的并行分块乘法的过程为(1)分块：将：Anxn与Bnxn分成P块A,j和B,j(0&i,j&JP-1),每块大小为(n/Jp)M(n/jp),并将它们分配给JpMjp个处理器(B,。，.，P0,«”.，Pqj,=)。开始时凡存放在块A,j与块B,j,然后计算块C,j。(2)通信：为了计算块C,j,总结需要所有块A,j与块B,j(0<k<7-1)，为19此在每行处理器之间要进行A矩阵块的多到多播送，已得到A,k；而在每列的处理器之间要进行B矩阵块的多到多播送，已得到Bk,j。(3