jacobi.c

这是一段使用并行方式运算的jacobi迭代

#include <stdio.h>
#include "mpi.h"

#define  TOTALSIZE 16
#define  MYSIZE 4

int step = 10;
int n, myid, numprocs, i, j;
double a[TOTALSIZE][MYSIZE + 2] = {0};//a[16][6]
double b[TOTALSIZE][MYSIZE + 2] = {0};//b[16][6]
double c[TOTALSIZE] = {0};//receive
double d[TOTALSIZE] = {0};//send

int main(int argc, char *argv[])
{
	int begin_col, end_col;
	int i, j;
	MPI_Status status;

	//初始化并行环境
	MPI_Init(&argc, &argv);
	MPI_Comm_rank(MPI_COMM_WORLD, &myid);
	MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
	fprintf(stderr, "Process %d of %d is alive!! \n", myid, numprocs);

	//MPI_Barrier(MPI_COMM_WORLD);

	//初始化数组
	//这是一个16*16的稀疏矩阵
	//这里按照f的习惯先列后行！
	for (j = 0; j < (MYSIZE + 2); j++  )
	{
		for (i = 0; i < TOTALSIZE; i++  )
		{
			a[i][j] = 0;
		}
	}
	if (myid == 0)
	{
		for (i = 0; i < TOTALSIZE; i++ )
		{
			a[i][1] = 8;
		}
	}
	if (myid == 3)
	{
		for (i = 0; i < TOTALSIZE; i++ )
		{
			a[i][MYSIZE] = 8;
		}
	} 
	for (j= 0; j < (MYSIZE + 2); j++ )
	{
		a[0][j] = 8;
		a[TOTALSIZE - 1][j] = 8;
	}

	/*for (i = 0; i < (TOTALSIZE ); i++)
	{
		printf("%d ", myid);
		for (j = 0; j < (MYSIZE+2); j++)
		{
			printf("  (%4.2f)  ", a[i][j]);
		}
		printf("\n");
	} */  	



	for (n = 1; n <= step; n++)
	{

		if (myid < 3)
		{
			MPI_Recv(c,TOTALSIZE, MPI_DOUBLE, (myid + 1), 
					10, MPI_COMM_WORLD, &status);
			//printf("%d -r: ", myid);
			for (i  = 0; i < TOTALSIZE; i++)
			{
				a[i][MYSIZE + 1] = c[i];
				//printf(" %d ", c[i]);
			}
			//printf("\n");

		}

		if (myid > 0)
		{
			for (i  = 0; i < TOTALSIZE; i++)
			{
				d[i] = a[i][1];
			}
			MPI_Send(d, TOTALSIZE, MPI_DOUBLE, (myid - 1), 
					10, MPI_COMM_WORLD);
		}

		if (myid < 3)
		{
			for (i  = 0; i < TOTALSIZE; i++)
			{
				d[i] = a[i][MYSIZE];
			}
			MPI_Send(d, TOTALSIZE, MPI_DOUBLE, (myid + 1), 
				10, MPI_COMM_WORLD);
		}

		if (myid > 0)
		{
			MPI_Recv(c, TOTALSIZE, MPI_DOUBLE, (myid - 1), 
				10, MPI_COMM_WORLD, &status);
			for (i  = 0; i < TOTALSIZE; i++)
			{
				a[i][0] = c[i];
			}

		} 

		begin_col = 1;
		end_col = MYSIZE ;

		if (myid == 0)
		{
			begin_col = 2;
		}
		if (myid == 3)
		{
			end_col = MYSIZE - 1;
		}

		for (j = begin_col; j <=end_col; j++)
		{
			for (i = 1; i <=(TOTALSIZE - 1) ;i++ )
			{
				b[i][j] = ( a[i][j+1] + a[i][j - 1] + a[i+1][j] + a[i-1][j] )  * 0.25 ;
			}
		}

		for (j = begin_col; j <=end_col; j++)
		{
			for (i = 1; i <= (TOTALSIZE - 1) ;i++ )
			{
				a[i][j]  = b[i][j];
			}
		} 
	} 
	
	for (i = 1; i < (TOTALSIZE - 1); i++)
	{
		printf("%d: ", myid);
		for (j = begin_col; j <= end_col; j++)
		{
			printf("(%f)  ", a[i][j]);
		}
		printf("\n");
	}   
	MPI_Finalize();
	return 0;
}