OpenMP 学习笔记

更新于 2025-01-04 Disqus：本文字数： 4k 阅读时长 ≈ 20 分钟

想学习多线程 programming 却不知道如何下手？不妨先看看较简单的 OpenMP 吧。

安装 OpenMP

一些不属于 C++ 标准库里的 lib 通常需要手动安装，但是 OpenMP 可不同，目前常见的编译器都自带了它。可以根据下表进行查询：

Compiler	OpenMP Version	Link
Clang	5.0/5.1	clang.llvm.org
GCC	5.0/5.1	gcc.gnu.org
MSVC	2.0 + 一部分 3.1 + 一部分 4.0	learn.microsoft.com
Intel® oneAPI DPC++/C++ Compiler	5.0/5.1/5.2	intel.com
NVIDIA HPC Compilers	5.0/5.1	docs.nvidia.com

当然如果你所使用的编译器不在列表里也不必更换，可以去 openmp.org 进行查询。

虽说大多数编译器已经自带 OpenMP，但一个巨大的例外是 MacOS 自带的 Clang，我们需要手动安装 libomp：

1	brew install libomp

之后还需要设定 OpenMP_ROOT：

1	echo 'export OpenMP_ROOT=$(brew --prefix)/opt/libomp' >> ~/.zshenv

启用 OpenMP

使用这个 section name 实在有些不妥，但实际上的意思是我们如何在拥有一些包含 OpenMP 的代码之后正确的编译并运行他们。

CMake

首先需要找到 OpenMP 这个 package，之后在 linking 的时候链接上就好了，默认链接的是 dynamic lib，也就是说程序的 executable 里并不包含运行所需要的全部文件，运行时还会在 LD_LIBRARY_PATH，rpath/runpath，系统默认的 lib 目录，以及 working directory 下寻找动态库。所以当我们需要把编译后的文件 copy 到其他机器上运行时，别忘了检查 target nodes 上是否具有所有的库文件。

find_package(OpenMP REQUIRED)

# rest code

target_link_libraries(<target> <PRIVATE|PUBLIC|INTERFACE> OpenMP::OpenMP_CXX)

Command-line

根据不同编译器的 syntax，会有些许的不同。

Compiler	Flag
Clang	`-fopenmp`
GCC	`-fopenmp`
MSVC	`/openmp`,`/openmp:experimental`,`/openmp:llvm`
Intel® oneAPI DPC++/C++ Compiler	Linux`-qopenmp`, Windows`/Qopenmp`
NVIDIA HPC Compilers	`-⁠mp`

使用 OpenMP 编写代码

由于 OpenMP 使用 directives 来指导与多线程有关的操作，并不包含如 MPI 中多进程互相发送数据的操作，所以实际上手更加简单。

基本结构

#include <omp.h>

// ...resting code...

#pragma omp construct [clause]

// ...resting code...

比如我们需要创建一个 parallel region 时，需要使用：

#pragma omp parallel
{
    // parallel region
}

当我们需要叠加 constructs 时，可以直接将他们写在同一行：

#pragma omp parallel
{
  #pragma omp for
    // ...
}

与 #pragma omp parallel for 等价。

常用环境变量&函数

Env variables	Explainations	Examples
OMP_NUM_THREADS	设置需要在 `parallel` 中使用的线程数量	`OMP_NUM_THREADS="4"`
OMP_PLACES	设置使用的线程应该分配的地方，是位于不同的 core 上，thread 上，还是不同的 socket 上，抑或是手动控制放置的位置	`OMP_PLACES="cores"`, `OMP_PLACES="{0, 1, 2, 3}"`
OMP_PROC_BIND	设置是否设置 CPU 亲和性，避免执行过程中切换核心	`OMP_PROC_BIND="true"`

1
2
3

void omp_set_num_threads (int num_threads); // 设置需要使用的线程数量，覆盖 OMP_NUM_THREADS env
int omp_get_num_threads (void); // 获取当前 team 下的线程总数
int omp_get_thread_num (void); // 获取当前线程的 id，从 0 开始，0 就是 master thread

`parallel` construct

#pragma omp parallel

在随后（紧跟着）的 block 中创建线程并且每个线程都独立的运行这个 block，在区块结束之后调用 implicit barrier。


