/*-------------------------------------------------------------------
 * This program implements the concurrent wave equation described  
 * in Chapter 5 of Fox et al., 1988, Solving Problems on Concurrent
 * Processors, vol 1.   
 *
 * A vibrating string is decomposed into points.  Each processor is  
 * responsible for updating the amplitude of a number of points over
 * time.
 *       
 * At each iteration, each processor exchanges boundary points with
 * nearest neighbors.  This version uses low level sends and receives
 * to exchange boundary points.
 *
 * Revised: 11/27/97 by Saleh Elmohamed
 *----------------------------------------------------------------- */

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

#define MASTER 0 
#define MAXPOINTS 1000        /* recommended max number of points */
#define MAXSTEPS 10000        /* recommended max number of steps  */
#define PI 3.14159265

int E_RtoL = 10;
int E_LtoR = 20;
int E_OUT1 = 30;
int E_OUT2 = 40;

void init_master(void);
void init_workers(void);
void init_line(void);
void update (int left, int right);
void output_master(void);
void output_workers(void);
MPI_Status status;
MPI_Request request;
int procid,              /* processor ID */
    nproc,               /* number of processes */
    nsteps,              /* number of time steps */
    tpoints,             /* total points along string */
    npoints,             /* number of points handled by this processor */
    first,               /* index of 1st point handled by this processor */
    rcode;               /* generic return code */
double values[MAXPOINTS+2],  /* values at time t */
       oldval[MAXPOINTS+2],  /* values at time (t-dt) */
       newval[MAXPOINTS+2];  /* values at time (t+dt) */

double begin_time, stop_time;

/* time_t  t1, t2; */
/*------------------------------------------------------
 *      Master obtains input values from user
 *------------------------------------------------------*/

void init_master(void)
{
  int buffer[2];
  
  /*==============================================*/
  /* Set number of points, number of iterations   */
  /* These numbers can be set by the user. It is  */
  /* recommended that tpoints not to exceed 1000  */
  /* and nsteps not to exceed 10000.              */
  /*==============================================*/

  tpoints = 100;
  nsteps = 1000;

#if 0
  while ((tpoints < nproc) || (tpoints > MAXPOINTS))
    {
      printf("Enter number of points along vibrating string\n");
      scanf("%s", tchar);
      tpoints = atoi(tchar);
      if ((tpoints < nproc) || (tpoints > MAXPOINTS))
	printf("enter value between %d and %d\n", nproc, MAXPOINTS);
    }
  while ((nsteps < 1) || (nsteps > MAXSTEPS))
    {
      printf("Enter number of time steps\n");
      scanf("%s", tchar);
      nsteps = atoi(tchar);
      if ((nsteps < 1) || (nsteps > MAXSTEPS))
	printf("enter value between 1 and %d\n", MAXSTEPS);
    }
#endif

  printf("Processor %d: points = %d, steps = %d\n", 
	 procid, tpoints, nsteps);

  /* then broadcast time advance parameter, total points, time steps */
  buffer[0] = tpoints;
  buffer[1] = nsteps;
  MPI_Bcast(&buffer, 2, MPI_INT, MASTER, MPI_COMM_WORLD);
}

/*-------------------------------------------------------
 *      Workers receive input values from master
 *------------------------------------------------------*/
  
void init_workers(void)
{
  int buffer[2];

  /* then receive total points, time steps */
  MPI_Bcast(&buffer, 2, MPI_INT, MASTER, MPI_COMM_WORLD);
  tpoints = buffer[0];
  nsteps = buffer[1];
}

/*---------------------------------------------------------
 *     All processes initialize points on line
 *------------------------------------------------------- */

void init_line(void)
{
  int nmin, nleft, npts, i, j, k;
  double x, fac;

  /* calculate initial values based on sine curve */
  nmin = tpoints/nproc;
  nleft = tpoints%nproc;
  fac = 2.0 * PI;
  for (i = 0, k = 0; i < nproc; i++)
    {
      npts = (i < nleft) ? nmin + 1 : nmin;
      if (procid == i)
	{
	  first = k + 1;
	  npoints = npts;
	  printf ("Processor %d: first = %d, npoints = %d\n", 
		  procid, first, npts);
	  for (j = 1; j <= npts; j++, k++)
	    {
	      x = (double)k/(double)(tpoints - 1);
	      values[j] = sin (fac * x);
	    }  
	}
      else 
	k += npts;
    }
  for (i = 1; i <= npoints; i++) 
    oldval[i] = values[i];
}

/*------------------------------------------------------------
 *      Calculate new values using wave equation
 *-----------------------------------------------------------*/
  
void do_math(int i)
{
  const double dtime = 0.3;
  const double c = 1.0;
  const double dx = 1.0; 
  double tau, sqtau;

  tau = (c * dtime / dx);
  sqtau = tau * tau;
  newval[i] = (2.0 * values[i]) - oldval[i]  
    + (sqtau * (values[i-1] - (2.0 * values[i]) + values[i+1]));
}

