テンポラルブロッキングを適用したステンシル計算コードのSIMD化とルーフラインモデルを用いた性能解析

全文

(1)Vol.2015-HPC-151 No.17 2015/10/1. 情報処理学会研究報告 IPSJ SIG Technical Report. テンポラルブロッキングを適用したステンシル計算コードの SIMD 化とルーフラインモデルを用いた性能解析佐藤真平1,2,a). 佐藤幸紀1,2,b). 遠藤敏夫1,2,c). 概要：大規模並列処理をおこなう場合においても CPU およびメモリに対するチューニングはアプリケーション全体の性能を決める重要な要素である．本稿では，ステンシル計算コードの SIMD 化に注目する．ステンシル計算において，メモリ参照の空間局所性を利用する最適化手法であるテンポラルブロッキングを適用したコードを対象に，性能モデルの 1 つであるルーフラインモデルを用いてステンシル計算の性能を見積もりを試みる．テンポラルブロッキングを適用したコードに SIMD 化に関するチューニングを施し，性能向上を確認する．また，その性能向上をルーフラインモデルを用いて視覚的にとらえられることを示す．. 1. はじめに. ウンまで性能向上を達成したことが報告されている．アプリケーションの性能チューニングを成功させるに. エクサスケール時代に向けた高性能計算システムにおい. は，アプリケーションを実行する計算機システムにおける. て，高速計算を実現するためには数千から数万に及ぶ並. ピーク性能を見積もることが重要である．性能の見積もり. 列性を抽出する大規模並列処理が必要である．一方で，1. には，一般に性能モデルが用いられる．ルーフラインモデ. ノードもしくは CPU 単体におけるアプリケーションの性. ル [2] は性能モデルの 1 つで，アプリケーションの潜在的. 能は，大規模並列化後の性能を決めるベースとなる要素で. なチューニングの可能性とボトルネック解析による性能限. ある．例えば，1 ノードあたりの性能を 1 割ほどしか引き. 界を視覚的に示す性能モデルである．このモデルでは，ア. 出せていないアプリケーションは，大規模並列化を行って. プリケーションの性能を計算量と DRAM へのアクセス量. も同程度しかシステム全体の性能を引き出せないという. で抽象化し，表現する．. ことになる．今後の高性能計算機システムにおいても，1 ノードあたりの性能チューニングは重要な課題である．. 本稿では，ルーフラインモデルによるアプリケーションの性能見積もりとそれに基づく性能チューニングの事例を. 近年の計算機システムにおける 1 ノードあたりの性能. 報告する．対象とするアプリケーションはステンシル計算. チューニングでは，メモリウォール問題のために複雑化す. コードとする．ステンシル計算コードを対象に，テンポラ. るメモリ階層への対応，コア数の増加や SIMD を用いたプ. ルブロッキングと呼ばれるメモリバンド幅の要求を改善す. ロセッサ内における並列性の抽出が求められる．Intel で. る手法 [3], [4], [5] を適用し，さらにコンパイラによる自動. は，C や C++ で単に書かれたコードとマルチコア，メ. SIMD 化を考慮したチューニングを施す．. ニーコアプロセッサ向けに最もチューニングされたコード. ステンシル計算は，流体シミュレーションなどの分野の. の性能差を「Ninja Gap」と称し，コンパイラの支援など. 重要カーネルである．一般に，シミュレーションする領域. を利用したユーザフレンドリーなチューニングでその性能. を格子で表し，それぞれの格子点の値を隣接する格子点の. 差を縮める取り組みを進めている [1]．この取り組みでは，. 値を用いて計算し更新する処理を時間ステップとして繰り. 最もチューニングされたコードと比較して平均で 24 倍の. 返し実行する．ステンシル計算はデータ並列性が高く容易. スローダウンであったコードが，平均で 1.3 倍のスローダ. に並列化が可能である．しかし，その性能はメモリバンド. 1. 幅に強く律速されることが知られており，並列化を行って. 2 a) b) c). 東京工業大学学術国際情報センター Global Scientific Information and Computing Center, Tokyo Institute of Technology JST CREST [email protected] [email protected] [email protected]. ⓒ 2015 Information Processing Society of Japan. も十分な性能向上が得られない場合がある．そこで，ルーフラインモデルを用いて計算機におけるステンシル計算コードのピーク性能を見積もることを試みる．同時に，最適化の適用前後のメモリバンド幅への要求. 1.

