next up previous contents
Next: 9 まとめ Up: 重力多体系の数値計算 Previous: 7 計算法 --- 空間領域

8 別のアプローチ --- 計算機を作る

ここでは、大規模な計算を可能にするためのこれまで述べてきたアルゴリズム の改善とは相補的なアプローチ、つまりハードウェアの計算速度そのものを速 くするという方向について簡単に述べたい。原稿依頼ではあんまりこの辺の話 は、、、とか書いてあったような気もするが、多体シミュレーションが数値実 験であるとすれば、ソフトウェアと同様にハードウェアも実験道具であり、ソ フトは自分で作るのにハードは買ってくるというのもある意味奇妙な話ではな いだろうか。

まあ、といっても、最近はソフト自体も大規模化し、個人、あるいは小さな研 究グループの手に負えなくなってきている。例えば、 FMM のような面倒なア ルゴリズムを、分散メモリ型並列計算機のようなプログラミングの手間が非常 に大きい計算機の上に実現し、しかもある程度まともな実行効率を得るのは容 易なことではない。

そういうことを考えると、ソフトだけでも大変なのに、ハードまで自分で作っ て、さらにそのハード用のソフトウェアまで自分で作るなんていうのは全く論 外である、、、と思われる方が多いのではないだろうか。

というわけで、以下必ずしもそうでもないのではないかという話をしたい。

8.1 基本的発想

多体問題シミュレーションにおいて、計算時間の大半というよりも 99% 以上 を占めるのは、粒子間の相互作用の計算、つまり運動方程式の右辺の評価であ る。従って、計算速度をあげるためには、何らかの方法で右辺を速く評価でき るようにするか、あるいはいろんな方法で右辺を評価する回数自体を減らすこ とになる。ここまででは、計算機は与えられているものとしてこれらをどうやっ て実現するかという話をしてきたわけである。

しかし、考えてみると粒子間の相互作用の計算というのは極めて単純な式の繰 り返し計算である。高級言語のプログラムで書けば数行だし、アセンブラでも 1ページくらいにしかならない。汎用計算機を使って多体シミュレーションを しているときは、ほとんどの時間汎用計算機をこの単純なプログラムを走らせ るためだけに使っているわけである。

このような状況では、粒子間の相互作用の計算だけを高速にやる補助回路を計 算機につけてやればシミュレーション全体を高速化できる。これは、例えばゲー ム用のコンピュータなどで、画像処理等のためにいくつかの専用回路をつけて そういった処理は CPU を使わないで済ませるというのと考え方としては同じ である。

しかし、多体問題用の専用回路には、実は例えば画像処理用回路に比べて原理 的に有利な点がある。画像処理の場合、3次元の特殊な処理でなければ処理の ための計算量はデータ量に比例する。従って、処理速度をあげるためには、デー タ転送速度も同時にあげる必要がある。このために、CPUやメインメモリと画 像処理/表示部分とを非常に高速なインターフェースでつなぐ必要が出てくる。 高速にしないといけないということは、現在の計算機の構成では CPU の「近 く」に持ってくる必要があるということであり、そのためには計算機の構成自 体がそういうものを考慮して作られている必要がある。例えば AGP バスがそ の例ということができる。

これに対し、粒子間の相互作用計算では、単純な計算法では計算量が粒子数 N の2乗に比例する。つまり、データ量の2乗に比例して計算量が増えること になる。もちろん、ツリー法などではそううまくはいかないが、それでもデー タ量に対して計算量が非常に多いということには変わりがない。このために、 それほど高速ではない汎用の I/O バスを通して接続してもまあまあの性能が でる。例えば、最近の大抵の計算機で採用されている PCI バスを使った接 続でもまあまあの性能を出せる。

そういうわけで、概念的には、粒子系シミュレーション用の専用計算機は図13 のような構成をとることになる。

  
図 13: 粒子系シミュレーション用の専用計算機の基本構造

専用回路の側は相互作用を計算するだけである。我々はこちらの部分に GRAPE (GRAvity PipE) という名前をつけている。とりあえず、もっとも単純な 直接計算で時間刻み一定の場合を考えると、専用計算機としては例えば系の全 粒子を受けとって自分のメモリにしまい、それから2重ループで全粒子への他 の全粒子からの力を計算する。最後に求まった力を送り返す。

とはいえ、これまで長々と話をしてきたように、直接計算で時間刻み一定しか 出来ないのではいかに専用回路が高速であったところで出来ることは制限され る。例えば、現在の汎用並列計算機で最高速のものは10Tflops 程度の性能を 持つが、これでも粒子数が になれば 1 ステップ 5 分、 なら 10時間となってすでに現実的とはいい難い。このために、専用計算機といって も実際に意味がある計算をするためにはこれまでに述べてきたような数値計算 法、特に独立時間刻みやツリー法との組み合わせが重要になる。これらについ ては次節以降で扱うとして、ここではとりあえず専用回路にすることでどの程 度高速になりえるかということを考える。

