/* trap.cpp -- Parallel Trapezoidal Rule, first version
 *
 * Input: None.
 * Output:  Estimate of the integral from a to b of f(x)
 *    using the trapezoidal rule and n trapezoids.
 *
 * Algorithm:
 *    1.  Each process calculates "its" interval of
 *        integration.
 *    2.  Each process estimates the integral of f(x)
 *        over its interval using the trapezoidal rule.
 *    3a. Each process != 0 sends its integral to 0.
 *    3b. Process 0 sums the calculations received from
 *        the individual processes and prints the result.
 *
 * Notes:  
 *    1.  f(x), a, b, and n are all hardwired.
 *    2.  The number of processes (p) should evenly divide
 *        the number of trapezoids (n = 1024)
 *
 * See Chap. 4, pp. 56 & ff. in PPMPI.
 */
#include <iostream.h>

/* We'll be using MPI routines, definitions, etc. */
#include <mpi++.h>


int main(int argc, char** argv) {
    int         my_rank;   /* My process rank           */
    int         p;         /* The number of processes   */
    float       a = 0.0;   /* Left endpoint             */
    float       b = 10.0;   /* Right endpoint            */
    int         n = 1024;  /* Number of trapezoids      */
    float       h;         /* Trapezoid base length     */
    float       local_a;   /* Left endpoint my process  */
    float       local_b;   /* Right endpoint my process */
    int         local_n;   /* Number of trapezoids for  */
                           /* my calculation            */
    float       integral;  /* Integral over my interval */
    float       total;     /* Total integral            */
    int         source;    /* Process sending integral  */
    int         dest = 0;  /* All messages go to 0      */
    int         tag = 0;
    MPI::Status  status;

    float Trap(float local_a, float local_b, int local_n,
              float h);    /* Calculate local integral  */

    /* Let the system do what it needs to start up MPI */
    MPI::Init(argc, argv);

    /* Get my process rank */
    my_rank = MPI::COMM_WORLD.Get_rank();

    /* Find out how many processes are being used */
    p = MPI::COMM_WORLD.Get_size();

    h = (b-a)/n;    /* h is the same for all processes */
    local_n = n/p;  /* So is the number of trapezoids */

    /* Length of each process' interval of
     * integration = local_n*h.  So my interval
     * starts at: */
    local_a = a + my_rank*local_n*h;
    local_b = local_a + local_n*h;
    integral = Trap(local_a, local_b, local_n, h);

    /* Add up the integrals calculated by each process */
    if (my_rank == 0) {
        total = integral;
        for (source = 1; source < p; source++) {
            MPI::COMM_WORLD.Recv(&integral, 1, MPI::FLOAT, source, tag, status);
            total = total + integral;
        }
    } else {  
        MPI::COMM_WORLD.Send(&integral, 1, MPI::FLOAT, dest, tag);
    }

    /* Print the result */
    if (my_rank == 0) {
        cout << "With n = " << n << " trapezoids, our estimate" << endl;
        cout << "of the integral from " << a << " to " << b 
             << " = " << total << endl;
    }

    /* Shut down MPI */
    MPI::Finalize();

    return 0;
} /*  main  */


float Trap(
          float  local_a   /* in */,
          float  local_b   /* in */,
          int    local_n   /* in */,
          float  h         /* in */) {

    float integral;   /* Store result in integral  */
    float x;
    int i;

    float f(float x); /* function we're integrating */

    integral = (f(local_a) + f(local_b))/2.0;
    x = local_a;
    for (i = 1; i <= local_n-1; i++) {
        x = x + h;
        integral = integral + f(x);
    }
    integral = integral*h;
    return integral;
} /*  Trap  */


float f(float x) {
    float return_val;
    /* Calculate f(x). */
    /* Store calculation in return_val. */
    return_val = x*x;
    return return_val;
} /* f */

/*
mpiCC -o trapCPP trap.cpp
export LAMRSH="ssh"
lamboot -v lamhosts

LAM 6.5.6/MPI 2 C++/ROMIO - University of Notre Dame

Executing hboot on n0 (helium.tjhsst.edu - 1 CPU)...
Executing hboot on n1 (sulfur.tjhsst.edu - 1 CPU)...
Executing hboot on n2 (phosphorus.tjhsst.edu - 1 CPU)...
Executing hboot on n3 (magnesium.tjhsst.edu - 1 CPU)...
topology done

mpirun N trapCPP
With n = 1024 trapezoids, our estimate
of the integral from 0 to 10 = 333.334

lamhalt

LAM 6.5.6/MPI 2 C++/ROMIO - University of Notre Dame

*/