/**********************************************************************
 *                      III TTTTTTT U     U                           *
 *                       I     T    U     U                           *
 *                       I     T    U     U                           *
 *                       I     T    U     U                           *
 *                      III    T    UUUUUUU                           *
 *                                                                    *
 *  ITU Informatics Institute                                         *
 *  Computational Science & Engineering Dept.                         *
 *  Istanbul / TURKIYE                                                *
 *                                                                    *
 *  By:H.Hakan GUREL                                                  *
 *  E-mail: gurel@be.itu.edu.tr                                       *
 *  Date: 13.10.2005                                                  *
 *                                                                    *
 *  Program: This is matrix matrix multiplication program. It         *
 *   calculates the matrix multiplication with defining matrix        *
 *   dimension in the code. A and B matrixes store the data which     *
 *   is generated with randomly in time. It can be seen the created   *
 *   matrixes and the result matrix using -DVERBOSE=3 flag while      *
 *   compiling. Program also considers the load balancing concept and *
 *   divide the data among processors equally. Master node only       * 
 *   creates the matrixes,also doing sending, receiving data and      *
 *   displaying the result so data is sharing among the total number  *
 *   of processor -1.                                                 *
 *                                                                    *
 *********************************************************************/

#include<stdio.h>

#include<math.h>
#include"mpi.h"
#define NRA 100            /* number of rows in matrix A */
#define NCA 100            /* number of columns in matrix A and rows in B*/
#define NCB 100            /* number of columns in matrix B */
#define MASTER 0           /* taskid of first task */
#define FROM_MASTER 1      /* setting a message type */
#define FROM_WORKER 2      /* setting a message type */

#ifndef VERBOSE
#define VERBOSE 0
#endif


MPI_Status status;
main(int argc, char **argv) 
{
int numtasks,           /* number of tasks in partition */
    taskid,             /* a task identifier */
    numworkers,         /* number of worker tasks */
    source,             /* task id of message source */
    dest,               /* task id of message destination */
    nbytes,             /* number of bytes in message */
    mtype,              /* message type */
    intsize,            /* size of an integer in bytes */
    dbsize,             /* size of a double float in bytes */
    rows,               /* rows of matrix A sent to each worker */
    averow, extra, offset,      /* used to determine rows sent to each worker */
    i, j, k,            /* misc */
    count; 
double a[NRA][NCA],     /* matrix A to be multiplied */
       b[NCA][NCB],     /* matrix B to be multiplied */
       c[NRA][NCB];     /* result matrix C */

double start, etime, 
       ttime, startcom, ecomtime, maxtime;
intsize = sizeof(int);
dbsize = sizeof(double);

MPI_Init(&argc, &argv);
start=MPI_Wtime();
MPI_Comm_rank(MPI_COMM_WORLD, &taskid);
MPI_Comm_size(MPI_COMM_WORLD, &numtasks);
numworkers = numtasks-1;
srand(time(NULL));
/**************************** master task ************************************/
if (taskid == MASTER) {
  printf("Number of worker tasks = %d\n",numworkers);

#if VERBOSE>2
  printf("\nMATRIX A\n"); 
#endif
  for (i=0; i<NRA; i++)
  {
    for (j=0; j<NCA; j++)
    {
      a[i][j]= rand()%100;
#if VERBOSE>2
      printf(" %6.2f ", a[i][j]);
#endif
    }
#if VERBOSE>2
   printf("\n");
#endif
  }

#if VERBOSE>2
  printf("\nMATRIX B\n"); 
#endif
  for (i=0; i<NCA; i++)
  {
    for (j=0; j<NCB; j++)
    {
      b[i][j]= rand()%100;
#if VERBOSE>2
      printf(" %6.2f ", b[i][j]);
#endif
    }
#if VERBOSE>2
    printf("\n");
#endif
  }
startcom=MPI_Wtime();
  /* send matrix data to the worker tasks */
  averow = NRA/numworkers;
  extra = NRA%numworkers;
  offset = 0;
  mtype = FROM_MASTER;
  for (dest=1; dest<=numworkers; dest++) 
	  {               
    rows = (dest <= extra) ? averow+1 : averow;     
    printf("******   sending %d rows to task %d\n",rows,dest);
    MPI_Send(&offset, 1, MPI_INT, dest, mtype, MPI_COMM_WORLD);
    MPI_Send(&rows, 1, MPI_INT, dest, mtype, MPI_COMM_WORLD);
    count = rows*NCA;
    MPI_Send(&a[offset][0], count, MPI_DOUBLE, dest, mtype, MPI_COMM_WORLD);
    count = NCA*NCB;
    MPI_Send(&b, count, MPI_DOUBLE, dest, mtype, MPI_COMM_WORLD);

    offset = offset + rows;
	    }

  /* wait for results from all worker tasks */
  mtype = FROM_WORKER;
  for (i=1; i<=numworkers; i++)      
	  {                 
    source = i;
    MPI_Recv(&offset, 1, MPI_INT, source, mtype, MPI_COMM_WORLD, &status);
    MPI_Recv(&rows, 1, MPI_INT, source, mtype, MPI_COMM_WORLD, &status);
    count = rows*NCB;
    MPI_Recv(&c[offset][0], count, MPI_DOUBLE, source, mtype, MPI_COMM_WORLD, 
               &status);
 
	    }
	ecomtime=(MPI_Wtime()-startcom);
//   print results 
#if VERBOSE>2
  printf("Here is the result matrix\n");
#endif
  for (i=0; i<NRA; i++) 
    { 
#if VERBOSE>2
    printf("\n"); 
#endif
    for (j=0; j<NCB; j++)
{ 
#if VERBOSE>2
      printf(" %6.2f ", c[i][j]);
#endif
}    

#if VERBOSE>2
  printf ("\n");
#endif
}
  }  /* end of master section */

/**************************** worker task ************************************/

if (taskid > MASTER) {
  mtype = FROM_MASTER;
  source = MASTER;
  MPI_Recv(&offset, 1, MPI_INT, source, mtype, MPI_COMM_WORLD, &status);
  MPI_Recv(&rows, 1, MPI_INT, source, mtype, MPI_COMM_WORLD, &status);
    count = rows*NCA;
  MPI_Recv(&a, count, MPI_DOUBLE, source, mtype, MPI_COMM_WORLD, &status);
    count = NCA*NCB;
  MPI_Recv(&b, count, MPI_DOUBLE, source, mtype, MPI_COMM_WORLD, &status);

  for (k=0; k<NCB; k++)
    for (i=0; i<rows; i++) {
      c[i][k] = 0.0;
      for (j=0; j<NCA; j++)
        c[i][k] = c[i][k] + a[i][j] * b[j][k];
      }

  mtype = FROM_WORKER;

  MPI_Send(&offset, 1, MPI_INT, MASTER, mtype, MPI_COMM_WORLD);
  MPI_Send(&rows, 1, MPI_INT, MASTER, mtype, MPI_COMM_WORLD);
  MPI_Send(&c, rows*NCB, MPI_DOUBLE, MASTER, mtype, MPI_COMM_WORLD);

  }  /* end of worker */
etime=(MPI_Wtime()-start);
/* compute total time*/

	if(taskid == MASTER)
{
	printf(" Calculation time %6.3f seconds\n", etime);

}

  MPI_Finalize();
} /* of main */