/* File:    mpi_simpson.c
 * Purpose: Calculate area using Simpson's rule and MPI
 *
 * Input:   a, b, n
 * Output:  Estimate of area between x-axis, x = a, x = b, and graph of f(x)
 *          using n subintervals Simpson's rule.  Run-time of parallel
 *          Simpson's rule code.
 *
 * Compile: mpicc -g -Wall -o mpi_simpson mpi_simpson.c -lm
 * Usage:   mpiexec -n <number of processes> ./mpi_simpson
 *
 * Notes:    
 * 1.  The function f(x) is hardwired.
 * 2.  The number of subintervals, n, and the number of processes, p,
 *     should satisfy n/p is even.
 */

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

void Get_data(int p, int my_rank, double* a_p, double* b_p, int* n_p);
double Simpson(double  a, double  b, int n, double h);
double f(double x);    /* Function we're integrating */
double Global_sum(double local_area, int my_rank, int p);
double Get_max_time(double par_elapsed, int my_rank, int p);

int main(void) {
   double  area;            /* Store result in area       */
   double  a, b;            /* Left and right endpoints   */
   int     n;               /* Number of subintervals     */
   double  h;               /* Length of subintervals     */
   int     p, my_rank;
   int     local_n;
   double  local_a, local_b;
   double  local_area;
   double  start, finish, my_elapsed, par_elapsed;
   
   MPI_Init(NULL, NULL);
   MPI_Comm_size(MPI_COMM_WORLD, &p);
   MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);

   Get_data(p, my_rank, &a, &b, &n);

   start = MPI_Wtime();
   h = (b-a)/n;
   local_n = n/p;
   local_a = a + my_rank*local_n*h;
   local_b = local_a + local_n*h;
   local_area = Simpson(local_a, local_b, local_n, h);
   area = Global_sum(local_area, my_rank, p);
   finish = MPI_Wtime();
   my_elapsed = finish - start;

   par_elapsed = Get_max_time(my_elapsed, my_rank, p);

   if (my_rank == 0) {
      printf("With n = %d subintervals, our estimates\n", n);
      printf("of the area from %f to %f is\n", a, b);
      printf("    %.15f with Simpson's rule\n\n", area);
      printf("Elapsed time = %e\n", par_elapsed);
   }

   MPI_Finalize();
   return 0;
}  /* main */

/*------------------------------------------------------------------
 * Function:     Get_data
 * Purpose:      Read in the data on process 0 and send to other
 *               processes
 * Input args:   p, my_rank
 * Output args:  a_p, b_p, n_p
 */
void Get_data(int p, int my_rank, double* a_p, double* b_p, int* n_p) {
   int        q;
   MPI_Status status;

   if (my_rank == 0) {
      printf("Enter a, b, and n\n");
      scanf("%lf %lf %d", a_p, b_p, n_p);

      for (q = 1; q < p; q++) {
         MPI_Send(a_p, 1, MPI_DOUBLE, q, 0, MPI_COMM_WORLD);
         MPI_Send(b_p, 1, MPI_DOUBLE, q, 0, MPI_COMM_WORLD);
         MPI_Send(n_p, 1, MPI_INT, q, 0, MPI_COMM_WORLD);
      }
   } else {
      MPI_Recv(a_p, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
      MPI_Recv(b_p, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
      MPI_Recv(n_p, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &status);
   }
}  /* Get_data */


/*------------------------------------------------------------------
 * Function:     Simpson
 * Purpose:      Estimate area using Simpson's rule
 * Input args:   a: left endpoint
 *               b: right endpoint
 *               n: number of subintervals
 *               h: stepsize = length of subintervals
 * Return val:   Simpson's rule estimate of area between x-axis,
 *               x = a, x = b, and graph of f(x)
 */
double Simpson(double  a, double  b, int n, double h) {
   double area, x, temp;
   int i;

   area = f(a) + f(b);

   temp = 0.0;
   for (i = 1; i <= n-1; i += 2) {
      x = a + i*h;
      temp += f(x);
   }
   area += 4*temp;
   
   temp = 0;
   for (i = 2; i <= n-2; i += 2) {
      x = a + i*h;
      temp += f(x);
   }
   area += 2*temp;

   area *= h;
   area /= 3;

   return area;
}  /* Simpson */


/*------------------------------------------------------------------
 * Function:    f
 * Purpose:     Compute value of function to be integrated
 * Input args:  x
 */
double f(double x) {
   double return_val;

   return_val = x*x + 1;
   return return_val;
}  /* f */


/*------------------------------------------------------------------
 * Function:   Global_sum
 * Purpose:    Find the sum of all the process' areas
 * In args:    my_rank:       calling process' rank
 *             p:             total number of processes
 *             local_area:    elapsed time on calling process
 * Ret val:    Process 0:     sum of all processes areas
 *             Other procs:   input value for area
 */
double Global_sum(double local_area, int my_rank, int p) {
   int source;
   MPI_Status status;
   double temp;
   double area = local_area;

   if (my_rank == 0) {
      for (source = 1; source < p; source++) {
         MPI_Recv(&temp, 1, MPI_DOUBLE, source, 0, MPI_COMM_WORLD, 
               &status);
         area += temp;
      }
   } else {
      MPI_Send(&local_area, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
   }
   return area;

}  /* Global_sum */


/*------------------------------------------------------------------
 * Function:   Get_max_time
 * Purpose:    Find the maximum elapsed time across the processes
 * In args:    my_rank:       calling process' rank
 *             p:             total number of processes
 *             par_elapsed:   elapsed time on calling process
 * Ret val:    Process 0:     max of all processes times
 *             Other procs:   input value for par_elapsed
 */
double Get_max_time(double par_elapsed, int my_rank, int p) {
   int source;
   MPI_Status status;
   double temp;

   if (my_rank == 0) {
      for (source = 1; source < p; source++) {
         MPI_Recv(&temp, 1, MPI_DOUBLE, source, 0, MPI_COMM_WORLD, 
               &status);
         if (temp > par_elapsed) par_elapsed = temp;
      }
   } else {
      MPI_Send(&par_elapsed, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
   }
   return par_elapsed;
}  /* Get_max_time */