在C++中使用openmp進行多執行緒程式設計
一、前言
多執行緒在實際的程式設計中的重要性不言而喻。對於C++而言，當我們需要使用多執行緒時，可以使用boost::thread庫或者自從C++ 11開始支援的std::thread，也可以使用作業系統相關的執行緒API，如在Linux上，可以使用pthread庫。除此之外，還可以使用omp來使用多執行緒。它的好處是跨平臺，使用簡單。
在Linux平臺上，如果需要使用omp，只需在編譯時使用"-fopenmp"指令。在Windows的visual studio開發環境中，開啟omp支援的步驟為“專案屬性 -> C/C++ -> 所有選項 -> openmp支援 -> 是(/openmp)”。
本文我們就介紹omp在C++中的使用方法。
二、c++ openmp入門簡介
openmp是由一系列#paragma指令組成，這些指令控制如何多執行緒的執行程式。另外，即使編譯器不支援omp，程式也也能夠正常執行，只是程式不會多執行緒並行執行。以下為使用omp的簡單的例子：
int main()
{
vector<int> vecInt(100);
#pragma omp parallel for
for (int i = 0; i < vecInt.size(); ++i)
{
vecInt[i] = i*i;
}
return 0;
}
12345678910
以上程式碼會自動以多執行緒的方式執行for迴圈中的內容。如果你刪除"#pragma omp parallel for"這行，程式依然能夠正常執行，唯一的區別在於程式是在單執行緒中執行。由於C和C++的標準規定，當編譯器遇到無法識別的"#pragma"指令時，編譯器自動忽略這條指令。所以即使編譯器不支援omp，也不會影響程式的編譯和執行。
三、omp語法
所有的omp指令都是以"#pragma omp“開頭，換行符結束。並且除了barrier和flush兩個指令不作用於程式碼以外，其他的指令都只與指令後面的那段程式碼相關，比如上面例子中的for迴圈。
四、parallel編譯指示
parallel告訴編譯器開始 一個並行塊，編譯器會建立一個包含N（在執行時決定，通常為伺服器的邏輯核數，在Linux上檢視機器的邏輯核數命令為：cat /proc/cpuinfo| grep "processor"| wc -l）個執行緒的執行緒組，所有執行緒都執行接下來的語句或者由”{…}"包含的程式碼塊，在這執行結束之後，又回到主執行緒，建立的這N個執行緒會被回收。
#pragma omp parallel
{
cout << "parallel run!!!\n";
}
1234
以上程式碼在邏輯核數為4的cpu的電腦上執行時，輸出了4行”parallel run!!!"。即編譯器建立了一個包含4個執行緒的執行緒組來執行這段程式碼。在這段程式碼執行結束後，程式執行回到主執行緒。GCC編譯器的實現方式是在內部建立一個函式，然後將相關的執行程式碼移至這個函式，這樣一來程式碼塊中定義的變數成為了執行緒的區域性變數，互不影響。而ICC的實現方式是利用fork()來實現。
執行緒之間共享的變數是通過傳遞引用或者利用register變數來實現同步的，其中register變數在程式碼執行結束之後或者在flush指令呼叫時進行同步。
我們也可以利用if條件判斷來決定是否對後續的程式碼採用並行方式執行，如：
externint parallelism_enabled;
#pragma omp parallel for if(parallelism_enabled)
for(int c=0; c<n;++c)
handle(c);
1234
在這個例子中，如果parallelism_enabled為false，那麼這個for迴圈只會由一個執行緒來執行。
五、for指令
omp中的for指令用於告訴編譯器，拆分接下來的for迴圈，並分別在不同的執行緒中執行不同的部分。如果for指令後沒有緊接著for迴圈，編譯器會報錯。例如，
#pragma omp parallel for
for (int i = 0; i < 10; ++i)
{
printf("%d ", i);
}
12345
以上的程式碼執行後，會打印出[0,9]這10個數字。但是它們的順序是隨機出現的，在我的電腦上，執行的輸出是"0 1 2 8 9 6 7 3 4 5"。事實上，輸出結果不會是完全隨機的，輸出的序列是區域性有序的，因為在編譯器對for迴圈的拆分相當於下面的程式碼：
int this_thread = omp_get_thread_num();
int num_threads = omp_get_num_threads();
int my_start = (this_thread)* 10 / num_threads;
int my_end = (this_thread + 1) * 10 / num_threads;
for (int n = my_start; n < my_end; ++n)
printf("%d ", n);
123456
以上程式碼中，omp_get_thread_num()用於獲取當前執行緒在當前執行緒組中的序號；omp_get_num_threads()用於獲取執行緒組中的執行緒數。所以執行緒組中每個執行緒都運行了for迴圈中的不同部分。
這裡提到for迴圈的不同部分線上程組中的不同執行緒中執行的，執行緒組是在程式遇到"#pragma omp parallel"時建立，在程式塊(parallel後的”{…}"或者語句)結束後，執行緒組中的只有一個主執行緒。因此上面示例中的程式碼事實上是以下程式碼的縮寫：
#pragma omp parallel
{
#pragma omp for
for (int i = 0; i < 10; ++i)
{
printf("%d ", i);
}
}
12345678
此處的"#pragma omp for"即使在“#pragma omp parallel”指令建立的執行緒組中執行的。加入此處沒有#pragma omp parallel指令，那麼for迴圈只會在主執行緒中執行。parallel指令所建立的執行緒組的執行緒數預設是有編譯器決定的，我們也可以通過num_threads指令來指定執行緒數，如#pragma omp parallel num_threads(3)即告訴編譯器，此處需要建立一個包含3個執行緒的執行緒組。
六、Schedule指令
Schedule指令提供對for指令中執行緒排程更多的控制能力。它有兩種排程方式：static和dynamic。
static：每個執行緒自行決定要執行哪個塊，即每個執行緒執行for迴圈中的一個子塊。
dynamic：一個執行緒並不是執行for迴圈的一個子塊，而是每次都向omp執行時庫索取一個for迴圈中的迭代值，然後執行這次迭代，在執行完之後再索取新的值。因此，執行緒有可能執行任意的迭代值，而不是一個子塊。
之前的”#pragma omp parallel for“實際上的效果是”#pragma omp parallel for schedule(static)"。如果我們將之前的示例採用dynamic排程方式，即：
#pragma omp parallel for schedule(dynamic)
for (int i = 0; i < 10; ++i)
{
printf("%d ", i);
}
12345
那麼列印的結果則有可能不是區域性有序的。
在dynamic排程方式中，還可以指定每次索取的迭代值數量。如
#pragma omp parallel for schedule(dynamic，3)
for (int i = 0; i < 10; ++i)
{
printf("%d ", i);
}
12345
在這個例子中，每個執行緒每次都索取3個迭代值。執行完之後，再拿3個迭代值，直到for迴圈所有迭代值都執行結束。在最後一次索取的結果有可能不足3個。在程式內部，上面的例子與下面的程式碼是等價的：
int a,b;if(GOMP_loop_dynamic_start(0,10,1,3,&a,&b)){do{for(int n=a; n<b;++n) printf(" %d", n);}while(GOMP_loop_dynamic_next(&a,&b));}1234567
七、ordered指令
ordered指令用於控制一段程式碼在for迴圈中的執行順序，它保證這段程式碼一定是按照for中的順序依次執行的。
#pragma omp parallel for ordered schedule(dynamic)for (int i = 0; i < 10; ++i){ Data data = ReadFile(files[i]);#pragma omp ordered PutDataToDataset(data);}123456
這個迴圈負責讀取10個檔案，然後將資料放入一個記憶體結構中。讀檔案的操作是並行的，但是將資料存入記憶體結構中則是嚴格序列的。即先存第一個檔案資料，然後第二個…，最後是第十個檔案。假設一個執行緒已經讀取了第七個檔案的，但是第六個檔案還沒有存入記憶體結構，那麼這個執行緒會阻塞，知道第六個檔案存入記憶體結構之後，執行緒才會繼續執行。
在每一個ordered for迴圈中，有且僅有一個“#pragma omp ordered"指令限定的程式碼塊。
八、sections指令
section指令用於指定哪些程式塊可以並行執行。一個section塊內的程式碼必須序列執行，而section塊之間是可以並行執行的。如，
#pragma omp parallel sections{{ Work1();}#pragma omp section{ Work2(); Work3();}#pragma omp section{ Work4();}}123456789
以上程式碼表明，Work1, Work2 + Work3 以及 Work4可以並行執行，但是work2和work3的執行必須是序列執行，並且每個work都只會被執行一次。
九、task 指令(OpenMP 3.0新增)
當覺得for和section指令用著不方便時，可以用task指令。它用於告訴編譯器其後續的指令可以並行執行。如OpenMP 3.0使用者手冊上的一個示例：
struct node { node *left,*right;};externvoid process(node*);void traverse(node* p){if(p->left)#pragma omp task // p is firstprivate by default traverse(p->left);if(p->right)#pragma omp task // p is firstprivate by default traverse(p->right); process(p);}123456789101112
上面的示例中，當處理當前節點process§時，並不能夠保證p的左右子樹已經處理完畢。為了保證在處理當前節點前，當前節點的左右子樹已經處理完成，可以使用taskwait指令，這個指令會保證前面的task都已經處理完成，然後才會繼續往下走，新增taskwait指令之後，以上示例程式碼變為：
struct node { node *left,*right;};externvoid process(node*);void postorder_traverse(node* p){if(p->left)#pragma omp task // p is firstprivate by default postorder_traverse(p->left);if(p->right)#pragma omp task // p is firstprivate by default postorder_traverse(p->right);#pragma omp taskwait process(p);}12345678910111213
以下示例演示瞭如何利用task指令來並行處理連結串列的元素，由於指標p預設是firstprivate方式共享，所以無需特別指定。
struct node {int data; node* next;};externvoid process(node*);void increment_list_items(node* head){#pragma omp parallel{#pragma omp single{for(node* p = head; p; p = p->next){#pragma omp task process(p);// p is firstprivate by default}}}}12345678910111213141516
十、atomic指令
atomic指令用於保證其後續的語句執行時原子性的。所謂原子性，即事務的概念，它的執行不可拆分，要麼執行成功，要麼什麼都沒有執行。例如，
#pragma omp atomic counter += value;12
以上程式碼中，atomic保證對counter的改變時原子性的，如果多個執行緒同時執行這句程式碼，也能夠保證counter最終擁有正確的值。
需要說明的是，atomic只能用於簡單的表示式，比如+=、-=、*=、&=等，它們通常能夠被編譯成一條指令。如果上面的示例改為"counter = counter + value"，那將無法通過編譯；atomic作用的表示式中也不能夠有函式呼叫、陣列索引等操作。另外，它只保證等號左邊變數的賦值操作的原子性，等號右邊的變數的取值並不是原子性的。這就意味著另外一個執行緒可能在賦值前改變等號右邊的變數。如果要保證更復雜的原子性可以參考後續的critical指令。
十一、critical 指令
critical指令用於保證其相關聯的程式碼只在一個執行緒中執行。另外，我們還可以給critical指令傳遞一個名稱，這個名稱是全域性性的，所有具有相同名字的critical相關聯的程式碼保證不會同時在多個執行緒中執行，同一時間最多隻會有一個程式碼塊在執行。如果沒有指定名稱，那系統會給定一個預設的名稱：
#pragma omp critical(dataupdate){ datastructure.reorganize();}...#pragma omp critical(dataupdate){ datastructure.reorganize_again();}123456789
以上程式碼展示的兩個名稱為"dataupdate"的critical程式碼一次只會有一個執行，即datastructure的reorganize()和reorganize_again()不會並行執行，一次最多隻會有一個線上程中執行。
十二、openmp中的鎖
omp執行庫提供了一種鎖：omp_lock_t，它定義在omp.h標頭檔案中。針對omp_lock_t有5中操作，它們分別是：
- omp_init_lock 初始化鎖，初始化後鎖處於未鎖定狀態.
- omp_destroy_lock 銷燬鎖，呼叫這個函式時，鎖必須是未鎖定狀態.
- omp_set_lock 嘗試獲取鎖，如果鎖已經被其他執行緒加鎖了，那當前執行緒進入阻塞狀態。
- omp_unset_lock 釋放鎖，呼叫這個方法的執行緒必須已經獲得了鎖，如果當前執行緒沒有獲得鎖，則會有未定義行為。
- omp_test_lock a嘗試獲取鎖，獲取鎖成功則返回1，否則返回0.
omp_lock_t相當於mutex,如果執行緒已經獲得了鎖，那在釋放鎖之前，當前執行緒不能對鎖進行上鎖。為了滿足這種遞迴鎖的需求，omp提供了omp_nest_lock_t，這種鎖相當於recursive_mutex可以遞迴上鎖，但是釋放操作必須與上鎖操作一一對應，否則鎖不會得到釋放。
十三、flush 指令
對於多執行緒之間共享的變數，編譯器有可能會將它們設為暫存器變數，意味著每個執行緒事實上都只是擁有這個變數的副本，導致變數值並沒有在多個執行緒之間共享。為了保證共享變數能夠線上程之間是真實共享，保證每個執行緒看到的值都是一致的，可以使用flush指令告訴編譯器我們需要哪些哪些共享變數。當我們要在多個執行緒中讀寫共同的變數時，我們都應該使用flush指令。
例如，
/* presumption: int a = 0, b = 0; *//* First thread */ /* Second thread */ b =1; a =1;#pragma omp flush(a,b) #pragma omp flush(a,b)if(a ==0)if(b ==0){{/* Critical section */ /* Critical section */}}12345678
在上面的例子中，變數a,b在兩個執行緒中是共享的，兩個執行緒任何時候看到的a,b的值都是一致的，即執行緒1所見的即使執行緒2所見的。
十四、private, firstprivate,lastprivate 及 shared指令控制變數共享方式
這些指令用於控制變數線上程組中多個執行緒之間的共享方式。其中private,firstprivate,lastprivate表示變數的共享方式是私有的，即每個執行緒都有一份自己的拷貝；而shared表示執行緒組的執行緒訪問的是同一個變數。
私有變數共享方式有三種指令，它們的區別在於：
private：每個執行緒都有一份自己的拷貝，但是這些變數並沒有拷貝值，即如果變數是int,long,double等這些內建型別，那麼這些變數在進入執行緒時時未初始化狀態的；如果變數是類的例項物件，那麼線上程中變數是通過預設構造得到的物件，假設類沒有預設構造，則編譯會報錯，告訴你類沒有可用的預設構造；
firstPrivate：每個執行緒有一份自己的拷貝，每個執行緒都會通過複製一份。如果變數是int,long,double等內建型別則直接複製，如果為類的例項物件，則會呼叫示例物件的拷貝建構函式，這就意味著，假如類是的拷貝構造不可訪問，則變數不能夠使用firstprivate方式共享；
lastprivate:變數在每個執行緒的共享方式與private一致，但不同的是，變數的最後一次迭代中的值會flush會主執行緒中的變數中。最後一次迭代的意思是，如果是for迴圈，則主執行緒的變數的值是最後一個迭代值那次迭代中賦的值；如果是section，則主執行緒的變數最終的值是最後一個section中賦的值。要注意的是，最終主執行緒的中變數的值並非通過拷貝構造賦值的，而是通過operator=操作符，所以如果類的賦值操作符不可訪問，那麼變數不能採用lastprivate方式共享。
十五、default 指令
default命令用於設定所有變數的預設的共享方式，如default(shared)表示所有變數預設共享方式為shared。除此之外，我們可以使用default(none)來檢查我們是否顯示設定了所有使用了的變數的共享方式，如：
int a, b=0;#pragma omp parallel default(none) shared(b){ b += a;}12345
以上程式碼無法通過編譯，因為在parallel的程式碼塊中使用了變數a和b，但是我們只設置了b的共享方式，而沒有設定變數a的共享方式。
另外需要注意的是，default中的引數不能使用private、firstprivate以及lastprivate。
十六、reduction 指令
reductino指令是private,shared及atomic的綜合體。它的語法是：
reduction(operator : list)
其中operator指操作符，list表示操作符要作用的列表，通常是一個共享變數名，之所以稱之為列表是因為執行緒組中的每個執行緒都有一份變數的拷貝，reduction即負責用給定的操作符將這些拷貝的區域性變數的值進行聚合，並設定回共享變數。
其中操作符可以是如下的操作符：
| Operator | Initialization value |
|---|---|
| +，-，|，^，|| | 0 |
| *，&& | 1 |
| & | ～0 |
以下為階乘的多執行緒的實現：
int factorial(int number){int fac =1;#pragma omp parallel for reduction(*:fac)for(int n=2; n<=number;++n) fac *= n;return fac;}12345678
開始，每個執行緒會拷貝一份fac；
parallel塊結束之後，每個執行緒中的fac會利用“*”進行聚合，並將聚合的結果設定回主執行緒中的fac中。
如果這裡我們不用reduction，那麼則需用適用atomic指令，程式碼如下：
int factorial(int number){int fac =1;#pragma omp parallel forfor(int n=2; n<=number;++n){#pragma omp atomic fac *= n;}return fac;}1234567891011
但是這樣一來，效能會大大的下降，因為這裡沒有使用區域性變數，每個執行緒對fac的操作都需要進行同步。所以在這個例子中，並不會從多執行緒中受益多少，因為atomic成為了效能瓶頸。
使用reduction指令的程式碼事實上類似於以下程式碼：
int factorial(int number){int fac =1;#pragma omp parallel{int fac_private =1;#pragma omp for nowaitfor(int n=2; n<=number;++n) fac_private *= n;#pragma omp atomic fac *= fac_private;}return fac;}1234567891011121314
注：最後的聚合實際是包括主執行緒中共享變數的初始值一起的，在階乘的例子中，如果fac的初始值不是1，而是10，則最終的結果會是實際階乘值的10倍！
十七、barrier和nowait指令
barrier指令是執行緒組中執行緒的一個同步點，只有執行緒組中的所有執行緒都到達這個位置之後，才會繼續往下執行。而在每個for、section以及後面要講到的single程式碼塊最後都隱式的設定了barrier指令。例如
#pragma omp parallel{/* All threads execute this. */ SomeCode();#pragma omp barrier/* All threads execute this, but not before * all threads have finished executing SomeCode(). */ SomeMoreCode();}12345678910
nowait指令用來告訴編譯器無需隱式呼叫barrier指令，因此如果為for、section、single設定了nowait標誌，則在它們最後不會隱式的呼叫barrier指令，例如：
#pragma omp parallel{#pragma omp forfor(int n=0; n<10;++n) Work();// This line is not reached before the for-loop is completely finished SomeMoreCode();} // This line is reached only after all threads from// the previous parallel block are finished. CodeContinues(); #pragma omp parallel{#pragma omp for nowaitfor(int n=0; n<10;++n) Work();// This line may be reached while some threads are still executing the for-loop. SomeMoreCode();} // This line is reached only after all threads from// the previous parallel block are finished. CodeContinues();1234567891011121314151617181920212223
十八、single 和 master 指令
single指令相關的程式碼塊只執行一個執行緒執行，但並不限定具體哪一個執行緒來執行，其它執行緒必須跳過這個程式碼塊，並在程式碼塊後wait，直到執行這段程式碼的執行緒完成。
#pragma omp parallel{ Work1();#pragma omp single{ Work2();} Work3();}123456789
以上程式碼中，work1()和work3()會線上程組中所有執行緒都 執行一遍，但是work2()只會在一個執行緒中執行，即只會執行一遍。
master指令則指定其相關的程式碼塊必須在主執行緒中執行，且其它執行緒不必在程式碼塊後阻塞。
#pragma omp parallel{ Work1();// This...#pragma omp master{ Work2();}// ...is practically identical to this:if(omp_get_thread_num()==0){ Work2();} Work3();}12345678910111213141516
十九、迴圈巢狀
如下程式碼並不會按照我們期望的方式執行：
#pragma omp parallel forfor(int y=0; y<25;++y){#pragma omp parallel for for(int x=0; x<80;++x) { tick(x,y); }}123456789
實際上內部的那個“parallel for"會被忽略，自始至終只建立了一個執行緒組。假如將上述程式碼改為如下所示，將無法通過編譯：
#pragma omp parallel forfor(int y=0; y<25;++y){#pragma omp for // ERROR, nesting like this is not allowed. for(int x=0; x<80;++x) { tick(x,y); }}123456789
在OpenMP 3.0中，可以利用collapse指令來解決迴圈巢狀問題，如：
#pragma omp parallel for collapse(2)for(int y=0; y<25;++y){ for(int x=0; x<80;++x) { tick(x,y); }}12345678
collapse指令傳遞的數字就代表了迴圈巢狀的深度，這裡為2層。
在OpenMP 2.5中，我們可以通過將多層迴圈改為單層迴圈的方法來達到目的，這樣便無需迴圈巢狀:
#pragma omp parallel forfor(int pos=0; pos<(25*80);++pos){ int x = pos%80; int y = pos/80; tick(x,y);}1234567
然而重寫這樣的程式碼也並非易事，另一個辦法是採用omp_set_nested方法開啟迴圈巢狀支援，預設是關閉的：
omp_set_nested(1);#pragma omp parallel forfor(int y=0; y<25;++y){#pragma omp parallel for for(int x=0; x<80;++x) { tick(x,y); }}12345678910
現在內層迴圈中，也會建立一個大小為N的執行緒組，因此實際上我們將得到N*N個執行緒，這便會導致頻繁的執行緒切換，會帶來較大的效能損失，這也就是為什麼迴圈巢狀預設是關閉的。也許最好的方法就是將外層迴圈的parallel指令刪除，只保留內層迴圈的parallel：
for(int y=0; y<25;++y){#pragma omp parallel for for(int x=0; x<80;++x) { tick(x,y); }}12345678
二十、取消執行緒（退出迴圈）
假如我們想要使用omp多執行緒來優化如下方法：
constchar* FindAnyNeedle(constchar* haystack, size_t size,char needle){for(size_t p =0; p < size;++p) if(haystack[p]== needle) { return haystack+p; //找到needle，直接退出函式 }return NULL;}123456789
我們最直觀的想法應該是在for迴圈外加上“#pragma omp parallel for"，但是讓人失望的是這將無法通過編譯。因為omp要求必須每個迴圈迭代都能得到處理，因此不允許直接退出迴圈，這也就是說在迴圈中不能使用return、break、goto、throw等能夠中斷迴圈的語句。為了能夠提前退出迴圈，我們需要退出時，通知執行緒組的其他執行緒，讓它們結束執行：
- 棄用omp，選擇其他多執行緒程式設計，如pthread
- 通過共享變數，通知執行緒組其他執行緒
以下為使用布林標記來通知其他執行緒的示例：
constchar* FindAnyNeedle(constchar* haystack, size_t size,char needle){constchar* result = NULL;bool done =false;#pragma omp parallel forfor(size_t p =0; p < size;++p){#pragma omp flush(done)if(!done){/* Do work only if no thread has found the needle yet. */if(haystack[p]== needle){/* Inform the other threads that we found the needle. */ done =true;#pragma omp flush(done) result = haystack+p;}}}return result;}12345678910111213141516171819202122
然而這樣寫有個缺點就是，即使done標記變為true了，其他執行緒仍然需要完成每次迭代，即使這些迭代是完全沒有意義的。. 當然，我們也可以不用上述的done標記：
constchar* FindAnyNeedle(constchar* haystack, size_t size,char needle){constchar* result = NULL;#pragma omp parallel forfor(size_t p =0; p < size;++p)if(haystack[p]== needle) result = haystack+p;return result;}123456789
但是這也並沒有完全解決問題，因為這樣當一個執行緒已經找到需要的結果是，也不能夠避免其他執行緒繼續執行，這也就造成了不必要的浪費。
事實上，omp針對這個問題並沒有很好的解決辦法，如果確實需要，那隻能求助於其他執行緒庫了。
參考資料：
https://blog.csdn.net/gengshenghong/article/details/7003110
https://www.cnblogs.com/liangliangh/p/3565136.html
https://blog.csdn.net/chen134225/article/details/107396923?utm_medium=distribute.pc_relevant.none-task-blog-2defaultbaidujs_title~default-0.essearch_pc_relevant&spm=1001.2101.3001.4242