インテル® Cilk™ Plus によるデータ並列処理とスレッド並列処理: ビジュアル・コンピューティングのケーススタディー

この記事は、インテル® デベロッパー・ゾーンに掲載されている「Using Intel® Cilk™ Plus to Achieve Data and Thread Parallelism: A Case Study for Visual Computing」の日本語参考訳です。

---

概要

高度なベクトル命令に対応したマルチコアシステムは、サーバー、ワークステーション、デスクトップからラップトップ、モバイルデバイスに至るまで、さまざまな分野に広がりを見せています。データ並列処理とスレッド並列処理の両方を利用することで、ユーザー・アプリケーションのパフォーマンスが大幅に向上する可能性があります。C/C++ 向けの言語拡張であるインテル® Cilk™ Plus では、データ並列とスレッド並列を簡単に表現することができます。簡潔で保守しやすく、将来のハードウェアでもスケーラブルなソースコードを維持しつつ、アルゴリズムを自然な方法で表現でき、信頼性が高く予測可能な並列処理を提供します。さらに、既存の C/C++ データ構造やコードともシームレスに動作します。

この記事では、ビジュアル・コンピューティング・ベンチマーク (AOBench – アンビエント・オクルージョン・ベンチマーク) を例に、インテル® Cilk™ Plus でデータ並列処理およびスレッド並列処理を表現し、シーケンシャルな C プログラムのパフォーマンスを大幅に向上させる方法を説明します。

1. はじめに

インテルでは、並列プログラミング向けの新しいプログラミング・モデルを提供しています。同時に、OpenMP* や MPI のような標準規格への取り組みも継続して行っています。インテル® Cilk™ Plus は、データ並列構造とタスク並列構造の両方を提供する言語拡張です。インテル® C/C++ コンパイラーを含む、インテル® Parallel Composer XE 2011 およびインテル® C++ Composer XE 2011^[3] 以降に含まれています。この記事では、最初にインテル® Cilk™ Plus で利用可能な並列言語拡張の概要、次にアンビエント・オクルージョン・ベンチマーク (AOBench) アプリケーションについて、最後にインテル® Cilk™ Plus への移植とその結果に関して説明します。

2. インテル® Cilk™ Plus

インテル® Cilk™ Plus は、ソースコードでタスク/データレベルの並列性を明示的に表現する言語拡張です。既存のコードに新しいキーワードを追加するだけでタスクを並列化できます。データ並列処理またはベクトル化に関しては、インテル® Cilk™ Plus の配列表記拡張を使ってアルゴリズムを記述することで、保守しやすく、ベクトル化が予測可能になります。さらに、スカラー関数を要素関数またはベクトル関数として宣言することもできます。これにより、ベクトル・コンテキストで呼び出された場合、コンパイラーはその関数のベクトル・バージョンを生成します。

2.1 タスク並列処理

インテル® Cilk™ Plus には、2 つのタスクレベルの並列化構造があります。1 つ目の構造は、_Cilk_spawn および _Cilk_sync の 2 つのキーワードを提供します。図 1 は、C で記述された有名なフィボナッチ・コードの再帰実装 (シリアルバージョン) と並列実装を比較したものです。

図 1. 左はフィボナッチのシリアル実装で、右はインテル® Cilk™ Plus キーワードを利用した (同一コードの) 並列実装。オリジナルコードにインテル® Cilk™ Plus キーワードを追加しただけであることが分かります。

シリアルプログラムから並列プログラムへの変更は、必要に応じて簡単なキーワードを挿入するだけで済みます。プログラム自体の変更は必要ありません。また、キーワードを挿入しても、プログラムの簡潔さや保守性には影響しません。オリジナルのシリアルコードと同じプログラムロジックの分かりやすい並列コードになります。インテル® Cilk™ Plus フレームワークは、再帰により発生する多くのタスクのロードバランスを動的に行います。これらは、既存のシリアル C/C++ コードを並列化するという、インテル® Cilk™ Plus の主要な設計方針の 1 つをよく表しています。インテル® Cilk™ Plus の利点は、シリアル等価性 (決定論性) を保証しつつ、最小限の変更でコードを並列化できることです。

