mpipar.c

在MPI上实现的矩阵相乘并行计算的源程序。

/*This is the parallel programm for the matrix multiplation */
/*programmer : Liu zhifeng 9500*/
#include "stdio.h"
#include "stdlib.h"
#include "mpi.h"
#define MAX 20
/* assume the maximum of the value of matrix*/
#define INTSIZE sizeof(int)
#define CHARSIZE sizeof(char)
#define A(x,y) A[x*K+y]
#define B(x,y) B[x*N+y]
#define a(x,y) a[x*K+y]
#define b(x,y) b[x*n+y]
#define c(l,x,y) c[x*(n*p)+y+l*n]
/*because the array A&B is inputed from outside files ,,it is hard
 to apply the array area   staticly ,so we need to apply array area
dynamicly.Both collum  and row of them maybe unknown,it is optimal
 to use pointee poiting to 1-demension array,but for the sake of
reading,we use macro to simulate 2-dimension array*/

/*use 1-dimension array to simulate 2-dimension*/
/*this is the parrallel program using pointee*/
/*-----------------------------------*/
/*  A :matrix  with row M,column K
    B :matrix  with row K,column N
    C = A * B
    C :matrix  with row M,column N
    Matrix A is divided by row
    Matrix B is divided by column
    just as the following:
    A= (
       a0
       a1
       .
       .
       .
       ap-1
       )
    B=(b0,b1,...,bp-1)
    C=A*B
     =(
       a0b0,  a0b1,  a0b2,........,a0bp-1
       a1b0,  a1b1,  a1b2,........,a1bp-1
         ..............................
         ..............................
       ap-1b0,ap-1b1,ap-1b2,......,ap-1bp-1
       )
     =
     (
     c0
     c1
     c2
     .
     .
     cp-1
     )
def: ci=(aib0,aib1,aib2,......aibp-1)
      p is number of processors.
     a is a part of matrix A
     b is a part of matrix B
     a: matrix  with row m,column K
     b: matrix  with row K,column n
   on the processor myid, a is in fact array amyid.
*/
int *A,*B;
int M,N,K ;
int m,n;
char fileC[40];/*store C int fileC*/ 
/*M is row-number of matrix A
  K is column-number of matrix A,row-number of matrix B
  N is column-number of matrix B
  m is row-number of matrix A
  K is column-number of matrix A,row-number of matrix B
  n  is column-number of matrix B
 */
int myid;
int p;
/*p is number of the processors
  inputed from keyboard */

FILE *fd;/*filepointer to fileC*/
MPI_Status status; /*MPI status*/
double time1,time2,time3;

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

void  read_A_and_B( )
{
char fileA[40],fileB[40];
 FILE *fdA,*fdB;
 int i,j;
 int RowB; /*Row of matrix B*/

printf("Processor %d:Read matrix A and B in\n",myid) ;
printf("First 2 integers of file should be the row and coloumn of the matrix\n");

printf("Filename of matrix A:");
/*scanf("%s",fileA);*/
strcpy(fileA,"liuA");
printf("Filename of matrix B:");
/*scanf("%s",fileB);*/
strcpy(fileB,"liuB");
time1=MPI_Wtime();
if ((fdA=fopen(fileA,"r"))==NULL)
  /*open the file for read only and set the file stream to the head of the file
*/
  {
  fatal("Open file A");
  }
  fread(&M,INTSIZE,1,fdA);
  fread(&K,INTSIZE,1,fdA);
  /*read M&K as row&column of matrix A*/


 if ((fdB=fopen(fileB,"r"))==NULL)
  /*open the file for read only and set the file stream to the head of the file*/
    {
    fclose(fdA);
    fatal("Open file B");/*free A & exit to the system*/
    }
  fread(&RowB,INTSIZE,1,fdB);
  /*RowB should be the same of K*/
  if (RowB!=K) 
  
{
fclose(fdA);fclose(fdB);
fatal("Column-number of A should equal to  Row-number of B,otherwise A CANNOT mutiple with B!");
}

  fread(&N,INTSIZE,1,fdB);
  /*read K&N as row&column of matrix B*/

m=M/p;if (M%p!=0) m++;  /*xiang4 shang4 qu3 zheng3 向上取整 */
n=N/p;if (N%p!=0) n++;


printf("Processor %d:M=%d,K=%d,N=%d,m=%d,n=%d\n",myid,M,K,N,m,n);

  A=(int *)malloc(INTSIZE*(m*p)*K);
 
  if (A==NULL) 
  {
  fclose(fdA);fclose(fdB);
  fatal("Dynamic allocation to matrix A fail!");
  }
  fread(A,INTSIZE,M*K,fdA);

 if (m*p!=M)
     for (i=M;i<m*p;i++)
       for (j=0;j<K;j++)
        A(i,j)=0; /*add some 0 to A*/
        /*for i j*/
   /*read matrix A in */

 
  B=(int *)malloc(INTSIZE*K*(n*p));
  if (B==NULL) 
 {
  fclose(fdA);fclose(fdB);
  fatal("Dynamic allocate space for matrix B fail!");
 }

 if (n*p==N)  fread(B,INTSIZE,K*N,fdB);
  else /*add some 0 to B*/
{
  for (i=0;i<K;i++)
   {
     fread(&( B(i,0) ),INTSIZE,N,fdB);
     for(j=N;j<n*p;j++)
        B(i,j)=0;
      /*for j*/
   } /*for i*/
    /*read matrix B in */
} /*else*/



fclose(fdA);
fclose(fdB);
/*close all the file pointer*/
free(fileA);
/*free the dynamilc allocated space*/

}  /*read_A_and_B()*/

