mat_inv.c

经典的并行计算程序用于大规模矩阵求逆问题

/*************************************************
 *			 Parallel Computation 
 *			 for Matrix Inversion
 *			 Using MPICH2 Library
 *        by Wei WANG, ISEE, ZJU
 *		   Date: July 19, 2006 Summer
 *************************************************
 *          IMPORTANT: Input file format 
 *           M -- dimension of a matrix   
 *           i1  j1  A(i1,j1)             
 *           i2  j2  A(i2,j2)				
 *           i3  j3  A(i3,j3)				
 *************************************************															
//////////////////// EXAMPLE /////////////////////////// 4 /////////////////////////
// matrix:   1.4  0.0  5.0  3.0		    \                0 0 1.4		2 2 5.9	////
//////////   0.0  3.0  0.0  0.0    ===== >  input file:  0 2 5.0		2 3 1.1	////
//////////   5.0  0.0  5.9  1.1		    /                0 3 3.0		3 3 -1	////
//////////	 3.0  0.0  1.1   -1	   ///////////////////   1 1 3.0	////////////////
////////////////////////////////////////////////////////////////////////////////////														 
														 
 * It is recommended that put your app.exe and input data file in the root directory.
 */
														 
#include<stdio.h>
#include<stdlib.h>
#include"mpi.h"
#include<time.h>

#define a(x,y) a[x*M + y]		// a is a M-by-N matrix
#define A(x,y) A[x*M + y] 		// A is a M-by-N matrix
#define doublesize sizeof(double)
#define intsize sizeof(int)

long M, N, i, j;
double *A;
int my_rank;
int p;							// size of the communicator
FILE *mat_f = 0;				// pointer to input data file
MPI_Status status;

clock_t start, stop;
double duration;

/*============warnning message==================*/
void fatal( char* message )
{
	printf( "%s\n", message );
	exit(1);
}

/*========= to finalize temp matrix ============*/
void env_finalize( double* a, double* f )
{
	free( a );
	free( f );
}

/*========= to print a matrix with M-by-N ======*/
void print(double *A)
{
	for(i=0; i < M; i++)
	{
		for(j=0; j<N; j++)
			printf("%lf\t", A(i,j));
		printf("\n");
	}

	return;
}

/*== to store the inversed matrix in "output.txt" ==*/
int store_res(double *A)
{
	FILE *mat_out_f = 0;

	if((mat_out_f = fopen("\\output.txt", "w")) == NULL)
	{
		fatal("Write file open ERROR!");
	}

	fprintf(mat_out_f, "The matrix AFTER inversion\n");
	for(i=0; i<M;i++)
	{
		for(j=0;j<N;j++)
		{
			if(j < N-1)
				fprintf(mat_out_f, "%lf ", A(i,j)); 
			else
				fprintf(mat_out_f, "%lf\n", A(i,j));
		}
	}
	fclose(mat_out_f);	

	return 0;
}