2.1.1 _Cilk_spawn キーワード

_Cilk_spawn キーワードは、図 1 の fib(n-1) のような関数呼び出しに直接利用します。呼び出し先と呼び出し元は同時に実行することができます。この並行性を除けば、通常の関数呼び出しと同じです。_Cilk_sync キーワードは、呼び出し関数 (親関数) に対しすべての子タスクが完了し、コントロールが返されるまで待機することを指示します。インテル® Cilk™ Plus の並列領域の概念は新しいものではなく、親から見た並列領域として、そしてスポーン可能な単位として、関数の概念を再利用しています。スポーンされる関数に渡す引数は、親タスクが子を切り離す前に評価されなければいけません。これは、C/C++ 言語において呼び出し前に引数が評価されなければいけないのと似ています。図 2 の例で、関数呼び出し g() と h() は f() の呼び出しがスポーンされる前に評価されます。そのため、g() と h() の呼び出しは、f() の呼び出しの前にシリアルに実行されます。

2.1.2 _Cilk_sync キーワード

図 1 の fib() 関数にある _Cilk_sync キーワードは、変数 x と y の値が参照される前に fib(n-1) と fib(n-2) の非同期実行により更新されることを保証するのに必要です。明示的な _Cilk_sync キーワードに加えて、すべてのスポーンを行う親関数の最後には暗黙的な _Cilk_sync があります。この暗黙的な _Cilk_sync により、スポーンを行う関数はそのすべての子タスクの完了を待機してからリターンします。そうすることで、呼び出された関数がリターンした後にダングリング・タスク (親なしの子タスク) が残りません。

2.1.3 _Cilk_for キーワード

2 番目の構造は _Cilk_for キーワードを提供します。このキーワードを指定して、C/C++ の for ループを並列化できます。_Cilk_for では、for ループを実行する前にその反復回数が判明している必要があります。図 3 の構文について考えてみます。

図 3: _Cilk_for キーワード

_Cilk_for 構造には次のような制限があります。

コンパイラーはコンパイル時に違反を検出し、これらの制限違反を通知します。ただし、比較値とインデックス変数がポインターにより更新される場合は除きます。

2.2 インテル® Cilk™ Plus の配列表記 (アレイ・ノーテーション) を利用するデータ並列処理

インテル® Cilk™ Plus で提供される新しい構文を使って、配列を簡単に操作することができます。配列表記により、プログラムで直接データ並列配列操作 (ベクトル) を表現できます。ほとんどの C/C++ 実装は、通常の C/C++ 配列 (例: float a[N];) に対するデータ並列操作を表現する手段を提供していません。代わりに、プログラマーはループを使って、配列の要素に対する操作として表現しなければいけないため、操作が連続したシリアル順序になります。シリアル順序のループには欠点があります。シリアルループをベクトルループに変換するため、コンパイラーはループのベクトル操作がスカラー操作と等価であり、言語の規格に準拠していることを証明しなければいけません。例えば、配列の参照にポインターを利用している場合、追加の属性やキーワードなしでは、コンパイラーは 2 つのポインターがメモリーのオーバーラップしない領域を指していることを証明できない可能性があります。また、保守を行う場合も、反復順序が重要でないことが分かりづらいでしょう。

配列表記 a[:] = b[:] + c[:] は、配列 b と c の要素を任意の順序で加算できることをコンパイラーに知らせます。配列表記のセマンティクスでは、右辺と左辺が部分的にオーバーラップする場合、動作は不定です。これらのセマンティクスは、詳細な解析や restrict のようなプラグマ、属性、キーワードがなくても、コンパイラーがベクトルコードを生成できるように支援します。図 4 は、配列表記で表現された倍精度実数 Alpha X Plus Y (daxpy) アルゴリズムの例です。x[:] と y[:] は倍精度の配列で、a は倍精度のスカラーです^[1]。y[:] と完全にオーバーラップしているため、+= 演算子を使えることが分かります。配列表記のスカラーオペランドは、部分配列の各要素に提供されます。

