并行程序设计导论学习笔记——OpenMP

Posted 2020-09-19 g63

 使用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。

参考书籍：并行程序设计导论（英文版），机械工业出版社。