OpenMP入門

前情提要：並行(parallel):需要多個運算核心同時完成

其中有多處理器和單處理器多核兩種實現方式，其中差異如下：

• 同一晶片上的多核通訊速度更快
• 同一晶片上的多核能耗更低

OpenMP初見

OpenMP環境設定

筆者當初剛進入的時候，走了許多彎路，這裡給出Dev-Cpp版本的OpenMP的環境設定

點選Tools，再點選裡面的Complier Options，在裡面的General下的Add the following commands when calling the compiler:加入這段控制語句-fopenmp就可以在Dev-Cpp中編寫omp程式

OpenMP作用

• 編譯器指令(#pragma)、庫函數、環境變數
• 極大地簡化了C/C++/Fortran多執行緒程式設計
• 並不是全自動並行程式語言(其並行行為仍需由使用者定義及控制)

#pragma

預處理指令

• 設定編譯器狀態或指定編譯器完成特定動作
  • 需要編譯器支援相應功能，否則將被忽略

#pragma once

該段程式碼實現的功能是標頭檔案只被編譯一次

#ifdef HEADER_H

#define HEADER_H

...

#endif

上述程式碼實現的功能也是標頭檔案只被編譯一次

#pragma once	#ifndef
需要編譯器支援	不需要特定編譯器
針對物理檔案	不針對物理檔案
需要使用者保證標頭檔案沒有多份拷貝	需要使用者保證不同檔案的宏名不重複

#pragma GCC poison printf

禁止使用printf函數，該段語句本質上就是申明printf被汙染了，緊跟在該條語句後面的程式碼段不能使用printf。

注意：對於形如#pragma的語句，其作用域一般只包含於該語句的下一個程式碼段，也就是一對大括號之內。

對於#pragma還有許多語法，類似

#pragma startup

#pragma exit

#pragma warning

這些語句對於gcc編譯器來說，不能進行正常的編譯，直接會被忽視，並不會執行，因此這裡不表

總之，對於編譯器來說如果不支援#pragma指令則會直接將其忽視，並不會執行

檢視OpenMP版本

• 使用 _OPENMP 宏定義

使用OpenMP編寫第一個程式

#include<stdio.h>
#include<omp.h>

int main() {
  #pragma omp parallel
  {
    printf("Hello World\n");
  }
  return 0;
}

我們一般通過 #pragma omp parallel 指明並行部分

執行緒的編號獲取是通過下列語句

omp_get_thread_num();

總共的執行緒總數獲取是通過下列語句

omp_get_max_threads();

同一執行緒的多個語句一般情況下並不是連續執行，其中可能有其他執行緒加入，導致該執行緒阻塞。

OpenMP執行機制

• 使用分叉(fork)與交匯(join)模型

fork : 由主執行緒(master thread)建立一組從執行緒(slave threads)

• 主執行緒編號永遠為0(thread 0)
• 不保證執行順序

join : 同步終止所有執行緒並將控制權轉移回至主執行緒

fork是各奔東西，join是百川入海

OpenMP語法

編譯器指令

#pragma omp parallel [clause[clause]...]{stuctured block}

• 指明並行區域及並行方式
• 這裡的clause是指明詳細的並行引數，其大致分為
  • 控制變數線上程間的作用域
  • 顯式指明執行緒數目
  • 條件並行

num_threads(int)

• 用於指明執行緒數目
• 當沒有指明時，將預設使用OMP_NUM_THREADS環境變數
  • 環境變數的值為系統運算核心數目(或超執行緒數目)
  • 可以使用 omp_set_num_threads(int) 修改全域性預設執行緒數，注意這需要在外面修改，線上程中修改筆者嘗試後沒有作用
  • 可以使用omp_get_num_threads()獲取當前設定的預設執行緒數
  • num_threads(int)優先順序高於環境變數，不保證建立指定數目的執行緒，會後到系統資源的限制

這裡做個區分

並行控制指令	功能
omp_get_num_threads	返回當前並行區域中的活動執行緒個數，如果在並行區域外部呼叫，返回1
omp_get_thread_num	返回當前的執行緒號，注意不要和之前的omp_get_num_threads混餚
omp_set_num_threads	設定進入並行區域時，將要建立的執行緒個數
omp_get_max_threads	用於獲得最大的執行緒數量，這個最大數量是指在不使用num_threads的情況下，OpenMP可以建立的最大執行緒數量。需要注意這個值是確定的，與它是否在並行區域呼叫沒有關係

並行for迴圈

#pragma omp parallel
{
  int n;
  for(n = 0; n < 4; n++) {
    int thread = omp_get_thread_num();
    printf("thread %d\n", thread);
  }
}

該段語句是每個執行緒都執行了一次該回圈體中的內容，也就是最後的結果會發現每個執行緒都執行了四次printf("thread %d\n", thread)語句。

如果想要將回圈中的迭代分配到多個執行緒並行，有兩種方式，

方式1

#pragma omp parallel
{
  int n;
  #pragma omp for
  for(n = 0; n < 4; n++) {
    int thread = omp_get_thread_num();
    printf("thread %d]n", thread);
  }
}

在並行區域內加入#pragma omp for
注意：在並行區域內，for迴圈外還可以加入其他並行程式碼

方式2

int n;
#pragma omp parallel for
for(n = 0; n < 4; n++) {
  int thread = omp_get_thread_num();
  printf("thread %d\n", thread);
}

合併為#pragma omp parallel for，這樣子寫法更加簡潔

巢狀並行

OpenMP中的每個執行緒同樣可以被並行化為一組執行緒

• OpenMP預設關閉巢狀，即預設情況下一個執行緒不能並行化為一組執行緒，需要使用omp_set_nested(1)開啟
• 巢狀並行仍然是fork-join模型，注意無論什麼時候，主執行緒的編號一定是0

語法限制

當迴圈分配並行的時候，有下列限制：

• 不能使用!=作為判斷條件(接受<,<=,>,>=)
• 迴圈必須為單入口單出口(也就是不能使用break，goto等跳轉語句)

資料依賴性

• 迴圈迭代相關(loop-carried dependence)

  • 依賴性與迴圈相關，去除迴圈則依賴性不存在

• 非迴圈迭代相關(loop-independent dependence)

  • 依賴性與迴圈無關，去除迴圈依賴性仍然存在

  如果資料與迴圈無關，那麼普通的迴圈並行是可以的，但是如果與迴圈有關，那麼並行就會出現問題

  #include<stdio.h>
  #include<omp.h>
  
  int main() {
    int a[4] = {1, 2, 3, 4};
    #pragma omp parallel for
    for(int i = 1; i < 4; i++) {
      a[i] = a[i-1] + 2;
    }//這裡執行出來的結果不一定是1,3,5,7，這就是迴圈迭代相關
    return 0;
  }
  

執行緒互動與同步

• OpenMP是多執行緒共用地址架構
  • 執行緒可通過共用變數通訊
• 執行緒及其語句執行具有不確定性
  • 共用資料可能造成競爭條件(race condition)
  • 競爭條件:程式執行的結果依賴於不可控的執行順序(此時結果是不可測的，即不同的執行緒執行順序會有不同的結果)
• 必須使用同步機制避免競爭條件的出現
  • 同步機制將帶來巨大開銷
  • 儘可能改變資料的存取方式減少必須的同步次數

高階語言的語句並非原子操作

比如count++並不是只有一個子操作，而是返回count值後再進行++操作。

對於count++其所包含的操作如下

LOAD Reg, count
ADD #1
STORE Reg, count

包含三個子操作，如果一個執行緒負責count++，另一個執行緒負責count--，那麼就有可能出現下面的情況(競爭條件)

        thread 1        thread 2
          Reg    count    Reg
  LOAD     10     10    
  ADD      11     10  
           11     10       10    LOAD
           11     10        9    SUB
           11      9        9    STORE
  STORE    11     11

同樣如果上述的模板中thread 1執行的是count--，而thread 2執行的是count++，那麼最終count儲存的是9，這樣子就會因為執行順序不同而造成不同的結果

臨界區(critical section)

#pragma omp critical

• 指的是一個存取共用資源(例如：共用裝置或是共用記憶體)的程式片段，而這些共用資源又無法同時被多個執行緒存取的特性

  • 同一時間內只有一個執行緒能執行臨界區內程式碼
  • 其他執行緒必須等待臨界區內執行緒執行完畢後才能進入臨界區
  • 常用來保證對共用資料的存取之間互斥

  如果具有競爭條件，那麼就可以使用臨界區來保證結果的正確，此時被#pragma omp critical修飾的程式碼段對於每個執行緒來說此時裡面的資源是它獨佔的，其他執行緒無權存取修改，因此可以消除競爭條件帶來的不確定性

  #pragma omp parrallel for
  for(int i = 0; i < 1000; i++) {
    int value = rand()%20;
    histogram[value]++;
  }
  int total = 0;
  for(int i = 0; i < 20; i++) {
    total +== histogram[i];
    cout << histogram[i] << endl;
  }
  cout << "total:" << total << endl;
  

  本段程式碼的輸出未必是1000，原因可以參考前面的count++，++並非原子操作，因此多執行緒使用的時候會形成競爭條件，導致最後的儲存除了問題

---

解決方案（新增臨界區）：

#pragma omp parallel for
for(int i = 0; i < 1000; i++) {
  int value = rand()%20;
  #pragma omp critical
  {
    histogram[value]++;
  }
}
int total = 0;
for(int i = 0; i < 20; i++) {
  total += histogram[i];
  cout << histogram[i] << endl;
}
cout << "total:" << total << endl;

原子(atomic)操作

#pragma omp atomic

• 保證對記憶體的讀寫更新等操作在同一時間只能被一個執行緒執行
  • 常用來做計數器、求和等
• 原子操作通常比臨界區執行更快
  • 不需要阻塞其他執行緒
• 臨界區的作用範圍更廣，能夠實現的功能更復雜

程式碼演示如下：

#pragma omp parallel for
for(int i = 0; i < 1000; i++) {
  int value = rand()%20;
  #pragma omp atomic
  histogram[value]++;
}

關於原子操作和臨界區的一些總結

指令	功能
臨界區	存取臨界資源的程式碼，一個執行緒進入了臨界區，如果第二個執行緒不會進入臨界區那麼第二個執行緒也可以執行
原子操作	原子操作獨佔處理器，其他執行緒必須等原子操作完了才可以執行

原子操作是山大王，臨界區是此路是我開，此樹是我栽

總的來說就是，原子操作時不允許其他執行緒搶佔資，也就是不能中斷，但是臨界區沒有這個硬性規定，只是相關的程式碼區段此時只能有一個執行緒執行。

柵障(barrier)

#pragma omp barrier

• 在柵障點處同步所有執行緒
  • 先執行至柵障點處的執行緒必須等待其他執行緒
  • 常用來等待某部分任務完成再開始下一部分任務
  • 每個並行區域的結束點預設自動同步執行緒

#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<time.h>
#include<omp.h>

int main() {
  srand(time(NULL));
  int total = 0;
  int histogram[20];
  memset(histogram, 0, sizeof(histogram));
  #pragma omp parallel num_threads(20)
  {
    for(int i = 0; i < 50; i++) {
      int value = rand()%20;
      #pragma omp atomic
      histogram[value]++;
    }

    int thread = omp_get_thread_num();
    #pragma omp atomic
    total += histogram[thread];
  }
  return 0;
}

這裡如果沒有加入柵障指令，那麼最後統計的total一般不是1000，造成資料丟失，因為不能保證執行到total += histogram[thread]對應的數位已經全部統計完畢，其他執行緒可能仍然在統計

修改方案如下：

int total = 0;

#pragma omp parallel num_threads(20)
{
  for(int i = 0; i < 50; i++) {
    int value = rand()%20;
    #pragma omp atomic
    histogram[value]++;
  }
  #pragma omp barrier//保證執行到下一個語句之前，所有的執行緒都已經統計完畢了亂數的次數
  int thread = omp_get_thread_num();
  #pragma omp atomic
  total += histogram[thread];
}

這裡給出#pragma omp for的一些相似之處

int total = 0;

#pragma omp parallel num_threads(20)
{
  #pragma omp for
  for(int i = 0; i < 1000; i++) {
    int value = rand()%20;
    #pragma omp atomic
    histogram[value]++;
  }
  int thread = omp_get_thread_num();
  #pragma omp atomic
  total += histogram[thread];
}

這裡的for分配出來的結果仍然是1000，也就是說#pragma omp for也是等到最後的執行緒執行完迴圈才進行下一步任務，前面先完成的執行緒會進入等待狀態。注意如果修改成#pragma omp for nowait，那麼結果就會不確定，因為此時先完成的執行緒不會再等待了(產生了一個隱式柵障，但是nowait可以去除)。

single & master

#pragma omp single{}

• 用於保證{}內的程式碼由一個執行緒完成
• 常用於輸入輸出或初始化
• 由第一個執行至此處的執行緒執行
• 同樣會產生一個隱式柵障
  • 可由#pragma omp single nowait去除

#pragma omp master{}

• 與single相似，但指明由主執行緒執行
• 與使用if的條件並行等價
  • #pragma omp parallel IF(omp_get_thread_num() == 0) nowait，這邊的話，筆者暫時沒有找到相關的資料，不過可以理解為只有主執行緒執行該程式碼塊，而其他執行緒直接跳過，可以執行後面的任務。
  • 預設不產生隱式柵障

#pragma omp master使用如下

int total = 0;
#pragma omp parrallel
{
  #pragma omp for
  for(int i = 0; i < 1000; i++) {
    int value = rand()%20;
    #pragma omp atomic
    histogram[value]++;
  }  
  #pragma omp master
  {
    for(int i = 0; i < 20; i++) {
      total += histogram[i];
    }
  }
}

注意single會產生柵障，但是master不會

single獨樂樂，master眾樂樂

並行Reduction

指明如何將執行緒區域性結果彙總

• #pragma omp for reduction(+:total)
• 支援的操作，八大金剛
  1. +
  2. -
  3. *
  4. &
  5. |
  6. &&
  7. ||
  8. ^

#include<string.h>
#include<stdio.h>
#include<stdlib.h>
#include<omp.h>
#include<time.h>

int main() {
  int total = 0;
  int histogram[20];
  memset(histogram, 0, sizeof(histogram));
  #pragma omp parallel
  {
    #pragma omp for
    for(int i = 0; i < 1000; i++) {
      int value = rand()%20;
      #pragma omp atomic
      histogram[value]++;
    }
    #pragma omp for reduction(+:total)
    for(int i = 0; i < 20; i++) {
      total += histogram[i];
    }
  }
  printf("%d\n", total);
  return 0;
}

reduction指令就是指明如何將執行緒區域性結果彙總，原先沒有reduction限制的時候，total是每個執行緒共用的，也就是所有的執行緒都可以存取修改total的值。但是加入reduction限制後，每個執行緒的total都是獨屬於該執行緒的全新total，執行緒1和執行緒2的total沒有關係，也就是執行緒1的total的值改變了，不會影響到執行緒2的total的值。這樣每個執行緒都有自己獨特的total，等到所有執行緒結束完成之後，再根據reduction後面所給的操作符，將所有執行緒的total經過該運算儲存到原來的total中。

變數作用域

OpenMP與序列程式的作用域不同

• OpenMP中必須指明變數為shared或private

Shared	Private
變數為所有執行緒所共用	變數為執行緒私有，其他執行緒無法存取
並行區域外定義的變數預設為shared	並行區域內定義的變數預設為private
	迴圈計數器預設為private

int histogram[20];//histogram是shared
memset(histogram, 0, sizeof(histogram));
int total = 0;//total是shared

int i,j;
#pragma omp parallel for
for(i = 0; i < 1000; i++) {//迴圈計數器i是private，雖然定義在並行區域外部
  int value = rand()%20;
  #pragma omp atomic
  histogram[value]++;
  for(j = 0; j < 1000; j++) {//迴圈計數器j也是private
    ...
  }
}

顯式作用域定義

• 顯式指明變數的作用域
• shared(var)
  • 指明變數var為shared
• default(none/shared/private)
  • 指明變數的預設作用域
  • 如果為none則必須指明並行區域內每一變數的作用域

  int a, b, c;
  ...//init a, b, c
  #pragma omp parallel default(shared) shared(a)
  {
    b += a;
  }
  

  這裡會報錯，因為b並沒有指定是private還是shared，注意只要default中選擇的是none，那麼並行體中出現的變數就都要說明是shared還是private型別。
  • 注意在C/C++中只存在none和shared兩種選擇，對於Fortran，可以多一個private引數選擇

  int sum = 0;
  #pragma omp parallel for
  for(int i = 0; i < 20; i++) {
    sum += i;
  }
  

  這裡連續和的答案不確定，也是因為sum是shared，shared本質上就是說所有執行緒的sum指向的都是同一塊記憶體空間
  • 注意：迴圈計數器變數無論前面是shared還是private都是private，否則會造成邏輯上的錯誤
• private(var)
  • 指明變數var為private

  int i = 10;
  #pragma omp parallel for
  for(int j = 0; j < 4; j++) {
    printf("Thread %d: i = %d\n", omp_get_thread_num(), i);
  }
  printf("i = %d\n", i);
  

  這裡執行緒中輸出的i的值都是未定義的，因為此時i是private，也就是說，對於每個執行緒來說，他們的i都是相互獨立的，擁有自己的空間。所以對於每個執行緒，他們的i都需要重新初始化。
• firstprivate(var)
  • 指明變數var為private，同時表明該變數使用master thread中變數值初始化

  int i = 10;
  #pragma omp parallel for firstprivate(i)
  for(int j = 0; j < 4; j++) {
    printf("Thread %d: i = %d\n", omp_get_thread_num(), i);
  }
  printf("i = %d\n", i);
  

  這裡輸出的結果i都是10，也就是說firstprivate(var)的作用就是仍然申明i是private，不過此時不需要我們顯式初始化，而是隱式的初始化，初始化的值是主執行緒中的i的值。
• lastprivate(var)
  • 指明變數為private，同時表明結束將最後一層結果賦予該變數

  int i = 10;
  #pragma omp parallel for lastprivate(i)
  for(int j = 0; j < 4; j++) {
    printf("Thread %d: i = %d\n", omp_get_thread_num(), i);
  }
  printf("i = %d\n", i);
  

  這邊會發現程式執行的結果，i的結果和最後一個執行緒，在這邊是執行緒3中的i一致。所以我們可以得出lastprivate(var)的功能就是將最後一個執行緒中的變數的值在join後賦值給原先主執行緒的。

資料並行和任務並行

資料並行

• 同樣指令作用在不同資料上
• 之前所舉的例子都是資料並行

任務並行

• 執行緒可能執行不同任務
  • #pragma omp section
• 每個section由一個執行緒完成
• 同樣有隱式柵障(可使用nowait去除)

方式1：

#paragma omp parallel

#pragma omp sections
{
  #pragma omp section
  task_A();
  #pragma omp section
  task_B();
  #pragma omp section
  task_C();
}

方式2：

#pragma omp parallel sections
{
  #pragma omp section
  task_A();
  #pragma omp section
  task_B();
  #pragma omp section
  task_C();
}

#pragma omp sections的語法如上，此時裡面的每個#pragma omp section都只會由一個執行緒來執行，如果沒有加nowait，那麼直到最後一個任務完成，所有的執行緒才能繼續往下，否則將在這裡等待，換句話說sections指令有隱式柵障。自測程式如下，加上nowait和不加輸出結果是不同的。

#include<stdio.h>
#include<omp.h>

int main() {
  #pragma omp parallel
  {
  	#pragma omp sections //nowait//這裡證明了sections自帶隱式柵障
  	{
  	  #pragma omp section
	  {
	    int thread = omp_get_thread_num();
		printf("Thread %d : tastk 1\n", thread);
	  }
	  #pragma omp section
	  {
	  	int thread = omp_get_thread_num();
	  	printf("Thread %d : task 2\n", thread);
	  }
	  #pragma omp section
	  {
	  	int thread = omp_get_thread_num();
	  	printf("Thread %d : task 3\n", thread);
	  }
	}
	int thread = omp_get_thread_num();
	printf("Thread %d : task 4\n", thread);
  }
  return 0;
}

執行緒排程

當迭代數多於執行緒數時，需要排程執行緒

• 某些執行緒將執行多個迭代

• #pragma omp parallel for schedule(type,[chunk size])

  • type包括static,dynamic,guided,auto,runtime
  • 預設為static

  排程方式	作用&缺陷
  static	執行緒的負載可能不均勻
  dynamic	在執行中動態分配任務；迭代任務依然根據chunk size劃分成塊；執行緒完成一個chunk後向系統請求下一個chunk
  guided	與dynamic類似；但分配的chunk大小在執行中遞減，最小不能小於chunk size引數
  auto&runtime	編譯器決定策略與環境變數指定策略

總的來說auto和runtime其實最後都是在static,dynamic,guided中挑選一個策略。所以理解這三種排程方式最為重要。

static

• #pragma omp parallel for schedule(static)
  • 沒有指定chunk size的引數，那麼預設就是迭代數除以執行緒總數作為size，然後給每個執行緒分配這麼多的任務
• #pragma omp parallel for schedule(static, 5)
  • 指定了size的大小，那麼就是每次給執行緒的任務大小就是size，注意這邊每個執行緒執行的總任務數仍然是相等的，只不過一次執行的塊的大小變化了(size小的時候)，如果size過大，那麼就會產生某個執行緒執行了非常多的任務，而其他執行緒執行很少或者不執行的現象。

dynamic

• #pragma omp parallel for schedule(dynamic)
  • 沒有指定size，size預設為1，也就是此時每個任務塊蘊含的任務數量為1，然後逐個分配給所有執行緒，執行快的執行緒可以先接收下一個任務，實現任務的動態分配。
• #pragma omp parallel for schedule(dynamic, 5)
  • 指定了size的大小，與不指定的區別只在於此時一個任務塊含有的任務數量為5，每次排程都是5，而不是1

guided

• #pragma omp parallel for schedule(guided)
  • 沒有指定size的情況下，預設為1，與 dynamic類似，不過此時分配的任務數量不是固定不變，而是一開始很多，然後逐漸減小，但是不會小於chunk size的數量
• #pragma omp parallel for schedule(guided, 5)
  • 指定了chunk size的數量，也就是說最後分配的任務數量不能小於這個值，其他的與guided是一樣的

OpenMP小結

軟硬體環境

• CPU多執行緒並行庫
  • 編譯器指令、庫函數、環境變數
• 共用記憶體的多核系統

基礎語法

老祖宗：#pragma omp construct [clause[clause]...]{structured block}

分類	語法
指明並行區域	#pragma omp parallel
迴圈	#pragma omp (parallel) for
巢狀	omp_set_nested(1)
常用函數	omp_get_thread_num(),num_threads(int)
同步	#pragma omp critical,#pragma omp atomic,#pragma omp barrier,nowait
變數作用域	default(none/shared/private),shared(),private(),firstprivate(),lastprivate()
排程	schedule(static/dynamic/guided,[chunk_size])

---

文末資源分享：

OpenMP常見的32個陷阱

OpenMP Reference Guide

禁止轉載，如需轉載請在評論區與筆者商量或者私信。