并行計算矩陣分塊乘法

1、目錄一、題目及要求11、題目12、要求1二、設(shè)計算法、算法原理1三、算法描述、設(shè)計流程23.1 算法描述23.2 設(shè)計流程4四、源程序代碼及運行結(jié)果61、超立方61.1 超立方的源程序代碼61.2 運行結(jié)果112、網(wǎng)孔連接112.1 源程序代碼112.2 運行結(jié)果183、在數(shù)學軟件中的計算結(jié)果19五、算法分析、優(yōu)缺點191、簡單的并行分塊乘法的過程為192、使用Cannon算法時的算法分析203、算法的優(yōu)缺點21六、總結(jié)22參考文獻23、題目及要求1、題目簡單并行分塊乘法：(1)情形1：超立方連接；(2)情形2:二維環(huán)繞網(wǎng)孔連接12334、,9-23-74)已知A=20121-143B=121

2、543,求C=AmB-271591013-5918-3-567>-11517172、要求(1)題目分析、查閱與題目相關(guān)的資料;(2)設(shè)計算法；(3)算法實現(xiàn)、代碼編寫；(4)結(jié)果分析和性能分析與改進；(5)論文撰寫、答辯；二、設(shè)計算法、算法原理要考慮的計算問題是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三、算法描述、設(shè)計流程3.1 算法描述超立方情形下矩陣的簡單并行分塊算法輸入：待選路的信包在源處理器中輸出：將原處理器中的信包送至其目的地Begin(1) fori=1tondoh=Sj-diendfor(2) i=1,V=S(3) whilei<ndo(3.(1) ifri=1then從當前節(jié)點V選路到節(jié)點為V1(3.(2) i=i+1endwhileEnd二維網(wǎng)孔情形下矩陣的簡單并行分塊算法輸入：待選路的信包處于源處理器中輸出：將各信包送至各自的目的地中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設(shè)計流程如圖3-1以下是二維網(wǎng)孔與超立方連接設(shè)計流程二維網(wǎng)孔步驟：(1)先進行行播送；(2)再同時進行列播送；144圖3-1二維網(wǎng)孔示意圖超立方步驟：依次從低維到高維播送,d-立方，d=0,1,2,3,4算法流程如圖所

5、示：圖3-2算法流程四、源程序代碼及運行結(jié)果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運行結(jié)果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個處理器的運行結(jié)果2、網(wǎng)孔連接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、區(qū)*/總矩陣,C=A*B*/、b、c表分塊，tmp_a、tmp_b表緩intdg,dl,dl2,p,sp;/*dg:總矩陣維數(shù);dl:矩陣塊維數(shù);dl2=dl*dl;p:處理器個數(shù);sp=sqrt(p)*/intmy_rank,my_row,my_col;/*my_rank:處理器ID;(my_row,my_col):處理器應(yīng)輯陣列坐溫*/一一一一11MPI_Statusstatus;/* 函數(shù)名：get_index* 功能：處理外邏輯陣列坐標至rank號的轉(zhuǎn)換* 輸入：坐標、邏輯陣列維數(shù)* 輸出：rank號* /intget_index(introw,intcol,intsp)(return

22、(row+sp)%sp)*sp+(col+sp)%sp;/*函數(shù)名：random_A_B*功能：隨機生成矩陣A和B*/voidrandom_A_B()(inti,j;floatm;/srand(unsignedint)time(NULL);/*設(shè)隨機數(shù)種子*/*隨機生成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;/*函數(shù)名：scatter_A_B12功能：rank為0的處理器向其他處理器發(fā)送A、B矩陣的相關(guān)塊*/voidscatter_A_B()(inti,j,k,l;intp_imin,p_imax,p_jmin,p_jmax;for(k=0;k<p;k+)(/*計算相應(yīng)處理器所分得的矩陣塊在總矩陣中的坐標范圍*/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中的相應(yīng)塊拷至tmp_a,tmp_b,準備向其

24、他處理器發(fā)送*/一一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的處理器直接將自己對應(yīng)的矩陣塊從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的處理器向其他處理器發(fā)送tmp_a,tmp_b中相關(guān)的矩陣塊*/(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* 函數(shù)名: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)向上循環(huán)移動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);./* 函數(shù)名：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);/*函數(shù)名：collect_c*功能：rank為0的處理器從其余處理器收集分塊矩陣c*/voidcollect_C()(.inti,j,i2,j2,k;intp_imin,p_imax,p_jmin,p_jmax;/*分塊矩陣在總矩陣中頂點邊界值*/*將rank為0的處理器中分塊矩陣c結(jié)果賦給總矩

29、陣C對應(yīng)位置*/for(i=0;i<dl;i+)for(j=0;j<dl;j+)Cij=ci*dl+j;for(k=1;k<p;k+)(/*將rank為0的處理器從其他處理器接收相應(yīng)的分塊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中的相應(yīng)位置，從而構(gòu)造出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/*函數(shù)名：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、數(shù)名：main*功能：主過程，Cannon算法，矩陣相乘*輸入：argc為命令行參數(shù)個數(shù)，argv為每個命令行參數(shù)組成的字符串數(shù)組*/intmain(intargc,char*argv)(inti;MPI_Init(&argc,&argv);/*啟動MPI計算*/MPI_Comm_size(MPI_COMM_WORLD,&p);/*確定處理器個數(shù)*/MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);/*確定各自的處理器標識符*/sp=sqrt(p);/*確保處理器個數(shù)是完全平方數(shù)，否則打印錯誤信息，程序退出*/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);/*總矩陣維數(shù)*/dl=dg/sp;/*計算分塊矩陣維數(shù)*/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的處理器向其他處理器發(fā)送A、B矩陣的相關(guān)塊*/else/*rank不為0的處理器接收來自rank為0的處理器的相應(yīng)矩陣分塊*/(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();/*結(jié)束MPI計算*/return0;2.2運行結(jié)果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個處理器的運行結(jié)果183、在數(shù)學軟件中的計算結(jié)果»尾口2,現(xiàn)£20121,-49;比75優(yōu)門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中的運行結(jié)果五、算法分析、優(yōu)缺點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,總結(jié)需要所有塊A,j與塊B,j(0<k<7-1)，為19此在每行處理器之間要進行A矩陣塊的多到多播送，已得到A,k；而在每列的處理器之間要進行B矩陣塊的多到多播送，已得到Bk,j。(3