さて、専用回路にすることで高速化出来るのはそもそもなぜであろうか?そん なのは当たり前であるとなんとなく思っているだけでは、実際にどういうふう に回路を作ればどれだけ高速化出来るか、またそもそも高速化の理論的な限界 はどこかといった問題に答えることが出来ない。それでは、何か作ったとして もそれに意味があるものかどうか評価出来ないのである。

というわけでまた余談になるが、並列計算の場合でもそうだが、多くの研究と いうか並列プログラムの開発で、そこで実装することにしたアルゴリズムがネッ トワークの速度やプロセッサ台数が与えられた時に、もっともうまくいってど の程度の性能がでるかということは普通は理論的に評価できる。理論的にどう 振舞うのかわからないほど複雑なアルゴリズムではうまくいく気遣いはないか らである。しかしながら、著者の知る限りではそのような分析をやってから仕 事を始めた、あるいは途中でもちゃんとそういう検討をして方針を考えたよう に見えるプロジェクトはあまり多くないというかほとんどない。このような状 況ではアルゴリズムの進歩はいいところメトロポリス・モンテカルロ法的なも のでしかないわけである。

もちろん、計算機科学の専門家が書いた論文の場合には、少なくとも実装の性 能についての評価はそれなりにちゃんとしている。しかし、その場合でもアル ゴリズムの性能の理論的な評価はほとんどされていないことが多い。なぜ並列 計算法の研究が、そのような理論不在の状況にあるのかということは著者には 大きな謎である。

と、話を戻す。高度なアルゴリズムと組み合わせた時の話は後でするので、こ こではまずは単純な直接計算で時間刻み一定の場合を考える。

普通の計算機、具体的には普通の PCとかそのクラスタの場合には、重力計算 は汎用 CPU で実行される。この時に実現可能な速度は、どんなに良くてもCPU の理論ピーク性能を超えることはあり得ない。例えば、 750 MHz の EV6 Alpha CPU では、あらゆることがうまくいったとしてクロックサイクル毎に実行でき るのは加算と乗算が一つづつである。最近の x86 系の CPU である AMD Athlon や Intel P4 でも事情は同じようなものである。さて、普通に2粒子 間の相互作用を計算するのに必要な演算の数を数えると、乗算と加算が合わせ て19(だったかな)、平方根が1、割算が1になる。平方根や割算にどれくらいの 時間がかかるかは CPU によるしまた演算ライブラリの出来によっても変わる (細かいことをいい出せば、要求する精度にもよる。例えばニュートン法反復 で平方根や割算を実現する場合には単精度なら倍精度に比べて一回反復数がす くなくてもいいはずである。)が、まあここではとりあえずそれぞれ 10 演算 程度としておこう。すると、あらゆることがうまくいったとして1ペアの重力 を計算するのにはこれらの CPU では30 クロックサイクル程度かかることにな る。

まあ、実際にはいろいろ無駄があるのでこの2-3倍はかかるが、とりあえずこ れよりも速くはなり得ないという上限はこのあたりになる。 1 GHz の CPU を 使った場合で 1 秒に3 千万ペアを評価できることになる。

さて、ここで、「上限」として、「CPUを作るのに使っているすべてのトラン ジスタを演算器に使うことができ、それらがすべて並列に動作したとしたらど れくらいのことが出来るか」というのを考えてみる。現在のデジタル回路を使っ た計算機である限り、これ以上には速くできないのは明らかであろう。

現在の典型的な CPU で使われているトランジスタ数は 1-2 で ある。これに対し、浮動小数点演算器を構成するのに必要なトランジスタの数 は、パイプライン段数とか細かいことで多少は変わるが、大雑把にいって表 1 の程度である。

  
表 1: 浮動小数点演算器に必要なトランジスタ数

すべての演算を倍精度で行なう意味はほとんどない。というのは、すでに見て きたように軌道の指数的発散の問題があるし、これに影響するのは近傍の粒子 との相互作用、すなわち桁落ちがもっとも大きく効くところである。従って、 最初の座標の引き算と、最後の力の累算は倍精度で行なうとしても、それ以外 は単精度で行なうというのがまあ妥当な選択であろう。これよりも精度を落し てもほとんどの問題であまり問題はないが、しかしそれほど劇的には演算器が小さくな らないからである。なお、理論的には演算器の大きさは乗算器の場合仮数のビッ ト数の2乗に比例するが、仮数が小さくなると指数のところの回路の大きさが 無視出来なくなる。加算器の場合は、高速な加算器やシフタを使った時に大き さがビット数 p に比例する。

