OpenMP/Tasks

如果你跟著上一章並做了練習，你就會發現我們開發的求和函式很不穩定。（公平地說，我們設定了資料來展示這種不穩定性。）實際上，簡單地求和的 $O(n)$ 舍入誤差：誤差與我們要累加的元素數量成正比。我們可以透過更改求和演算法來解決這個問題。有三個候選的穩定求和演算法

將數字從小到大排序，然後像以前一樣累加它們。問題：這將問題的複雜度從線性變為 $O(n\log n)$ .
使用Kahan 演算法。雖然該演算法需要線性時間並且非常穩定，但在實踐中它相當緩慢，並且比我們的第三種選擇更難並行化，第三種選擇是
分治遞迴。這種演算法的舍入誤差是 $O(\log n)$ ，即與元素數量的對數成正比。

基本的分治求和在 C 中很容易表達

float sum(const float *a, size_t n)
{
    // base cases
    if (n == 0) {
        return 0;
    }
    else if (n == 1) {
        return *a;
    }

    // recursive case
    size_t half = n / 2;
    return sum(a, half) + sum(a + half, n - half);
}

如果你在上一章為該程式開發的程式中使用這個 sum 的定義，你會發現它產生了完全預期的結果。但是這個演算法沒有迴圈，那麼我們如何使用 OpenMP 製作並行版本呢？

我們將使用 OpenMP 中的任務構造，將問題視為任務並行而不是資料並行。在第一個版本中，sum 的任務遞迴版本如下所示

float sum(const float *a, size_t n)
{
    // base cases
    if (n == 0) {
        return 0;
    }
    else if (n == 1) {
        return 1;
    }

    // recursive case
    size_t half = n / 2;
    float x, y;

    #pragma omp parallel
    #pragma omp single nowait
    {
        #pragma omp task shared(x)
        x = sum(a, half);
        #pragma omp task shared(y)
        y = sum(a + half, n - half);
        #pragma omp taskwait
        x += y;
    }
    return x;
}

我們引入了兩個任務，每個任務都設定一個被宣告為 shared 的變數，該變數與另一個任務共享。如果我們沒有宣告共享變數，每個任務都會設定自己的區域性變數，然後丟棄結果。然後我們使用 #pragma omp taskwait 等待任務完成，並結合遞迴結果。

你可能會對緊跟著 #pragma omp single nowait 的 #pragma omp parallel 感到驚訝。問題是第一個pragma會導致池中的所有執行緒執行下一段程式碼。single 指令會導致除了一個執行緒（通常是第一個遇到該程式碼塊的執行緒）以外的所有執行緒不執行它，而 nowait 會關閉 single 上的屏障；在封閉的 parallel 區域上已經有一個屏障，其他執行緒會蜂擁而至。

不幸的是，如果你真的嘗試執行這段程式碼，你會發現它仍然不是特別快。原因是任務太細粒度了：在遞迴樹的底部附近， ${\tfrac {n}{2}}$ 次呼叫將兩個元素的陣列分成每個處理一個元素的子任務。我們可以透過引入除了基本情況和遞迴情況之外的“中間情況”來解決這個問題，該情況是遞迴的，但不涉及設定並行任務：如果遞迴達到預先指定的截止值，它將不再嘗試為 OpenMP 執行緒池設定任務，而是隻進行遞迴求和。

練習：在遞迴中引入額外的情況，並衡量程式的速度。不要提前檢視下一個程式，因為它包含了這個練習的答案。

現在我們實際上有兩個遞迴合二為一：一個帶並行遞迴情況，另一個是序列遞迴。我們可以透過在每個級別進行更少的檢查來解開這兩個遞迴，以獲得更好的效能。我們還將並行設定程式碼分離到驅動函式中。

#include <stddef.h>

#define CUTOFF 100  // arbitrary

static float parallel_sum(const float *, size_t);
static float serial_sum(const float *, size_t);

float sum(const float *a, size_t n)
{
    float r;

    #pragma omp parallel
    #pragma omp single nowait
    r = parallel_sum(a, n);
    return r;
}

static float parallel_sum(const float *a, size_t n)
{
    // base case
    if (n <= CUTOFF) {
        return serial_sum(a, n);
    }

    // recursive case
    float x, y;
    size_t half = n / 2;

    #pragma omp task shared(x)
    x = parallel_sum(a, half);
    #pragma omp task shared(y)
    y = parallel_sum(a + half, n - half);
    #pragma omp taskwait
    x += y;

    return x;
}

static float serial_sum(const float *a, size_t n)
{
    // base cases
    if (n == 0) {
        return 0.;
    }
    else if (n == 1) {
        return a[0];
    }

    // recursive case
    size_t half = n / 2;
    return serial_sum(a, half) + serial_sum(a + half, n - half);
}

當並行任務內的程式碼花費更多時間進行計算，而花費更少時間進行記憶體訪問時，這種技術效果更好，因為這些記憶體訪問可能需要在處理器之間進行同步。