sorp.c

Parallel programming／Lou Baker, Bradley J.Smith ．—New York：McGraw-Hill Book Co.

/**********************************************************

	SPMD Poisson solver using MPI Message-Passing (WinMPI).
	Uses Succesive Over-Relaxation (SOR) Gauss-Seidel.
	L. Baker, Dagonet Software

***********************************************************/


/*#define WIN31*/
#include <stdio.h>
#include <stdlib.h>
/*below stmt for gcc compiler-non windows*/
/*
#define far
#define pascal
*/
#include "mpi.h"
#define DEBUG

#ifndef WIN31
#define Int int
#define LPPSTR char **
#endif

#define Index(i,j) [(j)+(i)*ny]
/*#define max(a,b) ((a)>(b)?(a):(b))*/
Int nproc,myid,nmax,ncol,nrow;
/*
   ncol,nrow= number of columns and rows of processors
nx,ny= number of grid lines on this processor
nmax= maximum of nx and ny
*/

struct neighborlist{
	Int procid,bindex;
	struct neighborlist *next;
	} *neighbors,*nn;/* send/recv to/from */

MPI_Status status;

/*double Iterate(Int nx,Int ny,float *u,float *b,float *rho,double omega);*/

double Iterate(Int nx,Int ny,float *u,float *b,float *rho,double omega)
{
Int i,j;
double om,resid,residg;
om=omega*.25;/* for 2-d*/
/*send interface data to neighbors*/
/* loop over neighbors, sending appropriate bndry of u array */
nn=neighbors;
while(nn){
	/*send to procid b Index(bindexs,...)*/
	switch(nn->bindex){
		case 0: for(i=0;i<nx;i++)b Index(nn->bindex,i)=u Index(i,0);break;
		case 1: for(i=0;i<nx;i++)b Index(nn->bindex,i)=u Index(i,ny-1);break;
		case 2: for(i=0;i<ny;i++)b Index(nn->bindex,i)=u Index(0,i);break;
		case 3: for(i=0;i<ny;i++)b Index(nn->bindex,i)=u Index(nx-1,i);break;
    }
#ifdef DEBUG
if(!myid)printf(" root Sending\n");
#endif 
	MPI_Send(&(b Index(nn->bindex,0)),nmax,MPI_FLOAT,nn->procid,1,MPI_COMM_WORLD);
	nn=nn->next;
	}
/* interior nodes*/
residg=0.;
for(i=1;i<nx-1;i++)
	for(j=1;j<ny-1;j++)
		{resid= rho Index(i,j) + u Index(i+1,j)+u Index(i-1,j)
				+ u Index(i,j+1)+ u Index(i,j-1)-4.* u Index(i,j);
		 /* immediate update= Gauss-Jordan */
		 u Index(i,j) +=(float)( om *resid);
		 residg+=resid*resid;
		 }
/*receive interface data*/
/* put into b array */
nn=neighbors;
while(nn){
	/*recv from procid b Index(bindexs,0...<ny)*/
	/*BLOCKING recv*/         
#ifdef DEBUG
if(!myid)printf(" root at Recv\n");
#endif 
	MPI_Recv(&(b Index(nn->bindex,0)),nmax,MPI_FLOAT,nn->procid,1,MPI_COMM_WORLD,&status);
	nn=nn->next;              
#ifdef DEBUG
if(!myid)printf(" root after Recv\n");
#endif 
	
    }

/* apply interface NEWS 2-d grid assumed*/
j=0;
for(i=1;i<nx-1;i++)
		{resid= rho Index(i,j) + u Index(i+1,j)+u Index(i-1,j)
				+ u Index(i,j+1)+ b Index(0,i)-4.* u Index(i,j);
		 /* immediate update= Gauss-Jordan */
		 u Index(i,j) += (float)(om *resid);
		 residg+=resid*resid;
         }
j=ny-1;
for(i=1;i<nx-1;i++)
		{resid= rho Index(i,j) + u Index(i+1,j)+u Index(i-1,j)
				+ u Index(i,j-1)+ b Index(1,i)-4.* u Index(i,j);
		 /* immediate update= Gauss-Jordan */
		 u Index(i,j) += (float)(om *resid);
		 residg+=resid*resid;
         }
i=0;
for(j=1;j<ny-1;j++)
		{resid= rho Index(i,j) + u Index(i+1,j)+u Index(i,j-1)
				+ u Index(i,j+1)+ b Index(2,j)-4.* u Index(i,j);
		 /* immediate update= Gauss-Jordan */
		 u Index(i,j) += (float)(om *resid);
		 residg+=resid*resid;
         }
i=nx-1;
for(j=1;j<ny-1;j++)
		{resid= rho Index(i,j) + u Index(i-1,j)+u Index(i,j-1)
				+ u Index(i,j+1)+ b Index(3,j)-4.* u Index(i,j);
		 /* immediate update= Gauss-Jordan */
		 u Index(i,j) += (float)(om *resid);
		 residg+=resid*resid;
         }
/* four corners: each with 2 b terms*/
		{resid= rho Index(0,0) + u Index(1,0)+b Index(0,0)
				+ u Index(0,1)+ b Index(2,0)-4.* u Index(0,0);
		 u Index(0,0) += (float)(om *resid);
		 residg+=resid*resid;
         }

		{resid= rho Index(0,ny-1) + u Index(1,ny-1)+b Index(1,0)
				+ u Index(0,ny-2)+ b Index(2,ny-1)-4.* u Index(0,ny-1);
		 u Index(0,ny-1) += (float)(om *resid);
		 residg+=resid*resid;
         }

		{resid= rho Index(nx-1,0) + u Index(nx-2,0)+b Index(0,nx-1)
				+ u Index(nx-1,1)+ b Index(3,0)-4.* u Index(nx-1,0);
		 u Index(nx-1,0) += (float)(om *resid);
		 residg+=resid*resid;
         }

        {resid= rho Index(nx-1,ny-1) + u Index(nx-2,ny-1)+b Index(1,nx-1)
                + u Index(nx-1,ny-2)+ b Index(3,ny-1)-4.* u Index(nx-1,ny-1);
         u Index(nx-1,ny-1) += (float)(om *resid);
         residg+=resid*resid;
         }
return residg;
}


