日本ソフトウェア科学会第 33 回大会 (2016 年度) 講演論文集

一般化 N 体問題アルゴリズムに対する局所性向上技法
佐藤 重幸

田浦 健次朗局所性向上は，ループ最適化の文脈で長く研究されている．多重ループによって添字アクセスによる規則的な配列走査を行う計算の場合，局所性向上は比較的容易である．一方，木データ構造を不規則に走査する計算を含む場合，局所性向上は容易ではない．近年，N 体問題アルゴリズム Barnes-Hut 法に代表される，木データ構造への独立な走査を反復するループに対する局所性改善技法 point blocking が開発された．しかし，Appel 法に代表される，より高速な N 体問題アルゴリズムでは，Barnes-Hut 法とモデル上で同等の計算が，木データ構造の相互再帰走査によって実装される．Point blocking は，この種の計算を扱わない．Appel 法の計算構造は，N 体問題アルゴリズムの一般形としての側面がある．したがって，その構造に対して機械的な局所性向上技法が開発できれば，N 体問題ア
ルゴリズム一般に対して潜在的に適用可能となる．本研究では，Appel 法の構造を対象として，機械的に適用可能な局所性向上技法を提案する．本研究は進行中である．

上技法は，ループブロッキングとループインターチェ

1 導入

ンジを組み合わせたループタイリングである．ただ参照局所性（以下，単に局所性）は，プロセッサ

し，タイリングは，任意のループに適用できるものでアーキテクチャが高速化のために仮定する計算の性質

はなく，実際上は，密行列積などの典型的な行列演算である．例えば，時間的局所性を仮定した高速化技法に適用範囲が限定される．一方，対象の計算を限定すがキャッシュであり，空間的局所性を仮定した技法がることで，より効果的なタイリングが行えることも知キャッシュのプリフェッチである．プロセッサが期待られている [8] [13] [12] ．これらの観測から，局所性向

する局所性が高いプログラムは，それらの高速化技法

上は，対象の問題や計算の構造に依存するものであの恩恵を最大限に受け，高い性能を達成する．それゆり，対象をある程度限定した技法の開発が有望である

えに，コンパイラ最適化，とりわけループ最適化の文と考える．そこで本研究では，局所性向上の対象とし脈では，局所性向上技法がよく研究されている [1] [5]．

て，N 体問題アルゴリズムに焦点を当てる．

ループ最適化における標準的な（時間的）局所性向

N 体問題を定義通りに計算する直接法は，直積をとる演算であり，2 重ループとして実現される．この˚ Locality Enhancement for Generalized N-body Al-

2 重ループには既存のループ最適化が適用できるが，gorithm.
This is an unrefereed paper. Copyrights belong to
the Author(s). This is ongoing work. If you would
like to cite this paper, please contact the first author
for auther’s version.
Shigeyuki Sato, 高知工科大学情報学群, School of Information, Kochi University of Technology.
Kenjiro Taura, 東京大学大学院情報理工学系研究科電子情報学専攻, Dept. of Information and Communication Engineering, Graduate School of Information
Science and Technology, University of Tokyo.

現実的には直接法が単純に使われることはない．空間中の粒子（相互作用を持つ体）の集合を再帰的に分割した空間分割木（3 次元空間では典型的には 8 分木）を用いた，近似に基づく高速解法の利用が一般的である．典型的な高速解法は，Barnes-Hut 法 [3] である．これは，空間分割木に対して 1 つの粒子への寄与を計算するクエリを反復するアルゴリズムである．

Barnes-Hut 法は，OpN log N q 時間アルゴリズムであるが，より高速な OpN q 時間アルゴリズムが存在す
る．その典型例は，Appel 法 [2] である．これは，複数の粒子への寄与を，空間分割木上で階層的に計算

def direct ( particles ):for pi in particles :するアルゴリズムである．直接法，Barnes-Hut 法，

for pj in particles :Appel 法は，見た目こそ異なるものの，本質的にはど

kernel ( pi , pj )れも直積をとる演算を行っている．とりわけ Appel

図1

直接法．

法は，容易に直接法や Barnes-Hut 法へ切り替えることができるため，その計算の構造が，N 体問題アルゴリズムとしての一般性を持っていると見做せる．

N 体問題ソルバは効率化は古くから広く研究されており，局所性に関しては分散メモリ環境向けの効率化

