「不⽼」における分⼦軌道計算プログ
ラムABINIT-MPの⾼速化⼿法の検討
満⽥ 晴紀, 星野 哲也（名⼤）, 望⽉ 祐志（⽴教⼤）,板倉 耕太（FOCUS）, ⽚桐 孝洋（名⼤）, ⼤島 聡史（九⼤）永井 亨, 河合 直聡（名⼤）2023/9/6

分⼦軌道計算プログラム ABINIT-MP• フラグメント分⼦軌道法(FMO)を⽤いた計算シミュレーションソフトウェア• 専⽤GUIによる可視化が容易• Fortranベース、MPIまたはOpenMPによる並列化• HPCI拠点で利⽤可能
ü A64FX：富岳(理研), 不⽼(名⼤), Wisteria(東⼤)ü TSUBASA：SQUID(阪⼤), AOBA-A(東北⼤)ü Intel Xeon系：ITO(九⼤), TSUBAME(東⼯⼤), Grand Chariot(北⼤)など• GPUへの移植も進んでいる• 性能解析のためのミニアプリが提供されている
ü 今回はこのミニアプリを⽤いる• 性能分析や最適化はまだ不⼗分

出典:計算⼯学ナビ ABINIT-MPOpen シリーズ (Ver.2 Rev.4)
http://www.cenav.org/?s=ABINITMP

2

研究⽬的
1. A64FXにおけるABINIT-MPミニアプリの性能評価
ü A64FXは省電⼒に重きを置いたプロセッサであり、省電⼒が重要な今後のプロセッサもA64FXに似た性能特性を持つと考えられるため、A64FXにおける性能評価は重要である
ü ABINIT-MPミニアプリではFMO計算に4つの⼿法を提供しているが、特にA64FXにおいてそれぞれの⽐較評価は⼗分ではない→ 本発表ではA64FXにおいてプロファイラを⽤いて4つの計算⼿法の性能分析を⾏う

2. ベクトル化を念頭においたコード最適化ü 分⼦軌道計算における軌道の組み合わせで81個の個別ルーチンが⽤意されており、ループ中でselect⽂により呼び分けられるため分岐コストが⼤きいü 個別ルーチン内のベクトル化対象ループ⻑が極端に⼩さいü ベクトル化対象ループ内のカットオフのための分岐がベクトル化を妨げている→ 本発表では分岐コストを減らすための実⾏順序ソート⼿法、ループ⻑の確保⼿法、最内ループ内分岐の削減⼿法について提案・評価

3. ⾃動チューニング(AT)適⽤に関する検討
ü 例えばベクトル化を促進するためのループブロッキング幅はプロセッサにより最適値が異なりうる→ 本発表ではATの適⽤可能性を検討。実際のAT適⽤は今後の課題
3

ABINIT-MPのFMO計算⼿法• フラグメント分⼦軌道(FMO)法ü ⽣体⾼分⼦のような巨⼤分⼦について⾮経験的分⼦軌道計算を⾏うことができる
ü 以下の⼿順をとる(右図参照)
1.
2.
3.

分⼦間のハミルトニアンから電⼦分を取り出す
電⼦に関するシュレディンガー⽅程式を⽴てる
シュレディンガー⽅程式をHatree-Fock(HF)近似で解いて系の電⼦エネルギーを求める• ABINIT-MPにおけるFMO計算
ü HF計算において最も重い⼆電⼦積分⽣成処理をObaraらの垂直漸化式
関係(VRR)法により⾏うØ VRR法は軌道⾓運動量のシフトアップによってs型(⾓運動量0), p型(同1),d型(同2)の軌道に関する積分を⽣成するアルゴリズム

ü s/p/d型関数の組み合わせによって34=81種類のルーチンを⽣成、select
⽂によって呼び分ける
Ø s型が4つの[ss-ss]の計算が基準
Ø ベクトル化等の最適化は81種類のルーチン全て書き換える必要がある
4

ABINIT-MPミニアプリの計算⼿法• ABINIT-MPの評価のため2電⼦積分⽣成とその積分を使ったHF計算のみ抜き出したもの
ü ABINIT-MPでは2電⼦積分の⽣成のみで実⾏時間の60%程度を占める• 以下の4⼿法を評価できる
ü ObaraØ ABINIT-MP本体同様、ObaraらのVRR法を⽤いている
Ø 81種類のサブルーチンを呼び分けて2電⼦積分⽣成及びHF計算を⾏うü Obara-vector
Ø NEC TSBASAを念頭におき、Obaraをベースに変更を加えたもの
Ø 多重ループを⼀重化して最内ループを⻑くしているØ 81種類のサブルーチンを呼び分けるのは同様

