逍遙文工作室

已經設定好OpenMP程式設計環境後，接下來就是實際來寫一則code，仔細分析thread們如何運作，在此我選擇積分這個主題，以梯形法則 (Trapezoidal Rule)來計算曲線下的面積。

簡單來講，即是計算梯形面積。為了讓曲線下面積的估量更加精確，a到b曲線下的面積可以細分成n塊，接著對這n塊梯形面積作估量。方便說明，我選y=x*x這條曲線，a=0 b=8 n=8，設定4條thread去計算8塊梯形面積，最後加總面積為172平方單位。

/**
	Theme: OpenMP Trapezoidal Rule
	Compiler: Visual Studio 2010 Professional
	Date: 100/10/25
	Author: HappyMan
	Blog: https://cg2010studio.wordpress.com/
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>

double f(double x);    /* Function we're integrating */
double Trap(double a, double b, int n);

double  global_result = 0.0;  /* Store result in global_result */
double  a, b;                 /* Left and right endpoints      */
int     n;                    /* Total number of trapezoids    */
int     thread_count;
int *thread_number;
double *index_area;
int i;

int main() {
   a=0,b=8,n=8,thread_count=4;
   global_result = Trap(a, b, n);

   printf("n = %d\n", n);
   printf("%.1f -> %.1f = %.3e\n", a, b, global_result);
   printf("#index sequence: \n");
   printf("index, thread_number, index_area\n");
   for(i=0; i<n; i++)
     printf("%5d, %13d, %10.3f\n",i,thread_number[i],index_area[i]);

   system("pause");
   return 0;
}  /* main */

/*------------------------------------------------------------------
 * Function:    f
 * Purpose:     Compute value of function to be integrated
 * Input arg:   x
 * Return val:  f(x)
 */
double f(double x) {
   double return_val;

   return_val = x*x;
   return return_val;
}  /* f */

/*------------------------------------------------------------------
 * Function:    Trap
 * Purpose:     Use trapezoidal rule to estimate definite integral
 * Input args:
 *    a: left endpoint
 *    b: right endpoint
 *    n: number of trapezoids
 * Return val:
 *    approx:  estimate of integral from a to b of f(x)
 */
double Trap(double a, double b, int n) {

   double h, global_area, local_area;
   int i;
   int my_rank;
   int x;

   thread_number=(int *)malloc(n * sizeof(int));
   index_area=(double *)malloc(n * sizeof(double));

   h = (b-a)/n;
   global_area = 0;

   printf("#execution sequence: \n");
   printf("index, thread_number, index_area\n");
#  pragma omp parallel for num_threads(thread_count) \
      reduction(+: global_area)
   for (i = 0; i < n; i++){
     my_rank = omp_get_thread_num();
	 local_area = (f(a + i*h) + f(a + (i+1)*h))*h/2;
     global_area += local_area;
	 *(thread_number+i) = my_rank;
	 *(index_area+i) = local_area;;
	 printf("%5d, %13d, %10.3f\n",i,thread_number[i],index_area[i]);
   }
   return global_area;
}  /* Trap */

程式跑出來的結果：

#execution sequence:
index, thread_number, index_area
    2,             1,      6.500
    4,             2,     20.500
    0,             0,      0.500
    6,             3,     42.500
    3,             1,     12.500
    5,             2,     30.500
    1,             0,      2.500
    7,             3,     56.500
n = 8
0.0 -> 8.0 = 1.720e+002
#index sequence:
index, thread_number, index_area
    0,             0,      0.500
    1,             0,      2.500
    2,             1,      6.500
    3,             1,     12.500
    4,             2,     20.500
    5,             2,     30.500
    6,             3,     42.500
    7,             3,     56.500
請按任意鍵繼續 . . .

概念可由這張圖來解釋：

參考：WiKiTrapezoidal rule。