斉次更新の実装 - 章 GPGPU 17 - JAIST Repository https://dspace.jaist.ac.jp/

第 3 章 GPGPU 17

4.2 斉次更新の実装

(

r!₁,⋯,r!_i′,⋯,r!_N

)

(

r!₁,⋯,r!_i,⋯,r!_N

)

CPU

GPU

Update Accept / Reject Next update r!₁,⋯,!

r_i,⋯,! r_N

( )

! ↓

r₁,⋯,r!_i′,⋯,! r_N

( )

r!₁,⋯,! r_i+1,⋯,!

r_N

( )

! ↓

r₁,⋯,r!_i+1′ ,⋯,! r_N

( )

^…^{Repeat N}_times

! ′

r

! ′

r

_i+1

φ₁( )!r_i′ ⁼ ^a1s s=1

∑

64 ^Θ^s^{( )}^r^!ⁱ^′

⋮ φL

!′ r_i

( )⁼ ^a^Ls

s=1

∑

64 ^Θ^s^{( )}^r^!ⁱ^′















φ₁( )!r_i+1′ ⁼ ^a1s s=1

∑

64 ^Θ^s^{( )}^!^rⁱ⁺¹^′

⋮ φL

!′ r_i+1

( )⁼ ^a^Ls

s=1

∑

64 ^Θ^s^{( )}^r^!ⁱ⁺¹^′















Fig. 4.2: 逐次更新アルゴリズム

単純実装（逐次更新の換装）段階における計算の流れを Fig.4.2に示す．任意の電子位置をCPU上で更新後(!ri →!r_i^$), 試行関数を構成するφj(!r_i^$) の再評価のために，それを得るのに必要なデータをGPUに転送し，φj(!r^$_i)をGPU上で構築する．構築されたφj(!r_i^$) はCPUへ戻し，試行関数を構成後，更新された電子位置の棄却/採択を行い，また次の電子位置を更新する，という逐次更新の流れになる．上述の流れにおいて，!riと !ri+1の更新はランダムで行なわれるため無関係である．したがって，!r₁から!r_N までの全電子の軌道関数φj(!r^$₁)からφj(!r^$_N) をGPU上で一度に処理する構造（斉次更新）に変更できれば，GPGPU並列数を逐次更新の電子数N 倍にすることが可能である．次節ではこの斉次更新の実装について述べる．

(

r!₁′,⋯,r!_i,⋯,r!_N

)

(

r!₁,⋯,r!_i,⋯,r!_N

)

⋮

(

r!₁,⋯,r!_i,⋯,r!_N′

)

(

r!₁,⋯,r!_i,⋯,r!_N

)

CPU

GPU

Update all

Accept / Reject

!′

r₁,⋯,r!_i,⋯,r!_N

( )

! ⋮

r₁,⋯,r!_i′,⋯,r!_N

( )

! ⋮ r₁,⋯,!

r_i,⋯,r!_N′

( )

! ′ r

φ1

!′ r₁

( )= a_1s

s=1

∑

64 ^Θ^s^{( )}^!^r¹^′

⋮ φL

!′ r₁

( )⁼ ^aLs s=1

∑

64 ^Θ^s^{( )}^r^!¹^′

⋯ φ1

!′ r_N

( )

⁼ ^a^1s

s=1

∑

64 ^Θ^s

^{( )}

^r^!^N^′

⋮ φL

!′ r_N

( )

⁼ ^a^Ls

s=1

∑

64 ^Θ^s

^{( )}

^!^r^N^′















Fig. 4.3: 斉次更新アルゴリズム

進めるために，まずは斉次更新のCPU処理版を構築し，統計誤差範囲での一致を確認の上，斉次更新CPU版を段階的にGPUに換装した．

4.2.2 _{斉次更新の} GPU _換装

