Cannon矩阵乘法的MPI实现及性能分析

#includestdio.h#includempi.h#includemath.h#includestdlib.h#includetime.h#includememory.hMPI_Statusstatus;double**A,**B,**C;//C=A*Bdouble*a,*b,*c;//各个进程的缓冲区intn;//矩阵的行列数intnp;//每个进程控制的小矩阵的行列数intp,rank;//进程个个数、当前进程的编号,笛卡尔进程编号double*tempa,*tempb;voidProduceABC();//在根处理器中生成矩阵AB,初始化矩阵CvoidPrintABC();//输出结果voidScatterAB();//分发矩阵AB中的元素到各个进程中voidMainProcess();//cannon算法的主过程voidcollectC();//收集结果矩阵CvoidMutiply();//矩阵相乘voidPrintab();voidPrintc();intmain(intargc,char*argv[]){inti;doublestarttime,endtime;MPI_Init(&argc,&argv);MPI_Comm_size(MPI_COMM_WORLD,&p);MPI_Comm_rank(MPI_COMM_WORLD,&rank);if(rank==0){printf(请输入矩阵的行列数n=);fflush(stdout);scanf_s(%d,&n);printf(\n);}MPI_Bcast(&n,1,MPI_DOUBLE,0,MPI_COMM_WORLD);//n=atoi(argv[1]);np=n/(int)sqrt(p);a=(double*)malloc(np*np*sizeof(double));b=(double*)malloc(np*np*sizeof(double));c=(double*)malloc(np*np*sizeof(double));memset(c,0,np*np*sizeof(double));tempa=(double*)malloc(np*np*sizeof(double));tempb=(double*)malloc(np*np*sizeof(double));if(rank==0){//在根处理器中为矩阵ABC分配空间A=(double**)malloc(n*sizeof(double*));B=(double**)malloc(n*sizeof(double*));C=(double**)malloc(n*sizeof(double*));for(i=0;in;i++){A[i]=(double*)malloc(n*sizeof(double));B[i]=(double*)malloc(n*sizeof(double));C[i]=(double*)malloc(n*sizeof(double));}ProduceABC();//在根处理器中随机生成矩阵AB，初始化矩阵CScatterAB();//分发矩阵AB中的元素到各个进程中}else{MPI_Recv(a,np*np,MPI_DOUBLE,0,1,MPI_COMM_WORLD,&status);MPI_Recv(b,np*np,MPI_DOUBLE,0,2,MPI_COMM_WORLD,&status);}starttime=MPI_Wtime();//开始时间MainProcess();//cannon算法的主过程if(rank==0){collectC();//收集结果矩阵CPrintABC();//输出结果endtime=MPI_Wtime();printf(timeused:%lf\n,endtime-starttime);for(i=0;in;i++){free(A[i]);free(B[i]);free(C[i]);}free(A);free(B);free(C);}else{MPI_Send(c,np*np,MPI_DOUBLE,0,1,MPI_COMM_WORLD);}free(a);free(b);free(c);free(tempa);free(tempb);MPI_Finalize();return0;}voidProduceABC()//在根处理器中生成矩阵AB{inti,j;for(i=0;in;i++)for(j=0;jn;j++){A[i][j]=1.0;B[i][j]=1.0;C[i][j]=0.0;}}voidPrintABC()//输出结果{printf(A[0][0]=%f\nB[0][0]=%f\nC[0][0]=%f\n,A[0][0],B[0][0],C[0][0]);}voidScatterAB()//分发矩阵AB中的元素到各个进程中{intimin,imax,jmin,jmax;intsp;inti,j,k,m;for(k=0;kp;k++){/*计算相应处理器所分得的矩阵块在总矩阵中的坐标范围*/sp=(int)sqrt(p);imin=(k/sp)*np;imax=imin+np-1;jmin=(k%sp)*np;jmax=jmin+np-1;/*rank=0的处理器将A,B中的相应块拷至tempa,tempb，准备向其他处理器发送*/m=0;for(i=imin;i=imax;i++){for(j=jmin;j=jmax;j++){tempa[m]=A[i][j];tempb[m]=B[j][i];//矩阵B按列优先存储m++;}}/*根处理器将自己对应的矩阵块从tempa,tempb拷至a,b*/if(k==0){memcpy(a,tempa,np*np*sizeof(double));memcpy(b,tempb,np*np*sizeof(double));}else/*根处理器向其他处理器发送tempa,tempb中相关的矩阵块*/{MPI_Send(tempa,np*np,MPI_DOUBLE,k,1,MPI_COMM_WORLD);MPI_Send(tempb,np*np,MPI_DOUBLE,k,2,MPI_COMM_WORLD);}}}voidMainProcess()//cannon算法的主过程{MPI_Commcomm;//笛卡尔结构通讯器intcrank;intdims[2],periods[2],coords[2];intsource,dest,up,down,right,left;inti;dims[0]=dims[1]=(int)sqrt(p);periods[0]=periods[1]=1;MPI_Cart_create(MPI_COMM_WORLD,2,dims,periods,1,&comm);MPI_Comm_rank(comm,&crank);MPI_Cart_coords(comm,crank,2,coords);MPI_Cart_shift(comm,1,-1,&right,&left);MPI_Cart_shift(comm,0,-1,&down,&up);MPI_Cart_shift(comm,1,-coords[0],&source,&dest);MPI_Sendrecv_replace(a,np*np,MPI_DOUBLE,dest,1,source,1,comm,&status);MPI_Cart_shift(comm,0,-coords[1],&source,&dest);MPI_Sendrecv_replace(b,np*np,MPI_DOUBLE,dest,1,source,1,comm,&status);Mutiply();//矩阵相乘for(i=1;idims[0];i++){MPI_Sendrecv_replace(a,np*np,MPI_DOUBLE,left,1,right,1,comm,&status);MPI_Sendrecv_replace(b,np*np,MPI_DOUBLE,up,1,down,1,comm,&status);Mutiply();//矩阵相乘}MPI_Comm_free(&comm);}voidcollectC()//收集结果矩阵C{inti,j,k,s,m;intimin,imax,jmin,jmax;intsp=(int)sqrt(p);/*根处理器中的c赋给总矩阵C*/for(i=0;inp;i++){for(j=0;jnp;j++)C[i][j]=c[i*np+j];}for(k=1;kp;k++){/*根处理器从其他处理器接收相应的分块c*/MPI_Recv(c,np*np,MPI_DOUBLE,k,1,MPI_COMM_WORLD,&status);//printf(rank=%d\n,k);Printc();imin=(k/sp)*np;imax=imin+np-1;jmin=(k%sp)*np;jmax=jmin+np-1;/*将接收到的c拷至C中的相应位置,从而构造出C*/for(i=imin,m=0;i=imax;i++,m++){for(j=jmin,s=0;j=jmax;j++,s++)C[i][j]=c[m*np+s];}}}voidMutiply()//矩阵相乘{inti,j,k;for(i=0;inp;i++)for(j=0;jnp;j++)for(k=0;knp;k++)c[i*np+j]+=a[i*np+k]*b[j*np+k];//b按列优先来搞}注意：1.以.c来编译2.realease需单独配环境，并以realease跑比较快3.64位的mpich2要以x64平台编译，否则通不过4.mpiexec中一定要勾选以单独窗口运行，否则跑不动5.矩阵乘法，第二个矩阵按列优先存储性能提升数倍。若以二维存储，先转置再乘，同样的效果。