まあ、ここでは、倍精度乗算器を使わないとすれば典型的な演算器の大きさは トランジスタと思っておいていいであろう。ということは、 今の典型的な CPU 程度の LSI には最大で 500 もの演算器が入るということ である。 これが 1 GHz のクロックで動くとすれば、チップ単体での性能はピー クで比べて汎用の CPU の 250 倍ということになる。この 250倍というところ に、専用回路の潜在的なアドバンテージがあるということになる。

さて、専用回路(LSI)を作った時の上限の性能がこれくらいであったとして、 このアプローチが現実的なものであると主張するためには以下のような疑問に 答える必要があろう。

以下、順に見ていくことにする。

まず、実際に500も演算器があるチップを作れるのかということだが、これは 以下の3つの手法の組み合わせによって可能になる。第一は、相互作用計算自 体をパイプライン化したハードウェアを作ることである。つまり、相互作用計 算の演算の手順通りに演算器を並べてしまうのである。概念を図 14に示す。このパイプラインは基本的にクロッ クサイクルごとに新しい粒子データを一つ受けとり、パイプラインをずーっと データが流れていって求まった相互作用を最終段のアキュムレータで集計して いく。各演算器が何をするかは始めから決まっているので、このパイプライン のなかで制御するべきものは最終段のアキュムレータだけである。これは、計 算を始める時にデータをクリアするのと、それから計算が終ったところで結果 を固定するというだけで、電子回路としては標準的なロードとクリア機能がつ いたフリップフロップでアキュムレータのレジスタを構成して、それを外から の制御信号で制御するというだけである。このパイプライン化によって、30個 程度の演算器を集積することができる。エルミート公式を使うために時間導関 数も計算することにすれば、50個程度の演算器が集積可能である。このパイプ ラインの基本的な使い方としては、外づけのメモリとメモリアドレス生成回路 (といっても、実質は単なるカウンタ)と組み合わせて全体回路を構成する。 計算の手順としては、最初にメモリに力を及ぼすほうの粒子を書き込み、それ からチップ内の力を受ける方の粒子座標レジスタに値を書き込み、カウンタを 回して計算させ、終ったところでアキュムレータから値を読み出す。これを、 力を受ける粒子の数だけ繰り返す。

  
図 14: 重力計算用パイプライン

第二は、複数パイプラインを同一チップに集積することである。複数のパイプ ラインは、それぞれ別の粒子への力を計算することにしておく。こうすれば、 計算中にチップに供給する必要があるデータ量はパイプラインの本数に無関係 になるので、沢山パイプラインを入れても効率が下がることはない。これは通 常の計算機では並列化が一般には難しく、そのために半導体技術が進んで大量 のトランジスタが集積可能になっても、チップに集積された演算器の数がほと んど増えていないのとは大きな違いである。もちろん、これは使い道を限って 並列性が高いものにしているからである。

第三は、我々が仮想マルチパイプラインと呼んでいる方法である。この方法で は、一つのパイプラインを時分割で使い、複数の粒子への力を計算する。なぜ こんなややこしいことをするかというと、メモリからのデータ転送の必要バン ド幅を減らすためである。つまり、一つのパイプラインが複数の粒子への力を 計算するということは、メモリからの粒子データはその分だけ同じものを使い 回せるということになる。バンド幅がどの程度になるかをここで紹介しておく と、位置を64ビット、質量を32ビットで表現するとして、パイプラインのクロッ ク速度が 100 MHz とすると、単純な設計では転送速度が 2.8 GBytes/s とな る。これは最新のマイクロプロセッサでのプロセッサとレベル2キャッシュの データ転送速度に相当し、実現するのは不可能ではないがかなり特殊な技術を 使った実装が必要になる。マイクロプロセッサのように大量に製造する場合に は量産効果でコストを下げられるが、数が少ないという前提の下ではあまり特 殊な技術を使うのは現実的ではない。ここで、仮想マルチパイプラインの多重 度を例えば 8 とすれば、データ転送速度は 350MB/s 程度となってそれほど難 しい話ではなくなる。パイプラインクロックが 100 MHz というのは最近のマ イクロプロセッサの 1.5 GHz とかいうのに比べてずいぶん低いではないかと 思うかもしれない。これは純粋にチップ設計にどれだけの時間、人とお金を掛 けられるかによっている。

次に、お金の話である。チップを作るのにはかなり莫大なお金が掛かるのは確 かである。が、上でもちょっと出たがあまり極限的な設計をしなければマイク ロプロセッサで掛かるような 100億円単位というわけではない。我々が作った ものの例をあげると、 GRAPE-4 の場合で試作コスト(サンプルチップまで) が2500万円、GRAPE-6 ではだいぶ高くなって1億ちょっとになっている。これ は半導体技術が GRAP-4 の から と進歩した 分初期費用が上がったためと、LSIの詳細設計を基本的には外注にしたためで ある。量産時のチップ単価は GRAPE-6 で3万円程度で、まあまあ我慢出来る範 囲であろう。

