/* 
Jonathan Hauenstein - University of Notre Dame
Center for Applied Mathematics Minicourse on Parallel Computing
April 24, May 1, May 8, 2008
Sample programs are designed for educational use 
*/

#include "trapezoid_dyn.h"
#include <mpi.h>

void create_trapInfo(MPI_Datatype *mpi_trapInfo)
{ // create the MPI datatype mpi_trapInfo
  MPI_Aint sAdd, add; trapInfo A;

  // arrays for length, displacement and datatypes
  int blockLen[3] = {1,1,1}; // 1 element in each block
  MPI_Datatype types[3] = {MPI_INT, MPI_DOUBLE, MPI_DOUBLE};
  MPI_Aint indices[3] = {0,0,0};

  // calculate the starting address
  MPI_Address(&A.n, &sAdd);
  // calculate the second address and compute the displacement
  MPI_Address(&A.a, &add); indices[1] = add - sAdd;
  // calculate the third address and compute the displacement
  MPI_Address(&A.b, &add); indices[2] = add - sAdd;

  // build the mpi datatype
  MPI_Type_struct(3, blockLen, indices, types, mpi_trapInfo);
  // commit the mpi datatype
  MPI_Type_commit(mpi_trapInfo);

  return;
}

double head_trapezoid(int totalTrap, double a, double b, int num_proc)
{ // calculate the trapezoid approximation using a dynamic distribution
  // n == 0 means that it is time to quit
  int i, active, count = 0, maxSize = 100000;
  double tempD, curr_a = a, stepSize = (b - a) / totalTrap, ans = 0;
  trapInfo info;
  MPI_Datatype mpi_trapInfo;
  MPI_Status status;

  // create mpi_trapInfo
  create_trapInfo(&mpi_trapInfo);

  // send out the first set of packets and count the number that are active
  active = 0;
  for (i = 1; i < num_proc; i++)
  { // determine the size of the packet
    if (count + maxSize < totalTrap)
    { // send a full packet
      count += info.n = maxSize;
    }
    else
    { // send out a partial packet
      count += info.n = totalTrap - count;
    }
    // setup a & b
    info.a = curr_a;
    curr_a = info.b = curr_a + stepSize * info.n;
    // send it off
    MPI_Send(&info, 1, mpi_trapInfo, i, 0, MPI_COMM_WORLD);
    // update active
    if (info.n > 0) active++;
  }

  // wait to recv answers back
  while (count < totalTrap)
  { // wait to recv
    MPI_Recv(&tempD, 1, MPI_DOUBLE, MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
    // update ans & active
    ans += tempD; active--;

    // send next packet
    if (count + maxSize < totalTrap)
    { // send a full packet
      count += info.n = maxSize;
    }
    else
    { // send out a partial packet
      count += info.n = totalTrap - count;
    }
    // setup a & b
    info.a = curr_a;
    curr_a = info.b = curr_a + stepSize * info.n;
    // send it off to where the last packet came from
    MPI_Send(&info, 1, mpi_trapInfo, status.MPI_SOURCE, 0, MPI_COMM_WORLD);
    // update active
    if (info.n > 0)
      active++;
  }

  // wait until all active has sent data back
  while (active > 0)
  { // wait to recv
    MPI_Recv(&tempD, 1, MPI_DOUBLE, MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
    // update ans
    ans += tempD;

    // send empty packet back
    info.n = 0;
    // send it off to where the last packet came from
    MPI_Send(&info, 1, mpi_trapInfo, status.MPI_SOURCE, 0, MPI_COMM_WORLD);

    // decrement active
    active--;
  }

  // clear mpi_trapInfo
  MPI_Type_free(&mpi_trapInfo);

  return ans;
}

void worker_trapezoid(int my_id, int num_proc)
{ // recv data, compute trapezoid approximation and send it back
  // n == 0 means that it is time to quit
  double ans;
  trapInfo info;
  MPI_Datatype mpi_trapInfo;
  MPI_Status status;

  // setup mpi_trapInof
  create_trapInfo(&mpi_trapInfo);

  do
  { // recv a packet
    MPI_Recv(&info, 1, mpi_trapInfo, MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &status);

    // see if we need to compute an approximation
    if (info.n > 0)
    { // perform trapezoid calculation
      ans = composite_trapApprox(&info);

      // return data back
      MPI_Send(&ans, 1, MPI_DOUBLE, status.MPI_SOURCE, 0, MPI_COMM_WORLD);
    }

  } while (info.n > 0);

  // clear mpi_trapInfo
  MPI_Type_free(&mpi_trapInfo);

  return;
}