OpenMP Tutorial

发表于 更新于 分类于 本文字数： 7.1k 阅读时长 ≈ 12 分钟

概述

   OpenMP 可以认为是用于共享内存型并行编程 的 API 接口，“MP” 是多处理（”multiprocessing”）的缩写。因此，OpenMP 被用于这样的系统：任何进程或线程都可以潜在地访问所有可用的内存资源。当我们使用 OpenMP 编程时，我们可以将系统看作 CPUs 的集合，这些 CPU 都可以访问内存。

参考 《Introduction to ParallelProgramming》 by Peter Pacheco.

OpenMP Vs Pthreads

   尽管 Pthreads 也用于共享内存式的编程，但它和 OpenMP 有些不同：

  • Pthreads 要求用户显示定义每个线程的行为；OpenMP 只用声明目标代码块需要并行，具体的行为交给编译器和运行时环境完成。
  • Pthreads 可以作为与主程序进行链接的库，使用 C 编译器编译即可；但 OpenMP 需要编译器的额外支持，因此使用不支持 OpenMP 的 C 编译器来编译 OpenMP 程序不会获得预想的并行执行行为。
  • Pthreads 更加底层，赋予用户控制线程行为细节的能力，但也给用户增加了编程负担；OpenMP 使得并行编程更加简单，但控制底层线程行为的能力较弱。

  对于大型的高性能程序，学术界和工业界的很多人都相信使用 Pthreads API 是非常困难的。OpenMP 被显式设计为一种更高层的 API，使得用户可以增量式地将串行程序并行化。

简单使用

   OpenMP 使用 “基于声明式” 的API，在C/C++语言中表现为使用预处理指令 #pragma。 #pragma 所定义的很多指示字是编译器特有的，在不同的编译器间是不可移植的。预处理期将忽略它不认识的 #pragma 指令，不同的编译器可能以不同的方式解释同一条 #pragma 指令。对于 OpenMP, 使用 #pragma omp 作为声明的前缀。

一个简单的 HelloWorld 程序 omp_hello.c：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
#include <stdio.h>
#include <stdlib.h>

#ifdef _OPENMP
    #include <omp.h>
#endif

void hello() {

#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

    printf("Hello from thread %d of %d\n", my_rank, thread_count);
}

int main(int argc, char *argv[])
{
    
    int thread_count = strtol(argv[1], NULL, 10);

    #pragma omp parallel num_threads(thread_count)
    hello();

    return 0;
}

带有 -fopenmp 选项进行编译：

$ gcc -g -Wall -fopenmp -o omp_hello omp_hello.c

运行（指定4个线程）：
$ ./omp_hello 4

输出结果：

1
2
3
4
Hello from thread 0 of 4
Hello from thread 3 of 4
Hello from thread 2 of 4
Hello from thread 1 of 4

由于 4 个线程竞争使用 stdout ，所以上述输出顺序不一定与线程的编号一致，当再次执行时，输出顺序就可能不同。

使用 single 声明可以指定代码只能被一个线程执行。

在上述代码中，#pragma omp声明中 parallel 指令表示紧跟的结构块（structured block）需要被多个线程执行。结构块可以是单个语句或者只有一个入口和出口组合语句（类似于编译理论中的“基本块”），即它不包括分支语句，也不是任何分支语句的目标。

num_threads 被称为子句（clause），它可以修改 parallel 指令的默认行为。

上述声明所代表的底层行为到底是怎么样的呢？
首先当主线程（Master）运行到 parallel 指令时，会生成 thread_cnt-1 个从线程（Slave），然后每个线程都会执行结构块中的代码，最后当某个线程执行完成后，会等待其他线程完成。当所有线程完成后，从线程会被终结，主线程继续执行结构块之后的代码。Master线程和所有Slave线程组成一个 team 。

一个复杂例子

   该节以计算某个函数 f 在区间 [a,b] 上的积分为例。我们使用梯形公式：
令 $ h = \frac{(b-a)}{n} $ （n 代表将定义域分为等长区间的个数）：
$$ S = h[\frac{f(x_0)}{2} + f(x_1) + f(x_2)+\cdots+\frac{f(x_n)}{2}] $$