ü Obara-generalØ 81種類のサブルーチンを汎⽤的な⼀つのサブルーチンとしてまとめたもの

ü Hgp-generalØ Head-Gordonらの⽔平漸化式関係(HRR)法を⽤いたもの
Ø ⾓運動量和が同じ組み合わせに対して⼀括⽣成を⾏うため、基底関数の短縮⻑が⻑い場合にはVRR法に⽐べて演算数の低減が期待できる5

Obaraのプログラム構造
do NF = 1,nf !フラグメントループ
MPI
do I = 1,nshell
MPI/OMPdo J = 1,nshell

ミニアプリ範囲
nshell : 30~100、FMOの場合はその倍

do K =sub.get_tei_pq_fix
1,nshelldo L = 1,nshellcall get_tei_pq

call suberi

index = calcId(I,J,K,L)
select (index)
case (000):
call sub_ssss
case (001):
call sub_sssp…
case (624)call sub_xxxx
end select

enddo
enddo
enddo
enddo
enddo

sub.suberi(ssss~dddd)I/J/K/L_ex: 1 or 3 or 6 Iに依存
do p = 1,I_ex
do q =*sub_xxxx
1,J_ex 一般式
do r = 1,K_ex
do s = 1,L_ex
if (カットオフ)
積分計算小原漸化式(VRR)
enddo
enddo
enddo
enddo

s,p,d,f,g軌道の組み合わせで54=625通りを呼び分け。ただし特別に関数として用意されているのはs,p,d軌道の組
み合わせの34=81通り。今回の入力で利用されるのもs,p,d軌道のみ。
6

Obaraのプログラム構造
do NF = 1,nf !フラグメントループ
MPI
do I = 1,nshell
MPI/OMPdo J = 1,nshell

ミニアプリ範囲
nshell : 30~100、FMOの場合はその倍

do K =sub.get_tei_pq_fix
1,nshelldo L = 1,nshellcall get_tei_pq

enddo
enddo
enddo

sub.suberi
index = calcId(I,J,K,L)(ssss~dddd)I/J/K/L_ex: 1 or 3 or 6 Iに依存
select (index)
do p = 1,I_ex
実質的なループ⻑はほとんど⼀桁（中央値8）case (000):
do q =*sub_xxxx
1,J_ex 一般式
でベクトル化は厳しい。最⼤は1296。
call sub_ssss
do r = 1,K_ex
case (001):
do s = 1,L_ex
call sub_sssp
if (カットオフ)
call suberi
積分計算小原漸化式(VRR)
…
enddo
enddo
case (624)enddocall sub_xxxx
enddo
end selectenddo
s,p,d,f,g軌道の組み合わせで54=625
enddo
通りを呼び分け。ただし特別に関数として用意されているのはs,p,d軌道の組
み合わせの34=81通り。今回の入力で利用されるのもs,p,d軌道のみ。
7

Obara-vectorのプログラム構造
do NF = 1,nf !フラグメントループ
ミニアプリ範囲
MPI
do I = 1,nshell
MPI/OMPnshell : 30~100、FMOの場合はその倍
do J = 1,nshellsub.get_tei_pq_fix
do KLtype= 0,8call get_tei_pq

call suberi

index = calcId(I,J,K,L)
select (index)
case (000):
call sub_ssss_vec
case (001):
call sub_sssp_vec…
case (624)call sub_xxxx_vec
end selectenddo
enddo
enddo
enddo

sub.suberi(ssss~dddd)

do p = 1,I_ex*sub_xxxx 一般式 I/J_ex: 1 or 3 or 6 Iに依存
KL*rs: 150程度
do q = 1,J_ex
do klrs = KL*rs
if (カットオフ)
積分計算小原漸化式(VRR)
enddo
enddo
enddo

s,p,d,f,g軌道の組み合わせで54=625通りを呼び分け。ただし特別に関数として用意されているのはs,p,d軌道の組
み合わせの34=81通り。今回の入力で利用されるのもs,p,d軌道のみ。
8