図 4. 配列表記で表現された Daxpy/saxpy アルゴリズム

2.2.1 部分配列演算子

インテル® Cilk™ Plus 言語拡張は、次の構文を持つ部分配列演算子を定義します: array_base[begin:length:stride]
各要素の説明は次のとおりです。

• array_base は、配列、ポインター、C99 の可変長配列を含む C/C++ で許可されている配列表現です。
• begin は、部分配列の最初の要素のインデックスです。
• length は、部分配列中の要素の数です。
• stride は、要素間のインクリメントです。stride はオプションで指定できます。省略した場合、デフォルト値 1 が使用されます。

セクション演算子の省略形 [:] は、部分配列中のすべての要素を意味します。

セクション演算子は、例えば a[0:5][3][0:8:2] のように、多次元配列に適用したり、配列インデックスと組み合わせることもできます。この例で、a は 3 次元配列または配列へのポインターです。ランクは、セクション演算子が使われる次元数であって、添字ではありません。この例では、ランクは 2 です。つまり、このセクションは 3 次元配列の 2 次元のスライスです。

2.2.2 組込み関数

配列の要素の和のように、よく使われる操作の組込み関数^[2] が提供されています。例えば、ドット積演算は次のように表現できます。

図 5.   _sec_reduce_add で表現されたドット積演算

配列表記は、配列操作レベルでアルゴリズムを表現できる自然な方法です。配列表記では、ソースコードがアーキテクチャーに依存しないため、簡潔で保守しやすく、将来のハードウェアでもスケーラブルなソースコードを記述できます。

2.2.3 インテル® Cilk™ Plus の要素関数

プログラマーは C/C++ の標準表記でスカラー関数を記述し、要素関数として宣言することができます。要素関数は、通常、配列の複数の要素を任意の順序で処理する際に使用されます。図 6 は、要素関数を用いて配列の要素単位の加算を実行する単純な例です。

図 6.   要素関数で配列の要素単位の加算を実行する例

__declspec(vector) を追加すると、コンパイラーは関数「v_add」を要素関数として使用します。そして、スカラー・コンテキストで呼び出されるバージョンとベクトル・コンテキストで呼び出されるバージョン (つまり、ベクトル引数を使って、関数の本体でスカラーではなくベクトルで処理する) を生成します。そのため、プログラマーがアセンブリー・コードで組込み関数の呼び出しを記述したり、アセンブリー命令でベクトル・バージョンを記述しなくてもパフォーマンスが向上する可能性があります。

「for」ループで関数を呼び出す代わりに、プログラマーは「_Cilk_for」を指定することができます。その場合、コンパイラーはベクトル・バージョンを呼び出し、スレッドを活用して複数のベクトルを並列に実行します。この言語構造を利用することで、マルチコアとベクトルレベルの並列処理の両方の利点が得られます。

3. ベンチマークについて