calc 函数代表每个线程需要执行的结构块，它会根据当前线程的编号执行对应的子任务，其结果被写入 group_result_ptr 变量。group_result_ptr 变量必须进行互斥访问，可以使用 #pragma omp critical 来声明。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
int f(int x) {
    return x* x;
}
void calc (double a, double b, int n, double *global_result_p){
    double h, x, my_result;
    double local_a, local_b;
    int i, local_n;

    int my_rank = omp_get_thread_num();
    int thread_count = omp_get_num_threads();

    h = (b-a) / n;
    local_n = n / thread_count;
    local_a = a + my_rank*local_n*h;
    local_b = local_a + local_n*h;
    
    my_result = (f(local_a) + f(local_b)) / 2.0;
    
    for(i = 1; i <= local_n - 1; i++) {
        x = local_a  + i * h;
        my_result += f(x);
    }

    my_result = my_result*h;

#pragma omp critical
    *global_result_p += my_result;
}

我们定义 f(x) = x*x 来作为待积分的函数，main 函数如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
int main(int argc, char *argv[])
{
    double global_result = 0.0;
    double a, b;
    int n;
    int thread_count;
    
    thread_count = strtol(argv[1], NULL, 10);
    printf("Enter a, b, and n: \n");
    scanf("%lf %lf %d", &a, &b, &n);

#pragma omp parallel num_threads(thread_count)
    calc(a, b, n, &global_result);

    printf("With n = %d trapezoids, our estimate\n", n);
    printf("of the integral from %f to %f = %.14e\n",a, b, global_result);
    
    return 0;
}

main 函数与第2节中的代码基本相同，这里calc函数被声明需要并行执行。

编译后执行：

1
2
3
4
5
6
$ ./trap 4
Enter a, b, and n: 
3 6 100
With n = 100 trapezoids, our estimate
of the integral from 3.000000 to 6.000000 = 6.30004500000000e+01

精确值为 63，上面的输出结果是合适的，如果增加n值，会得到更精确的积分结果。

变量的作用域

还原（reduction）、还原变量及还原操作符（reduction operator）
还原子句： reduction(<operator> : <variable list>)

使用 private(var1,var2,...) 子句可以声明某些变量位于私有域中。
使用 shared(var1, var2,...) 子句可以声明某些变量位于共享域中。

使用parallel for 指令

   parallel for指令可用于声明 for 循环为并行执行，上文的积分程序可以不使用 parallel 指令，而使用 parallel for指令，如下所示：

1
2
3
4
5
6
7
h = (b-a)/n;
approx = (f(a) + f(b))/2.0;
#pragma omp parallel for num_threads(thread_count) reduction(+:approx)
for (i=1; i<= n-1; i++) {
    approx += f(a+i*h);
}
approx += h* approx;

值的注意的是，approx 必须被声明为还原变量，+为还原运算符。另外，在使用parallel声明时，每个线程所要完成的子任务由线程自己分配，但使用parallel for声明时，任务的分配由系统决定。

for 循环中的 i 位于私有域中，即每个线程有一个 i 的拷贝，如果位于共享域，i++语句的执行 会导致数据竞争。对于被parallel for声明的for循环，其中的所有变量都位于共享域，但循环变量除外，例如这里的 i 。

可以使用parallel for 声明的 for 循环是有限制的，它需要满足一定的规范形式。否则的话，运行时系统如何在执行之前知道for循环的迭代次数呢？

另外，非常重要的两点是：

  • 对于使用 parallel for声明的for循环，OpenMP 编译器不会检测循环中不同迭代间的依赖关系，这些依赖必须由编程人员自己来检查。
  • 所以，对于不同迭代间存在依赖的 for循环，通常不能被 OpenMP 并行化。

使用 default(none) 子句可以强制用户为每个变量声明域，例如循环变量，在前文我们知道，循环遍历默认是私有的，但使用default(none)子句后，我们必须使用private显式声明它是私有的。

并行化排序程序

