mpipoisson.cpp

Jacobi迭代法求解Poisson方程的MPI并行编程。可以接受参数作为问题剖分大小

//*************************************************************************
//
//							并行计算程序设计代码
//									之
//				Jacobi迭代法求解Poisson方程-Δu(x,y)=f(x,y)
//
//
//				设计：廖能想
//
//				日期：2008年4月13日
//
//				Emal：[e-mail]
//
//	简介：
//
//		1、可以在命令行中使用参数来自定义问题大小（Mx×My），如，
//			mpiexec.exe -np 4 xxx.exe 128 128	 (Window环境)。
//			mpirun -np 4 xxx 128 128  (Linux环境)
//
//		2、程序可以根据进程数自动调整并行任务划分。
//
//		3、计算结果将输出到两个目标：控制台和文件。
//			控制台：	内容为任务信息及计算终止代码
//			文件：	内容为任务大小、剖分步长、计算解和真解（用于测试比较）
//
//

#define MPICH_SKIP_MPICXX	// 用于VC6.0+WinXP环境的编译
#define DA	3.0
#define DB	3.0				// 区域Ω边界
#define DEF_NX	64
#define DEF_NY	64			// 默认规模
#define MAX_NUM_XY	100000	// 问题规模限制
#define NORM	1			// 收敛范数：1 无穷范数 0 二范数
#define Max_Count	9999999	// 最大迭代次数
#define EP	1e-14			// 误差下限
#define PI	3.141592653579
#define MAXDOUBLE	1.797693E+308
#include "mpi.h"
#include <stdio.h>
#include <math.h>
#include <memory.h>
#pragma comment(lib,"mpi.lib")

// 原函数
double uxy(double x,double y) {
	x-=1.5,y-=1.5;
	return  -cos( pow(x,3)+pow(y,3) );
}

// 右端函数
double fxy(double x,double y) {
	double temp;
	x-=1.5, y-=1.5;
	temp =9.0*cos( pow(x,3)+pow(y,3) )*(pow(x,4)+pow(y,4));
	temp+=6.0*sin( pow(x,3)+pow(y,3) )*(x+y);
	return -temp;
}

// 边界
double gxy(double x,double y) {
	return uxy(x,y);
}

/*	任务自适应划
*
*	·只是进程布局的划分·
*	·不一定能使得问题网格均匀划分·
*
*	要求：	mx<my
*	原理：
*		首先选取长度不占优的方向（x）做基准，按两个
*	方向上的长度比例来计算该方向上的进程划分（px），
*	并以此值为起点向数值1依次搜索直至能找到满足进程
*	总数size被px整除，这时px和py=size/px即为最佳划分。
*/
void _AutoCarve(
				int mx,int my,int size,
				int &px,int &py	)
{
	px=int(size/(1.0*my/mx+1))+1;
	for(py=size/px;py*px!=size;py=size/(--px));
}