def bh ( t ):for p in t . particles ():bh_q (p , t )（すなわち通信最適化）がよく研究されている [14] [7]．
ただし，大抵はソルバのアルゴリズムないし実装のdef bh_q ( pi , tj ):改善という形式であり，ループタイリングのような機

if faraway ( pi , tj ):械的なプログラム変換の形式ではない．その数少な

kernel ( pi , tj )い例外は，Jo と Kulkarni によって提案された point

elif tj . is_leaf ():blocking [10] である．これは，ループタイリングをfor pj in tj . particles ():kernel ( pi , pj )Barnes-Hut 法の構造に適応させたものであり，その構造に強く依存している．したがって，Appel 法の構造への適用は困難であり，N 体問題アルゴリズム一般

else :
for cj in tj . children ():bh_q ( pi , cj )に対する局所性向上技法には至っていない．
そこで本研究では，Appel 法の構造に適用可能で，

図2

Barnes-Hut 法．

潜在的に N 体問題アルゴリズム一般を対象とする局
所性向上技法を提案する．これは，point blocking のアイデアを拡張したものであり，point blocking 自体のある種の一般化として見做すこともできる．

ルゴリズムである．
次に，Barnes-Hut 法 [3] を図 2 に示す．これは，粒子集合を空間分割木，典型的には d 次元空間では 2d的なアイデアを述べる．そして，現状得ている予備的

分木を用いて保持する．ここでは，各葉ノードが少数の粒子を集合として直接保持する †1 ．この入力木な結果と今後の展望を述べて，本稿の結びとする．

（t）に対して，粒子 1 つ（p）への寄与を計算するク

本研究は進行中である．本稿では，提案手法の基本

エリ（bh_q(p,t)）を反復するのが，アルゴリズムの2 N 体問題アルゴリズム概観である（bh）．寄与を計算するクエリにおいて，

本節では，N 体問題アルゴリズムの代表例である，

被作用側の粒子（pi）と作用側の木（tj）が表す部分

直接法，Barnes-Hut 法，Appel 法を簡単に紹介する．

空間が十分に離れている（すなわち faraway(pi,tj)

アルゴリズムは Python の構文に従って記述する．

が真）ならば，その部分空間を 1 つの体と近似して寄まず，直接法を図 1 に示す．全粒子について力を受与を計算する．近似できない場合は，木を降下して，ける側（被作用側）pi と力を与える側（作用側）pj

十分離れている部分空間を探す．葉ノードに到達したを列挙して，その力の寄与を計算して作用側 pi を更ら，それが保持する粒子について直接法を適用する．

新する．ここでは，計算カーネル kernel は，第 1 引
数を書き換えるものとし，pi と pj が等しいときには，相互作用は定義されない．自明に OpN 2 q 時間ア†1 元論文 [3] では葉ノードが保持する粒子数は 1 だが，
今日ではそれをチューニングパラメタとするのが標準的である．

全ての粒子同士は十分離れていると判定されることを

1
2

def appel ( t ):34

仮定している．8 行目では，コレクション c の k-組合

if t . is_leaf ():わせを生成する関数 comb(c,k) を利用している．こpassの走査によって十分離れた組を全て見つける計算は，
OpN q 時間で行える †2 ．

else

5

cs = t . children ()6

for c in cs :純な変更で対応できる．1 つは，appel で入力木の再7

appel ( c )帰的な分解で葉ノードに到達したとき（3 行目）に，葉ノードが複数の粒子を保持する場合は，2 つの単8for ci , cj in comb ( cs ,2):

直接法を適用すること．もう 1 つは，appel_dual で，

9appel_dual ( ci , cj )作用側と被作用側の両方を走査するときに，片方が10葉ノードに到達したら，それ以降は Barnes-Hut 法

11

def appel_dual ( ti , tj ):と同様の計算を適用することである．ただし，その12

if faraway ( ti , tj ):Barnes-Hut 法では，作用側と被作用側の両方を更新

13

kernel ( ti , tj )14kernel ( tj , ti )15else :

するものとする．

Appel 法は，このように直接法や Barnes-Hut 法へと容易に切り替わることができる枠組みになっており，

16

if tj < ti : tS , tL = tj , ti

それらよりも一般的な計算の構造である．実際，木の17

else :tS , tL = ti , tj