次に斉次更新処理を以下のようにGPU上に換装した：この段階では，§4.1で行なったのと同様に，相異なる格子点sでの積算をスレッドにアサインする（スレッド内の演算は唯一つの積算のみ）．斉次更新への変更に伴い，軌道のインデックスjに加え，電子のインデックスiに関する並列性が加わる．これを(j, i)という二次元構造のブロックで処理する（→Fig.4.4）．

結果がTable 5.2に示されている．スレッド数は64で，逐次更新の場合と変わらないが

ブロック数が（電子数N×軌道数L）に増加したことで，並列度が大幅に向上し，216電子の場合，6.39倍，1536電子の場合には5.61倍の高速化が達成された．

φj

!′ r₁

( )

⁼^a^j1^Θ¹

( )

^r^!¹^′ ⁺^a^j2^Θ²

( )

^r^!¹^′^+,⋯,^+a^j64^Θ⁶⁴

( )

^r^!¹^′

φk

!′ r₁

( )

⁼^a^k1^Θ¹

( )

^r^!¹^′ ⁺^a^k2^Θ²

( )

^r^!¹^′^+,^⋯^,^+a^k64^Θ⁶⁴

( )

^r^!¹^′

各スレッドで処理

ブロック

…

に分割

φ_L

( )

r^!₁′ ⁼^aL1Θ1

!′ r₁

( )

⁺^aL2Θ2

!′ r₁

( )

^+,⋯,^+aL64Θ64

!′ r₁

( )

φj

!′ r_N

( )

⁼^a^j1^Θ¹

( )

^!^r^N^′ ⁺^a^j2^Θ²

( )

^r^!^N^′ ^+,⋯,^+a^j64^Θ⁶⁴

( )

^r^!^N^′

φk

!′ r_N

( )

⁼^a^k1^Θ¹

( )

^r^!^N^′ ⁺^a^k2^Θ

…

( )

^!^r^N^′ ^+,^⋯^,^+a^k64^Θ⁶⁴

( )

^r^!^N^′

φL

!′ r_N

( )

⁼^a^L1^Θ¹

( )

^!^r^N^′ ⁺^a^L2^Θ²

( )

^r^!^N^′ ^+,⋯,^+a^L64^Θ⁶⁴

( )

^r^!^N^′

Fig. 4.4: 斉次更新におけるブロックとスレッドの割り当て

cavc ( j, x, y, z, spin)! M×M×M! 4×4×4=64

a_js

Fig. 4.5: cavcとajsの関係

4.2.3 スレッディング構造の改変と最適化

前節までの逐次更新および斉次更新では，積和算

"64 s=1

ajsΘs(!ri) (4.2)

において，各スレッドでa_jsΘ_s(!r_i)を計算し，スレッド間のリダクションにより総和が計算された．係数ajsのデータ構造に着目し，ランダムアクセスを解消するようなスレッディング構造の組み換えを試みたところ，大きな性能向上を達成することが出来た：係数ajs

はcavcという配列名で，cavc(j, x, y, z, spin)の5次元配列に格納されている．jは軌道のインデックス，x, y, zは格子点のインデックス，spinは電子スピンに関するインデックスである．本研究では非偏極系を取り扱うため，スピンのインデックスは退化・消失する．

x, y, zに関する要素数はM ×M ×M で，216電子系の場合M = 60，1536電子系では

M = 50である．

グローバルメモリ上に置かれたcavc(j, x, y, z, spin)に対し，各電子位置!r_iを担当するスレッドは，!ri近傍の4×4×4 = 64点の格子点でのajsの値にアクセスする（→Fig.4.5）．

前節までの逐次更新，および，斉次更新の実装では相異なるsがスレッドにアサインされているため，sを足とする連続メモリアクセスが発生するが，cavc配列はメモリ構造上，

