Chiselの最新版を追いかける、あるいは、ソースをバグ修正する方法について。

ソースコードはgithubで公開されているので自由にcloneできる。ただ、sbtの仕組みがよくわからないので、ScalaやChisel自体がどこにインストールされているのか不明で、そのため、仮にChiselのソースコードを自分で修正したとしてそれを使ってプロジェクトを実行するにはどうすればよいのか、やってみるまでよくわからなかった。実際には、やってみると非常に簡単であった。

1. ソースをcloneする。必要であればgithubでforkしてそれをcloneする
2. ソースを修正する。Chselの実装は”chisel/src/scala/main”にある。
3. cloneしたソースのトップディレクトリにはMakefileが用意されていて、これで修正したソースをベースにbuildが行われdeployされる。
4. Chiselを利用しているプロジェクトのMakefileを適切に修正する。

一連の流れをコマンドで書くと：

1
2
3
4
5
6

git clone https://github.com/ucb-bar/chisel.git
cd chisel/src/mains/scala/
emacs Tester.scala
git commit -am .... (必要な場合のみ)
cd ../../..
make publish-local

buildされたJavaのclassファイルは、デフォルトでは”~.$hostname/local/edu.berkeley.cs/chisel_2.10/2.3-SNAPSHOT”にdeploy(コピー)される。”2.3-SNAPSHOT”は、執筆時の場合で、バージョンによっては違うかもしれない。あとは利用してるプロジェクトのMakefileを修正すればよい。変数”CHISEL_DEFAULT_VERSION”等が利用しているclassファイルの場所を指していて、デフォルトではこれが”latest.release”にセットされているので、これをmakeのオプションなりで上書きすればよい。

gitのコマンドについては、forkをして修正をフィードバックする予定ならばbranchを切っておく必要がある。これは別の話。負の値pokeすると動作しない問題は修正して、報告すると無事mergeされた。

Scalaの埋め込みDSL(emmbeded DSL;eDSL)でHDLを設計するChisel https://chisel.eecs.berkeley.edu/ について。

VHDLは冗長だ。色々と面倒くさい決まりごとがたくさんある。 特に困惑するのは、Nビット加算器の記述のような「パラメータ化された回路」を設計する時だった。 VHDLはAdaをベースにしているだけあり、非常に豊富な機能があり、 パラメータ化された回路を記述することなどわけない。 え、Adaって何？という困惑を除けばだ。

それを解決する方法として高位合成という技術がある。 これはC言語などのふつーのプログラミング言語のコードから、回路として実装可能なHDLを生成する技術だ。 これはこれで、我々の目的であった「パラメータ化された浮動小数点演算回路」を実現するには、 ハードルが高かった。そもそも、高位合成ツールはHDLでの実装を細かい粒度でコントロールしたいという、 我々の目的にはあまり向いてないような気がする。実際には、真面目に使うことを検討したことがない。 そこで、2006年頃、パラメータ化されたVHDL生成する手法として、 埋め込みRuby(eRuby)によるメタプログラミングを採用した。 その詳細は説明を省略し、その例だけを直接に示す：http://galaxy.u-aizu.ac.jp/note/wiki/UnconventionalHDLProgramming

要は、VHDLをべたに手書きするのは大変につらい(EmacsのVHDL modeは挙動がおかしいし)。 Verilog HDLなら無駄に冗長なところは取り除かれていると思うが、 それでも「パラメータ化」をするのは大変そうである。 これも真面目に使うことを検討したことがない。 なぜなら、Chisel があるから。論文や スライドを見てみると、 eRubyを使った我々の方法よりエレガントで、 なによりHDLシミュレータを使わずに回路記述のシミュレーションできるという点が素晴らしい。 テストベンチもこのeDSLで書くことができる。 eDSLで記述することの利点は、少ないハードウエア記述で回路を生成できるということである。 上の論文によると、3020行のVerilog HDLコードと同等のものを記述するのにこのeDSLでは1046行ですんだという。 そして、eDSLからはVerlog HDLが生成されるだけではなく、シミュレーションのためのC++コードが生成される。 この実行速度は商用シミュレータより7倍以上高速(約5400秒 vs. 約580秒)だという。 なお、FPGAによる同じ回路のエミュレーションは約80秒であったそうだ。 一方、コンパイル時間は、シミュレータ, Chisel, FPGAの順に、約20秒, 約120秒, 約3700秒である。 綺麗なトレードオフ関係。とはいえ、FPGAでのエミュレーションは面倒が多いので、実用的かどうかは疑問が残る。

Chiselの導入方法についてはそのWebサイトに説明があるのだが、いまいちはっきりとはしない。 チュートリアルとGetting Start Guideとマニュアルと。試行錯誤した結果をメモしておく。

Ubuntu 14.04LTSの場合

Chiselを利用するのに明示的にインストールが必要なものはJRE, JDK, sbt(Scala用のbuildシステムの一種。RubyのRakeみたいなもの)。 シミュレータを使う場合にはg++も必要。

Ubuntuの場合JREとJDKはシステムのパッケージをインストールすればよい。

1

sudo aptitude install openjdk-7-jre

sbtは、sbtのサイトからDebian用のパッケージをDLして、dpkgでインストールした

1
2

wget https://dl.bintray.com/sbt/debian/sbt-0.13.6.deb
sudo dpkg -i sbt-0.13.6.deb

これで準備が整う。ScalaやChiselは明示的にインストールする必要はなく、 chisel-tutorialをレポジトリをcloneし、その中のサンプル用ディレクトリでmakeすると、 sbtが必要なもののインストールなどは面倒見てくれる。

1
2
3

git clone https://github.com/ucb-bar/chisel-tutorial.git
cd chisel-tutiroial/hello
make

以上で、”hello”というサンプルが実行できるだろう。 ただし、Ubuntu 14.04の場合、”/bin/sh”が”/bin/dash”へシンボリックリンクになっており、 makeコマンドが失敗をする。これはMakefileがshellとしてbashを想定しているためで、 おそらく実害はないはずだが、”/bin/sh”を”/bin/bash”へのシンボリックリンクにする必要があった。 dashや古いzshが受け付けないコマンドオプションを使っているため。

MacOSの場合

JREとJDKはシステムのものでOK。g++もXcodeのコマンドラインツールをインストールすればいい。 sbtはMac用のuniversal package(zip)をDLして展開し、 “bin/sbt”, “bin/sbt-launch-bash”, “bin/sbt-launch.jar”をpathの通ったところに置けばよい。 詳しくはここを参照する：http://www.scala-sbt.org/0.13/tutorial/ja/Manual-Installation.html

あとはUbuntuと同じ。ただし、shellの調整は必要ない。

1
2
3

git clone https://github.com/ucb-bar/chisel-tutorial.git
cd chisel-tutiroial/hello
make

自分用のプロジェクトを作る

sbtを直接利用する場合にはこの限りではないが、tutorialと同じ方法でChiselを使うには、以下の2つのファイルが必要。 tutorialのファイルをベースに編集修正したもの。