次に、誰がそんなチップを設計してくれるのか?ということだが、これは、結 局、予算さえあればどこでもやってくれる。仕様書レベルのものから論理設計 を全部やってもらうとしても、論理設計と検証は優秀なエンジニアが一人でやっ たとして一年程度で出来る。優秀なエンジニア一人、1年間のコストというの がどれくらいかというのは評価が難しいところだが、大企業に頼むとまあ 1 億かもうちょっと高いあたりになろう。これはもちろんそのエンジニアが そんな給料をもらっているというわけでは(日本の企業では)ない。 が、まあ、優秀なエンジニアあたりの企業のトータルコストというのは、良く てこれくらいであろう。悪いともっとかかる。

とはいえ、これは、仕様書と検証用テストデータおよび出力はこちらで準備す るとしての話である。つまり、論理設計そのものは自分ではやらないとしても、 パイプラインやそのまわりの制御回路の動作仕様はこちらで決める必要がある。 とはいえ、高エネルギー実験や天文観測の研究グループが日常的にやっている レベルの回路設計・開発に比べればたいして難しいものではないので、そのあ たりの知り合いとかがいれば教わりながら、、、というので出来る程度の話で はある。

もちろん、これはある程度というよりは巨額な予算をとってきた場合の話であ る。読者のなかでこの記事を読んでそういうことをやってみるのも悪くはない と思う人がひとりくらいはいたとして、予算が1億円からないと始まらないと いうのではもちろん話にならない。

我々が専用計算機を作ろうというのを始めたのは 1988年で既に 13 年も前の ことになる(なんか嫌な感じ)。当時は、半導体技術、基本的には集積度のレ ベルが今とは2桁以上違っており、例えば完全パイプライン動作の 64ビットの 浮動小数点演算器を1チップで実現したLSIは存在しなかった。この時点では、 市販の浮動小数点演算LSI、あるいは固定小数点のALU や ROM をつかったテー ブルルックアップを組み合わせて演算パイプラインを組み上げて、 ラッピングボード1枚で当時の汎用マイクロプロセッサの 100 倍以上の性能を 実現することができた。当時の汎用マイクロプロセッサの 100 倍以上の性能 というのは、いうまでもなく現在のマイクロプロセッサ程度の性能ということ であり、このような市販の演算用LSIを組み合わせる方法ではそれ以上の速度 を実現するのは難しい。これは単にボード上にLSIを並べた時にその間のデー タ転送を速くするのは大変だからである。また、それ以前に汎用の演算 LSIというもの自体が軍用の特殊なものなどを除いてほとんど存在しなくなっ ている。もともとは演算LSIを使って作られていた回路のほとんどが、高速を 必要とする場合には専用のカスタムLSI、半導体技術の進歩によってそれほど の速度を必要としなくなった場合には汎用のCPU や DSP によって置き換えら れたからである。もともと、このような演算専用LSIのかなり大きなマーケッ トはマイクロプロセッサに付加する浮動小数点演算ユニットであった(WTL 1264 など)。マイクロプロセッサが浮動小数点演算回路も集積した時点で、 マーケットのほとんどが失われたのである。

従って、現在において小規模の専用計算回路を設計・開発し、しかも意味があ る程度の性能を実現するのは 10 年前に比べてある意味でははるかに困難になっ ている。 10 年前は、一つの Window of opportunity であったことは間違い ない。というのは、さらにその 10 年まえ、つまり 20 年前を考えると、演算 器を並べて専用演算回路を作る以前に、4ビットとか 8 ビットの乗算器 LSI を並べて乗算器を作る必要があったわけであり、ある程度複雑な専用演算回路 を作り上げることは不可能ではないもののかなり大がかりで費用も人手もかか るプロジェクトとなったからである。

さて、では現在ではもうどうしようもないのかというと、これは実は必ずしも そうではない。 FPGA (Field programmable gate array)というものが半導 体技術の進歩につれて順調にその規模、速度を拡大してきているからである。 FPGA というのはその名前の通り、基本的には後から論理回路を書き込む(プ ログラムする)ことのできる汎用のLSIといってよい。どうやってプログラム 可能にしているかというと、基本回路には小さな SRAM を使い、テーブル参照 によって任意の論理式を実現する。さらにその基本回路間をまたプログラム可 能なスイッチと配線マトリックスでつなぐことで、小さな基本回路では実現不 可能な大規模・複雑な論理を実現するわけである。

