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

/* Global storage for the integral: */
double sum = 0.0; 
pthread_mutex_t sum_mutex = PTHREAD_MUTEX_INITIALIZER; 
/* Global counter of thread number: */ 
int thread_id = 0; 

/* 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 ) {
    /* free( inPck ); */ 
    pthread_mutex_lock(&sum_mutex);
    sum += aL+aR;
    pthread_mutex_unlock(&sum_mutex); 
    return; 
  }
  else{
    /* Create two threads to do this subintegral: */ 
    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;
  }
}

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

  double tol,xL,fL,xR,fR,area0;
  package_t *istruct;
  
  /* Argument checking */
  if( argc != 4 ){
    fprintf( stderr, "Expecting three 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]);

  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); 

}