simpson.c

// This program seems to work.
#include "stdio.h"
#include "math.h"
#include "assert.h"

#if defined(_OPENMP)
	#include "omp.h"
#endif

#define EPS 0.001
#define NUM_THREADS 1

typedef double f_t(double x);
typedef struct {
	double so,se;
	double xo,xe; 
} calc_t;

double sim(f_t * f, double a, double b, int n) {
	double par = f(a)+f(b);
	calc_t s[NUM_THREADS];  // s means "state"
	double h=(b-a)/n;
	int twot = 2*NUM_THREADS; 
	double in = twot*h;
	assert (n % (twot) == 0);
	int k = n / twot, i=0, j=0, t;


#if defined(_OPENMP)
	omp_set_num_threads(NUM_THREADS);
	#pragma omp parallel private(t,i,j)
        {
		t=omp_get_thread_num();
#else
		t=0;
		//j=0
		//printf("Thread %d\n",t);
#endif	
		s[t].xe=a+h*(t+1);
		s[t].xo=s[t].xe+h;
		s[t].se=s[t].so=0.0;		
		
		for (i=0;i<k;i++) {
			s[t].se+=f(s[t].xe);
			s[t].so+=f(s[t].xo);
			s[t].xe=s[t].xe+in;
			s[t].xo=s[t].xe+h;
		//	j++;
		}
#if defined(_OPENMP)
		// printf("Thread end %d\n",t);
		// printf("Thread end %d %d\n",t,j);
	}
#endif
	double so=0.0,se=0.0;
  	for (j=0;j<NUM_THREADS;j++) {
		so+=s[j].so;
		se+=s[j].se;
	}
	return (b-a)/(3*n)*(par+4*so+2*se);
}

double sin_(double x) {
	return sin(x);
}

void test_sim_spec() {
	double r = sim(sin_, 0, 2*M_PI, 8*800);
	printf("test_sim_spec..");
	printf("int = %f\n",r);
	assert (fabs(r)<EPS);
	printf("DONE\n");
}
void test_sim_non_null() {
	double r = sim(sin_, 0, M_PI/2, 8*800);
	printf("test_sim_non_null..");
	printf(" r=%f \n", r);
	assert (fabs(r-1)<EPS);
	printf("DONE\n");
}

#define ITERS 1

void timing() {
	int i;
	double start = omp_get_wtime( );  
	for (i=0;i<ITERS;i++) {
		sim(sin_, 0, M_PI*2, 64000000);
	}
    	double end = omp_get_wtime( );  
	printf("start = %.16g\nend = %.16g\ndiff = %.16g\n",  
             start, end, end - start); 
}


int main() {
	test_sim_spec();
	test_sim_non_null();
	timing();
}