＊ここに掲載するのは、STFC Daresbury研究所のLaurence J. Ellison博士によるHECToRレポート「A load balancing algorithm for Molecular Dynamics simulations, Laurence J. Ellison, 10th May 2011」を要約したものです。
概要
領域分割（DD）法を用いて並列化された分子動力学（MD）シミュレーションにおいては、シミュレーション領域（MDセル）は同サイズの空間領域に分割され、シミュレーションで用いるマシン上の特定のコアに配置されます。各コアの計算負荷は次の2種類です：(i)領域内の原子に働く力の計算、(ii)力の計算に必要となる原子座標を隣接コアから転送するメッセージ通信。
MDセル内の原子の分布が一様であれば、各コアは大まかにいえば同じ計算負荷を受け持ちます：計算負荷がバランスしている、という言い方をします。もし原子の分布が大きく偏っている場合は、ある領域は原子を多く含み、その他は比較的疎らに原子が分布します。この結果、高密度な原子を持つコアは多くの計算負荷を担い、他方疎らな密度の原子を担うコアはアイドルする時間が長くなります。ここには大きな欠点・不利益が2つあります。第一に、MDシミュレーションは最も実行に時間のかかるコアで速度が決まるため、平均ステップ時間、つまり解を得る時間は最も負荷の高いコアで決まります。第二に、多くのコアが少量の作業しか行わないため、ハードウェア効率が悪化します。
高負荷なコアから低負荷なコアへの原子の再配置により、この状況を改善する負荷バランスプログラムのプロトタイプを開発しました。このために、各領域を複数の同サイズのセルに分割しました。繰り返し、MD時間ステップ内で、各コアは自身の計算負荷を現状の原子配置の関数として推定します。この計算負荷は集団通信を通じて蓄積され、計算負荷が分散した時点でセルは高負荷コアからより低負荷なコアへ転送されます。セルの転送により、他のコアに配置されていた対象セル内の原子が、この時点で対象セル内の原子に働く力の計算対象になることを意味します。セルの転送は、指定の許容範囲内で負荷がバランスするまで続けられます。こうしてシミュレーションの高速化が達成されるはずです。即ち、負荷バランスを行わない場合に比較して、平均時間ステップ長が削減されるでしょう。負荷分散がうまくいけば、より短い時間で結果を得られ、ハードウェアをより効率的に用いることが可能になります。
当初の目的は、(i)並列MDシミュレーションの性能改善のための負荷バランサープロトタイプの開発、(ii)負荷バランサーの徹底的なテストと性能に影響するパラメータの理解、(iii)並列領域分割を用いるが負荷バランス機能を持ち合わせていない最新のDL_POLY4への負荷バランサーの組み込み。最終的にロバストかつ柔軟なアルゴリズムを開発し、特にセルと領域の次元を原子間相互作用のカットオフ距離と完全に切り離して選択可能にしました。コードは誤差に関して多くの試験が行われ、性能テストの結果は良好でした。多くのテストの中で、負荷バランサーは理論的な最高値に近い高速化が達成されましたが、極めて稀にコントロールテストよりも遅い場合がありました。
言い換えれば、ほとんどの場合には負荷バランサー開発で失うものは僅かで、効率で得をする場合のほうが多いでしょう。残念ながら時間の関係で、DL_POLY4への組み込みは達成しませんでした。これらのコードは共に独自の複雑な設計がされており、一方を他方に組み込むのは容易ではありません。さらに注意深く行っていた負荷バランサー開発の後半に予期せぬ設計上の問題が生じました。この問題解決には、早期にDL_POLY4へ組み込みを急ぐよりも、プロトタイプのバグを取り、可能な限り効率的に、かつ確実なものにするほうが先決でした。
現状の計算手法は大規模原子集団の運動力学、あるいはマクロな流体の振る舞いのいずれかに独立してフォーカスされています。最小単位では分子動力学（MD)シミュレーションが、閉じた領域内での個々の相互作用に従う個別の原子の運動を追跡します。この手法は液体の輸送やせん断性の決定に成功を収めてきました[5,6]。他方、直接数値シミュレーション（DNS)は連続体の仮定を採用し、大きな計算領域内の各点における流速に対する非線形偏微分方程式を解きます。例えば、ZakiらによるＤNSコード Transflow[13,14,15]は、乱流における渦モード相互作用のような純粋な流体現象に用いられてきました。しかしながら現状では、これら方法論を一つのフレームワークへ統合する、特に大規模並列科学計算の観点で幾つかの試みが成されています。
このプロジェクトの目標は、様々な物理現象の効率的な連成シミュレーションの最先端の計算フレームワークを開発、検証、最適化することです。特にこのプロジェクトでは、現実に近いシミュレーションのために高忠実度ＤNSをＭDシミュレーションと連成することに焦点を当てています。このＭD-DNS連成コードは、ミクロな表面やコーティングのような分子スケールの現象がどのようにマクロな流体運動に影響を及ぼすか、という点に関して新たな能力と洞察を提供します。この全スケールを解像する能力を用いることにより、簡略化され、時に不正確な分子スケールの構成的なモデルを不要となります。
コード概要
負荷バランサーのソースコードはFortran90で記述され、メッセージ通信にはMPIを用います。F90標準に準拠しているか、また不適切なプログラミングを除去するために、ツールForcheck[softeng.cse.rl.ac.uk]を用いて、シリアル版のソースコードを検証しました。トップレベルのプログラムコードは3つのステージで構成されています。
ステージ1では、多くのシミュレーションコードに共通の、静的配列の割り当てや、定数の定義、入力ファイルからのユーザ定義指示の読み込み、初期原子座標の設定などが実行されます。
ステージ3では、主に診断データの出力が実行されます。診断には、全セルと原子数のチェックや、主要サブルーチンのローカルおよびグローバル平均実行時間（作業負荷のボトルネックや不均衡の推定に便利です）、プログラム内で使用された主要な動的配列に割り当てられたメモリー量などが含まれます。
コードの主要な部分はステージ2に存在します。これは通信と計算の多くを実行する時間ステップのループです。
·初期の負荷分散決定のためのセル分布計算
以下の式から各コア上の相対的作業負荷Wを推定します：
W = W1 + ρW2
ここでW1はこのコアで評価すべき2体相互作用数で、W2はこの2体相互作用を完結するのに必要となる、他のコアから転送すべき原子座標数です。パラメータρは、1つの相互作用計算の平均コストに対する1つの原子の輸送にかかる平均コストです。負荷バランスサブルーチンの最初にWは各コアで独立に計算され、その結果が蓄積された後に全作業の分布が求まります。
ここで作業負荷の計算を実行するためには、ローカルコアは必要となるコアの分布を知る必要があります。ここで、セルAとBにそれぞれ原子iとjがカットオフ距離以内に置かれていると考えましょう。
iがjに及ぼす力は、その反対のjがiに及ぼす力と同じです。よって一つの力の成分fijのみを計算して、jにかかる力の計算へはその逆符号fji=-fijを加えればよいことになります。このようにして、与えられたセル内の原子の力の計算には、カットオフ距離内のすべてのセルに在る原子座標すべてが必要なのでなく、その内の“required cell”と呼ぶ「部分的なハロ」セルにあるもののみを計算すれば可能です。
あるコアに割り当てられた任意のセルに対する“required cell”は、それと異なるコアに割り当てられている場合があります、これを以降“オフコア”と呼びます。“required cell”内の原子座標がオフコアである場合は、転送が必要になります。負荷バランス処理ではコアからコアへセルを移動させるため、特定のセルがどのコアにあるかの検出と、その後のセル内の原子座標の検出には特別な困難さがあります。この問題の解決方法は、セルが物理的に配置される領域を担当するコア（“ホームコア”と呼びます）を、コアの領域上に物理的に配置されるセル（ネイティブセルと呼びます）の情報と共に周期的に更新することです。ホームコアに格納すべき情報は、そのネイティブセルの現在のホスト（セルが現在配置されているコア）とその分布です。オフコアのセル情報を必要とする如何なるコアも、そのセルが何処に配置されていようとも、その必要とするセルのホームコアにアクセスします。この“forwarding address”パラダイムは、任意のセルはシステムの如何なるコアにも配置は可能ですが、ユニークなグローバルインデックスで指定されるセルや、その物理的配置、つまりそのホームコアは常に決定可能であるという事実に基づきます。これは必要とするセルの分布を求める際に用いされます。このやり方の注意すべき点は、情報の転送が1対1メッセージ通信で行われることです。この種の効率は、負荷バランス処理において、データをコア外に移動させる際の通信オーバーヘッドが目立たない場合に効果があります。同様にスケーラビリティを維持するために、集団通信は最小限にする必要があります。
·負荷バランス処理
“required cell”の分布が得られたら、各コアは作業負荷の計算を行います。全コアの作業負荷はグローバル総和で蓄積され、各コアは作業負荷分布の平均と標準偏差を計算します。標準偏差と平均の比が望ましい閾値以上であれば負荷バランス処理が実行されます。負荷バランス処理の最初のステップは、作業負荷の降順にコアをランク付けすることです。次に最高負荷コアと最低負荷コアをペアリングします。この関係の中で高負荷コアは一つのセルを低負荷コアへ送ります。各コアは自身の負荷を再評価し、コアのペアはこの情報を共有します。コア間の負荷の差と平均の比が指定された閾値よりも大きければ、さらにセルを転送します。この手順を負荷バランス処理ループで要求されたバランスに達するまで繰り返します。さらに別のペアについて負荷バランス処理を行っていけば、作業負荷の分布は全体としてバランスしていきます。
·力の計算に必要となるセルの決定
力の計算ループの前に、各コアは全ての相互作用計算に必要となるセルデータを転送しなくてはなりません。通常の負荷バランスを考慮しない領域分割ベースのMDシミュレーションではコアが、自身の領域に隣接した領域のコアからハロデータを受け取るように実行されます。負荷バランス処理が行われる際は、このままではセルをそのホームコア外から転送できません。コアは、所望の原子がハロデータに存在すると仮定できなくなり、要求されるセルが現在どのコアにあるのかを見つけてそのセルの原子座標を要求する必要があります。この手法の実装はセル分布を求めるやり方と似たものです。この処理の中でホームコアは、ホストコアがセルに対してどのくらいの量の要求を受信するかを見積もります。コアが受信するメッセージ数の予測はプログラム設計の最大の課題です。
·ローカルでの力の計算に必要な原子の転送
ここでは、要求元コアは必要なセルをどこへ要求するか、また、要求先コアは必要なセルを要求する他のコアからどのくらいの要求量があるかを知っています。要求された原子座標は転送されて、力の計算ループが処理されます。
· 原子に働く総力の計算
力の計算ループ内では、12-6レナード-ジョーンズポテンシャル形式の2体相互作用のみが計算されます。明らかに、DL_POLY_4のような分子間および分子内相互作用を扱う多目的MDコードでは機能不足です。ここでは負荷分散アルゴリズムをスクラッチから開発することを主眼としてしるため、簡潔な相互作用のみを扱っています。よってこのプロトタイプは、(i)アルゴリズムが技術的にロバストで、(ii)より現実的なMDシミュレーションの効率を改善することを実証するために十分なできる限りの機能のみを実装しています。基本的なレベルの力の計算ループですが、MDシミュレーションで扱われるコアの計算負荷レベルを十分に表現しています。
· オフコアの原子に働く逆符号の力成分の送信
力の計算では、原子に働く力を計算後に、転送された座標の原子上に働く逆符号の力の成分の計算が必要で、この結果は元のコアへ送り返されます。こうして系の各原子に働く総力が評価されます。
· 原子の移動
MDシミュレーションでは一般に、力の計算ループ後に運動方程式を解いて次の原子位置を計算します。ここではこの計算はせず、テストベッドとして、ランダムな方向にランダムな移動をさせるか、指定された方向に決まった距離移動させるかのどちらかを行います。
· セル境界を越えて新たなセルへ原子を移動する
通常の領域分割ベースのMDシミュレーションでの原子移動は、領域境界を跨る際のみ原子は異なるコアに割り当てられます。そのコアはその対象となる領域を担当するコアです。しかしながら、もしコアがそのホームコア外へ転送される場合、負荷バランス処理が実装されている場合は、セル境界を跨ぐ原子は再割り当て対象の候補としなくてはなりません。これを考慮するために、再度'forwarding address'パラダイムの方法を採用しました。ここで、ホームコアはそのネイティブセルの境界を跨る原子の座標を要求します。
·実測値からのρ計測
ρは、一つの原子転送の平均コストと一つの相互作用計算の平均コストの比です。これは、これらを実行するサブルーチンの実測値をベースにしています。
性能テスト結果
負荷バランサーの性能には様々なパラメータが影響します。一方、MDシステムは、原子密度の不均一の程度や、原子総数、原子間相互作用の計算内容、相互作用カットオフ距離などの物理的特性を持ちます。また、ハードウェア条件としては、プロセッサ数や、プロセッサの特性として、クロック周波数や、バンド幅、レイテンシーやプロセッサ当たりのコア数などがあります。さらにコンパイラやそのフラグも重要な要素です。また、セル分割や、ρの値、負荷のバランスをどこまで求めるかなどの負荷バランサ自身のパラメータも重要な指標です。通常のMDシミュレーションではカットオフ値などのようにあらかじめ決まったものもありますが、テストのためには任意に選択するものとします。これはパラメータ空間を広く扱い、膨大な仕事になるため、ここでは数個の簡潔な原子配置のみを用い、コア数を変えて、異なるカットオフ値をそれぞれで用いることとします。
全ての試験で、負荷バランス有り/無しの場合の平均速度比を計測します。テストは、2011年3月14-18日の間、HECToRフェーズ2b（XE6）を用いました。XE6はプロセッサ当たり24コアを搭載しています。使用したコンパイラはPGIで、”-O3 -fastsse”フラグを用いました。ρは常に25.0としました。負荷バランス条件：平均負荷に対する負荷バランスサブルーチン内の2つの負荷の差の比、の閾値は0.05としました。全てのテストは100時間ステップ計算しました。
·システム1
最初のテストは、立方体のMDセルの8分割分に原子が配置されて残りが真空のものが格子上に並べられた例です。これは人工的ですが、結果の解釈が容易なものです。MDセルの体積Vはコア数Pに比例して配置してテストを実行します。ここで、領域分割には、4x4x4=64か5x5x5=125分割を選択するものとします。セルメッシュが荒過ぎると、負荷バランサーによるコアからコアへの原子データ移動量が大きくなり過ぎます。逆に細か過ぎるとセルの再割り当てが多くなり過ぎ、副作用として原子データは多くの小さなパケットに分割されて、マシン全体に細かく分散されて通信オーバーヘッドが大きくなり、性能が劣化します。
レナード―ジョーンズ相互作用の様々なカットオフ値rcutに対しても同じテストを行いました。最少カットオフは0.449auで、原子間距離0.2auの2.5倍弱です。この選択理由は、通常のMDシミュレーションではレナード―ジョーンズポテンシャルが2.5σで打ち切られるからです。平衡距離21/6σに対応するポテンシャル極小は、ここでは格子間隔と等しくしています。しかしながら金属相互作用や静電相互作用においてはより大きなカットオフ値が用いられるため、大きなカットオフ値を用いたテストも行いました。その場合には原子の移動は行わず、与えられた領域分割とカットオフにおいて、負荷バランサがどのように不均衡な原子配置を処理するかを調べました。こうした静的な配置を用いると、テストの結果解釈が容易になります。
ほぼ全てのシステム1テストにおいて、6.0から8.0の高速化を観測しました。8.0倍はMDセルの1/8に原子が配置された際の理論的な最高性能です。相対的に加速が低い例(1.68倍)は、最大コア数と最大カットオフを用いた実行の場合です。この理由は、占有された領域が分割されたセルの総数が最大になり、多くのセルが他のコアに割り当てられて、原子データは分散されているためです。これは大きなカットオフと組み合わされて、要求セルの受信と力の送信に大きな通信オーバーヘッドが発生します。力の計算ループ前後の、原子データ受信と力の送信での通信負荷が大きくなって性能が劣化すると疑うことができます。
この推測は、同数のコアでプロセッサコアを飽和(24コア)させずに12コアを用いた場合に高速化率が1.68倍から7.12倍へ上昇することからも確かと考えられます。プロセッサ当たりの利用コア数の減少はバンド幅を拡大させ、通信コストを削減します。ここで、すべてのカットオフに対して、高速化率はP=216でピークに達しますが、コアを増やしていくとゆっくりと減少していきます。これも上記と同様な理由によるものと考えられます。
·システム2
このテストでは、原子を配置した一つの領域のみで構成され、原子間距離は同じく0.2auで、静的に配置されます。領域は一つのため、P（コア数）とそれによる分割領域数の増加に伴い、原子を含むMDセルの比率は比例して小さくなります。領域は以前と同様に125セルに分割されます。
高速化率のピークはP=125かP=216にあり、その後はコア数の増加と共に劣化していきます。
·システム3
システム1と同じ設定を用いますが、今回は系の物理サイズを固定して、コア数を増加させます。領域当たりのセル分割数は125で、静的な配置とします。
この設定ではコア数の増加に伴い、コアやセルあたりの原子数は減少します。これは負荷バランスを仮定せずともハードウェアを浪費する状態です。コア数を増やすとコア当たりの計算負荷は減りますが、同時に通信オーバーヘッドは増加し、カットオフが効いてきます。これは負荷バランスの有無にかかわらず生じ、平均ステップ時間は10以上で大きな減少をしなくなります。
·システム4
最後に現実的なMDシミュレーションにより近いテストを行いました。この系は半径1.25auの球領域に高密度に原子を配置し、これを4.0auの立方体MDセルに入れて、領域サイズ1.0auの4x4x4=64個の領域に分割しました。よってテストは64コアで行いました。この球とMDセルの体積比は1/8で、システム1や3と同じ比率です。セルはさらに4x4x4=64に分割されました。
球1として、原子間距離を0.1auとして球内に8144原子を配置しました。これは64コアで実行するのに比較的小さな数字です。カットオフ距離は0.249auとしました。次に球2として、原子間距離を半分の0.05au、球内原子数を72327と設定しました。これはさらに多くの計算負荷がかかるケースです。この際はカットオフ距離を0.1249auとしました。この他、下記の様々な条件を課しました：
1.原子位置を固定しました。
2.球外は真空にしたまま、時間ステップ当たり0.01au距離をX方向に原子を移動させました。これにより、系の原子が大挙してセルおよび領域境界を横切る際に、負荷バランサがどのように対処するかについて見ることができます。
3.球を固定して球外は真空にしたまま、球内の原子をランダムに動かします。ただし、時間ステップ当たりの最大移動量は球1では0.10au、球2では0.05auとします。
4.球内の原子密度の約1%の密度の原子を、ランダムに球外へ配置します。これにより、ガスが蒸気で囲まれた液滴へ凝集するシミュレーションの単純な模擬を表現します。
5.上記の場合に、すべての原子を項目3と同様にランダムに移動させます。
6.上記の場合に、原子分布が一様になるまで長時間にわたり時間ステップを実行させます。
テスト結果から、球1ケースのほうが球2ケースよりも高速化率が低いことが判りました。これはコア当たりの原子数が球1のほうが少ないことによるものです。負荷バランサはコア当たりの計算負荷が多いほうが効果的です。また、原子の移動も希薄な原子ガスの存在も、全ての場合において特に高速化の効果がないことが示されました。よって負荷バランスアルゴリズムが実シミュレーションにおいても効果があることが予想されます。
ケース6では負荷バランスの有用性が最もよく表されます。不均一な原子分布が一様になったシミュレーションの最後では、負荷バランスを行った場合のステップ時間のほうが僅かに大きいです。初期の不均一な分布では大きな効果がありますが、一様分布になった場合でもバランスさせる処理が含まれるためです。しかしながら差は僅かで、一般的にいえば負荷バランスアルゴリズムが有効であるといえます。
性能テスト結果
ここで開発された負荷バランスアルゴリズムには、処理の初期化とリファインを止めるタイミングを決める作業負荷の推定など、いくつかの欠点があります。より最適な計算負荷の分散のためには、より正確な作業負荷推定が必要です。そこでは、オフコアへ転送されるセルがどれか、またそれが割り当てられる新しいホストコアがどれなのかを選択するアルゴリズムも改善されるでしょう。その最大の利得は、データ分散の最小化と、負荷分散処理による通信オーバーヘッドの削減です。
それでも性能評価の結果は、このコードが良好な結果を生成することを示しています。比較的高計算負荷の系では、負荷バランス処理は良好な高速化を示しました。稀に劣化する場合がありますが大きな減速はありませんでした。
並列計算機においては、コアのクロック周波数はこれ以上大きく増えることはなく、プロセッサー当たりのコア数の増大が見込まれるため、コア間の通信に改善の努力が必要となる傾向にあります。クロック周波数が向上しない中では、コア間の作業負荷が分散する場合に、その負荷を再配置することが重要になります。一方、メッセージ通信がより有効になるに従い、データ再配置の負荷を小さくせねばなりません。こうした開発は、ここで開発された負荷バランスアルゴリズムを有利に導きます。
いずれにせよMDプログラムの設計や、負荷バランス処理を用いるかどうかの如何に関わらず、系列計算リソースの効率的な利用は多くの部分がユーザに依存します。分子系はそれ自身の特性を持ち、これらは科学的な見地のみならず実用性からも、シミュレーションの最適なパラメータ選択といった考察が必要です。
謝辞
このプロジェクトは、NAG Ltd.が運営するHECToRの分散計算科学および工学(CSE)サービスの基に実行されました。英国の国立スーパーコンピューティング・サービスである、HECToR：英国リサーチ・カウンシル・ハイエンド計算サービスは、リサーチ・カウンシルを代行するEPSRCが管理しています。そのミッションは英国学術界の科学および工学の研究支援です。HECToRスーパーコンピューターは、UoE HPCx Ltd.およびNAG Ltd.のCSEサポートサービスにより管理運営されています。