并行计算矩阵分块乘法

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：二维环绕网孔连接已知求。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从当前节点V选路到节点为V1 （3.2） endwhileEnd二维网孔情形下矩阵的简单并行分块算法输入：待选路的信包处于源处理器中输出：将各信包送至各自的目的地中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设计流程以下是二维网孔与超立方连接设计流程。 如图3-1二维网孔步骤：(1)先进行行播送;(2)再同时进行列播送;37625149804444444442233331111213141510 图3-1 二维网孔示意图超立方步骤：依次从低维到高维播送, d-立方, d=0,1,2,3,4;算法流程如图所示：包括mpi的头文件相关变量声明MPI_INIT()MPI_COMM_RANK()MPI_COMM_SIZE()进入MPI系统矩阵内部的通信应用控制实体：矩阵内部的计算程序MPI_FINALIZE()退出MPI系统结束开始循环直至结束图3-

5、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(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运行结果图4.1 4个处理器的运行结果2、网孔连接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表缓冲区 */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 stat

18、us;/* *函数名: get_index *功能：处理器逻辑阵列坐标至rank号的转换 *输入：坐标、逻辑阵列维数 *输出：rank号 */int get_index(int row, int col, 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; 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; /* 函数名：scatter_A_B * 功能：rank为0的处理器向其他处理器发送A、B矩阵的相关块 */void scatter_A_B() int i,j,k,l; int p_imin,p_imax,p_jmin,p_jmax; for(k=0; kp; k+) /*计算相应处理器所分得的矩阵块在总矩阵中的坐标范围*/ 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中的相应块拷至tmp_a,tmp_b，准备向其他处理器发送*/ 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=

21、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_WORLD); MPI_Send(tmp_b, dl2, MPI_FLOAT, k, 2, MPI_COMM_WORLD); /* *函数名:init_alignment *功能:矩阵A和B初始对准 */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中坐标为(i,j)的分块B(i,j)向上循环移动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) );/* *函数名：main_shift *功能：分块矩阵左移和上移，并计算分块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); /* *函数名：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;idl;i+) for(j=0;jdl;j+) Ci

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

27、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;idg;i+) for(j=0;jdg;j+) printf(%15.0f ,mij); printf(n); printf(n);/* *函数名：main *功能：主过程，Cannon算法，矩阵相乘 *输入：argc为命令行参数个数，argv为每个命

28、令行参数组成的字符串数组 */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); /* 确保处理器个数是完全平方数，否则打印错误信息，程序退出 */ 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); /* 总矩阵维数 */ dl = dg / sp; /* 计算分块矩阵维数 */ dl2 = dl * dl; /* 计算处理器在逻辑阵列中的坐标 */ 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的处理器随机化生成A、B矩阵 */ scatter_A_B(); /* 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_

33、COMM_WORLD, &status); init_alignment(); /* A、B矩阵的初始对准 */ main_shift(); /* 分块矩阵左移、上移, cannon算法的主过程 */ 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(); /* 结束MPI计算 */ return 0;2.2运行结果图4.2 4个处理器的运行结果3、在数学软件中的计算结果图4.3 在MATLAB中的运行结果五、算法分析、优缺点1、简单的并行分块乘法的过程为（1）分块：将: Ann与 Bnn分成p块Ai,j和Bi,j(0i,j ),每块大小为，并将它们分配给个处理器 （ ）。开始时Pi,j存放在块Ai,j与块Bi,j，然后计算块Ci,j。（2）通信：为了计算块