コンテンツにスキップ

第05章 集合（後）

5.4　スキャン

スキャンは、入力シーケンスの部分リダクションを全て作成して、新しい出力シーケンスにする。包括的スキャンと排他的スキャンの2種類がある。

スキャンの図

・包括的スキャン：n番目の出力値は最初のnの入力値に対するリダクションである。

template<typename T, typename C>
void inclusion_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;
    }
}

包括的スキャンの例

・排他的スキャン：n番目の出力値が除外されたリダクション。

template<typename T, typename C>
void inclusion_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;
    }
}

排他的スキャンの例

２つの例には初期値がある。これらには2つの理由がある。

  • 排他的スキャンの最初の出力値を演算するときに単位元を持つ必要性を回避するため
  • タイルされた並列スキャンを記述するためにシリアルスキャンを有用なビルディング・ブロックとするため

スキャンは、ループ伝播依存があるが、並列化することが出来る。また、レデュースの並列化と似ているため、操作の順序を変更するための集約関数の結合性を活用することが出来るが、レデュースとは異なり、冗長なコスト演算コストがかかる。スパンがO(n)からO(logn)になる代わりに、多くのアルゴリズムではワークがほぼ倍になる。

スキャンを並列化する効率のいいアプローチは、フォーク・ジョイン・パターンに基づく手法である。

OpenMP

OpenMPには、組込み並列スキャンは用意されていない。8.11節に記載してあるフォーク・ジョインCilk Plusコードに似たツリーベースのスキャンは記載することが出来る。

３フェーズのスキャンの漸近的実行時間は、$T_P=O(N/P+P)$となる。P≪Nの場合、P項を無視して、N/P項のみ考慮すればよく、スピードアップは線形になる。Nが固定の場合、$N/P+P$の値は$P=\sqrt{N}$のときに最小化される。最大漸近的スピードアップは次のようになる。

T1TP=O(NN/P+P)=O(NN/N+N)=O(N) \frac{T_1}{T_P}=O(\frac{N}{N/P+P})=O(\frac{N}{N/\sqrt{N}+\sqrt{N}})=O(\sqrt{N})

これは$O(logN)$時間となるツリーベースのスキャン程漸近的に優れていないが、定数係数により魅力的になるかもしれない。しかし、他の並列スキャンのように、３フェーズのスキャンでは集約関数の呼び出しがシリアルスキャンの約２倍になる。

どのくらい時間が違うのか

・３フェーズのスキャン(よく利用されるらしい)

  • 最後を除いて、等しいサイズのタイルに入力を分割し、各タイルを並列にレデュースする。
  • リダクション値の排他的スキャンを実行する。並列スキャン全体が包括的であってもスキャンは常に排他的である。
  • 各タイルのスキャンを実行する。各タイルの初期値はフェーズ２の排他的スキャンの結果である。並列スキャンが包括的か排他的かによって、並列スキャンもそれに対応する。

図5.6 フェーズの図

OpenMP実装

template<typename T, typename R, typename C,typename S>
void openmp_scan(
    size_t n,　//インデックスの空間サイズ
    T initial,　//スキャンの初期値
    size_t tilesize,　//反復空間の各タイルのサイズ
    R reduce,　//ファンクター。reduce(i,size)は、[i,size)のインデックスに対するリダクション値を返す
    C combine,　//ファンクター。combine(x,y)はx⊕yを返す
    S scan　　　//ファンクター。scan(i,size,initial)は指定された初期値で開始する[i,i+size)のインデックスに対してスキャンを行う。
) {
    if(n > 0){
        //使用するタイルの数にtをセット。最大で1スレッドあたり1タイルで
        //リクエストされたtilesizeより小さいタイルはない
        size_t t = std::min(size_t(omp_get_max_threads()), (n - 1)/tilesize + 1);
        //各タイルのリダクション値を保持する空間を割り当てる
        temp_space<T> r(t);
        //1タイルあたり1つのスレッドをリクエスト
    }
#pragma omp parallel num_threads(t)
    {
        //実際に許可されたスレッド数を確認(要確認！！)
        size_t p = omp_get_num_threads();
        //スレッドごとに1つのタイルになるようにtilesizeを再計算
        tilesize = (n+p-1)/p;
        //最後のタイルのインデックスを示す。
        size_t m = p - 1;
#pragma omp for
        //ステップ1のリダクション
        for(size_t i = 0;i <= m; ++i){
            r[i] = reduce(i*tilesize, i==m ? n-m*tilesize : tilesize);
        }
#pragma omp single
        //ステップ2の排他的スキャン
        for(size_t i = 0;i <= m; ++i){
            T tmp = r[i];
            r[i] = initial;
            initial = combine(initial,tmp);
        }
#pragma omp for
        //ステップ3の包括的スキャン
        for(size_t i = 0;i <= m; ++i){
            scan(i*tilesize, i==m ? n-m*tilesize : tilesize, r[i]);
        }
    }
}

このコードは利用可能なスレッド数を確認し、1スレッドあたり１つのタイルを扱うようにしている。(要確認！！)の場所で、確認しなかった場合、コードは正しいままだが、いくつかのスレッドが他のスレッドよりも多くタイルを実行する場合があり、その場合ロード・インバランスが発生しパフォーマンスに影響を及ぼす。

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

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

図5.7 マップとスキャンの融合

左の図は、マップ→フェーズ(1.2.3)→マップの順で行われている。最初のマップのタイルはスキャンの最初のフェーズのシリアル・リダクションと組み合わせることが出来る。また、３番目のフェーズのスキャンされたタイルは、次のタイルされたマップと組み合わせることが出来る。融合した結果が右の図。 　この結果、算術的なコードブロックが作成され、メモリー・トラフィックおよび同期のオーバーヘッドが大幅に削減される。

