本論文では、衝突系N体計算のためのいくつかの高速化手法について述べる。

衝突系の計算では、粒子の散乱や二体緩和の効果が系の進化に重要な影響を与える。衝突系の例としては、惑星系形成、球状星団、銀河中心などがあげられる。衝突系の計算では幅広いタイムスケールに対応すること、近接遭遇を正しく取り扱うこと、粒子の軌道を正確に積分することが要求され、このため、近似を用いない直接計算法と粒子ごとに独立の時間刻みを与える「独立時間刻み法(Aarseth1963)」の組み合わせが主に用いられる。本論文ではこの独立時間刻み法の高速化、高精度化について述べる。

まずはじめに、6次精度及び8次精度のHermite積分法について述べる。これまで衝突系の計算でよく用いられてきたのは、加速度の1階微分(jerk)までを解析的に計算して用いる4次のHermite積分法(Makino and Aarseth, 1992)である。我々は、2階微分(snap)もしくは3階微分(crackle)程度までの計算であれば現実的な計算コストで実現できることに着目し、6次/8次精度の積分法を構築した。加速度のみの場合から3階微分までを計算する場合の相互作用あたりの浮動小数点演算数を列挙してみると、38、60、97、144となる。6次の積分法は、典型的な要求精度(エネルギーの相対誤差が10(-8)程度)であれば、4次のものの3倍程度の刻み幅をとることができる。演算数と照らし合わせて、2倍程度の高速化が実現できたことになる。エネルギーの相対誤差で10(-10)以上の高い精度が要求される場合、8次の積分法が有効である。

次に、独立時間刻み法の新しい並列化手法について述べる。独立時間刻み法の並列化は、これまで十分に研究し尽くされてきたとはいい難い。効率のよい並列化の手法が確立されてこなかったため、近年急速に普及が進む超並列計算機の衝突系計算への応用も進んでいなかった。独立時間刻み法の並列化が難しい理由は、1、そもそも一度に積分される粒子が全体のごく一部であるため、計算の粒度が小さく低レイテンシなネットワークが要求されること、2、ある粒子への力を計算するためには他の全粒子の情報が必要となるため、アルゴリズム上の工夫がない限りネットワークのバンド幅への要求も高いものとなること、である。我々は新しい2次元並列アルゴリズムを開発し、少なくとも後者のバンド幅の問題は解決した。

これまでに試みられてきた独立時間刻み法の並列実装は、いずれも1次元並列のものである。N体計算そのものには、力を受ける粒子についての並列度であるi並列と、力そのものについての並列度であるj並列とがある。1次元並列のアルゴリズムでは、このいずれかを用いる。1次元並列アルゴリズムでは、プロセッサ数をpとしたときプロセッサあたりの計算量O(N2/p)に対して通信量はO(N)となり、結局のところネットワークのバンド幅が与えられた問題規模に対して用いることのできるプロセッサ数を制限してしまう。2次元並列アルゴリズムでは両方の並列度を同時に用いることで、通信量をO(Nl√P)にまで軽減できる。

2次元並列アルゴリズムそのものは、Makino2002で提案されている。しかしここで提案された手法には、1、平方数のプロセッサ数しか使えない、2、独立時間刻み法との組み合わせで負荷分散が不均等になる、といった問題点があった。我々はこの問題を同時に解決した、新しい2次元並列アルゴリズムを開発した。このことにより、問題規模が比較的小さいところから、良好なスケーラビリティを得た。1次元並列アルゴリズムでも、問題規模を際限なく大きくすれば一見良好なスケーラビリティは得られる。しかしながら衝突系の計算では長時間積分を要求するものが多く、小さな問題規模で高い性能を得ることが重要になってくる。

800CPUコアのCrayXD1を用いたシミュレーションでは、64kの粒子数に対し2.03Tflops、ピーク性能の57.7%を実現している。

最後に、GPU (Graphics Processing Unit)を用いた衝突系N体計算について述べる。GPUは低価格ながら高い浮動小数点演算性能を持ち、近年科学技術計算への注目を集めている。例えば、今回使用したGPU、GeForce 9800 GTX+は2万2千円程度から入手できるが約500Gflopsの単精度浮動小数点演算性能を持つ。GPUを用いた数値計算はGPGPU(General Purpose GPU)とも呼ばれ、GPUの数値計算への応用が広まった背景としては近年GPUがグラフィックス専用のハードウェアから徐々に汎用性を高めてきたことがあげられる。

最近になって、GPUを用いたN体計算でも数百Gflopsという高い実行性能を持つ研究結果が複数発表されている(Hamada and Iitaka, 2007; Belleman et al., 2008; Nyaland et al., 2007, Schive et al., 2008)。しかしながら、これらの結果は全て単精度で計算されたもので、そのまま衝突系の計算へと応用できるものではない。高い精度での軌道の積分を要求する衝突系の計算では、少なくとも「座標の差分」と「力の積算」を倍精度で計算する必要がある。前者が倍精度を要求するのは、近接遭遇時の桁落ちを防ぐためで、後者が倍精度を要求するのはたくさんの小さな力を足し合わせて全体の力とするためである。

我々はこの二カ所を、2語の単精度浮動小数点語と複数の単精度浮動小数点演算の組み合わせで置き換える「疑似倍精度法」を開発し、GPUで衝突系の計算に耐えうる精度をもつ実装を実現した。実装には4次精度のHermite積分法のものと6次精度のものとがあるが、両者ともGflops値でGRAPE-6A(Fukushige et al. 2005)の2倍程度の性能を実現している。また、6次精度の実装ではN<64kにおいて1TflopsのGRAPE-6(Makino et al. 2003)よりも高速な計算を実現している。SSEとOpenMPを用いてクアッドコアのCPU用に最適化した実装に比べても、8~9倍高速である。

