/*---------------------------------------------------------------------
 * 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 task is  
 * responsible for updating the amplitude of a number of points over
 * time.
 *       
 * At each iteration, each task exchanges boundary points with
 * nearest neighbors.  This version uses low level sends and receives
 * to exchange boundary points.
 *
 * This version illustrates use of persistent communication.
 *
 * REVISED 4/15/96 R. Leibensperger   revised timing
 * REVISED 8/30/96 R. Leibensperger   use of MPI_PROC_NULL
 * Revised 11/15/97 Saleh Elmohamed to run on SunHPC MPI
 *---------------------------------------------------------------------*/
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <sys/time.h>
#include <time.h>
#include "mpi.h"
#define MASTER 0 
#define MAXPROC 8
#define MAXPOINTS 1000
#define MAXSTEPS 1000000
#define PI 3.14159265

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

void get_data(void);
void init_line(void);
void pcomm_create (int right, int left);
void update (int left, int right);
void output_master(void);
void output_workers(void);
MPI_Request request;
MPI_Status status;
MPI_Request Prequest[4];
MPI_Status Pstatus[4];
int rank,                  /* rank */
  ntask,                   /* number of tasks */
  nsteps,                  /* number of time steps */
  tpoints,                 /* total points along string */
  npoints,                 /* number of points handled by this task */
  first,                   /* index of 1st point handled by this task */
  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) */

/*  ------------------------------------------------------------------------
 *      Obtain input values from user
 *  ------------------------------------------------------------------------ */

