xu.cpp

vc并行算法 论文

// Canno's MM algorithm.
#include<stdlib.h>
#include<mpi.h>
//将处理器的网格坐标转换为处理器号

int CarToRank(int row,int col,int npdim)
{
	return((row+npdim)%npdim)*npdim+(col+npdim)%npdim;
}
//块矩阵相乘
 void MatrixMultiply(int nlocal,double *a,double *b,double *c)
 {
	 int i,j,k;
	 for(i=0;i<nlocal;i++)
		 for(j=0;j<nlocal;j++)
			 for(k=0;k<nlocal;k++)
			 {
				 c[i,j]+=a[i,j]*b[k,j];
			 }
 }
 void main(int argc,char**argv)
 {
	 int i;
	 int n,nlocal;
	 double *a,*b,*c,*wspace;
	 int npes,npdim;
	 int myrank,myrow,mycol;
	 MPI_Status status;
	 /*初始化，得到系统信息*/
	 MPI_Init(&argc,&argv);
	 MPI_Comm_size(MPI_COMM_WORLD,&npes);
     MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
	 if(argc!=2)
	 {//需要给出矩阵的维数
		 if(myrank==0)
	printf("Usage:%s<the dimension of the matrix>\n",argv[0]);
		 MPI_Finalize();
		 exit(0);
	 }
	 //矩阵的维数
	 n=atoi(argv[1]);
	 //处理器网格每一维的处理器数P的平方根

	 npdim=sqrt(npes);
	 if(npdim*npdim!=npes)
	 {
		 处理器数目必须是平方数
       if(myrank==0)
    printf("The number of processes must be a perfect square.\n");
	   MPI_Finalize();
	   exit(0);
	 }
	 //本地的块矩阵的维数
	 nlocal=n/npdim;
	 //根据处理器号得到自己在处理器网格的坐标（myrow,mycol)
	 myrow=myrank/npdim;
	 mycol=myrank%npdim;
//分配本地计算所需的空间
	 a=(double*)malloc(nlocal*nlocal*sizeof(double));
     b=(double*)malloc(nlocal*nlocal*sizeof(double));
     c=(double*)malloc(nlocal*nlocal*sizeof(double));
     wspace=(double*)malloc(nlocal*nlocal*sizeof(double));
	 //经矩阵赋初值（a,b的内容实际上没有实际意义）。
	 srand48((long)myrank);
	 for(i=0;i<nlocal*nlocal;i++)
	 {a[i]=b[i]=drand48();
	 c[i]=0.0;
	 }
	 /*Perform the initial matrix alignment.First for a and then for b*/
	/*进行初始对齐，先a再b。设处理器坐标为(i,j)，则本地的a将传给p(i,j-i)，
	 即左移i列，b将传给p(i-j,j)，即上移j列*/
	 /*MPI_Sendrecv的形式
	 int MPI_Sendrecv(void*sendbuf,int sendcount,MPI_Datatype senddatatype,
	 int dest,int sendtag,void*recvbuf,int recvcount,
	 MPI_Datatype recvdatatype, int source, int recvtag,
	 MPI_Comm comm,MPI_Status *status)
	 */

	MPI_Sendrecv(a,nlocal*nlocal,MPI_DOUBLE,
		 CartToRank(myrow,mycol-myrow,npdim),
		 1,wspace,nlocal*nlocal,MPI_DOUBLE,
		 CartToRank(myrow,mycol+myrow,npdim),
		 1,MPI_COMM_WORLD,&status);

	 memcpy(a,wspace,nlocal*nlocal*sizeof(double));

	 MPI_Sendrecv(b,nlocal*nlocal,MPI_DOUBLE,
		 CartToRank(myrow-mycol,mycol,npdim),
		 1,wspace,nlocal*nlocal,MPI_DOUBLE,
		 CartToRank(myrow+mycol,mycol,npdim),
		 1,MPI_COMM_WORLD,&status);
memcpy(b,wspace,nlocal*nlocal*sizeof(double));
/*主循环*/
   for(i=0;i<npdim;i++)
   {MatrixMultiply(nlocal,a,b,c);  /*c=c+a*b*c*/
   /*a左移*/
   MPI_Sendrecv(a,nlocal*nlocal,MPI_DOUBLE,
		 CartToRank(myrow,mycol-1,npdim),
		 1,wspace,nlocal*nlocal,MPI_DOUBLE,
		 CartToRank(myrow,mycol+1,npdim),
		 1,MPI_COMM_WORLD,&status);

	 memcpy(a,wspace,nlocal*nlocal*sizeof(double));
/*b上移1*/
   MPI_Sendrecv(b,nlocal*nlocal,MPI_DOUBLE,
		 CartToRank(myrow-1,mycol,npdim),
		 1,wspace,nlocal*nlocal,MPI_DOUBLE,
		 CartToRank(myrow+1,mycol,npdim),
		 1,MPI_COMM_WORLD,&status);
   memcpy(b,wspace,nlocal*nlocal*sizeof(double));
  
   /*得到原来的a,b，不是必须，但可以检验得法实现*/
 MPI_Sendrecv(a,nlocal*nlocal,MPI_DOUBLE,
		 CartToRank(myrow,mycol+myrow,npdim),
		 1,wspace,nlocal*nlocal,MPI_DOUBLE,
		 CartToRank(myrow,mycol-myrow,npdim),
		 1,MPI_COMM_WORLD,&status);

	 memcpy(a,wspace,nlocal*nlocal*sizeof(double));

	 MPI_Sendrecv(b,nlocal*nlocal,MPI_DOUBLE,
		 CartToRank(myrow+mycol,mycol,npdim),
		 1,wspace,nlocal*nlocal,MPI_DOUBLE,
		 CartToRank(myrow-mycol,mycol,npdim),
		 1,MPI_COMM_WORLD,&status);
memcpy(b,wspace,nlocal*nlocal*sizeof(double));
//算法完成

 if(myrank==0)
	 printf("Done multiplying the matrices!\n");
 free(a);
 free(b);
 free(c);
 free(wspace);
 MPI_Finalize();
   }