/*--------------------------------------------------------------------
 *      All processes update their points a specified number of times  
 *-------------------------------------------------------------------*/
  
void update(int left, int right)
{
  int i, j, id_rtol, id_ltor;

  /* update values for each point along string */
  for (i = 1; i <= nsteps; i++)
    {
      /* exchange data with "left-hand" neighbor */
      if (first != 1)  
	{
	  MPI_Isend(&values[1], 1, MPI_DOUBLE, left, E_RtoL,
		    MPI_COMM_WORLD, &request);
	  MPI_Wait(&request, &status);
	  MPI_Recv(&values[0], 1, MPI_DOUBLE, left, E_LtoR,
                   MPI_COMM_WORLD, &status);
	}
      /* exchange data with "right-hand" neighbor */
      if (first + npoints -1 != tpoints)
	{
	  MPI_Isend(&values[npoints], 1, MPI_DOUBLE, right,
		    E_LtoR, MPI_COMM_WORLD, &request);
	  MPI_Wait(&request, &status);
	  MPI_Recv(&values[npoints+1], 1, MPI_DOUBLE, right,
		   E_RtoL, MPI_COMM_WORLD, &status);
	}
      /* update points along line */
      for (j = 1; j <= npoints; j++)
	{
	  /* global endpoints */
	  if ((first + j - 1 == 1) || (first + j - 1 == tpoints))
	    newval[j] = 0.0;
	  else
	    do_math(j);
	}
      for (j = 1; j <= npoints; j++)
	{
	  oldval[j] = values[j];
	  values[j] = newval[j];
	}
    }
}

/*-----------------------------------------------------------
 *      Master receives results from workers and prints
 *---------------------------------------------------------- */

void output_master(void)  
{
  int i, start, npts, tpts, buffer[2];
  double results[MAXPOINTS];
  
  /* store worker's results in results array */
  for (i = 1; i < nproc; i++)
    {
      /* receive number of points and first point */
      MPI_Recv(&buffer, 2, MPI_INT, MPI_ANY_SOURCE, E_OUT1,
	       MPI_COMM_WORLD, &status);
      
      start = buffer[0];
      npts = buffer[1];
	 
      /* receive results */
      MPI_Recv(&results[start-1], npts, MPI_DOUBLE, MPI_ANY_SOURCE,
	       E_OUT2, MPI_COMM_WORLD, &status);
    }

  /* store master's results in results array */
  for (i = first; i < first + npoints; i++)
    results[i-1] = values[i];

  tpts = (tpoints < 10) ? tpoints: 10;  
  printf("The first %d points (just for validation):\n", tpts);
  for (i = 0; i < tpts; i++) 
    printf("%4.2f  ", results[i]);
  printf("\n");
}

/*-----------------------------------------------------------
 *      Workers send the updated values to the master
 *----------------------------------------------------------*/
  
void output_workers(void)
{
  int buffer[2];

  /* send first point and number of points handled to master */
  buffer[0] = first;
  buffer[1] = npoints;
  MPI_Isend(&buffer, 2, MPI_INT, 0, E_OUT1, MPI_COMM_WORLD, &request);
  MPI_Wait(&request, &status);

  /* send results to master */
  MPI_Isend(&values[1], npoints, MPI_DOUBLE, 0, E_OUT2, MPI_COMM_WORLD,
	    &request);
  MPI_Wait(&request, &status);
}

/*  -----------------------------------------------
 *      Main program
 *  --------------------------------------------- */

main(int argc, char **argv)
{
  int left, right;

  /* learn number of tasks and task ID */
  rcode = MPI_Init(&argc, &argv);
  begin_time = MPI_Wtime();
  MPI_Comm_rank(MPI_COMM_WORLD, &procid);
  MPI_Comm_size(MPI_COMM_WORLD, &nproc);

  if (rcode >= 0) 
    printf ("Processor %d: Wave Program running ... \n", 
	    procid);

  /* determine left and right neighbors */
  if (procid == nproc-1) 
    right = 0;
  else 
    right = procid + 1;
  if (procid == 0) 
    left = nproc - 1;
  else 
    left = procid - 1;

  /* get program parameters and initialize wave values */
  if (procid == MASTER) 
    init_master();
  else 
    init_workers();
  init_line();

  /*  t1 = time(NULL); */
  /* update values along the line for nstep time steps */
  update(left, right);


  /* collect results and print */
  if (procid == MASTER) 
    output_master();
  else 
    output_workers();

  stop_time = (MPI_Wtime() - begin_time);
  printf("... it took %6.3f seconds\n", stop_time);

  /*  t2 = time(NULL); */

  /*
  printf("The whole thing took about %4.2f seconds\n", 
         difftime(t2, t1));
  */

  MPI_Finalize();
}