プログラムするべきハードウェア回路は、通常のカスタムLSIと同じように VHDL などのハードウェア記述言語によって設計できる。多数ビットの加算器 や乗算器は、多くの場合にFPGAメーカーが提供する VHDL コンパイラによって 自動的に最適な回路構造に変換される。 (速度優先か回路サイズ最小化優先 かをある程度の範囲で指定できるものが多い)

回路規模でいうと、現在はメーカー公称で100 万ゲートを超えるものが入手可 能(高価ではあるが)になっている。ここでゲートはカスタム LSI や FPGA の大きさを表示する時の伝統的な単位で、NAND 回路一つ分を 1 ゲートと数え る。トランジスタ数としては CMOS の場合4トランジスタにおおむね相当する。 動作速度も、設計に注意を払えば 100 MHz 近くまで到達できなくもない。

まあ、100万ゲートというのはメーカーさんの主張で、実際にはその半分程度 の規模の回路が入れば万歳だが、それでも 200 万トランジスタ、したがって 単精度なら浮動小数点演算器が100個程度入ることになる。これは上に書いた カスタムLSIの数字に比べるとかなり小さいが、これはもちろんプログラム可 能性のために無駄が発生するからである。逆に、その差が5倍程度なのは、一 つにはカスタムLSIの数字が少し古いからではある。この原稿を書いているの は 2001年2月である。 FPGA の数字は現時点のものだが、カスタム LSI の数 字は我々が実際に LSI を作った 1998年時点のものになっている。そのために FPGA のほうが技術進歩の分の数倍の得になっている。もう一つは、カスタム LSI のほうは技術的な極限まで大きいものを考えてはいないからである。

100 MHz で演算器が 100 個だと 10 Gflops にしかならないわけで、現在の最 高速のマイクロプロセッサに比べて10倍しか速くないということにはなる。し かし、それでも以下に述べる理由からマイクロプロセッサを 10 個並べるより は有利になる可能性が高い。一つは、カスタム LSI の場合と同様に、 100 個 の演算器を使ってパイプラインを構成し、繰り返し演算の間はそのすべてを有 効に使うというのが比較的容易に出来るからである。さらに、チップ1-2個の システムで実証が出来れば数百万円の予算で100チップ程のシステムが構成出 来るであろう。これはマイクロプロセッサ 1000台に相当し、例えば通常の大 学の研究室レベルでは購入予算もそうだがスペース・消費電力からも実現困難な システムになる。

さて、もちろん、FPGAでは試作コストは掛からないといっても、中身を設計し ないといけないのはカスタム LSI と同じである。固定小数点の演算器程度な らばコンパイラが作ってくれるし、ソフトウェアを買ってくれば浮動小数点演 算器もコンパイラが作ってくれるようになる。従って、どういう計算をさせた いのかこちらにわかっていれば、中身の設計は(面倒ではあるが)それほど難 しいという話ではない。

これはもちろんカスタム LSI の場合でも全く同じことであるが、むしろ面倒 なのは次の2つ、つまり一つは動作検証、もう一つは周辺回路の設計である。

動作検証は、要するに「この回路でいい」かどうかを確認するわけである。例 えば、データ形式を普通の 64 ビット浮動小数点と決めてしまえば、回路規模 は大きくなるが検証は簡単になる。普通にPCとかで計算した答と合わせるだけ だからである。しかし、これは多くの場合に回路が大きくなりすぎるので、や はり普通とは違うデータ形式を使いたくなる。 例えば 32ビット浮動小数点演算を使うなら、64ビット型との変換をちゃんと やればやはり普通にPCとかで計算した答と合わせるだけで済み、検証は比較的 簡単である。

とはいえ、実際にはいろいろ悩ましいことはある。大抵の場合に全演算に32ビッ ト精度がいるわけではないし、またそもそも浮動小数点がいるところは少ない。 こういったところを、必要に応じて最小限のデータ長や形式を選んで実現して いくことによって、全体を 64 ビットでやったりあるいは 32ビットしか使わ なかったりする場合に比べてかなり回路規模を小さく出来る。しかし、そのた めにはそのような演算器を作らないといけないし、もっと大変なことにそのよ うな演算器を使って作ったパイプラインからでてくるべき答を知らないといけ ない。 さらに、その精度が必要を満たしているかどうか、つまり、そのパイプライン でシミュレーションしてまともな答が出るかどうかもチェックする必要がある。

とはいえ、FPGA の場合は、これらすべてを実際に回路を作ってからやるとい うのも可能ではある。カスタム LSI の場合には作り直しには莫大なお金と時 間が掛かるので、十分に時間をかけ、まず汎用計算機上でパイプラインのシミュ レータをつくってそれを単体テストし、さらに実際の N 体シミュレーショ ンプログラムに組み込んで評価する。その評価が終った後で、実際のハードウェ アを設計しその動作がシミュレータと同じかどうかをテストデータで検証する というかなり回りくどい手順が必要になる。この検証に、設計よりもはるかに 長い時間がかかるといっても過言ではない。 FPGA ではこのへんはある程度い い加減にやれる。

