jacobi_mpi.cpp

Jacobi_MPI. ln2 & pi 的计算,基于Romberg and Simpson. MPI_File.

#define MPICH_SKIP_MPICXX
#include "mpi.h"
#include "stdio.h"
#include "math.h"
#include "stdlib.h"
#include <MEMORY.H>
#include "dos.h"
#define eps 0.0000001
#define n 4
#pragma comment(lib,"mpi.lib")

int iterprocess(double *x0,double *ai,double *bi,double *new_x)
{	
	int rank=0,size,ndiag=rank,i,j;
	double sum;
	MPI_Comm_rank(MPI_COMM_WORLD,&rank);         //得到当前正在运行的进程标识号，放在rank中
	MPI_Comm_size(MPI_COMM_WORLD,&size);         //得到所有参加运算的进程的个数，放在size中
	//for(int p=0;p<size;p++){					 //p:进程序号
	int num=n/size; // if n%size ==0
	for(j=0;j<num;j++){
		ndiag=rank*num+j; // line		
		for(i=0,sum=0;i<n;i++){
			sum+=(i==ndiag?0:-ai[i+j*n])*x0[i];
		}
		sum+=bi[ndiag];
		new_x[j]=sum/ai[ndiag+j*n];
	}
	return 0;
}


int getresidul(double *a,double *x0,double *b,double *res)
{
	for(int i=0;i<n;i++){
		res[i]=0;
		for(int j=0;j<n;j++)
			res[i]+=a[i+j*n]*x0[j];
	}
	for(i=0;i<n;i++)  
		res[i]=res[i]-b[i];                                 //残量r=Ax-b
	return 0;
}


int getl2norm(double *res,double *rnorm)          
{	 
	int i;
	for(i=0,*rnorm=0.0;i<n;i++)
		*rnorm+=res[i]*res[i];       
	*rnorm=sqrt(*rnorm);                                   //残量的二范数
	return 0;
}


void main(int argc,char *argv[]){
	int rank,size,niter;
	double a[n*n],b[n];
	double new_x[n],x0[n],ai[n*n],res[n],bi[n],rnorm,rnorm0;
	MPI_Init(&argc,&argv);
	MPI_Comm_rank(MPI_COMM_WORLD,&rank);         //得到当前正在运行的进程标识号，放在rank中
	MPI_Comm_size(MPI_COMM_WORLD,&size);		 //得到所有参加运算的进程的个数，放在size中
	if(rank==0){
	//read the data A[] []  ,b[]      
	//MPI_File_open(MPI_Comm comm, char * filename, int amode, MPI_Info info,MPI_File * fh)
		{
		double temp[]={4,0,1,1,0,2,0,1,1,0,4,0,0,1,0,2};
		memcpy(a,temp,sizeof(temp));
		}
		b[0]=6.0;
		b[1]=3.0;
		b[2]=5.0;
		b[3]=3.0;
		getl2norm(b,&rnorm0);
		x0[0]=0.0;
		x0[1]=0.0;
		x0[2]=0.0;
		x0[3]=0.0;
	} 
	
	MPI_Scatter(a,n*(n/size),MPI_DOUBLE,ai,n*(n/size),MPI_DOUBLE,0,MPI_COMM_WORLD);
	MPI_Scatter(b,n/size,MPI_DOUBLE,bi,n/size,MPI_DOUBLE,0,MPI_COMM_WORLD);
		{
			//测试发送数据是否成功
			/*printf("\n%d",rank);
			for(int i=0;i<n*(n/size);printf("(%lf) ",ai[i++]));
			for(int i=0;i<(n/size);printf("(%lf) ",bi[i++]));*/
		}
	niter=0; 		//迭代次数	
	do{
		MPI_Bcast(x0,n,MPI_DOUBLE,0,MPI_COMM_WORLD);
				{
					//测试发送数据是否成功
					//printf("\n%d",rank);
					//for(int i=0;i<n;i++,printf("Process %d got %lf,rank,*x0++));
				}
		
		//各进程完迭代任务
		iterprocess(x0,ai,bi,new_x); 
		MPI_Gather(new_x,n/size,MPI_DOUBLE,x0,n/size,MPI_DOUBLE,0,MPI_COMM_WORLD);
		if(rank==0){
			getresidul(a,x0,b,res);
			getl2norm(res,&rnorm);
			rnorm/=rnorm0;
			printf("第 %d 步rnorm的值: %lf\n",niter,rnorm);
		}
		MPI_Bcast(&rnorm,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
		niter=niter+1;
		//printf("%.20lf ",eps); 
	}while(rnorm>eps);
	
	if(rank==0){
		printf("原线性方程组的近似解为:\n");
		for(int i=0;i<n;i++)
			printf("x[%d]=%lf\n",i+1,x0[i]);
		printf("迭代次数:",niter);
	}
	MPI_Finalize();
}