build.sbt: https://gist.github.com/dadeba/97e5b03a49e0133bed59

Makefile: https://gist.github.com/dadeba/45631a382c6fb12a2a64

makeコマンドで”*.scala”に記述された回路がコンパイルされ、シミュレーションが実行される。

Webブラウザを使ったチャットがあると、ネットをつかった会議や議論の時に便利だろう。 SkypeやGoogle Hangoutは、チャットがbuilt-inされているので、それを併用すればよいのだけど、 いわゆる「テレビ会議」システムを使った、他拠点での会議をおこなう時に、 言葉で説明するより文字で説明した方がよいことがある。 たとえば、定量的に精密に議論をしたい時に、「テレビ会議」では誤って伝わる可能性があり、 文字を使ったコミュケーションを併用したいことがある。 それと議論のためにファイル(おもにグラフのファイル)やURLの共有をしたい時にもチャットは役に立つ。

巷にはそのためだけの商用サービスなどもあり、各所で利用されているとは思うのだが、 研究プロジェクト内でのみ共有するような議論は、外部に保存したくないなどの微妙な問題があり、 あまり使われるとは思えない。研究費からそのような費用を出すのは、不可能ではないが、 研究室レベルでやるのは簡単ではないという事情もあり。 それで、オープンソースでなんとかならないかというと、あまりよい選択肢はない。 往年の掲示板CGIを使うというのもあると思うし、Wikiベースのシステムでも、 あとで議事録を残したり、議論内のファイルを保持記録しておくにはよいかもしれない。

だからChatonというWebベースの(ほぼ)リアルタイムチャットソフトがある。

http://chaton.practical-scheme.net/

ソースコードは https://github.com/shirok/Chaton で公開されており、 自分の使えるLinux機があれば、Gaucheが動くことを前提として、手軽に設置利用できる。 ただ、一つ面倒なのは、ChatonはCOMETというプロトコルを採用していて、 そのためにhttpポートに加えて、ある特定のポートをコンテンツのポーリングのために必要とする。 githubで公開されているコードは、このポートが任意に公開できることを前提としていて、 80と443しかポートが空いていない場合には、reverse proxyを設定する必要があり、 そのためにはソースコードを一部修正する必要があった。

ので自分用に修正をしたコードをフォークした。https://github.com/dadeba/Chaton

以下、そのセットアップ方法を記録しておく。

ソースをcloneする

1

git clone git@github.com:dadeba/Chaton.git

“conf/site.conf”を適切に編集する。

最低限必要なのはgoshのパスと、COMET serverを保持するディレクトリ(下記4.参照)。

設置するチャット用の設定ファイルを作る。

Chatonオリジナルものは”sample.conf”になる。 このforkバージョンでは追加の設定が必要なのでそれについて。 COMETのポートをApacheのreverse proxyを使って外部に公開するため、 それ用のpathを(comet-proxy-path "comet1")で指定する。これはこの行を追記する。 また、(comet-port "9997")も適当なポート番号を指定する。このポート番号はあとで利用する。

Chaton用ファイルを生成

オリジナルのドキュメントにあるように

1

gosh ./build-site mychat.conf

としてChatonを設置する。この結果、server-htdocs-dirとserver-data-dirに必要なファイルが生成される。 server-htdocs-dirのほうはhttpdからアクセスできる場所ではないといけないし、 server-data-dirはhttpdが書き込みできる必要がある。一番簡単には、どちらもhttpdが動くユーザーの権限、 Ubuntuの場合”www-data:www-data”がowner:groupとなる場所を指定した方がよい。 httpdによる、server-data-dirへの適切な書き込みやアクセス権限がないと、原因不明で動作しない。 面倒を避けるために、COMET serverを保持するディレクトリもwww-dataがownerとした方がよい。

