POPプロジェクト事例：一般個体輸送モデルによる粒子トラッキングコードGITMの性能分析

１．背景

お客様：Johan Van Der Molen博士
アプリケーション：GITM
言語：Fortran90
プログラミングモデル：OpenMP, netCDF形式

２．アプリケーション構造

 個体行動モデル（IBM）GITM（一般個体輸送モデル）コードは、オフラインの粒子トラッキングモデルであり、流体力学的モデルからの速度場を用います。これは、物理的な遺留、拡散、および生物学的な発展と行動を含んでいます。この手法は、粒子が流線に正確に従うことを保証します。さらに、拡散を模擬するために、移流補正を伴うランダムウォーク手法も用いています。GITMの生物学的な発展と行動モジュールにより、ユーザが定義した卵数と幼虫成長ステージ数を用いて粒子を発展させることが出来ます。

 移流・拡散コードであるGITMコードは、初期条件をNetCDFファイルから読み込み、解の時間発展を計算します。時間積分の間、各時間ステップでの解はNetCDFファイルに追加されていきます。このコードはOpenMP並列化されており、インテルコンパイラ16.0.2でコンパイルされています。計算に用いるマシンは、Intel Xeon E5-2670 CPU 2.60 GHzです。ハイパースレッディング機構が利用可能な状態にセットされていましたが、コンパイルフラグでは抑制されていました：

-par-affinity="verbose,granularity=fine,compact,1,0"

あるいは、環境変数をコード実行前に設定して、同等の設定を得ることもできます：

export KMP_AFFINITY="verbose,granularity=fine,compact,1,0"

ここではCPUの論理コア0-15を用います。本レポートでは、コードの入力ファイルに、プロファイリング目的のためシミュレーション時間を短めの1時間に設定します。16スレッドで実行させるため、環境変数で、export OMP NUM THREADS=16と設定しました。

３．関心領域(ROI)

 関心領域は、ディスクから入力条件を読み込んだ後の時間積分であることから、この計算のみに焦点を当てます。OpenMPモデルではI/Oは並列化出来ないため、ここでは扱いません。Intel VTuneで計測したところ、I/O時間は0.0165%程度であったため、ボトルネックではないことが分かりました。

 図1に、VTuneにより示された最も計算負荷の高いルーチンが示してあります。この図から、サブルーチンinterpolate_m3d(), do_advection()がホットスポットとして、並列化対象に特定されました。この他のものは、組み込み関数と同期に用いられるOpenMPランタイムです。

図1：GITMプロファイル

図1では、4つのOpenMP並列領域が示されており、これらが本レポートの中心になります。ここでのCPU時間は、16スレッド全てに関して累積された値です。

４．スケーラビリティ

 GITMのparallel scalabilityを図2に示します。ここで、左のy軸は高速化率、右はシリアル実行されているNetCDFのI/Oを含むparallel efficiencyです。

図2：GITMのスケーラビリティ

図2からわかる通り、コードのparallel scalabilityは、2から4スレッドへ急激に落ち込み、その後徐々に悪化します。16スレッドを用いた際の高速化率は約5倍、効率は31%にとどまり、効率的とは全く言えません。

５．効率

 特定された4つのOpenMP領域は、それぞれ同期時間を持ちます。同期時間は、他のスレッドが処理を完了してコード内のバリアのようなポイントに到着するまでの待ち時間です。Potential gainは、完璧な負荷バランスとオーバーヘッドが無いことが仮定されたもので、図3に示されています。これは実行時間Wall clockの削減につながります。

図3：GITMのPotential gain

この図から、OpenMP領域の最適化に余地があることが分かります。表1は、関心領域におけるGITMの同期時間を示しています。これは全16スレッドの合計です。CPI（cycles per instructions）の値は、小さければ高性能、大きければ低性能であることを示します。

CPU Total Time	801.530 seconds
CPU Time in user code	601.904 seconds
CPU Synchronisation Time	199.118 seconds
CPI Rate	2.761

表1：16スレッドのGITMの同期時間とCPI

この表から、CPU時間の25%程度が同期時間にとられていることが分かります。CPIは2.761で、VTuneは改善対象としてハイライトしました。しかしながら、CPIの算定ではベクトル命令は1命令とカウントされており（この実行では、ベクトルユニットは4つの倍精度演算を1演算とすることが出来ます）、理想的なCPI値が明確でありません。さらに、CPIを改善する方法は非常に曖昧です。

