并行计算矩阵分块乘法

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 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系统矩阵内部的通信应用控制实体：矩阵内部的计算程序MPI_FINALIZE

5、()退出MPI系统结束开始循环直至结束图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表缓冲区 */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的处理器向其他处理器发送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，准备向其他处理器发送*/ 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的处理器向其他处理器发送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)向上循环移动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的处理器向其他处理器发送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中的运行结果五、算法分析、优缺点1、简单的并行分块乘法的过程为（1）分块：将: An×n与 B