并行程序设计导论学习笔记——OpenMP（1）
使用OpenMP需要在编译器上打开OpenMP开关，并包含omp.h文件。我使用的是在Windows下的Visual Studio 2015，只需在工程选项中打开OpenMP支持就可以了。按照书上的说法，GCC增加参数-fopenmp就可以了。
OpenMP有两个重要的函数：
omp_get_thread_num()
omp_get_num_threads()
他们的返回值都是无符号整数，第一个用来返回当前执行的线程编号，而第二个返回一共有的线程数量。和C语言的数组下标类似，当前执行的线程编号的值从0开始，最大值总是线程数量减一。至于这两个函数的用途，是进行线程的任务分配。
#pragma omp为OpenMP的预处理器指令，他将告诉编译器的OpenMP的具体行为。
在5.1节中主要涉及到了OpenMP的#pragma omp parallel指令，用来将当前指令并行化。注意这个指令并不包括任何的任务分配，只是简单的将他下面的指令送给不同线程去执行而已。所以，简单的一条#pragma omp parallel指令下面接一个for循环，而for循环的起止点都是常数而不是根据当前线程分配的，只会把这个for循环用线程数量执行一遍，与我们预期的“将内容拆开，分配给不同线程”的想法相悖，无法起到并行化的效果。
书后面提到#pragma omp parallel for就可以自动分配线程任务，但是具体的分配方式还是取决于具体的控制语句。
对于OpenMP指令，还可以跟上一些从句(clause)。比如，在上面的#pragma omp parallel，可以跟上num_threads(thread_count)，表示用thread_count的线程数量来运行。除此之外还有一条#pragma omp critical指令，用来在不同线程访问或者更新所有线程共享的资源（比如全局变量），用来防止竞争的情况(race condition)。
书上还说了一些概念性的东西：不同的线程是从进程forked出来的，共享大部分进程的资源（比如stdin，stdout，等），但是不同的线程享有不同的栈(stack)空间和程序计数器(program counter)。当一个线程执行完毕，将返回进程（的主线程）。不同的栈空间意味着局部变量是独立的，尽管可能是使用了同一个函数，但是他们的局部变量在OpenMP创建了子线程后依然是相互独立不相干的。
例题：用OpenMP实现定义法计算曲边梯形（定积分几何意义）的并行算法。
定义法计算定积分包括分割，求和，取极限三个步骤。串行算法比较容易实现，用循环累加分割的小区间的值就可以了。并行的时候，因为各个小区间求和互不相干（没有依赖性(dependence)和相互通信(communication)），因此只需要将不同小区间分配给不同线程就可以了。最后#pragma omp critical指令强调了汇总操作时只能由一个线程更新其共享资源（也相当于锁）。
例子的完整的实现如下：
#include<stdio.h>
#include<time.h> /* clock() function */
#include<stdlib.h>
#include<omp.h> /* OpenMP header */
#include<math.h>
double result_parallel;
#define N 400000000
double get_approx_serial(double(*func)(double x), double interval_a, double interval_b, long long int n)
{
double h;
double result = 0;
double x_i;
int i;
h = (interval_b - interval_a) / n;
for (i = 0; i < n; i++)
{
x_i = interval_a + h*i;
result += func(x_i)*h;
}
return result;
}
void get_approx_parallel(double(*func)(double x), double interval_a, double interval_b, long long int n)
{
double h;
double result = 0;
double x_i;
long long int local_i_start, local_i_end;
int thread_id, total_threads;
int i;
/* Get thread num */
total_threads = omp_get_num_threads();
thread_id = omp_get_thread_num();
/* Task Arrangment between threads */
local_i_start = (n / total_threads)*thread_id;
if (thread_id == total_threads - 1)
local_i_end = n;
else
local_i_end = (n / total_threads)*(thread_id + 1);
/* approx Algorithm */
h = (interval_b - interval_a) / n;
for (i = local_i_start; i < local_i_end; i++)
{
x_i = interval_a + h*i;
result += func(x_i)*h;
}
#pragma omp critical
result_parallel += result;
return;
}
int main()
{
clock_t start_clock, end_clock;
int total_thread_num;
/* Get total threads to display */
#pragma omp parallel
total_thread_num = omp_get_num_threads();
printf("This computer Has %d Threads in total\n", total_thread_num);
/* Do serial algorithm and get the used time. */
start_clock = clock();
printf("The area between y=sin(x) , y=0 , x=0 ,x=pi is %lf\n", get_approx_serial(&sin, 0, 3.141593, N));
end_clock = clock();
printf("Serial total use Time %d ms", end_clock - start_clock);
/* Do parallel algorithm and get the used time. */
start_clock = clock();
#pragma omp parallel
get_approx_parallel(&sin, 0, 3.141593, N);
printf("\nThe area between y=sin(x) , y=0 , x=0 ,x=pi is %lf\n", result_parallel);
end_clock = clock();
printf("Parallel total use Time %d ms", end_clock - start_clock);
system("pause");
return 0;
}
critical语句的反汇编情况如下，证明了编译器在解释这条语句时转换成了两个函数调用，来处理线程的竞争冲突。
程序的运行结果如下所示：我的机器是6核12线程的，所以一共有12条线程可以使用，并行算法比串行算法实际加速比(Speedup)$S=T_s/T_p=9912/984=10.07$，效率(efficiency)$E=(1/p)S=0.839$（计算方法和定义见书58页）。根据阿姆达尔定律(Amdahl's law) ，随着数据规模$N$的增大，串行代码占比几乎为零（只有#pragma omp critical子句声明的为串行代码），这个算法可以达到理论上几乎线性的加速比。实际上因为SMT超线程效率有限，加速比为并未达到12。
参考书籍：并行程序设计导论（英文版），机械工业出版社。