循环调度

  在使用parallel for声明 for 循环时，任务分配是由系统决定的，系统可以使用简单的块划分方法：前n/thread_count次迭代分配给线程0，接着的 n/thread_count 次迭代分配给线程1，依次类推，直到分配完所有迭代。但有时候这并不是最优的分配方案，比如当不同迭代的负载并不相同的时候。

OpenMP 的 schedule子句可以用来调度迭代，来获得性能提升。格式为：

schedule(<type> [, <chunksize>])

<type>可以是下面的任意一种：

  • static：在循环执行前分配迭代。
  • dynamic 或 guided：循环执行期间分配迭代，即一个线程完成它的任务后，还可以从运行时系统获取更多的任务。
  • auto：由编译器或者运行时系统决定调度。
  • runtime：由运行时系统决定调度。

对于 static 调度，使用 round-robin 分配方法，如果我们有12个迭代和3个线程，schedule(static, 1)的分配方案如下：

  • Thread0 ：0，3，6，9
  • Thread1 ：1，4，7，10
  • Thread2 ：2，5，8，11

<chunksize>默认为total_iterations/thread_count，在这种情况下的分配方案几乎等价于默认的分配方案（即简单的块划分）。

对于 dynamic调度，<chunksize>默认为 1 。对于 guided调度，当某个线程完成它的任务时，重新获得的迭代个数会递减。

对于 runtime 调度，具体的调度策略由环境变量 OMP_SCHEDULE 决定。

那么应该使用哪一种调度策略呢？

schedule子句是存在开销（overhead）的，且guided调度的开销最大，dynamic次之，static最小。如果不使用 schedule子句时的程序性能是可以接受的，就不要使用它。下面是一些总结：

  • 如果每次迭代的负载基本相同，就不要使用schedule子句。
  • 如果每次迭代的负载是线性递增的，使用schedule(static,1) 或许可以获得更好的性能。
  • 如果迭代的负载不确定，可以使用不同的调度选项进行探测，这是可以使用 schedule(runtime)，每次探测可以对 OMP_SCHEDULE环境变量设定不同的值。

同步

   pragma omp barrier 可以在 parallel结构块的中间显式声明同步语义，即 team 中的任一线程必须在屏障点等待其他所有线程也到达屏障点之后才能继续执行。另外，parallel、parallel for 和 for 指令在结构块的结尾加入了一个隐含的屏障。对于一些 OpenMP 实现，支持使用nowait声明来取消这个隐式的屏障。

   pragma omp atomic 指令类似 critical 但拥有更好的性能。使用 atomic声明的临界区只能是单语句：

1
2
3
4
5
x <op>= <expression>
x++;
++x;
x--;
--x;

另外，critical(name) 与 critical不同，它可以声明不同的临界区（以不同的 name 命名），这些临界区的代码可以并行执行。

锁（locks）：被用于数据结构而非代码块的互斥访问。

1
2
3
4
void omp_init_lock(omp_lock_t *);
void omp_set_lock(omp_lock_t *);
void omp_unset_lock(omp_lock_t *);
void omp_destory_lock(omp_lock_t *);

另外需要注意的是：

  • 对于单个临界区不应该混合使用不同类型的同步原语。
  • 不保证线程一定不会饥饿。
  • 嵌套使用同步原语是“危险”的，有可能造成死锁。

Cache一致性

线程安全

   在结构块中调用外部函数时，要清楚地知道该函数是否线程安全（thread-safe），如果不是线程安全的，则该结构块不能被并行化，除非我们可以使用线程安全的函数替换它。

其他 directive 和 clause

directive

  • atomic : 内存位置将会原子更新（Specifies that a memory location that will be updated atomically）。
  • flush : 所有线程对所有共享对象具有相同的内存视图（view of memory）。
  • master : 指定由主线程来执行接下来的程式。
  • ordered: 指定在接下来的代码块中，被并行化的 for 循环将依序执行（sequential loop）。
  • sections : 将接下来的代码块包含将被并行执行的section块。
  • single : 之后的程式将只会在一个线程（未必是主线程）中被执行，不会被并行执行。
  • threadprivate : 指定一个变量是线程局部存储（thread local storage）。

clause