graph TD
A([Sequential]) -->|omp parallel| B(Spawn threads)
B -->|Thread 0| C[Block]
B -->|Thread 1| D[Block]
B -->|Other threads| E@{ shape: procs, label: "Block"}
C --> F(omp barrier)
D --> F(omp barrier)
E --> F(omp barrier)
F --> G([Sequential])

#pragma omp parallel // << 在此处生成 OMP_NUM_THREADS 个线程，所有线程均会单独执行下面 block 里的内容
  {
    int totalThreads = omp_get_num_threads();
    int threadId = omp_get_thread_num();
    std::cout << std::format("This is thread {0} out of {1} threads\n",
                             threadId, totalThreads);
  } // << 在 block 结束之后，会等待所有线程执行完毕，再继续执行接下来的部份

运行(OMP_NUM_THREADS=4)后的输出结果（顺序可能会不一样，因为每个线程执行到输出的时间不一定按顺序排列）为:

This is thread 0 out of 4 threads
This is thread 1 out of 4 threads
This is thread 2 out of 4 threads
This is thread 3 out of 4 threads

可以看到，在 block 里的 variables 为每个线程所私有的，因为 threadId 在不同线程中的值是不同的，但是当我们把变量的 declaration 放到 block 外面时，情况有所不同：

  int totalThreads, threadId;
#pragma omp parallel
  {
    totalThreads = omp_get_num_threads();
    threadId = omp_get_thread_num();
    std::cout << std::format("This is thread {0} out of {1} threads\n",
                             threadId, totalThreads);
  }

此时的输出结果顺序则更为混乱：

This is thread 2 out of 4 threads
This is thread 1 out of 4 threads
This is thread 2 out of 4 threads
This is thread 3 out of 4 threads

Master thread 似乎并没有输出？其实不然，此时所有的输出语句对应的 thread 并不一定是 threadId 中所指示的线程，此时我们将变量的 declaration 放到 block 外面之后，他们在 block 创建线程之后为所有线程所共享的，所以一个线程中的 threadId 可以被其他线程所随意修改。

所以我们在设计并行区域时，需要注意 declaration 的位置，需要在多个进程间共享的变量置于外边，将每个线程所私有的变量置于 block 内部。

`for` construct

#pragma omp for

我们最常用的 directives 出现了！他被用在 for-loop 前，用于将一整个循环的工作量分配到不同的线程之中去。


graph TD
A([Sequential]) -->|omp parallel for| B(Spawn threads)
B -->|Thread 0| C[Work 1]
B -->|Thread 1| D[Work 2]
B -->|Other threads| E@{ shape: procs, label: "Work..."}
C --> F(omp barrier)
D --> F(omp barrier)
E --> F(omp barrier)
F --> G([Sequential])

  constexpr size_t ARRAY_SIZE = 1000000;
  auto arr = std::make_unique<int32_t[]>(ARRAY_SIZE);

  for (int32_t i = 0; i < ARRAY_SIZE; ++i) { // 初始化数组
    arr[i] = i;
  }

#pragma omp parallel for // spawn threads
  for (size_t i = 0; i < ARRAY_SIZE; ++i) {
    arr[i] *= 2;
  } // implicit barrier

有关如何在多个线程中分配工作量： Scheduling

`schedule` clause

#pragma omp for schedule([modifier: ]<kind>[, chunk])

这个 directive 是用来辅助 for clause 的，它让我们指定该如何将整个 for loop 的工作量分配到各个线程。

Static

使用其的语法为 static, <chunkSize>。

在这种 scheduler 下，整个 for-loop 的 iterations 将被分为 chunkSize 大小。如此时 iteration 的总数为 $1000$，chunkSize = 100，那么整个循环将被分为 $10$ 个部分（若不能被整除，那么会尽量保证每个部分的大小相近），之后将各部分逐个按顺序分配到线程中。这样的分配过程在循环实际执行前执行。若未指定 chunkSize，此时将整个循环分为大小相近的部分，将每个部分分配到线程中（可以认为此时 chunkSize 为 omp 计算而得来）。

#pragma omp parallel for schedule(static)
  for (size_t i = 0; i < 8; ++i) {
    std::cout << std::format("Thread: {} - i: {}\n", omp_get_thread_num(), i);
  }

对应的输出为：

Thread: 0 - i: 0
Thread: 1 - i: 2
Thread: 2 - i: 4
Thread: 3 - i: 6
Thread: 3 - i: 7
Thread: 1 - i: 3
Thread: 2 - i: 5
Thread: 0 - i: 1

