CUDA を用いた Pure SU(3)gauge theory のための Twisted
Gradient Flow 計算の加速石川健一（広島大学理学研究科）

K.-I.I., I.Kanamori, Y.Murakami,A.Nakamura,M.Okawa,R.Ueno,”Non-perturbative determination of the Λ-parameter in the pureSU(3) gauge theory from the twisted gradient flow coupling,” JHEP
12(2017)067 [arXiv:1702.06289]
での計算に用いたプログラムの加速について紹介します。
2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光1

もくじ
1. 格子QCD（SU(3)ゲージ理論）のTwisted Gradient
Flow 計算とは
2. 計算すべき式とどういうふうに加速するか
3. 大体の結果
4. まとめ2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光2１．格子QCD（SU(3)ゲージ理論）のTwisted Gradient Flow 計算とは• 格子QCDでも純SU(3)ゲージ理論という、クォークのいない世界のQCDでの計算です。

• グルーオンだけの世界です。
• 電磁気学で言うところの光子だけの世界です。（電荷を帯びた粒子＝電子や陽子がいない世界）古典場の理論だと真空中のMaxwell 方程式で線形微分方程式なので解析的に解けます。
• QCDの場合、グルーオンだけでもずいぶん違います。
グルーオン自身が色電荷（カラーチャージ）を帯びているからです。QCDのMaxwell 方程式は非線形になります。古典場の理論の範囲でも解析的に解ける訳ではないです。
• 純SU(3)ゲージ理論の量子理論は更に複雑ですが、通常のQCD(クォークのいる世界)の性質の重要なもの（漸近自由性、閉じ込め）を持っているようなので、調べる価値があります。
2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光

水素原子

光子

電子

陽子

3

１．格子QCD（SU(3)ゲージ理論）のTwisted Gradient Flow 計算とは• 格子QCDでも純SU(3)ゲージ理論という、クォークのいない世界のQCDでの計算です。

グルーボール• グルーオンだけの世界です。
• 電磁気学で言うところの光子だけの世界です。（電荷を帯びた粒子＝電子や陽子がいない世界）古典場の理論だと真空中のMaxwell 方程式で線形微分方程式なので解析的に解けます。
• QCDの場合、グルーオンだけでもずいぶん違います。
グルーオン自身が色電荷（カラーチャージ）を帯びているからです。QCDのMaxwell 方程式は非線形になります。古典場の理論の範囲でも解析的に解ける訳では
ないです。

グルーオン

光子

• 純SU(3)ゲージ理論の量子理論は更に複雑ですが、通常のQCD(クォークのいる世界)の性質の重要なもの（漸近自由性、閉じ込め）を持っているようなので、調べる価値があります。
2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光4１．格子QCD（SU(3)ゲージ理論）のTwisted Gradient Flow 計算とは
• 純SU(3)ゲージ理論の量子理論は更に複雑ですが、通常のQCD(クォークのいる世界)の性質の重要なもの（漸近自由性、閉じ込め）を持っているようなので、調べる価値があります。

𝑔• グルーオン同士の相互作用の強さ g (=クォークとグルーオンの相互作用の強さと等しいけど今はクォークがいないので忘れてください。でも基本的に電子と光子の相互作用の強さ e 電荷のQCD版
です)• g は量子力学の効果で、距離に応じてわずか（対数的）に変化します。（距離＝エネルギーの逆数）• QCDでは短距離で弱くなります。長距離で強くなります。（量子電気力学はこの逆の振る舞い）

2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光5１．格子QCD（SU(3)ゲージ理論）のTwisted Gradient Flow 計算とは
• 純SU(3)ゲージ理論の量子理論は更に複雑ですが、通常のQCD(クォークのいる世界)の性質の重要なもの（漸近自由性、閉じ込め）を持っているようなので、調べる価値があります。

𝑔• QCDでは短距離で弱くなります。長距離で強くなります。（量子電気力学はこの逆の振る舞い）• この近いところで力が弱くなる振る舞いは「漸近自由性」と呼ばれます。
• 結合定数 g QCDの相互作用の強さを表す量（無次元量）𝑔2
𝛼𝑠 =
4𝜋
𝛼𝑠 (𝑞)のエネルギースケール依存性の図。 𝑞はエネルギース
ケール。 𝑞の大きさで対数的に変化する。この図は純SU(3)計算で求められた最初の理論計算の例[Luscher,Sommer,Weisz,Wolff NPB413(1994)481.]
2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光6１．格子QCD（SU(3)ゲージ理論）のTwisted Gradient Flow 計算とは
• 純SU(3)ゲージ理論の量子理論は更に複雑ですが、通常のQCD(クォークのいる世界)の性質の重要なもの（漸近自由性、閉じ込め）を持っているようなので、調べる価値があります。

𝑔• QCDでは短距離で弱くなります。長距離で強くなります。（量子電気力学はこの逆の振る舞い）• この近いところで力が弱くなる振る舞いは「漸近自由性」と呼ばれます。

• 結合定数 g QCDの相互作用の強さを表す量
図の横軸の絶対値（GeV単位）を決める基本パラメータがΛ「ラムダ」（無次元量）
𝑔2
パラメータと呼ばれる。QCDの基本定数。 𝛼𝑠 = 4𝜋
すべての次元を持つ物理量はこのラムダパラメータのベキに定数をかけ
𝛼𝑠 (𝑞)のエネルギースケール依存性の図。 𝑞はエネルギース
たものになっている。低エネルギーでの物理観測量をラムダパラメータケール。 𝑞の大きさで対数的に変化する。この図は純SU(3)計を単位として測ったものを正確に決めることは理論的に需要算で求められた最初の理論計算の例[Luscher,Sommer,Weisz,Wolff NPB413(1994)481.]
2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光7

