/* 
   Compute integral from 0  to 1  4/(1+x*x). 
   This shoudl be equal to PI. 
   

   Second attempt at parallelization: idea: make i, x, pi, sum private. 
   This works. 

[ltoma@lobster /Volumes/CS/ltoma/public_html/teaching/cs3225-GIS/fall16/Code/OpenMP]$./a.out 
Start(2)
Start(0)
Start(4)
Start(6)
Start(1)
Start(5)
Start(7)
Start(3)
Done(0): pi = 3.142
Done(1): pi = 3.142
Done(4): pi = 3.142
Done(2): pi = 3.142
Done(5): pi = 3.142
Done(3): pi = 3.142
Done(7): pi = 3.142
Done(6): pi = 3.142
Done. time = 0.1


   But every thread does the whole work.  We need to split the loop
   among the threads.  OpenMP has special support to parallelize loops
   and makes it about as easy as it could be.

*/ 


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



int main(int argc, char** argv) {

  long num_steps = 10000000; 
  double step, x, pi, sum; 
  int i; 
  double t1, t2; 
  
  step = 1/(double)num_steps; 
  
  t1 = omp_get_wtime(); 
  

#pragma omp parallel  private (i, x, pi, sum)
  //another way would be to declare i,x,pi,sum locally inside the
  //block, then they are private
   { 
     sum = 0;  //<----- this is important to be inside
     int thread_id = omp_get_thread_num(); 
     printf("Start(%d)\n", thread_id); 
     
     for (i=0; i<num_steps; i++) {
       //value in the middle of the rectangle
       x = (i+.5)*step; 
       sum = sum + 4/(1+x*x); 
     }
     pi = step * sum; 
          
     printf("Done(%d): pi = %.3f\n", thread_id, pi); 
  }

  t2 = omp_get_wtime();   
  printf("Done. time = %.1f\n", t2-t1);
  
  



  return 0; 
}