Obara-vectorのプログラム構造
do NF = 1,nf !フラグメントループ
ミニアプリ範囲
MPI
doルーチン
I = 1,nshellsub_ddddに対する富⼠通コンパイラ最適化出⼒
MPI/OMPnshell : 30~100、FMOの場合はその倍
do J = 1,nshellsub.get_tei_pq_fix
do KLtype= 0,8call get_tei_pq

call suberi

index = calcId(I,J,K,L)
select (index)
case (000):
call sub_ssss_vec
case (001):
call sub_sssp_vec…
case (624)call sub_xxxx_vec
end selectsub.suberi(ssss~dddd)

do p = 1,I_ex*sub_xxxx 一般式 I/J_ex: 1 or 3 or 6 Iに依存
KL*rs: 150程度
do q = 1,J_ex
do klrs = KL*rs
if (カットオフ)
積分計算小原漸化式(VRR)
enddo
enddo
enddo

最内ループ長を確保するためにさらに内側の細かいループを全てunrollingしており、
s,p,d,f,g軌道の組み合わせで54=625
enddo
その結果として大量の変数、ループボディ1万行超。富士通コンパイラの出力では、レ通りを呼び分け。ただし特別に関数としジスタスピル14,392となってベクトル化されない。
enddoて用意されているのはs,p,d軌道の組

enddo
enddo

み合わせの34=81通り。今回の入力で利用されるのもs,p,d軌道のみ。
9ミニアプリの性能評価
A64FX及びIntel Xeon (CascadeLake)における実⾏時間(秒)5回実⾏した際の中央値

データセットObara
A64FX (2.2GHz)

OpenMP 48コア使用
3,379.2 GFLOPS
富士通コンパイラ4.9.0–Kopenmp,fast,nosimd,ocl –Nlst=tIntel Xeon Platinum
8280 (CascadeLake)

OpenMP 1ソケット28コア使用
2,419.2 GFLOPSIntel compiler 2021.7.1 –O3 –fopenmp

6-31G

6-31G*

cc-pVDZ

2.1

5.4

12.6

Obara-vector

5.9

15.8

45.1

Obara-general

185.8

346.4

638.1

Hgp-general

361.8

656.9

1208.0Obara1.2

3.1

7.6

Obara-vector

3.2

11.6

37.0

Obara-general

5.1

13.9

32.9

Hgp-general

3.9

11.7

27.2

• いずれもアダマンタン（C10H16）分⼦に対する計算• 基底関数が6-31G, 6-31G*, cc-pVDZで異なるü 6-31Gではd関数がないので、呼び出されるサブルーチンは24=16パターンのみ
ü 6-31G*, cc-pVDZでは34=81通りのルーチンが呼び出されるü cc-pVDZは⽔分⼦にもp型の分極関数を加えたもので、6-31G系より短縮⻑が⻑く、カットオフが多い• A64FXが全般的に遅いü Obara, Obara-vectorを⽐較すると、1.22-1.85倍ほど遅い
ü いずれもベクトル化率が低い
ü *-generalが特に遅いのは、critical指⽰⽂を使っているため
Ø 最適化の進度がObara, Obara-vectorより遅れていることによるものと思われる
10

