并行計(jì)算矩陣分塊乘法

1、目錄一、題目及要求11、題目12、要求1二、設(shè)計(jì)算法、算法原理1三、算法描述、設(shè)計(jì)流程23.1算法描述23.2設(shè)計(jì)流程4四、源程序代碼及運(yùn)行結(jié)果61、超立方61.1超立方的源程序代碼61.2運(yùn)行結(jié)果112、網(wǎng)孔連接112.1源程序代碼112.2運(yùn)行結(jié)果183、在數(shù)學(xué)軟件中的計(jì)算結(jié)果19五、算法分析、優(yōu)缺點(diǎn)191、簡(jiǎn)單的并行分塊乘法的過(guò)程為192、使用Cannon算法時(shí)的算法分析203、算法的優(yōu)缺點(diǎn)21六、總結(jié)22參考文獻(xiàn)23一、題目及要求1、題目簡(jiǎn)單并行分塊乘法：（1）情形1: 超立方連接；（2）情形2：二維環(huán)繞網(wǎng)孔連接已知求。2、要求（1）題目分析、查閱與題目相關(guān)的資料；（2）設(shè)計(jì)算法；（3

2、）算法實(shí)現(xiàn)、代碼編寫；（4）結(jié)果分析和性能分析與改進(jìn)；（5）論文撰寫、答辯；二、設(shè)計(jì)算法、算法原理要考慮的計(jì)算問(wèn)題是C=AB,其中A與B分別是矩陣。A、B和C分成的方塊陣,和,大小均為 ，p個(gè)處理器編號(hào)為, 存放,和。通訊:每行處理器進(jìn)行A矩陣塊的多到多播送(得到, k=0 )每列處理器進(jìn)行B矩陣塊的多到多播送(得到, k=0 )乘-加運(yùn)算: 做三、算法描述、設(shè)計(jì)流程3.1算法描述超立方情形下矩陣的簡(jiǎn)單并行分塊算法輸入：待選路的信包在源處理器中輸出：將原處理器中的信包送至其目的地Begin（1） for i=1 to n doendfor (2) (3) while do (3.1)if th

3、en從當(dāng)前節(jié)點(diǎn)V選路到節(jié)點(diǎn)為V1 （3.2） endwhileEnd二維網(wǎng)孔情形下矩陣的簡(jiǎn)單并行分塊算法輸入：待選路的信包處于源處理器中輸出：將各信包送至各自的目的地中Begin（1） 沿x維將信包向左或向右選路至目的地的處理器所在的列（2） 沿y維將信包向上或向下選路至目的地的處理器所在的行分塊乘法算法 /輸入: , ; 子快大小均為 輸出: n Begin (1)for i=0 to do for all par-do if ik then endif if jk then B(i+1)mod , j endif endfor endfor for i=0 to do for all pa

4、r-do =+ endfor Endfor End3.2設(shè)計(jì)流程以下是二維網(wǎng)孔與超立方連接設(shè)計(jì)流程。 如圖3-1二維網(wǎng)孔步驟：(1)先進(jìn)行行播送;(2)再同時(shí)進(jìn)行列播送;37625149804444444442233331111213141510 圖3-1 二維網(wǎng)孔示意圖超立方步驟：依次從低維到高維播送, d-立方, d=0,1,2,3,4;算法流程如圖所示：包括mpi的頭文件相關(guān)變量聲明MPI_INIT()MPI_COMM_RANK()MPI_COMM_SIZE()進(jìn)入MPI系統(tǒng)矩陣內(nèi)部的通信應(yīng)用控制實(shí)體：矩陣內(nèi)部的計(jì)算程序MPI_FINALIZE()退出MPI系統(tǒng)結(jié)束開始循環(huán)直至結(jié)束圖3-

5、2 算法流程四、源程序代碼及運(yùn)行結(jié)果1、超立方1.1超立方的源程序代碼#include stdio.h#include stdlib.h#include mpi.h#define intsize sizeof(int)#define floatsize sizeof(float)#define charsize sizeof(char)#define A(x,y) Ax*K+y#define B(x,y) Bx*N+y#define C(x,y) Cx*N+y#define a(x,y) ax*K+y#define b(x,y) bx*n+y#define buffer(x,y) buffer

