我們目前的計算機都是基於馮偌伊曼結構的，在MIMD作為主要研究對象的系統中，分為兩種類型：共享內存系統和分布式內存系統，之前我們介紹的基於MPI方式的並行計算編程是屬於分布式內存系統的方式，現在我們研究一種基於OpenMP的共享內存系統的並行編程方法。OpenMP是一個什么東東？首先我們來看看來之百度百科中的定義：OpenMp是由OpenMP Architecture Review Board牽頭提出的，並已被廣泛接受的，用於共享內存並行系統的多處理器程序設計的一套指導性的編譯處理方案(Compiler Directive)。OpenMP支持的編程語言包括C語言、C++和Fortran；而支持OpenMp的編譯器包括Visual studio，Sun Compiler，GNU Compiler和Intel Compiler等。OpenMP功能中最強大的一個功能是：在我們之前串行程序的源碼基礎上，只要進行少量的改動，就可以並行化許多串行的for循環，達到明顯提高性能的效果。
1.OpenMP環境配置與例子
廢話不多少，首先上一個例子，由於本文以Vs 2015作為IDE，C++作為開發語言，在正式進行OpenMP編碼之前，需要對編譯器稍微配置一下。啟動VS2015新建一個C++的控制台應用程序，如下圖所示：
然后在項目解決方案資源管理器上選擇項目名稱，點擊右鍵，選擇“屬性”，如下圖所示：
然后在屬性頁上左側選擇“配置屬性”——“C/C++”——“語言”，然后在右側“OpenMP支持”后選擇“是（/openmp）”，如下圖所示：
這樣就能在我們VS2015中進行OpenMP的程序開發了。將如下串行程序運行一下
//數字累加 void blog4::SumForNumber() { int a = 0; for (int i = 0; i<100000000; i++) a++; } //測試串行程序 void blog4::Test1(int argc, char* argv[]) { clock_t t1 = clock(); for (int i = 0; i<100; i++) SumForNumber(); clock_t t2 = clock(); std::cout << "Serial time: " << t2 - t1 << std::endl; }對上面的串行程序增加一句話：#pragma omp parallel for。變為OpenMP的並行，繼續運行。
void blog4::Test2(int argc, char* argv[]) { clock_t t1 = clock(); #pragma omp parallel for for (int i = 0; i<100; i++) SumForNumber(); clock_t t2 = clock(); std::cout << "Parallel time: " << t2 - t1 << std::endl; }同兩次運行時間的對比圖如下，可以看出，只要對程序做很少的改動，程序的運行時間就會減少為之前的1/3，性能提升有3倍左右，效果還是蠻好的。
2.基礎理論
通過上面的例子，相信大部分同學已經了解了怎么寫一個簡單的OpenMP的並行計算程序，已經注意到OpenMP的並行計算程序總是以
#pragma omp作為開始的。
在pragma opm后面是一條parallel指令，用來表明之后是一個結構化代碼塊，最基本的parallel指令可以用如下的形式表示:
#pragma omp parallel如果使用上面這條簡單的指令去運行並行計算的話，程序的線程數將由運行時系統決定（這里使用的算法十分復雜），典型的情況下，系統將在每一個核上運行一個線程。如果需要執行使用多少個線程來執行我們的並行程序，就得為parallel指令增加num_threads子句，這樣的話，就允許程序員指定執行后代碼塊的線程數。
#pragma omp parallel num_threads(thread_count)在程序中，為了能夠使用OpenMP函數，還需要在程序包含omp.h頭文件。
#include <omp.h>在這里，我們還是用之前介紹過的求積分的函數作為例子，
積分函數
//積分函數 double blog4::Trap(double a, double b, int n) { 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; return my_result; }
數學函數
//數學函數 double blog4::f(double x) { return pow(x, 3); }
OpenMP代碼
//第二個案例 void blog4::Test2(int argc, char* argv[]) { int thread_count = strtol(argv[1], NULL, 10); double global_result = 0.0; #pragma omp parallel num_threads(thread_count) { double my_result = 0.0; my_result = Trap(0, 3, 1024); #pragma omp critical global_result += my_result; } printf("%f\n", global_result); }
歸約
void blog4::Test2(int argc, char* argv[]) { int thread_count = strtol(argv[1], NULL, 10); double global_result = 0.0; #pragma omp parallel num_threads(thread_count) reduction(+: global_result) { global_result += Trap(0, 3, 1024); } printf("%f\n", global_result); }
運行：
后面的2指定使用2個線程運行程序。
通過以上的例子，大概知道OpenMP只要添加一條很簡單的parallel for指令，就能夠並行化大量的for循環所組成的串行程序。但使用它有幾天限制條件：
（1）、OpenMP只能並行化for循環，它不會並行while和do-while循環，而且只能並行循環次數在for循環外面就確定了的for循環。
（2）、循環變量只能是整型和指針類型（不能是浮點型）
（3）、循環語句只能是單入口單出口的。循環內部不能改變index，而且里面不能有goto、break、return。但是可以使用continue，因為它並不會減少循環次數。另外exit語句也是可以用的，因為它的能力太大，他一來，程序就結束了。
3.數據依賴問題
使用OpenMP在明面上有上面的這些限制規則，而需要關注的是一個更隱匿的問題：數據依賴問題。這其中涉及兩種依賴情況：數據依賴和循環依賴。為了弄清楚這兩種依賴的具體是什么，我們先來看一個例子。
計算前n個斐波那契（fibonacci）數：
然后加上OpenMP並行
通過測試，我們有時會得到：
1 1 2 3 5 8 13 21 34 55（正確）
然后也可能偶爾會得到：
1 1 2 3 5 8 0 0 0 0.（錯誤）
這究竟是發生了什么？似乎運行時，系統將fibo[2], fibo[3], fibo[4], 和 fibo[5]的計算分配給了一個線程，將fibo[6], fibo[7],fibo[8], 和 fibo[9]分配給了另外一個線程。所以在計算fibo[6]的時候，如果第一個線程在此線程之前已經完成，入口初始化為fibo[5]=8，得到的計算結果就是正確的；而如果在計算fibo[6]的時候，如果第一個線程在此線程之前還沒有完成，入口初始化為fibo[5]=0，則得到的計算結果就是錯誤的。
從這里，我們看到兩個要點：
（1）OpenMP編譯器不檢查被parallel for指令並行化的循環所包含的迭代間的依賴關系，而是由程序員來識別這些依賴關系。
（2）一個或更多個迭代結果依賴於其它迭代的循環，一般不能被OpenMP正確的並行化。
從求斐波那契（fibonacci）數的這個例子中，我們發現fibo[3], fibo[4]計算間的依賴關系稱之為數據依賴。而fibo[5]的值在一個迭代中計算，它的結果是在之后的迭代中使用，這樣的依賴關系稱之為循環依賴。
4.怎么尋找循環依賴及解決辦法
當我們試圖使用一個parallel for指令時，需要小心循環依賴關系，而不用擔心數據依賴。例如，在下列循環中：
int i = 0; for (i = 0; i < n; i++) { x[i] = a + i*h; y[i] = exp(x[i]); }
在這個程序中y[i]數據依賴於x[i]的，這種依賴使用OpenMP並行化的話，是沒有關系的。
在看一個例子，如下公式計算π：
使用能夠在串行代碼中實現的公式
//第5個案例 void blog4::Test5(int argc, char* argv[]) { double factor = 1; double sum = 0.0; int n = 1000; for (int i = 0; i < n; i++) { sum += factor / (2 * i + 1); factor = -factor; } double pi_approx = 4.0*sum; printf("%f", pi_approx); }
得到的結果3.140593，符合預期。我們在上面這個段代碼中加上OpenMP的並行化指令
//第5個案例 void blog4::Test5(int argc, char* argv[]) { double factor = 1; double sum = 0.0; int n = 1000; int thread_count = strtol(argv[1], NULL, 10);//omp_get_num_threads(); cout << thread_count << endl; #pragma omp parallel for num_threads(thread_count) reduction(+:sum) for (int i = 0; i < n; i++) { sum += factor / (2 * i + 1); factor = -factor; //cout << i << "-" << factor << endl; } double pi_approx = 4.0*sum; printf("%f", pi_approx); }
運行程序
得到結果2.730014，不符合預期。這是存在一個循環依賴，不能保證第k次迭代中對factor的更新對第k+1次迭代有影響，所以導致結果不正確。通過分析，對代碼進行相應的改動。
if (i % 2 == 0) factor = 1.0; else factor = -1.0; sum += factor / (2 * i + 1);
運行結果為：2.976587，依舊不正確，為什么呢？
缺省情況下，任何在循環前聲明的變量，在線程間都是共享的，。為了消除這種循環依賴關系之外，還需要保證每個線程有它自己的factor副本，OpenMP的語句為我們考慮了，通過添加private子句到parallel指令中來實現這一目標。
//第5個案例 void blog4::Test5(int argc, char* argv[]) { double factor = 1; double sum = 0.0; int n = 1000; int thread_count = strtol(argv[1], NULL, 10);//omp_get_num_threads(); cout << thread_count << endl; #pragma omp parallel for num_threads(thread_count) reduction(+:sum) private(factor) for (int i = 0; i < n; i++) { if (i % 2 == 0) factor = 1.0; else factor = -1.0; sum += factor / (2 * i + 1); //factor = -factor; //cout << i << "-" << factor << endl; //SumForNumber(); } double pi_approx = 4.0*sum; printf("%f", pi_approx); }
這下終於正確了。