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

int main(  int argc, char *argv[] )
{ 

   double *A, *B, *C, *tC;
   int L, M, N, L0, L1;
   double tempo2, tempo1, tempro;

   int i, j, k, t, st;
   double v, rmin, rmax;

/*  MPI variables */
   char server[MPI_MAX_PROCESSOR_NAME];
   int ierr, my_rank, numprocs, ls, dest, tag;
   MPI_Status status;
   int LetsGo;

   ierr = MPI_Init(&argc,&argv);
   ierr = MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
   ierr = MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
   ierr = MPI_Get_processor_name(server,&ls);

   LetsGo = 1;
   if ( my_rank == 0 ) {
      /* Master process looks for input data */
      fprintf(stdout,"Input Matrices dimensions L, M, N:\n");
      fscanf(stdin,"%d, %d, %d",&L, &M, &N);
      fprintf(stdout," Matrices dimensions are: %d, %d, %d\n",L,M,N);
      if ( LetsGo == 1 ) {
         /*  Master allocates entire matrices */
         A = (double*) malloc(sizeof((double)1.0)*L*M);
         B = (double*) malloc(sizeof((double)1.0)*N*M);
         C = (double*) malloc(sizeof((double)1.0)*L*N);
         if ( A == NULL || B == NULL || C == NULL ) {
            fprintf(stderr,"Error allocating matrices \n");
      	     LetsGo = 0;
         }
      }
   }
   /* Every process synchronizes */
   ierr = MPI_Barrier(MPI_COMM_WORLD);
   /* Send start flag */
   ierr = MPI_Bcast(&LetsGo,1,MPI_INT,0,MPI_COMM_WORLD);
   if ( LetsGo == 0 ) ierr = Terminate(my_rank);
   /* Send data to processes */
   ierr = MPI_Bcast(&L,1,MPI_INT,0,MPI_COMM_WORLD);
   ierr = MPI_Bcast(&M,1,MPI_INT,0,MPI_COMM_WORLD);
   ierr = MPI_Bcast(&N,1,MPI_INT,0,MPI_COMM_WORLD);
   if ( my_rank == 0 ) {
      fprintf(stdout," Initializing matrices on server %s. Input data are: %d, %d, %d\n",server, L, M, N);
   }

   /* Compute local rows and allocate local memory */
   ierr = estremi(L,&L0,&L1,numprocs,my_rank);
   if ( my_rank != 0 ) {
         A = (double*) malloc(sizeof((double)1.0)*(L1-L0)*M);
         B = (double*) malloc(sizeof((double)1.0)*N*M);
         C = (double*) malloc(sizeof((double)1.0)*(L1-L0)*N);
         if ( A == NULL || B == NULL || C == NULL ) {
            fprintf(stderr,"Error allocating matrices process %d\n",my_rank);
            ierr = Terminate(my_rank);
         }
   }

   /* Every process executes */
   for ( i = 0; i < (L1-L0)*M; i++ ) A[i] = 0.0;
   for ( i = 0; i < N*M; i++ ) B[i] = 0.0;
   for ( i = 0; i < (L1-L0)*N; i++ ) C[i] = 0.0;

   /*  Compute multiplication */
   for ( i = L0; i < L1; i++ ) {
      for ( k = 0; k < M; k++ ) {
         A[ (i-L0)*M + k] = (double)(i-k) / (double)(i+k+1);
      }
   }

   for ( j = 0; j < N; j++ ) {
      for ( k = 0; k < M; k++ ) {
         B[ k*N + j ] = 0.01 * (double)(k-j) / (double)(k+j+1);
      }
   }

   ierr = MPI_Barrier(MPI_COMM_WORLD);
   tempo1 = MPI_Wtime();
   /*  Compute multiplication */
   fprintf(stdout,"Process %d computes C(:,%d:%d) \n",my_rank,L0,L1-1);
   for ( j = 0; j < N; j++ ) {
      for ( i = L0; i < L1; i++ ) {
         C[ (i-L0)*N + j ] = 0.0;
         for ( k = 0; k < M; k++ ) {
            C[ (i-L0)*N + j ] = C[ (i-L0)*N + j ] + A[ (i-L0)*M + k ] * B[ k*N + j ];
         }
      }
   }
    
   if ( my_rank != 0 ) {
   /*   Send data to master */
      dest = 0; tag = 100;
      ierr = MPI_Send(&C[0], (N*(L1-L0)), MPI_DOUBLE, dest, 
                       tag, MPI_COMM_WORLD);
      fprintf(stderr," Process %d sends %d elements from %lf \n",my_rank,(N*(L1-L0)),C[0]);
   }
	
   if ( my_rank == 0 ) {
      
      /* Master process gathers results */
      for ( t = 1; t < numprocs; t++ ) {
         tag = 100;
         ierr = estremi(L,&L0,&L1,numprocs,t);
         ierr = MPI_Recv(&C[L0*N], (N*(L1-L0)), MPI_DOUBLE, t, tag, MPI_COMM_WORLD, &status);
         fprintf(stderr," Process %d receives  %d elements from %lf \n",my_rank,(N*(L1-L0)),C[L0*N]);
      }
   
      tempo2 = MPI_Wtime();
   
      tempro =  (tempo2 - tempo1);

      
      fprintf(stdout,"Elapsed time is %f = %f - %f \n",tempro, tempo2, tempo1);
      	   	   
      /* Compare results */ 
      tC = (double*) malloc(sizeof((double)1.0)*L*N);
      if ( tC == NULL ) {
         fprintf(stderr," Process %d: Error allocating matrix tC \n",my_rank);
      }
      /*  Define whole A matrix  */
      for ( i = 0; i < L; i++ ) {
         for ( k = 0; k < M; k++ ) {
            A[ i*M + k] = (double)(i-k) / (double)(i+k+1);
         }
      }

      for ( j = 0; j < N; j++ ) {
         for ( i = 0; i < L; i++ ) {
            tC[ i*N + j ] = 0.0;
            for ( k = 0; k < M; k++ ) {
               tC[ i*N + j ] = tC[ i*N + j ] + A[ i*M + k ] * B[ k*N + j ];
            }
         }
      }

      rmin = rmax = tC[0] - C[0];
      for ( j = 0; j < N; j++ ) {
         for ( i = 0; i < L; i++ ) {
            v = tC[ i*N + j ] - C[ i*N + j ];
            if ( v < rmin ) {
               rmin = v;
            }
            if ( v > rmax ) {
               rmax = v;
            }
         }
      }

      fprintf(stderr," Check results: rmin, rmax = %f, %f\n",rmin,rmax);
      
   }
   ierr = MPI_Barrier(MPI_COMM_WORLD);
   ierr = Terminate(my_rank);
   return(0);
}

int estremi(int length, int* start, int* end, int intervals, int n)
{
    int delta;

    delta = ( length + (intervals - 1) ) / intervals;
    *start = delta * n; *end = delta * ( n + 1);
    if ( *end > length ) *end = length;
   
   return(0);
}

int Terminate(int p)
{
  int ierr;
 
  fprintf(stderr," Process %d terminates\n",p);
  ierr = MPI_Finalize();
 
   return(0);
}