例えば”/var/www/chaton”を作成し、“`chown www-data:www-data /var/www/chaton”として、 COMET serverのpathは”/var/www/chaton/bin”として、各チャットのserver-data-dirは “/var/www/chaton/mychat/logs”のようにディレクトリを指定するのが一番簡単。 オリジナルのドキュメントに従って、各自のホームディレクトリ下に作成すると、 権限の設定が面倒なことになる場合がある。 この場合、build-siteの実行も、www-dataの権限で行うか、 rootで実行後、生成されたファイルのownerとgroupをwww-data:www-dataとする必要がある。

Apacheのreverse proxyの設定

reverse proxy moduleは読み込まれているという前提で 以下の設定をApacheの設定ファイルに追記する。

1
2
3

Alias /chaton/mychat /var/www/chaton/mychat/htdocs
ProxyPass        /comet1 http://localhost:9997
ProxyPassReverse /comet1 http://localhost:9997

最初の行はserver-htdocs-dirで指定したディレクトリと、ApacheでのURLとの対応を指定する。 残りの2行でreverse proxyを設定する。COMET serverは9997ポートで動作するので、 それを”http://example.com/comet1“のURLでアクセスできるように指定する。どちらの行も必要。

COMET serverを起動

1

sudo -u www-data gosh /var/www/chaton/bin/chaton-viewer-mychat

Apacheを再起動する。

http://example.com/chaton/mychat

にアクセスすると利用開始できると思う。

「Coders at Work」から引用するページ。なお、読み終わるまで更新される。 各インタビューごとに最低ひとつ引用する。

36ページ

私は普通C++のテンプレートが大好きというような人にはあまり近づかないことにしています。

第1章ジェイミー・ザウィンスキー(LISP、Netscape開発など)より。人外界にも派閥があるもので。

76ページ

最適化が楽しいのは必要ではないからです。

第2章ブラッド・フィッツパトリック(LiveJournal)より。目標が定量的に評価できて明確なのが最適化のよいところである。速くて正しいことは大正義なので。

104ページ

「これは私たちのやり方ですが、皆さんはほかのやり方をしてかまいません」と言うのは、人類全体に大きな損害をもたらしたし、今後もずっと続くことでしょう。

第3章ダグラス・クロッフォード(ECMAScript(Javascript)の規格を批判)より。この部分はK&RのC言語ではソースコードの整形を規定しなかったため、様々な流儀のソースコードがのさばってしまったことを指している。コーディングスタイルという本来はプログラミングでは本質的でないところに知的リソースは使わないのが正しいという。

133ページ

結局は実数を有限精度の浮動小数点数で表現することの頭がおかしくなりそうなトレードオフの話になって、大変なだけです。

第4章ブレンダン・アイク(Mozilla CTO)より。これは、数値計算には興味を持てなかったという理由について。

199ページ

なぜなら知性というのはスカラー値ではなく、ベクター値だからです。共感や感情的知性を欠いている人は、APIやGUIや言語の設計をすべきではないのです。

第5章ジョシュア・ブロック(Google チーフJaveアーキテクト)より。プログラミング能力の高い人材は、大抵その組織で最も知的な人たちではあるが、コーディング以外のこと組織や経営についての判断まで、彼らにまかせると(略)という話。

219ページ

私はどちらかというとCを手で書くよりはプログラムで自動生成しますね。そのほうが簡単ですから。

第6章ジョー・アームストロング(Erlang開発者)より。全てはメタになる。ベビメタDeath！

240ページ

自分のプログラムを800桁の数字として入力するわけです。

第7章サイモン・ペイトン・ジョーンズ(Haskell開発者)より。学生時代にプログラミングした、手作りみたいな計算機でのプログラミングについて。往時のマイコンのように、ブートコードのマシン語を直接入力していたのだろう。

313ページ

一時期私はあの本をモニタスタンドにしていました。私の持っている本で一番大きな本だったし、ちょうどいい高さになったのです。これは好都合でした。いつも目の前にあったので、リファレンスとしてよく使うようになりました。

第8章ピーター・ノーヴィック(Googleの研究部門責任者。元NASA計算科学部門)より。

351ページ

違っていると感じるのは、文章の主たる読者は、コンピュータとは随分違った種類のプロセッサを持っているということです。

第9章ガイ・スティール(Schemeとか色々)より。一般人は逸般人とは思考回路が違うよね。

。。。

佐藤文隆編「林忠四郎の全仕事」(京都大学学術出版会)より引用するページ。

多くの理論天文学関係者の大ボス(私からするとボスのボスのボス)である、 有名科学者の関わった文章や対談などをまとめた大著。 いわゆる追悼文集に相当するものであるが、第一章は自叙伝であり、 自らが晩年に関係者に配布していた手作り？の自伝を編集したもの。 残りの章は日本語論文や解説、対談、弟子達による文章、追悼文を集積して、 全781ページハードカバーという、故人の句作や詩作を集めたような追悼文集では決してない。 自叙伝部分は、生い立ちから始まって、旧制高校、大学入学、第二次大戦時の従軍を経て、 研究者となりという折々の出来事が事細かく記述されている。 研究者となって以後は、研究上の重要なアイデアや論文だけでなく、海外出張の顛末、 学生とのやりとりや、他大学の先生との会議なども含めて、とにかく細かい。 詳細な日記をつけていないとここまでの記録ができないのではと思うが、 博覧強記すぎるのかもしれない。以下、面白いと思ったところを抜粋。

48ページ「電子計算機の初めての使用(1959)」 1959-1960にNASAゴダード研究所に滞在した時に。

さて、ここNASAでのゼロックスによる便利な複写やIBMの高速電子計算機の使用は私の初めての経験であって、帰国後は、教室にゼロックスやカード・パンチ機を導入するとともに、共同の計算室を設置するように努力した。また、東大教養の小野周教授らと協力して、わが国の共同利用の計算機センターを七箇所に新設し、その一つを強大に設置することに数年間尽力した。

林忠四郎とその門下生たちの重要な仕事は、偏微分方程式の数値解に依存していて、日本で計算機を正しく使い始めた最初の研究者の一人であり、その影響は今でも非常に大きい。天体物理学研究では数値計算とコンピュータシミュレーションは欠かせない研究手法であり、そのため計算機の利用が得意な天文学者が多いのは偶然ではなく必然なのだろう。計算機科学や計算機工学が確立する前、研究者の必要性から計算科学は存在していた。

81ページ「パソコンとソフトの勉強」 1984年に京都大学を定年退官した後のこと。

退官の時に贈られたパソコンとソフトは富士通製であり、これは、富士通のハードウエアは信用できるとの松田君たちの推薦によるものであった。しかし、実際はソフトが駄目で、大原謙一君の約二ヶ月の援助の後に、やっとパソコンは動き出した。東大物理に勤務していた佐藤勝彦君が後で言うところによると、当時、MSDOSのソフトを用いていた同君のNEC製のパソコンは簡単に動いていたのである。つまり、UNIXを簡略化したMSDOSの有用性を軽視したのは、富士通の技術者の大きな失敗であった。

自分の弟子達なのだから君づけなのは当然ではあるが、著名な先生の名前が頻出する。

ところで、私は大型計算機の設置の運動をしていた頃から、計算機のハードウエアとソフトの両面の機構の奥底を一度は知りたい（ブラック・ボックスの中身を一度は覗きたい）と思っていたので、この機会に長期にわたる勉強を始めた。

改めて念の為、定年退官後のことです。

パソコンの買い換えは、1986年にはFM11に代えてFM16βを購入。1996年には富士通製のMSウィンドウ95を購入、1999年にはNEC製のウィンドウMS98を購入、2003年にはウィンドウMSXPを孫の浩平に頼んで組み立て、2006年にはデル製のウィンドウMSXP(3.2MHz)を購入した。

最後のクロック周波数はあきらかに誤植。この時期は定年退官後の70代後半から80代にかけての時期なのだから、やはりメモ魔だったのか。

93ページ「放送大学の視聴」

1999年には、36吋の大型テレビを購入し、パラボラ・アンテナを二つ設置して、BS(衛星放送)とCS(通信衛星)の両方の受信が可能になった。有益なテレビ番組がないときに、CSによって無料の放送大学の番組を見ることができるようになったのは、予想外に幸いであった。

その頃、弟子の一人であるこ杉本大一郎先生(ボスのボス)が、東大退官後に放送大学で教鞭をとっていた。

97ページ

さて、私のこれまでの経験によると、GoogleのWikipediaの記事は、かなり信用がおけるものである。

2000年代にはパソコンを研究の計算や作図に使うだけでなく、情報取得にも利用するようになった、その経験から。

99ページ 2004年の開催された「星形成と太陽系起源」シンポジウムでの研究発表で、自らがおこなったN体計算について報告。この時84歳ですね。

このN(=30〜1000)体の運動の計算は、倍精度と4倍精度の、4次と6次のルンゲ・クッタ法を用いて行った。重力のカットオフはしなかった。計算の結果、初期が非回転の場合は、自由落下時間の3倍以上の時間がたつと、粒子の存在領域は、球状のコア・楕円体状のハロー・拡がったエスケープ領域の3領域に分かれることを見いだした。(中略) さて、精度や次数の異なった計算で得られた値を比較することによって、このN体問題の本質はカオス的であるという重大な発見をした。すなわち、決定論的な予測ができる時間には限界があり、また初期条件への強い依存性があることがわかったのである。

個人的には、大学院の研究テーマとしてN体計算(SPH法)を選んだのは1995年ころのことであるが、理研に在職した2004年頃に演算器のハードウエア実装をやりはじめ、その関連で4倍精度演算について教示をうけたのは2006年のことであった。その後、4倍精度計算の高速化手法の研究をソフト・ハードとも色々とやってきたが、すでにその2年前「4倍精度でN体計算」をおこない、自明ながらもそのカオスの性質を計算機上で発見し報告されていたことには、複数の意味で驚いた。科学の世界は巨人の肩に乗って発展するものではあるが、孫悟空がお釈迦様の手の上から逃れることも難しい！

1章の自叙伝部分は他の部分も読み物として面白い。決して天体核、天体物理学という専門分野だけでなく、日本の物理学の歴史にも関わりが深い。うまく編集すれば新書として売れそう。

1章の最後には、付録として年表の他、論文業績リスト、指導学生のリスト、そして媒酌人をしたカップルのリストがある。

テストするコードをテストしてみようと思い立ち、Google Test https://code.google.com/p/googletest/ を利用した記録。

数値シミュレーションのコード、特にN体シミュレーションのコード開発は、 自動テストをしながら開発するのには本質的にあまりむいていない。 これは、全体テストにかかる時間のほうがシミュレーションの計算より時間がかかりかねないし、 MPIで並列化されているコードで大規模なテストはバッチジョブとでしか実行できない場合、 自動テスト化するにはそれ自体が別の一仕事になってしまう。 個人でおこなうには多少無理がある。

それでも、個々の関数や小さなクラスにユニットテストをおこなう価値はあるかもしれない。 必要があって試してみようかと思いつく例として、 既存のコア関数を最適化してSSE/AVX化するという場合がある。 AVX化して4倍高速にするぜといきっていても、ベクトル長に関わるループの間違いが混入したりするのはありがちだし、 intrinsicを使うことで計算自体を間違う可能性は大きい。 これは既存のコードをCUDA/OpenCL化するときも同様で、 参照実装とOpenCLコードの結果を比較せずにコードを書くことはできない。 このような場合、経験的にはその度に自分で小さなテスト用プログラムを書いて、 適当なセットの入力値に対して同じ結果になるか、 あるいは誤差がulpだとか期待される相対誤差範囲にあるかを調べて「できた！」と思って終わりにすることはよくある。 ほとんどの場合にはそれで問題なく、広い意味でのSIMD化は成功といえるのだが。しかし。

以下のURLに、International HPC Summer School 2014の講義スライドがアップロードされている。 その中に”HPC Software Engineering”と題した講義があった。

Internation HPC Summer School 2014 Schedule https://www.xsede.org/web/international-hpc-summer-school/2014-wiki/-/wikid/9khB/2014+Wiki/Full+Schedule

スライド：lindahl_hpc_software_engineering_20140605.pdf

ローカルコピー：http://galaxy.u-aizu.ac.jp/permanent/lindahl_hpc_software_engineering_20140605.pdf

ここで紹介されているGROMACSというオープンソースで開発されている分子動力学コードの開発では、 様々なsoftware engineeringの手法が利用されている。 具体的に利用しているソフトウエアは、Git(ソースコード管理), CMake(build system), Doxygen(C++クラス文書生成), Google Test(テスト), gerrit(コードレビュー), Jenkins(継続的インテグレーション), Redmine(バグ・プロジェクト管理), copernicus(高レベルでのコード実行の制御)になる。 最後のものはsoftware engineeringとは直接関係ない。詳しくはスライドを参照して欲しい。 これらの手法を一度に全て採用するわけにも行かないが、それぞれは以前より見聞きしたことはあったので、 GoogleTestとCMakeを使って、最近必要になったコード書き換えをおこなった。

このコード書き換え(リファクタリング)は、ある関数を動作を変えずに速度について最適化したい、という目的がある。 入力のパターンは膨大なため、全てをテストすることは不可能であるから、 これまで似たようなリファクタリングをやった時には、乱数をベースに入力を発生させて、 元の関数の出力と比較するということを100万パターン行うプログラムを作った。 そして、そのプログラム自体を、shell script等で複数回実行するということをやってきた。 よくある「車輪の再発明」で、これまでこういうことを違う問題に対して何度も行ったことがある。 そのたびに乱数で入力を発生させて処理を繰り返すコードを、単純なことでもあるので、 昔のソースを探すよりは手っ取り早くその場でアドホックに作っていた。

このような時にGoogleTestを使うとどうなるか。結果からいうと、このような無限に近いパターンを検証したいなら、 非常に簡単だし役に立つ。そのためにはGoogleTestの利用方法を学ぶ必要があるが、それは 「Google Test ドキュメント日本語訳」 http://opencv.jp/googletestdocs/index.html を参照した。 さらに、Google Test自体をより簡単に利用するために、githubで公開されている https://github.com/jpilet/gtestcmake CMakeスクリプトを利用した。ただし、このオリジナルは今回利用した環境では、pathの問題があったので、 修正を加えた。一行ではあるが修正したものは https://github.com/dadeba/gtestcmake にforkした。 このCMakeスクリプトを使うと、最初にconfigurationした時に Google Testの最新コードを自動的にcheckoutしてくれるので非常に便利であった。

1
2
3
4
5
6
7

cd project_directory_somewhere
cp -r /elsewhere/gtestcmake/ .
vim CMakeLists.txt
mkdir build
cd build
cmake ../
make

CMakeLists.txtにGoogle Testに対応したテスト用のコードの依存関係を書く。 cmakeコマンドを実行するとMakefileがbuildディレクトリに生成されるので、 あとはそれを使ってテスト用の実行ファイルが作ることができる。 上記ドキュメント日本語訳には、様々なテスト記述方法や　オプションが説明されている。 その中で特に今回の目的に役立ったのは「テストを繰り返す」 http://opencv.jp/googletestdocs/advancedguide.html#adv-repeating-the-tests の部分。 Google Testのライブラリとリンクすると、ここで説明されている実行時引数を指定することができる。 今回のように検証すべき入力が乱数で生成可能で、しかも無限のようにある場合、

1

./my_test_code --gtest_repeat=1000 --gtest_break_on_failure

とするだけで、このテストを1000回繰り返し、 途中1つでも失敗したらそこで実行が止まる。別途スクリプトを書いたりする必要がない！ これは予期していなかったよい点であった。 この方法を使って一応十分に検証したコードの速度は、リファクタリング前より3倍高速になった。めでたし。

Blaauw&Brooks “Computer Architecture: Concepts and Evolution”(1997)から、 コンピュータアーキテクチャの系統樹。

1997年の出版のため、だいたい1980年代までの話題についてのみ。 ペーパーバックで1000ページを超える大著であるが、ペーパーバック版は2分冊になっている。 著者らはIBM 360の設計者・マネージャーであっため、 前半は360を開発した時の経験に基づいてコンピュータアーキテクチャという単語の発明者(Brooks)が、 コンピュータアーキテクチャを講義するという形になっている。 アーキテクチャを説明する手段として、APL言語(!)を採用していて、それによりISAの詳細を紹介している。

後半は「A Computer Zoo」と題して、系統樹にある過去のエポックメイキングなアーキテクチャをすべて「定義」するという、大変なことをおこなっている。 この動物園の最初のエントリーは「Difference Engines of Babbage and Scheutz」である。 そのあと、Harvard Mark I、Zeus Z4と続く。 CDC 6600やCray 1など、いわゆるスパコンのISAが説明されているだけでなく、 それらの浮動小数点演算やベクトル演算についても記述がある。 この2つの計算機はその子孫がまがりなりにも残っているのに対して、 いきなり系譜が途切れてしまうIlliac IVについては記述がないし、Connection Machineはそもそも載っていない。。。 一方で、VAX 11やPDP8などのUNIXに深く関わる計算機や、 最後にはいくつかマイクロプロセッサについても。 とわいえ、この系統樹に載っていて現代でも広く使われているISAは、SPARC, MIPS, 8080A, AS400とその子孫のみだろうか。 当たり前だが、ARMやGPUのようなマルチコアについては一切触れられてはいない。

なお、IlliacについてはHord “The Illiac IV”(1982)という文献があり、 Connection MachineについてはHillis “The Connection Machine”(1985)がある。 前者をAmazonで購入してみるとプリントオンデマンドのものが送られてきた。 中身は二昔前の論文のようにタイプライタで組版された、非常に読みにくいものであるが、内容は面白いところもある。 後者は日本語訳(＋日本語版独自内容の第二部あり)が1990年に出版されていて、すごく時々大きな書店で在庫をみかけることがある。 見かけたら、ぜひカバーにある訳者紹介を見て欲しい。ニヤリと出来るから。

距離の計算は粒子が関わるシミュレーションでは必ず必要な演算で、それを高速化することが結構クリティカルである。伝統的なGRAPE式の演算の換算方法では、1相互作用(加速度とポテンシャル)の計算に必要な演算数を38演算としている。実際に必要な演算数を数えると、平方根を除いて20演算しかない。そのうち1つは除算であるが、平方根と除算は結果が出るまで複数サイクルかかるので、それらを10演算とすると39演算と換算できる。1演算の誤差は無視して欲しい。

この「換算」は最初のGRAPEが開発された90年代ころには、それなりに根拠のある換算であって、実際に必要なサイクル数を計測して決められている(Warren 1997)。ポイントはこの時に演算精度は指定していないことで、GRAPEが非常に効率が良く高速であった理由の1つは演算に単精度や倍精度のような標準形式を使わなかったからである。重力相互作用計算に最適に各演算に必要な演算精度を評価し、演算形式もradix 2の浮動小数点にこだわらずに、固定小数点や対数形式を組み合わせている。

2014年は最初のGRAPE-1が開発されてから25年周年とのことである。現代ではこの換算はあまりにも過大評価すぎる。いや、90年代の終わりくらいからすでにそうだったろう。それでもいまだに除算や平方根は加算や乗算のように1サイクルでは計算できないので、速く計算できるに超したことはない。逆数平方根が高速に計算できると、このような場合の除算を省略できるし、単位ベクトルの計算にも便利であるので、x86の命令には12ビット精度で逆数平方根をかえす命令が追加されている。より精度が必要であれば、それを初期値としてニュートン収束を適用すればよい。

そのような命令が追加される前、PCゲームの世界では有名なトリックがあった。

Lomont 2003より

この”magic number”の謎についてはWikipedia http://en.wikipedia.org/wiki/Fast_inverse_square_root の記事が詳しい。さらには逆数平方根についての学位論文(学部) Robertson 2012 http://shelfflag.com/rsqrt.pdf があるとは！ 起源を調べるとPCのゲームではなく、1986年のdraft論文までさかのぼることができるらしい。

http://www.netlib.org/fdlibm/e_sqrt.c

このソースコードのコメント部分に、シフトと整数加算で平方根/逆数平方根が約8bit精度で得られるとある(draftの著者はあのIEEE standard 754のKahanとK.C.Ng)。このコメントの”magic number”(逆数平方根)は0x5fe80000である。上の数字と少し違う。というのもこの”magic number”は必ずしも最適ではないからだ。

上記Wikipedia記事から参照されているLomont 2003とMcEniry 2007は、最適な”magic number”を見つける手法についてまとめている。どちらも浮動小数点フォーマットの説明からはじめて、なぜ2つの整数演算で逆数平方根が得られるのかを証明したうえで、誤差評価をおこない、最適な値を計算機で検索するという方法を使っている。後者のほうがより一般性があり詳しい。結果、得られた最適な”magic number”は0x5f375a86 or 0x5f375a85であるそうだ。Lomontの方法で前者が計算され(x86と関係あるんですかねぇ…)、McEniryはそれに加えて85でも同じくほぼ最適であるとした。Robertson 2012は、McEniryや他の仕事を参照して、0x5f375a86が最適値であることの証明をしている。

さらに、2014年のblog記事 http://multigrad.blogspot.jp/2014/04/math-evolution-and-dirty-tricks.html では、遺伝的プログラミングを使ってこの数字を発見しようと試みている。しかも、”magic number”だけではなく、組み合わせる整数演算もパラメータとして、遺伝的プログラミングで探索・最適化してみたという。結果は0x5f33f870。得られた数式は”(a+a) - (x >> 1)“だそうだ。惜しい。

得られる演算精度は約8bitだから、符号部指数部とつなぐと1+8+7(hidden bitを除いて) = 16bitになり、数値シミュレーションのようにそのあとにニュートン収束を使うことが前提ならば、一番最初の”0x5fe80000”で実は十分である。Quake 3などのゲーム用CGで、単位ベクトルの計算等に使うのであれば8bit精度の中で最適なほうがよいのだろう。

2009年の最近の話題 http://www.geocities.jp/andosprocinfo/wadai09/20090314.htm エントリによると、HPC-ACEでは8 bit精度の逆数平方根(と逆数)、PowerISA V2.06では5 bit精度の命令が追加されている。とすると、テーブルを使うまでもなく、この方法を使っているのかもしれない。一応、1e-12から1e12の間の全ての単精度変数について、誤差の頻度分布を計算してみた(単精度変数を網羅してテストする手法の参考： http://randomascii.wordpress.com/2014/01/27/theres-only-four-billion-floatsso-test-them-all/ )。ただし、引き算は上位16bitだけ。”Newton N”はニュートン収束をN回したもの。SSE2のrsqrtps命令は比較のために。

rsqrtps命令と比べると、誤差の頻度分布は変な形をしている。ニュートン収束を2回してもまだ完全に単精度にはならない。単精度でニュートン収束が3回以上必要だとすると、なんということでしょう、「GRAPE換算」よりサイクルがかかってしまうではありませんか。biasも心配ではあるし、安全を取るならば5bitでもテーブルを使った近似なのかなあ。

なおGRAPEシリーズでの平方根と除算の取り扱いは世代によって異なる。最初のGRAPE-1は、そもそも全ての演算がROM化されていたので、平方根もROMをひくだけだったろう。その後継である低精度バージョンのGRAPE-3やGRAPE-5は、入出力は整数でパイプライン途中は対数形式の浮動小数点を採用していて、この場合平方根は1bitシフトするだけという夢のような高速化が可能であった。一方いわゆる高精度なGRAPE-2にはじまる偶数シリーズは、世代ごとに詳細が異なるし全てを把握しているわけではないが、基本的には内部は通常の浮動小数点演算形式のため、テーブル補間をベースにしたものであろう。分子動力学用のMD-GRAPEシリーズは、距離を引数とした任意のテーブルを設定することができるので、明示的に平方根は計算していないと思う(要確認)。最後に、我々が設計した高精度演算でプログラマブルなGRAPE-MP( http://dx.doi.org/10.1016/j.procs.2011.04.093 )は、128 bit浮動小数点演算をハードとして実装しているが、加算器乗算器に加えて単精度相当の逆数平方根のユニットを持つ。これに採用したのは2次のチェビシフ補間である。この場合三回ニュートン収束すると必要な精度が得られる。

2012年ころから、GPUクラスターでの並列ツリーコードがいくつか発表されている。

私の関係者や知るところでは、Ogiya etal. (doi:10.1088/1742-6596/454/1/012014) と Bedorf etal. BONSAI (doi:10.1016/j.jcp.2011.12.024) がある。 二つのコードは、グローバルな並列化方針は似通っていて、Makino 2004ベースであるが、 GPUの利用方法が異なる。 Ogiyaらのコードはツリー構築はホストでおこないツリー走査と相互作用計算をGPUでおこなう (Nakasato 2011 の方法 http://dx.doi.org/10.1016/j.jocs.2011.01.006)。 それに対して、BONSAIはツリー構築、ツリー走査、相互作用計算を全てGPUでおこなう。 どちらもツリー構築のアルゴリズムはMortonキー＋ソートを使っている。

以下の図はそれぞれのコードの性能を筑波大学計算科学センターのHA-PACSで評価したものである。

Ogiya etal.のstrong scaling

N = 32Mの場合

BONSAIのstrong scaling

http://www.ccs.tsukuba.ac.jp/files/slides/sympo2013/24-fujii.pdf より

加速度の演算精度パラメータが不明であるのと、粒子分布が異なるため単純には比較できない。 どちらも、HA-PACSでGPUあたりの粒子数が多いならば、よくscalingしている。

最近、我々は(Watanabe & Nakasato HEART’14発表 プレプリント http://arxiv.org/abs/1406.6158 )、 また別の並列ツリーコードを実装し、HA-PACSで性能評価をした。 我々の手法は、PPPT法(Oshino, Funato & Makino 2011)にもとづいて、 短距離部分と長距離部分で力を分割して、それぞ別のツリーによって計算をおこなう。 元々のPPPT法は衝突系のために提案され、短距離部分は四次のHermite法で高精度に積分をし、 長距離部分はツリーを使って演算量を減らしまた時間積分も二次のLeapfrogにすることで、 トータルの計算時間を減らすことを主眼としている。 我々の方法は、計算時間を減らすと共にノード間通信量を減らすことで、 よりよいスケーリングを達成することを目的としている。 短距離部分のツリーは局所的なツリーとなるため、プロセッサ数が多い(かつ最適に領域分割されている)ならば、 近接ノード間での通信のみで更新できる。 長距離部分のツリーは全対全の通信が必要になるが、 PPPT法と同様に二つの成分毎に時間積分の頻度を変えているため、頻繁には発生しない。 これによりさらなる演算量と通信量の削減を達成できる。 GreeMなどのTreePM法と比べると、FFTで計算していたPM部分を、 独立したグローバルなLETツリーとして計算するものである。 どちらのツリー計算も、Nakasato 2011の方法でOpenCLによりGPUで加速化している。 領域分割はMakino 2004の方法であり、負荷分散はIshiyama etal. 2009の手法である。

我々のHybrid Tree法のHA-PACSでのstrong scalingの評価は以下の図になる。 ただし、これは「短距離のKickを3回、長距離＋短距離のKickを1回、領域分割を1回」を1サイクルとして、 このサイクルあたりの計算時間である。通常のツリー法は毎step長距離＋短距離のKickをしていることに相当する。 粒子はPlummer分布。短距離力のカーネルはDLLである。 他パラメーターがたくさんあるが、説明はとりあえず省略。

N = 32Mの時のbreakdownは以下の通り。 ただし上記の説明の通りT_total = 3*T_hard(短距離)+ T_soft(長距離＋短距離) + T_DD(領域分割)である

短距離部分の計算はグローバルなLETツリーの部分よりも高速であるしscalingも若干よい。 領域分割の時間がflatになってしまうのは、他のコードと同じ。 問題によっては長距離部分の計算をおこなう頻度をもっと減らしても問題ないであろう。 渡辺さんの修士論文では、それをパラメータとして、エネルギー保存を比較し、 様々な場合の比較をおこなっている。近いうちにTechnical Reportとして公開予定。

概要

重力多体問題のために実装されて発表された、並列ツリーコードについてまとめる。 できるだけそれまでになかった新しいアイデアが導入された論文を古い順にとりあげているので、 この分野の全てのコード・実装を取り上げているわけではない。 並列ツリーコードの特徴は、以下の凡例にある三つのポイントで分類できるので、 各項、それを最初に提示した。

凡例

Parallel Hierarchical N-body Methods (Salmon 1990, PhD thesis)

論文がどこにあるのか検索中。「Parallel Computing Works!」に書いてあったかもしれない。

Parallel Hashed Oct-Tree (Warren & Salmon 1993)

ツリーのデータ構造として、Barnes&Hut(1986)のFortranでの実装では、 八分木をポインタ(のようなもの)で表現していた。この素直な実装は、 ある節(以下cell/セル)の粒子が少ない場合には、効率が悪い。 そこで、fully threaded treeという方法では、八分木の代わりに2つのポインタでツリーを表現できる。 というのは、この論文とは関係ない背景知識で、 この論文ではhashを使って直接ツリー構造を保持するという手法を提案している。ハッシュに使うキーはMortonキーと呼ばれる、三次元の粒子座標を一次元の値にマップしたものである。Mortonキーは、3 bitごとに区切られていて、八分木のある階層のどの三次元象限にいるかに対応しており、このキーが計算できると、簡単にツリー上の上の階層や下の階層のキーを計算できる。粒子の分布を一次元であらわすことができるのなら、これほどよいことはないが、問題はこの一次元のキーが疎であること。そこでWarren&Salmonは、配列をhashで置き換えてこの問題を解決した。このhashはMortonキーを入力として、粒子のindexを返すようなものである。これにより、あるセルの子セルを全て数え上げる、などの典型的なツリーの処理は、子セルのキーを計算して(自分のキーに3bitをつけたす)、hashから子セルのindexを得ることで実現できる。さらに、並列化でもMortonキーの性質を利用すると、全プロセッサがひとつのglobalな木を構築すると仮定した場合、個々のプロセッサは自分の粒子のMortonキーを計算してローカル木を作るだけでよい。Mortonキーは全プロセッサで一元化されているので、あるセルのデータがどのプロセッサにあるのかは、簡単に求めることができる。これを利用して、彼らはLocally Essential Tree(LET)という方法を提案した。これは、あるローカル木の粒子と相互作用する他プロセッサの粒子とセルだけを含むような木である。LETが構築できれば、その後のツリーに関わる演算は通信なしでできる。Mortonキーを使えば、LETを得ることは原理的には難しくない(実際はきっと面倒)。Mortonキーのよい性質は、2つの粒子のキーの値が近い場合には、粒子同士が空間的に近いところにマッピングされることである(このようなキーに基づいて粒子の座標をつないで得られる曲線は、その意味で空間充填曲線とよばれる。空間充填曲線の考え方は、メモリアクセスの局所性を高めるのに有用なので、HPCではよく使われるアイデアのひとつ)。ただし、Mortonキーでソートした場合、一部空間でのジャンプが含まれるため、より空間的な局所性の高いマッピングがPeano-Hilbert(PH)キーとよばれる。これも空間充填曲線の一つである。Warren&Salmonは、粒子の領域分割にはPHキーを使っている。

Implications of Tree methods for multiprocessor architectures (Singh, Hennessy & Guputa 1995)

http://dl.acm.org/citation.cfm?doid=201045.201050

この論文だけは、他のエントリと異なり、 コンピュータアーキテクチャの研究者(へね)がツリー法をマルチプロセッサで性能評価した論文である。 第一著者であるSinghの博士論文 “Parallel hierarchical N-body methods and their implications for multiprocessors” (Doctoral Dissertation, Stanford Univesity, 1993)を投稿したものと思われる。

すごく長い。以下、Section 9のメモを。

ここでは、SalmonのLETを使ったツリー法について色々と議論している。 論調は、LETとメッセージパッシングは使えねー、ってこと。 ただし、これは90年代前半の仕事なので、今とは対象としている問題のサイズが圧倒的に小さい。 この部分である意味かなりdisられている、Salmonの結果でも512プロセッサで10万粒子とからしい。

ちょっとおもしろかったので引用。

The difference is also demonstrated by our experience with having groups of graduate students implement the application on SAS-CC and message-passing machines, in a ten-week parallel programming project course at Stanford University. The SAS versions were produced very quickly, and yield very good speedups on the DASH multiprocessor. A message-passing version, however, is yet to be completed within the time allotted for the project

この部分の前に、Salmonの提案によるMPを使ったLETのツリー法は、 SASでの実装より数倍時間がかかったと書いてある(誰が？)。 その上で、ツリー法の実装をStanfordの大学院生に10週のプロジェクトとしてやらせてみたら、 共有メモリ(SAS-CC)とMessage-paasing(MP)の生産性の「違い」が示されたという。 SAS-CCとはshared address space-cache coherentの略。著者らはSAS-CCな計算を推しているのは言うまでもない。

まず、LETではメモリ使用量が非常に多くなる。 これはLETでは保守的にTreeを構築するので、 結果的には無駄な粒子とノードをLETに追加してしまうため。 よって、メモリ利用量も多いし、そのためのデータの交換に必要なデータ量も多いし、 それを実現するためのプログラミングも非常に複雑になる。ボロクソである。

一方SAS-CCでのツリーの実装では、データ交換は暗黙的のためプログラミングが簡単であり、 また、本当に必要最低限のデータ交換しか発生しない。 それは、SASでの実装では、LETを構築するのではなく、一つの共有されたtreeを実装するためで、 “Costzone”と呼ばれるロードバランスの手法と対になっている。 利用するメモリ量は、シリアルのコードとほとんど変わらない程度になる(はずである)。 ただし、彼らが議論している粒子数はたかだか16000粒子程度であり(Table 2参照)、 それのための必要なキャッシュ量は16KBとあり、現代には通用しない。 そもそも、現在でも128プロセッサを超えるような共有メモリ機を、 ハード的に実現するのは現実的には非常に難しい(性能が出ない)。

それでも、無理やり彼らの結論を敷衍すると、ハード的に共有メモリを実現することは難しくても、 N-bodyに専用の共有メモリフレームワークを改めて考えて、 まともに実装することには意味があるかもしれない、ということになる。 実際、彼らは仮想ポインタを利用したツリーの実装と、 ハッシュ値を使ったツリーの実装があり得ると提案していて、 これらはある種N-Bodyのための共有メモリフレームワークになっている。 脚注には、ハッシュ値を使ったツリーの実装を「この論文」を書いている時に、 Salmonらがおこなっているとある。これは上のHOT論文のことである。

Parallel Tree Code (Dubinski 1996)

J.SalmonのPh.D Thesis(1990)で提案されている、ORBによる領域分割を利用した実装。ORBでは、粒子を含んだ領域を再帰的に二分割する。単なる八分木のように等体積になるように分割するのではなく、ニ分割する際には両側で必要な計算量が等しくなるように分割する。分割する面の向きは、原理的には任意の方向にできるが、このコードではxyz軸に垂直な面でcyclicに分割をする。ORBで分割された各プロセッサの粒子は、個々にローカルな木を構築しする。そのあとLETを作成するという部分は、Warren&Salmonと同じである。”with Linked List”と書いた意味は、上記”fully threaded tree”のこと。

Forest Method (Yahagi, Mori & Yoshii 1999)

この論文のポイントは、”Forest Method”と名付けられた、新たな領域分割法を提案していることである。これはボロノイ分割による領域分割と、そのボロノイ領域を活用したロードバランシングがみそである。ボロノイ領域により粒子を分割すると、個々のプロセッサの担当する粒子分布は球状に近くなる。一方で、ORBのような分割方法では、粒子分布はどうしても直方体になってしまう。三次元の立体で表面積と体積の比が一番小さいのは球であり、領域分割された並列計算においては、この比によって通信量が決定されるため、可能であれば領域分割は球に近いほうが有利である。また、ボロノイ分割では、他の領域との境界は、それぞれの領域の中心点を結んだ直線を垂直に分割する面である。本来のボロノイ分割では、分割面はこの直線の中点を通るようなものである。Forest Methodでは、この分割面をロードバランスに応じて、中点からずらすことで、各プロセッサの領域の体積を調整し、それによってロードバランスを実現する。この時に「計算負荷(ロード)」の平均化に、拡散方程式を利用していることも特徴である。領域分割が面でおこなわれているので、LETを構築することは、比較的簡単になる。

Gadget-1 (Springel, Yoshida & White 2001)

天体物理学シミュレーション業界では最も広く使われている実装の初期のもの。ポイントは、ローカルツリーを作った後で、LETを作ることをしないで、力を計算する粒子(GRAPEの用語でいうとi粒子)の情報を通信する。この時に各プロセッサが、リング上に接続されていることを仮定して、バケツリレー式にi粒子のデータを隣に渡す。受け取ったi粒子に対して、ローカルツリーからの加速度(重力とSPH)を計算し、次はi粒子のデータとともにこの加速度も通信する。このような手法を採用している理由は、Gadgetではではblocked individual timestep schemeを使っていて、つまりある時間に積分する粒子は何らかの条件で選んだ粒子のみである。この選ばれた粒子の数は、通常全粒子数よりかなり少ないので、LETを作るコストが馬鹿にならない。欠点は、この方法ではプロセッサが多くなると性能が落ちることである。大きな利点は、プログラムが非常に単純になること。Gadgetでは、ツリー構築の時間を最小限にするため、粒子が動いた場合にも、全体の木構造を作り直すのではなく、局所的な変更のみを反映させるなど、詳細な最適化をおこなっている。これは実際には結構大変なので(セルの重心や四重極モーメントの予測子を計算する！)、LETを作ったほうが実は簡単かもしれない。宇宙論的シミュレーションで必要な周期境界条件は、Ewald和によって得られた加速度場を利用している。

Gasoline (Wadsley, Stadel, Quinn 2004)

GasolineはPKDGRAV(Stadel 2002 Phd http://adsabs.harvard.edu/abs/2001PhDT……..21S )を拡張し、元々は重力のみだったものに、SPHを加えたコードである。このコードが、他の実装と最も違うところは、基本的なデータ構造として二分木を採用していることである。ツリーが出来てしまえば、二分木か八分木に本質的な違いはないが、ツリー構築の方法がbottom upかtop downかという違いがある。一方で、今では一番下の階層にある節の粒子数が1個になるような純粋な八分木を構築することは無駄であるとわかっているので、Gasolineでも一番下の階層の粒子数は8個以下とかになっている。その意味では、Gasolineで採用しているものは、PKDGRAVの論文のタイトル通り”k-d”ツリーである(ほぼ二分木。他のコードはほぼ八分木)。領域分割はORBである。並列化の方法は、ORBツリーにたいして直接ツリーウォークを行う方法である。つまり、各プロセッサはORBツリー(その下の階層は各プロセッサのローカルツリー)のツリーの根から直接ツリーをたどっていき、ローカルではないセルにたどり着いたら、暗黙的に（？）通信をおこなうという方法である。周期境界条件の扱いはGadget-1と同じ。Gasolineのその他の部分には目新しい部分はないので省略。

PKDGRAVについてはwebページがあったので、 詳しくはこちら：http://hpcc.astro.washington.edu/faculty/trq/brandon/pkdgrav.html

Fast Parallel Treecode with GRAPE (Makino 2004)

ORB分割では、並列プロセッサの数が2のべき乗であることを仮定していて、ある軸の方向には必ず二分割しかしないため、任意のプロセッサ数では利用できない。この実装では、その制限を排除し、分割数を任意として、より幅広い並列プロセッサでORB likeな手法を使う。もう一つORBと異なるのは、LETの構築の仕方である。ORBでは、各プロセッサの持つローカルな粒子群が、ORBによって作られたツリーとなっているので、他のプロセッサからは必要最低限に枝切りしたツリー(LET)をもらえばよい。結果、グローバルなツリーは、ORBツリーとローカルな八分木ツリーの組み合わせとなる。この実装では、ローカルなツリーを作るとのは同じで、そのプロセッサに属する粒子と相互作用する粒子とセルのリストを作り、それをローカルなツリーに追加するというアイデアを採用している。これにより、個々のプロセッサの持つツリーはほぼ最適となる(ORB+ローカルツリーでは、幾何学的に近い粒子が別のノードに割り当てられてしまう)。

Gadget-2 (Springel 2005)

Gadget-2は、Gadget-1の後継であるが、実装の詳細は大きく変わっている。この業界のスタンダードなコード。TreePM法とは、重力のような長距離力を、短距離成分と長距離成分に分離して、それぞれを別の重力ソルバーで解く手法であり、名前が表すとおり、短距離にはツリーを、長距離にはPM法特にFFTベースでPoisson方程式を解く。この方法は本質的にはParticle-Particle Particle-Mesh(PPPM)の拡張であり、短距離成分の計算をツリーで加速するということがポイントになる。なおTreePM法自体はSpringelの提案ではなく、Xu 1995によって提案された。Gadget-2ではより精密化された定式化によるSPHを採用しているだけでなく、SPHと組み合わせた多相星間物質のモデルや、星形成モデル、冷却関数等々、様々な物理過程も含んでいる。N体シミュレーションで最も大規模な計算のひとつはGadget-2(の特殊な最適化バージョン)でおこなわれた。TreePM法なので、周期境界条件はPMの部分をFFTで解くことで、自然と満たされている。

PPPM and TreePM Methods on GRAPE Systems (Yoshikawa & Fukushige 2005)

GRAPE用に最適化されたTreePMコード。TreePMでは、短距離成分の計算はカットオフをもったようなカーネルによる和に帰着するので、ほとんどSPHの計算と変わらない。そのため、そもそもグローバルなツリーを作る必要はない。この実装では、領域分割としては1次元のslab分割を採用している。slab分割は一般的には、あまりスケーラブルではないが、この実装でslab分割を採用している理由は、個々のノードがGRAPE-6Aをもつクラスターをターゲットとしているからである。GRAPE-6Aをもつクラスターでは、ノードあたりの重力演算能力が普通のx86_64コア1つより10倍以上高速である。これは何を意味するかというと、同じ計算をするのにノード数を十分の一以下にできるということである。MPPで1024コア必要な計算を、GRAPEのクラスターであれば32ノードとかで実現できる。その意味では、GRAPEやGPUのような加速演算装置を持ったクラスターは、必要なMPIプロセスの数を減らすことができるので、並列ツリー法には適したアーキテクチャといえる。この論文は、ノード数が数十であれば、PMのメッシュ数が512以上もあれば、slab分割であっても十分活用できることを示している。

GreeM (Ishiyama, Fukushige & Makino 2009)

Yoshikawa&Fukushige(2005, YF)のTreePMコードに、Makino(2004, M)提案のツリー法を組み合わせたもの。M2004, YF2005ともに、GRAPEと組み合わせたツリーコードであるが、GreeMはGRAPEボードをターゲットとはしていない。国立天文台のCRAY XT-4をターゲットとして開発された(のだろう)。そのため、ノードあたりの演算性能が落ちるため、YF2005とは異なり、領域分割はORBベースのものに変更されている。ただし、GRAPEの代わりに、Nitadori(2006)による、SSE2を利用した高速なGRAPE互換x86_64プロセッサ用ライブラリPhantomを利用している。Phantomを使うことで、1コアあたりの重力演算性能として10 - 30 Gflopsを得ることができる。CRAY XT-4では、NP = 1024コアまでよくスケーリングしていて、これはGadget-2よりもよい。

京コンピュータ用に最適化されたバージョンで2012年のGordon Bell Prize(Ishiyama, Nitadori & Makino)を受賞した。

HACC (Habib etal. 2012/2013)

2012年のGB prizeにて、GreeM(京)と直接対決(BlueGeneQ)して、敗退したコード。 2013年もGB filnalistであったが、残念ながら再び受賞を逃している。

2HOT : Improved Parallel Hashed Oct-Tree (Warren 2013)

http://dl.acm.org/citation.cfm?id=2503220

HOTが帰ってきた。熱すぎるということで。

最近流行のTreePMではなく、あくまでもツリー法にこだわり高速化している。 そのために、遠距離の成分は高次の多重極展開を利用することで効率化をはかっている。 全体的な並列化の手法は元のHOT論文からほぼ変わってはいない。 宇宙論計算に特化することで、一様密度の成分を差し引いたり、 宇宙論にあわせた積分手法を採用するなど、問題に特化した最適化が行われている。