/*******************************************************************************
 * This is a simple program that shows the numerical integration example.
 * It computes pi by approximating the area under the curve:
 *      f(x) = 4 / (1+x*x) between 0 and 1.
 * To do this intergration numerically, the interval from 0 to 1 is divided into
 * some number (num_sub_intervals) subintervals and added up the area of
 * rectangles
 * The larger the value of the num_sub_interval the more accurate your result
 * will be.
 *
 * The program first asks the user to input a value for subintervals, it
 * computes the approximation for pi, and then compares it to a more 
 * accurate aproximate value of pi in the math.h library.
 *
 * Hwo to Compile:
 *  $ cc omp_pi_integration.c -o omp_pi_integration
 * 
 * How to Setup Runtime Environment:
 *  $ export OMP_NUM_THREADS=8
 * 
 * How to Run:
 *  $ aprun -n 1 -d 8 ./omp_pi_integration
 * 
 * Run the program multiple times. Each time Increase the number of threads.
 * Also when it ask you for the number of iterations, you can start with 1000
 * and each time you run add a zero, and see what have changed.
 * 
 *
 ******************************************************************************/
// The following C libraries are needed.
#include <stdio.h>	 // Is needed for printing the final results
#include <stdlib.h>  // Is needed for exiting early if an error occurs
#include <math.h>   // Is needed for fabs()/absolute value of floating point numbers
#include <omp.h>   // Is needed for MPI parallelization

int main(int argc, char *argv[]) {
    int num_sub_intervals = 0;   // number of sub intervals
    double start_time, end_time, time_diff;  //for timeing most computation intese part of the program
    double x, pi;
    double sum = 0.0; // for partial sum
    double step;
	int i;

    printf("Please enter the number of iterations used to compute pi: \n");
    scanf("%d",&num_sub_intervals);

    step = 1.0/(double) num_sub_intervals;

    // Record the start time from here on:
    start_time = omp_get_wtime();

// The next block of code in {} will be run in parallel by all the threads
// Since its a conditional parallelism, if the number of sub intervals are 
// more than 100, it will fork the number of threads setup by OMP_NUM_THREADS
// the dynamic scheduling will do an equal distribution of loop itterations or
// sub_intervals to each thread in the pool.  
#pragma omp parallel if(num_sub_intervals > 100) reduction(+:sum) private(x)
{   
	#pragma omp for schedule(dynamic) 
    for(i=0; i < num_sub_intervals; i++){
        x = (i+0.5)*step;
        sum += 4.0/(1.0+x*x);
    }
}
	// once all the threads compute their partial sums of each sub intervals
	// the thread 0, will finally put together all the partial sums.
    pi = step * sum;

    //End recording the time here.
    end_time = omp_get_wtime();
    time_diff = end_time - start_time;

    // print the result here: 
    printf("computed pi value is = %g (%17.15f)\n\n", pi,pi);
    printf("M_PI accurate value from math.h is: %17.15f \n\n", M_PI);
    printf("Difference between computed pi and math.h M_PI = %17.15f\n\n",
            fabs(pi - M_PI));
    printf("Time to compute = %g seconds\n\n", time_diff);

    return EXIT_SUCCESS;
}