用openmp实现梯形积分法

/**
* Inplement trapezoidal integration with OPENMP
* @filename: 
* 	trapezoidal_integration.c
* @compile: 
* 	gcc trapezoidal_integration.c -o trapezoidal_integration.out -fopenmp
* @note：
* 	n must be divided by thread_num with no remainder
* 	n means the number of tiny trapezoidal
*/

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

#ifdef _OPENMP
	#pragma message "Complier did support OPENMP"
	#include <omp.h>
#else
	#pragma message "Complier did NOT support OPENMP"
#endif

double f(double x){
	double y = 2*x; 
	return y;
} 

double area(
		double 	a				/*in*/,
		double 	b				/*in*/,
		int 	n				/*in*/,
		double 	*global_area	/*out*/){
	
	#ifdef _OPENMP
		int my_rank = omp_get_thread_num();
		int thread_count = omp_get_num_threads();
	#else
		int my_rank = 0;
		int thread_count = 1;
	#endif
	
	//The height of each tiny trapezoid
	double h = (b-a)/n; 
	//Number of trapezoids processed by each thread
	int local_n = n/thread_count;
	
	double local_start = a+my_rank*local_n*h;
	double local_per_area = 0;
	double x1 = local_start,x2=local_start+h;
	double local_total_area = 0;
	for(int i=0;i<local_n;++i){
		local_per_area = (f(x1)+f(x2))*h/2.0;
		x1 = x2;
		x2 = x1+h;
		local_total_area += local_per_area;
	}
	#pragma omp critical
	{
	printf("thread %d is going to enter the critical section\n",my_rank); 
	*global_area+=local_total_area;
	printf("thread %d has gone out the critical section\n",my_rank);
	} 
	
	return 0;
}

int main(int argc,char* argv[]){
	int thread_count = strtol(argv[1],NULL,10);
	double a,b;
	int n;
	
	scanf("%lf%lf%d",&a,&b,&n);
	//To ensure b>a
	if(a>b){
		int temp = a;
		a = b;
		b = temp;
	}
	double global_area = 0.0;

	#ifdef _OPENMP
	#pragma omp parallel num_threads(thread_count)
	#endif
	area(a,b,n,&global_area);

	printf("global_area=%lf\n",global_area);
	return 0;
}