分解戦略を抽象化するだけで，直接法や Barnes-Hut

18

for c in tL . children ():法を包含した構造になる．被作用側を粒子になるまで

19

appel_dual (c , tS )分解してから，作用側を分解すると，Barnes-Hut 法

図3

Appel 法．

と同等になる．作用側と被作用側の両方（つまり入力木）を同時に粒子になるまで分解するのが，直接法に

粒子 1 つのクエリは，木の高さの定数倍程度のノー相当する．したがって，Appel 法に適用可能であり，ドしか訪問せず，木は通常 Oplg N q の高さになるよ

かつ「部分空間として大きい方を分解する」戦略に依うに構築するので，1 つの粒子への寄与の計算には存しない技法であれば，N 体問題アルゴリズム一般

Oplg N q 時間掛かる．それを N 回反復するので，全に対して，広く適用可能であることが期待できる．本体として OpN lg N q 時間アルゴリズムである．

論文では，分割戦略に依存せずに Appel 法に適用可

最後に，Appel 法 [2] を図 3 に示す．これは，Barnes-

Hut 法と同様に，粒子集合を空間分割木の形で保持能であることを以って，提案する局所性向上技法の一般性を主張する．

し，遠方の寄与を近似する．しかし，Barnes-Hut ア
ルゴリズムとは違って，全体として入力木をトップダ

3 既存の局所性向上技法

ウンに分解しながら（appel），部分木の 2 つを作用N 体問題アルゴリズムを念頭において，既存の局側ないし被作用側と見做して，両方の木を走査して寄所性向上技法を紹介する．

与を計算する（appel_dual）．これは，十分離れた 2
つの部分空間を，それぞれ 1 つの体に近似して，寄3. 1 ループタイリング

与をまとめて計算することを意味する．したがって，直接法は，同一要素へのアクセスが重複する 2 重探索において被作用側と作用側は同様に扱われ，部分

ループである．このときループタイリングは効果的空間として大きい方を分解する（16–19 行目）というに機能する．図 4 にループタイリングの適用結果を戦略が取られる．ここでは，記述の簡潔さのために，†2 厳密には fair split という条件で空間分割がなされている前提を要する．詳細は [6] を参照．

入力の空間分割木は粒子 1 つになるまで分割を行い，

def tiled_direct ( particles ):def bh ( t ):n = len ( particles )ps = t . particles ()# range(0,n,k) ” [0,k,2*k,...,n-1]

for i in range (0 , len ( ps ) , k ):bh_block_q ( ps [ i : i + k ] , t )for i in range (0 ,n , k ):for j in range (0 ,n , k ):for pi in particles [ i : i + k ]:def bh_block_q ( pi_block , tj ):new_b = []for pj in particles [ j : j + k ]:for pi in pi_block :kernel ( pi , pj )if faraway ( pi , tj ):図4

ループタイリングを適用した直接法．

kernel ( pi , tj )else :

示す．

new_b . append ( pi )被作用側と作用側を表す添字 i と j が，定数 k 分

if tj . is_leaf ():だけブロック化され，それらに対応する粒子ブロッfor pi in new_b :ク particles[i:i+k] と particles[j:j+k] のルー

for pj in tj . particles ():プが導入されている．そして i についての粒子ブロッkernel ( pi , pj )ク particles[i:i+k] のループが，j のループの内else :

側に入っている．すなわち，ループネストがインターチェンジしている．このとき，2k 個分の粒子がキャッfor cj in tj . children ():シュに載りきるならば，内側の 2 重ループは，局所性

bh_block_q ( new_b , cj )が高い．なぜなら，一度粒子ブロックがキャッシュに

図5

Point blocking を適用した Barnes-Hut 法．

載りきった後は，それが 2 重ループ内で再利用され続けるからである．

セスの重複が生じる計算なので，1 つ目の前提条件はこのループタイリングという変換には，時間的局所満たす．したがって，N 体問題アルゴリズムの時間的性向上の基本的なアイデアが凝縮されている．重要な局所性向上での主な課題は，2 つ目と 3 つ目の実現でのは次の 3 点である．

ある．

‚ 同一データ要素へのアクセスが重複する計算を対象とする．

‚ キャッシュに載りきるサイズにデータを分割し，それに対する部分計算を導入する．

