Ref: https://wdxtub.com/2016/03/20/openmp-guide/
簡介
這門課作為 ECE 中少有的跟計算機科學相關的課，自然是必上不可。不過無論是 OpenMP 還是 CUDA，對於平時極少接觸並行編程的我來說，都是十分吃力的，第一次作業的 OpenMP 編程已經讓意識到了個中的差別，當然，在單個核心的計算速度基本達到極致的現在，掌握並行編程可以算是程序員的基本素養，而 OpenMP 其實是一個非常好的開始，簡單，易懂，見效飛快。所以我們的旅程，就從這里開始吧。
Hello OpenMP
OpenMP是一種面向共享內存以及分布式共享內存的多處理器多線程並行編程語言。一段簡單的代碼如下：
1 #include <omp.h>
2 #include <iostream>
3 using namespace std; 4
5 int main(){ 6 #pragma omp parallel for
7 for (int i = 0; i < 10; ++i) 8 { 9 cout << i; 10 } 11 cout << endl; 12 return 0; 13 }
通過#pragma omp預處理指示符指定要采用OpenMP
通過#pragma omp parallel for來指定下方的for循環采用多線程執行，此時編譯器會根據CPU的個數來創建線程數，對於雙核系統，編譯器會默認創建兩個線程執行並行區域的代碼。
這段程序的輸入如下（省略前面的終端信息）
dawang$ ./a.out
3680479152
dawang$ ./a.out
8603971425
dawang$ ./a.out
3086419752
dawang$ ./a.out
6038714925
常用的庫函數
函數原型 / 功能
返回當前可用的處理器個數
int omp_get_num_procs(void)
返回當前並行區域中的活動線程個數，如果在並行區域外部調用，返回1
int omp_get_num_threads(void)
返回當前的線程號（個人感覺這里為omp_get_thread_ID好一些）
int omp_get_thread_num(void)
設置進入並行區域時，將要創建的線程個數
int omp_set_num_threads(void)
下面的這個例子演示了四個庫函數
1 #include <iostream>
2 #include <omp.h>
3 using namespace std; 4
5 int main(){ 6 cout << "CPU number: " << omp_get_num_procs() << endl; 7
8 cout << "Parallel area 1: " << endl; 9
10 #pragma omp parallel //下面大括號內部為並行區域
11 { 12 cout << "Num of threads is: " << omp_get_num_threads(); 13 cout << "; This thread ID is " << omp_get_thread_num() << endl; 14 } 15
16 cout << "Parallel area 2:" << endl; 17 omp_set_num_threads(4); // 設置為並行區域創建4個線程
18 #pragma omp parallel //下面大括號內部為並行區域
19 { 20 cout << "Num of threads is: " << omp_get_num_threads(); 21 cout << "; This thread ID is " << omp_get_thread_num() << endl; 22 } 23
24 return 0; 25 }
大家可以自己運行一次看看自己的輸出
數據相關性
在循環並行化時，由於多個線程同時執行循環，迭代的順序是不確定的。如果是數據不相關的，則可以采用基本的#pragma omp parallel for預處理器指示符。
如果語句S2與語句S1相關，那么必然存在以下兩種情況之一：
1. 語句S1在一次迭代中訪問存儲單元L，而S2在隨后的一次迭代中訪問同一存儲單元，稱之為循環迭代相關（Loop-Carried Dependence）；
2. S1和S2在同一循環迭代中訪問統一存儲單元L，但S1的執行在S2之前，稱之為非循環迭代相關（Loop-Independent Dependence）。
for 循環並行化的聲明形式
1 #include <iostream>
2 #include <omp.h>
3 using namespace std; 4
5 int main(){ 6 // for 循環並行化聲明形式1
7 #pragma omp parallel
8 { 9 #pragma omp for
10 for (int i = 0; i < 10; ++i){ 11 cout << i << endl; 12 } 13 } 14
15 // for 循環並行化聲明形式2
16 #pragma omp parallel for
17 for (int j = 0; j < 10; ++j){ 18 cout << j << endl; 19 } 20 return 0; 21 }
上邊代碼的兩種聲明形式是一樣的，很顯然第二種聲明形式更為簡潔緊湊。但是第一種聲明形式有一個好處，即可以在並行區域內、for循環以外寫其他並行代碼。
for 循環並行化的約束條件
盡管OpenMP可以方便地對for循環進行並行化，但並不是所有的for循環都可以進行並行化。以下幾種情況不能進行並行化：
1. for循環中的循環變量必須是有符號整形。例如，for (unsigned int i = 0; i < 10; ++i){}會編譯不通過；
2. for循環中比較操作符必須是<, <=, >, >=。例如for (int i = 0; i != 10; ++i){}會編譯不通過；
3. for循環中的第三個表達式，必須是整數的加減，並且加減的值必須是一個循環不變量。例如for (int i = 0; i != 10; i = i + 1){}會編譯不通過；感覺只能++i; i++; –i; 或i–；
4. 如果for循環中的比較操作為<或<=，那么循環變量只能增加；反之亦然。例如for (int i = 0; i != 10; –i)會編譯不通過；
5. 循環必須是單入口、單出口，也就是說循環內部不允許能夠達到循環以外的跳轉語句，exit除外。異常的處理也必須在循環體內處理。例如：若循環體內的break或goto會跳轉到循環體外，那么會編譯不通過。
基本 for 循環並行化舉例
1 #include <iostream>
2 #include <omp.h>
3
4 int main(){ 5 int a[10] = {1}; 6 int b[10] = {2}; 7 int c[10] = {0}; 8
9 #pragma omp parallel
10 { 11 #pragma omp for
12 for (int i = 0; i < 10; ++i){ 13 // c[i] 只跟 a[i] 和 b[i] 有關
14 c[i] = a[i] + b[i]; 15 } 16 } 17
18 return 0; 19 }
嵌套 for 循環並行化舉例
1 #include <omp.h>
2
3 int main(){ 4 int a[10][5] = {1}; 5 int b[10][5] = {2}; 6 int c[10][5] = {3}; 7
8 #pragma omp parallel
9 { 10 #pragma omp for
11 for (int i = 0; i < 10; ++i){ 12 for (int j = 0; j < 5; ++j){ 13 // c[i][j] 只跟 a[i][j] 和 b[i][j] 有關
14 c[i][j] = a[i][j] + b[i][j]; 15 } 16 } 17 } 18
19 return 0; 20 } 21
22 -------------------------------------------------------
23
24 對於雙核 CPU 來說，編譯器會讓第一個cpu完成： 25 for (int i = 0; i < 5; ++i){ 26 for (int j = 0; j < 5; ++j){ 27 // c[i][j] 只跟 a[i][j] 和 b[i][j] 有關
28 c[i][j] = a[i][j] + b[i][j]; 29 } 30 } 31
32 會讓第二個 cpu 完成： 33 for (int i = 5; i < 10; ++i){ 34 for (int j = 0; j < 5; ++j){ 35 // c[i][j] 只跟 a[i][j] 和 b[i][j] 有關
36 c[i][j] = a[i][j] + b[i][j]; 37 } 38 }
數據的共享與私有化
在並行區域中，若多個線程共同訪問同一存儲單元，並且至少會有一個線程更新數據單元中的內容時，會發送數據競爭。本節的數據共享與私有化對數據競爭做一個初步的探討，后續會在同步、互斥相關章節中進行進一步描述。
除了以下三種情況外，並行區域中的所有變量都是共享的：
1. 並行區域中定義的變量
2. 多個線程用來完成循環的循環變量
3. private、firstprivate、lastprivate或reduction字句修飾的變量
例如：
1 #include <iostream>
2 #include <omp.h>
3 using namespace std; 4
5 int main(){ 6 int share_a = 0; // 共享變量
7 int share_to_private_b = 1; // 通過 private 子句修飾該變量之后在並行區域內變為私有變量
8
9 #pragma omp parallel
10 { 11 int private_c = 2; 12
13 #pragma omp for private(share_to_private_b)
14 for (int i = 0; i < 10; ++i) //該循環變量是私有的，若為兩個線程，則一個線程執行0~4，另一個執行5~9
15 cout << i << endl; 16
17 } 18
19 return 0; 20 }
聲明方法 / 功能
1 並行區域中變量val是私有的，即每個線程擁有該變量的一個拷貝 2 private(val1, val2, ...) 3
4 與private不同的是，每個線程在開始的時候都會對該變量進行一次初始化。 5 first_private(val1, val2, ...) 6
7 與private不同的是，並發執行的最后一次循環的私有變量將會拷貝到val 8 last_private(val1, val2, ...) 9
10 聲明val是共享的 11 shared(val1, val2, ...)
如果使用private，無論該變量在並行區域外是否初始化，在進入並行區域后，該變量均不會初始化。
Reduction 的用法
直接上例子
1 #include <iostream>
2 #include <stdio.h>
3 #include <omp.h>
4 using namespace std; 5
6 int main(){ 7 int sum = 0; 8 cout << "Before: " << sum << endl; 9
10 #pragma omp parallel for reduction(+:sum)
11 for (int i = 0; i < 10; ++i){ 12 sum = sum + i; 13 printf("%d\n", sum); 14 } 15
16 cout << "After: " << sum << endl; 17
18 return 0; 19 }
其中sum是共享的，采用reduction之后，每個線程根據reduction（+: sum）的聲明算出自己的sum，然后再將每個線程的sum加起來。
reduction聲明可以看作：
1. 保證了對sum的原則操作
2. 多個線程的執行結果通過reduction中聲明的操作符進行計算，以加法操作符為例：
假設sum的初始值為10，reduction（+: sum）聲明的並行區域中每個線程的sum初始值為0（規定），並行處理結束之后，會將sum的初始化值10以及每個線程所計算的sum值相加。
我們在上邊已經了解了reduction的聲明形式，其具體如下：
reduction (operator: var1, val2, …)
其中operator以及約定變量的初始值如下：
1 運算符 數據類型 默認初始值 2 + 整數、浮點 0
3 - 整數、浮點 0
4 * 整數、浮點 1
5 & 整數 所有位均為1 6 | 整數 0
7 ^ 整數 0
8 && 整數 1
9 || 整數 0
線程同步之 atomic
在OpenMP中，線程同步機制包括互斥鎖同步機制和事件同步機制。互斥鎖同步的概念類似於Windows中的臨界區（CriticalSection）以及Windows和Linux中的Mutex，以及VxWorks中的SemTake何SemGive（初始化時信號量為滿），即對某一塊代碼操作進行保護，以保證同時只能有一個線程執行該段代碼。
atomic（原子）操作語法
1 #pragma omp atomic
2 x< + or * or - or * or / or & or | or << or >> >=expr 3 (例如x <<= 1; or x *=2;) 4
5 或 6
7 #pragma omp atomic
8 x++ //or x--, --x, ++x
可以看到atomic的操作僅適用於兩種情況：
1. 自加減操作
2. x<上述列出的操作符>=expr
例如
1 #include <iostream>
2 #include <omp.h>
3 using namespace std; 4
5 int main(){ 6 int sum = 0; 7 cout << "Before: " << sum << endl; 8
9 #pragma omp parallel for
10 for (int i = 0; i < 20000; ++i){ 11 #pragma omp atomic
12 sum++; 13 } 14 cout << "Atomic-After: " << sum << endl; 15
16 sum = 0; 17 #pragma omp parallel for
18 for (int i = 0; i < 20000; ++i){ 19 sum++; 20 } 21 cout << "None-atomic-After: " << sum << endl; 22 return 0; 23 }
輸出20000。如果將#pragma omp atomic聲明去掉，則輸出值不確定。
線程同步之 critical
這里的臨界區與Windows下的CriticalSection類似。
臨界區聲明方法
1 #pragma omp critical [(name)] //[]表示名字可選
2 { 3 //並行程序塊，同時只能有一個線程能訪問該並行程序塊
4 }
例如
1 #include <iostream>
2 #include <omp.h>
3 using namespace std; 4
5 int main(){ 6 int sum = 0; 7 cout << "Before: " << sum << endl; 8
9 #pragma omp parallel for
10 for (int i = 0; i < 100; ++i){ 11 #pragma omp critical(a)
12 { 13 sum = sum + i; 14 sum = sum + i * 2; 15 } 16 } 17
18 cout << "After: " << sum << endl; 19
20 return 0; 21 }
critical 與 atomic 的區別在於，atomic 僅適用於上一節規定的兩種類型操作，而且 atomic 所防護的僅為一句代碼。critical 可以對某個並行程序塊進行防護。
For a simple increment to a shared variable, atomic and critical are semantically equivalent, but atomic allows the compiler more opportunities for optimisation (using hardware instructions, for example).
In other cases, there are differences. If incrementing array elements (e.g. a[i]++ ), atomic allows different threads to update different elements of the array concurrently whereas critical does not. If there is a more complicated expression on the RHS (e.g. a+=foo() ) then the evaluation of foo() is protected from concurrent execution with critical but not with atomic.
Using a critical section is a legitimate way of implementing atomics inside the compiler/runtime, but most current OpenMP compilers do a better job than this.
線程同步之事件同步機制
互斥鎖同步包括atomic、critical、mutex函數，其機制與普通多線程同步的機制類似。而事件同步則通過nowait、sections、single、master等預處理器指示符聲明來完成。
1. 隱式柵障
barrier為隱式柵障，即並行區域中所有線程執行完畢之后，主線程才繼續執行。
2. nowait 用來取消柵障
其用法如下：
1 #pragma omp for nowait //不能用#pragma omp parallel for nowait
2
3 或 4
5 #pragma omp single nowait
例如
1 #include <stdio.h> 2 #include <omp.h> 3 4 int main(){ 5 #pragma omp parallel 6 { 7 #pragma omp for nowait 8 for (int i = 0; i < 20; ++i){ 9 printf("%d+\n", i); 10 } 11 12 #pragma omp for 13 for (int j = 0; j < 10; ++j){ 14 printf("%d-\n", j); 15 } 16 17 for (int j = 0; j < 10; ++j){ 18 printf("%dx\n", j); 19 } 20 } 21 return 0; 22 }
第一個 for 循環的兩個線程中的一個執行完之后，繼續往下執行，因此同時打印出了第一個循環的 + 和第一個循環的 - 。
可以看到，第二個 for 循環的兩個線程都執行完之后，才開始同時執行第三個 for 循環，並沒有交叉。也就是說，通過 #pragma omp for 聲明的 for 循環結束時有一個默認的柵障。
3. 顯式同步柵障 #pragma omp barrier
1 #include <stdio.h>
2 #include <omp.h>
3
4 int main(){ 5 #pragma omp parallel
6 { 7 for (int i = 0; i < 100; ++i){ 8 printf("%d+\n", i); 9 } 10 #pragma omp barrier
11 for (int j = 0; j < 10; ++j){ 12 printf("%d-\n", j); 13 } 14 } 15 }
兩個線程(具體數目不同 CPU 不同)執行了第一個for循環，當兩個線程同時執行完第一個for循環之后，在barrier處進行了同步，然后執行后邊的for循環。
4. master 通過#pragma omp mater來聲明對應的並行程序塊只由主線程完成
1 #include <stdio.h>
2 #include <omp.h>
3
4 int main(){ 5 #pragma omp parallel
6 { 7 #pragma omp master
8 { 9 for (int j = 0; j < 10; ++j){ 10 printf("%d-\n", j); 11 } 12 } 13
14 printf("This will be shown two or more times\n"); 15 } 16 return 0; 17 }
進入 parallel 聲明的並行區域之后，創建了兩個(或更多)線程，主線程執行了 for 循環，而另一個線程沒有執行 for 循環，而直接進入了 for 循環之后的打印語句，然后執行 for 循環的線程隨后還會再執行一次后邊的打印語句。
5. section 用來指定不同的線程執行不同的部分
通過一個示例說明其使用方法：
1 #include <stdio.h>
2 #include <omp.h>
3
4 int main(){ 5 #pragma omp parallel sections // 聲明該區域分為若干個 section, section 之間的運行順序為並行的關系
6 { 7 #pragma omp section // 第一個 section, 由某個線程單獨完成
8 for (int i = 0; i < 5; ++i){ 9 printf("%d+\n", i); 10 } 11
12 #pragma omp section // 另一個 section, 由某個線程單獨完成
13 for (int j = 0; j < 5; ++j){ 14 printf("%d-\n", j); 15 } 16 } 17 return 0; 18 }
因為並行區域中有兩個線程，所以兩個section同時執行。
線程的調度優化
通過前邊的介紹，知道了並行區域，默認情況下會自動生成與CPU個數相等的線程，然后並行執行並行區域中的代碼，對於並行區域中的for循環，有特殊的聲明方式，這樣不同的線程可以分別運行for循環變量的不同部分。通過鎖同步（atomic、critical、mutex函數）或事件同步（nowait、signal、section、master）來實現並行區域的同步控制。
具體的調度策略均由底層完成，本節介紹幾種可以在上層對for循環進行控制的調度策略。
1 determines which iterations are executed by each thread 2
3 STATIC 4 The iteration space is broken in chunks of approximately size N/(num of threads). Then these chunks are assigned to the threads in a Round-Robin fashion. 5 STATIC, CHUNK 6 The iteration space is broken in chunks of size N. Then these chunks are assigned to the threads in a Round-Robin fashion. 7 Characteristics of static schedules 8 Low overhead 9 Good locality (usually) 10 Can have load imbalance problems 11 DYNAMIC[,chunk] 12 Threads dynamically grab chunks of N iterations until all iterations have been executed. If no chunk is specified, N = 1
13 GUIDED[,chunk] 14 Variant of dynamic. The size of the chunks deceases as the threads grab iterations, but it is at least of size N. If no chunk is specified, N = 1. 15 Characteristics of static schedules 16 Higher overhead 17 Not very good locality (usually) 18 Can solve imbalance problems 19 AUTO 20 The implementation is allowed to do whatever it wishes. (Do not expect much of it as of now) 21 RUNTIME 22 The decision is delayed until the program is run through the sched-nvar ICV. It can be set with: 23 The OMP_SCHEDULE environment variable 24 The omp_set_schedule() API call
能看到這里，如果都跑過一遍的話，應該也就差不多了。上課過程中有啥想法再追加吧。我要去改代碼了再見。