(2) Vol.2015-HPC-151 No.17 2015/10/1. 情報処理学会研究報告 IPSJ SIG Technical Report 表 1. の変化，SIMD 化による性能向上を視覚的にとらえられることを実証する．. 2. ルーフラインモデルを用いたステンシル計算コードの性能見積もり. 実験に使用する計算機環境. プロセッサ. Intel Xeon E5-2650 v2 × 2. コア数. 8×2. 動作周波数. 2.6 GHz. 単精度浮動小数点演算. 83.2 GFlops/s × 2. ピーク性能（SSE 命令）. ステンシル計算コードについて，我々が使用している計. L1 データキャッシュ. 32 KB, 8-way / core. 算機環境における性能見積もりを行う．まず我々の計算機. L2 キャッシュ. 256 KB, 8-way / core. 環境におけるルーフラインモデルを示し，次にそのモデル. L3 キャッシュ. 20 MB, 20-way / CPU. 上でステンシル計算コードの性能を見積もる．. メモリ. DDR3 1866 MHz 64 GB. コンパイラ. icc 14.0.2. HW カウンタ取得ツール. likwid 3.1.3. 2.1 ルーフラインモデルメモリバンド幅を考慮した性能モデルとしてルーフライ. Roofline model (Xeon E5-2650 V2, 2 CPU). ンモデル [2] が知られている．このモデルでは，演算性能. ^^ƉĞĂŬƉĞƌĨŽƌŵĂĐĞ͕ϭϲϲ͘ϰ'&ůŽƉƐͬƐ. とメモリバンド幅のどちらか一方がアプリケーションの能，演算強度（Operational Intensity）を用いてアプリケーション性能を見積もる．演算強度は，アプリケーションにおける浮動小数点演算量と DRAM アクセスのデータ量の. ŶŽͲǀĞĐƉĞĂŬƉĞƌĨŽƌŵĂĐĞ͕ϰϭ͘ϲ'&ůŽƉƐͬƐ. Performance (GFlops/s). 性能を律速すると仮定し，メモリバンド幅，ピーク演算性. 128. 32. 8. 2. 比（Flops/Byte）である．仮定から，ルーフラインモデル. 0.5. roofline, no-vec roofline, SSE. による演算性能は次の式で表される．. 0.125 0.0078125 0.03125. 図 1. 0.125 0.5 2 8 Operational Intensity (Flops/DRAM Byte). 32. 128. 実験に使用する計算機環境におけるルーフラインモデル. 演算性能 = min{ ピーク演算性能, メモリバンド幅×演算強度 }. (1) この式は，横軸に演算強度，縦軸に演算性能として図示. キャッシュは L1，L2 がプライベート，L3 が共有となっている．メモリは，DDR3 1866 MHz の 8 GB のモジュール. すると，メモリバンド幅に律速され演算強度が高くなるに. を各 CPU に 4 枚ずつ搭載し，4 チャネルすべてを利用し. つれて演算性能が向上する線と CPU ピーク性能に律速さ. ている．コンパイラはインテルコンパイラ 14.0.2 を使用. れ，演算強度にかかわらず演算性能が一定になる線で表現. する．最適化オプションはすべての実験において O3 を使. される．. 用する．実験において，アプリケーション実行時のメモリ. ルーフラインモデルにおけるメモリバンド幅とピーク演. バンド幅などはハードウェアカウンタを用いて計測する．. 算性能は適切に設定する必要がある．例えばメモリバンド. ハードウェアカウンタの取得には likwid 3.1.3 [6] を利用. 幅については，実行するスレッド数やメモリ構成によって. する．. 変わる．ピーク演算性能についても，使用する計算機環境や SIMD の使用，不使用などによって変化する．. likwid は，ハードウェアカウンタの取得のみならず，スレッド・アフィニティーの制御も行うことができるツールである．以降のすべての実験において，スレッドはコアに. 2.2 実験環境におけるルーフラインモデル. 固定で割り当てて実行する．. ここでは，本稿でアプリケーションを実行する計算機環. 図 1 に，実験で使用する計算機環境におけるルーフライ. 境についてまとめ，その環境におけるルーフラインモデル. ンモデルを示す．モデルは，2 CPU（16 スレッド）実行時. を示す．. の SSE および SIMD 不使用（no-vec）の 2 つの場合につ. 表 1 に実験で使用する計算機環境をまとめる．この計算. いて図示している．ピーク演算性能は，単精度浮動小数点. 機は Ivy Bridge 世代の CPU，Intel Xeon E5-2650 v2 を. 演算のピーク性能を用い，SSE の場合は 166.4 GFlops/s，. 2 基搭載する NUMA 型のマシンである．コア数は 1 CPU. no-vec の場合は 41.6 GFlops/s とする．式 1 より，図中. あたり物理 8 コア，論理 16 コアである．実験では，ハイ. の斜線の傾きはメモリバンド幅となる．メモリバンド幅は. パースレッディングは使用せずに，1 コアに 1 スレッド. STREAM ベンチマーク [7] で計測した Triad の結果から. を割り当ててアプリケーションを実行する．動作周波数は. 81.7 GB/s とする．. 2.6 GHz で，ターボブーストは不使用とする．単精度の浮. 以降の実験において，スレッド数は 2 CPU を使用し 16. 動小数点演算性能は，SSE 命令の使用で， CPU あたり. スレッド固定とする．浮動小数点演算については，SIMD. 4 F lops × 8 cores × 2.6 GHz = 83.2 GF lops/s となる．. は SSE 命令のみを用いる．使用する計算機では浮動小数. ⓒ 2015 Information Processing Society of Japan. 2.