则每个线程均被分配到 $2$ 次循环，且此时每个线程所执行的为 adjacent iterations，也就是 $0$ 号线程执行 $0-1$，$1$ 号线程执行 $2-3$，以此类推。

当我们指定了 chunkSize 后，分配的方法会发生变化，若 $\text{iterations} / \text{chunkSize} > \text{threadNum}$，那么将会将剩下的任务重新从 $0$ 号线程开始进行分配（就是 round-robin fashion，相信这个词可以更好的表达意思🤓）。

#pragma omp parallel for schedule(static, 1)
  for (size_t i = 0; i < 8; ++i) {
    std::cout << std::format("Thread: {} - i: {}\n", omp_get_thread_num(), i);
  }

可以发现此时 $0$ 号线程执行的是 $0\text{th}$ 与 $4\text{th}$ 循环：

Thread: 0 - i: 0
Thread: 2 - i: 2
Thread: 1 - i: 1
Thread: 3 - i: 3
Thread: 0 - i: 4
Thread: 2 - i: 6
Thread: 1 - i: 5
Thread: 3 - i: 7

请不要将 schedule(static) 视为 schedule(static, 1)，它们对任务的分配并不相同！

若我们执行下段代码：

#pragma omp parallel for schedule(static, 4)
  for (size_t i = 0; i < 8; ++i) {
    std::cout << std::format("Thread: {} - i: {}\n", omp_get_thread_num(), i);
  }

则会发现只有 $0$ 与 $1$ 号线程各自接收到 $4$ 次循环：

Thread: 0 - i: 0
Thread: 1 - i: 4
Thread: 0 - i: 1
Thread: 0 - i: 2
Thread: 0 - i: 3
Thread: 1 - i: 5
Thread: 1 - i: 6
Thread: 1 - i: 7

由此我们可以发现，使用 static scheduler 的优势是具有 predictability，可以让我们提前知道每个线程将会被分配到多少工作量以及每个线程执行的循环具体为哪一些。其次他的分配过程开销较小，因为此过程在循环进行前就执行了。但此调度的缺点也显而易见，，如果每次 iteration 所需要的时间差距较大，就会出现“一核有难，九核围观”的名场面（毕竟整个循环结束后会有 implicit barrier）。

Dynamic

使用其的语法为 dynamic, <chunkSize>。

与 static 类似，此调度也是将所有 iterations 分配为大小相近的几个部分。但不同的是，他并不会在循环执行前将不同部分按顺序分配给不同线程，而是根据线程的运行状况来进行分配，一个线程结束当前所分配的任务后，他就会被分配新任务以执行（当一个线程 idle 时，它会请求新的大小为 chunkSize 的任务）。

#pragma omp parallel for schedule(dynamic)
  for (size_t i = 0; i < 8; ++i) {
    std::cout << std::format("Thread: {} - i: {}\n", omp_get_thread_num(), i);
  }

Thread: 0 - i: 0
Thread: 1 - i: 2
Thread: 0 - i: 1
Thread: 3 - i: 6
Thread: 1 - i: 3
Thread: 2 - i: 4
Thread: 0 - i: 5
Thread: 3 - i: 7

则说明 $0$ 号线程执行之间更快，所以它执行了 $3$ 次循环。

让我们在循环中添加一些小把戏来更显著的查看此差异：

#pragma omp parallel for schedule(dynamic)
  for (size_t i = 0; i < 8; ++i) {
    if (omp_get_thread_num() == 0) {
      std::this_thread::sleep_for(std::chrono::seconds(1));
    }

    std::cout << std::format("Thread: {} - i: {}\n", omp_get_thread_num(), i);
  }

Thread: 1 - i: 2
Thread: 2 - i: 4
Thread: 3 - i: 6
Thread: 1 - i: 3
Thread: 2 - i: 5
Thread: 3 - i: 7
Thread: 1 - i: 1
Thread: 0 - i: 0

此时由于 $0$ 号线程的执行时间过长，导致其他线程把剩下的循环给解决了。

使用这种调度方式的优势是 better load balancing，每个线程可以在执行完当前任务后执行新的任务，从而降低某些耗时循环在同一线程上执行的概率。但此种方法由于是在循环运行过程中动态决定任务分配，具有较大的 overhead，同时我们也不可预测每个线程所执行的任务。