3. 2 Point Blocking
端的に言えば，Barnes-Hut 法とは，直接法における内側ループが，木走査の再帰関数に置き換わった‚ 分割されたデータが再利用されるように制御をものである．この観察に基づいてループタイリング

変えて，部分計算をスケジューリングする．

を Barnes-Hut 法に適応させたのが，point blocking

1 つ目は，時間的局所性向上の前提である．同一デーである．具体的には，外側の粒子ループをブロック化

タ要素を一度しかアクセスしない計算（所謂ストリー

し，粒子ブロックのループを木走査の内側に入れるとム処理）では，キャッシュされたデータの再利用が生いう変換である（図 5）．

じる余地はない．つまり，時間的局所性は向上し得な

Point blocking のアイデアは単純だが，これが局

い．2 つ目はループブロッキングに，3 つ目はループ所性向上を達成するためには，重要な仮定を 2 つ要

インターチェンジに，それぞれ対応する．N 体問題

する．

の定義から，N 体問題アルゴリズムはどれも，本質的には直積を取る演算を行う．直積演算は，要素アク‚ 同一粒子ブロック内の粒子が，空間モデルの上で近傍にある．

‚ 粒子 1 つのクエリでアクセスされるデータが，キャッシュに載りきらない．

性が悪化する事態が生じる．つまり，point blocking
は，ブロック導入の空間オーバーヘッド分だけ，局所

第 1 の仮定は，粒子集合 ps が，空間インデクスに性へのオーバーヘッドを潜在的に伴っている．この実沿って整列されることで，自然と達成される．この粒装技法は，point blokcing の空間オーバーヘッドによ子集合の整列自体は，Barnes-Hut 法では常識的な処理である †3 ．この仮定が満たされることで，同一ブる局所性への影響を小さく抑える効果がある．

ロック内の粒子が類似した木走査を行うことになり，関する 1 つの仮定を導入することで，point blocking

その結果，pi_block が k 個に近い要素数を保持したは，元論文 [10] で示されるように，局所性向上技法とまま，bh_block_q による木走査が行われる．これは，して成功裡に機能する．

上記のように，入力に関する 2 つの仮定と実装に

ループタイリングが適用された直接法における内側の 2 重ループと，時間的局所性の観点で類似する．

4 提案技法

第 2 の仮定は，入力木が十分に大きいことを仮定す

本節では，point blocking のアイデアを拡張した時ることで達成される．この仮定は，point blocking の間的局所性向上技法を提案する．まず，point blocking

利得に強く関わる．第 1 の仮定により，粒子集合の上が Appel 法に適用しがたい原因をまとめ，次にそれをのループにおいて隣接する粒子は，類似した木ノード

解決するアイデアと提案技法の適用例の概観を示す．

へのアクセスを行うことが期待される．したがって，
第 2 の仮定が満たされてしまうと，隣接する粒子の4. 1 Appel 法の難しさ

クエリによって，キャッシュ上のデータの再利用が自Point blocking は単純でありながら，入力と実装技

然に生じる．つまり，元の計算の構造において時間的術に関する仮定を導入することで上手く機能すること

局所性が十分向上してしまう．第 2 の仮定が満たさは，前節で述べた．この point blocking の単純さは，れないと，計算の構造を複雑化させる point blocking

外側の粒子集合の上のループに依るところが大きい．

のオーバーヘッドだけが影響して，時間的局所性向上

まず，ループタイリングの既存研究から，外側のルーの利得が生じない．

プについてブロッキングして，ブロックのループを内以上の時間的局所性の本質に関わる仮定が満たさ側にインターチェンジすれば局所性が向上することがれた上で，高い実験的性能を達成するためには，実装即座に予想できる．そして，外側のループのブロッキ技法において重要な工夫を 1 点要する．それは，ングは自明であり，インターチェンジも内側の再帰関‚ 粒子ブロック用の領域を，木走査の途中で新た数の本体にループを導入するだけで済む．これに加えに確保せずに，入力木の深さ毎にブロック 1 つて，内側の再帰関数の振る舞いが単純であることも重分の領域を確保して再利用する要である．具体的には，引数の粒子だけに依存して，ことである．これは，入力木の高さが小さい定数であ単一の入力木を線型に走査する点である．これによっり，木走査が深さ優先探索であることを仮定した実装て，類似した木走査（すなわち重複したノードアクセ

技法である．

