/* matrix4.c - naive parallel transpose & dot product */

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

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

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,j;

   for(i=0;i<m;i++)
   {  for(j=0;j<n;j++)
       fprintf(fout,"%d ",M[i][j]);
      putc('\n',fout);
   }
}

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)
int ** M1;
int ** M2;
int m,n;
{  int i,j;
   int ** M3;
   MPI_Status status;

   M2 = transpose(M2,m,n);
   M3 = makeMatrix(m,m);
   for(i=0;i<m;i++)
   { for(j=0;j<m;j++)
     {  MPI_Send(M1[i],n,MPI_INT,j+1,0,MPI_COMM_WORLD);
        MPI_Send(M2[j],n,MPI_INT,j+1,0,MPI_COMM_WORLD);
     }
     for(j=0;j<m;j++)
      MPI_Recv(&(M3[i][j]),1,MPI_INT,j+1,0,MPI_COMM_WORLD,&status);
   }
   return M3;
}

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

   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);
      if(p!=m+1)
      {  printf("must have %d processors\n",m+1);
         MPI_Finalize();
         exit(0);
      }
      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);
      writeMatrix(stdout,m3,m,m);
      putchar('\n');
        
   }
   else
   {  int m,n,i;
      MPI_Status status;
      int *R,*C;
      int dp;

      MPI_Recv(&m,1,MPI_INT,0,0,MPI_COMM_WORLD,&status);
      if(p!=m+1)
      {  printf("must have %d processors\n",m+1);
         MPI_Finalize();
         exit(0);
      }
      MPI_Recv(&n,1,MPI_INT,0,0,MPI_COMM_WORLD,&status);

      R = makeVector(n);
      C = makeVector(n);

      for(i=0;i<m;i++)
      {  MPI_Recv(C,n,MPI_INT,0,0,MPI_COMM_WORLD,&status);
         MPI_Recv(R,n,MPI_INT,0,0,MPI_COMM_WORLD,&status);
         dp = dotProd(C,R,n);
         MPI_Send(&dp,1,MPI_INT,0,0,MPI_COMM_WORLD);
      }

   }
   MPI_Finalize();
}