/*
 *  OpenMP lecture exercises
 *  Copyright (C) 2011 by Christian Terboven <terboven@rz.rwth-aachen.de>
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
 *
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <omp.h>
#include <unistd.h>


double CalcPi(int n);

int main(int argc, char **argv)
{
    int n = 150000000;
    const double fPi25DT = 3.141592653589793238462643;
    double fPi;
    double fTimeStart, fTimeEnd;
    
#ifdef READ_INPUT  
    printf("Enter the number of intervals: ");
    scanf("%d",&n);
#endif

    if (n <= 0 || n > 2147483647 ) 
    {
        printf("\ngiven value has to be between 0 and 2147483647\n");
        return 1;
    }
    fTimeStart = omp_get_wtime();

	
    /* the calculation is done here*/
    fPi = CalcPi(n);
	

    fTimeEnd = omp_get_wtime();
    printf("\npi is approximately = %.20f \nError               = %.20f\n",
           fPi, fabs(fPi - fPi25DT));
    printf("  wall clock time     = %.20f\n", fTimeEnd - fTimeStart);

#   pragma omp parallel 
   {
#   pragma omp for
	for (int i=0; i < omp_get_num_threads(); i++) 
	printf("Hello from thread %d \n", omp_get_thread_num());
#pragma omp critical
	{
	sleep(1);
 	}    
	printf("sleep %d\n", omp_get_thread_num());
   }	
    return 0;
}


double f(double a)
{
#pragma pomp inst begin(foo)
    double z=(4.0 / (1.0 + a*a));
#pragma pomp inst end(foo)
    return z;
}


double CalcPi(int n)
{
    const double fH   = 1.0 / (double) n;
    double fSum = 0.0;
    double fX;
    int i;
#pragma omp parallel for private(i) reduction(+:fSum)
    for (i = 0; i < n; i += 1)
    {
        fX = fH * ((double)i + 0.5);
        fSum += f(fX);
    }
    return fH * fSum;
}