ス）をスケジュールする方法が単純になる．そして類

この実装技法が重要なのは，粒子ブロックの表現の似した木走査が連続することを仮定できるなら，走査空間コストが定数ではないことに由来する．新たに作中のノードの深さに注目することで，配列を添字アクられるブロックのために確保される領域が，キャッシュ

セスするループのように，木走査をモデル化できる．

に載っていないデータ領域であった場合，時間的局所

前述の point blocking の実装技術は，このモデル化を暗黙に利用している．

†3 ただし，実際に整列されているかを機械的に判定するのは困難である．そのため，Jo と Kulkarni は，point
blocking に合わせて，クエリの実行途中に粒子集合を整列させる手法も開発している [11]．

上記のような Barnes-Hut 法の構造上の単純さを，Appel 法の構造は持たない．それゆえに，下の点が非自明になる．

‚ 木走査における重複アクセスのモデル化．
‚ ブロッキングの対象．
‚ 計算のスケジューリング．
これはそれぞれ，3. 1 節で述べた時間的局所性向上の基本的アイデアの 3 点に対応する．
前述のように，3. 1 節で挙げた第 1 の前提条件は，‚ 再帰関数が木を κ 段分（κ は小さい定数）進める計算には重複アクセスを含む．

‚ 再帰関数の継続集合をブロッキングする．
‚ 木ペアの深さから，継続の深さを定義して，継続の制限付き幅優先スケジューリングを行う．
以下では，それぞれの意味を順に説明する．

N 体問題アルゴリズムでは自然と満たされる．Appel

第 1 の解の背後にある観察は，Appel 法の再帰関法でも，全体としては重複アクセスが生じている．第数の内側に直積演算が存在する（9 行目と 19 行目）

1 の問題点は，どのように重複アクセスが生じるかをことと，その再帰呼出し毎に木を 1 段階分解してい把握しにくい点である．実行時の完全な振る舞いは，ることである．これらの観察から，あるノード（もしプログラム変換時には利用できない．機械的なプログ

くはノードのペア）から深さ κ までノード上の再帰

ラム変換として局所性向上を実現するためには，対象関数では，近い親族同士（例えば兄弟や従兄弟）の直の計算をモデル化して，計算を静的に予測可能にし積が取られるので，そのノードアクセスが重複するとなければならない．実際，point blocking では，類似仮定できる．

した木走査の連続であるとモデル化することで，2 重

第 2 の解は，第 1 の解から自然に導出される．再帰

ループと同列の計算と見做している．Appel 法では，関数の再帰を，木の深さ κ に至る所で止めると，深同様のモデル化は適用できない．

