mpieqsolver.cpp

一个使用MPI的并行Equation Solver。采用高斯迭代。在双核和4核机器上加速比可到2和4。并行的效果非常好。采用ghost rows。可操作于任意的process数

#include "mpi.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
int N=500;  

#define AR_INDEX(ar, a, b) (ar[((a) * (N+2)) + (b)])

int main(int argc, char** argv)
{
	double tol=0.075;
    int i, j, flag=0; 
	double diff = 0.0;
	double totaldiff=0.0;
    double* temp;
    int iteration = 0;
    double* grid1;
    double* grid2;
	clock_t start;
	clock_t finish;
	int rank, numtasks, namelen;
	int n;
	MPI_Status Stat1,Stat2;
	char processor_name[20];

	/* Allocate and initialize matrices */
	start=clock();
	MPI_Init(&argc,&argv);
	MPI_Comm_rank(MPI_COMM_WORLD,&rank);
	MPI_Comm_size(MPI_COMM_WORLD,&numtasks);
	MPI_Get_processor_name(processor_name,&namelen);
	fprintf(stderr,"Process %d on %s\n", rank, processor_name);
	fflush(stderr);

	n=N/numtasks;
	if(rank==numtasks-1)
		n=N/numtasks+N%numtasks;

	grid1 = (double*)malloc((n+2) * (N+2) * sizeof(double));
    grid2 = (double*)malloc((n+2) * (N+2) * sizeof(double));
	  
	for(i=0; i<=n+1; i++)
	{
		for(j=0; j<=N+1; j++)
		{
			if(i==0 || i==(n+1) || j==0 || j==(N+1)) {
				AR_INDEX(grid1, i, j) = 0.0;
				AR_INDEX(grid2, i, j) = 0.0;
			}
			else
				AR_INDEX(grid1, i, j) = 10.0*(i+rank*n) + j;
		}
	}
    
	while(flag==0){	
		iteration++;

		if(rank<numtasks-1)
			MPI_Send(&AR_INDEX(grid1,n,0), N+2, MPI_DOUBLE, rank+1, 0, MPI_COMM_WORLD);
		if(rank>0)
			MPI_Send(&AR_INDEX(grid1,1,0), N+2, MPI_DOUBLE, rank-1, 1, MPI_COMM_WORLD);

	    if(rank>0)
			MPI_Recv(&AR_INDEX(grid1,0,0), N+2, MPI_DOUBLE, rank-1, 0, MPI_COMM_WORLD, &Stat1);	
		if(rank<numtasks-1)
			MPI_Recv(&AR_INDEX(grid1,n+1,0), N+2, MPI_DOUBLE, rank+1, 1, MPI_COMM_WORLD, &Stat2);

        for (i = 1; i <= n; i++) { // row
            for (j = 1; j <= N; j++) { // column
                AR_INDEX(grid2, i, j) = 0.2 * (AR_INDEX(grid1, i, j)
                                               + AR_INDEX(grid1, i-1, j) + AR_INDEX(grid1, i+1, j)
                                               + AR_INDEX(grid1, i, j-1) + AR_INDEX(grid1, i, j+1));
                diff += fabs(AR_INDEX(grid2, i, j)  - AR_INDEX(grid1, i, j) );
            }
        }
		MPI_Reduce(&diff, &totaldiff, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
		if(rank==0){
		//	printf("Iteration %d: diff=%f\n", iteration, totaldiff/(N*N));
			if(totaldiff/(N*N) <= tol) {
				flag=1;	
			}
		}
	  temp = grid1;
      grid1 = grid2;
      grid2 = temp;
	  diff=0;
      MPI_Bcast(&flag, 1, MPI_INT, 0, MPI_COMM_WORLD);
    }; 

	free(grid1); /* free grid memory */
	free(grid2); /* free grid memory */
		
	finish=clock();
	
	if(rank==0){
		printf("Iteration count %d: final diff=%f\n", iteration, totaldiff/(N*N));
	    printf("The duration time is: %lf\n", (double)(finish-start)/CLOCKS_PER_SEC);
	}
	MPI_Finalize();
	//getchar();
}