(3) Vol.2015-HPC-151 No.17 2015/10/1. 情報処理学会研究報告 IPSJ SIG Technical Report. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47. for ( int t = 0; t < TMAX ; ++ t ) { // time step // update for ( int z = 1; z < ZMAX -1; ++ z ) { for ( int y = 1; y < YMAX -1; ++ y ) { for ( int x = 1; x < XMAX -1; ++ x ) { # if defined POINTS_27 // 27 point stencil * APTR ( dst , x , y , z ) = cc * * APTR ( src , x , y , z ) + cx * (* APTR ( src , x -1 , y , z ) + * APTR ( src , x +1 , y , z )) + cy * (* APTR ( src , x , y -1 , z ) + * APTR ( src , x , y +1 , z )) + cz * (* APTR ( src , x , y , z -1) + * APTR ( src , x , y , z +1)) + cxy * (* APTR ( src , x -1 , y -1 , z ) + * APTR ( src , x -1 , y +1 , z ) + * APTR ( src , x +1 , y -1 , z ) + * APTR ( src , x +1 , y +1 , z )) + cyz * (* APTR ( src , x , y -1 , z -1) + * APTR ( src , x , y -1 , z +1) + * APTR ( src , x , y +1 , z -1) + * APTR ( src , x , y +1 , z +1)) + czx * (* APTR ( src , x -1 , y , z -1) + * APTR ( src , x -1 , y , z +1) + * APTR ( src , x +1 , y , z -1) + * APTR ( src , x +1 , y , z +1)) + cxyz * (* APTR ( src , x -1 , y -1 , z -1) + * APTR ( src , x -1 , y -1 , z +1) + * APTR ( src , x -1 , y +1 , z -1) + * APTR ( src , x -1 , y +1 , z +1) + * APTR ( src , x +1 , y -1 , z -1) + * APTR ( src , x +1 , y -1 , z +1) + * APTR ( src , x +1 , y +1 , z -1) + * APTR ( src , x +1 , y +1 , z +1)); # elif defined POINTS_19 // 19 point stencil * APTR ( dst , x , y , z ) = cc * * APTR ( src , x , y , z ) + cx * (* APTR ( src , x -1 , y , z ) + * APTR ( src , x +1 , y , z )) + cy * (* APTR ( src , x , y -1 , z ) + * APTR ( src , x , y +1 , z )) + cz * (* APTR ( src , x , y , z -1) + * APTR ( src , x , y , z +1)) + cxy * (* APTR ( src , x -1 , y -1 , z ) + * APTR ( src , x -1 , y +1 , z ) + * APTR ( src , x +1 , y -1 , z ) + * APTR ( src , x +1 , y +1 , z )) + cyz * (* APTR ( src , x , y -1 , z -1) + * APTR ( src , x , y -1 , z +1) + * APTR ( src , x , y +1 , z -1) + * APTR ( src , x , y +1 , z +1)) + czx * (* APTR ( src , x -1 , y , z -1) + * APTR ( src , x -1 , y , z +1) + * APTR ( src , x +1 , y , z -1) + * APTR ( src , x +1 , y , z +1)); # else // 7 point stencil * APTR ( dst , x , y , z ) = ( cc * * APTR ( src , x , y , z ) + cx * (* APTR ( src , x -1 , y , z ) + * APTR ( src , x +1 , y , z ))) + ( cy * (* APTR ( src , x , y -1 , z ) + * APTR ( src , x , y +1 , z )) + cz * (* APTR ( src , x , y , z -1) + * APTR ( src , x , y , z +1))); # endif } } } }. 図 2. 実験に用いるステンシル計算の疑似コード. 表 2. 点命令を乗算と加算に限り 2 命令同時実行することができるが，本稿では 1 命令実行の場合をピーク性能として議論. ステンシル計算の 16 スレッド実行時の性能 GFlops/s GB/s. 7 point stencil. を進める．. 2.3 ステンシル計算コード. SSE. 35.0. 68.7. no-vec. 32.7. 64.6. 19 point stencil. 図 2 に実験で用いるステンシル計算の疑似コードを示. SSE. 84.2. 66.6. す．このステンシル計算は，単精度浮動小数点のサイズ. no-vec. 39.8. 31.5. XM AX × Y M AX × ZM AX の 3 次元配列に対して，7. 27 point stencil SSE. 106.7. 61.1. no-vec. 37.3. 22.2. 点，19 点，27 点のいずれかのステンシル計算を TMAX ステップ行う．コードの記述言語は C 言語である．配列 src および dst は実行時に動的に Linear Array*1 で確保する．. Roofline model (Xeon E5-2650 V2, 16 threads). 配列の要素へのアクセスは，インデックス計算により最内. ^^ƉĞĂŬƉĞƌĨŽƌŵĂĐĞ͕ϭϲϲ͘ϰ'&ůŽƉƐͬƐ. ループの次元である X 次元で連続になるようしている．チマーク [8] では，要素値が格納されている配列の他に係数のために複数の配列が宣言されており，合計で 10 本ほどの配列が用いられている．これにより，キャッシュにお. ŶŽͲǀĞĐƉĞĂŬƉĞƌĨŽƌŵĂĐĞ͕ϰϭ͘ϲ'&ůŽƉƐͬƐ. Performance (GFlops/s). ステンシル計算のベンチマークとして知られる姫野ベン. 128. 32. 8. 7pt no-vec. では，配列は要素値が格納されている 2 つの 3 次元配列のみである．したがって，キャッシュにおける競合ミスを考. 19pt SSE. 2. 19pt no-vec. ける競合ミスが発生し十分な性能を達成するには適切なパディングが必要となる [9]．本稿で用いるステンシル計算. 7pt SSE. 27pt SSE. 0.5. 27pt no-vec 0.125 0.0078125 0.03125. 0.125 0.5 2 8 Operational Intensity (Flops/DRAM Byte). 32. 128. 図 3 ステンシル計算の 16 スレッド実行時のルーフラインモデル. 慮する必要がなく，パディングは実施していない．点，19 点，27 点ステンシルそれぞれについて行う．また，. 2.4 ルーフラインモデルによるステンシル計算の性能見積もり前述の計算機を用いて，ステンシル計算を実行しルーフラインモデルに基づき性能を見積もる．性能の計測は，7. 比較のために SIMD を用いた場合と用いない場合の 2 通りの計測を行う．サイズは 512 × 512 × 512 とする．以降の実験でも，ステンシル計算のサイズはこのサイズを用いる．表 2 に 7 点，19 点，27 点ステンシルの 16 スレッド実行. *1. Linear Array は，1 次元配列を多次元配列とみなして扱う配列．本稿のステンシル計算では malloc(sizeof(float) ×XM AX × Y M AX × ZM AX) のように配列を確保している．. ⓒ 2015 Information Processing Society of Japan. 時の性能とメモリバンド幅の計測結果を示す．また，図 3 にそれぞれの場合についてルーフラインモデルにプロット. 3.

(4) Vol.2015-HPC-151 No.17 2015/10/1. 情報処理学会研究報告 IPSJ SIG Technical Report. する．ルーフラインは，SSE を使用した場合のピーク性能と SSE を使用しない場合のピーク性能を描画している．7. ĐŽƉǇ. 点ステンシルは黒，19 点ステンシルは青，27 点ステンシ ĐŽƉǇ. ルは赤のマーカーで性能を示している．また，SSE を使用した時の性能（SSE）を四角，SSE を使用しない時の性能（no-vec）を三角のマーカーで示している．. 7 点ステンシルにおいて，演算強度が小さく，SSE の使用にかかわらず性能がメモリバンド幅律速になり 35. GFlops/s 程度となることがわかる．19 点，27 点ステンシルにおいて，SSE を使用しない場合は演算性能律速とな ǆсϬ. りピーク演算性能のあたりの性能となっていることがわか. ƚŝŵĞďůŽĐŬƐƚĞƉƐ. る．これら 2 つについても，SSE を使用した場合はメモリバンド幅律速となり，演算強度に応じた性能となっている．. 図 4. ƚŝŵĞďůŽĐŬƐƚĞƉƐ. 冗長な計算が発生するテンポラルブロッキングの実装例. SIMD の使用による演算強度の変化は小さく，図 3 で示されている通り SSE 命令を使用した場合でもメモリバン. ĐŽƉǇ. ド幅律速となり十分に性能向上しないことがわかる．さらなる高速化として AVX 命令を使用したとしても，SSE と. ĐŽƉǇ. 同様にメモリバンド幅律速となり十分な性能向上は達成し ĐŽƉǇ. ないことが推測できる．このステンシル計算を高速化するには，メモリバンド幅への要求を減らすアルゴリズムの変更が必要である．. 3. SIMD 化を考慮したテンポラルブロッキング. ǆсϬ ƚŝŵĞďůŽĐŬƐƚĞƉƐ. 3.1 テンポラルブロッキングテンポラルブロッキングは，ステンシル計算のメモリバ. 図 5. ƚŝŵĞďůŽĐŬƐƚĞƉƐ. 冗長な計算がないテンポラルブロッキングの実装例. ンド幅要求を減らす最適化手法 [3], [4], [5] である．メモリアクセスの局所性を利用し，キャッシュの利用率を高める. ブロックのステップにおける計算結果を再利用しつつ計算. 手法として空間ブロック分割が知られている．テンポラル. 範囲をずらしながら時間ブロックの計算を進める．図の例. ブロッキングは，さらにメモリアクセスの局所性を利用す. では，全体の配列から必要な格子点をコピーした後の時間. るために空間ブロックに分割した領域ごとに時間ステップ. ステップの計算において，前の空間ブロックの計算の結果. を数ステップ分計算する．以降では，テンポラルブロッキ. から次のステップのブロックの計算に必要な格子点（水色. ングで進める時間ステップを時間ブロックと呼ぶ．. で示した格子点）をコピーしながら計算を進める．これに. テンポラルブロッキングの実装は 2 種類ある．時間ステップを進めるにあたって冗長な計算が発生する場合 [3] としない場合 [5] である．図 4 に冗長な計算が発生する実装例，図 5 に冗長な計算がない実装例を示す．. より，計算する格子点の数はテンポラルブロッキングを用いない場合と同じになる．冗長な計算が発生する実装は，計算する格子点が多くなるが，格子点の計算がはじめに確保した空間ブロックの範. 図 4 は，冗長な計算が発生する実装において，1 次元の. 囲から外れることはないためメモリアクセスの局所性が高. 配列に対して時間ブロックの計算を進める様子を示してい. い．一方で，冗長な計算がない実装は，計算する格子点の. る．緑で示した空間ブロックの結果を得るためには，その. 数はテンポラルブロッキングを用いない場合と変わらない. 空間ブロックサイズより大きな空間からピラミッド状に時. が，計算範囲がずれていくためにメモリアクセスの局所性. 間ブロックの計算を進める．このとき，前の空間ブロック. が冗長な計算が発生する実装ほど高くない．しかし，テン. のステップにおいて同じ計算を行っている格子点が存在す. ポラルブロッキングを用いない場合と比較すれば十分にメ. る．この実装では時間ブロックサイズを大きくすると，冗. モリアクセスの局所性は高く，本稿では冗長な計算がない. 長に計算する格子点が多くなる．. 実装を採用する．. 図 5 は，冗長な計算がない実装において，1 次元の配列に対して時間ブロックの計算を進める様子を示している．緑で示した空間ブロックの結果を得るためには，前の時間. ⓒ 2015 Information Processing Society of Japan. 3.2 SIMD 化を考慮したテンポラルブロッキング SIMD を使用することで浮動小数点演算を 1 度にメモリ. 4.

(5) Vol.2015-HPC-151 No.17 2015/10/1. 情報処理学会研究報告 IPSJ SIG Technical Report. 表 3 7 点ステンシル計算のテンポラルブロッキングの評価. ĐŽƉǇ. Spatial Block Temporal GFlops/s GB/s Temporal SIMD GFlops/s GB/s Temporal Proposal GFlops/s GB/s. ĐŽƉǇ. ĐŽƉǇ. (512, 64, 64). (256, 256, 32). (128, 128, 128). 44.7 8.9. 37.9 20.8. 22.0 40.5. 94.9 16.7. 51.0 29.9. 24.2 42.9. 90.2 18.4. 52.7 28.9. 25.6 43.4. Roofline model (Xeon E5-2650 V2, 16 threads) ^^ƉĞĂŬƉĞƌĨŽƌŵĂĐĞ͕ϭϲϲ͘ϰ'&ůŽƉƐͬƐ. Performance (GFlops/s). 128. ǆсϬ ƚŝŵĞďůŽĐŬƐƚĞƉƐ. 図 6. ƚŝŵĞďůŽĐŬƐƚĞƉƐ. 冗長な計算がないテンポラルブロッキングで 4 要素ずつずら. ŶŽͲǀĞĐƉĞĂŬƉĞƌĨŽƌŵĂĐĞ͕ϰϭ͘ϲ'&ůŽƉƐͬƐ. 32. 8 Temporal 2. Temporal SIMD Temporal Proposal. 0.5. す実装例. Non Temporal Non Temporal SIMD. 0.125 0.0078125 0.03125. 上に連続に並んだ複数のデータに対して実行することができる．例えば SSE 命令を使用すると，ステンシル計算においては 4 つの格子点の計算を 1 度に行うことができる．. 図7. 0.125 0.5 2 8 32 Operational Intensity (Flops/DRAM Byte). 128. テンポラルブロッキングを適用した 7 点ステンシル計算のルーフラインモデル. 冗長な計算がないテンポラルブロッキングの実装では，計算範囲を 1 要素ずつずらしながら時間ブロックの計算を進める．これに対して SIMD 化を行うと配列の両端のブロックにおいて SIMD で計算できない格子点ができてしまう．. 4. 評価ここでは，ステンシル計算へのテンポラルブロッキング. 高速化には，なるべく多くの格子点を SIMD で計算するこ. の適用および SIMD 化についてルーフラインモデルに基. とが必要である．. づきその性能を評価する．評価対象は以下の通りである．. SIMD を使用した浮動小数点演算では，1 度に計算する. • Temporal：テンポラルブロッキングを適用したコー. データのメモリアライメントが整っていると効率よくメモ. ド．SIMD 化のためのディレクティブの挿入は行なわ. リ参照でき高い演算性能を達成することができる [10]．例. ずにコンパイルする．. えば，SSE 命令では 16 バイト境界にデータが整列してい. • Temporal SIMD：Temporal のコードにコンパイラ. るとより演算性能が高くなる．テンポラルブロッキングに. による SIMD 化のためのディレクティブを挿入した. おいて 1 要素ずつずらしたとき，空間ブロックの計算を. コード．malloc によるメモリアライメントを調整し，. SIMD で実行してもデータのメモリアライメントが 16 バ. 計算する格子点を 16 バイト境界に整列．. イト境界に整列していない演算になってしまう．高速化に. • Temporal Proposal：SSE による SIMD 化を考慮. は，なるべく多くの SIMD 演算のメモリアライメントが. し，4 要素ずつずらすテンポラルブロッキングを適用. 整っていることが必要である．. したコード．. 以上のことを踏まえて，SIMD 化を考慮したテンポラル. ステンシル計算を行う全体のサイズは 512 × 512 × 512 と. ブロッキングとして，ずらす要素数を SIMD の演算数に合. する．テンポラルブロッキングにおける空間ブロックサイズ. わせるチューニングを実施する．この実装では，SSE 命令. は，(x, y, z) = (512, 64, 64)，(256, 256, 32)，(128, 128, 128). の場合は 4 要素，AVX 命令の場合は 8 要素ずらしながら. の 3 通り，時間ブロックサイズは，128 とする．空間ブ. 時間ブロックの計算を進める．. ロックサイズは，すべて要素数が等しくなっている．. 図 6 に SSE 命令の使用を考慮して，冗長な計算のない. すべての評価は 16 スレッドで行う．実験に使用する計. テンポラルブロッキングにおいて，ずらす要素数を 4 とし. 算機は NUMA 型のマシンのため，スレッドは 2 つの CPU. たときの実装例を示す．4 要素ずらすことで，SSE 命令を. に 8 スレッドずつ均等に分配し，コアに固定して割り当て. 使用したときに配列の終端以外すべて SIMD で格子点の. る．また，メモリを確保した配列に対してファーストタッ. 計算を行うことができる．また，配列の先頭においてメモ. チを実施し，ロードバランスを均等にしている．. リアライメントを整えると，すべての SIMD による計算がアライメントが整ったデータに対して行われる．この実装. 4.1 7 点ステンシル計算の評価. は，SIMD による演算性能の向上が見込まれるが，前のス. 表 3 に，7 点ステンシル計算における Temporal，Tem-. テップでの計算結果をより多く再利用する必要があり，メ. poral SIMD，Temporal Proposal の性能をまとめる．表 2. モリアクセスの局所性が低くなるという欠点がある．. において，SIMD 化なしで 32.7 GFlops/s の性能であった. ⓒ 2015 Information Processing Society of Japan. 5.

(6) Vol.2015-HPC-151 No.17 2015/10/1. 情報処理学会研究報告 IPSJ SIG Technical Report 表 4 19 点ステンシル計算のテンポラルブロッキングの評価 (512, 64, 64). (256, 256, 32). (128, 128, 128). 37.6 3.4. 35.8 8.7. 30.8 24.5. 119.0 9.2. 96.4 21.9. 51.2 39.6. 120.3 10.6. 96.8 22.7. 50.6 38.9. Roofline model (Xeon E5-2650 V2, 16 threads) ^^ƉĞĂŬƉĞƌĨŽƌŵĂĐĞ͕ϭϲϲ͘ϰ'&ůŽƉƐͬƐ. 128 Performance (GFlops/s). Spatial Block Temporal GFlops/s GB/s Temporal SIMD GFlops/s GB/s Temporal Proposal GFlops/s GB/s. 表 5 27 点ステンシル計算のテンポラルブロッキングの評価 Spatial Block Temporal GFlops/s GB/s Temporal SIMD GFlops/s GB/s Temporal Proposal GFlops/s GB/s. (512, 64, 64). (256, 256, 32). (128, 128, 128). 36.2 2.5. 35.3 3.0. 32.4 19.4. 114.5 6.6. 102.4 17.1. 65.7 37.6. 116.0 7.6. 104.2 18.2. 68.2 39.2. ŶŽͲǀĞĐƉĞĂŬƉĞƌĨŽƌŵĂĐĞ͕ϰϭ͘ϲ'&ůŽƉƐͬƐ. 32. 8 Temporal. 2. Temporal SIMD Temporal Proposal. 0.5. Non Temporal Non Temporal SIMD. 0.125 0.0078125 0.03125. 図 8. 0.125 0.5 2 8 32 Operational Intensity (Flops/DRAM Byte). 128. テンポラルブロッキングを適用した 19 点ステンシル計算のルーフラインモデル Roofline model (Xeon E5-2650 V2, 16 threads) ^^ƉĞĂŬƉĞƌĨŽƌŵĂĐĞ͕ϭϲϲ͘ϰ'&ůŽƉƐͬƐ. 128 Performance (GFlops/s). ステンシル計算が，空間ブロックサイズが (512, 64, 64) の時に最も良い性能を示し，SIMD 化を含めて約 3 倍の性能向上を達成している．図 7 に空間ブロックサイズ (512, 64, 64) の時のそれぞれの性能をルーフラインモデルに図示する．また，図中に赤. Temporal. 2. Temporal SIMD Temporal Proposal Non Temporal Non Temporal SIMD. 0.125 0.0078125 0.03125. Temporal，Non Temporal SIMD）を示す．図から，テンポリバンド幅への要求が改善し，演算強度が高くなっている. 8. 0.5. のマーカーでテンポラルブロッキング適用前の性能（Non ラルブロッキングを適用することでステンシル計算のメモ. ŶŽͲǀĞĐƉĞĂŬƉĞƌĨŽƌŵĂĐĞ͕ϰϭ͘ϲ'&ůŽƉƐͬƐ. 32. 図 9. 0.125 0.5 2 8 32 Operational Intensity (Flops/DRAM Byte). 128. テンポラルブロッキングを適用した 27 点ステンシル計算のルーフラインモデル. ことがわかる．テンポラルブロッキング適用前では SIMD 化をしても性能向上が望めなかったが，適用後は SIMD 化により性能向上が望める演算強度となっている．. 4.2 19 点，27 点ステンシル計算の評価表 4，表 5 に，19 点，27 点ステンシル計算における. SIMD を考慮したテンポラルブロッキング（Temporal. Temporal，Temporal SIMD，Temporal Proposal の性能を. Proposal）については，空間ブロックサイズ (512, 64, 64). まとめる．表 2 において，SIMD 化なしで 39.8 GFlops/s. では性能が向上しなかった．これは，4 要素ずらす X の. の性能であった 19 点ステンシル計算が，空間ブロックサ. 次元のサイズが 512 で，空間分割されていないことが原. イズが (512, 64, 64) の時に最も良い性能を示し，SIMD 化. 因としてあげられる．Temporal SIMD の実装においても，. を含めて約 3 倍の性能向上を達成している．同様に，表 2. 先頭の空間ブロックはメモリアライメントを 16 バイト境. において，SIMD 化なしで 37.3 GFlops/s の性能であった. 界に整列させている．そのため，X の次元を分割しない場. 27 点ステンシル計算が，SIMD 化を含めて約 3 倍の性能. 合，先頭の空間ブロックで計算する格子点が多く，多くの. 向上を達成している．. SIMD 演算がメモリアライメントの整ったデータに対して. 図 8，図 9 に空間ブロックサイズ (512, 64, 64) の時の 19. 行われる．Temporal Proposal の実装では，より多くの格. 点，27 点それぞれの性能をルーフラインモデルに図示す. 子点について SIMD で計算をすることができ，またメモ. る．また，図中に赤のマーカーでテンポラルブロッキング. リアライメントが整っているデータが多い．しかし，メモ. 適用前の性能（Non Temporal，Non Temporal SIMD）を. リバンド幅への要求が悪化することによる性能低下が上回. 示す．図から，テンポラルブロッキングを適用することで. り，空間ブロックサイズ (512, 64, 64) では性能が向上しな. ステンシル計算のメモリバンド幅への要求が改善し，演算. かったと考えられる．. 強度が高くなっていることがわかる．19 点，27 点ステンシ. 表 3 において，X の次元を分割しており，SIMD 化に. ル計算については，SIMD 化によりある程度の性能向上が. よる性能向上が見込める空間ブロックサイズ (256, 256, 32). 望めたが，メモリバンド幅律速により性能向上が制限され. では，Temporal Proposal が Temporal SIMD より高い性. ていた．テンポラルブロッキングを適用することで SIMD. 能を示している．これは，Temporal Proposal によるアラ. 化により十分な性能向上が望める演算強度となっている．. イメントが整ったデータの SIMD 演算による性能向上が，. 19 点，27 点ステンシル計算では，SIMD を考慮したテ. メモリバンド幅への要求が悪化することによる性能低下を. ンポラルブロッキング（Temporal Proposal）が空間ブロッ. 上回ったと考えられる．. クサイズ (512, 64, 64) において Temporal SIMD と比較し. ⓒ 2015 Information Processing Society of Japan. 6.

(7) Vol.2015-HPC-151 No.17 2015/10/1. 情報処理学会研究報告 IPSJ SIG Technical Report. て若干高い性能を示している．これは，7 点ステンシルと. [5]. 比べて格子点の計算が多く，アライメントが整ったデータの SIMD 演算による性能向上が，メモリバンド幅への要求が悪化することによる性能低下を上回ったと考えられる．. 5. まとめと今後の課題. [6] [7] [8] [9]. エクサスケール時代に向けた高性能計算システムにおいて，高速計算を実現するために大規模並列処理を行う場合においても，1 ノードあたりの性能チューニングはアプリケーション全体の性能を決める重要な要素である．本稿では，ルーフラインモデルによるステンシル計算. [10]. 南武志，岩下武史，中島浩：冗長な計算を伴わない 3 次元 FDTD 法の時空間タイリング，情報処理学会論文誌コンピューティングシステム（ACS），Vol. 6, No. 1, pp. 56–65 (2013). Likwid: https://github.com/rrze-likwid/likwid. STREAM Benchmark: https://www.cs.virginia.edu/stream/. Himeno Benchmark: http://accc.riken.jp/2444.htm. 佐藤真平，佐藤幸紀，遠藤敏夫：ルーフラインモデルによる性能幅推定とステンシル計算コードにおけるメモリレイアウト最適化による性能最大化，情報処理学会研究報告 2015-ARC-216， No. 32, pp. 1 – 6 (2015). A Guide to Vectorization with Intel C++ Compilers: https://software.intel.com/en-us/articles/a-guide-toauto-vectorization-with-intel-c-compilers?language=en.. コードの性能見積もりとそれに基づく性能チューニングの事例を報告した．ステンシル計算に，テンポラルブロッキングと呼ばれるメモリバンド幅への要求を改善する手法を適用し，さらにコンパイラによる SIMD 化を考慮したチューニングを施した．まずはじめに，ステンシル計算についてテンポラルブロッキング適用前は SIMD 化による性能向上が十分に得られないことをルーフラインモデルに基づき確認した．次に評価では，テンポラルブロッキングおよび SIMD 化を考慮したチューニングを施したコードが，SIMD 化前と比較して約 3 倍の性能向上を達成することを示した．また，ルーフラインモデルからテンポラルブロッキングによるメモリバンド幅への要求が改善していることを確認した．SIMD 化を考慮したチューニングでは，演算量が多い 19 点，27 点ステンシルにおいて性能が向上していることを確認した．今後の課題としては，AVX による SIMD 化，テンポラルブロッキングの姫野ベンチマークへの適用，SIMD 化を考慮したテンポラルブロッキングのプリフェッチを用いた高速化，などが挙げられる．参考文献 [1]. [2]. [3]. [4]. Satish, N., Kim, C., Chhugani, J., Saito, H., Krishnaiyer, R., Smelyanskiy, M., Girkar, M. and Dubey, P.: Can Traditional Programming Bridge the Ninja Performance Gap for Parallel Computing Applications?, In Proceedings of the 39th Annual International Symposium on Computer Architecture (ISCA ’12), pp. 440–451 (2012). Williams, S., Waterman, A. and Patterson, D.: Roofline: An Insightful Visual Performance Model for Multicore Architectures, Communications of the ACM, Vol. 52, No. 4, p. 65 (2009). Wolfe, M.: More Iteration Space Tiling, Proceedings of the 1989 ACM/IEEE Conference on Supercomputing, Supercomputing ’89, New York, NY, USA, ACM, pp. 655–664 (online), DOI: 10.1145/76263.76337 (1989). Nguyen, A., Satish, N., Chhugani, J., Kim, C. and Dubey, P.: 3.5-D Blocking Optimization for Stencil Computations on Modern CPUs and GPUs, In Proceedings of 2010 International Conference for High Performance Computing, Networking, Storage and Analysis (SC ’10), pp. 1–13 (online), DOI: 10.1109/SC.2010.2 (2010).. ⓒ 2015 Information Processing Society of Japan. 7.

(8)