さ κ ` 1 以降について再帰関数の継続が複数生成さ第 2 の問題点は，計算の全体構造が単純なループをれる．これを再帰的に適用すれば，継続集合から（よ成していないことに起因する．前述のように，ループり大きい）継続集合が構成される．κ が小さければ，が最外に存在すれば，point blockig のように，単純に生成される継続も少なく，その継続集合から，次の継それをブロッキングの対象にできる．しかし，Appel

続集合を作る計算の中で同じデータがアクセスする法における最外に存在するのは，入力木の再帰的分解

可能性が高い．しかし，継続が増えすぎると，アクセである．これを 3. 1 節で挙げた第 2 の点が満たされスされるデータが散らばり，キャッシュ内に残っているようにブロッキングする方法は，自明ではない．

る可能性が低下する．したがって，継続集合のブロッ

第 3 の問題点は，第 2 の問題点と密接に関係する．
ブロッキングする方法が非自明であるということは，

キングとは，計算対象の継続集合をある程度小さく抑えることを意味する．

ブロックに対する部分計算も非自明になり，それを第 3 の解は，第 2 の解に対して詳細な実装を与えスケジュールする方法は，更に非自明になる．Point

るものである．Appel 法の再帰関数が生成する継続

blocking の実装技法から分かるように，局所性向上

は，木のペアで表現できる．一方の木だけを基準にし技法の実装上の問題は，最終的には，スケジューリンてスケジューリングするのは，木の分解戦略に依存しグの実装技法へと集約される．スケジューリングが非自明ということは，実装技法も非自明であることを即た技法となり，望ましくない．そこで，そのペアの入力木における深さの和 †4 を，継続の「深さ」として座に意味する．

扱う．同一の深さの継続が，入力木の上で同じ深さのAppel 法を含んだ N 体問題アルゴリズム一般を対ノードにアクセスすることは保証されない．しかし，象する局所性向上を設計する上で，上述の 3 点は，避各継続における木の分解戦略が同一であると期待でけられない問題である．

きることから，この深さの定義は，継続がアクセスするノードの深さの近似として妥当である．継続の深さ

4. 2 提案技法のアイデアと適用例
提案技法のアイデアは，前述の 3 つの問題点に次のような解を与えることに集約される．

†4 深さの平均としてもよい．しかし，和にすると，継続を進めると必ず次の深さに進むことが保証されるので，スケジューリングを理解しやすい．

を利用することで，継続の深さ優先スケジューリング
と幅優先スケジューリングが定義できるようになる．

5 現状と今後

幅優先スケジューリングは，次の深さの継続ブロック

本節では，まず現状行った実装と実験について述べを大きくする一方，深さ優先スケジューリング小さくる．そして今後の課題と展望について述べて，本稿の抑える．継続のブロッキングは，幅優先に深さ優先を結びとする．

組み合わせて，次の深さの継続ブロックに制限を加えた幅優先スケジューリングによって実現できる．

5. 1 実装と予備実験

図 6 に，提案技法を Appel 法に適用した結果を示す．

2 節で述べたように，Appel 法は，入力木の葉ノー

threshold は，継続ブロックのサイズに関する閾値をドが複数の粒子を保持する場合，直接法と Barnes-

表すパラメータである．kappa は，前述の κ に相当す

Hut 法を自然に内包する．この混成版の Appel 法を，る．appel_dual_bounded は，図 3 の appel_dual のC++14 で実装した．この参照実装を元に，図 6 と同計算を木の深さ minpbi, biq 分で打ちきって，blocks

様に，提案技法を適用した Appel 法を手で実装した．

に継続を入れる関数である．継続の深さの計算には，

手実装ではあるが，本研究の目的に沿って，Appel 法

木 ti と tj の深さ di と dj を使う．木の深さは，継の分解戦略に依存した技法は導入していない．ただ

続として木と共に保存する．appel_bounded も同様

し，機械的に実現可能な変換による効率化は施した．

である．appel_block が，継続の制限付き幅優先スケ

ここで機械的に実現可能という意味は，C++テンプ

ジューリングを実装している関数である．peek_size

レートメタプログラミングによるコンパイル時コード

は，次の深さの継続ブロックのサイズをみる関数である †5 ．深さ毎にブロックの領域を確保して再利用す

生成によって実現可能という意味である．典型例は，

るのは，point blocking と共通である．

や appel_dual_bounded）から，インライン展開を用

定数深さ分で再帰を打ちきる計算（appel_bounded

提案技法は，一見すると point blocking と全く異いて実際の再帰呼出しを除去したことである．つまなるように見える．しかし，非作用側の木を常に優り，κ に関わる計算は全て展開され，オーバーヘッド

先する分解戦略を Appel 法に適用した場合，アルゴは極小化さている．また，各種パラメータは，テンプ

リズムは Barnes-Hut 法に類似し，提案技法は point

レートパラメータの形式となっており，コンパイル時

blocking に類似する．入力木の高さを h としたとき，定数である．このような C++テンプレートメタプロ

深さ h の継続ブロックが，point blocking における粒グラミングによって，提案技法の適用したことによる子集合を分割したブロックに相当する．そして，深さ

コードの複雑化によって，コンパイラ最適化が適用さh ` 1 から 2h まに継続ブロックが，point blocking

れないという事態を極力回避している．

の再帰呼出しで利用される粒子ブロックに相当する．

この 2 つの Appel 法の実装を，2 次元空間上にラpoint blocking が入力に仮定する，同一ブロック内のンダム生成された粒子（粒子数 N “ 107 ）から作ら粒子が空間デル上で近傍にあることは，提案技法の場

れた 4 分木（ノード数 1.4 ˆ 106 ，高さ 12）を与えて合，入力木の葉に近い位置から生成さた継続ブロック

