并行計算矩陣分塊乘法

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

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

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

4、r i=0 to do for all par-do =+ endfor Endfor End3.2設計流程以下是二維網孔與超立方連接設計流程。 如圖3-1二維網孔步驟：(1)先進行行播送;(2)再同時進行列播送;37625149804444444442233331111213141510 圖3-1 二維網孔示意圖超立方步驟：依次從低維到高維播送, d-立方, d=0,1,2,3,4;算法流程如圖所示：包括mpi的頭文件相關變量聲明MPI_INIT()MPI_COMM_RANK()MPI_COMM_SIZE()進入MPI系統(tǒng)矩陣內部的通信應用控制實體：矩陣內部的計算程序MPI_FINALIZE

5、()退出MPI系統(tǒng)結束開始循環(huán)直至結束圖3-2 算法流程四、源程序代碼及運行結果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(

6、x,y) ax*K+y#define b(x,y) bx*n+y#define buffer(x,y) bufferx*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(); dataF

7、ile=fopen("yin.txt","r"); fscanf(dataFile,"%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); fscanf(dataFile,"%d%d", &P, &N); if (K!=P) printf("the

8、 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); printf("Input of file "yin.txt"n"); printf("%dt %dn",M, K); for(i=0;i<M;i+) for(j=0;j<

9、K;j+) printf("%ft",A(i,j); printf("n"); printf("%dt %dn",K, N); for(i=0;i<K;i+) for(j=0;j<N;j+) printf("%ft",B(i,j); printf("n"); C=(float *)malloc(floatsize*M*N); int gcd(int M,int N,int group_size) int i; for(i=M; i>0; i-) if(M%i=0)&&a

10、mp;(N%i=0)&&(i<=group_size) return i; return 1;void printResult() int i,j; printf("nOutput of Matrix C = ABn"); for(i=0;i<M;i+) for(j=0;j<N;j+) printf("%ft",C(i,j); printf("n"); endtime=MPI_Wtime(); printf("n"); printf("Whole running time

11、 = %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_size,mp1,mm1; MPI_Init(&argc,&argv); MPI_Comm_size(MPI_COMM_

12、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); 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

13、_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(myid<p) a=(float *)malloc(floatsize*m*K); b=(float *)malloc(floatsize*K*n); c=(float *)mallo

14、c(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;i<m;i+) for (j=0;j<K;j+) a(i,j)=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<

15、;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); else MPI_Recv(a,K*m,MPI_FLOAT,0,myid,MPI_COMM_WORLD,&status); for (j=0;j<K;j+) MPI_Recv(&b(j,0),n,MPI_FLOAT,0,myid,MPI_COMM_WORLD,

16、&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,MPI_FLOAT,mm1,mm1,MPI_COMM_WORLD); MPI_Recv(b,K*n,M

17、PI_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+) for(j=0;j<N;j+) C(i,j)=*(c+i*N+j); if

18、(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(b); free(c); if(myid=0) fre

19、e(C); if(myid%2!=0) free(buffer); return (0);1.2運行結果圖4.1 4個處理器的運行結果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; /* 總矩陣,C = A * B */float *a, *b, *c, *tmp_a, *t

20、mp_b; /* a、b、c表分塊，tmp_a、tmp_b表緩沖區(qū) */int dg, dl, dl2,p, sp; /* dg:總矩陣維數;dl:矩陣塊維數;dl2=dl*dl;p:處理器個數;spsqrt(p) */int my_rank, my_row, my_col; /* my_rank:處理器ID;(my_row,my_col):處理器邏輯陣列坐標 */MPI_Status status;/* *函數名: get_index *功能：處理器邏輯陣列坐標至rank號的轉換 *輸入：坐標、邏輯陣列維數 *輸出：rank號 */int get_index(int row, int col

21、, int sp) return (row+sp)%sp)*sp + (col+sp)%sp;/* *函數名：random_A_B *功能：隨機生成矩陣A和B */void random_A_B() int i,j; float m; /srand(unsigned int)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+)

22、for(j=0; j<dg ; j+) scanf("%f",&m); Bij = m;m=0; /* 函數名：scatter_A_B * 功能：rank為0的處理器向其他處理器發(fā)送A、B矩陣的相關塊 */void scatter_A_B() int i,j,k,l; int p_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 = (

23、k - (k % sp)/sp * dl; p_imax = (k - (k % sp)/sp +1) *dl -1; l = 0; /*rank=0的處理器將A,B中的相應塊拷至tmp_a,tmp_b，準備向其他處理器發(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的處理器直接將自己對應的矩陣塊從tmp_a,tmp_b拷至a,b*/ if(k=0) memcpy(a, tmp_a, dl2 * sizeof(float);

24、memcpy(b, tmp_b, dl2 * sizeof(float); else /*rank=0的處理器向其他處理器發(fā)送tmp_a,tmp_b中相關的矩陣塊*/ MPI_Send(tmp_a, dl2, MPI_FLOAT, k, 1, MPI_COMM_WORLD); MPI_Send(tmp_b, dl2, MPI_FLOAT, k, 2, MPI_COMM_WORLD); /* *函數名:init_alignment *功能:矩陣A和B初始對準 */void init_alignment() MPI_Sendrecv(a, dl2, MPI_FLOAT, get_index(my_

25、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中坐標為(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_ro

26、w+my_col,my_col,sp), 1, MPI_COMM_WORLD, &status); memcpy(b, tmp_b, dl2 * sizeof(float) );/* *函數名：main_shift *功能：分塊矩陣左移和上移，并計算分塊c */void main_shift() int i,j,k,l; for(l=0; l<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;

27、/* 將分塊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); MPI_Recv(b

28、, dl2, MPI_FLOAT, get_index(my_row+1, my_col, sp), 1, MPI_COMM_WORLD, &status); /* *函數名：collect_c *功能：rank為0的處理器從其余處理器收集分塊矩陣c */void collect_C() int i,j,i2,j2,k; int p_imin,p_imax,p_jmin,p_jmax; /* 分塊矩陣在總矩陣中頂點邊界值 */ /* 將rank為0的處理器中分塊矩陣c結果賦給總矩陣C對應位置 */ for (i=0;i<dl;i+) for(j=0;j<dl;j+) Cij

29、=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<=p

30、_imax; i+) j2=0; for(j=p_jmin; j<=p_jmax; j+) Cij=ci2*dl+j2; j2+; i2+; /*函數名：print *功能：打印矩陣 *輸入：指向矩陣指針的指針，字符串 */void print(float *m,char *str) int i,j; printf("%s",str); /*打印矩陣m*/ for(i=0;i<dg;i+) for(j=0;j<dg;j+) printf("%15.0f ",mij); printf("n"); printf(&quo

31、t;n");/* *函數名：main *功能：主過程，Cannon算法，矩陣相乘 *輸入：argc為命令行參數個數，argv為每個命令行參數組成的字符串數組 */int main(int argc, char *argv) int i; MPI_Init(&argc, &argv); /* 啟動MPI計算 */ MPI_Comm_size(MPI_COMM_WORLD, &p); /* 確定處理器個數 */ MPI_Comm_rank(MPI_COMM_WORLD, &my_rank); /* 確定各自的處理器標識符 */ sp = sqrt(p);

32、/* 確保處理器個數是完全平方數，否則打印錯誤信息，程序退出 */ if (sp*sp != p) if (my_rank = 0) printf("Number of processors 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(

33、argv1); /* 總矩陣維數 */ dl = dg / sp; /* 計算分塊矩陣維數 */ dl2 = dl * dl; /* 計算處理器在邏輯陣列中的坐標 */ my_col = my_rank % 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 *

34、/ 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) /* rank為0的處理器為A、B、C分配空間 */ A = (float *)malloc( dg * sizeof(float*) ); B = (float *)malloc( dg * sizeof(float*) ); C = (float *)

35、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的處理器隨機化生成A、B矩陣 */ scatter_A_B(); /* rank為0的處理器向其他處理器發(fā)送A、B矩陣的相關塊 */ else /* rank不為0的處理器接收來自ra

36、nk為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); init_alignment(); /* A、B矩陣的初始對準 */ main_shift(); /* 分塊矩陣左移、上移, cannon算法的主過程 */ if(my_rank = 0) collect_C(); /* rank為0的處理器從其余處理器收集分塊矩陣c */ print(A,"

37、;random matrix A : n"); /* 打印矩陣A */ print(B,"random matrix B : n"); /* 打印矩陣B */ print(C,"Matrix C = A * B : n"); /* 打印矩陣C */ else MPI_Send(c,dl2,MPI_FLOAT,0,1,MPI_COMM_WORLD); MPI_Barrier(MPI_COMM_WORLD); /* 同步所有處理器 */ MPI_Finalize(); /* 結束MPI計算 */ return 0;2.2運行結果圖4.2 4個處理器的運行結果3、在數學軟件中的計算結果圖4.3 在MATLAB中的運行結果五、算法分析、優(yōu)缺點1、簡單的并行分塊乘法的過程為（1）分塊：將: An×n與 B