#ifdef WIN31
int MPI_main (int argc, LPPSTR argv)
#else
int main (int argc, char *argv[])
#endif
/*
int MPI_main (int argc, LPPSTR argv)
*/
{
	double residual,localresid,abserr,/*relerr,*/omega,start,fini;
	float *u,*b,*rho;
	Int iterno,nx,ny,nxtot,nytot,i,j,k,ibtm,jbtm,myrow,mycol;
	Int sendbuf[2],recbuf[2];
#ifdef WIN31
	LPPSTR far *argmnt;
#else 
	char **argmnt[];
#endif	
	/*init*/
#ifdef DEBUG
if(!myid)printf(" before MPI calls\n");
#endif	

	MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD,&myid);
    MPI_Comm_size(MPI_COMM_WORLD,&nproc);
#ifdef DEBUG
if(!myid)printf(" after MPI calls\n");
#endif
    /*define pblm */

	nxtot=nytot=10;/*full grid size nytot>=nxtot*/
	switch(nproc){
		case 1: nx=nxtot;ny=nytot;ncol=nrow=1;break;
		case 2: nx=nxtot/2;ny=nytot;ncol=2;nrow=1;break;
		case 4: nx=nxtot/2;ny=nytot/2;ncol=nrow=2;break;
		default:
			printf(" allowed proc number 1,2,4\n");
			exit(1);
    }
	nmax= max(nx,ny);
	abserr=1.e-6;
	abserr=1.e-1;
	/*relerr=1.;*/
	omega=1.6;/*relaxation parameter. Should be >=1. */

	/*printf(" enter omega, abserr\n");scanf("%lf%lf",&omega,&abserr);
	printf(" echo omega, abserr %lf %lf\n",omega, abserr);
	*/
	if(nproc==1)myid=0;
#ifdef DEBUG
//if(!myid)printf(" before MPI timer\n");
#endif

	if(!myid)start=MPI_Wtime();
#ifdef DEBUG
//if(!myid)printf(" after MPI timer\n");
#endif
    u= (float *)malloc(sizeof(float)*((int)(nx*ny)));/*fill with 0 */
	rho= (float *)malloc(sizeof(float)*((int)(nx*ny)));/* fill with 0*/
	for(i=0;i<nx;i++)
		for(j=0;j<ny;j++){u Index(i,j)= rho Index(i,j) =(float)0.0;}
	b= (float *)malloc(sizeof(float)*4*(int)(ny));/* fill with x*y*/