６．シリアル実行性能

 GITMはOpenMP並列モデルを採用していますが、各スレッドを更にベクトル化することで性能を上げることは可能です。本レポートでは、図3の4つのOpenMP領域のベクトル化に焦点を当てます。インテルコンパイラの最適化メッセージ/レポートから、どの演算がベクトル化されたのかを見ることが出来ます（コンパイルフラグ：-qopt-report=2 -qopt-report-phase=vecを用います。）。

 GITMコードは多くのFortran90配列演算を含んでおり、インテルコンパイラによるベクトル化が期待できます。ベクトル化は、コンパイル時の自動ベクトル化や、OpenMPのSIMD指示で可能です。ここでは自動ベクトル化について調査します。配列要素がもしメモリー境界に整列されていない場合は、ベクトル演算は「ピール」、「ループ」、「リマインダー」の順で処理されます。これらフェーズが少ないほど、すなわち、「ループ」フェーズで使用される配列境界がメモリ境界上に完全に整列されていると、演算はより高速に実行されます。

６．１　706:interpolate_m3d()のベクトル化

 記法706:interpolate_m3d()の意味は、OpenMP領域がサブルーチンinterpolate_m3d()の706行目から始まることを意味します。このセクションは、ベクトル化可能なFortran90配列演算を多く含みます。例えば以下のような文です：

uu(:,inij:endj,:) = (uu1(:,inij:endj,:)*t2+uu2(:,inij:endj,:)*t1)/ttplus

 Fortranはカラムメジャー順であるため、最初のインデックスがベクトル化されて、後ろのインデックスはループとなります。これは、インテルコンパイラにより自動で処理されるFortran90言語機能の一つです。この領域に含まれるこの他の演算も同じ構造です。インテルコンパイラのメッセージは、こうした配列演算が3つのフェーズでベクトル化されていることを示していました。

６．２　748:interpolate_m3d()のベクトル化

 この領域は、上記と同じ配列演算を持ちますが、下記のようなケースも含まれていました：

where (hun(:,inij:endj,:) > 0) u_lag(:,inij:endj,:) = &
& uu(:,inij:endj,:)/hun(:,inij:endj,:)/dx

 この配列演算は、where文が真の場合のみ実行されますが、これは実行時に決まる内容です。コンパイラは、これがどの程度の頻度で実行されるかがわからないことから、コンパイル時にベクトル化をしません。ベクトルレジスタをパッキングするためにはスキャッター/ギャザー操作が必要になるため、オーバーヘッドが存在します。コンパイラは、静的解析によるヒューリスティックスに基づいて、これをベクトル化しないことに決めました。これは、-vec-threshold0コンパイルフラグを使用すると、強制的にベクトル化することができます。詳細については、8.1.6項の推奨事項で説明します。

６．３　732:interpolate_m3d()のベクトル化

 この領域には、下記の構造の2つのFortran90配列演算が有ります：

hun(i,j,1:kmax) = 0.5*(hn(i+1,j,1:kmax)+hn(i,j,1:kmax))

この配列演算は、「ループ」と「リマインダー」フェーズでベクトル化されていました（「ピール」は除かれました）。つまり、より高速に実行されます。

６．４　273:do_advection()のベクトル化

 この領域は、ベクトル化候補であることが直ぐにはわからない様々な演算を含むブロックを含んでいます。ベクトル化には、OpenMP SIMD指示をマニュアルで追加する必要があります。

７．　負荷バランス

 706:interpolate_m3d()の負荷バランスを図4に示します。横軸は、同時に動いているコア数、縦軸は経過Wall clock時間です。全てが16コアで動く場合とシリアル部分が1コアで動作している場合が、完全な負荷バランスを持つときに示されます。

図4：706:interpolate_m3d()のOpenMP領域の負荷バランス

 図5は、748:interpolate_m3d()の負荷バランスです。

図5：706:interpolate_m3d()のOpenMP領域の負荷バランス

 図6は、732:interpolate_m3d()の負荷バランスです。

図6：706:interpolate_m3d()のOpenMP領域の負荷バランス

 図7は、273:do_advection()の負荷バランスです。

図7：706:interpolate_m3d()のOpenMP領域の負荷バランス

図4では非常に大きな負荷分散が発生しています。図5-7ではそれと比較してより小さな負荷分散です。

８．　観察結果

 このレポートでは、コード内の計算負荷の大きな部分であるOpenMP並列領域に焦点を当てました。負荷分散が大きな問題であることが判明しました。以下に今後の改善に関する提言を記載します。

８．１　短期的な提言

 このセクションでは早急な実装が可能な推奨事項を記します。

８．１．１　コンパイラ

 より多くの配列演算をベクトル化できるよう、最新バージョンのコンパイラを利用すべきです。古いバージョンは同じ最適化を導く可能性は低いです。最新のバージョンは実利用前に適切なテストが必要です。今回用いたインテルコンパイラのバージョンは16.0.2です。

