| Previous | ToC | Next |
GRAPE (GRAvity PipE) は、私達が 1988 年頃から開発を続けてきている一連 の重力多体問題専用計算機です。どういうもので、なんのためにやっているか、 ということを、私の個人的な視点からまとめます。ここは事実を客観的に書く べきところですが、当事者である以上そんなのは無理なので書く本人の意志と しては客観的であろうとしていますが読む人はそんなことはないかもと思って 読んで下さい、というものです。
もうひとつは、これも相互作用が重力であることによるのですが、断熱壁や等 温壁で囲まれた密度一様なものを考えるわけではなくて、自分自身の重力でま とまっているものを考える必要があるということです。
この話と計算機との関係は、こういった系の研究ではなるべく沢山の粒子を使 いたい、ということと、でも粒子の数を少し増やすと計算量がどーんと増える ので計算は急速に大変になる、ということです。
重力多体系ではそう簡単ではないところがあります。この辺は少し難しい話に なりますが、大雑把にでも説明しないと何故そもそも速い計算機が欲しいか、 というのがわからないと思いますので少しお付き合い下さい。
杉本はこの振動に「重力熱力学的振動」という名前をつけていたのですが、こ の研究会をきっかけにこの言葉は世界に受け入れられていくことになりました。
以下に当時書いた Cyber 205 用のコードを示します。
C***
do 40 i = 1, 3
maxnod2=maxnode
dxs(1,i;npnnul) = poscp(1,i;npnnul)
& - q8vgathr(pos(1,i;maxnod2), pnodes(1;npnnul); npnnul)
40 continue
stols(1;npnnul) = q8vgathr(sizetol2(1;maxnod2),
& pnodes(1;npnnul); stols(1;npnnul))
masss(1;npnnul) = q8vgathr(mass(1;maxnod2),
& pnodes(1;npnnul); masss(1;npnnul))
nexts(1;npnnul) = q8vgathr(nxtpos(1;maxnod2),
$ pnodes(1;npnnul);nexts(1;npnnul))
childs(1;npnnul) = q8vgathr(fstchd(1;maxnod2),
$ pnodes(1;npnnul);childs(1;npnnul))
rsqs(1;npnnul) = dxs(1,1;npnnul)**2 + dxs(1,2;npnnul)**2
$ + dxs(1,3;npnnul)**2
bstom(1;npnnul) = (rsqs(1;npnnul) .ge. stols(1;npnnul))
where (bstom(1;npnnul))
npnods(1;npnnul)=nexts(1;npnnul)
otherwise
npnods(1;npnnul)=childs(1;npnnul)
endwhere
where ((.not. bstom(1;npnnul)) .or.
& (ptclcp(1;npnnul) .eq. pnodes(1;npnnul)))
& masss(1;npnnul) = 0.0
rsqs(1;npnnul) = 1.0/(rsqs(1;npnnul) + eps2)
phims(1;npnnul) = masss(1;npnnul) *
& vsqrt(rsqs(1;npnnul); npnnul)
phicp(1;npnnul) = phicp(1;npnnul) - phims(1;npnnul)
phims(1;npnnul) = phims(1;npnnul) * rsqs(1;npnnul)
do 90 i = 1, 3
acccp(1,i;npnnul) = acccp(1,i;npnnul) -
$ dxs(1,i;npnnul) * phims(1;npnnul)
90 continue
bnewcp(1;npnnul)= (npnods(1;npnnul) .ne. null)
j=q8scnt(bnewcp(1;npnnul))
if(j .lt. npnnul) then
nend=npnnul-j
endlst(1;nend)=q8vcmprs(pindex(1;npnnul),
$ .not. bnewcp(1;npnnul);endlst(1;nend))
do 230 i=1,nend
k=endlst(i)
pbody=ptclcp(k)
phi(pbody)=phicp(k)
termcnt(pbody)=termcp(k)
acc(pbody,1)=acccp(k,1)
acc(pbody,2)=acccp(k,2)
acc(pbody,3)=acccp(k,3)
230 continue
if(j .ne. 0) then
ptclcp(1;j)=
$ q8vcmprs(ptclcp(1;npnnul),bnewcp(1;npnnul);j)
phicp(1;j)=
$ q8vcmprs(phicp(1;npnnul),bnewcp(1;npnnul);j)
termcp(1;j)=
$ q8vcmprs(termcp(1;npnnul),bnewcp(1;npnnul);j)
do 240 k=1,ndim
poscp(1,k;j)=
$ q8vcmprs(poscp(1,k;npnnul),bnewcp(1;npnnul);j)
acccp(1,k;j)=
$ q8vcmprs(acccp(1,k;npnnul),bnewcp(1;npnnul);j)
240 continue
endif
endif
if(j .ne. 0) then
pnodes(1;j)=q8vcmprs(npnods(1;npnnul),
$ bnewcp(1;npnnul);pnodes(1;j))
endif
npnnul=j
こちらは、アルゴリズムの基本は同じですがあんまり邪悪なチューニングをす
る前の基本的なプログラムです。ループの入れ替えをした後の最内側ループを
示しています。
C***
do 40 i = 1, pcount
pbody = pfirst + i - 1
if (pnode(i) .ne. null) then
p=pnode(i)
dx = pos(pbody, 1) - pos(p, 1)
dy = pos(pbody, 2) - pos(p, 2)
dz = pos(pbody, 3) - pos(p, 3)
rsq = dx**2 + dy**2 + dz**2
if (rsq .ge. sizetol2(p) .or. p .eq. pbody) then
if(p .ne. pbody) then
rsq = rsq + eps2
r2inv = 1.0/rsq
rinv = sqrt(r2inv)
phim = mass(p)*rinv
phi(pbody) = phi(pbody) - phim
phim = phim * r2inv
acc(pbody, 1) = acc(pbody, 1) - dx*phim
acc(pbody, 2) = acc(pbody, 2) - dy*phim
acc(pbody, 3) = acc(pbody, 3) - dz*phim
endif
pnode(i)=nxtpos(p)
else
nodediv = nodediv + 1
pnode(i)=fstchd(p)
endif
endif
40 continue
Cyber 用のコードは3倍くらいの長さになった上、なんだか意味がわからない
ものである、というのはわかるかと思います。これはコンパイラの能力の問題、
ということもできなくもないのですが、下のプログラムから上のプログラムに
対応するようなアセンブラ、というか上のプログラムは q8xxxx というのはそ
のままベクトル命令に変換されるのでほぼアセンブラなのですが、を出すこと
をコンパラに期待するのはなかなか現実的ではない、というかコンパイラを書
く人にそんなことをやってというのには無理があるかもしれません。
まあ、とにかく、そういうわけでツリーコードのベクトル化も結構上手くでき ました。そこで、 CM でもやってみよう、という話になって、 87年の秋から 冬にかけて 3ヶ月ほどボストンの MIT AI Lab の近くにある TMC で机をもらっ てました。
TMC は、数値計算機のための計算機に何をもたせるべきか、というのは全然考 えずに人工知能屋さんがハードウェアとソフトを作っていた会社でした。 86 年にもバーンズが CM 上にツリーコードの実装をしたのですが、その時に使え た言語は *Lisp という、 Lisp に CM 上のアセンブラ呼び出し機能をつけた だけのものでした。
但し、元々 Hillis は人工知能屋で、言語設計には G. L. Steele 等の人が入っ ていたこともあり、アセンブラといっても上の FTN200 なんかに比べると全く 次元が違うものでした。CM のハードウェア、ソフトウェアの画期的な点は、 仮想プロセッサの概念を導入したことです。それまでの SIMD 並列機や、ある いはベクトル機をアセンブラでプログラムする場合もそうですが、ハードウェ アの並列性、つまり何個プロセッサがあるか、というのがプログラマに見えま す。このため、例えば 8192 粒子のシミュレーションをするプログラムを書く 時、使えるプロセッサが 256個だったとすると、 1 プロセッサが32個の粒子 を分担することになって、そのための配列宣言とかループ構造をももったプロ グラムを書くとかそういうことをする必要が起きます。 これに対して、 CMは 1つの物理プロセッサが複数の仮想プロセッサの実行をエミュレーションする、 というのを SIMD の実行レベルで行う、という仕掛けがありました。このため、 上の例では、アセンブラで書く時でも 8192 個のプロセッサを使うことができ ます。これはプログラムの生産性を非常に上げるものです。
なんだか話がそれていきますが、とにかくこの時は CM 上の C* という、C++ のクラスの機能を使って C にSIMD 並列実行機能をつけた言語でツリーコード とかを書いてました。これは、上の仮想プロセッサのおかげでかなり複雑な処 理が比較的簡単なプログラムで実現でき、まあまあの性能がでました。 ここでまあまあというのは理論ピークの10パーセントとかそれくらいです。ベ クトル機でもそんな程度なので、価格性能比の高い CM では十分競争力があり ました。
但し、これは性能評価くらいにおわって、あんまり実際のシミュレーションに は使いませんでした。シミュレーションには、プリンストンにはいった ETA-10 とか、当時リクルートがもっていたスーパーコンピューター研究所の Cray XMP や NEC SX-1 を使っていました。このスーパーコンピューター研究 所はなかなか不思議なところで、とても採算がとれるあてはなかったのではと 思うのですが、 Cray XMP、NEC SX-1、 富士通 VP-400 等の大型スパコンの他、 研究所に Alliant FX-8、 Ardent Titan、PAX-64J、Intel iPSC/2 等の並列機 がありました。大型スパコンは計算サービス部門のものですが、結構研究所の 関係者には使わせてくれました。当時東大の大型計算機センターのスーパーコ ンピューターを使うと1時間 1 万円ほどしたので、これは極めてありがたいも のでした。
と、そんなことをしていた 1988年の秋、指導教官の杉本から4ページほどの研 究会集録原稿のコピーみたいなものを渡されました。当時野辺山宇宙電波観測 所にいた近田の、「計算機を使う、計算機の技術を使う」という原稿です。 近田がその年の天文・天体物理若手夏の学校という、大学院生が中心となっ て組織する勉強会みたいなものに講演者として招待されて話をしたものをまと めたものです。
野辺山宇宙電波観測所は世界で初めて 5 素子のミリ波干渉計を動かしたとこ ろです。干渉計というのは、複数の電波望遠鏡からの信号を計算機の上で合成 することで口径の大きな電波望遠鏡と同じ分解能を実現しようとするものです。 そのための計算機の能力は望遠鏡の台数が多いほど、また観測する波長が短い ほど高いものが必要になり、野辺山の場合は 100Gflops ほどが必要でした。 ベクトル機の能力が 1 Gflops なかった時代の話です。
まともに計算機を買ったのではとても無理なのですが、近田たちは信号処 理用の専用計算機を作ることでこの処理速度を実現しました。出来たものは FX という名前で、基本的には望遠鏡からの信号を FFT (高速フーリエ変換) する専用回路です。データのビット長を FFT の部分は 6 ビット程度とするこ とでハードウェアの規模を小さくし、さらに完全専用回路とすることで 100Gops (浮動小数点ではないので fl をつけない)の速度を普通のメインフレー ム程度の規模の回路で実現しました。
近田の原稿の主旨は、そのようなやり方をシミュレーションでも使えばよ い、というものでした。見出しだけを引用すると
「昨日、今日、明日」
「性能は値段に大体比例する」
「実はオーバーヘッドの多少によって性能/価格は3桁違う」
「効率の良い implement しやすいアルゴリズムを選ぶことが大事」
となっていて、基本的な思想はこれだけでわかるかと思います。但し、杉本に とって大事だったのは、近田が単に例として「こんなハードウェアができます」 とブロック図を書いたものが粒子間の重力を計算する専用回路で、まさに私達 が必要としていたものだったことです。杉本は早速近田と接触し、PAX の慶応 の川合、その時大学院合格が決まっていた伊藤、私の3人と一緒に野辺山に出 かけて近田の話を聞きました。
その時にどういう話をしたか、というのは私は実は殆ど憶えていないのですが、 基本的には伊藤が修士の間に設計する、いろいろわからないことは近田に聞く、 というような話をしたはずです。近田の原稿を引用すると
具体的イメージをつかんでもらうために、例として図2にパソコンにつない でN 体問題の重力を計算する回路の概念を示します。 10MHz クロックで間 口がパラレルのパイプライン式で32bit固定小数点で計算させるとすると、 市販の IC で 200個、0.1 立方米くらいの嵩で、250MOPS のスーパーコン ピューター並みの性能が得られるでしょう。値段としては自作で 400万円 位でしょう。もし億円オーダーの金をかけて、大型計算機につなぐような バンとしたものを作る場合は、高性能を狙って専用 IC を作り、 IC 5000 個、 4 byte 浮動小数点演算、 20MHz クロック、20立方米、 100Gflops、 10万体問題を4万ステップ 1 日に計算できる。とあります。
| Previous | ToC | Next |