我々の研究成果をもういちどまとめてみる。

1.従来用いられてきたの4次精度のHermite積分法に対して、演算数あたりの精度/速度に優れる6次/8次のHermite積分法を開発した。

2.新しい並列化アルゴリズムを開発し、超並列計算機が衝突系のN体計算に有用であることを世界で初めて示した。

3.疑似倍精度法を用い、単精度のみをサポートしたGPUを衝突系のN体計算に応用することに成功した。低価格なGPUを用いて高速かつ高精度の計算を実現した。

高密度恒星系は恒星どうしの近接衝突(大角散乱) が頻繁におこる系として特徴づけられる。このような衝突系の進化の数値シミュレーションは、粒子(恒星)の軌道を正確に積分する必要があること、そのために粒子ごとに時間積分のきざみを決める必要があるなど、無衝突恒星系の進化シミュレーションの数値手法とは大きな違いがある。論文提出者は、高密度恒星系の数値シミュレーションを大幅に加速するための3 つの新しい手法を提案した。

本論文は4 章からなる。論文は、共著者との連名ですでに出版あるいは出版予定であるが、そのすべてが論文提出者の似鳥啓吾が筆頭著者であるだけでなく、彼の主導で研究が進められたものであることを論文審査において確認した。なお、その論文の内容を主論文のなかに含めることについては、共著者の承諾書が得られている。

本論文第1 章は序論であり、研究の背景や従来の研究の問題点をまとめ、本研究の目的と意義、およびその概観を述べている。

第2 章では、高密度恒星系のN 体シミュレーションを効率的に行うために、論文提出者が提案した6 次精度エルミート積分法および8 次精度エルミート積分法について述べられている。この積分法の特徴は、今までこの分野のスタンダードとして使われて来た4 次精度エルミート積分法に対し、1 ステップあたりの計算量は、6 次精度エルミート積分法では約1.5 倍、8 次精度エルミート積分法では約2 倍に増えるものの、1 ステップあたりの時間きざみ幅をそれぞれ、1.5 倍以上、2 倍以上とることができれば、実計算時間は大幅に短縮できることである。論文提出者は、実際の系にこれらの手法を応用し、実計算時間において、2 倍程度の短縮をはかれることを実証した。また、全エネルギーを10(-8) の精度で保存させなければならないような系においては6 次精度エルミート法が、また、10(-11) の精度で保存させなければならないような系においては8 次精度エルミート法が非常に有利であることも示した。この研究成果はすでに発表されており、この手法を使った計算も実際に行われていて、実計算時間の大幅な短縮に貢献している。すでにこの手法が高密度恒星系の進化計算に大きな寄与をしていることからも、その成果は高く評価できる。

第3 章では、高密度恒星系の進化シミュレーションに必要な粒子ごとの独立時間きざみ法の並列計算の新しい手法について述べられている。重力計算においては、重力を受ける粒子をi 粒子とし、重力を及ぼす粒子をj 粒子とする場合に、重力計算の並列化を、i 粒子について行うか、j 粒子について行うかを、それぞれ、i 並列、j 並列として区別する。これまで行われて来た並列化は、どちらか一方に関してのみの、1 次元並列化がほとんどであった。1 次元並列化の場合には、プロセッサの数をp とすると、プロセッサあたりの重力の計算量がO(N2=p) となるのに対し、プロセッサ間の通信量がO(N) となるので、結局のところ与えられた問題の粒子数が有効なプロセッサ数を決めてしまう。これに対して、2 次元アルゴリズムは、i とj の両方の並列度を同時に用いることで、プロセッサ間の通信量をO(N=pp) まで軽減できる。論文提出者は、i とj の両方に任意のプロセッサ数を割り当てることと、各粒子は独立時間きざみでありながら、十分均等に負荷分散が達成されるアルゴリズムを開発し、その有効性を実際の大規模並列計算機上の数値シミュレーションにおいて実証した。この成果が今後大規模並列計算機上で大いに使われて、この分野の発展に寄与すると期待される。

第4 章では、最近大いにその性能をあげてきたSIMD ハードウェアの一種であるグラフィク描画専用プロセッサ(GPU) を高密度恒星系の進化計算に使うために、論文提出者が新しく開発した専用計算ライブラリについて述べられている。これらの汎用GPU の弱点は、単精度浮動小数点計算については、非常に高い性能を持つものの、高精度の計算が必要とされる高密度恒星系の進化計算に必要な、倍精度浮動小数点計算の能力は10 分の1 以下と低くなってしまうことである。論文提出者は、2ワードの単精度浮動小数点数を使うことにより、倍精度と同じような計算精度を達成できる手法を開発し、その実用性を実際の数値計算において実証した。特に、市販の安価なグラフィックボードを用いて、専用計算機GRAPE-6A を凌駕するような計算性能を達成したことは、高く評価できる。

最後の章は、6 次および8 次精度エルミート積分法に関する詳細を述べたAppendixである。

論文提出者は高密度恒星系の進化シミュレーション用に、新しく3つの手法を開発し、実用的な計算においてその有効性を実証した。また、その一部はすでに他でも使われ、この分野の発展に貢献している。これらの結果は高密度恒星系の進化シミュレーションを大きく進展させる画期的なものである。

以上を要するに、本論文は恒星系天文学の分野において、新しい知見をもたらすとともに、新しい発展の可能性を開くものである。よって本論文は博士(理学)の学位論文としてふさわしいものであると、審査委員会は認める。