5.4スキャン

スキャン(scan)収集操作は、入力シーケンスの部分リダクションをすべて作成して、新しい出力シーケンスに操作。スキャンには以下の２つのバリエーションがある。

• 包括的スキャン: n番目の出力値は最初のn個の入力値に対するリダクションである。
• 排他的スキャン: n番目の出力値は最初のn-1個の入力値に対するリダクションである。

包括的スキャンのシリアル実装（リスト5.14）

template<typename T, typename C>
void inclusive_scan(
    size_t n,          // 要素の数 [# シーケンスのサイズを指定]
    const T a[],       // 入力コレクション [# 処理対象となる配列]
    T A[],             // 出力コレクション [# 結果を格納する配列]
    C combine,         // 集約ファンクター [# 値を結合する関数を指定]
    T initial          // 初期値 [# 計算のスタートとなる値]
) {
    for (size_t i = 0; i < n; ++i) {
        initial = combine(initial, a[i]); // [# 現在の初期値と入力値を結合]
        A[i] = initial;                  // [# 結果を出力配列に格納]
    }
}

- 実行結果

排他的スキャンのシリアル実装（リスト5.15）

template<typename T, typename C>
void exclusive_scan(
    size_t n,          // 要素の数 [# シーケンスのサイズを指定]
    const T a[],       // 入力コレクション [# 処理対象となる配列]
    T A[],             // 出力コレクション [# 結果を格納する配列]
    C combine,         // 集約ファンクター [# 値を結合する関数を指定]
    T initial          // 初期値 [# 計算のスタートとなる値]
) {
    if (n > 0) {
        for (size_t i = 0; i < n - 1; ++i) {
            A[i] = initial;              // [# 現在の初期値を出力配列に格納]
            initial = combine(initial, a[i]); // [# 初期値を更新]
        }
        A[n - 1] = initial;             // [# 最後の要素を設定]
    }
}

• 実行結果

OpenMP

OpenMPには組み込みの並列スキャンは用意されていないが、OpenMP 3.0でサポートされたタスクを使用して、フォーク・ジョインCilk Plusコードに似たツリーベースのスキャンを実装できる(8.11節に記載)。以下のような3フェーズのアプローチが一般的である。

3フェーズのスキャンの実行概要

• フェーズ1: 入力を等しいサイズのタイルに分割し、各タイルを並列にリダクションする。
• フェーズ2: リダクション値の排他的スキャンを実行する。このスキャンは常に排他的である。
• フェーズ3: 各タイルのスキャンを実行する。初期値はフェーズ2の排他的スキャンの結果から取得する。

3フェーズのスキャンの漸近的実行時間は $T_p = \Theta(P + \frac{N}{P})$ $P \ll N$: スレッドの数 $P$ がデータの総量 $N$ に比べて非常に少ない場合、$P$ の影響を無視でき、実行時間は $T_p \approx \frac{N}{P}$ と見なせる。 スピードアップは、シリアル実行時間 $T_1$ に対する並列実行時間 $T_p$ の比として定義される。 固定 $N$ の場合:$P + \frac{N}{P}$の値は$P = \sqrt{N}$のとき最小化される。その結果以下のようにスピードアップを計算できる $$\frac{T_1}{T_p} = \Theta\left(\frac{N}{\frac{N}{\sqrt{N}} + \sqrt{N}}\right) = \Theta(\sqrt{N})$$

OpenMPによる3フェーズのスキャンのタイル実装

template<typename T, typename R, typename C, typename S>
void openmp_scan(
    size_t n,          // 要素の数 [# シーケンスのサイズを指定]
    T initial,         // 初期値 [# 計算のスタートとなる値]
    size_t tilesize,   // タイルサイズ [# 分割単位のサイズを指定]
    R reduce,          // リダクション関数 [# 範囲をリダクション]
    C combine,         // 結合関数 [# 値を結合]
    S scan             // スキャン関数 [# 個別スキャン操作]
) {
    if (n > 0) {
        size_t t = std::min(size_t(omp_get_max_threads()), (n - 1) / tilesize + 1);
        temp_space<T> r(t); // [# 各タイルのリダクション値を格納]

        #pragma omp parallel num_threads(t)
        {
            size_t p = omp_get_thread_num();
            tilesize = (n + p - 1) / p;

            #pragma omp for
            for (size_t i = 0; i < t; ++i) {
                r[i] = reduce(i * tilesize, tilesize); // [# タイルごとにリダクション]
            }

            #pragma omp single
            {
                for (size_t i = 0; i < t; ++i) {
                    T temp = r[i];
                    r[i] = initial;                   // [# スキャン結果を初期化]
                    initial = combine(initial, temp); // [# スキャン結果を更新]
                }
            }

            #pragma omp for
            for (size_t i = 0; i < t; ++i) {
                scan(i * tilesize, tilesize, r[i]);   // [# タイルごとにスキャン]
            }
        }
    }
}

