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

#define DIM 8

int printMat(float *mat, int numRow, int numCol, int mype, int nproc, char *str){

    float *buf;
    int i, j, k;
    MPI_Status status;

    if (mype == 0)
    {   
		printf("%s\n",str);
		for (i=0; i<numRow; i++)
		{
			for (j=0; j<numCol; j++)
			printf("%.2f\t", mat[(i*numCol)+j]);
			printf("\n");
		} 

		buf = (float *)malloc(numCol*numRow*sizeof(float));

		for (k=1; k < nproc; k++)
		{
			MPI_Recv(buf, numRow*numCol, MPI_FLOAT, k, 0, MPI_COMM_WORLD, &status);
			for (i=0; i < numRow; i++)
			{
				for (j=0; j<numCol; j++)
					printf("%.2f\t", buf[(i*numCol)+j]);
				printf("\n");
			} 

		}

		printf("\n\n\n");
    }
    else
		MPI_Send(mat, numRow*numCol, MPI_FLOAT, 0, 0, MPI_COMM_WORLD);

    return 0;
}



int main (int argc, char *argv[])
{
  FILE *fp;
  int count = 0, ind;
  float *A, *B, *C, *bufRecv, *bufSend, somma;
  int numRow, numCol, startInd;
  int nproc, mype; 
  int i, j, y, k, sup, r;

  MPI_Init(&argc,&argv);
  MPI_Comm_size(MPI_COMM_WORLD, &nproc);
  MPI_Comm_rank(MPI_COMM_WORLD, &mype);

  /* Get local matrix dimensions */
  numRow = DIM / nproc;
  numCol = DIM;
  /* Avoid remainders */
  if (r = DIM % nproc)
    {
      if (mype == 0) printf("The number of process must divide %d exactly.\n",DIM);
      MPI_Finalize();
      return -1;      
    }

  /* Allocate workspace */
  A = (float *)malloc(numCol*numRow*sizeof(float));
  B = (float *)malloc(numCol*numRow*sizeof(float));
  C = (float *)malloc(numCol*numRow*sizeof(float));

  bufRecv = (float *)malloc(numCol*numRow*sizeof(float));
  bufSend = (float *)malloc(numRow*numRow*sizeof(float));

  /* Initialize matrices */
  for(i=0;i<numRow;++i)
  {
    for(j=0;j<numCol;++j)
      {
			A[(i*numCol)+j] = (((numRow * mype) + (i+1)) * (j+1));
			B[(i*numCol)+j] = 1 / A[(i*numCol)+j];
			C[(i*numCol)+j] = 0;
      }
  }
  /* Perform multiplication */
  for(count=0;count<nproc;++count)
  {
      startInd = count * numRow;
      for(i=0;i<numRow;++i)
		for(j=0;j<numRow;j++)
		{
			bufSend[j+(i*numRow)] = B[(i*numCol)+startInd+j];
		}
/* - use MPI_Allgather to distribute columns of B[][] matrix: */
			---
		for(k=0;k<numRow;++k)
		{
			for(i=0, somma=0.0;i<numRow;++i)
			{
				somma=0.0;
				for(j=0;j<numCol;++j)
				{
					somma += A[(k*numCol)+j] * bufRecv[(numRow*j)+i];
				}
				C[(count*numRow)+(k*numCol)+i] = somma;
			}
		}

   }

  /* Print matrices */
  printMat(A, numRow, numCol, mype, nproc,"Matrix A:");
  printMat(B, numRow, numCol, mype, nproc,"Matrix B:");
  printMat(C, numRow, numCol, mype, nproc,"Matrix C:"); 

  /* Free workspace */
  free(A);
  free(B);
  free(C);
  free(bufSend);
  free(bufRecv);

  MPI_Finalize();
  return 0;

}