6、x*n+y #define c(l,x,y) cx*N+y+l*nfloat *a,*b,*c,*buffer;int s;float *A,*B,*C; int M,N,K,P ;int m,n;int myid;int p; FILE *dataFile; MPI_Status status;double time1;double starttime,endtime;void readData() int i,j; starttime = MPI_Wtime(); dataFile=fopen(yin.txt,r); fscanf(dataFile,%d%d, &M, &K); A=(fl

7、oat *)malloc(floatsize*M*K); for(i = 0; i M; i+) for(j = 0; j K; j+) fscanf(dataFile,%f, A+i*K+j); fscanf(dataFile,%d%d, &P, &N); if (K!=P) printf(the input is wrongn); exit(1); B=(float *)malloc(floatsize*K*N); for(i = 0; i K; i+) for(j = 0; j N; j+) fscanf(dataFile,%f, B+i*N+j); fclose(dataFile);

8、printf(Input of file yin.txtn); printf(%dt %dn,M, K); for(i=0;iM;i+) for(j=0;jK;j+) printf(%ft,A(i,j); printf(n); printf(%dt %dn,K, N); for(i=0;iK;i+) for(j=0;j0; i-) if(M%i=0)&(N%i=0)&(i=group_size) return i; return 1;void printResult() int i,j; printf(nOutput of Matrix C = ABn); for(i=0;iM;i+) for

9、(j=0;jN;j+) printf(%ft,C(i,j); printf(n); endtime=MPI_Wtime(); printf(n); printf(Whole running time = %f secondsn,endtime-starttime); printf(Distribute data time = %f secondsn,time1-starttime); printf(Parallel compute time = %f secondsn,endtime-time1);int main(int argc, char *argv) int i,j,k,l,group

10、_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;ip;i+) MPI_Send(&M,1,MPI_INT,i,i,MPI_COMM_WORLD); MPI_Send(&K,1,MPI_INT,i,i,MPI_COMM_WORLD); MPI_Send(&N,1,MPI_INT,i,i,MPI_CO

11、MM_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); m=M/p; n=N/p; if(myidp) a=(float *)malloc(floatsize*m*K); b=(float *)malloc(floatsize*K*n); c=(floa

12、t *)malloc(floatsize*m*N); if (myid%2!=0) buffer=(float *)malloc(K*n*floatsize); if (a=NULL|b=NULL|c=NULL) printf(Allocate space for a,b or c fail!); if (myid=0) for (i=0;im;i+) for (j=0;jK;j+) a(i,j)=A(i,j); for (i=0;iK;i+) for (j=0;jn;j+) b(i,j)=B(i,j); if (myid=0) for (i=1;ip;i+) MPI_Send(&A(m*i,

13、0),K*m,MPI_FLOAT,i,i,MPI_COMM_WORLD); for (j=0;jK;j+) MPI_Send(&B(j,n*i),n,MPI_FLOAT,i,i,MPI_COMM_WORLD); free(A); free(B); else MPI_Recv(a,K*m,MPI_FLOAT,0,myid,MPI_COMM_WORLD,&status); for (j=0;jK;j+) MPI_Recv(&b(j,0),n,MPI_FLOAT,0,myid,MPI_COMM_WORLD,&status); if (myid=0) time1=MPI_Wtime(); for (i

14、=0;ip;i+) l=(i+myid)%p; for (k=0;km;k+) for (j=0;jn;j+) for (c(l,k,j)=0,s=0;sK;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,MPI_FLOAT,mm1,mm1,MPI_COMM_WORLD); MPI_Recv(b,K*n,MPI_FLOAT,mp1,myid,MPI_COMM_WORLD,&status); else for(k=0;kK;k+) for(

15、j=0;jn;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;im;i+) for(j=0;jN;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;kp;k+) MPI_Recv(c,m

16、*N,MPI_FLOAT,k,k,MPI_COMM_WORLD,&status); for(i=0;im;i+) for(j=0;jN;j+) C(k*m+i),j)=*(c+i*N+j); if(myid=0) printResult(); MPI_Finalize(); if(myidp) free(a); free(b); free(c); if(myid=0) free(C); if(myid%2!=0) free(buffer); return (0);1.2運(yùn)行結(jié)果圖4.1 4個(gè)處理器的運(yùn)行結(jié)果2、網(wǎng)孔連接2.1源程序代碼#include #include #include #in

