Estimate of the integral f(x) using MPI parallel programming 2

/*
* Author : Emin Yüce
* Input:
* n: number of trapezoids.
* Output: Estimate of the integral from 1 to 500 of f(x) = 6x^6-7x^5+3x^2+11
* using the trapezoidal rule and n trapezoids.
*
*/
#include
#include "mpi.h"
/*Method declerations used in program*/
int Ceiling_log2(int x);
int I_receive( int stage, int my_rank, int* source_ptr);
int I_send(int stage, int my_rank, int p, int* dest_ptr);
void Send_Trapezoid(int* n_ptr, int my_rank,int p);
float Calculate_Integral(float local_a, float local_b, int local_n, float h); /* Calculate local integral */
void Collect_Local_Integral(float *,int,int);
void Send_Global_Integral(float* globalSum ,int my_rank /* in */,int p /* in */);
int Ceiling_log2(int x);
int I_receive_2(int stage ,int my_rank,int p, int * source_ptr);
int I_send_2(int stage, int my_rank, int* dest_ptr);
void Send2(float partial_sum , int dest) ;
void Receive2(float* partial_sum, int source);


main(int argc, char** argv) {

const float a = 1; /* Left endpoint */
const float b = 500; /* Right endpoint */
int my_rank; /* My process rank */
int p; /* The number of processes */
int n; /* 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;




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

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

/* Find out how many processes are being used */
MPI_Comm_size(MPI_COMM_WORLD, &p);

/*
*PE0 sends the input data to all other PEs
*PEs other than PE0 receive the input data that is sent
*/
Send_Trapezoid(&n, my_rank, 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' 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;

/*Partial calculation of integral*/
integral = Calculate_Integral(local_a, local_b, local_n, h);

/*Collect partial integral value on PE0 using tree strategy*/
/*Add up the integrals calculated by each process */
Collect_Local_Integral(&integral,my_rank, p);

// PE0 sends the global sum to other PEs
Send_Global_Integral(&integral,my_rank,p);


/*PE0 displays the global sums on the screen*/
if (my_rank == 0) {
printf("\nn = %d trapezoids\n", n);
printf("p = %d processor number\n",p);
printf("Process 0 sends the area of %d trapezoids to other processes at Log(p) times\n",n/p);
printf("Collect local computations at process 0 at Log(p) times \n");
printf("Inform global sum the all other processes with Log(p) times \n");
printf("The integral of 6x^6-7x^5+3x^2+11 function from %f to %f is %f\n", a, b, integral);

}
/* Shut down MPI */
MPI_Finalize();
} /* main */

/********************************************************************/
void Send_Trapezoid(int* n_ptr ,int my_rank /* in */,int p /* in */) {
/*
PE 0 sends number of trapezoidal to all other PEs
All PEs except PE0 receive the sending data
*/
int source;
int dest;
MPI_Status status;

if (my_rank == 0){
printf("The Program is written to calculate integral of f(x) = 6x^6-7x^5+3x^2+11 by using parallel computing \n");
printf("Enter number of trapezoids\n");
scanf("%d", n_ptr);
for (dest = 1; dest < p; dest++){
MPI_Send(n_ptr, 1, MPI_INT, dest, 0, MPI_COMM_WORLD);
}
}else{
MPI_Recv(n_ptr, 1, MPI_INT, MPI_ANY_SOURCE, 0, MPI_COMM_WORLD, &status);
}
}
/*
Determines Which PE is the destination for receive operation
Used for collecting partial sums
*/
int I_receive_2(int stage /* in */,int my_rank /* in */,int p /* in */,int * source_ptr /* out */) {

int power_2_stage;

/* 2^stage = 1 << stage */
power_2_stage = 1 << stage;
if (my_rank < power_2_stage){
*source_ptr = my_rank + power_2_stage;
if (*source_ptr >= p) return 0;
else return 1;
} else return 0;
} /* I_receive */

/*
Determines Which PE is the target for receive operation
Used for collecting partial sums
*/
int I_send_2(int stage /* in */,int my_rank /* in */, int* dest_ptr /* out */) {

int power_2_stage;

/* 2^stage = 1 << stage */
power_2_stage = 1 << stage;
if ((power_2_stage <= my_rank) && (my_rank < 2*power_2_stage)){
* dest_ptr = my_rank - power_2_stage;
return 1;
} else return 0;
}
/*
- all PE has a partial sum value
- collect sum with tree strategy in pe 0
*/
void Collect_Local_Integral(float * partial_sum,int my_rank /* in */,int p/* in */) {

int source;
int dest;
int stage;
float k=0; // temp


for (stage = Ceiling_log2(p) -1; stage >= 0 ; stage--){
if (I_receive_2(stage, my_rank,p, &source)){
k = *partial_sum;
MPI_Status status;
MPI_Recv(partial_sum, 1, MPI_FLOAT, source, 0,
MPI_COMM_WORLD, &status);
*partial_sum = k + *partial_sum;

/*PE0 displays the partial sums on the screen in the form.*/
if(my_rank == 0){
printf("Sender process: %d Partial sum: %f \n",source,*partial_sum);
}
}
else if (I_send_2(stage, my_rank, &dest)){
MPI_Send(partial_sum, 1, MPI_FLOAT, dest, 0, MPI_COMM_WORLD);
}
}
}
/********************************************************************/
int I_receive(int stage /* in */,int my_rank /* in */,int* source_ptr /* out */) {
int power_2_stage;

/* 2^stage = 1 << stage */
power_2_stage = 1 << stage;
if ((power_2_stage <= my_rank) &&
(my_rank < 2*power_2_stage)){
*source_ptr = my_rank - power_2_stage;
return 1;
} else return 0;
}
int I_send(int stage /* in */,int my_rank /* in */,int p /* in */,int* dest_ptr /* out */) {
int power_2_stage;

/* 2^stage = 1 << stage */
power_2_stage = 1 << stage;
if (my_rank < power_2_stage){
*dest_ptr = my_rank + power_2_stage;
if (*dest_ptr >= p) return 0;
else return 1;
} else return 0;
}
/*
PE 0 sends the global sum to all other PEs using tree strategy
*/
void Send_Global_Integral(float* globalSum ,int my_rank /* in */,int p /* in */) {

int source;
int dest;
int stage;
MPI_Status status;


for (stage = 0; stage < Ceiling_log2(p); stage++){
if (I_receive(stage, my_rank, &source)){
MPI_Recv(globalSum, 1, MPI_INT, source, 0,
MPI_COMM_WORLD, &status);
}
else if (I_send(stage, my_rank, p, &dest)){
MPI_Send(globalSum, 1, MPI_INT, dest, 0, MPI_COMM_WORLD);

}
}
}

/*
calculates the integral of function using trapezoidal numeric method
*/
float Calculate_Integral(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;
}
/*
return the value of the function for given x parameter
*/
float f(float x) {


float return_val;
float power (float base, int n); // the function that calculate power of given number

/* Calculate f(x). */
/* Store calculation in return_val. */
return_val = 6 * power(x,6) - 7 * power(x,5) + 3 * power(x,2) + 11;
return return_val;
} /* f */

/*
calculate the n th power of base number
*/
float power(float base, int n) {

int i;
float p;
p = 1;
for (i = 1; i <= n; ++i){
p *= base;
}
return p;
}

/********************************************************************/
int Ceiling_log2(int x /* in */) {
/* Ceiling of log_2(x) is just the number of times
* times x-1 can be divided by 2 until the quotient
* is 0. Dividing by 2 is the same as right shift.
*/

/* Use unsigned so that right shift will fill
* leftmost bit with 0
*/
unsigned temp = (unsigned) x - 1;
int result = 0;

while (temp != 0) {
temp = temp >> 1;
result = result + 1 ;
}
return result;
}