jを足にとる１次元配列である．したがって，cavcへのアクセスはメモリアドレス上，飛び飛びのランダムアクセスとなってしまう．この場合，§3.4で述べたとおり，レイテンシの大きいグローバルメモリへのアクセスが大量に発生し，高速化を妨げる．これを避けるには，各スレッドが，軌道のインデックスjを足とした１次元配列cavcへの連続アクセスを構成するようにスレッディングを組み替え，コアレッシングを行なう必要がある．この場合，相異なる軌道jがスレッドにアサインされることになる．スレッド内の演算は，ある一つの軌道関数を構築する(4.2)式の64項の積和算となる．残されたインデックスであ

るi（電子インデックス）についてブロック化を行なった．すなわち，同一の電子インデッ

クスi を有する複数のスレッド（L本）がブロックにまとめられ処理される（→Fig.4.6）．

前節までの実装ではスレッドで並列処理されていたs=1〜64の積和算は，本実装では，

各スレッドのループ内で処理される．この構造により，従来ではスレッド間で行なっていたリダクション処理の必要がなくなった．jをスレッド並列にしたことで，以前の実装でのΘs(!ri)に関する連続アクセスが失われる．しかしながら，一つのブロック（電子インデックスiが固定）は，高々，64要素の同一のΘ_s(!r_i)を参照するため，再利用性があり，

キャッシュによる効果が期待される：あるWarpが一度，グローバルメモリ上にあるΘs(!ri) にアクセスすると，その近傍の配列は，L1やL2キャッシュに保持され，キャッシュヒットが起きやすく，Θs(!r_i)に関する連続アクセスの喪失はメモリアクセス遅延上の問題にはならないと期待される．

結果がTable 5.2に示されているが，この実装で最終的に216電子の場合，16.58倍，1536 電子の場合には30.67倍の高速化が達成された．

φ_j

( )

r!₁′ ⁼^aj1Θ₁

( )

r!₁′ ⁺^aj2Θ₂

( )

r!₁′^+,^⋯^,⁺^aj64Θ₆₄

( )

r!₁′

φ_k

( )

!r₁′ ⁼^ak1Θ₁

( )

r!₁′⁺^ak2Θ₂

( )

!r₁′^+,⋯,^+ak64Θ₆₄

( )

r!₁′ 各スレッドで処理

ブロックに分割

…

φ_L

( )

r^!₁′ ⁼^aL1Θ₁

( )

!r₁′ ⁺^aL2Θ₂

( )

r!₁′^+,⋯,^+aL64Θ₆₄

( )

r!₁′

φj

!′ r_N

( )

⁼^aj1Θ₁

( )

r!_N′ ⁺^aj2Θ₂

( )

r!_N′ ^+,^⋯^,^+aj64Θ₆₄

( )

r!_N′ φ_k

( )

r!_N′ ⁼^a^k1^Θ¹

( )

^!^r^N^′ ⁺^a^k2^Θ

…

( )

^r^!^N^′ ^+,⋯,^+a^k64^Θ⁶⁴

( )

^r^!^N^′

φL

!′ r_N

( )

⁼^aL1Θ₁

( )

r!_N′ ⁺^aL2Θ₂

( )

r!_N′ ^+,⋯,^+aL64Θ₆₄

( )

!r_N′

Fig. 4.6: 最適化によるブロックとスレッドの割り当て

ドキュメント内 JAIST Repository https://dspace.jaist.ac.jp/ (ページ 36-41)

斉次更新の実装

第 3 章 GPGPU 17

4.2 斉次更新の実装

(

)

(

)

( )

( )

( )

( )

! ′

r

! ′

r

∑

∑

∑

∑

(

)

(

)

(

)

(

)

( )

( )

( )

! ′ r

∑

∑

( )

∑

( )

( )

∑

( )

4.2.2 斉次更新の GPU 換装

( )

( )

( )

( )

( )

( )

( )

( )

…

( )

( )

( )

( )

( )

( )

( )

( )

( )

( )

…

( )

( )

( )

( )

( )

( )

4.2.3 スレッディング構造の改変と最適化

( )

( )

( )

( )

( )

( )

( )

( )

…

( )

( )

( )

( )

^{( )}

^{( )}

4.2.2 _{斉次更新の} GPU _換装