周辺回路のほうは、設計自体が面倒ではある。例えば、 FPGA の石を買ってき たとして、それに少なくともメモリをつけないといけないしホスト計算機に組 み込むためのインターフェースもいる。これらをまとめてプリント基板を設計 し、どこか業者に製造を依頼しないといけない。

とはいえ、例えば標準的な PCI インターフェースと外づけメモリに大規模 FPGA を組み合わせたようなものならば、国内でもいくつかの企業が商品化し ているものがある。そのようなものを使ってみるのがもっともハードルが低 い方法であろう。

最後に、並列化について簡単にまとめる。LSI一つでマイクロプロセッサの100 倍以上の性能が実現できるとしても、一つしか使えないならマイクロプロセッ サ100個、というよりはPC を 100 台並べたクラスターの方が速いことになる。 もちろん、既に見たように重力多体問題に対して 100台を超える規模の PC ク ラスタで十分な性能を出した例というのはあまりないが、まあ、原理的にはそ ういうことではある。

重力多体問題専用計算機、というより相互作用計算専用回路の利点は、かなり 大規模な並列性が比較的容易に実現出来るということである。もちろん、これ はやりたい計算自体に高い並列性があるからである。アルゴリズムによっては 次節以降に述べるようにいろいろややこしい話がでてくるが、単純な直接計算 で一定時間刻みの場合には、数万から100万以上に及ぶ粒子から、またそれだ けの数の粒子への力を計算するわけで、力を受けるほうはすべて独立、つまり 並列に計算出来るしまた力を及ぼすほうも、最終的に合計するところの計算時 間が問題にならない範囲では並列に計算できる。

この並列性は既に単一のチップの話をする時に利用しているわけで、パイプラ イン化できるのは力を及ぼす方がある程度多いからであり、また複数パイプラ インを並列動作させられるのは力を受ける側が複数あるからである。パイプラ イン段数はせいぜい 100 段程度なので、粒子数が100万とかいっているときに はほとんど無視できる。従って、一つの粒子への力をさらに 100 とか 1000 チップにわけで、並列に計算してもパイプライン遅延は問題にならない。もっ とも、1000 もパイプラインがあるとそれを順番に合計していったのでは時間 が掛かり過ぎるが、これはトーナメント方式で合計すれば2入力の加算器を使っ たとして 10 段、 4 入力のものを特別に設計すれば 5 段のツリー構造ネット ワークで済む。さらにこの全体がパイプライン化できるので、段数が多くても 性能にはそれほど影響しない。 こういった回路は、規模的には現在の FPGA に容易に収まる。

複数粒子への力を並列に計算するほうはもっと簡単で、単にメモリとアドレス 生成ユニットの回路に複数のチップをぶら下げるだけである。

適当な大きさのプリント基板を作れば、チップ数十個くらいは載せられる。こ れで、普通のPCクラスタとしては 1000 台程度の規模(大規模な計算機センター にある程度)とほぼ同等の性能が実現できる。それ以上の数を並べたければ、 プリント基板を複数並べればいい。

並べたものをどうやって使うかというような話になると、段々問題というより もアルゴリズム固有の側面が大きくなる。これらは以下の節でもうちょっと議 論する。

8.2 独立時間刻みと専用計算機

すでに繰り返し述べてきたように、専用計算機がいかに高速であっても、ある 程度賢いアルゴリズムが使えなければ実用的な意味はほとんどない。まず、こ こでは、独立時間刻みについて考える。

独立時間刻みのために必要な機能を専用計算機に作り込むのは実はそれほど大 したことではない。普通の場合となにが違うかというと、単に重力を計算する 時に位置がそのまま使えるわけではなくて、予測子を評価してから重力を計算 するというだけである。予測子は単なる時間の多項式なので、汎用計算機で計 算するのと全く同様にホーナーの公式を使って高い次数から順に計算していけ ば良い。

とはいえ、一つ問題なのは同時に積分できる粒子の数である。もともとの独立 時間刻みでは一度に一つの粒子しか積分しない。これでは、多数パイプライン があってもそれらを一つの粒子への重力を計算するのにしか使えないことにな る。但し、すでに並列化のところで述べたような、ブロックステップ法を使え ばそこそこの数の粒子を並列に積分出来るので、これはあまり大きな問題では ない。[Mak91a]

我々が開発したシステムの中では、 GRAPE-4[MTES97]および GRAPE-6 が独立時間刻み用(といっても、もちろんそうでない計算コードでも 使えるが)に設計されている。 GRAPE-4 は 1992 年に開発を始めて 1995 年 に完成したシステムであり、カスタム LSI 一つ(HARP chip)に重力とその時間 導関数を計算するパイプラインを入れ、もう一つのカスタムLSI(PROMETHEUS chip)に予測子計算を入れた。HARP チップは2粒子への力を仮想マルチパイプ ラインにより並列に計算する。