１．格子QCD（SU(3)ゲージ理論）のTwisted Gradient Flow 計算とは
• Twisted Gradient Flow 法（以下TGF法）という方法で𝛼𝑠 (𝑞)の𝑞依存性を精度よく計算できることがわかってきています（2006年以降）• 特に低エネルギー側での𝛼𝑠 (𝑞)の値と、そのエネルギースケール𝑞を単位として測った物理量が精度良く計算できると良いです。
エネルギースケール(質量次元1):𝑞

物理量(質量次元１とする)：𝑀𝑞での結合定数の値：𝛼𝑠 (𝑞)𝑞を単位として測った物理量の値：𝑀/𝑞
結合定数の𝑞依存性：𝛼𝑠 (𝑞)ラムダパラメータを単位として測った物理量の値：𝑀/Λ

格子QCDの中でのいろいろな計算方法のうちTGF法という方法を用いてこれらを精度よく計算した
2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光8２．計算すべき式とどういうふうに加速するか• 格子QCDの説明はざっくりです• 4次元正方格子空間のイジング模型のイジングスピン（𝜎 (𝑛) =+ 1, −1）• 4次元正方格子空間の純SU(3)ゲージ理論のグルーオン変数（𝑈𝜇 (𝑛):3x3 ユニタリ行列の４次元ベクトル）• イジング模型分配関数
𝑍𝑖𝑠𝑖𝑛𝑔 =ෑ෍

exp −𝛽𝑆𝑖𝑠𝑖𝑛𝑔 [𝜎]𝑛∈4𝐷 𝐿𝑎𝑡𝑡𝑖𝑐𝑒 𝜎 𝑛 =±1• 純SU(3)ゲージ理論分配関数

𝑍𝑃𝑢𝑟𝑒𝑆𝑈(3) =ෑෑ න𝑑𝑈𝜇 𝑛 exp[−𝛽𝑆𝑃𝑢𝑟𝑒𝑆𝑈(3) [𝑈]]𝑛∈4𝐷 𝐿𝑎𝑡𝑡𝑖𝑐𝑒 𝜇=1,2,3,4 𝑆𝑈 3 𝑚𝑎𝑡𝑟𝑖𝑥• モンテカルロ法による重み付け積分の評価（統計力学と同じ手法が使える）
2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光9２．計算すべき式とどういうふうに加速するか• TGF法の結合定数の定義• モンテカルロ法でグルーオン変数のサンプルが得られたら、各サンプ
ルごとに次の計算を行う。それらのサンプル平均が量子力学的な期待値。すなわち結合定数（観測値、物理量に対応）となる。
モンテカルロ法で得られたグルーオン変数のあるサンプル：{𝑈𝜇 𝑛 }Gradient Flow （勾配流）計算によるグルーオン変数の平滑化（平滑度 𝑡）：
早い。
計算量は
少ない。

𝑑𝑉𝜇 𝑛, 𝑡
2 × 3 𝜕𝑆𝑝𝑢𝑟𝑒𝑆𝑈 3 𝑉
=−
𝑑𝑡
𝛽
𝜕𝑈𝜇 𝑛

𝑉𝜇 (𝑛, 𝑡)を初期値𝑉𝜇 𝑛, 0 = 𝑈𝜇 (𝑛)から平滑度が𝑡 = 𝑐 2 𝐿2 /8になるまで計算する。
平滑化されたグルーオン変数：{𝑉𝜇 𝑛, 𝑡 }遅い。計算量が多い。

2{𝑉𝜇 𝑛, 𝑡 }を用いて TGF法結合定数𝑔𝑇𝐺𝐹
1/𝐿 を求める。統計平均を取る
2
𝑔𝑇𝐺𝐹
1/𝐿 = 𝐶𝑜𝑛𝑠𝑡 × 𝐸 𝑉 ቚ
2018/12/1

𝑡=𝑐 2 𝐿2 /8

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光

𝑐 = 0.3, 𝐿 : 4次元格子空間の全体の一辺の大きさ10２．計算すべき式とどういうふうに加速するか• TGF法の結合定数の定義• モンテカルロ法でグルーオン変数のサンプルが得られたら、各サンプ
ルごとに次の計算を行う。それらのサンプル平均が量子力学的な期待値。すなわち結合定数（観測値、物理量に対応）となる。
モンテカルロ法で得られたグルーオン変数のあるサンプル：{𝑈𝜇 𝑛 }Gradient Flow （勾配流）計算によるグルーオン変数の平滑化（平滑度 𝑡）：
早い。
計算量は
少ない。

𝑑𝑉𝜇 𝑛, 𝑡
2 × 3 𝜕𝑆𝑝𝑢𝑟𝑒𝑆𝑈 3 𝑉
=−
𝑑𝑡
𝛽
𝜕𝑈𝜇 𝑛

𝑉𝜇 (𝑛, 𝑡)

こいつを速くする

を初期値𝑉𝜇 𝑛, 0 = 𝑈𝜇 (𝑛)から平滑度が𝑡 = 𝑐 2 𝐿2 /8になるまで計算する。
平滑化されたグルーオン変数：{𝑉𝜇 𝑛, 𝑡 }遅い。計算量が多い。

2{𝑉𝜇 𝑛, 𝑡 }を用いて TGF法結合定数𝑔𝑇𝐺𝐹
1/𝐿 を求める。統計平均を取る
2
𝑔𝑇𝐺𝐹
1/𝐿 = 𝐶𝑜𝑛𝑠𝑡 × 𝐸 𝑉 ቚ
2018/12/1

𝑡=𝑐 2 𝐿2 /8

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光