17、clude #include #include /* 全局變量聲明 */float *A, *B, *C; /* 總矩陣,C = A * B */float *a, *b, *c, *tmp_a, *tmp_b; /* a、b、c表分塊，tmp_a、tmp_b表緩沖區(qū) */int dg, dl, dl2,p, sp; /* dg:總矩陣維數(shù);dl:矩陣塊維數(shù);dl2=dl*dl;p:處理器個(gè)數(shù);spsqrt(p) */int my_rank, my_row, my_col; /* my_rank:處理器ID;(my_row,my_col):處理器邏輯陣列坐標(biāo) */MPI_Status stat

18、us;/* *函數(shù)名: get_index *功能：處理器邏輯陣列坐標(biāo)至rank號(hào)的轉(zhuǎn)換 *輸入：坐標(biāo)、邏輯陣列維數(shù) *輸出：rank號(hào) */int get_index(int row, int col, int sp) return (row+sp)%sp)*sp + (col+sp)%sp;/* *函數(shù)名：random_A_B *功能：隨機(jī)生成矩陣A和B */void random_A_B() int i,j; float m; /srand(unsigned int)time(NULL); /*設(shè)隨機(jī)數(shù)種子*/*隨機(jī)生成A,B,并初始化C*/ for(i=0; idg ; i+) for

19、(j=0; jdg ; j+) scanf(%f,&m); Aij = m; Cij = 0.0;m=0; for(i=0; idg ; i+) for(j=0; jdg ; j+) scanf(%f,&m); Bij = m;m=0; /* 函數(shù)名：scatter_A_B * 功能：rank為0的處理器向其他處理器發(fā)送A、B矩陣的相關(guān)塊 */void scatter_A_B() int i,j,k,l; int p_imin,p_imax,p_jmin,p_jmax; for(k=0; kp; k+) /*計(jì)算相應(yīng)處理器所分得的矩陣塊在總矩陣中的坐標(biāo)范圍*/ p_jmin = (k % sp

20、 ) * 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，準(zhǔn)備向其他處理器發(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的處理器直接將自己對(duì)應(yīng)的矩陣塊從tmp_a,tmp_b拷至a,b*/ if(k=

21、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_WORLD); MPI_Send(tmp_b, dl2, MPI_FLOAT, k, 2, MPI_COMM_WORLD); /* *函數(shù)名:init_alignment *功能:矩陣A和B初始對(duì)準(zhǔn) */void init_alignment() MPI_

22、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(float) ); /*將B中坐標(biāo)為(i,j)的分塊B(i,j)向上循環(huán)移動(dòng)j步*/ MPI_Sendrecv(b, dl2, MPI_FLOAT, get_index(my_row-my_col,my_col,sp), 1, t

23、mp_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 *功能：分塊矩陣左移和上移，并計(jì)算分塊c */void main_shift() int i,j,k,l; for(l=0; lsp; l+) /*矩陣塊相乘，c+=a*b */ for(i=0; idl; i+) for(j=0; jdl; j+) for(k=0; kdl; k+) ci*dl+j +=

24、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(b , dl2, MPI_FLOAT, get_index(my_row-1, my_col, sp), 1, MPI_COMM_WORLD)

25、; MPI_Recv(b , dl2, MPI_FLOAT, get_index(my_row+1, my_col, sp), 1, MPI_COMM_WORLD, &status); /* *函數(shù)名：collect_c *功能：rank為0的處理器從其余處理器收集分塊矩陣c */void collect_C() int i,j,i2,j2,k; int p_imin,p_imax,p_jmin,p_jmax; /* 分塊矩陣在總矩陣中頂點(diǎn)邊界值 */ /* 將rank為0的處理器中分塊矩陣c結(jié)果賦給總矩陣C對(duì)應(yīng)位置 */ for (i=0;idl;i+) for(j=0;jdl;j+) Ci

26、j=ci*dl+j; for (k=1;kp;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=p_imax; i+)

27、j2=0; for(j=p_jmin; j=p_jmax; j+) Cij=ci2*dl+j2; j2+; i2+; /*函數(shù)名：print *功能：打印矩陣 *輸入：指向矩陣指針的指針，字符串 */void print(float *m,char *str) int i,j; printf(%s,str); /*打印矩陣m*/ for(i=0;idg;i+) for(j=0;jdg;j+) printf(%15.0f ,mij); printf(n); printf(n);/* *函數(shù)名：main *功能：主過(guò)程，Cannon算法，矩陣相乘 *輸入：argc為命令行參數(shù)個(gè)數(shù)，argv為每個(gè)命

28、令行參數(shù)組成的字符串?dāng)?shù)組 */int main(int argc, char *argv) int i; MPI_Init(&argc, &argv); /* 啟動(dòng)MPI計(jì)算 */ MPI_Comm_size(MPI_COMM_WORLD, &p); /* 確定處理器個(gè)數(shù) */ MPI_Comm_rank(MPI_COMM_WORLD, &my_rank); /* 確定各自的處理器標(biāo)識(shí)符 */ sp = sqrt(p); /* 確保處理器個(gè)數(shù)是完全平方數(shù)，否則打印錯(cuò)誤信息，程序退出 */ if (sp*sp != p) if (my_rank = 0) printf(Number of pro

29、cessors is not a quadratic number!n); MPI_Finalize(); exit(1); if (argc != 2) if (my_rank = 0) printf(usage: mpirun -np ProcNum cannon MatrixDimensionn); MPI_Finalize(); exit(1); dg = atoi(argv1); /* 總矩陣維數(shù) */ dl = dg / sp; /* 計(jì)算分塊矩陣維數(shù) */ dl2 = dl * dl; /* 計(jì)算處理器在邏輯陣列中的坐標(biāo) */ my_col = my_rank % sp ; my

30、_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; idl2 ; i+) ci = 0.0; /* 為tmp_a、tmp_b分配空間 */ tmp_a = (float *)malloc( dl2 * sizeof(float) ); tmp_b =

