// This is an alternative parallel version of quadrature_sequential.c 
//
// It differs from quadratur_parallel.c by the way work is distributed.
// This one uses a 'round robin' strategy.
// 
// To compile: mpicc quadrature_parallel.c -lm
//
// This code also exhibits linear speedup. See listed output at the 
// bottom of this file.
//
// Ole Nielsen, SUT 2003

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


//Define the circle arc to be integrated.
long double fct(long double x) {
  return sqrt(1.0-x*x);
}


int main(argc,argv)
     int argc;
     char *argv[];
{
  int N, i, lo, hi;
  long double result=0.0, my_result=0.0;
  long double reference = 3.1415926535897932; //Reference value, 17 digits

  
  double t0;
  int p;  // Myid
  int P;  // Number of processors 

  int namelen, tag=0, j, k;
  char processor_name[MPI_MAX_PROCESSOR_NAME];  
  int *C;

  MPI_Status status;
    
  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &p);
  MPI_Comm_size(MPI_COMM_WORLD, &P);
  MPI_Get_processor_name(processor_name, &namelen);  
   
  printf("This is process %d of %d running on %s\n", 
	 p, P, processor_name);

  N = 100000000;   // Resolution  
  //N = 29;   // Resolution  
  
  if (p == 0) {
    printf("Resolution = %d\n", N);
  }

  // Begin computation
  t0 = MPI_Wtime();

  
  for (k=0; k<15; k++) {
    //Each process p handles own stride of points
    my_result = 0.0;
    for (i=p; i<N; i+=P) {
      my_result += fct( (i+0.5)/N );
    } 

    my_result = 4*my_result/N;
  }
	 
  printf("P%d: Partial result = %.16Lf\n", p, my_result);   

  
  // Communication phase
  if (p == 0) {
  
    //Receive results from others 
    result = my_result;  
    for (j=1; j<P; j++) {          
      printf("P%d: receiving from P%d\n", p, j);
      MPI_Recv(&my_result, 1, MPI_LONG_DOUBLE, j, tag,
	       MPI_COMM_WORLD, &status);
      result += my_result;
    }
     
    printf("Error  = %.2e\n", fabs(result-reference)); 
    printf("Time   = %.2f sec\n", MPI_Wtime() - t0);
    
  } else {
    // send my_result to intended proc 0
    MPI_Send(&my_result, 1, MPI_LONG_DOUBLE, 0, tag,
	     MPI_COMM_WORLD);     
  }  

  MPI_Finalize();
  return 0;
}

//Observed output, Resolution = 100000000
//P = 1
//Error = 3.44e-13
//Time = 13.52 sec
//
//P = 2
//Error = 3.44e-13
//Time = 6.76 sec
//
//P = 4
//Error = 3.44e-13
//Time = 3.39 sec
//
//P = 8
//Error = 3.44e-13
//Time = 1.69 sec
//
//
//SPEEDUP & EFFICIENCY
//t = array([13.52, 6.76, 3.39, 1.69])
//
//t[0]/t
//
//array([ 1.        ,  2.        ,  3.98820059,  8.        ])
//
//p = [1, 2, 4, 8]
//
//t[0]/t/p
//array([ 1.        ,  1.        ,  0.99705015,  1.        ])