OpenMP併發程式設計快速入門
OpenMP是目前被廣泛接受的，用於共享記憶體並行系統的多處理器程式設計的一套指導性的編譯處理方案。它提供了對並行演算法的高層的抽象描述，程式設計師通過在原始碼中加入專用的pragma來指明自己的意圖，由此編譯器可以自動將程式進行並行化，並在必要之處加入同步互斥以及通訊。本文是OpenMP使用的一個初步介紹，期望能引領讀者進入併發程式設計的世界。
下面這個兩篇文章介紹了OpenMP環境配置方法的詳細步驟，對此不甚瞭解的朋友可以先參閱它們以做為必要之準備：
歡迎關注白馬負金羈的部落格 http://blog.csdn.net/baimafujinji，為保證公式、圖表得以正確顯示，強烈建議你從該地址上檢視原版博文。本部落格主要關注方向包括：數字影象處理、演算法設計與分析、資料結構、機器學習、資料探勘、統計分析方法、自然語言處理。
本文中的示例程式主要圍繞一個簡單的“利用積分法求
#include "stdio.h"
#include "omp.h"
static long num_steps = 100000000;
double step;
int main()
{
int i;
double 執行上述程式碼所得之結果如下：
pi = 3.141593, time = 0.307644 s
下面我們來寫第一個基於OpenMP的並行程式。其實OpenMP併發程式設計的本質就是自動生成若干個執行緒，以期充分利用現代多核處理器的平行計算能力。現在很多程式語言（例如Java）中也提供有多執行緒程式設計的介面或者函式，但是開發人員必須顯式地控制這些併發的執行緒。而OpenMP則相當於提供了一套自動的管理方案，大大簡化了開發的難度。
- 函式omp_set_num_threads()用來指定你準備建立的執行緒數；
- 函式omp_get_thread_num()用來獲取已經建立的執行緒編號；
- 最後把你準備併發實現的部分放進#pragma omp parallel{}標記的大括號內即可。
下面給出示例程式碼：
#include "stdio.h"
#include "omp.h"
#define NUM_THREADS 4
static long num_steps = 100000000;
double step;
int main()
{
int i, nthreads;
double x, pi, sum[NUM_THREADS];
step = 1.0/(double)num_steps;
omp_set_num_threads(NUM_THREADS);
double start = omp_get_wtime();
#pragma omp parallel
{
int i, id, nthrds;
double x;
id = omp_get_thread_num();
nthrds = omp_get_num_threads();
if(id == 0) nthreads = nthrds;
for(i=id, sum[id]=0.0; i<num_steps; i=i+nthrds){
x = (i+0.5)*step;
sum[id] += 4.0/(1.0+x*x);
}
}
for (i = 0, pi=0.0; i < nthreads; i++)
{
pi += step * sum[i];
}
double end = omp_get_wtime();
printf("pi = %f, number of threads = %d, time = %f s\n", pi, NUM_THREADS, end-start);
return 0;
}然後我們使用不同的執行緒數1~4來觀測程式的輸出如下：
pi = 3.141593, number of threads = 1, time = 0.444881 s
pi = 3.141593, number of threads = 2, time = 0.551260 s
pi = 3.141593, number of threads = 3, time = 0.564319 s
pi = 3.141593, number of threads = 4, time = 0.477291 s
這個結果值得我們討論的地方有：1）首先當執行緒數是1的時候，我們的併發程式就會退化成一個序列程式，但是由於我們增加了很多分發和收集的工作（有時我們說這是communication成本），所以它會比我們的序列程式其實更加耗時；2）提高執行緒的數量，似乎也沒有提高程式執行的速度。這是由於出現了“false sharing”現象，即如果獨立的資料元素碰巧位於相同的cache line上，那麼每次update都會導致cache lines與執行緒之間“slosh back and forth”（來來回回）。也就是說在序列時的記憶體訪問次數因為比並行時訪問的次數少，所以反而序列的程式執行得更快。
為了解決上述這個問題，我們就指定一下cache line的大小，這樣單獨一個執行緒每次訪問記憶體時，會一次性帶走它所能帶走的最大資料量，而且我們還會通過調整程式碼使得這些資料都會被訪問它的執行緒所使用。從而減少因為“slosh back and forth”而浪費的時間。來看下面這段示例程式碼：
#include "stdio.h"
#include "omp.h"
#define NUM_THREADS 4
#define PAD 8 //assume 64 byte L1 cache line size
static long num_steps = 100000000;
double step;
int main()
{
int i, nthreads;
double x, pi, sum[NUM_THREADS][PAD];
step = 1.0/(double)num_steps;
omp_set_num_threads(NUM_THREADS);
double start = omp_get_wtime();
#pragma omp parallel
{
int i, id, nthrds;
double x;
id = omp_get_thread_num();
nthrds = omp_get_num_threads();
if(id == 0) nthreads = nthrds;
for(i=id, sum[id][0]=0.0; i<num_steps; i=i+nthrds){
x = (i+0.5)*step;
sum[id][0] += 4.0/(1.0+x*x);
}
}
for (i = 0, pi=0.0; i < nthreads; i++)
{
pi += step * sum[i][0];
}
double end = omp_get_wtime();
printf("pi = %f, number of threads = %d, time = %f s\n", pi, NUM_THREADS, end-start);
return 0;
}同樣我們使用不同的執行緒數1~4來觀測程式的輸出如下：
pi = 3.141593, number of threads = 1, time = 0.446988 s
pi = 3.141593, number of threads = 2, time = 0.231959 s
pi = 3.141593, number of threads = 3, time = 0.224227 s
pi = 3.141593, number of threads = 4, time = 0.211466 s
顯然這一次我們的程式就執行得快了很多！但是要知道我們所使用的機器的cache line的大小到底是多少，這顯然有點為難開發人員了，我們能不能有一種更加優雅的方式來做上面那些事情呢？OpenMP確實也為我們提供了支援。
下面程式碼演示了OpenMP中提供的一種叫做 Mutual Exclusion的同步模式，在程式碼中由#pragma omp critical來指示：
#include "stdio.h"
#include "omp.h"
#define NUM_THREADS 2
static long num_steps = 100000000;
double step;
int main()
{
int nthreads;
double pi = 0.0;
step = 1.0/(double)num_steps;
omp_set_num_threads(NUM_THREADS);
double start = omp_get_wtime();
#pragma omp parallel
{
int i, id, nthrds;
double x, sum;
id = omp_get_thread_num();
nthrds = omp_get_num_threads();
if(id == 0) nthreads = nthrds;
for(i=id, sum=0.0; i<num_steps; i=i+nthreads){
x = (i+0.5)*step;
sum += 4.0/(1.0+x*x);
}
#pragma omp critical
pi += step * sum;
}
double end = omp_get_wtime();
printf("pi = %f, number of threads = %d, time = %f s\n", pi, NUM_THREADS, end-start);
return 0;
}同樣我們使用不同的執行緒數1~4來觀測程式的輸出如下：
pi = 3.141593, number of threads = 1, time = 0.367218 s
pi = 3.141593, number of threads = 2, time = 0.184200 s
pi = 3.141593, number of threads = 3, time = 0.183366 s
pi = 3.141593, number of threads = 4, time = 0.178658 s
恰當地運用OpenMP提供的各種同步模式（例如Atomic, Barrier, Critical等）能夠對程式執行效率的提升起到相當大的幫助。OpenMP中還有很多話題值得探討，有興趣的讀者可以參閱相關資料以瞭解更多。
等等，這樣我就算OpenMP入門了嗎？如果你翻看其他OpenMP的資料，下面我要介紹的這部分內容才往往是其他資料的開篇內容！這些方法和技巧會讓你頓時感覺OpenMP真的很容易上手。那前面的內容到底算什麼？按照Tim Mattson的話來說，前面這些其實是幫助你更好更深刻的理解OpenMP！
平行計算中最值得我們優化的地方就是迴圈。計算機最擅長做的事情就是重複大量地簡單計算，而這個“重複”計算就是通過迴圈來實現的。如果迴圈裡面的內容不存在相互依賴關係，也就是迴圈體中語句之間的順序是可以任意調整的，那麼你可以採用下面這種語法來讓OpenMP協助進行迴圈的並行優化：
#pragma omp parallel
{
#pragma omp for
for(i = 0; i < N; i++){
do_something;
}
}特別地，如果#pragma omp parallel{}的大括號中唯一的內容就是一個for迴圈，那麼上面的程式碼還可以採用下面這種簡單的寫法，它們的作用是完全等價的：
#pragma omp parallel for
for(i = 0; i < N; i++){
do_something;
}但是如果迴圈體中的內容之間有一定的依賴關係，我們該怎麼做呢？舉個例子：
int i, j, A[MAX];
j = 5;
for (i = 0; i < MAX; i++){
j+=2;
A[i]=big(j);
}這個程式碼中，迴圈體裡的兩條語句之間就具有一定的依賴性，每次big(j)操作都依賴於本輪迴圈最新得到之j值，如果你貿然使用並行方法，那麼整個迴圈中的順序將被徹底打亂，那程式最終所得之結果就難以保證了！
一個解決方案就是改寫程式碼消除依賴性，例如：
int i, A[MAX];
#pragma omp parallel for
for (i = 0; i < MAX; i++){
int j = 5 + 2*(i+1);
A[i]=big(j);
}上述程式碼的思路就是每輪迴圈都重新建立一個變數j，這個j的值由當前的i直接算得，那麼下面的big(j)就只能依賴於本輪的變數。因為沒有全域性變數，只有區域性變數，如果j還沒被建立，那麼big(j)也不會執行，這就是與之前程式最大不同的地方。
另外一個比較特殊，但也更為常見的情況是類似下面這種的：
double ave = 0.0, A[MAX];
int i;
for(i=0;i<MAX;i++){
ave+=A[i];
}
ave = ave/MAX;注意累加操作必須依賴於之前（也就是上一輪迴圈）所得之ave才能算得本輪之ave。這種“累加”操作非常常見，但由於存在依賴關係，我們之前的做法並不能奏效。類似的還有“累乘”等等。OpenMP把這類情況統稱為Reduction，並提供了很好的支援。
reduction(op:list)
A local copy of each list variables is made and initialized depending on the “op”。例如如果是累加那麼op就被初始為0，如果是累乘op就被初始為1。
Then updates occur on the local copy. Local copies are reduced into a single value and combined with the original global value.
下面這個例子演示了Reduction的使用
double ave = 0.0, A[MAX];
int i;
#pragma omp parallel for reduction (+:ave)
for(i = 0; i < MAX; i++){
ave += A[i];
}
ave = ave/MAX;最後我們改寫本文最初給出的那個序列程式，只作簡單改動：
#include "stdio.h"
#include "omp.h"
static long num_steps = 100000000;
double step;
int main()
{
int i;
double pi, sum = 0.0;
step = 1.0/(double)num_steps;
double start = omp_get_wtime();
#pragma omp parallel
{
double x;
#pragma omp for reduction(+:sum)
for(i=0; i<num_steps; i++){
x = (i+0.5)*step;
sum = sum + 4.0/(1.0+x*x);
}
}
pi = step * sum;
double end = omp_get_wtime();
printf("pi = %f, time = %f s\n", pi, end-start);
return 0;
}請讀者自行編譯並執行上述程式，正常情況下，這個程式應該會比使用critical的程式略慢一些，但速度提升仍然非常顯效。到此為止，你已經算是基本掌握了OpenMP並行程式設計的一些方法和思路了，對於有深入學習需求的朋友，可以參閱一下文獻【1】中後半程（14~27）視訊課程的學習。本文中的程式主要來自（1~13）視訊教學片中的示例。
參考文獻
[2] 左飛，程式碼揭祕——從C/C++的角度探祕計算機系統，電子工業出版社