【並行計算】基於OpenMP的並行編程

    咱們目前的計算機都是基於馮偌伊曼結構的，在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);
}

    這下終於正確了。