実行し，実行時間を比較した．近似度のパラメータはが，深さ h に至ることから，自然と達成される．こθ “ 0.4 とし，葉ノードが持つ最大粒子数 Nl を 20 とのことか分かるように，提案技法は，Appel 法の構造RRした．実験環境は，Intel⃝Xeon⃝CPU
E5-2699 v3の中で point blocking を潜在的に含した，より一般(2.3GHz) と DDR4-2133 768 GB を搭載したサーバ的な技法となっている．

である．このプロセッサは，45 MB Smart Cache を備えている．実験対象はどちらも逐次実装なので，こ†5 ここで，「次の深さ」の定義には選択肢がある．文字通りに深さを 1 つ先のブロックのサイズとしても良い
し，κ 分だけ先にあるブロックの総和とみてもよい．

のキャッシュがほぼ全て利用されたと考えられる．
実験した範囲での最良のパラメータは，κ “ 2 で，threshold “ 64 であった．このとき，提案技法による速度向上は，参照実装比で 4.1%であった．これは，トメタプログラミングが駆使されてライブラリが実

5 回実行した中央値である．N と Nl を半分にして，装されている．ここに提案技法を上手く組み込むこ

ほぼ同じ木のサイズで実験した場合は，12.3%にまでとは，可能であるが一仕事である．組み込んだ後の展速度向上が増えた．しかし，これは Nl が小さすぎて，望として，我々のライブラリのユーザレベル API を

Appel 法の実行時間が現実的でないほど大きくなっ

解釈することで，より効率的なスケジューリングのた．この性質は point blocking にも共通する．元論

実装に切り替えることを想定している．例えば，修文 [10] で示される実験結果は Nl “ 1 の場合であり，正 Barnes-Hut 法 [4] は，より効率的なスケジューリ

低次元の N 体問題として現実的なパラメータ設定でングを実現できる見通しがある．我々のパターンライはない．2 次元 N 体問題という設定では，提案技法

ブラリは，現状共有メモリ環境向けだが，分散メモリの速度向上を 4.1%と見做すのが妥当である．

環境に適応させる予定である．その過程で，提案技法

4.1%という結果は，Barnes-Hut 法に対する point
blocking と比べると，僅かな効果にも見える．しかの分散メモリ環境での有効性も調べることも想定している．

し，OpN log N q 時間アルゴリズムの Barnes-Hut 法と比べて，OpN q 時間アルゴリズムの Appel 法は，そ参 考 文 献もそもデータアクセスの重複が定数回である．提案技[ 1 ] Allen, R. and Kennedy, K.: Optimizing Compilers for Modern Architectures: A Dependence-BasedApproach, Morgan Kaufmann Publishers, 2001.
[ 2 ] Appel, A. W.: An Eﬃcient Program for ManyBody Simulation, SIAM J. ScI. Stat. Comput.,Vol. 6, No. 1(1985), pp. 85–103.
[ 3 ] Barnes, J. and Hut, P.:
A hierarchical
OpN log N q force-calculation algorithm, Nature,
Vol. 324(1986), pp. 446–449.
[ 4 ] Barnes, J. E.: A Modified Tree Code: Don’t
Laugh; It Runs, J. Comput. Phys., Vol. 87,
No. 1(1990), pp. 161–170.
[ 5 ] Bondhugula, U., Hartono, A., Ramanujam, J.,
and Sadayappan, P.: A Practical Automatic Polyhedral Program Optimization System, Proc. PLDI
’08, ACM, 2008, pp. 101–113.
[ 6 ] Callahan, P. B. and Kosaraju, S. R.: A Decomposition of Multidimensional Point Sets with Applications to k-Nearest-Neighbors and n-Body Potential Fields, J. ACM, Vol. 42, No. 1(1995), pp. 67–90.
[ 7 ] Driscoll, M. B., Georganas, E., Koanantakool, P., Solomonik, E., and Yelick, K. A.: A
Communication-Optimal N-Body Algorithm for Direct Interactions, Proc. IPDPS ’13, IEEE, 2013,
pp. 1075–1084.
[ 8 ] Franchetti, F., Voronenko, Y., and Püschel, M.:
Formal Loop Merging for Signal Transforms, Proc.
PLDI ’05, ACM, 2005, pp. 315–326.
[ 9 ] Gray, A. G. and Moore, A. W.: ‘N-Body’ Problems in Statistical Learning, Proc. NIPS ’00, 2000,
pp. 521–527.
[10] Jo, Y. and Kulkarni, M.: Enhancing Locality for
Recursive Traversals of Recursive Structures, Proc.
OOPSLA ’11, ACM, 2011, pp. 463–482.
[11] Jo, Y. and Kulkarni, M.: Automatically Enhancing Locality for Tree Traversals with Traversal Splicing, Proc. OOPSLA ’12, ACM, 2012, pp. 355–374.
[12] Krishnamoorthy, S., Baskaran, M., Bondhugula,法によって計算の構造が複雑化し，オーバーヘッドが導入された上でも，4.1%の有意な性能向上があった
ことは，局所性向上技法として有意義な結果であると評価している．