/*================ main program =================*/
int main( int argc, char** argv )
{
	long k;
	int my_rank;
	int group_size;
	int i1, i2;
	long v, w;
	long m;
	double *a, *f;
	char ch;

	MPI_Init (&argc, &argv);							// parallel computation initialize
	MPI_Comm_size (MPI_COMM_WORLD, &group_size);		// get the number of units to join computation
	MPI_Comm_rank (MPI_COMM_WORLD, &my_rank);			// get current rank in communicator
	p = group_size;										// p is the number of units which join computation
	
	/* ------------------------------------------------------------------------------------------------- */
	/* -------------------------------------- INPUT DATA FROM FILE ------------------------------------- */
	
	if(my_rank == 0)
	{
		mat_f = fopen( "\\input.txt", "r" );			// open data file from root directory
		if(mat_f == NULL)								// exception dealing
		{
			fatal("Can not open input file!");
		}
	
		fscanf(mat_f, "%ld", &M);						// read dimension of matrix

		N = M;

		A = (double*) malloc ( doublesize * M * N );	// allocate memory space to store matrix

		for(i=0; i<M; i++)
			for(j=0; j<N; j++)
				A(i,j) = 0;								// each element to be zero in advance
				
		while( (ch=fgetc(mat_f)) != EOF )				// read data file
		{	
			fscanf(mat_f, "%ld", &i);
			fscanf(mat_f, "%ld", &j);
			fscanf(mat_f, "%lf", A + i*M + j);
			A(j,i) = A(i,j);							// a symmetrical matrix
		}

		fclose( mat_f );								// close the data file
		
//		printf("The matrix BEFORE inversion\n");		// print the matrix to be inversed
//		print( A );										// I dont recommend print it when the matrix is large:)
	}
	
	/* --------------------------------------------------------------------------------------------------*/
	/* --------------------------------------------------------------------------------------------------*/
	start = clock();									// start to time, just to have a test

	N = M;
		
	MPI_Bcast(&M, 1, MPI_LONG, 0, MPI_COMM_WORLD);		// processor '0' broadcast value of matrix's dimension

	m = M / p;											// number of rows each processor will receive
	if(M % p != 0)
		m++;

	f = (double*) malloc (doublesize * M);				// to create a buffer for send & receive data
	a = (double*) malloc (doublesize * m * M);

	if( a == NULL || f == NULL )						// exception dealing
		fatal("allocate buffer error!\nMaybe your PC is busy now or memory available is too small!\n");
	
	/* ---------------------------------------------------------------------------------------------------*/
	/* ---------------------------------- ASSIGN MISSION TO EACH PROCESSOR -------------------------------*/

	if(my_rank == 0)									// processor '0' divides matrix A into p small
	{													// m-by-M matrice. Then send each "m rows" to processor
		for(i = 0; i < m; i++)							// '1', '2' ... 'p-1'
			for(j = 0; j < M; j++)
				a(i, j) = A(i * p, j);

		for(i = 0; i < M; i++)
			if((i % p) != 0)
			{
				i1 = i % p;
				i2 = i/p + 1;
				MPI_Send(&A(i, 0), M, MPI_DOUBLE, i1, i2, MPI_COMM_WORLD);
			}
	}
	else												// if not processor '0', just receive data from
	{													// processor '0'
		for(i = 0; i < m; i++)
			MPI_Recv(&a(i, 0), M, MPI_DOUBLE, 0, i + 1, MPI_COMM_WORLD, &status);
	}
	
	/* -------------------------------------------------------------------------------------------------- */
	/* ----------------------------------- DO PARALLEL COMPUTATION -------------------------------------- */
	
	for(i = 0; i < m; i++)
	{
		for(j = 0; j < p; j++)
		{
			if(my_rank == j)							// if processor 'j', broadcast the elements in main row
			{
				v = i * p + j;
				if(a(i,v) != 0)////////////////////////////////////////////////////////////////////////////////
					a(i, v) = 1 / a(i, v);
				for(k = 0; k < M; k++)
					if(k != v)
						a(i, k) = a(i, k) * a(i, v);
				for(k = 0; k < M; k++)
					f[k] = a(i, k);
				MPI_Bcast(f, M, MPI_DOUBLE, my_rank, MPI_COMM_WORLD);
			}
			else
			{
				v = i * p + j;
				MPI_Bcast(f, M, MPI_DOUBLE, j, MPI_COMM_WORLD);
			}
			
			if(my_rank != j)							// if not processor 'j', use the main rows to do basic
			{											// conversion with other m rows.
				for(k = 0; k < m; k++)
					for(w = 0; w < M; w++)
						if(w != v)
							a(k, w) = a(k, w) - f[w] * a(k, v);
				for(k = 0; k < m; k++)
					a(k, v) = -f[v] * a(k, v);
			}
			
			if(my_rank == j)							// the processors, which send the main rows, do basic 
			{											// conversion with other m-1 rows which are not main rows
				for(k = 0; k<m; k++)
					if(k != i)
					{
						for(w = 0; w<M; w++)
							if(w != v)
								a(k, w) = a(k, w) - f[w] * a(k, v);
					}
				for(k = 0; k<m; k++)
					if(k != i)
						a(k,v) = -f[v] * a(k,v);
			}
		}
	}

	/* -------------------------------------------------------------------------------------------------- */
	/* --------------------------------	RECEIVE DATA, INVERSION COMPLETE -------------------------------- */

	if(my_rank == 0)									// if processor '0', fix the main m rows						
	{
		for(i=0; i<m; i++)
			for(j=0; j<M; j++)
				A(i * p, j) = a(i,j);
	}
	
	if(my_rank != 0)									// if not processor '0', send the inversed rows to '0'
	{
		for(i=0; i<m; i++)
			for(j=0; j<M; j++)
				MPI_Send(&a(i,j), 1, MPI_DOUBLE, 0, my_rank, MPI_COMM_WORLD);
	}
	else												// if processor '0', receive the inversed rows from processors
	{
		for(i=1; i<p; i++)
			for(j=0; j<m; j++)
				for(k=0; k<M; k++)
				{
					MPI_Recv(&a(j,k),1,MPI_DOUBLE,i,i,MPI_COMM_WORLD,&status);
				}
	}

	/* -------------------------------------------------------------------------------------------------- */
	/* -------------------------------- PRINT THE INVERSED MATRIX --------------------------------------- */
//	if(my_rank == 0)
//	{
//		printf("\n\nThe matrix AFTER inversion\n");
//		print( A );										// well, you can print it on screen, it would last a long time
//	}													
	
	/* -------------------------------------------------------------------------------------------------- */
	/* --------------------------------- STORE THE INVERSED MATRIX IN A FILE ---------------------------- */
	store_res( A );
		
	/* -------------------------------------------------------------------------------------------------- */
	/* ------------------------------------ END OF MPI PROCESS ------------------------------------------ */
	MPI_Finalize();
	env_finalize(a,f);									// free temporary storage space
	
	stop = clock();
	duration = ((double)(stop - start)) / CLK_TCK;
	printf("The duration time for this process: %lf seconds.\n", duration);
	
	free(A);
	
	return 0;
}