./note033.html

 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 日に計算できる。

とあります。

ToC