８．１．２　配列アラインメント

 配列がメモリー境界に整列されてない場合はベクトル化では3フェーズが必要になります。フェーズが少なければオーバーヘッドも小さくなります。配列を出来る限りメモリー境界に割当てるには、コンパイラフラグ-align array8byteを用います。これは、配列要素を8バイト境界、すなわちFortran倍精度値に割当てます。これはインテルコンパイラのデフォルトとなっていますが、これが常に設定されるようにMakefile内に直接定義するほうが無難です。

８．１．３　浮動小数点除算

 浮動小数点の除算は乗算よりも遥かにCPUサイクルを多く消費します。これはセクション6.1と6.2で用いられています。分母の逆数を掛けることで除算と置き換えることが出来ます。

８．１．４　OpenMPの呼出し

 この推奨事項は性能ではなく、シリアル実行を行う場合のものです。OpenMPサブルーチンの呼出しは、例えば以下の様にプレフィックス!$を用いるべきです。これは、OpenMPを用いない指定で、すなわちフラグ-qopenmpを用いてコンパイルした場合、呼び出されません。

!$　　　　NTHREADS = OMP_GET_NUM_THREADS( ) !$　　　　TID = OMP_GET_THREAD_NUM( )

８．１．５　706:interpolate_m3d()について

 このOpenMP領域はFortran90配列演算を含みます。ここでは、OpenMP並列領域内で2番目の配列インデックスにスレッドIDを使用します。Fortran90配列演算は洗練され使い易いものですが、必ずしもキャッシュ効率が良いとは限りません。多くのキャッシュミスL1D.REPLACEMENT=10,476,015,714（レベル1データキャッシュ置換）が存在し、図6.1で見たような大きな負荷分散を引き起こしています。2番目と3番目の次元（ランク）の配列操作を2つの外側DOループとして、最初の次元をFortran 90配列演算として残すことが推奨されます。これはIntel compilerでベクトル化可能です。この場合、これらのDOループをOpenMPを用いて並列化することで、負荷分散を低減させるように様々なスケジューリングポリシーを適用することが可能になります：

1. schedule(static, chunk)：異なるchunkサイズを試す。
2. schedule(dynamic, chunk)：異なるchunkサイズを試す。これはより良好なスレッド間負荷バランスに効果的
3. schedule(guided, chunk)：dynamicを似ているが、作業ブロックがさらにスレッド間で可能な限り均等に割り当てられる。

 ここでは、最適なchunkサイズを持つOpenMP DOループがより良いキャッシュ効率を実現するであろうことが考えられます。最適なchunkサイズは発見的に得られるケースがほとんどです。

８．１．６　748:interpolate_m3d()について

 このOpenMP領域の負荷バランスも改善可能です。負荷分散は上記と同じ配列演算の不足から生じています。コンパイルフラグ-vec-threshold0を用いると、where文が常に真になるとしてベクトル化されます。-vec-thresholdnについて、nに対してどの値でコード性能が高いかを調べる必要があります。ここで、nは1から99の間の値を指定します。

８．１．７　732:interpolate_m3d() について

  このOpenMP領域は負荷分散がありますが、OpenMP DOループで囲まれているため、OpenMPスケジューリングを以下のようにすることで簡単に改善できます：

!$　　　　OMP PARALLEL DEFAULT(SHARED) PRIVATE(i, j) !$　　　　OMP DO SCHDULE(TYPE, CHUNK)

ここで、TYPEはセクション8.1．5でのスケジューリング型です。

８．１．８　273:do_advection()について

 このOpenMP領域の改善も可能ですが、実行人に占める割合が小さく、並列領域が大きなコードブロックを持つため大きな修正が必要になります。この領域は、OpenMP SIMD指示行による直接的なベクトル化が可能と、より細かい粒度の並列処理となるようにOpenMP領域サイズを縮小することで改善できます。

８．２　将来的な推奨事項

 将来的な性能改善のための推奨事項は以下のとおりです：

1. MPI並列化を行う：これにより、複数の計算ノードでの実行が可能になります。現状は一つのノード上でしか実行できません。大きな仮想(共有)メモリーマシンで行っても、一般的にOpenMPのスケール性を上げるのは困難です。また、こうした仮想メモリーマシンは通常の計算ノードを集めるよりも高価です。
2. 並列NetCDFライブラリを使用する：並列ファイルシステムに対してMPI-IOを用いたI/Oを行います。I/Oはボトルネックとして識別されていませんが、並列I/Oを使用することでより大きな問題に対処できるようになります。

---

補足事項：この作業はPOPプロジェクト(*)の任務として担当したNAGのスタッフが実施しました。