5.6 積分

関数の累積として知られる集計関数のスキャンには、いくつかの応用性がある。演算を行えば、一定時間での任意の区間に対する近似の積分を行える。

5.6.1 問題の説明

関数fと区間[a,b]に対し、[a,b]の任意の部分区間でfの定積分を迅速に計算できるテーブルをあらかじめ計算する。$\Delta x=(b-a)/(n-1)$であるとし、テーブルはfのサンプルの累積和で、スケールは$\Delta x$です。

tablei=Δx∑0if(a+iΔx) table_i = \Delta x \sum_{0}^{i} f(a+i\Delta x)

区間[c,d]のfの積分は次のように推定できる。

∫cdf(x)dx≈interp(d)−interp(c) \int_{c}^{d} f(x)dx \approx interp(d) - interp(c)

interp(x)は、テーブルの線形補間を示す。

線形補間：既知のデータ点を使って未知のデータ点を推定する方法。既知の２つの点を直線で結び、その直線上の任意の点を求めることにより補間を行う。

どんな風に積分するのか

与えられた区間[a,b]がある。

その区間をn等分する。そうしたら、aから一個目の青い線の面積を求める。次にaから2個目・・・のように求めていく。これがスキャン操作にあたる。

与えられた区間[c,d]を求めたいので、そのc,dを線形補間する。

与えられた区間[c,d]を上記のように、[a,d]から[a,c]の面積を引く形にすると求めたい[c,d]が求まる。

5.6.2 シリアル実装

次のコード5.17は、サンプリングおよび和計算のシリアル実装を示している。ループ反復がそれぞれの前の反復に依存するため、ループ伝播依存がある。

リスト5.17

```C++　 template void serial_prepare_integral_table( X a,　//サンプリングの開始位置 X b,　//サンプリングの終了位置 size_t n,　//処理するサンプルの数 Y table[],　//テーブルサンプルのデスティネーション F f　//関数パラメーター ) { // 空のリクエストを制御、0の場合何もせず終了 if(n==0) return; 

// Δxの計算。サンプルの間隔を計算する。aからbまでの間をn-1区間に分割する。
const X dx = (b-a)/X(n-1);

//積分スキャン
Y sum = Y(0);
for(size_t i = 0; i < n; ++i){
    sum += f(a+dx*i);
    table[i] = sum * dx;
}

} リスト5.20は２つのサンプルから定積分を求める方法を示している。配列をリニアに補間してルックアップするヘイパー関数serial_sampleを定義する。 リスト5.20C++ tempalte Y serial_sample{ size_t n, // テーブルの数 Y table[],　// テーブル(面積が入る) X x　//　サンプルの位置 } { //サンプル位置の整数部を計算 X i = floor(x); n=3,i=1 //y0とy1はtable配列の値が入る。 Y y0 = i < X(0) ? Y(0) : table[i < X(n) ? size_t(i) : n - 1]; y0=table[1] Y y1 = i+1 < X(0) ? Y(0) : table[i+1 < X(n) ? size_t(i+1) : n - 1]; y1=table[2] // サンプル間をリニアに補間 return y0+(y1-y0)(x-i); } tempalte Y serial_integrate( size_t n,　//テーブルのサンプルの数 Y table[], // サンプルの累積 X a, // 区間 a X b, // 区間 b X x0, //積分の下限 x0 X x1　//積分の上限 x1 ) { // 1/Δxの計算。 X scale = X(n-1)/(b-a); // x0での補間値の計算。interp(c)の計算 Y y0 = serial_sample(n, table, scale(x0-a)); // x1での補間値の計算。interp(d)の計算 Y y1 = seiarl_sample(n, table, scale*(x1-a)); //積分計算 return y1-y0; } ```

視覚的に表してみる。

5.6.4 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にマップする関数
) {
    if (n == 0) return;

    const X dx = (b - a) / (n - 1);

    openmp_scan(
        n, Y(0), //インデックスの空間サイズ,スキャンの初期値
        1024,    //反復空間の各タイルのサイズ
        [=,&table](size_t i, size_t m) -> Y {
            Y sum = Y(0);
            for (; m > 0; --m, ++i) {
                sum += (table[i] = f(a + dx * i));
            }
            return sum;
        }, //指定された範囲のサンプリング点の部分和の計算
        std::plus<Y>(), //各タイルで計算された部分和を結合し、累積和にする
        [=,&table](size_t i, size_t m, Y initial) {
            for (; m > 0; --m, ++i) {
                initial += table[i];
                table[i] = initial * dx;
            }
        }
    ); //累積和を計算し、それにスケールΔxをかけて積分値を計算
}

OpenMPのこのコードは、上のserial_prepare_integral_tableに対応している。実際の[c,d]の計算は、また別の関数などになる。

・このコードの配列表記はループ構造に置換することができる。

OpenMPの配列表記には、スレッド並列とベクトル並列を活用する効果的な方法であるため、上の変更はオプションである。もちろん、この変更を行うには両方の並列化モデルをサポートするコンパイラーが必要になる。

スキャンはOpenMPの組込みではないが、5.4節の３フェーズのアプローチを使用する。

5.7　まとめ

この章では、レデュースとスキャンの２つの集合パターンの概要と実装におけるオプションなどを説明した。より詳細なものは後の章でやる。

レデュースとスキャンはマップとともによく使用され、隣接するマップと実装(の一部)を融合、あるいはマップとレデュース/スキャンの順序を逆にして、最適化することが出来る。