void send_a_and_b( )
{
int i,j; /*loop variable*/



for (i=1;i<p;i++)
{

/*send(a_i,i)*/
MPI_Send(&(A(m*i,0) ),K*m,MPI_INT,i,i,MPI_COMM_WORLD);
/*send(b_i,i)*/
 for (j=0;j<K;j++)
  MPI_Send(&( B(j,n*i) ),n,MPI_INT,i,i,MPI_COMM_WORLD);
    /*for j*/

}/*for  i*/
printf("Processor %d:already send_a_and_b\n",myid);
} /*send_a_and_b*/

void get_a_and_b(int *a,int *b)
{
int i,j;/*loop variable*/

if (myid==0)
  /*get a and b by directly copying data from A and B*/
  {
   for (i=0;i<m;i++)
    for (j=0;j<K;j++)
       a(i,j)=A(i,j);
   for (i=0;i<K;i++)
    for (j=0;j<n;j++)
     b(i,j)=B(i,j);
  } /*if*/
else /*get a and b from data sent by Process 0 */
{
/*recv(a,0)*/
MPI_Recv(a,K*m,MPI_INT,0,myid,MPI_COMM_WORLD,&status);

/*recv(b,0)*/
for (j=0;j<K;j++)
 MPI_Recv(&( b(j,0) ),n,MPI_INT,0,myid,MPI_COMM_WORLD,&status);
}
printf("Processor %d:already get_a_and_b\n",myid);


} /*get_a_and_b*/





void matrix_multiplation(int l,int *a,int *b,int *c)
/* a*b-->c[l] */
{ int i,j,s;
  for (i=0;i<m;i++)
    for (j=0;j<n;j++)
      for (c(l,i,j)=0,s=0;s<K;s++)
          c(l,i,j)+=a(i,s)*b(s,j);
} /* matrix_multiplation*/



void parallel_computation(int *a,int *b,int *c)
{
int i,j;/*loop variable*/
int l;
int *buffer;
int mm1,mp1; /*its usage will be declared later*/
for (i=0;i<p;i++)
{
 l=(i+myid)%p;

/*Now compute a_myid*b_l*/
matrix_multiplation(l,a,b,c) ;
   /* multiplate matrix A with matrix B,and save the result into matrix C[l]*/
mm1=(p+myid-1)%p;
/*Process mm1 will use b_myid next */
mp1=(myid+1)%p;
/*Process myid will use b_mp1 next*/

  if (i!=p-1)
   /*send(B,mm1), recv(B,mp1)*/
  {
   if(myid%2==0)/*myid is even,first send,then recv */
  {
    MPI_Send(b,K*n,MPI_INT,mm1,mm1,MPI_COMM_WORLD);
    MPI_Recv(b,K*n,MPI_INT,mp1,myid,MPI_COMM_WORLD,&status);
    } /*if*/
   else
   /*myid is odd,recv before send*/
   { 
   if( (buffer=(int *)malloc(K*n*INTSIZE) )==NULL)
     fatal("error in allocating buffer!");
   for(j=0;j<K*n;j++) buffer[j]=b[j];
   /*store b in buffer to avoid it cleared by message sending from others*/
   MPI_Recv(b,K*n,MPI_INT,mp1,myid,MPI_COMM_WORLD,&status);
   MPI_Send(buffer,K*n,MPI_INT,mm1,mm1,MPI_COMM_WORLD);
    free(buffer);
   }
}  /*if*/
printf("Processor %d:Complete compute  i=%d,l=%d\n",myid,i,l);
}/*for  i */
printf("Processor %d:All parallel computation Complete\n",myid);
} /*parralel_computation*/