void get_data(void)
{
  char tchar[8];
  int buffer[2];

  if (rank == MASTER) {
    /* get number of points, number of iterations from user */
    tpoints = 1000;
    nsteps = 10000;
#if 0
    while ((tpoints < ntask) || (tpoints > MAXPOINTS))
      {
	printf("Enter number of points along vibrating string\n");
	scanf("%s", tchar);
	tpoints = atoi(tchar);
	if ((tpoints < ntask) || (tpoints > MAXPOINTS))
	  printf("enter value between %d and %d\n", ntask, 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("%d: points = %d, steps = %d\n", rank, tpoints, nsteps);
    buffer[0] = tpoints;
    buffer[1] = nsteps;
  }

  /* then broadcast number of points and number of iterations */
  MPI_Bcast(buffer, 2, MPI_INT, MASTER, MPI_COMM_WORLD);
  tpoints = buffer[0];
  nsteps = buffer[1];
}

/*  ------------------------------------------------------------------------
 *     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/ntask;
  nleft = tpoints%ntask;
  fac = 2.0 * PI;
  for (i = 0, k = 0; i < ntask; i++)
    {
      npts = (i < nleft) ? nmin + 1 : nmin;
      if (rank == i)
	{
	  first = k + 1;
	  npoints = npts;
	  printf ("%d: first = %d, npoints = %d\n", rank, 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];
}

/*  -------------------------------------------------------------------------
 *      Create persistent communication requests
 *  -------------------------------------------------------------------------*/

void pcomm_create (int right, int left)
{
  int i;
  /*  Initialize all requests to NULL
   *  All tasks wait for all four requests, even though some tasks
   *  do not use all four requests */

  for (i=0;i<4;i++) Prequest[i] = MPI_REQUEST_NULL;

  /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
   * Step 1: Create persistent communication requests for:           *
   *         sending to left-hand neighbor        Prequest[0]        *
   *         receiving from left-hand neighbor    Prequest[1]        *
   *         sending to right-hand neighbor       Prequest[2]        *
   *         receiving from right-hand neighbor   Prequest[3]        *
   * Remember: the arguments to the initialization calls are the     *
   * same as for the non-blocking communications you are replacing   *
   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

  /* create requests for communication with left-hand neighbor */
  if (left != MPI_PROC_NULL)
    {
      MPI_Send_init(&values[1], 1, MPI_DOUBLE, left, E_RtoL,
		    MPI_COMM_WORLD, &Prequest[0]);
      MPI_Recv_init(&values[0], 1, MPI_DOUBLE, left, E_LtoR,
		    MPI_COMM_WORLD, &Prequest[1]);
    }
  /* create requests for communication with right-hand neighbor */
  if (right != MPI_PROC_NULL)
    {
      MPI_Send_init (&values[npoints], 1, MPI_DOUBLE, right,
		     E_LtoR, MPI_COMM_WORLD, &Prequest[2]);
      MPI_Recv_init (&values[npoints+1], 1, MPI_DOUBLE, right,
		     E_RtoL, MPI_COMM_WORLD, &Prequest[3]);
    }
}


/*---------------------------------------------------------------------
 *      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 tasks update their points a specified number of times  
 *---------------------------------------------------------------------*/
  
void update(int left, int right)
{
  int i, j;
  struct timeval tv1, tv2;
  int dt=0;

  /* update values for each point along string */
  for (i = 1; i <= nsteps; i++)
    {
      /*  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
       * Step 2:                                                      *
       * Replace the four non-blocking communication calls with       *
       * calls that "start" the four persistent communications.       *
       *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  */

      /* time the communication */
      gettimeofday(&tv1, (struct timeval*)0);

      /* exchange data with "left-hand" neighbor */
      if (left != MPI_PROC_NULL)  
	{
	  MPI_Start(&Prequest[0]);
	  MPI_Start(&Prequest[1]);
	}
      /* exchange data with "right-hand" neighbor */
      if (right != MPI_PROC_NULL)
	{
	  MPI_Start(&Prequest[2]);
	  MPI_Start(&Prequest[3]);
	}
      MPI_Waitall(4, Prequest, Pstatus);
      /*  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  */

      /* complete the timing */
      gettimeofday(&tv2, (struct timeval*)0);
      dt = dt + ((tv2.tv_sec - tv1.tv_sec) * 1000000 + tv2.tv_usec 
		 - tv1.tv_usec);
      
      /* 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];
	}
    }
  printf("time for communication in update loop   uSec = %8d\n", dt);
}

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

void output_master(void)  
{
  int i, start, npts, buffer[2], istep, source;
  double results[MAXPOINTS];
  
  /* store worker's results in results array */
  for (i = 1; i < ntask; 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 */
      source = status.MPI_SOURCE;
      MPI_Recv(&results[start-1], npts, MPI_DOUBLE, 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];

  istep = (tpoints <= 10) ? 1: tpoints/10;
  printf ("Points for validation:\n");
  for (i = 0; i < tpoints; i+=istep)
    printf ("%d:%4.2f  ", i, results[i]);
  if (i-istep != tpoints - 1) 
    printf ("%d:%4.2f  ", tpoints-1, results[tpoints-1]);
  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
 *------------------------------------------------------------------*/

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

  /* learn number of tasks and rank in MPI_COMM_WORLD */

  rcode = MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  MPI_Comm_size(MPI_COMM_WORLD, &ntask);
  if (rcode >= 0) printf ("%d: Wave Program running\n", rank);

  /* determine left and right neighbors */
  if (rank == ntask-1) right = MPI_PROC_NULL;
  else right = rank + 1;
  if (rank == 0) left = MPI_PROC_NULL;
  else left = rank - 1;

  /* get program parameters and initialize wave values */
  get_data();
  init_line();

  /* create persistent communication requests */
  pcomm_create(right, left);

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

  /*  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
   * Step 3:                                                   *
   *--- deallocate persistent communication requests here      */

  for (i=0; i<4; i++) 
    if (Prequest[i] != MPI_REQUEST_NULL) MPI_Request_free (&Prequest[i]);

  /*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*/
  /* collect results and print */

  if (rank == MASTER) output_master();
  else output_workers();

  MPI_Finalize();
  return(0);
}