此 scheduler 更适合 unpredictable workload per iteration。

Guided

有没有一种折中的方案呢？像 TCP congestion control 那样可以动态调整的 chunkSize 有木有🤔？

🤓☝️还真有，下面隆重介绍 schedule(guided, <chunkSize>)。:D

“~~Slow~~Fast start”：在一开始时内部的 chunkSize 将会是一个较大的值，来降低 overhead，但是不会小于手动设定的的大小（除非手动设定的太大了）。

“Adaptability”：随着线程们执行完各自的任务，新的 chunkSize 将会慢慢变小，这更有效的保证了每个线程都有事可做。

这种调度方法更适合特别需求每个线程都保持 busy 的环境，但其代价则为最高的 overhead，因为不但要在循环进行过程中进行调度，还要计算调度所需要的信息。

Auto & runtime

设置为 auto 则表明让编译器/运行时自行决定使用何种调度方法。

runtime 则会获取环境变量 OMP_SCHEDULE 的值（例如 OMP_SCHEDULE="nonmonotonic:dynamic,4"）并将其作为调度器。

Modifier

在类型之前还可以增加 modifier 来控制 iterations 是如何在线程间进行分配的。

Modifier	Explaination
`monotonic`	Default behavior of `static`，对于分配到每个线程的 chunk，均会按照先后顺序执行，例如 $0-100$ 分配给 $0$ 号线程后，会按照 $0, 1, 2, \cdots$ 的顺序执行
`nonmonotonic`	Default behavior but `static`，对于分配到的 chunk，执行的顺序并不保证，例如 $0-100$ 分配给 $0$ 号线程后，可能会按照 $0, 10, 4, \cdots$ 的顺序执行
`simd`	当此 modifier 被运用在 SIMD consurtct 里时，omp 会执行类似 SIMD 的 fashion，将循环展开而更好的利用处理器的特性，`chunkSize` 会因此而改变为 $\lceil \text{chunkSize}/\text{simdWidth} \rceil \times \text{simdWidth}$。当此 modifier 并没有被运用在 SIMD construct 里时会被忽略

`simd` construct

`#pragma omp simd`

常被用来表明接下来的循环可以被 vectorize，用一次向量化操作同时执行多个 iterations。

例如下段代码：

constexpr size_t N = 1e9;

auto a = std::make_unique<double[]>(N);
auto b = std::make_unique<double[]>(N);
auto c = std::make_unique<double[]>(N);

for (size_t i = 0; i < N; ++i) { // initialization
  a[i] = 1.0;
  b[i] = 2.0;
}

for (size_t i = 0; i < N; ++i) { // calculation
  c[i] = a[i] * b[i];
}

我们使用类似 VTune 这样的 perf analysis 工具进行 profile 之后，可以发现默认情况下编译器生成的为 scalar 指令：

Scalar program

而当我们加入 simd construct 之后，编译器就会生成使用向量指令的程序：

  constexpr size_t N = 1e9;

  auto a = std::make_unique<double[]>(N);
  auto b = std::make_unique<double[]>(N);
  auto c = std::make_unique<double[]>(N);

#pragma omp simd
  for (size_t i = 0; i < N; ++i) { // initialization
    a[i] = 1.0;
    b[i] = 2.0;
  }

#pragma omp simd
  for (size_t i = 0; i < N; ++i) { // calculation
    c[i] = a[i] * b[i];
  }

Vector program

当我们在编译指令中添加一些更高级的向量化命令（如 -mavx2）后，omp 会尝试使用更高阶的向量化指令，若我们什么都不添加，即使 CPU 支持 AVX2 或 AVX512， omp 也只会采用 128bit 向量。

需要注意的是，只使用 simd construct 时，仅会有单个线程来执行整个循环，若需要多线程同时执行且包含 SIMD 指令时，请查看 omp for simd

`#pragma omp for simd`

当我们了解了 for 与 simd construct 之后，两者合并起来便成为了在多个线程中均使用 SIMD 指令。

  constexpr size_t N = 1e9;

  auto a = std::make_unique<double[]>(N);
  auto b = std::make_unique<double[]>(N);
  auto c = std::make_unique<double[]>(N);

#pragma omp parallel for simd
  for (size_t i = 0; i < N; ++i) { // initialization
    a[i] = 1.0;
    b[i] = 2.0;
  }

