/*
 * Parallel Trapezoidal Rule
 *
 * 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 rules.
 * 	3a.	Each process != 0 sends its integral to 0.
 * 	3b.	process 0 sums the calculations received from the individual
 * 		processes and prints the result.
 *
 * Note: f(x), a, b, and n are hardwired.
 *
 * Note: This example from the book "Parallel Programming with MPI" by Peter
 * 	 S. Pacheco.
 */

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

#include <mpi.h>

float Trap(float, float, int, float);
float f(float);

int
main (int argc, char **argv)
{
	int my_rank;		/* process rank */
	int p;			/* number of processes */
	float a = 0.0;		/* left endpoint */
	float b = 1.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 calculations */
	float integral;		/* integral over my interval */
	float total;		/* total integral */
	int source;		/* process sending integral */
	int dest = 0;		/* all messages go to process 0 */
	int tag = 0;		/* tag for messages */
	MPI_Status status;	/* return status for receive */

	/* Initialize MPI */
	MPI_Init(&argc, &argv);

	/* Get my process rank */
	MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);

	/* How many processes are being used? */
	MPI_Comm_size(MPI_COMM_WORLD, &p);

	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's 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_Recv(&integral, 1, MPI_FLOAT, source, tag,
					MPI_COMM_WORLD, &status);
			total += integral;
		}
	} else {
		MPI_Send(&integral, 1, MPI_FLOAT, dest, tag, MPI_COMM_WORLD);
	}

	/* print the result */
	if (my_rank == 0) {
		printf("With n = %d trapezoids, our estimate\n", n);
		printf("of the integral from %f to %f = %f\n", a, b, total);
	}

	/* Shut down MPI */
	MPI_Finalize();

	return 0;
}

float
Trap(float local_a, float local_b, int local_n, float h)
{
	float integral;
	float x;
	int i;

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

float
f(float x)
{
	/* simplistic example - perform other calculations here */
	return sin(x);
}