アンビエント・オクルージョン・ベンチマーク (AOBench) は、有名なビジュアル・コンピューティング・ベンチマークで、アンビエント・オクルージョンのレンダリングに用いられます。アンビエント・オクルージョンとは、非反射面を作成することで、3D グラフィックの局部的反射モデルにリアルさを追加するシェーディング手法です。AOBench は、当初は Processing (http://processing.org/ (英語)) という言語でコーディングされていました。その後、ハイエンド CPU から GPU、スマートフォンまで、多種多様なプラットフォームをサポートするべく、C、Java*、Objective-C、CUDA、OpenCL*、Go などのさまざまな言語に移植されました (http://code.google.com/p/aobench (英語))。

AOBench は、複数の言語とプラットフォームにおける浮動小数点パフォーマンスのベンチマークを提供しています。僅か 400 行のコードで構成されており、標準の算術関数のみを必要とするため、その分かりやすさから、さまざまな言語およびプラットフォームへの移植が促進されました。AOBench をインテル® Cilk™ Plus へ移植することで、インテル® Cilk™ Plus とその他の実装によるパフォーマンスを比較できます。AOBench は、3 つの球体と 1 つの平面からなる決められたテストシーンにレイ・トレーシング・カーネルを適用します。各シーンの光線とオブジェクトの交点ごとに、シーンへランダムに光線を投射してアンビエント・オクルージョンの近似値を計算します。

ここでは、次の手法を使って AOBench の計算リソースの使用率を向上させます。

1. 構造体配列データ構造
2. ネストされた計算の反復
3. データ依存性を持つ分岐
4. C のショート関数
5. ビルトインの超越関数 (sin、cos、sqrt、rand)

AOBench のデータ構造は struct データ型を使用しています。3 つの球体と 1 つの平面からなるシーンは、この struct データ型で定義されています。出力として、AOBench は、球体と平面が交差するグレースケールのイメージを生成します。

以下は AOBench の擬似コードです。

上記の関数はすべて、ドット積、外積、正規化のような小さなベクトル関数を呼び出します。また、平方根、余弦、正弦のような超越関数も利用します。最も多くの計算が行われるのは compute_ambient_occlusion() 関数です。この関数には、乱数ジェネレーターを呼び出し、ランダムな光散乱を作成する小さな入れ子のループがあります。作成される光散乱ごとに、シーン・オブジェクトとの交点を計算する関数が呼び出されます。合計で 12 個の関数があり、実行時に入れ子の深さは最大で 6 段になります。

4. インテル® Cilk™ Plus による変更

このセクションでは、AOBench のシリアル・アルゴリズムを、データレベルの並列性 (DLP) とタスクレベルの並列性 (TLP) の両方を利用する並列実装に変更する基本ステップを説明します。最初にインテル® Cilk™ Plus の配列表記を利用して、シリアル・アルゴリズム内の DLP をシングルコアの VPU (ベクトル・プロセシング・ユニット、例: インテル® ストリーミング SIMD 拡張命令 (インテル® SSE) など) にマップし、次にインテル® Cilk™ Plus のタスクメカニズムを利用して、アルゴリズムの高度な TLP を複数のコアにマップする方法を示します。

VPU を利用するには、まずアルゴリズム中の SIMD (Single-Instruction, Multiple-Data) 演算を見つけます。多くの計算集約型アプリケーションと同様に、AOBench は浮動小数点型のデータセットの走査にループを使っています。各反復で実行される演算が同じで、オペランドが独立している場合 (または独立させることが可能な場合)、これらのループは VPU に適した有力な DLP 並列の候補となります。DLP 候補が見つかったら、それらを配列表記で表現して、コンパイラーが SIMD 演算を認識できるようにします。そして、コンパイラーが配列表記を解析し、大きな配列操作を VPU のネイティブな幅の小さな SIMD 演算にストリップマイニングするのに任せます。ここでは、4 つの 32 ビット浮動小数点データ要素を処理するベクトル命令に対応したインテル® Core™ i7 プロセッサー上でベンチマーク・テストを実施します。

また、複数のコア (それぞれに VPU がある) の利点をさらに活用するため、インテル® Cilk™ Plus のタスクメカニズムでタスクレベルの並列性を表します。。図 7 は、DLP (データレベルの並列性) と TLP (タスクレベルの並列性) を実装する AOBench の変換プロセスです。ここからは、AOBench をインテル® Cilk™ Plus に移植するプロセスについて述べます。

図 7: インテル® Cilk™ Plus により DLP を最初に実装し、次に TLP を実装する AOBench の変換

以下は、AOBench をインテル® Cilk™ Plus に変換するステップです。

1. 実行時間のほとんどが費やされている候補関数またはコードブロックを見つけます。
2. ループ伝播の依存性がない (または依存性を排除できる) 候補ループを見つけます。
3. 変換後も既存のデータ依存性を保持できる安全な命令を選択します。
4. 安全な命令にある定義関数が配列形式でない場合は、データレイアウトを変換し、使用されている箇所をすべて置換します。
5. 配列表記を利用して、安全な命令にある単一のデータ参照を複数のデータ参照に拡張します。
6. インテル® Cilk™ Plus のタスクメカニズムでタスクレベルの並列性を表現し、配列表記を処理する複数の独立した反復を並列に実行します。

アルゴリズム 1: 関数の配列表記の生成

hotspot 解析 (インテル® VTune Amplifier^[3] のようなパフォーマンス・アナライザーを使用するか、コードにインストルメントを行うプローブを挿入する) から、ambient_occlusion() で実行時間の 90% が費やされており、関数内に候補となる内部ループがあるため、この関数が候補関数であることが分かります。候補ループの最初の 8 つの命令は、8 つのスカラー変数の定義で、ループ内でのみ使用されます。theta 変数と phi 変数は標準の算術ライブラリーを用いて計算され、その他の変数はこれらの変数から計算されます。8 つの命令にはループ伝播の依存性がないため (前後のループ反復の結果に依存していないため)、安全に変換できます。これらの命令により処理されるデータは配列形式ではないため、命令で定義されているスカラー変数ごとにローカル配列を作成し、配列表記を適用する必要があります。図 8 のシリアルコードに適用されるアルゴリズム 1 に従って、ループの外側で 8 つの配列変数を宣言し、対応するスカラー変数の定義と使用箇所をループの反復に応じた特定の配列要素に置換します。

ここで、配列サイズはループの反復回数によって決定されます。単一のデータ操作を複数のデータ操作に変換するには、配列全体にわたり配列表記でデータ並列性を表現します。コンパイラーは、VPU ベクトル命令に渡されるデータのチャンク (塊) を作成することで配列全体を VPU にマップするため、命令はループに依存しません。単純に命令をループの外側に移動することができます。

図 9 は等価な配列表記コードです。theta と phi は乱数から計算されるため、シリアルバージョンの結果と一致するように使用前に rand1 と rand2 に格納しておきます。

図 8: ambient_occlusion() のシリアルコード

次に、構造体変数と反復ごとに異なる値で構造体のスカラーメンバーを定義する 6 つの命令を処理します。配列表記にするためスカラー変数を配列変数に変換したように、構造体内部のスカラーメンバー変数を配列に変換する必要があります。そのためには、データ構造をスカラー変数の構造体から配列構造体 (SOA) に変換します。これにより、インテル® Cilk™ Plus の配列表記を利用して、1 つのメンバー変数を定義する 1 つの命令を、一度に複数のメンバー変数を定義する命令に効率良く置換できます。正しく変換するため、構造体メンバーが使われている箇所すべてを、新しい構造体の対応する配列メンバーの適切な要素に置換します。現在、この SOA への変換はプログラマーが手動で行う必要がありますが、インテル® Cilk™ Plus の将来のバージョンでは自動的に処理されるようになる予定です。

図 9: 図 8 のスカラー変数と等価な配列表記


図 10: 構造体変数の配列表記への変換

ここまでの作業で、アンビエント・オクルージョンのほとんどの計算にデータ並列性が適用されました。ループで処理される次の命令は ray_sphere_intersect() 関数の呼び出しです。ループの内側にある内部関数も、関数内の命令がすべて安全な命令であれば、配列表記により並列化の可能性がもたらされます。関数内の命令の変換は、ループの変換とほぼ同じステップで行えます。関数の引数を配列参照に変換し、関数の本体で部分配列表記を利用して、複数の配列要素が並列に処理されるようにします。呼び出し位置は適宜調整します。図 11 は、ray_sphere_intersect() とその呼び出し位置の一部を変換したものです。この場合、関数の引数はスカラー変数からなる構造体で、すでに SOA に変換されています。

すでに ambient_occlusion() の最内ループの命令は、1 つの累積命令を除いて、すべて配列表記に変換されています。この累積命令は以前の反復からのデータに依存しているため、そのままにします。ここで実行時間を測定すると、このアプリケーションはシリアル実行よりも 2.49 倍速いことが分かりました (表 1 を参照)。さらに、一部の関数をインライン展開すると (冗長な命令が削除され)、3.21 倍のスピードアップ (表 1 を参照) が得られます。例えば、ループの反復全体で ray->org と sphere->center は変わらないため、ray_sphere_intersect() で rs を nphi 回計算する代わりに、1 回計算するだけで済みます。

ここまでの作業では、データ並列性のみを適用し、シングルコアのベクトル命令だけを利用しました。4 VPU を並列に実行するクアッドコア/8 スレッドのインテル® Core™ i7 プロセッサーの性能を最大限に活用するため、ここでタスクレベルの並列性にも注目してみましょう。タスクレベルの並列性はインテル® Cilk™ Plus により簡単に実装できます。AOBench の大量の計算はループで行われるため、タスクレベルの並列性を表現するには、候補ループを見つけます。一般に、スレッドの粒度が大きいほど、タスク並列処理によりパフォーマンスが向上します。この点を考慮して、render() の最外ループでタスクをスポーンします。このループには前述の _Cilk_for の制限がなく、このタスクメカニズムを実装しただけでも 6.7 倍のスピードアップが達成されたため、ここでは _Cilk_for キーワードを選択しています。図 12 は _Cilk_for によるコードの変換を示しています。

図 12: インテル® Cilk™ Plus による変換でタスクレベルの並列性をサポート

5. パフォーマンス結果

このセクションでは、インテル® Parallel Studio XE 2011 SP1 による測定結果とシリアルバージョンの比較結果を示します。測定に使われた AOBench のシリアル・ソース・コードは、次の Google* のプロジェクト・ホスティング・サイトから取得しました:

http://code.google.com/p/aobench/ (英語)。

インテル® Cilk™ Plus への変換には、インテル® Parallel Studio XE 2011 SP1 に含まれるインテル® C++ コンパイラーを使用しました。測定したシステムは、インテル® Core™ プロセッサー 2.80GHz、4GB RAM を搭載するものです。表 1 は、AOBench のさまざまな実装 (シリアル、配列表記、インテル® Cilk™ Plus のタスク並列処理、インテル® Cilk™ Plus をフル活用した場合など) の実行時間 (秒単位) とスピードアップを示します。インテル® Cilk™ Plus でデータ並列処理とスレッド並列処理を実装することにより、平均で 16.47 倍のスピードアップが得られました。

 

表 1: AOBench の実行時間とスピードアップ

6. まとめ

この記事では、有名な AOBench ベンチマークを使って、インテル® Cilk™ Plus の設計および実験に基づく評価について述べました。インテル® Cilk™ Plus により、プログラマーはコードの並列ブロックを特定することに集中し、複雑な並列リソースへの割り当てはコンパイラーに任せることで、並列プログラミングを簡単に実現できることを示しました。測定結果から、インテル® Cilk™ Plus はワークスチール・スケジューラーによりマルチコアを効率良く利用し、16 倍のスピードアップを達成できることが分かりました。

参考資料

http://software.intel.com/en-us/intel-cilk-plus (英語)

注

1. このコード例では倍精度を使っていますが、整数および単精度でも同様の方法で表すことができます (Saxpy)。

2. 組込み関数は、コンパイラーにより実装および処理されます。使用法はユーザー定義関数に似ていますが、インライン展開によりパフォーマンスが向上する可能性があります。

3. インテル® VTune™ Amplifier およびインテル® Composer の詳細については、次の Web サイトを参照してください:
http://www.xlsoft.com/jp/products/intel/

* その他の社名、製品名などは、一般に各社の表示、商標または登録商標です。

コンパイラーの最適化に関する詳細は、最適化に関する注意事項を参照してください。