2019 Volume 18 Issue 4 Pages 169-175
分子動力学プログラムLAMMPSは，分散並列処理により高い効率で実装されており，スーパーコンピュータ「京」でも，大規模なノードを使ってもよくスケールし，高性能を発揮している．LAMMPSの さらなる高速化を目指し，Mod -FixLan機能で利用されている乱数ルーチンのSingle Instruction Multi Data (SIMD/ベクトル) 化及びOpenMPでのスレッド並列化によるさらなる高速化の実現を試みた．乱数生成は逐次処理のアルゴリズムが基本でありSIMD化及びOpenMPによる並列化の難しい部分の高速化を新たに実装した．特に乱数ルーチンの実装の改良によって，全体で46%程度の性能向上が観測された．まだまだ，LAMMPSには高速化の余地がある．
多くの実用プログラムがスーパーコンピュータ「京」[1] の本格利用を目的に，クラスタータイプの高性能並列計算機用に並列化されている．
材料設計シミュレーションに広く用いられているのは，Large-scale Atomic/Molecular Massively Parallel Simulator. LAMMPS[2]である．LAMMPSは，古典分子動力学(MD)計算プログラムで，金属や半導体などの固体，生体分子やポリマーなどのソフトマターやメゾスコピック系物質などに対応した多彩なポテンシャルが用意されている．
そのため，並列処理や高速数値演算装置のGraphics Processing Unit: GPU[3]やXEON Phi[4]を利用して容易に高速化が可能である．もともと，逐次計算用フリーソフトウェアであり，プログラムの把握が比較的容易である．またGNU GPLライセンス(https://ja.wikipedia.org/wiki/GNU_General_Public_License)[5] に基づいて自由に改変可能である．
高速計算に大きな影響を与える高速演算処理機構の，General-purpose computing on graphics processing units; GPUによる汎用計算:GPGPU[6]やXEON Phiは，プログラムを組んでみると性能上昇効果が容易に得られ，従来のベクトルスーパーコンピュータのように，計算機特有の最適化によって性能上昇を見る事ができる
LAMMPSもGPGPU対応ライブラリとモジュールがあり，GPGPUが実装されている計算機で高性能を発揮する．
本報告は，LAMMPSを更に高速化するために，Message Passing Interface: MPI[7]による分散一定のスケーラビリティを確認した後，さらにプログラムの高速化を目指し，single instruction multi data (SIMD)[7]化及びスレッド並列化されてない部分(本報告では乱数を生成している部分)を最適化した．これによりMPI分散並列と化及びスレッド並列化によるハイブリッド並列に高速化されプログラム全体で46%の高速化を実現した．
プログラムのSIMD化及びOpenMP化は，GPGPUを使わない場合でも，大きな効果があるということも改めて確認できた．MPIによる大規模分散並列化とSIMD化及びOpenMP化は両立する．
性能評価の条件と環境をまとめて示す．Figure 1に計算時間をまとめた．
Wall clock time of LAMMPS with a test datum (0.13*109 atoms * 10000 steps) on the K computer
測定内容は，以下の通りである．
・測定計算機:京 786ノード(+1536ノード)
・プログラム: MPI[7], OpenMP[8]を用いた並列化LAMMPS
・性能の測定方法:MPI 並列 とOpenMP並列 を組み合わせたハイブリッド並列．(「京」1ノード は8コア構成．ノード はMPI並列に割付．コアは MPI 並列 and/or OpenMP並列に割付) 786ノードでのハイブリッド並列の組み合わせは3種類で以下の通りである． 1) 786MPI x 8 OpenMP(768 × 8 =6144) 2) 1536MPI x 4 OpenMP(1536 × 4 =6144) 3) 6144MPI x 1 OpenMP(6144 × 1 =6144)
ここで，pair, bond, neigh,等は，LAMMPSの各計算機能を示す．測定値は経過(wall-clock)時間で，5回の測定で最短時間を採用した．
Figure 1に，ノード数が768の時の，どの3つの組み合わせのハイブリッド並列での実行時間は，同じコア数での計算となり，どれも同じ時間になると期待される．事実，pairやbondはどのハイブリッド並列の組み合わせでは同じ実行時間となっている．しかし，modifyだけは，3)番目のカラムに示す条件のような場合(すなわち，MPI並列の数が大きくOpenMP並列の割合が少ない場合)，実行時間が短くなっている．
この原因を調査するとmodifyの並列化は，MPI 分散並列化のみであり，LAMMPSの当該部分は，OpenMPによる並列化がされてない．このためMPI並列の割合が大きい3)の場合，性能が向上しているように見える．
Important program step 1: Step of random number generation. It looks hard to apply the SIMD technique. void FixLangevin: post_force_untemplated ... for (int i= 0; i< nlocal; i++) { .... fran[0] = gamma2*(random->uniform()-0.5); fran[1] = gamma2*(random->uniform()-0.5); fran[2] = gamma2*(random->uniform()-0.5); .... fsum[0] += fran[0]; fsum[1] += fran[1]; fsum[2] += fran[2]; ... }
modifyの実行時間はFixLangevin (温度効果を乱数で緩和する手続き)で消費されており，この手続きで，乱数が生成されるため，この処理に実行時間が集中していた．そのためFixLangevin から乱数を呼び出すオーバヘッドの削減やSIMD並列化，スレッド化等々の高速化策を試み，これらの最適化によりさらなる性能向上の実現を目指したプログラムを作成し性能を評価する事とした．
FixLangevin で乱数を呼び出す箇所をImportant program step 1: に示した．用いたテストでの計算コスト比は，乱数発生4，計算本体1である．このようにFixLangevinは ，乱数発生が処理時間を占めている．LAMMPSで用意された乱数発生のアルゴリズムは，逐次に発生させている．プログラム中，乱数の逐次発生では，最適化を行えないため，乱数を複数生成し最適化を促進する．プログラムは，nlocal個×3 の乱数を作成するrand[][]関数とし，SIMD化などの高速化技術の適用を検討する．
Important program step 2で示すように，乱数を複数同時生成することでベクトル化(SIMD化) 及びスレッド並列化(OpenMP化)の有効性を試し，FixLangevinの大幅な性能改善を実現する．
Important program step 2: Step calling random number generation. This is the most time consuming step for generation of random numbers. It looks hard to apply the SIMD technique, too.
void
FixLangevin post_force_untemplated ... for (int i = 0; i < nlocal ; i ++) .... fran[0] = gamma2*(rand[1][i] −0.5); fran[1] = gamma2*(rand[2][i] −0.5); fran[2] = gamma2*(rand[3][i] −0.5); .... fsum[0] += fran [0]; fsum[1] += fran [1]; fsum[2] += fran [2]; ... }
乱数のアルゴリズムをImportant program step 3. にプログラムで示す．本プログラムを元に，乱数生成でSIMD化と並列化が可能性を検証する．
Important program step 3: Program for confirmation of speed up by parallel generation of random numbers double RanMars::uniform() { double uni = u[i97] − u[j97]; if (uni < 0.0) uni += 1.0; u[i97] =uni; i97–; if (i97 == 0) i97 = 97; j97–; if (j97 == 0) j97 = 97; c -= cd; if (c < 0.0) c += cm; uni -=c; if (uni < 0.0) uni +=1.0; return uni
このプログラムは，乱数を逐次に生成しuniに格納して呼び出し元に戻る．配列uは97個の乱数のseedである．
i97, j97はポインタである．再び乱数を生成する時は，前ステップで計算されたseed uを用いて，u, i97, j97が逐次に更新される．
ただし，この関数のプログラムは，SIMD化及びOpenMP化を意図したコードになっていない．このため乱数をn個同時発生するように変更し，実行効率の良いSIMD化及びOpenMP化が可能な部分がないか検討する．
Important program step 4でC++からfortranにコードを変更し，乱数を複数生成するようにプログラムを修正した手続きを示す(実際のプログラムでは，異なる言語は使わないほうが良い．アドレスの計算方法や演算順序が異なる事が多く，プログラムのデバグに労力がかかることがある)．
本プログラムでは，乱数をn個生成し配列 uni に格納後，呼び出し元に戻る．
Important program step 4: Confirmation of balky generation of random numbers. Improve index reference of balky generation of random numbers are also shown i97 and j97 at after “enddo” do i=1,n uni(i) = u(i97) − u(j97) if (uni(i) < 0.0) uni(i) = uni(i) + 1.0 u(i97) = uni(i) i97=i97-1 if (i97 == 0) i97 = 97 j97=j97-1 if (j97 == 0) j97 = 97 c = c-cd if (c < 0.0) c = c+cm uni(i) = uni(i)-c if (uni(i) < 0.0) uni(i) = uni(i) + 1.0 enddo
まずは，cの依存性を解消する．cはループでcd回引かれるため，c-cd*iとすることで，ループ更新によるcの依存性は解消される(Important program step 5)．この結果，cの更新はOpenMPにより並列化できる．(Important program step 6)．
Important program step 5: generation of random numbers. do i=1,n uni(i) = u(i97) − u(j97) if (uni(i) < 0.0) uni(i) = uni(i) + 1.0 u(i97) = uni(i) i97=i97-1 if (i97 == 0) i97 = 97 j97=j97-1 if (j97 == 0) j97 = 97 ctmp = c-cd*i if(ctmp < 0.0) ctmp = ctmp & + (int(-ctmp/cm)+1)*cm uni(i) = u(ii)-ctmp if (uni(i) < 0.0) uni(i) = uni(i) + 1.0 enddo c = c-cd*n if(c < 0.0) c = c + (int(-c/cm)+1)*cm
Important program step 6: Thread Parallelization by loop division of “c” do i=1,n uni(i) = u(i97) − u(j97) if (uni(i) < 0.0) uni(i) = uni(i) + 1.0 u(i97) = uni(i) i97=i97-1 if (i97 == 0) i97 = 97 j97=j97-1 if (j97 == 0) j97 = 97 enddo Comp parallel do private(i,c_tmp) do i=1,n ctmp = c-cd*i if(ctmp < 0.0) ctmp = ctmp + & (int(-ctmp/cm)+1)*cm uni(i) = uni(i)-ctmp if (uni(i) < 0.0) uni(i) = uni(i) + 1.0 enddo
Important program step 7: Program for consuming index generation do i=1,n ii=mod(i97-i+1,97) if(ii<=0) ii=ii+97 jj=mod(j97-i+1,97) if(jj<=0) jj=jj+97 u(ii) = u(ii) − u(jj) if (u(ii) < 0.0) u(ii) = u(ii) + 1.0 uni(i)=u(ii) enddo i97=mod(i97-n,97) if(i97<=0) i97=i97+97 j97=mod(j97-n,97) if(j97<=0) j97=j97+97
Important program step 8: Confirmation of generation of index i, ii, jj.
i ii jj
1 97 33
2 96 32
…
33 65 1
34 64 97
35 63 96
…
97 1 34
98 97 33
99 96 32
次に配列uの更新について検討する．Important program step 7にある配列uの更新のインデックスi,ii,jjを確認すると，Important program step 8から配列uのインデックスは97で折返されていることが分かる．ここで，インデックスを折り返さず，十分長い配列uを使い，配列uの更新において，iiをi-97とし， jjをi-33とする．このuの更新は等価である．すなわち，この乱数生成のアルゴリズムは，Lagged Fibonacciである．n個の配列を作るのに，前後97個は，別な処理が必要であるが単純なループとなる．これをプログラムで書くと次のようになる． do i =,,, u(i)= u(i-97) − u(i-33) enddo 配列uから長い配列uniに97個セットした後，計算を続行するプログラムをImportant program step 9に示した．なお，33以上離れるとuniの要素間に依存関係があるためスレッド並列化は不可である．しかし，SIMDの処理要素数(ベクトル長)が33以下ならSIMD化は可能である．
Important program step 9: Confirmation of efficiency of random number generation uni(i),u(...)=u(...),... i=1,97 do i=98, n-97 tmp=uni(i-97)-uni(i-33) if (tmp < 0.0) tmp = tmp + 1.0 uni(i) = tmp enddo u(...),uni(i)=uni(...),... i=n-96,n
LAMMPSで使われていた乱数は，M系列と同じで，Lagged Fibonacci のアルゴリズムである．ラグは，97と33の組み合わせで，周期は297である．LAMMPSの当該部分がMPIで並列化されていたのは，乱数の初期値をMPIプロセス番号から誘導しているためである．厳密に言うと相関がない初期のseedとは言えない．それでもFixLangevinで生成された乱数を使う事では問題ない．逆に，そのような乱数の初期値を選んだことで乱数がMPIで並列化できているように見える．
乱数の性能を評価する．性能測定条件は，12万原子の時の乱数発生(11,616*3個を1,000回)をシミュレートするテストプログラムを「京」にて実行する．プログラムはFortranコードからC++言語に書き戻し評価する．データの収集と解析は「京」で利用できるfipp/fappo (富士通製性能評価ツール[9]を使用する．
Table 1 に結果を示す．乱数発生は，上から1) 1個ずつ生成(オリジナル)，2)複数生成，3) OpenMP化，4) SIMD化と，順次実装した．結果的に3)と4)のスレッド並列化及びベクトル化の効果である．細かな点は，性能解析ツールを使ってホットスポットを特定でき，性能のボトルネックであることを評価し，当該部分をOpenMP化やSIMD化することでボトルネックが解消することを確認する．これらの最適化により乱数生成の性能向上を確認できる．
| Elapsed | MFLOPS | MFLOPS/Peak | |
| (sec.) | (%) | ||
| 1.249 | 207.9 | 0.162 | sequential gen. |
| 1.032 | 251.7 | 0.196 | parallel gen. |
| 0.658 | 926.6 | 0.723 | OpenMP |
| 0.104 | 6329.0 | 4.944 | SIMD |
最終ページに掲載のFigure 2 に性能評価の結果を示した．測定条件は，2章に記したものと同じである．
Confirmation of improvement of LAMMPS's efficiency by the random number generation.
Figure 2の左側は最適化前で，MPIとOpenMPによるハイブリッド並列という観点から見て，乱数生成の並列化は不十分であった．Sandia National Laboratoriesの開発者はハイブリッド並列のことは考慮していない．Figure 2の右側に示したように，ハイブリッド並列化による並列処理の最適化後では，同数のコアで分散並列化及びスレッド並列化のどの組み合わせでも，ほぼ同等の実行時間となっている．
左右の結果を見ると，ハイブリッド並列化とSIMD化の実現によりmodifyの性能向上がなされたことが分かる．また，右の図で768ノードと1536ノードで良くスケールしていることが分かり，MPIとOpenMPのハイブリッド並列化が実装されていてそれらがバランスよく働いている事を示す結果となっている．
全実行時間で，再左端のオリジナルの並列化プログラムと最右端の比較では，ハイブリッド並列化と，乱数生成のSIMD化による高速計算の実現で，並列不能部分の並列化により46%の性能向上となった．
本研究で開発したLAMMPSの高速化アルゴリズムは，2019年8月で運用が終わる「京」に続く「富岳」でも，高性能を発揮する．富岳は基本的なチップ構成が旧京と同じであるが，メモリアクセスの性能が良いので，より一層の力を発揮する事が期待できる． 本研究結果はすでにLAMMPSの管理機関であるSandia National Laboratoriesに報告済みである．本研究の成果が公開版LAMMPSに組み込まれるかどうかは，論文執筆時点(2019/05/01)では，不明である．
本研究で開発・実装した並列乱数生成アルゴリズムに関し，片桐孝弘教授(名古屋大)，青山智夫研究員(江戸川大)，小柳義夫サイエンスアドバイザー(RIST)の励ましに感謝する．
愛媛大学教育学部の原本博史准教授，理学部数学科の土屋卓也教授，山崎義徳准教授達との深い議論を持てた．感謝する．
本研究は愛媛大学理学部(長岡伸一教授)とRIST (小久保達信次長)の共同研究契約による研究の一部である．
