OpenMP初探

 OpenMP 提供了对并行算法的高层的抽象描述，程序员通过在源代码中加入各种专用的pragma 指令（directive，指示/命令）来指明自己的意图，由此编译器可以自动将程序进行并行化，并在必要之处加入同步互斥以及通信。当选择忽略这些pragma，或者编译器不支持OpenMP 时，程序又可退化为通常的（串行）程序，代码仍然可以正常运作，只是不能利用多线程来加速程序执行。OpenMP 提供的这种对于并行描述的高层抽象降低了并行编程的难度和复杂度，这样程序员可以把更多的精力投入到并行算法本身，而非其具体实现细节。对基于数据分集的多线程程序设计，OpenMP 是一个很好的选择。同时，使用OpenMP 也提供了更强的灵活性，可以较容易的适应不同的并行系统配置。线程粒度和负载平衡等是传统多线程程序设计中的难题，但在OpenMP 中，OpenMP 库从程序员手中接管了部分这两方面的工作。但是，作为高层抽象，OpenMP 并不适合需要复杂的线程间同步和互斥的场合。OpenMP的另一个缺点是不能在非共享内存系统（如计算机集群）上使用。在这样的系统上，MPI使用较多。

OpenMP 编程技术

1）循环并行化
循环并行化（loop parallelize）是指使用OpenMP 的parallel for 指导语句将C/C++的for循环并行化，即将循环中的迭代（平均）分配给线程组中的各个线程分别执行后再汇总。

2）并行区域
并行区域（parallel region）指由OpenMP 的parallel 语句定义的并行控制语句块, 并行区域中的代码被所有的线程重复执行。

3）工作分区
工作分区（sections）是指利用OpenMP 的sections 编译指导语句，将用section 语句指定的内部代码，划分成多个工作区分配给线程组中的各个线程，不同的section 由不同的线程执行。各线程工作量的总合等于原来的工作量。这显然比并行区域所做的重复性劳动更有意义。

4）归约操作
归约操作是OpenMP 编程方式给同步编程带来的特殊的编程功能，该操作会反复将一个二元运算符应用在一个变量和另一个值上，并把结果保存在原变量中。该操作通过
reduction 语句提供。

下面我们用矩形法则的数值积分方法来估算π的值：

1）串行程序
static long num_steps = 100000;
double x, pi, sum = 0.0, step = 1.0/(double) num_steps;
for (long i = 1; i <= num_steps; i++){
x = (i + 0.5)*step;
sum += 4.0/(1.0 + x*x);
}

 

2）并行程序1(并行区域)
#define NUM_THREADS 2
73
static long num_steps = 100000;
double x, pi = 0.0, sum[NUM_THREADS], step = 1.0/(double)num_steps;
omp_set_num_threads(NUM_THREADS);
#pragma omp parallel
{
int id = omp_get_thread_num();
sum[id] = 0.0;
for (long i = id; i < num_steps; i += NUM_THREADS) {
x = (i + 0.5)*step;
sum[id] += 4.0/(1.0 + x*x); // 根据ID 号分配任务
}
}
for(int i = 0; i < NUM_THREADS; i++) pi += sum[i] * step;
printf("Pi = %f\n", pi);

 

3）并行程序2(循环并行化)
#define NUM_THREADS 2
static long num_steps = 100000;
double x, pi = 0.0, sum[NUM_THREADS], step = 1.0/(double)num_steps;
omp_set_num_threads(NUM_THREADS);
#pragma omp parallel
{
int id = omp_get_thread_num();
sum[id] = 0.0;
#pragma omp for
for (long i = 0; i < num_steps; i++)
{
x = (i + 0.5)*step;
sum[id] += 4.0/(1.0 + x*x); // 存储每次的结果
}
}
for(int i = 0; i < NUM_THREADS; i++) pi += sum[i] * step;
printf("Pi = %f\n", pi);

 

4）并行程序3(私有化和临界区)
#define NUM_THREADS 2
static long num_steps = 100000;
double x, pi = 0.0, sum, step = 1.0/(double)num_steps;
omp_set_num_threads(NUM_THREADS);
#pragma omp parallel private(x, sum)
{
int id = omp_get_thread_num();
sum = 0.0;
for (long i = id; i < num_steps; i += NUM_THREADS)
{
x = (i + 0.5)*step;
sum += 4.0/(1.0 + x*x);
}
#pragma omp critical
pi += sum*step;
}
printf("Pi = %f\n", pi);

 

5）并行程序4(归约操作)
#define NUM_THREADS 2
static long num_steps = 100000;
double x, pi = 0.0, sum = 0.0, step = 1.0/(double)num_steps;
omp_set_num_threads(NUM_THREADS);
#pragma omp parallel for reduction(+:sum) private(x)
for (long i = 0; i < num_steps; i++)
{
x = (i + 0.5)*step;
sum += 4.0/(1.0 + x*x);
75
}
pi += sum*step;
printf("Pi = %f\n", pi);