5. 2 今後の課題と展望
本研究の喫緊の課題としては，複数の例について提案技法を適用して，その効果を調べることである．
（低次の）k 次元空間については，現状の実装をそのまま再利用して実験が可能である．Appel 法と同様の計算の構造を持つが，中間ノード（部分空間）におけるデータと計算量が大きい FMM というアルゴリ
ズムについても提案技法を適用し，その効果を調べる．高次元空間においても Appel 法と同等のアルゴ
リズムは存在する [9]．高次元用の空間分割木と共に，
そのアルゴリズムを実装し，提案技法を適用して効果を調べる．高次元の空間分割木は，大抵は 2 分木であり，葉ノードの粒子数が少ないので，提案技法の効果が大きいと期待される．
我々は，一般 N 体問題アルゴリズムパターン cosplit
を C++14 ライブラリとして開発している [15]．これは，2 節で述べたように，木の分解戦略を抽象化する
ことで，多様な N 体問題アルゴリズムの記述を可能にしている．総称性，モジュール性，低オーバーヘッ
ドの 3 つを実現するために，複雑な C++テンプレー

U., Ramanujam, J., Rountev, A., and Sadayappan,
P.: Eﬀective Automatic Parallelization of Stencil
Computations, Proceedings of the 2007 ACM SIGPLAN Conference on Programming Language Design and Implementation (PLDI ’07), ACM, 2007,
pp. 235–244.
[13] Sahoo, S. K., Krishnamoorthy, S., Panuganti,R., and Sadayappan, P.: Integrated Loop Optimizations for Data Locality Enhancement of Tensor

Contraction Expressions, Proc. SC ’05, IEEE, 2005,
pp. 13.
[14] Warren, M. S. and Salmon, J. K.: AstrophysicalN-body Simulations Using Hierarchical Tree Data
Structures, Proc. SC ’92, IEEE, 1992, pp. 570–576.
[15] 佐藤重幸, 田浦健次朗: N 体問題のためのデータ並
列スケルトンの設計と実装, 第 18 回プログラミングおよびプログラミング言語ワークショップ PPL 2016 予
稿集, 2016.

blocks = Nonedef appel ( t ):h = height ( t )blocks = [[] for _ in range ( h *2)]
blocks [0]. append (( t ,0 , None , None ))for d_sum in range (0 , h *2):
appel_rec ( d_sum )def appel_block ( d_sum ):while blocks [ d_sum ] != []:ti , di , tj , dj = blocks [ d_sum ]. pop ()if tj is None :appel_bounded ( ti , di , kappa )else :
appe l_ du al_ bo un de d ( ti , di , tj , dj , kappa )if peek_size ( blocks , d_sum ) >= threshold :appel_block ( d_sum )def appel_bounded (t , d , b ):if t . is_leaf ():passelif b == 0:
blocks [ d *2]. append (( t ,d , None , None ))else :
cs = t . children ()for c in cs :appel_bounded (c , d + 1 , b - 1)
for ci , cj in comb ( cs ,2):appe l_ du al_ boun de d ( ci , d + 1 , b - 1 , cj , d + 1 , b - 1)

def ap pel_d ual_bo un ded ( ti , di , bi , tj , dj , bj ):if faraway ( ti , tj ):kernel ( ti , tj )kernel ( tj , ti )elif bi == 0 or bj == 0:
blocks [ di + dj ]. append (( ti , di , tj , dj ))else :
if tj < ti : tS , dS , bS , tL , dL , bL = tj , dj , bj , ti , di , bielse :tS , dS , bS , tL , dL , bL = ti , di , bi , tj , dj , bj

for c in tL . children ():appe l_ du al_ boun de d (c , dL + 1 , bL - 1 , tS , dS , bS )図6

提案技法を適用した Appel 法．