𝑐 = 0.3, 𝐿 : 4次元格子空間の全体の一辺の大きさ11２．計算すべき式とどういうふうに加速するか• TGF法の結合定数の定義• モンテカルロ法でグルーオン変数のサンプルが得られたら、各サンプ
ルごとに次の計算を行う。それらのサンプル平均が量子力学的な期待値。すなわち結合定数（観測値、物理量に対応）となる。
モンテカルロ法で得られたグルーオン変数のあるサンプル：{𝑈𝜇 𝑛 }Gradient Flow （勾配流）計算によるグルーオン変数の平滑化（平滑度 𝑡）：
早い。
計算量は
少ない。

𝑑𝑉𝜇 𝑛, 𝑡
2 × 3 𝜕𝑆𝑝𝑢𝑟𝑒𝑆𝑈 3 𝑉
=−
𝑑𝑡
𝛽
𝜕𝑈𝜇 𝑛

𝑉𝜇 (𝑛, 𝑡)

こいつをCUDA
で速くする

を初期値𝑉𝜇 𝑛, 0 = 𝑈𝜇 (𝑛)から平滑度が𝑡 = 𝑐 2 𝐿2 /8になるまで計算する。
平滑化されたグルーオン変数：{𝑉𝜇 𝑛, 𝑡 }遅い。計算量が多い。

2{𝑉𝜇 𝑛, 𝑡 }を用いて TGF法結合定数𝑔𝑇𝐺𝐹
1/𝐿 を求める。統計平均を取る
2
𝑔𝑇𝐺𝐹
1/𝐿 = 𝐶𝑜𝑛𝑠𝑡 × 𝐸 𝑉 ቚ
2018/12/1

𝑡=𝑐 2 𝐿2 /8

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光

𝑐 = 0.3, 𝐿 : 4次元格子空間の全体の一辺の大きさ12２．計算すべき式とどういうふうに加速するか• TGF法の結合定数の定義
ホストCPUが遊ぶので• モンテカルロ法でグルーオン変数のサンプルが得られたら、各サンプ
モンテカルロと物理量
ルごとに次の計算を行う。それらのサンプル平均が量子力学的な期待計算とGPU（CUDA）
値。すなわち結合定数（観測値、物理量に対応）となる。
計算を同時に行う。
モンテカルロ法で得られたグルーオン変数のあるサンプル：{𝑈𝜇 𝑛 }Gradient Flow （勾配流）計算によるグルーオン変数の平滑化（平滑度 𝑡）：
早い。
計算量は
少ない。

𝑑𝑉𝜇 𝑛, 𝑡
2 × 3 𝜕𝑆𝑝𝑢𝑟𝑒𝑆𝑈 3 𝑉
=−
𝑑𝑡
𝛽
𝜕𝑈𝜇 𝑛

𝑉𝜇 (𝑛, 𝑡)

こいつをCUDA
で速くする

を初期値𝑉𝜇 𝑛, 0 = 𝑈𝜇 (𝑛)から平滑度が𝑡 = 𝑐 2 𝐿2 /8になるまで計算する。
平滑化されたグルーオン変数：{𝑉𝜇 𝑛, 𝑡 }遅い。計算量が多い。

2{𝑉𝜇 𝑛, 𝑡 }を用いて TGF法結合定数𝑔𝑇𝐺𝐹
1/𝐿 を求める。統計平均を取る
2
𝑔𝑇𝐺𝐹
1/𝐿 = 𝐶𝑜𝑛𝑠𝑡 × 𝐸 𝑉 ቚ
2018/12/1

𝑡=𝑐 2 𝐿2 /8

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光

𝑐 = 0.3, 𝐿 : 4次元格子空間の全体の一辺の大きさ13２．計算すべき式とどういうふうに加速するか• TGF法の結合定数計算の流れ（time line）• 初期バージョン：CPU Fortran バージョン• 格子サイズ 36^4 の場合• 計算機：九州大学 tatara （2017年：今もう無いです。）
１６コアOpenMP

１ノード

モンテカルロ法計算：初期配位{𝑈𝜇 𝑛 }から新しい配位{𝑈𝜇 𝑛 }を作る繰りGradient Flow （勾配流）計算により{𝑈𝜇 𝑛 }から 𝑉𝜇 𝑛, 𝑡 を作る

返す

モンテカルロ計算時間 1 : 5 勾配流計算時間2{𝑉𝜇 𝑛, 𝑡 }を用いて TGF法結合定数𝑔𝑇𝐺𝐹
1/𝐿 計算パラメータスキャンと高統計による高精度計算をしたいので 36^4格子で１パラメータで最大４万回これを繰り返す必要が出てきた
2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光14

２．計算すべき式とどういうふうに加速するか• TGF法の結合定数計算の流れ（time line）• バージョン2：CPU Fortran + 1GPU CUDA (v7.5)逐次バージョン• 格子サイズ 36^4 の場合• 計算機：九州大学 tatara （2017年：今もう無いです。） １ノード
１６コアOpenMP Tesla K20m （SM3.5 : Kepler 世代）

モンテカルロ法計算：初期配位{𝑈𝜇 𝑛 }から新しい配位{𝑈𝜇 𝑛 }を作る
繰

Gradient Flow （勾配流）計算により{𝑈𝜇 𝑛 }から 𝑉𝜇 𝑛, 𝑡 を作るり2{𝑉𝜇 𝑛, 𝑡 }を用いて TGF法結合定数𝑔𝑇𝐺𝐹
1/𝐿 計算CUDA加速返
モンテカルロ計算時間 1 : 1.08 勾配流計算時間

すパラメータスキャンと高統計による高精度計算をしたいので 36^4格子で１パラメータで最大４万回これを繰り返す必要が出てきた
2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光15

