/* 
% 
% Written by:
% -- 
% John L. Weatherwax                2005-05-26
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

Iterative trapezoidal rule (parallel).

*/ 

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <pthread.h>

/* Global storage for the integral: */
double sum = 0.0; 

/* Protect this integral with a mutex lock: */
pthread_mutex_t sum_lock = PTHREAD_MUTEX_INITIALIZER;

/* Storage for anything I need to pass to the worker threads */ 
typedef struct {
  int thread_number; 
  double xL, xR; 
} package_t; 

/* The function to integrate: */ 
double f( double x ){
  return 4/(1+x*x);
  /* return sin(x)*exp(x); */
  /* return 1.0; */
}

/* Worker routine for each thread: 

   This rouine will integrate the function f(x) 
   over ONE step.
*/ 
void *mult_worker( void *arg ){
  
  int thread_number,status;
  double xL,xR,tmpSum;

  package_t *p = (package_t*) arg; 

  /* If this thread finished release its resources to the "main" routine */
  status = pthread_detach( pthread_self() ); 
  if( status != 0 ) {
    fprintf(stderr,"pthread_detach failed\n"); 
    exit(2);
  }

  /* Extract this threads "package" */ 
  thread_number = p->thread_number; 
  xL = p->xL;
  xR = p->xR;

  printf( "I AM thread # %4d ... let me work!\n", thread_number );
   
  /* Do all work I can without locking anything: */
  /* Can I do function evalutations atomically? */ 
  tmpSum = 0.5*(xR-xL)*( f(xR) + f(xL) );
  /* Lock before updating sum: */
  pthread_mutex_lock(&sum_lock); 
  sum += tmpSum; 
  pthread_mutex_unlock(&sum_lock); 

  free(p); 

  printf( "I AM thread # %4d ... just finished!\n", thread_number );
}

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

  int n,ti;
  double xL,xR,tXR,tXL,dx;
  pthread_t *threads;
  package_t *p;                  /* Argument list each thread will get */ 
  
  /* Argument checking */
  if( argc != 4 ){
    fprintf( stderr, "Expecting three input arguments\n" ); 
    exit(1); 
  }
  /* Input the left/right end points of the inteval */ 
  xL = atof(argv[1]); 
  xR = atof(argv[2]);
  /* Input the number of intervals: */ 
  n = atoi(argv[3]);

  /* Create a number of pthread_t objects: */ 
  threads = (pthread_t *) malloc( n*sizeof(pthread_t) ); 
  if( !threads ) 
    perror( "malloc" ); 

  /* Divide work; creating threads to do each subinterval */ 
  dx = (xR-xL)/n;
  for( ti=0; ti < n; ti++ ){
    tXL = xL + dx*ti; 
    tXR = tXL + dx; 
    printf( "Thread %3d will integrate over [%6f,%6f]\n", ti, tXL, tXR ); 
    
    /* Create and fill the args needed for the worker threads: */ 
    p = (package_t*) malloc( sizeof( package_t ) ); 
    p->thread_number = ti; 
    p->xL = tXL; 
    p->xR = tXR; 

    if( pthread_create( &threads[ti], NULL, 
			mult_worker, (void *) p ) != 0 ){
      perror( "pthread_create" ); 
      exit(2); 
    }
  }

  /* Synchronize the completion of all threads */
  for( ti=0; ti < n; ti++ ){
    pthread_join(threads[ti],NULL); 
    printf( "Thread %d is joining the fray\n", ti ); 
  }
  /* Clear the thread memory ... no longer needed */ 
  free(threads);

  /* Display the resulting integral: */ 
  printf("Integral = %f\n", sum); 
  


}