/* 
% 
% 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.
% 
%-----

Recursive trapezoidal rule (parallel).

*/ 

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

typedef struct {
  double xL, fL, xR, fR, area0, tol; 
} package_t;

int max_number_of_threads; 
/* Global storage for the thead count and a protective mutex: */ 
int thread_count = 0; 
pthread_mutex_t thread_count_mutex = PTHREAD_MUTEX_INITIALIZER;

/* Global storage for the integral and a protective mutex: */
double sum = 0.0; 
pthread_mutex_t sum_mutex = PTHREAD_MUTEX_INITIALIZER; 

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

void * 
rec_quad( void *argIn ){

  /* Local variables: */ 
  double tol,xL,xR,fL,fR,area0;
  double xM,fM,aL,aR;
  pthread_t *threadL, *threadR;
  package_t *lPackage, *rPackage; 

  /* When run with a VERY strict tolerance this code runs out of memory.
     To alleviate this I placed the following code in this subroutine, 
     but this seemed to screw things up further ??? 
  : If this thread finishs, release its resources to the "main" routine :
  if( pthread_detach( pthread_self() ) != 0 ) {
    fprintf(stderr,"pthread_detach failed\n"); 
    exit(2);
  }
  */ 

  package_t *inPck = (package_t *) argIn; 
  /* Extract into local variables: */ 
  tol = inPck->tol; 
  xL = inPck->xL; xR = inPck->xR; 
  fL = inPck->fL; fR = inPck->fR; 
  area0 = inPck->area0; 

  /* 
     1) Split the interval into left and right subareas
     2) compute the approximate area of each subareas using the trapezoidal rule
     3) Stop if the relative tolerance is met
  */
  xM = 0.5*(xL+xR); fM = f(xM); 
  aL = 0.5*(xM-xL)*(fM+fL); aR = 0.5*(xR-xM)*(fM+fR); 
  printf( "rec_quad #(): aL=%f, aR=%f, aL+aR=%f, area0=%f\n", aL, aR, aL+aR, area0 );
  printf( "            : fabs(area0-(aL+aR))=%f; fabs(area0-(aL+aR))/area0=%f\n", 
	  fabs(area0-(aL+aR)), fabs(area0-(aL+aR))/area0 ); 

  if( fabs(area0-(aL+aR))/area0 <= tol ) {

    printf( "I've converged, adding in my bit...\n" ); 
    pthread_mutex_lock(&sum_mutex);
    sum += aL+aR;
    pthread_mutex_unlock(&sum_mutex); 
    return; 

  }else{
    
    /* Create two threads IF I can still create them.
       (limited by max_number_of_threads): */
    pthread_mutex_lock(&thread_count_mutex);  
    if( thread_count+2 <= max_number_of_threads ){ /* Do work with threads... */ 
      thread_count += 2; 
      pthread_mutex_unlock(&thread_count_mutex); 
      /* Create two threads to do this subintegral: */ 
      printf( "I've not converged yet, lets thread out this work...\n" ); 
      lPackage = (package_t *) malloc( sizeof(package_t) ); 
      if( !lPackage ){
	fprintf(stderr,"No memory\n"); 
	exit(2); 
      }else{
	lPackage->tol = tol; 
	lPackage->xL = xL; lPackage->fL = fL; 
	lPackage->xR = xM; lPackage->fR = fM; 
	lPackage->area0 = aL; 
      }
      rPackage = (package_t *) malloc( sizeof(package_t) ); 
      if( !rPackage ){
	fprintf(stderr,"No memory\n"); 
	exit(2); 
      }else{
	rPackage->tol = tol; 
	rPackage->xL = xM; rPackage->fL = fM; 
	rPackage->xR = xR; rPackage->fR = fR; 
	rPackage->area0 = aR; 
      }
      threadL = (pthread_t *) malloc( sizeof(pthread_t) ); 
      /* Spawn two new threads to do each subinterval: */ 
      if( pthread_create( threadL, NULL, rec_quad, (void*) lPackage ) != 0 ){
	perror( "pthread_create" ); 
	exit(2); 
      }
      threadR = (pthread_t *) malloc( sizeof(pthread_t) ); 
      if( pthread_create( threadR, NULL, rec_quad, (void*) rPackage ) != 0 ){
	perror( "pthread_create" ); 
	exit(2);
      }
      /* Wait for the threads to join and clear used memory */ 
      pthread_join(*threadL,NULL);
      free( lPackage ); free( threadL ); 
      pthread_join(*threadR,NULL); 
      free( rPackage); free( threadR ); 
      return;

    }else{                                      /* Do work serially... */ 
      pthread_mutex_unlock(&thread_count_mutex); 
      printf( "I've not converged yet, lets serialize the remaining work...\n" ); 
      /* 
	 Rather than a recursive solution at this point I have choosen
	 an iterative one, I felt programming the recursive solution
	 would have been easier.  For instance it would be:

	 threadSum = rec_quad_nonthread( tol, xL, fL, xR, fR, area0 ); 

	 with the definition of rec_quad_nonthread is given in the function 
	 prob_1_7_a_ser_rec.c file.

	 For the algebra for this section see the handwritten notes
	 that go with this problem.
      */ 
      double threadSum, dx, xL_i, xR_i; 
      int niters, sc; 

      threadSum = 0.0; 
      niters = (int) ( ( (xR-xL)/sqrt(tol) ) + 1 );
      dx = (xR-xL)/niters;
      for( sc=0; sc < niters; sc++ ){
	xL_i = xL + dx*sc;
	xR_i = xL_i + dx; 
	threadSum += 0.5*dx*(f(xL_i)+f(xR_i)); 
      }
      /* Finished, add to global sum... */ 
      pthread_mutex_lock(&sum_mutex); 
      sum += threadSum; 
      pthread_mutex_unlock(&sum_mutex);
      return; 

    } /* end thread_count < max_number_of_threads "else" clause */ 

  } /* end if( fabs(...) <= tol */ 
} /* end rec_quad */ 

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

  double tol,xL,fL,xR,fR,area0;
  package_t *istruct;
  
  /* Argument checking */
  if( argc != 5 ){
    fprintf( stderr, "Expecting four input arguments\n" ); 
    exit(1); 
  }
  /* Input the left/right end points of the interval */ 
  xL = atof(argv[1]); 
  xR = atof(argv[2]);
  /* Input the tolerance required: */ 
  tol = atof(argv[3]);
  /* Input the max number of threads to create: */ 
  max_number_of_threads = atoi(argv[4]);

  fL = f(xL); fR = f(xR); 
  area0 = 0.5*(xR-xL)*(fL+fR); 
  printf( "main: fL=%f, fR=%f, area0=%f\n", fL, fR, area0 );

  /* Package this information into a package type: */ 
  istruct = (package_t *) malloc( sizeof(package_t) ); 
  if( !istruct ){
    fprintf(stderr,"No memory\n"); 
    exit(2); 
  }
  istruct->tol = tol; 
  istruct->xL = xL; istruct->fL = fL; 
  istruct->xR = xR; istruct->fR = fR; 
  istruct->area0 = area0; 

  printf( "main: calling rec_quad...\n" ); 
  rec_quad( istruct ); 

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

}