#pragma omp parallel for simd
  for (size_t i = 0; i < N; ++i) { // calculation
    c[i] = a[i] * b[i];
  }

这样就可以将整个 for loop 的工作量分给不同线程的同时，利用 CPU 的向量化能力，进一步提升程序性能（其实还是挺明显的，具体可以查看 SIMD 的学习笔记）。

`sections` construct

#pragma omp sections

当我们需要在非循环中分配任务时应该如何？sections construct 给予我们更大的权限来指定每个线程所做的事：

#pragma omp parallel sections
  {
#pragma omp section
    {
      // Code block for section 1
    }
#pragma omp section
    {
      // Code block for section 2
    }
    // Additional sections can be added here
  }

我们可以在每个 omp section 后的 block 中编写希望某个线程所作的事，在实际运行过程中，每个 block 都会被运行一次，不同的 block 同时运行在不同的 thread 之上。

`single` construct

#pragma omp single

设置随后的 block 仅由其中任意一个 thread 来执行，但在 block 结束之后会有 implicit barrier。当我们的 parallel block 中出现一些不应该被许多线程同时执行，或不能被执行多次的 section 时（例如 I/O operations），就可以使用这个 construct。

  int sharedData = 0;

#pragma omp parallel
  {
#pragma omp single
    { // only one thread will execute this block
      std::this_thread::sleep_for(std::chrono::seconds(1));

      sharedData = 42;
      std::cout << std::format("Thread {} sets shared_resource = {}\n",
                               omp_get_thread_num(), sharedData);
    } // << implicit barrier

    // all threads will execute this block
    std::cout << std::format("Thread {} reads shared_resource = {}\n",
                             omp_get_thread_num(), sharedData);
  }

Thread 0 sets shared_resource = 42
Thread 0 reads shared_resource = 42
Thread 3 reads shared_resource = 42
Thread 2 reads shared_resource = 42
Thread 1 reads shared_resource = 42

`masked` construct

#pragma omp masked

当我们想指定某个 section 在一些特定的 thread （默认为 master 也就是 thread $0$ ）上运行时，可以使用 masked construct。对于 initialization 与 finalization 有所帮助。但在 masked 的 block 之后并没有 implicit barrier。

  int sharedData = 0;

#pragma omp parallel
  {
#pragma omp masked
    {
      std::this_thread::sleep_for(std::chrono::seconds(1));

      sharedData = 42;
      std::cout << std::format("Thread {} sets shared_resource = {}\n",
                               omp_get_thread_num(), sharedData);
    } // << no barrier

    // all threads will execute this block without waiting for masked block
    std::cout << std::format("Thread {} reads shared_resource = {}\n",
                             omp_get_thread_num(), sharedData);
  }

Thread 1 reads shared_resource = 0
Thread 2 reads shared_resource = 0
Thread 3 reads shared_resource = 0
Thread 0 sets shared_resource = 42
Thread 0 reads shared_resource = 42

`critical` construct

#pragma omp critical

当我们不想让某 section 同时被多个线程所执行，但又想让所有线程均可以执行该部分，就可以使用 critical construct，保证当前最多有一个 thread 在执行这个 block。

  int sharedData = 0;

#pragma omp parallel
  {
    for (int i = 0; i < 1000; ++i) {
#pragma omp critical
      {
        // Only one thread will execute this block at a time
        sharedData++;
      }
    }
  }

  std::cout << std::format("sharedData = {}\n", sharedData);

1	sharedData = 4000

`atomic` construct

`task` construct

`taskloop` construct

`target` construct

`teams` construct

`distribute` construct

`barrier` construct

`flush` construct

`ordered` construct

`cancel` construct

References

MacOS 安装 OpenMP
IBM
OpenMP Tutorial
OpenMP - Scheduling

安装 OpenMP

启用 OpenMP

CMake

Command-line

使用 OpenMP 编写代码

基本结构

常用环境变量&函数

parallel construct

for construct

schedule clause

Static

Dynamic

Guided

Auto & runtime

Modifier

simd construct

#pragma omp simd

#pragma omp for simd

sections construct

single construct

masked construct

critical construct

atomic construct

task construct

taskloop construct

target construct

teams construct

distribute construct

barrier construct

flush construct

ordered construct

cancel construct