31、(float *)malloc( dl2 * sizeof(float) ); if (my_rank = 0) /* rank為0的處理器為A、B、C分配空間 */ A = (float *)malloc( dg * sizeof(float*) ); B = (float *)malloc( dg * sizeof(float*) ); C = (float *)malloc( dg * sizeof(float*) ); for(i=0; idg; i+) Ai = (float *)malloc( dg * sizeof(float) ); Bi = (float *)malloc(

32、dg * sizeof(float) ); Ci = (float *)malloc( dg * sizeof(float) ); random_A_B(); /* rank為0的處理器隨機(jī)化生成A、B矩陣 */ scatter_A_B(); /* rank為0的處理器向其他處理器發(fā)送A、B矩陣的相關(guān)塊 */ else /* rank不為0的處理器接收來(lái)自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_

33、COMM_WORLD, &status); init_alignment(); /* A、B矩陣的初始對(duì)準(zhǔn) */ main_shift(); /* 分塊矩陣左移、上移, cannon算法的主過(guò)程 */ if(my_rank = 0) collect_C(); /* rank為0的處理器從其余處理器收集分塊矩陣c */ print(A,random matrix A : n); /* 打印矩陣A */ print(B,random matrix B : n); /* 打印矩陣B */ print(C,Matrix C = A * B : n); /* 打印矩陣C */ else MPI_Send

34、(c,dl2,MPI_FLOAT,0,1,MPI_COMM_WORLD); MPI_Barrier(MPI_COMM_WORLD); /* 同步所有處理器 */ MPI_Finalize(); /* 結(jié)束MPI計(jì)算 */ return 0;2.2運(yùn)行結(jié)果圖4.2 4個(gè)處理器的運(yùn)行結(jié)果3、在數(shù)學(xué)軟件中的計(jì)算結(jié)果圖4.3 在MATLAB中的運(yùn)行結(jié)果五、算法分析、優(yōu)缺點(diǎn)1、簡(jiǎn)單的并行分塊乘法的過(guò)程為（1）分塊：將: Ann與 Bnn分成p塊Ai,j和Bi,j(0i,j ),每塊大小為，并將它們分配給個(gè)處理器 （ ）。開始時(shí)Pi,j存放在塊Ai,j與塊Bi,j，然后計(jì)算塊Ci,j。（2）通信：為了計(jì)算塊