２．計算すべき式とどういうふうに加速するか• TGF法の結合定数計算の流れ（time line）• バージョン3：CPU Fortran + 1GPU CUDA 並行バージョン• 格子サイズ 36^4 の場合• 計算機：九州大学 tatara （2017年：今もう無いです。） １ノード
１６コアOpenMP Tesla K20m （SM3.5 : Kepler 世代）
モンテカルロ法計算
繰

Gradient Flow （勾配流）計算2
TGF法結合定数𝑔𝑇𝐺𝐹
1/𝐿 計算CUDA側が空いていたら配位{𝑈𝜇 𝑛 }を送りつけて計算させる
CUDA側が空くまで待つり返すGPU並行CUDA加速

勾配流計算時間にモンテカルロ計算が隠蔽されたパラメータスキャンと高統計による高精度計算をしたいので 36^4格子で１パラメータで最大４万回これを繰り返す必要が出てきた
2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光16

２．計算すべき式とどういうふうに加速するか• 並行計算のやり方• モンテカルロ計算はOpenMP スレッド並列しつつ、勾配流計算をGPUに行わせたい。
• さらにGPUは単独では動かなくCPUから制御する必要がある。
• CUDA の stream の概念である並行性の概念があるがstream の管理はCPUが担うのでようわからん。
• OpenMPのスレッド並列をタスク並列と混ぜて使う使い方は良うわか
らん。
• MPI並列でマスター（モンテカルロ計算）・スレーブ（GPUを制御して勾配流と結合定数を計算）形式で実装する。（MPMD (Multiple
Program Multiple Data モード）
モンテカルロ法計算
繰
CUDA側が空いていたら配位{𝑈𝜇 𝑛 }を送りつけて計算させる
CUDA側が空くまで待つり返す2018/12/1

Gradient Flow （勾配流）計算2
TGF法結合定数𝑔𝑇𝐺𝐹
1/𝐿 計算GPU並行CUDA加速

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光
勾配流計算時間にモンテカルロ計算が隠蔽された

17

２．計算すべき式とどういうふうに加速するか• CUDAプログラミングについては省略（前半の講義）• Fortran で正しく書かれている関数を逐次ほぼそのまま CUDA に移植。
OpenMP+ベクトル化の経験があれば、おんなじような感じでCUDAのthread –block 構造へデータレイアウトの変更を行う。
• Fotran 2003 の ISO_C_FOTRAN module により C の関数をFortran側で呼べるので混ぜるのは簡単。CUDA はC++なので公開する関数はすべて extern “C” つけとく。もちろんFortranなので引数は pointer
渡しとしておく（値渡しもできるけど）。

モンテカルロ法計算
繰
CUDA側が空いていたら配位{𝑈𝜇 𝑛 }を送りつけて計算させる
CUDA側が空くまで待つり返す2018/12/1

Gradient Flow （勾配流）計算2
TGF法結合定数𝑔𝑇𝐺𝐹
1/𝐿 計算GPU並行CUDA加速

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光
勾配流計算時間にモンテカルロ計算が隠蔽された

18

２．計算すべき式とどういうふうに加速するか• MPIによる並行計算のやり方• コミュニケータMPI_COMM_WORLD を切り分けてタスクごとに通信コ
ミュニケータを作る。今の場合3つ必要。
1. モンテカルロ法の複数ノード間MPI並列のための通信コミュニケータと勾配流計
算の複数ノード間MPI並列のための通信コミュニケータ
2. モンテカルロ法計算プロセス（マスター）と勾配流計算プロセス（スレーブ）間
でコマンド送受信とグルオン配位転送のための通信コミュニケータ
• 論文のデータを計算するときは１ノード計算をパラメータや乱数を変えて複数走らせたので上記 1. と 2. は結局使いませんでした。プログラムは一応動くようにはしましたが。

• １ノードに２プロセス載せて、一つはマスター（モンテカルロ計算）もう一つはスレーブ（勾配流計算）とする。OpenMP thread の割り当て
はマスターは 16 OMPthread, スレーブは 1 OMP thread 。
• スレーブのCPUはほとんど計算しなくてコマンド受付・配位受付・GPU制御を行う。こちらで結合定数まで計算してファイルに書き込む。
• マスターはthread 並列でモンテカルロ計算を行い、新しいグルーオン
配位ができたら、スレーブに問い合わせて配位を送って勾配流を計算させる。
2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光19

２．計算すべき式とどういうふうに加速するか• MPIによる並行計算のやり方
モンテカルロ法計算
HOST #1
CPU

HOST #2CPU

HOST #3CPU

HOST #3
CPU

2018/12/1

勾配流+結合定数計算

GPU

GPU

GPU

GPU

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光20

２．計算すべき式とどういうふうに加速するか• MPIによる並行計算のやり方
モンテカルロ法計算
HOST #1GPU

CPU

RANK #0

RANK #1

HOST #2

RANK #2

RANK #3

HOST #3

今の場合はランク番号の連続する偶奇が近くなるように割り当てた

GPU

CPU

RANK #4

RANK #5

HOST #3GPU

RANK #6

2018/12/1

最初に各MPIランクをどのようにホストに乗せていくかのランク
マップは hostfile とか使用しているシステムのマニュアルに従う

GPU

CPU

CPU

勾配流+結合定数計算

RANK #7

MPI_COMM_WORLD

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光21

２．計算すべき式とどういうふうに加速するか• MPIによる並行計算のやり方
モンテカルロ法計算
ｍ_host_commｍ_host_comm
HOST #1GPU

CPU

RANK #0

RANK #0

HOST #2

GPU

CPU

Color 0RANK #1

RANK #1

HOST #3

RANK #2

RANK #2

HOST #3GPU

RANK #3

2018/12/1

Color 1

GPU

CPU

CPU

勾配流+結合定数計算

RANK #3

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光

MPI_COM_SPLIT で
MPI_COMM_WORLDを２つに分割
この割り方はホスト間通信。Color はSPLIT
用語
通信はcolor 毎に閉じる。
22

２．計算すべき式とどういうふうに加速するか• MPIによる並行計算のやり方
モンテカルロ法計算
ｍ_h2g_comm

HOST #1

Color 0GPUCPURANK #1

RANK #0

ｍ_h2g_comm

HOST #2

Color 1GPUCPURANK #1

RANK #0

ｍ_h2g_comm

HOST #3

Color 2勾配流+結合定数計算

GPU

CPU

RANK #1

RANK #0

ｍ_h2g_comm
HOST #3

Color 3
2018/12/1

CPU

GPU

RANK #0

RANK #1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光

MPI_COM_SPLIT で
MPI_COMM_WORLDを4つに分割
この割り方はマスター
スレーブ間通信。
Color はSPLIT 用語
通信はcolor 毎に閉じる。
23

２．計算すべき式とどういうふうに加速するか• MPIによる並行計算のやり方
MPI_Send モンテカルロ法計算
で送る場合
が多い
いろいろな関数• スレーブ状態問い合わせ• スレーブ初期化命令転送したあと関連するデータ転送• 勾配流・結合定数計算命令を転送した
あと関連するデータ転送• スレーブ終了化命令転送• 全体終了命令転送

勾配流+結合定数計算 MPI_Recv

で待つ場合
こちらの状態は空が多い[コマンド受付ループ]｛[m_h2g_comm でマスターからのコマンドを受け付けるループ]{コマンドが来るまで待つ}[受け付けたコマンドに応じて仕事を分岐するスイッチ]{
状態問い合わせ{こちらの状態を返す}
初期化コマンド｛こちらを初期化する｝
勾配流・結合定数計算コマンド｛計算する。計
算が終わったらこちらの状態は空とする｝

これらの関数を適宜発行終了化コマンド｛こちらを終了化する｝
全体終了コマンド｛MPI_Finalize してコマンド
受付ループを抜ける｝}｝

m_h2g_comm コミュニケータでの RANK 0 と RANK 1 の間での MPI_Send,
MPI_Recv を用いたコマンドのやり取りやデータ転送を記述。モンテカルロ法
計算の実行ファイルと勾配流+結合定数計算の実行ファイルは異なる MPMD 形式の並列計算。場面毎に何を送って何を受けるかがわかっているので場面毎にそれぞれ Send と Recv を時系列に合わせて記述しておく。
2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光24

３．大体の結果（計算性能の）• 九大の TATARA システムは2017年度で終わっていて、現在はITO システムになっています。

• 現在のITO システムBにもGPUがついていてそれはPascal なのでめっちゃ速いと思うのですが、すでに物理の計算は終わっているので特にITOでの性能測定はやっていません。
• ですので2017年の6月頃に計測したデータをもとにお話します。
• 2017年の6月にはまだ 筑波大CCSのHA-PACS(TCA)が使えたのでそこでの結果と• 広島大学にあるGPUワークステーションの結果も合わせて紹介
します。
システム名

GPUの数

GPUの種類

CUDA バージョン

Tatara 九大

1 GPU/NODE

Tesla K20m (SM3.5:Kepler)

CUDA 7.5

HA-PACS(TC) 筑波大

4 GPU/NODE

Tesla K20X (SM5.3:Kepler)

CUDA 8.0

広大ワークステーション

3 GPU/NODE

Tesla C2075 (SM2:Fermi)

CUDA 8.0

2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光25

３．大体の結果（計算性能の）• CUDA計算のカーネル
カーネル名

役割

cu_force_u

勾配を計算する

cu_update_u

グルーオン場を更新する

cu_update_force

勾配を更新する

cu_energy

結合定数を計算する

reduce_double_kernel

結合定数を計算するときに総和計算で使う• 勾配流の微分方程式は拡散方程式の形。SU(3)群上のルンゲ
クッタの4次で表の上位3つの関数を組み合わせたループで平滑化していく。

2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光26

３．大体の結果（計算性能の）• CUDA計算のカーネルの演算量 [ 単位： GFLOP ]カーネル名

nv_prof による計測

理論見積り

cu_force_u

18.2

19.376

cu_update_u

4.925

4.683

cu_update_force

0.16120.1612

cu_energy

20.89921.415

reduce_double_kernel

0.000065408---• 実行しての測り方
nvprof --kernels cu_force –m ¥flop_count_dp,flop_count_dp_add,flop_count_dp_mul,flop_count_dp_fma ¥–o prof.out ./main_tgf_coupling
flop_count_dp で倍精度計算数を数える。この例はさらに、足し算と掛け算と積和の内訳も数えている。
計測したいカーネル関数名を指定して実行すれば良い。ここではそれぞれに計測したので5回行っている。
2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光27

３．大体の結果（計算性能の）• CUDA計算のコンパイル時のレジスタ割り当て状況
カーネル名

SM2

SM3.5

cu_force_u

264 bytes stack frame, 624 bytes spill stores,
612 bytes
spill loads Used 63 registers, 56 bytes cmem[0],28 bytes
cmem[16]

0 bytes stack frame, 0 bytes spill stores, 0
bytes spillloads Used 128 registers, 344 bytes cmem[0],40 bytes
cmem[2]

cu_update_u

80 bytes stack frame, 96 bytes spill stores, 96
bytes
spill loads Used 63 registers, 56 bytes cmem[0],68 bytes
cmem[16]

0 bytes stack frame, 0 bytes spill stores, 0
bytes spillloads Used 80 registers, 344 bytes cmem[0],44 bytes
cmem[2]

cu_update_force

0 bytes stack frame, 0 bytes spill stores, 0 bytes
spill
loads Used 37 registers, 72 bytes cmem[0]0 bytes stack frame, 0 bytes spill stores, 0
bytes spillloads Used 44 registers, 360 bytes cmem[0]cu_energy

768 bytes stack frame, 1592 bytes spill stores,
1748 bytes
spill loads Used 63 registers, 144 bytes smem,
64 bytes
cmem[0], 28 bytes cmem[16]48 bytes stack frame, 0 bytes spill stores, 0
bytes
spill loads Used 170 registers, 144 bytes smem,
352 bytes
cmem[0], 28 bytes cmem[2]• レジスタスピルすると性能劣化する。GPU世代ごとにレジスタ数が異なるのでレジ
スタ割り当てとかが異なる。コンパイル時に確認することができる。

• SM2 は1スレッド 63 レジスタまでしか使えないため、スピルが多く発生。
2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光28

３．大体の結果（計算性能の）• 計算性能（広島大学GPUワークステーション、SM2 Fermi世代）• ノード内GPU並列も実装している。スレーブプロセス内で GPU 並列。

• ３GPUまでよくスケールしている。単位は秒とGLOPS
2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光

格子サイズを 36^4 に固定した強スケーリング
29

３．大体の結果（計算性能の）• 計算性能（TATARA、SM3.5 Kepler世代）

格子サイズは 36^4

• 単位は秒とGLOPS• Fermi 3 GPU に比べてすでに、Kepler 1 GPU が1.6倍速い[1824 sec VS 823 sec]
• Fermi 1GPU と Kepler 1GPU では

3.8 倍Kepler が速い [ 3144 sec VS 823 sec]

• force と energy 計算が特にFermiは遅い（レジスタスピル多い）

2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光30

３．大体の結果（計算性能の）• 計算性能（HA-PACS(TCA)、SM3.5 Kepler世代）• スレーブプロセス内で GPU 並列。
• 計算カーネル性能はよくスケールしているが、袖通信(copy_boundary)のスケー
リングが悪い。

• ２GPUの計算時間はすでにモンテカルロ計算より速い
2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光

格子サイズを 36^4 に固定した強スケーリング 31

３．大体の結果（計算性能の）• TATARAでの並行計算• TATARA は 1GPU/NODE であるが、1GPU の勾配流結合定数
計算の時間とモンテカルロ計算の時間はだいたい同じぐらいとなった。
• 逐次実行の場合• 一つのグルオン配位を作成するのにかかるモンテカルロ計算の時間：
175 秒• 勾配流計算して結合定数を計測する時間：198秒• 合計：373秒• 並行実行の場合• マスター：一つのグルオン配位を作成するのにかかるモンテカルロ計
算の時間： 183 秒• マスター：からみた勾配流計算の時間: 16秒• マスター合計：199秒• スレーブ：勾配流計算して結合定数を計測する時間：199秒• 並行実行で 373秒=> 199 秒へ高速化できた（1.9倍）
2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光32

４．まとめ• TGF法の結合定数計算をCUDA で加速した。
• クォークのいないQCD（純SU(3)ゲージ理論）なのでグルーオ
ン配位のモンテカルロ計算は比較的軽量。
• 一方TGF結合定数計算のための勾配流計算は計算量が多い。
• GPUを用いて加速したが、まだモンテカルロ計算と同じ程度の経過時間であった。
• GPUとCPUを並行に使うことで計算を加速した。
• タスク並列（マスター・スレーブ）でモンテカルロ計算と勾配流計算を切り分けた• MPIを用いてMPMD 形式でマスタースレーブ関係を実装した。
• 並行実行によりCUDA逐次よりも 1.9 倍速くなった。

• 以上です。ありがとうございました。
2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光33

コード（CUDAFortran混成）• データ構造• ゲージリンク変数（グルーオン変数）を成分とする4次元ベクトル

3x3 行列（特殊ユニタリー行列）

𝑈𝜇𝑎𝑏 𝑛 , 𝑎, 𝑏 = 1,2,3, 𝜇 = 1,2,3,4, 𝑛:4次元格子座標 n=(ix,iy,iz,it)• 格子サイズ：NX,NY,NZ,NT （4次元の各一辺の長さ）• 次元: NDIM=4• Fortran:

complex(8) :: ue(0:NX/2,3,3,0:NY+1,0:NZ+1,0:NT+1,NDIM)
complex(8) :: uo(0:NX/2,3,3,0:NY+1,0:NZ+1,0:NT+1,NDIM)
X方向を2で割って4次元格子点の偶奇((nx+ny+nz+nt % 2))に応じて配列ue,uo に分けて保存している。

• 偶奇（Red-Black）色分けしている。CPUでの計算に適した配列構造。
• 格子点番号は１が原点。NX が端っこ。0と NX+1 は境界条件用の袖領域。
今は周期境界条件。
2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光34

コード（CUDAFortran混成）

ISO_C_BINDING モジュールを使ってFortran構造体とC言語構
造体でどちらでも使えるデータ構造を定義した （BIND(C)キーワード）• データ構造• NVIDIA CUDA の計算に適したデータ構造にする必要がある。
• Thread 長は長い方が大体よい。（32の倍数だとなおよい）• Fortran:• 配列の連続する部分はこのままだと NX/2 なので短い。したがって、ここでは X 方向と Y 方向を1次元化するような感じで長い連続する添え字を一番内側に持ってくるようにデータのレ
イアウトを変更する。

!
! CUDA gauge field!
type, BIND(C) :: cuDComplex_threadcomplex(DP) :: d(0:NTHREADS-1)
end type• GPU並列も考慮してT 方向をGPU分割する。

type, BIND(C) :: cu_thread_gauge
type(cuDComplex_thread) :: u(COL,COL)
end type
type, BIND(C) :: cu_sfield_gauge_gtype(cu_thread_gauge) :: v(NDIM)
end type

type, BIND(C) :: cu_vfield_gauge_gtype(cu_sfield_gauge_g) :: s(0:NBLOCKS_ARRAY-1)
integer(C_INT) :: gid
integer(C_INT) :: idummy(3)end type!
! fortran gauge field!
type, BIND(C) :: f_vfield_gauge
type(cu_vfield_gauge_g) :: g(0:NGPUS-1)
end type
2018/12/1

複素数が X方向連続＊Y方向の1部部分の連続する格子点分集まった集合
が 3x3 こ（3x3特殊ユニタリ行列）連続して集まった集合が ４こ(4次元ベクトル）連続して集まった集合
が Y方向ブロック番号＊Z方向＊（T方向/GPUあたり）格子点分連続して集まった集合

が GPU個、集まった集合＝グルーオン場全体データ

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光35

コード（CUDAFortran混成）

ISO_C_BINDING モジュールを使ってFortran構造体とC言語構
造体でどちらでも使えるデータ構造を定義した （BIND(C)キーワード）• データ構造• NVIDIA CUDA の計算に適したデータ構造にする必要がある。

• Fortran:• 古いデータ構造と CUDA用のデータ構造変更プログラムを書いてよくチェックする。
! subroutine convert_field_gauge(ue,uo,cu)!
! convert even/odd orderd gauge field to CUDA array format.
!
! ue : even site gauge field! uo : odd site gauge field! cu : gauge field in CUDA array format.
!
use lattice_mod
implicit none
complex(DP), intent(in) :: ue(0:NXH,NG,NG,0:NY1,0:NZ1,0:NT1,NDIM)
complex(DP), intent(in) :: uo(0:NXH,NG,NG,0:NY1,0:NZ1,0:NT1,NDIM)
type(f_vfield_gauge), intent(out) :: cu
integer :: mu,nu,it,iz,iy,ix,ixb,ieo,ic,jc
integer :: iblk,igpu,igt,iby,ithy,ithd!$OMP PARALLEL DO PRIVATE(mu,it,iz,iy,ieo,ixb,ix,igpu,igt,iby,ithy,ithd,iblk,ic,jc) COLLAPSE(4)do mu = 1,NDIM
do it = 1,NTdo iz = 1,NZ
do iy = 1,NY
ieo = mod(iy+iz+it+1+ieopu,2)
do ixb=ieo,NXH+ieo-1
ix = 2*ixb+1-ieo
igpu =((it-1)/NGT )igt = mod((it-1),NGT )iby =((iy-1)/NTHY)
ithy = mod((iy-1),NTHY)
ithd = ix-1 + (NTHX)*(ithy)
iblk = iby + (NBY)*( iz-1+1 + (NZ+2)*(igt+1) )2018/12/1• 配列添え字の変換

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光36

コード（CUDAFortran混成）

ISO_C_BINDING モジュールを使ってFortran構造体とC言語構
造体でどちらでも使えるデータ構造を定義した （BIND(C)キーワード）• データ構造• NVIDIA CUDA の計算に適したデータ構造にする必要がある。

• Fortran:• 古いデータ構造と CUDA用のデータ構造変更プログラムを書いてよくチェックする。
igpu =((it-1)/NGT )igt = mod((it-1),NGT )iby =((iy-1)/NTHY)
ithy = mod((iy-1),NTHY)
ithd = ix-1 + (NTHX)*(ithy)
iblk = iby + (NBY)*( iz-1+1 + (NZ+2)*(igt+1) )!DIR$ UNROLL(3)do jc = 1,NG
!DIR$ UNROLL(3)do ic = 1,NG
cu%g(igpu)%s(iblk)%v(mu)%u(ic,jc)%d(ithd) = ue(ixb,ic,jc,iy,iz,it,mu)
enddo
enddo

• データの詰め替え

enddo
ieo = mod(iy+iz+it+0+ieopu,2)
do ixb=ieo,NXH+ieo-1
ix = 2*ixb+1-ieo
igpu =((it-1)/NGT )igt = mod((it-1),NGT )iby =((iy-1)/NTHY)
ithy = mod((iy-1),NTHY)
ithd = ix-1 + (NTHX)*(ithy)
iblk = iby + (NBY)*( iz-1+1 + (NZ+2)*(igt+1) )!DIR$ UNROLL(3)do jc = 1,NG
!DIR$ UNROLL(3)do ic = 1,NG
cu%g(igpu)%s(iblk)%v(mu)%u(ic,jc)%d(ithd) = uo(ixb,ic,jc,iy,iz,it,mu)
enddo
enddo
第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光2018/12/137

コード（CUDAFortran混成）

ISO_C_BINDING モジュールを使ってFortran構造体とC言語構
造体でどちらでも使えるデータ構造を定義した （BIND(C)キーワード）• データ構造• NVIDIA CUDA の計算に適したデータ構造にする必要がある。
• C（CUDA）側のデータ構造体://!
! CUDA gauge field!
type, BIND(C) :: cuDComplex_threadcomplex(DP) :: d(0:NTHREADS-1)
end type
type, BIND(C) :: cu_thread_gauge
type(cuDComplex_thread) :: u(COL,COL)
end type
type, BIND(C) :: cu_sfield_gauge_gtype(cu_thread_gauge) :: v(NDIM)
end type// CUDA gauge field//
typedef struct _ALIG_ {cuDoubleComplex d[_NTHREADS];} cuDComplex_thread;

CUDA

thread で並列化

typedef struct _ALIG_ {cuDComplex_thread u[COL*COL];} cu_thread_gauge;
typedef struct _ALIG_ {cu_thread_gauge v[NDIM];} cu_sfield_gauge_g;typedef struct _ALIG_ {cu_sfield_gauge_g s[_NBLOCKS_ARRAY];
type, BIND(C) :: cu_vfield_gauge_g
int gid;type(cu_sfield_gauge_g) :: s(0:NBLOCKS_ARRAY-1)
CUDA block で並列化
int idummy[3];integer(C_INT) :: gid} cu_vfield_gauge_g;
integer(C_INT) :: idummy(3)CUDAでは各GPU毎
end type
CとFortranで同じ構造体にGPU上にメモリをtypedef struct _ALIG_ {メモリ上のデータが両言確保しないといけな
cu_vfield_gauge_g *g[_NGPUS];
語で同様にアクセスして!
いのでポインタ型も}cu_vfield_gauge;
同じ値が保障される! fortran gauge field
用意しておく!
//type, BIND(C) :: f_vfield_gauge// fortran gauge field
type(cu_vfield_gauge_g) :: g(0:NGPUS-1)//
end type
typedef struct _ALIG_ {cu_vfield_gauge_g g[_NGPUS];}
f_vfield_gauge;
第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光2018/12/138

コード（CUDAFortran混成）• 呼び出しFortran
うこと

ISO_C_BINDING モジュールを使• CUDA言語（またはC言語）のルーチンを Fortran から呼び出す。 Fortranmodule で interface 宣言で CUDA ルーチンの引数並びを宣言しておけばよい。
CUDA側（C++言語ベース）
extern “C” で公開

Fortran module 内で引用仕様(interface)宣言
interface!
! interfaces to CUDA routines!
function get_num_gpu_devices() RESULT(numgpus) BIND(C)importinteger(C_INT) :: numgpus
end function
subroutine new_cuda_tgf(this,z) BIND(C)importtype(cuda_tgf),
intent(inout) :: thistype(f_tfield_twist), intent(in):: z
end subroutine
subroutine delete_cuda_tgf(this) BIND(C)importtype(cuda_tgf), intent(inout) :: thisend subroutine

extern "C" void
cuda_evolve_gdf(cuda_tgf *ctgf, const int *const itraj, const int*const ntau, const double *const dtau){
startTimer(m_etime_evolve_gdf);
double dtau4 = 0.25*(*dtau);#pragma omp parallel forfor ( int igpu = 0; igpu < _NGPUS; ++igpu) {cuda_bind_gpu(igpu);
CUDA_SAFE_CALL(cudaMemcpy(ctgf->d_dt[ igpu],
dtau ,sizeof(double),cudaMemcpyHostToDevice));
CUDA_SAFE_CALL(cudaMemcpy(ctgf>d_dt4[igpu],&dtau4,sizeof(double),cudaMemcpyHostToDevice));}double esym;
double tau;
int imd;

subroutine cuda_evolve_gdf(this,itraj,ntau,dtau) BIND(C)imd = 0;importcuda_copy_boundary_gauge(ctgf->d_u);
CUDAでガンガンtype(cuda_tgf), intent(inout) :: this
計算していくcuda_energy(&esym,ctgf->d_u);
integer(C_INT), intent(in):: itraj… skip …
integer(C_INT), intent(in):: ntau
Fortran側:
real(DP),intent(in):: dtau
Bind(C) をつけてC言語の関数であることを指定。Fortran の引数はC側ではポインタ型であること。
end subroutine
Fortran は名前(Entity)に大文字小文字の区別がないので C 側で大文字と小文字を混ぜた関数名を使っ…skip…
ているときには2018/12/1

Bind(C,name=“CudaEvolveTGF_HOGE_v2”) とか直接 C側の名前を書いておくこと（エイリアス)
第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光39

コード（CUDAFortran混成）

Fortran
うこと

ISO_C_BINDING モジュールを使• 呼び出し• Fortran側で interface 宣言されているCUDAで実装されたサブルーチンを呼べばよい。

subroutine evolve_gdf_cuda_tgf_f(this)!
! cuda side code
! evolve tgf!
type(cuda_tgf), intent(inout) :: thisinteger(C_INT) :: itraj
integer(C_INT) :: ntau
real(DP):: dtautype(f_vfield_gauge), allocatable :: cutype(C_PTR) :: cptr!----------------------------------! set gauge and twist field on GPU!----------------------------------allocate(cu)cptr = C_LOC(cu)call comlib_recv_host(cptr,C_SIZEOF(cu))
call set_cuda_tgf(this,cu)
deallocate(cu)if ( 1 == m_is_initialized ) then!==============================! receive data from MC code!==============================
cptr = C_LOC(itraj)
call comlib_recv_host(cptr,C_SIZEOF(itraj))!==============================! start CUDA code!==============================
call cuda_evolve_gdf(this,itraj,ntau,dtau)

endif
returnend subroutine

GPUでの計算のためGPU側でメモリを確保してグルーオン
場を転送、GPUj上での計算のための初期化をしておく

CUDAでガンガン
計算していくcptr = C_LOC(ntau)
call comlib_recv_host(cptr,C_SIZEOF(ntau))
cptr = C_LOC(dtau)
call comlib_recv_host(cptr,C_SIZEOF(dtau))!----------------------------------! set gauge and twist field on GPU!-----------------------------------

大体こんな感じです。

CPU/GPU 並列のマスター・スレーブモデルのスレーブ側のデータ受信部
2018/12/1

第2回HPC-Phys 勉強会(2018/12/1)@理化学研究所和光40