GRAPE-4 では、プリント基板には PROMETHEUS チップ 1 つと HARP チップ 48 個を載せた。計算中は HARP チップは PROMETHEUS チップからすべて同じデー タを受けとるので、 個の粒子への力が並列に計算される ことになる。 GRAPE-4 全体はこのボード(プロセッサボード) 36 枚からな る。これを 9 枚ずつラックに入れ、それらをコントロールボードと呼ばれる 別の基板から制御する。プロセッサボードで求まった力はコントロールボード 上でハードウェアロジックで加算され、4枚のコントロールボードで求まった 力は最終的にホスト計算機上のソフトウェアで合計する。GRAPE-4 のピーク速 度はおよそ 1.1 Tflops に相当する。なお、動作クロックは 32 MHz である。

GRAPE-6 は 1997 年に開発を始め、現時点で小規模なシステム(4 Tflops 相 当)が動作しており、大規模システムの組み立て中である。GRAPE-6 では、 GRAPE-4 に比べてはるかに大きな LSI が利用できるようになったため、 1 チッ プに重力計算パイプライン 6 本と予測子パイプラインを組み込んだ。 一本の重力計算パイプラインは 8本の仮想パイプラインとして動く。これによ り必要なメモリバンド幅を減らしている。その代わりに、1チップですでに48 粒子への力が同時に計算されてしまう。複数チップでこれがさらに増えるのは 嫌らしいので、チップごとにメモリをつけて、ボード上でも複数ボード間でも 違うチップからの力は合計する。但し、現在のところ普通に合計することを考 えているのは最大 512 チップまでで、それ以上はチップ数に比例して並列に 計算される粒子数が増えるような使い方をする予定である。

なお、予定としては 4096 チップを並列動作させ、 130 Tflops を実現したい なあということになっているが、予算的にちょっと苦しい。

さて、独立時間刻みができれば次はネイバースキームも、、、というのが人情 であるが、いまのところ良い方法は見つかっていない。問題なのは、ネイバー リストが粒子ごとに違うということである。このために、ネイバーからの力を 並列計算するためには違う粒子への力を計算するパイプラインは、力を及ぼす 粒子も違うということになる。このために、ネイバーからの力を高速に計算す ることが出来るハードウェアは、どのように設計しても本質的にメモリバンド 幅が大きくなる。ところが、メモリバンド幅は製造の困難さに大きく影響する ので、そういうものをつくるのは難しいわけである。

8.3 ツリー・FMMと専用計算機

ツリー法や FMM でも、計算量の多くを占めるのは粒子同士の相互作用の計算 である。まあ、高次の多重極展開を使うと単純には粒子同士の相互作用という ことにはならないが、既に述べた疑似粒子法を使えばツリー法の場合は結局計 算するのは粒子からの重力になる。FMMの場合でも、もっとも計算量の多い多 重極展開を局所展開に変換するところは、疑似粒子法とアンダーソン法を組み 合わせることで粒子からの重力を計算するだけで済むようになる。以下、主に ツリー法の場合の話をする。FMMは専用計算機がどうこうという以前にそもそ も使われていないからである。

実際にはツリー法と専用計算機を組み合わせて高い性能を出すのに はもうちょっといろいろ工夫が必要である。ネイバースキームの場合と同様に、 ツリー法では重力を受ける粒子ごとにすべて違った計算をする。しかも、ある 意味ネイバースキームよりもさらに面倒なことに、どのような粒子(節点) からの力を受けるかはあらかじめ決まっているわけではなくて、実行時にツリー 構造をたどっていくことできまる。そういうものをそのままハードウェア化す るのは、不可能ではないがやはりパイプライン1本で粒子一つしか計算出来な いとか、そういう問題が起きる。

実はベクトル化や並列化の場合でも、規模は小さいながら類似の問題が起きて いて、その問題に対応するようなアルゴリズムも提案・実現されている。ここ での問題は、ベクトル化は出来るにしてもベクトル化されたループの中に条件 分岐(ツリーを降りるか降りないか)があり、またメモリアクセスも連続では ない(ツリー自体はポインタを通してアクセスされる)ためにどうしてもハー ドウェアのピーク性能からみると低い性能になってしまうということである。 典型的には、ベクトル計算機の場合理論ピーク性能の 10% 程度の性能になっ てしまう。