/*	主函数
*
*	·包括所有的运算，虽然没有进行拆解成函数模块，不过已分层次·
*
*	可接受参数 (Mx,My)或使用默认值 (DEF_NX, DEF_NY)
*
*/
int main(int argc,char *argv[])
{	
	//	┟┈┈← nx+2 →┈┈┧
	//	┃              　  ┊
	//	┃   (*,?,?,…,?)   ↑
	//	┃        ...      ny+2
	//	┃   (*,?,?,…,?)   ↓
	//	┃                  ┊
	//	┞┈┈<包括边界>┈┈┦
	
	int rank,size;			// 进程信息
	
	int MX,MY;				// 网格大小（问题规模）
	
	int px,py;				// 数据矩阵的进程分块
	
	int nx,ny,sx,sy;		// 进程任务块信息：n 大小 s 起点
	
	double dx,dy;			// 剖分步长
	
	int i,o,j,k,l;			// 辅助变量
	
	int MSize;				// 矩阵大小
	
	double *u,*u0,*u1,*f;	// 任务数据·矩阵

	MPI_Status sta;			// 状态，用于消息接收等
	
	// 数据类型，用于数据传输
	
	MPI_Datatype vstype;	// 纵向
	
	MPI_Datatype hstype;	// 横向
	
	int nSideL, nSideR, nSideU, nSideD;		// 周边进程
	
	// MPI 初始化
	
	MPI_Init(&argc,&argv);
	
	MPI_Comm_rank(MPI_COMM_WORLD,&rank);
	
	MPI_Comm_size(MPI_COMM_WORLD,&size);
	
	{	//*********************** 数据初始化 ********************
		
		{	// 确定区域网格大小 如参数合法则使用参数的设置，否则取默认值
			
			if( argc==3 ) {
				
				i=0;
				
				if(sscanf(argv[1],"%d",&MX)<1) i++;
				
				if(i || sscanf(argv[2],"%d",&MY)<1) i++;
				
			} else i=1;
			
			if( i || !(	MX>2 && MX<MAX_NUM_XY && 
				MY>2 && MY<MAX_NUM_XY) )
			{
				MX=DEF_NX,MY= DEF_NY;
			}
		}
		
		{	// 按进程数划分任务块
			
			py=1,px=1;
			
			if(size>1) {
				if(MX<MY) _AutoCarve(MX,MY,size,px,py);
				else _AutoCarve(MY,MX,size,py,px);
			}
		}
		
		// 矩阵块大小
		nx=MX/px;if(nx*px<MX) nx++;
		ny=MY/py;if(ny*py<MY) ny++;
		
		// 修正剖分数(也可以划分成不均匀的块，那么就不必修改剖分值)
		MX=nx*px, MY=ny*py;
		
		// 等分步长
		dx=double(DA)/(MX+1.0), dy=double(DB)/(MY+1.0);
		
		// 确定任务范围
		sx=rank%px*nx+1, sy=rank/px*ny+1;
		MSize=(nx+2)*(ny+2);		// 包含边界 : +2
		
		MPI_Type_vector( ny, 1, nx+2, MPI_DOUBLE, &hstype );
		MPI_Type_commit(&hstype);	// 横向数据传输类型提交
		
		MPI_Type_vector( 1, nx, nx+2, MPI_DOUBLE, &vstype );
		MPI_Type_commit(&vstype);	// 纵向数据传输类型提交
		
		// 空间分配
		u=new double[2*MSize], u0=u, u1=u0+MSize;
		f=new double[nx*ny];	// 由于不必考虑f的边界，故可少分配空间
		memset(u,0,sizeof(double)*2*MSize);
		memset(f,0,sizeof(double)*(nx*ny));
		
		{	// 周边进程
			if( rank-px>=0 ) nSideU =rank-px;
			else nSideU=MPI_PROC_NULL;	// 上
			
			if( rank+px<size ) nSideD=rank+px;	
			else nSideD=MPI_PROC_NULL;	// 下
			
			if( rank>0 && (rank-1)/px==rank/px ) nSideL=rank-1;
			else nSideL =MPI_PROC_NULL;	// 左
			
			if( (rank+1)/px==rank/px ) nSideR =rank+1;
			else nSideR =MPI_PROC_NULL;	// 右
		}
		{	// Dirichlet 边界
			if(MPI_PROC_NULL== nSideU) {	//  上
				j=1,k=rank%px*nx+1;
				for(i=0;i<nx;i++) *(u1+i+j)=*(u0+i+j)=gxy(dx*(k+i),0);
			}
			
			if(MPI_PROC_NULL== nSideD) {	//  下
				j=MSize-2,k=rank%px*nx+nx;
				for(i=0;i<nx;i++) *(u1+j-i)=*(u0+j-i)=gxy(dx*(k-i),DB);
			}
			
			if(MPI_PROC_NULL== nSideL) {	//  左
				j=nx+2,k=rank/px*ny+1;
				for(i=0;i<ny;i++) *(u1+j+i*j)=*(u0+j+i*j)=gxy(0,dy*(k+i));
			}
			
			if(MPI_PROC_NULL== nSideR) {	//  右
				j=nx+2,o=j+j-1,k=rank/px*ny+1;
				for(i=0;i<ny;i++) *(u1+o+i*j)=*(u0+o+i*j)=gxy(DA,dy*(k+i));
			}
		}
		{	// 右端函数 f(x,y)
			for(j=0;j<ny;j++) {
				for(i=0;i<nx;i++) {
					f[i+j*nx]=fxy((sx+i)*dx,(sy+j)*dy);
				}
			}
		}
	}	//*********************** 数据初始化结束 ***********************
	
	{	// Jacobi 迭代求解
		
		int niter=0;		// 消息序号标志
		
		double err=MAXDOUBLE,errT,errA,*temp;	// 误差
		
		double t1=dx*dx,t2=dy*dy,t3=t1*t2,t4=0.5/(t1+t2);	// 迭代式中的常数
		
		int count=0; 	// 计数
		
		double time_a=MPI_Wtime(),time_b;	// 计时
		
		int nExitCode=0;	// 循环跳出标志
		
		MPI_Op iOp=NORM==1?MPI_MAX:MPI_SUM;	// 归约操作类型
		
		while(true && !nExitCode) {	// 迭代循环
			
			niter=niter<MPI_TAG_UB?niter+1:1;
			
			for(j=0,errT=0;j<ny;j++) {
				
				for(i=0;i<nx;i++) {
					
					u0[i+1+(j+1)*(nx+2)]=  (	t3 * f[i+j*nx] + \
						t1 * (u1[i+1+j*(nx+2)] + u1[i+1+(j+2)*(nx+2)]) + \
						t2 * (u1[i+(j+1)*(nx+2)] + u1[i+2+(j+1)*(nx+2)]) \
						)*t4;
					
					errA=u0[i+1+(j+1)*(nx+2)]-u1[i+1+(j+1)*(nx+2)];
					
					if (NORM==1) {	// 无穷范数
						if( fabs(errA)>errT ) errT=fabs(errA);
					} else errT+=errA*errA;
				}
			}
			
			if(size>0) {
				MPI_Allreduce(	&errT,&errA,1,MPI_DOUBLE,iOp,MPI_COMM_WORLD );
				errT=errA;
			}
			
			if (NORM==1) {	// 无穷范数
				errA=errT;
			} else  errA=sqrt(errT);
			
			if( errA>err) {
				fprintf(stderr,"\n开始发散！");
				nExitCode|=0x4;
			}
			
			err=errA, temp=u1, u1=u0, u0=temp;
			
			if(err<EP) {	// 收敛
				nExitCode|=0x1;
			}
			
			if( count>Max_Count ) {	// 迭代次数不足
				nExitCode|=0x2;
			} else count++;
			
			{	// 边界交互信息·常规消息模式
				
				// 上
				MPI_Send( u1+nx+2+1, 1, vstype, nSideU, niter, MPI_COMM_WORLD);
				
				// 下
				MPI_Send( u1+(ny+0)*(nx+2)+1, 1, vstype, nSideD,niter, MPI_COMM_WORLD);
				
				// 左
				MPI_Send( u1+nx+2+1, 1, hstype, nSideL,niter, MPI_COMM_WORLD);
				
				// 右
				MPI_Send( u1+nx+2+nx, 1, hstype, nSideR, niter, MPI_COMM_WORLD);
				
				// 上
				MPI_Recv( u1+1, 1, vstype, nSideU, niter, MPI_COMM_WORLD, &sta);
				
				// 下
				MPI_Recv( u1+(ny+1)*(nx+2)+1, 1, vstype, nSideD, niter, MPI_COMM_WORLD, &sta);
				
				// 左
				MPI_Recv( u1+nx+2, 1, hstype, nSideL,  niter, MPI_COMM_WORLD, &sta);
				
				// 右
				MPI_Recv(	u1+nx+2+nx+1, 1, hstype, nSideR, niter, MPI_COMM_WORLD, &sta);
			}
		}
		time_b=MPI_Wtime();	// 计算模块结束时间
		
		MPI_Type_free(&hstype);
		
		MPI_Type_free(&vstype);
		
		if(true && rank==0) {	// 输出任务信息
			
			fprintf(stdout,"\n%d<P : %d > - %.18lf : %lf %lf",count,rank,err,dx,dy);
			
			fprintf(stdout,"\n%d %d %d %d --->>> size:%d %d",nx,ny,px,py,MX,MY);
		}
		
		fprintf( stdout, "\nProcesser<%d-(%d,%d)> -- %s:", \
			rank, sx, sy, (NORM==1?"无穷范数":"二范数") );
		
		fprintf( stdout,"\n   Error:%.18lf Time %.18lf Count: %d Exit: %d",\
			err, time_b-time_a, count, nExitCode );
	}

	{	// 真解 u(x,y) 用于分析误差

		for(j=1;j<ny+1;j++) {
			for(i=1;i<nx+1;i++) {
				u0[i+j*(nx+2)]=uxy((sx+i-1)*dx,(sy+j-1)*dy);
			}
		}
	}
	{	//	输出到文件
		//
		//	注 : 使用默认文件格式 native
		//
		//	结构 : double(Mx,My,dx,dy),u(MPI_DOUBLE MX*MY)
		
		char filename[256]={0};
		
		sprintf(filename,"Out_NO2_%d_%d.dat",MX,MY);
		
		MPI_Datatype filetype,wtype;
		
		MPI_Type_vector(ny,nx,MX,MPI_DOUBLE,&filetype);
		MPI_Type_commit(&filetype);
		
		MPI_Type_vector(ny,nx,nx+2,MPI_DOUBLE,&wtype);
		MPI_Type_commit(&wtype);
		
		MPI_Offset disp;
		
		MPI_Type_size(MPI_DOUBLE,&i);
		
		disp=(rank/px*ny*MX+rank%px*nx+4)*i;
		
		MPI_File fh;
		
		MPI_File_open (	MPI_COMM_WORLD, filename, MPI_MODE_CREATE+\
			MPI_MODE_WRONLY, MPI_INFO_NULL, &fh );
		
		if(0==rank) {
			
			double M[4]={(double)MX,(double)MY,dx,dy};
			
			MPI_File_write( fh, M, 4, MPI_DOUBLE, &sta);
		}
		
		// 文件窗口
		
		MPI_File_set_view(fh,disp,MPI_DOUBLE,filetype,"native",MPI_INFO_NULL);
		
		MPI_File_write(fh,u1+nx+3,1,wtype,&sta);
		
		MPI_File_seek(fh,(MX*MY-nx*ny)*sizeof(MPI_DOUBLE),MPI_SEEK_CUR);
		
		MPI_File_write(fh,u0+nx+3,1,wtype,&sta);
		
		MPI_File_close(&fh);
		
		MPI_Type_free(&wtype);
		
		MPI_Type_free(&filetype);
	}
	
	delete []f;
	
	delete []u;
	
	MPI_Finalize();
	
	return 1;
}