提案①：ベクトル化を念頭においた実⾏順ソート
do NF = 1,nf !フラグメントループ
ミニアプリ範囲
MPI!-- ループ前にあらかじめ実行順をソート
Select文のキーとなるindex (s, p, dの組み合わせによって決まる数値）をあらかじめ参照し、その際のループインデックスを覚えておく。
!$omp parallel
do index = 0, 624!$omp doselect (index)
case (000):
call sub_ssss
case (001):
call sub_sssp…
case (624)call sub_xxxx
end selectenddo!$omp end parallel
enddo

外側の4重ループを一重化。(30-100)4 なのでそれなりのサイズのループを
do IJKL = 1, num_xxxx
I = index_to_i(IJKL)
selectによる分岐後に実行でき、ベクトル化を行いやすくなった。
J = index_to_j(IJKL)
発想としてはObara-vectorに似ているが、このループを最内ループ化はしない点で異なる。
K = index_to_k(IJKL)
L = index_to_l(IJKL)
do p = 1,I_ex
do q = 1,J_ex
Sedコマンドによって元ソースから切り貼りを行い、81個の関数を生
do r = 1,K_ex成している
do s = 1,L_ex
if (カットオフ)
積分計算小原漸化式(VRR)
enddo
enddo
enddo
ループの計算コストは一定ではなく負荷不均衡が生じる
enddo
ため、nowaitにより同期コストを減らす。それぞれのサブ
enddo
ルーチンは独立実行可。
!$omp end do nowait
11

提案①のメリット・デメリット• メリットü分岐・関数呼び出しが減るü分岐後のループ⻑が⻑く、ベクトル化などの最適化を適⽤しやすいØ これが今回の実装を⾏なった主⽬的
Ø GPUなどにおいても有効と考えられる• デメリットüソートのオーバーヘッド
Ø ただし、電⼦状態が安定するまでの繰り返し計算の間、ソートは⼀回のみで良いü冗⻑計算の増加
Ø 外側の4重ループを⼀重化した影響

üキャッシュ効率の低下
üメモリアクセス効率の低下
Ø データの格納順を考えることで影響を減らせるが、今後の課題
12提案①実装の性能評価
A64FX及びにおける実⾏時間(秒)6-31G

6-31G*

cc-pVDZ

Obara2.1

5.4

12.6

Obara-vector

5.9

15.8

45.1

提案手法(括弧内はソートの内時間）

2.5 (0.2)

5.8 (0.3)

17.1 (0.6)

データセットA64FX (2.2GHz)

OpenMP 48コア使用
3,379.2 GFLOPS
富士通コンパイラ4.9.0–Kopenmp,fast,nosimd,ocl –Nlst=t• いずれもアダマンタン（C10H16）分⼦に対する計算• 基底関数が6-31G, 6-31G*, cc-pVDZで異なる
ü 6-31Gではd関数がないので、呼び出されるサブルーチンは24=16パターンのみ
ü 6-31G*, cc-pVDZでは34=81通りのルーチンが呼び出されるü cc-pVDZは⽔分⼦にもp型の分極関数を加えたもので、6-31G系より短縮⻑が⻑く、カットオフが多い• 提案⼿法はObaraと⽐較すると少し遅いが、Obara-vectorより速いü ソートは1コアのみで実⾏している。並列化により⾼速化の⾒込みはあるü 少しのオーバーヘッドでベクトル⻑を稼げたことが成果
Ø ただしこの段階では全くベクトル化はされていない

13

提案②：カットオフ分岐の排除及びブロッキング• 分岐排除基本⽅針
ü 条件判定部と計算部にループを分割する
ü 実⾏されるループボディの回数を数え、計算部のループ⻑とする
ü 実⾏される時のループインデックスを記憶しておき、計算ループで参照する
do p = 1,I_ex
do q = 1,J_ex
do r = 1,K_ex
do s = 1,L_ex
if (カットオフ)
積分計算小原漸化式(VRR)
enddo
enddo
enddo
enddo

これで計算ループからは分岐がなくなり、SIMD化が促進されるはずだが、そのままではループ⻑が短すぎる（中央値8）ため、効果が薄い

ix = 0do p = 1,I_ex
do q = 1,J_ex
do r = 1,K_ex
do s = 1,L_ex
if (カットオフ) thenix + 1
p_array(ix) = p
q_array(ix) = q
r_array(ix) = r
s_array(ix) = s
end if
enddo
enddo
enddo
enddo
do i = 1, ix
p = p_array(ix)q = q_array(ix)r = r_array(ix)s = s_array(ix)
積分計算小原漸化式(VRR)
enddo

条件判定ループ

計算ループ

14

提案②：カットオフ分岐の排除及びブロッキング• ブロッキングの⽅針
ü提案①で⼀重化したループに対してブロック化を⾏うüある程度まとめて条件判定を⾏うØこの際、内側４重ループが取りうる⻑さ（=64=1296)分のバッファを⽤意しておく必要がある

ü計算ループの⻑さの中央値が8だったので、BLK=32あたりで決め打ち
Øこのパラメータは最適化の余地があるが、今後の課題!$omp dodo IJKL = 1, num_xxxx, BLK
ix_blk = 0
do b = 0, BLK-1
I = index_to_i(IJKL+b)
J = index_to_j(IJKL+b)
K = index_to_k(IJKL+b)
L = index_to_l(IJKL+b)do p = 1,I_ex
do q = 1,J_ex
do r = 1,K_ex
do s = 1,L_ex
if (カットオフ)
ix_blk = ix_blk + 1
enddo
enddo
enddo
enddo
end dodo i = 1, ix_blk!積分計算小原漸化式(VRR)
end doend do!$omp end do nowait

条件判定ループ

計算ループ

15

提案②実装の性能評価• 今回はsub_ssssルーチンのみの適⽤（グラフもsub_ssssのみ）
ü81個の個別関数全てへの適⽤は今後の課題

16

提案②実装の性能評価• 今回はsub_ssssルーチンのみの適⽤（グラフもsub_ssssのみ）
ü81個の個別関数全てへの適⽤は今後の課題
全体的に1.2倍程度⾼速化

パイプライン並列が促進されたことで、浮動⼩数点演算待ち時間が減少している
インデックスの保存などを⾏なった影響
で、整数演算待ちが増えている
17

ATの適⽤に関する検討• 今回⾏ったループのブロッキングパラメータはプロセッサによって最適な値が変わりうる• 81個の個別関数でも最適値が変わりうる• ⼊⼒データによって、カットオフが⾏われる割合に変化が出ることがわかっているü6-31G*基底ではカットオフされるのが68%程度に対し、cc-pVDZでは82%程度カットオフされ、実際に実⾏されるのは2割以下となる• これらの組み合わせを⼈⼿で検証するのは困難→ AT適⽤の効果があると考えられ、⾃動チューニング⾔語であるppOpen-ATによる⾃動チューニングを今後⾏う予定である
18

まとめ• 分⼦軌道計算プログラムABINIT-MPに対し、性能分析を⾏った• ループ構造が複雑であること、計算内容の異なる81種類のルーチンを呼び分けること、分岐後のループ⻑が短いこと、最内ループにカットオフ分岐があることから、ベクトル化が課題であることが判明した
üベクトル化向けのコードであるObara-vectorでも、実際にはほとんどベクトル化されていなかった• ベクトル化を念頭におき、実⾏順序のソート⼿法、ループ⻑の確保⼿法、ブロッキングによるベクトル化促進を⾏った
üSub_ssssルーチンに適⽤した結果、1.2倍程度の⾼速化を確認した

19

今後の課題• 今回適⽤した実⾏順序の変更⼿法は、GPUなどのプロセッサでも効果的であると考えられ、GPU化が今後の課題である• 最適化の適⽤はsub_ssssのみであり、他の80種のルーチンに対する適⽤も今後の課題である• ブロッキングパラメータのチューニング等を81種類のルーチンや複数のプロセッサに対し適⽤することは⼈⼿は困難であり、今後ppOpen-ATを⽤いた⾃動最適化が課題である

20

21

提案①実装の性能評価
A64FX及びIntel Xeon (CascadeLake)における実⾏時間(秒)5回実⾏した際の中央値

6-31G

6-31G*

cc-pVDZ

Obara

2.1

5.4

12.6

Obara-vector

5.9

15.8

45.1

提案手法(括弧内はソートの内時間）

2.5 (0.2)

5.8 (0.3)

17.1 (0.6)

Obara

1.2

3.1

7.6

Obara-vector

3.211.637.0

データセットA64FX (2.2GHz)

OpenMP 48コア使用
3,379.2 GFLOPS
富士通コンパイラ4.9.0–Kopenmp,fast,nosimd,ocl –Nlst=tIntel Xeon Platinum
8280 (CascadeLake)

OpenMP 1ソケット28コア使用
2,419.2 GFLOPSIntel compiler 2021.7.1 –O3 –fopenmp• いずれもアダマンタン（C10H16）分⼦に対する計算提案手法(括弧内はソートの内時間）• 基底関数が6-31G, 6-31G*, cc-pVDZで異なる
ü 6-31Gではd関数がないので、呼び出されるサブルーチンは24=16パターンのみ
ü 6-31G*, cc-pVDZでは34=81通りのルーチンが呼び出されるü cc-pVDZは⽔分⼦にもp型の分極関数を加えたもので、6-31G系より短縮⻑が⻑く、カットオフが多い• 提案⼿法はObaraと⽐較すると少し遅いが、Obara-vectorより速いü ソートは1コアのみで実⾏している。並列化により⾼速化の⾒込みはあるü 少しのオーバーヘッドでベクトル⻑を稼げたことが成果
Ø ただしこの段階では全くベクトル化はされていない
22