　多くのOpenMPと同様に、このコードは利用可能なスレッド数を確認し、スレッドあたり１つのタイルを扱うようにする。実際のスレッド数に基づいてタイルサイズと数を再計算する。もし再計算しないと、いくつかのスレッドがほかのスレッドよりも多くのタイルを実行すると、 ロード・インバランス（負荷不均衡） が発生しパフォーマンスに影響を及ぼす。

• 長所：メモリトラフィックが最小限でスレッドの再利用が可能。
• 短所：並列領域内でスレッドが待機する場合があり、効率が低下する可能性がある。

5.5 マップとスキャンの融合

レデュースのように、スキャンも隣接する操作と融合することで最適化できる。

最初のマップタイルは、スキャンの最初のフェーズのシリアル・リダクションと組み合わせ可能。3番目のフェーズのスキャンされたタイルは、次のタイルされたマップと組み合わせることが可能。 これにより 算術的なコードブロックが作成され、メモリー・トラフィックおよび同期のオーバーヘッドが大幅に軽減される。

5.6 積分

スキャンは並列化が困難なアルゴリズムを並列化する要因となる場合がある。

問題の説明

関数 $f$ と区間 $[a, b]$ が与えられたとき、区間内の任意の部分区間で迅速に定積分を計算するテーブルを構築する。

テーブル構築式: $$\text{table}_i = \Delta x \cdot \sum_{j=0}^i f(a + j \cdot \Delta x)$$ ここで、$\Delta x = \frac{b-a}{n-1}$。

区間 $[c, d]$ の積分推定: $$\int_c^d f(x) dx \approx \text{interp}(d) - \text{interp}(c)$$ $\text{interp}(x)$ は線形補間関数。

$$f(x) = x^3 - 3x + 1$$

---

OpenMP

template<typename X, typename Y, typename F>
void openmp_prepare_integral_table(
    X a, // サンプリングの開始位置
    X b, // サンプリングの終了位置
    size_t n, // サンプルの数
    Y table[], // テーブル内サンプルの累積ディスティネーション
    F f // X -> Y にマップする関数
) {
    // サンプル数が0の場合は計算を行わずに終了
    if (n == 0) return;

    // サンプル間隔を計算 (全体の範囲をサンプル数で分割)
    const X dx = (b - a) / (n - 1);

    // 部分和を格納するベクターを初期化
    std::vector<Y> partial_sums(n, Y(0));

    // 各サンプルポイントで関数fを評価して部分和を計算
    #pragma omp parallel for
    for (size_t i = 0; i < n; ++i) {
        partial_sums[i] = f(a + dx * i); // 関数fの評価結果を格納
    }

    // 累積和の計算 (スキャン操作)
    for (size_t step = 1; step < n; step *= 2) { // ステップサイズを倍にして繰り返す
        #pragma omp parallel for
        for (size_t i = step; i < n; ++i) {
            // ステップサイズ分だけ前の値を加算
            partial_sums[i] += partial_sums[i - step];
        }
    }

    // スケーリングして結果をテーブルに格納
    #pragma omp parallel for
    for (size_t i = 0; i < n; ++i) {
        table[i] = partial_sums[i] * dx; // スケールされた積分値を格納
    }
}

5.7 まとめ

この章ではレデュースとスキャンの2つの集合パターンの概要と実装におけるさまざまなオプションの説明。 一般的にTBBやCilk Plusを使用すると、この章のレデュースの実装について考慮する必要はない。OpneMPによる3フェーズのスキャンはこのあとのフォーク・ジョインの実装よりはスケーラブルではない。