仮にある程度計算量が増えるにしても、力を計算するループのなかから条件分 岐や間接アクセスを追放できれば理論ピーク性能の半分以上は出せるので、そ ういう方法はないかというのが問題であった。これに対する一応の回答は、 「あらかじめある粒子に力を及ぼす粒子とセルのリストを作っておいて、それ を連続な配列にしまっておく」というものである。こうすれば、少なくとも重 力を計算している間はハードウェアのピークに近い性能がでる。

とはいえ、これは問題をすり替えただけなわけで、じゃあリストは誰がいつ作 るんだというのがもちろん問題である。リストを各粒子用に作るわけなので、 そのための計算量は重力計算自体の計算量とたいして変わらない。こちらが今 度は時間が掛かってしまう。

ここで Barnes が考えたのは、リストのほうを複数の粒子で使い回そうという ことである。もちろんこのためにはリストが複数の粒子のどれにとっても良い ものでないといけない。このためには、リストを作る時のツリーをたどる判定 条件をすこし厳しくして、複数の粒子のどれか一つについてでも条件が満たさ れなければ下におりるということにすればよい。

実際には物理的に近くにある複数の粒子は、ツリー構造を使って適当な大きさ のセルの中にある粒子ということにするのが簡単である。この時には、中の粒 子を見なくてもセルの幾何学的な形状から判定条件を決めればいい。つまり、 普通に相手との距離を計算する代わりにセル内の任意の点から相手への最短距 離を求めればいい。

セルの中にある粒子数を大きくとると、ツリーをたどる回数はいくらでも減ら すことができる。しかし、こんどはリストの長さのほうが伸びていくので、重 力計算の計算量がだんだん増えてしまう。従って、どこかに計算時間が最小に なるようなセルの大きさがあることになる。この最適点はもちろんツリーをた どる計算の速さと粒子間の重力を計算する速さの比によっていて、粒子間の重 力計算が速いと粒子数を大きくしてツリーをたどる回数を減らすことができる。 つまり、 GRAPE のほうが速ければある程度まではホスト計算機側の仕事を減 らせるということになる。

実際にどの程度の性能かというと、少し古いデータで恐縮だが GRAPE-5[KFMT00] 2枚にホスト CPU として EV6 600 MHz の Alpha を組み合わせたもの[KFM99]で、 Warren と Salmon によ る並列ツリーコードを 140 台の EV5 600 MHz Alpha 上で実行したもの [WGL+98]の半分程度の性能が出ている。言い替えると、 GRAPE-5 2枚で Alpha 70 台分の性能ということになる。 600 MHz EV5 Alpha というと若干旧式だが、名目ピーク性能では 1.2 Gflops で最近の Intel プ ロセッサとさして変わらない。さらに、Warren と Salmon のコードは非常に 高速なものであり、他の同様なコードに比べて 2-3倍は速く、理論ピークから そんなに遠いわけではない性能が出ている。

従来は、GRAPE では粒子からの力しか計算出来ないということで、モノポール (重心において、重力には符号がないので精度としてはダイポール)しか扱え ず、高精度の計算では計算量が急速に増大するという問題があったわけだが、 疑似粒子法によってその問題もほぼ解消した。

一層の高速化のためにはホスト計算機側を並列化する必要がでてくるが、これ は原理的には既存の並列プログラムを手直しして GRAPE が使えるようにする という作業になる。これまで、そのような並列実行環境が実際にはなかったた めに現時点ではホストが並列でGRAPEを使うプログラムというものは存在しな いが、この1月に国立天文台に16台の Alpha クラスタのそれぞれに GRAPE-5 をつける並列システムが導入されたことを契機に開発が進むことを期待してい る。これにより、従来は1000台規模の超並列システムでのみ可能であった1億 粒子を超える大規模シミュレーションが、比較的小規模なクラスターで実行可 能になる。

8.4 すべてを組み合わせる?

さて、専用計算機 GRAPE と独立時間刻みをうまく組み合わせることが出来る ことはわかった。また、 GRAPE とツリー法もまあまあうまく組み合わせてそ れなりの性能を出すことができる。となると、このすべてを組み合わせること はできないかという希望が出てくることになる。ツリー法の項で述べたように、 ツリー法と独立時間刻みを組み合わせるというのは一応出来ている話である。 さらに GRAPE を使うのも、それぞれに対しては出来ている。ので、全部組み 合わせるというのもなんとかなりそうだが、今のところうまくいっていない。 まあ、単に誰もそんな方法のことを真面目に検討していないからというのが最 大の理由なので、いろいろ考えればなんとかなるのではないかとは思う。まだ、 ここで紹介出来るほどまとまってはいないが、基本的には実装できそうなアイ ディアはあるので1-2年のうちにはなんとか出来ると個人的には思っている。



next up previous contents
Next: 9 まとめ Up: 重力多体系の数値計算 Previous: 7 計算法 --- 空間領域



Jun Makino
Sat Feb 24 22:13:53 JST 2001