paraw.c

在MPI上用并行LU分解求解线性方程组的源程序。

/*This is the parallel programm for the matrix multiplation */
#include "stdio.h"
#include "stdlib.h"
#include "mpi.h"

/* assume the maximum of the value of matrix*/
#define a(x,y) a[x*M+y]
#define b(x) b[x]
#define HOSTNAME "pine"
#define DIM 50

int port;
int M ;
int size;
int m;
float B;
float value;
double time1;
double time2;
float A[DIM];
float *B1;
/*char *fileX;*//*store X int fileX*/
/*M is size of matrix A and B
  m is size of matrix A and B in each point
 */
int my_rank;
int p;
/*p is number of the processors
  inputed from keyboard */

MPI_Status status; /*MPI status*/

void fatal(char *message)
/*print error message and exit system*/
{
printf("%s\n",message);
exit(0);/*exit system*/
}




void Environment_Finalize(float *a,float *b,float *x)
/*before exiting system,free dynamilc-allocated space*/
{
 free(a);
 free(b);
 free(x);
 
}/*Environment_Finalize*/



/*---------------- main ---------------------*/
int main(int argc, char **argv)
{
 
 int i,j,k,l,my_rank,group_size,mp1,mm1;
 /* i,j,k: loop variable
    l: is  the  column of matrix
    my_rank: the number of cuurrent processor
  group_size:number of process
 */
int i1,i2;
int icol,jcol;
float s; /*exchange register*/
float d; /*middle add register*/
float *f;
float maxa;
float *Afinal,*Bfinal;
float *a;  /*matrix Amy_rank*/
float *b;  /*matrix Bmy_rank*/
float *x;  /*matrix Xmy_rank*/
float nb;
char *fileA,*fileB;
char buffer[100];
FILE *fdA,*fdB,*fdB1;
int sizeB; /*size of matrix B*/



MPI_Init(&argc,&argv);/*MPI initialization*/
MPI_Comm_size(MPI_COMM_WORLD,&group_size);/*size of the process group*/
MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
p=group_size;



 /*get the number standing for the current processor*/

/* read matrix A&B from outside files*/


/*--------------------read A and B--------------------*/


if(my_rank==0)
{

 fileA="a";
 fdA=fopen(fileA,"r");
 fileB="b";
 fdB=fopen(fileB,"r");
 fdB1=fopen(fileB,"r");
if ((fdA=fopen(fileA,"r"))==NULL) printf("reada error");
	  /*open the file for read only and set the file stream to the head of the file
		*/
/* fatal("Open file A");
  fread(&M,sizeof(int),1,fdA) ;*/
  /*read Size as row&column of matrix A*/


 if ((fdB=fopen(fileB,"r"))==NULL)  printf("readb error");
  /*open the file for read only and set the file stream to the head of the file*/
/*  fread(&sizeB,sizeof(int),1,fdB);*/
  /*SizeB should be the same of M*/
 /* if (sizeB!=M) fatal("Column-number of A should equal to  Row-number of B,otherwise A CANNOT mutiple with B!");
 *//*read size of matrix B*/

 if ((fdB=fopen(fileB,"r"))==NULL)  printf("readb error");

  M=DIM;

  B1=(float *)malloc(sizeof(float)*M);
  if (B1==NULL) fatal("allocate space A fail");

}





 /*------------------read A and B-----------------*/

MPI_Bcast(&M,1,MPI_INT,0,MPI_COMM_WORLD);

m=M/p;if (M%p!=0) m++;/*the space allocated to each point*/
f=(float *)malloc(sizeof(float)*M);
a=(float *)malloc(sizeof(float)*m*M);
b=(float *)malloc(sizeof(float)*m);
if (a==NULL||b==NULL||f==NULL) fatal("allocate space a or b or f fail!");



if (my_rank==0)
     {
	
	printf("\n the coefficient matrix \n");
	for(i=0;i<M;i++)
	       
		       if ((i%p)!=0)
			      {
				i1=i%p;
				i2=i/p+1;
				fread(A,sizeof(float),M,fdA);
				for (j=0;j<M;j++)
					printf("   %f",A[j]);
				printf("\n");
				MPI_Send(A,M,MPI_FLOAT,i1,i2,MPI_COMM_WORLD);
				fread(&B,sizeof(float),1,fdB);
                                MPI_Send(&B,1,MPI_FLOAT,i1,i2,MPI_COMM_WORLD);
                              }

			else
			      {
			
				fread(A,sizeof(float),M,fdA);

				for(j=0;j<M;j++)
                 			{ printf("   %f",A[j]);
					  a[(i/p)*M+j]=A[j];
					}
				printf("\n");
				fread(&B,sizeof(float),1,fdB);
				b[i/p]=B;
			      }

	fread(B1,sizeof(float),M,fdB);
	printf("\nthe result matrix\n");
	for (i=0;i<M;i++)
		printf(" b[%d]=%f\n",i,B1[i]);

        free(B1);
     }
else
     {
	for(i=0;i<m;i++)
	     if ((i*p+my_rank)<=(M-1))
		{
		  MPI_Recv(&a[i*M],M,MPI_FLOAT,0,i+1,MPI_COMM_WORLD,&status);
                  MPI_Recv(&b(i),1,MPI_FLOAT,0,i+1,MPI_COMM_WORLD,&status);
		}

     }



time1=MPI_Wtime();
icol=0;
	for (j=0;j<=M-2;j++)
	{
	   if (my_rank==(j%p))
	       {
		maxa=a(icol,j);
		l=j;
		for (i=j;i<=M-1;i++) if (a(icol,i)>maxa)
				{maxa=a(icol,i);l=i;}
		if (l!=j)
		     for(i=0;i<m;i++)
			{s=a(i,j);a(i,j)=a(i,l);a(i,l)=s;}
		if (a(icol,j)==0)
			{printf("A is wrong");
			 return(0);
			}
		for(i=j+1;i<=M-1;i++)
			a(icol,i)=a(icol,i)/a(icol,j);
		b(icol)=b(icol)/a(icol,j);
		for (i=j+1;i<=M-1;i++)
			f[i-j-1]=a(icol,i);
		f[M-1]=b[icol];
		
		MPI_Send(&l,1,MPI_FLOAT,(my_rank+1)%p,1,MPI_COMM_WORLD);
		MPI_Send(&f[0],M,MPI_FLOAT,(my_rank+1)%p,1,MPI_COMM_WORLD);
		icol=icol+1;
	       }
	    else
	       {
	        MPI_Recv(&l,1,MPI_INT,(p+my_rank-1)%p,1,MPI_COMM_WORLD,&status);
		MPI_Recv(&f[0],M,MPI_FLOAT,(p+my_rank-1)%p,1,MPI_COMM_WORLD,&status);
		if (((my_rank+1)%p)!=(j%p))
			{
			 MPI_Send(&l,1,MPI_FLOAT,(my_rank+1)%p,1,MPI_COMM_WORLD);
			 MPI_Send(&f[0
],M,MPI_FLOAT,(my_rank+1)%p,1,MPI_COMM_WORLD);
			}
	       }
	    if (my_rank!=(j%p))
	     if (l!=j)
	       for(i=0;i<=m-1;i++)
		    {
		     s=a(i,j);a(i,j)=a(i,l);a(i,l)=s;}

	    
	    for (k=icol;k<=m-1;k++)
		{    for (i=j+1;i<=M-1;i++)
			   {a(k,i)=a(k,i)-f[i-j-1]*a(k,j);}
		     b(k)=b(k)-f[M-1]*a(k,j);
		}/*for*/

     } /*for j finish computing matrix*/






Afinal=(float *)malloc(sizeof(float)*m*p*M);
Bfinal=(float *)malloc(sizeof(float)*m*p);
x=(float *)malloc(sizeof(float)*M);



if (my_rank!=0)
    {
     
     for (i=0;i<m;i++)
	 {	  MPI_Send(&a(i,0),M,MPI_FLOAT,0,i+1,MPI_COMM_WORLD);
	  MPI_Send(&b(i),1,MPI_FLOAT,0,i+1,MPI_COMM_WORLD);
	 }
    }
else
    {
     
     for(i=0;i<m;i++)
	{
	  Bfinal[i*p]=b(i);
	  for(k=0;k<M;k++)
		Afinal[i*p*M+k]=a(i,k);
	  for (j=1;j<p;j++)
		{
		 MPI_Recv(&Afinal[(i*p+j)*M],M,MPI_FLOAT,j,i+1,MPI_COMM_WORLD,&status);
		 MPI_Recv(&Bfinal[i*p+j],1,MPI_FLOAT,j,i+1,MPI_COMM_WORLD,&status);
		}
	}
     }
 



if (my_rank==0)
     {
	  if(abs(Afinal[M*M-1])==0)
		{printf("wrong");return(0);}/*wrong*/
	  x[M-1]=Bfinal[M-1]/Afinal[M*M-1];
	  for(i=M-2;i>=0;i--)
	     {d=0;
	      for(j=i+1;j<M;j++) d=d+Afinal[i*M+j]*x[j];
	      x[i]=Bfinal[i]-d;
	     }
    }

free(f);



time2=MPI_Wtime();
time2=time2-time1;
if (my_rank==0)
	printf("the result matrix is \n");
if(my_rank==0)
	for(i=0;i<M;i++)
		printf("x[%d]=%f\n",i,x[i]);

if (my_rank==0)
	printf("\nrunning time=%f\n",time2);


MPI_Barrier(MPI_COMM_WORLD);
MPI_Finalize();

Environment_Finalize(a,b,x);
return (0);
/*Pay attention to TAG of MPI*/
}
/*main*/