#ifdef DEBUG
if(!myid)printf(" after malloc\n");
#endif
    if(!myid){
		start=MPI_Wtime();
#ifdef DEBUG
if(!myid)printf(" after MPI_Wtime\n");
#endif
		/* assign processors their part of the problem */
		if(nproc==1){ibtm=jbtm=0;}
		else{
			for(i=0;i<ncol;i++){
				ibtm=i*nx;
				for(j=0;j<nrow;j++){
					if(!i&&!j)continue;/*self*/
					jbtm=j*ny;
					/* pack ibtm,jbtm into buffer and send*/
					sendbuf[0]=ibtm;
					sendbuf[1]=jbtm;
					MPI_Send(sendbuf,(Int)2,MPI_INT,(Int)(i+j*ncol),(Int)0,MPI_COMM_WORLD);
#ifdef DEBUG
if(!myid)printf(" prc=%ld i,jbtm=%ld %ld\n",i+j*ncol,ibtm,jbtm);
#endif
                    }
				}
			}
	ibtm=jbtm=0;/*me*/
	}/* root process*/
	else{
		/*receive my part of pblm: ibtm,jbtm*/
		MPI_Recv(recbuf,2,MPI_INT,0,0,MPI_COMM_WORLD,&status);
		/*ibtm, jbtm from root*/
		ibtm=recbuf[0];
		jbtm=recbuf[1];
	}/* slave*/
	neighbors=NULL;
	/* determine who my neighbors are and what portion of the b array
	   they send/recv */
	for(i=0;i<ncol;i++)
		for(j=0;j<nrow;j++){
			k=i+j*ncol;
			if(k==myid)
				{
				myrow=j;
				mycol=i;
				break;
				}
		}
	for(i=0;i<ncol;i++)
		for(j=0;j<nrow;j++){
			k=i+j*ncol;
			if(k==myid)continue; /*self*/
			if(i==mycol && abs(j-myrow)==1)
				{
				if(j>myrow){
					/* other guy is to my top(north) */
					nn= /*(struct *neighborlist)*/malloc(sizeof(struct neighborlist));
                    nn->next=neighbors;
                    nn->procid=k;
					nn->bindex=1; 
					neighbors=nn;
#ifdef DEBUG
if(!myid)printf(" root has neighbor to north\n");
#endif 
                    }
                else{
					/* other guy is to my btm(south) */
					nn= /*(struct *neighborlist)*/malloc(sizeof(struct neighborlist));
                    nn->next=neighbors;
                    nn->procid=k;
					nn->bindex=0;          
					neighbors=nn;
#ifdef DEBUG
if(!myid)printf(" root has neighbor to south\n");
#endif 
					
                    }
            }
			if(j==myrow && abs(i-mycol)==1)
				{
				if(i>mycol){
					/* other guy is to my right(east) */
					nn= /*(struct *neighborlist)*/malloc(sizeof(struct neighborlist));
                    nn->next=neighbors;
                    nn->procid=k;
					nn->bindex=3;
					neighbors=nn;
#ifdef DEBUG
if(!myid)printf(" root has neighbor to east\n");
#endif 
                    }
                else{
                    /* other guy is to my left(west) */
					nn= /*(struct *neighborlist)*/malloc(sizeof(struct neighborlist));
                    nn->next=neighbors;
                    nn->procid=k;
					nn->bindex=2;
					neighbors=nn;
#ifdef DEBUG
if(!myid)printf(" root has neighbor to west\n");
#endif 
                    }
			}
        }
	/*initialize*/
	for(i=0;i<nx;i++)
		{b Index(0,i)=(float)0.0;
		 b Index(1,i)=(float) (i+ibtm);
		 }
	for(j=0;j<ny;j++)
		{b Index(2,j)=(float)0.0;
		 b Index(3,j)=(float) (j+jbtm);
         }
	/*calculate*/
	if(!myid){
		iterno=0;
		do{
#ifdef DEBUG
if(!myid)printf(" root calling Iterate\n");
#endif 
		    /*broadcast to all to iterate*/             
		    k=0;
		    MPI_Bcast(&k,1,MPI_INT,0,MPI_COMM_WORLD);
			localresid=Iterate(nx,ny,u,b,rho,omega);
#ifdef DEBUG
if(!myid)printf(" root returned Iterate\n");
#endif 
			iterno++;
			if(nproc==1)residual=localresid;
			else
				{MPI_Reduce(&localresid,&residual,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);}
		printf(" iterno=%ld residual=%lf\n",iterno,residual);
		if((residual<abserr)){
		    k=1;
		    MPI_Bcast(&k,1,MPI_INT,0,MPI_COMM_WORLD);
			break;            
			}
		}while (1);
    }
	else
		{/*slave*/
			do{
			MPI_Bcast(&k,1,MPI_INT,0,MPI_COMM_WORLD);
			if(k)break;/* done if k nonzero*/
			localresid=Iterate(nx,ny,u,b,rho,omega);
			MPI_Reduce(&localresid,&residual,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
			}while(1);
		}	
	/* print solution*/
   if(!myid)
		{fini=MPI_Wtime();
		for(i=0;i<nx;i++)
			for(j=0;j<ny;j++){
				printf(" i=%ld,j=%ld u=%lf\n",i,j, u Index(i,j) );
			}
			for(k=1;k<nproc;k++){
			MPI_Recv(u,(Int)(nx*ny),MPI_FLOAT,(Int)k,2,MPI_COMM_WORLD,&status);
			for(i=0;i<nx;i++)
				for(j=0;j<ny;j++){
					printf(" i=%ld,j=%ld u=%lf\n",i+nx*(k%ncol),j+ny*(k/ncol), u Index(i,j) );
				}

			}
		printf(" number of iterations=%ld time=%lf\n",iterno,fini-start);
		}
   else{
		MPI_Send(u,nx*ny,MPI_FLOAT,0,2,MPI_COMM_WORLD);
	}
MPI_Finalize();
return 0;
}