/* matrix7.c - more rows than processors */

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

int * makeVector(n)
int n;
{  return (int *)malloc(n*sizeof(int));  }

writeVector(fout,v,n)
FILE * fout;
int * v;
int n;
{  int i;
   for(i=0;i<n;i++)
    fprintf(fout,"%d ",v[i]);
   putc('\n',fout);
}

showVector(v,n,m,p1,p2)
int * v;
int n;
char * m;
int p1,p2;
{  printf("%d %s %d: ",p1,m,p2);
   writeVector(stdout,v,n);
}

int ** makeMatrix(m,n)
int m,n;
{  int ** newm = (int **)malloc(m*sizeof(int *));
   int i;
   for(i=0;i<m;i++)
    newm[i] = makeVector(n);
   return newm;
}

readMatrix(fin,M,m,n)
FILE * fin;
int ** M;
int m,n;
{  int i,j;

   for(i=0;i<m;i++)
    for(j=0;j<n;j++)
     fscanf(fin,"%d",&(M[i][j]));
}

writeMatrix(fout,M,m,n)
FILE * fout;
int ** M;
int m,n;
{  int i;

   for(i=0;i<m;i++)
    writeVector(fout,M[i],n);
}

int ** transpose(M1,m,n)
int ** M1;
int m,n;
{  int ** M2;
   int i,j;

   M2 = makeMatrix(n,m);
   for(i=0;i<m;i++)
    for(j=0;j<n;j++)
     M2[j][i] = M1[i][j];
   return M2;
}

int dotProd(V1,V2,n)
int *V1,*V2;
int n;
{  int dp;
   int i;

   dp = 0;
   for(i=0;i<n;i++)
    dp = dp+V1[i]*V2[i];
   return dp;
}

int ** matrixProd(M1,M2,m,n,p)
int ** M1;
int ** M2;
int m,n,p;
{  int i,j;
   int r,s;
   int ** M3;
   MPI_Status status;

   M2 = transpose(M2,n,m);
   M3 = makeMatrix(m,m);

   s = m/(p-1);
   r = m%(p-1);

   for(i=0;i<r;i++)
    for(j=0;j<s+1;j++)
     MPI_Send(M1[(s+1)*i+j],n,MPI_INT,i+1,0,MPI_COMM_WORLD);
   for(i=0;i<p-r-1;i++)
    for(j=0;j<s;j++)
     MPI_Send(M1[(s+1)*r+s*i+j],n,MPI_INT,i+r+1,0,MPI_COMM_WORLD);

   for(i=0;i<m;i++)
    for(j=1;j<p;j++)
     MPI_Send(M2[i],n,MPI_INT,j,0,MPI_COMM_WORLD);

   for(i=0;i<r;i++)
    for(j=0;j<s+1;j++)
     MPI_Recv(M3[(s+1)*i+j],m,MPI_INT,i+1,0,MPI_COMM_WORLD,&status);
   for(i=0;i<p-r-1;i++)
    for(j=0;j<s;j++)
     MPI_Recv(M3[(s+1)*r+s*i+j],m,MPI_INT,i+r+1,0,MPI_COMM_WORLD,&status);

   return M3;
}

main(argc,argv)
int argc;
char ** argv;
{  int id,p;
   int s;

   MPI_Init(&argc,&argv);
   MPI_Comm_rank(MPI_COMM_WORLD,&id);
   MPI_Comm_size(MPI_COMM_WORLD,&p);

   if(id==0)
   {  FILE * fin;
      int ** m1;
      int ** m2;
      int ** m3;
      int m,n,i;

      fin = fopen(argv[1],"r");
      fscanf(fin,"%d %d",&m,&n);
      for(i=1;i<p;i++)
       MPI_Send(&m,1,MPI_INT,i,0,MPI_COMM_WORLD);
      for(i=1;i<p;i++)
       MPI_Send(&n,1,MPI_INT,i,0,MPI_COMM_WORLD);

      m1 = makeMatrix(m,n);
      readMatrix(fin,m1,m,n);
      m2 = makeMatrix(n,m);
      readMatrix(fin,m2,n,m);
      fclose(fin);
      m3 = matrixProd(m1,m2,m,n,p);
      writeMatrix(stdout,m3,m,m);
      putchar('\n');
        
   }
   else
   {  int m,n,i,j;
      MPI_Status status;
      int **R,*C,**RR;
      int r,s;

      MPI_Recv(&m,1,MPI_INT,0,0,MPI_COMM_WORLD,&status);
      MPI_Recv(&n,1,MPI_INT,0,0,MPI_COMM_WORLD,&status);

      s = m/(p-1);
      r = m%(p-1);
      if(id<=r)
       s = s+1;
      R = makeMatrix(s,n);
      C = makeVector(n);
      RR = makeMatrix(s,m);

      for(i=0;i<s;i++)
       MPI_Recv(R[i],n,MPI_INT,0,0,MPI_COMM_WORLD,&status);

      for(j=0;j<m;j++)
      {  MPI_Recv(C,n,MPI_INT,0,0,MPI_COMM_WORLD,&status);
         for(i=0;i<s;i++)
          RR[i][j] = dotProd(R[i],C,n);
      }
      
      for(i=0;i<s;i++)
       MPI_Send(RR[i],m,MPI_INT,0,0,MPI_COMM_WORLD);
   }
   MPI_Finalize();
}