概要
以前WSL2上でCUDA環境を構築し，テストプログラムを実行する手続きについて解説した記事を書きました。本稿では，この時作った環境においてGPUで動作するプログラムを実際にビルドし，実行してみる試みを紹介したいと思います。対象のプログラムは，古典分子動力学シミュレーターLAMMPSです。おおよそ以下のような内容を予定しています。
- コンパイルする方法
- GPUインターフェースライブラリーのビルド
- LAMMPSロードモジュールのビルド
- 使い方
- GPU上で走らせる方法
- 複数GPUで実行する方法
- テスト計算
- 用いるポテンシャル EAMポテンシャル，TIP4Pポテンシャル
- 環境 CPU : Intel core i7-12700 2.1 GHz, GPU : GeForce GTX 1650 super
- むすび
用いるGPUはGPGPU専用製品と比較すると安価なグラフィックス描画用のものですが，このような環境で実アプリケーションにおいてどの程度の性能が見込めるのか確認していきたいと思います。
コンパイルする方法
GPUインターフェースライブラリー
LAMMPSインストールディレクトリーの下の lib/gpu にGPUのインターフェースライブラリーのソースコードが存在するので，これをビルドしていきます。
用いるGPUは倍精度の性能が単精度の1/32と低いため，GPUにおいてはできれば単精度で計算を実行したいです。そこで，「混合精度版」のインターフェースをビルドします。この版では，GPUの処理は単精度，CPUでの処理は倍精度で行います。最終的に重要なエネルギーや原子間力は倍精度で表すため，単精度を用いつつ計算精度の劣化を比較的小さくすることができます。
デフォルトが「混合精度版」なので，その実現のために特別なことを行う必要はありません。Makefile.cuda_mps において以下の行を編集することによって，このふるまいを制御することができます。
CUDA_PRECISION = -D_SINGLE_DOUBLE
_SINGLE_DOUBLE のかわりに _SINGLE_SINGLE とすればすべて単精度， _DOUBLE_DOUBLE とすればすべて倍精度で計算を行うようになります。
以前作った環境において以下のコマンドを実行し，NVIDIA HPC SDKをMPI込みで有効にします。
module load nvhpc-hpcx-cuda12
Loading nvhpc-hpcx-cuda12/23.9
Loading requirement: hpcx
用いるコンパイル環境はNVIDIA HPC SDKなので，それに対応した Makefile.cuda_mps というMakefileを用います。すなわち，以下の要領でビルドします。
cd LAMMPS_INSTALL_DIR/lib/gpu
make -f Makefile.cuda_mps
nvcc -I/usr/local/cuda/include -DUNIX -O3 --use_fast_math -DLAMMPS_SMALLBIG -Xcompiler -fPIC -gencode arch=compute_50,code=[sm_50,compute_50] -gencode arch=compute_52,code=[sm_52,compute_52] -gencode arch=compute_60,code=[sm_60,compute_60] -gencode arch=compute_61,code=[sm_61,compute_61] -gencode arch=compute_70,code=[sm_70,compute_70] -gencode arch=compute_75,code=[sm_75,compute_75] -gencode arch=compute_80,code=[sm_80,compute_80] -gencode arch=compute_86,code=[sm_86,compute_86] -D_SINGLE_DOUBLE --fatbin -DNV_KERNEL -o .//amoeba.cubin .//lal_amoeba.cu
...
...
このようなコマンドによってlibgpu.a というライブラリーが生成されれれば成功です。
LAMMPS本体のビルド
LAMMPS本体をコンパイルし，libgpu.a にリンクしてロードモジュールを作成します。そのためパッケージ gpu を以下の要領で有効にします。その他，テスト計算で用いる予定のポテンシャルのパッケージを有効にしています。作業はLAMMPSインストールディレクトリーの下のsrcディレクトリーにおいて行います。
cd ../../src
make yes-gpu
make yes-manybody
make yes-kspace
make yes-rigid
make mpi
Gathering installed package information (may take a little while)
...
...
ロードモジュールlmp_mpi が生成されれば無事コンパイルできたことになります。
使い方
コンパイルしたロードモジュールをGPU上で実行するには二つの方法があります。
- ポテンシャル指定の最後に
/gpuをつける。例：lj/cut/gpu -
-sf gpuをつけて実行する。例：mpirun -np 2 lmp_mpi -in in -sf gpu；ここで-in inはインプットファイル指定
複数のGPUが搭載されているようなシステムの場合，以下の要領で複数のGPU上で動作させることも可能です。
mpirun -np 2 lmp_mpi -in in -sf gpu -pk gpu 2
-pk gpu 2によって2基のGPUを用いる設定になります。多くの場合MPIプロセス数とGPU数を一致させた方が良好なパフォーマンスを得ることができると思いますが，複数のGPUを搭載したマシンを使っている場合は試行錯誤して最良の組み合わせを見出すのがよいのではないかと思います。
テスト計算
以下のような問題を用いてテスト計算を実行してみました。
- EAMポテンシャルを用いたbcc-Feスーパーセル ($100 \times 100\times 100$) , 200万原子の20ステップのMD計算
- TIP4Pポテンシャルを用いた2484分子の1000ステップのMD計算
計算負荷の観点からは，二つの例題は長距離相互作用が含まれるかどうかという点が最も大きな違いです。すなわち，EAMポテンシャルは短距離相互作用のみなのに対し，TIP4Pポテンシャルは長距離相互作用であるクーロン相互作用を含みます。クーロン相互作用はPPPM法によって評価します。なお，LAMMPSのマニュアルによると汎用のPPPM法は一部の演算をGPUにおいて計算することが可能なようなのですが，TIP4Pで用いる専用PPPM法 (pppm/tip4p) はGPU非対応なのでCPUのみで動作します。
動作環境は以下の通りです。
- CPU : Intel core i7-12700 2.1 GHz
- GPU : GeForce GTX 1650 super
用いたCPUには8つのperformance coreが搭載されています。そこで，8並列までは性能向上が期待できると考えMPIプロセス数は非並列と8並列にしました。GPUは1基しか搭載されてないので，MPI非並列実行のみです。
計算に要した時間は以下の通り（単位は秒）。
| EAM bccFe 200万原子 | TIP4P water 2484分子 | |
|---|---|---|
| GPU | 1.9 | 7.0 |
| CPU 非並列 | 23.2 | 37.7 |
| CPU 8並列 | 6.4 | 9.9 |
表の数値から明らかなように，GPU上でCPUよりも高速に動作することが確認できました。
- EAMポテンシャルの場合は非並列と比較して20倍以上，8並列と比較しても3倍以上の高速化が達成できています
- TIP4Pポテンシャルの場合差はそこまで大きくなく，非並列の場合は5倍以上，8並列の場合は1.4倍となりました
TIP4Pポテンシャルの場合相対的にGPUが低速な理由は，PPPM法によるクーロン相互作用の計算がGPUに対応していないことです。GPU化されたPPPM法を用いることができれば改善する可能性があります。ただしマニュアルによるとGPU化がなされているのは一部で，アルゴリズム上必須の高速フーリエ変換はCPUで行われるとのことなので，どこまで効果があるかは試してみないと分かりません。
LAMMPSは，ログファイルの最後の方でどの演算にどのくらいの時間がかかったかを記したレポートを出力してくれます。TIP4Pの場合にその内容を紹介したいと思います。
GPU：
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.83107 | 0.83107 | 0.83107 | 0.0 | 11.90
Bond | 0.12234 | 0.12234 | 0.12234 | 0.0 | 1.75
Kspace | 3.7883 | 3.7883 | 3.7883 | 0.0 | 54.25
Neigh | 0.00084586 | 0.00084586 | 0.00084586 | 0.0 | 0.01
Comm | 0.075925 | 0.075925 | 0.075925 | 0.0 | 1.09
Output | 0.0044897 | 0.0044897 | 0.0044897 | 0.0 | 0.06
Modify | 2.1266 | 2.1266 | 2.1266 | 0.0 | 30.46
Other | | 0.0328 | | | 0.47
CPU 非並列：
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 28.519 | 28.519 | 28.519 | 0.0 | 75.66
Bond | 0.11389 | 0.11389 | 0.11389 | 0.0 | 0.30
Kspace | 3.8728 | 3.8728 | 3.8728 | 0.0 | 10.27
Neigh | 2.891 | 2.891 | 2.891 | 0.0 | 7.67
Comm | 0.079032 | 0.079032 | 0.079032 | 0.0 | 0.21
Output | 0.004175 | 0.004175 | 0.004175 | 0.0 | 0.01
Modify | 2.1671 | 2.1671 | 2.1671 | 0.0 | 5.75
Other | | 0.04474 | | | 0.12
CPU 8並列：
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 4.801 | 5.0749 | 5.4564 | 8.2 | 53.23
Bond | 0.018259 | 0.019488 | 0.02023 | 0.4 | 0.20
Kspace | 1.9346 | 2.3158 | 2.5893 | 12.1 | 24.29
Neigh | 0.56025 | 0.56047 | 0.56073 | 0.0 | 5.88
Comm | 0.15301 | 0.15364 | 0.15407 | 0.1 | 1.61
Output | 0.0017384 | 0.0019381 | 0.0026595 | 0.6 | 0.02
Modify | 1.2191 | 1.2713 | 1.302 | 2.4 | 13.33
Other | | 0.1373 | | | 1.44
この結果について何点かコメントしたいと思います。
- GPUでは
PairおよびNeighの計算時間が際立って短いです。その他の部分はCPU非並列の場合とほぼ同等で，GPU化がなされていないと考えられます - CPU非並列と8並列を比較すると，
Pair,Neighなどの演算は5～6倍程度高速化されています。Kspaceも速くなってはいますが，上であげた二つほどではなく2倍程度です。これは高効率の並列化が比較的難しい高速フーリエ変換を含む処理だからではないかと考えられます
むすび
本稿では，以前構築したGPU環境を用いて古典分子動力学シミュレーターLAMMPSによるテスト計算を行ってみました。コンパイルや実行の仕方を説明し，EAMポテンシャルを用いたbccFeとTIP4Pポテンシャルを用いた水のシミュレーションの計算時間を計測してみました。その結果，グラフィックス処理用でミドルクラスの比較的安価なGPUでも高速な計算を行えることが確認できました。特に，長距離相互作用を伴わない場合に大きな効果を発揮しました。今後機会があればGPGPU専用製品を用いた検証も行っていきたいと思います。