void output(int *c)
{
/*int i,j;
printf("Processor%d:\n",myid);
for (i=0;i<m;i++)
  {for(j=0;j<N;j++)
   printf("%d ", *(c+i*N+j)    );
   printf("\n");
  }for i*/
}/*output*/


void write_c_to_file(int i,int *c)
/*write c(i) to the outside file*/
{

int RowC,ColumnC;

if (i==0)
{
  strcpy(fileC,"liuC1");
  fd=fopen(fileC,"w");
  /*open a file for write only,if it already exists,truncate it to zero length*/
   if (fd==NULL) 
   {
  
   fatal("Error in open file C");
   }/*if*/

  RowC=   p*m;
  ColumnC=p*n;
  

  fwrite(&RowC,INTSIZE,1,fd);

  fwrite(&ColumnC,INTSIZE,1,fd);


}/*i=0,write the beginning of the file*/
 
 fwrite(c,INTSIZE,(p*m)*n,fd);

if (i==p-1)
  {fclose(fd);}

}/*write c(i) to the file*/





void Environment_Finalize(int *a, int *b,int *c)
/*before exiting system,free dynamilc-allocated space*/
{
 free(a);
 free(b);
 free(c);
 if (myid==0)
 { free(A);
 free(B);
/*free dynamil-allocated space*/
}
}/*Environment_Finalize*/



/*---------------- main ---------------------*/
int main(int argc, char **argv)
{
 int i,j,l,group_size,mp1,mm1;
 /* i: loop variable
    l: is  the current column Pmyid should  calculate
 myid: the number of cuurrent processor
  group_size:number of process
 */
int *a;  /*matrix Amyid*/
int *b;/* matrix Bmyid*/
int *c; /*cmyid with Row m,column N*/
double starttime,endtime;

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,&myid);
/*get the number standing for the current processor*/

  if(myid==0) starttime = MPI_Wtime(); 
    
p=group_size;
printf("Processor %d:p=%d\n",myid,p);

/* read matrix A&B from outside files*/
if (myid==0)
{  read_A_and_B( );
  /*get matrix A and B from outside files*/
   for(i=1;i<p;i++)
      /*send M,N & K to Processor i*/
  {
   MPI_Send(&M,1,MPI_INT,i,i,MPI_COMM_WORLD);
   MPI_Send(&K,1,MPI_INT,i,i,MPI_COMM_WORLD);
   MPI_Send(&N,1,MPI_INT,i,i,MPI_COMM_WORLD);
  }/*for i*/

} /*if*/
else /*myid<>=0*/
/*wait for Processor 0 to read matrix A&B in and send M,N,K*/
{
  MPI_Recv(&M,1,MPI_INT,0,myid,MPI_COMM_WORLD,&status);
  MPI_Recv(&K,1,MPI_INT,0,myid,MPI_COMM_WORLD,&status);
  MPI_Recv(&N,1,MPI_INT,0,myid,MPI_COMM_WORLD,&status);
 
 m=M/p;if (M%p!=0) m++;  /*xiang4 shang4 qu3 zheng3 向上取整 */
 n=N/p;if (N%p!=0) n++;

 

}/*else*/




a=(int *)malloc(INTSIZE*m*K);
b=(int *)malloc(INTSIZE*K*n);
c=(int *)malloc(INTSIZE*m*N);
if (a==NULL||b==NULL||c==NULL) fatal("Allocate space for a,b or c fail!");


 if (myid==0) 
{get_a_and_b(a,b); 
send_a_and_b(   );
}
else
 get_a_and_b(a,b);

if (myid==0) time2=MPI_Wtime();

/* parallel computation*/
parallel_computation(a,b,c);
/*output(c);*/

if (myid==0)
{
 time3=MPI_Wtime();
write_c_to_file(0,c);

/*write c(0) to outside file*/
  for(i=1;i<p;i++)
   {

   MPI_Recv(c,(p*m)*n,MPI_INT,i,i,MPI_COMM_WORLD,&status);
   write_c_to_file(i,c);
   printf("Processor%d:already  receive c[%d]\n",myid,i);
   }/*for i*/
}/*if*/
else
/*send(c_myid,0)*/
{

MPI_Send(c,(p*m)*n,MPI_INT,0,myid,MPI_COMM_WORLD);
printf("Processor%d:already send c[%d] to 0\n",myid,myid);
}
if(myid==0){
           endtime=MPI_Wtime();
           printf("Parallel whole running time=%f secends\n",endtime-starttime);
           printf("Time to distribute data:%f secends\n",time2-time1);
           printf("Time to parallel compute:%f secends\n",time3-time2);
           printf("Time to write output file:%f secends\n",endtime-time3); 
           }
MPI_Finalize();
Environment_Finalize(a,b,c);
return (0);
}
/*main*/