The Influence of the Los Alamos and Livermore National Laboratories on the Development of Supercomputing by Mackenzie (1991)

​http://dx.doi.org/10.1109/MAHC.1991.10014

を読みながらのメモ。以下、論文からの引用とコメント：

The failure of suppliers to produce the required integrated circuits, anti-war demonstrations, riots, and fire bombing on the campus of the University of Illinois in 1970, and a range of other circumstances prevented the smooth development of ILLIAC IV.

ILLIAC IVについて。

The requisite calculations for one time-step for one cell might amount to 200 floating-point operations for the early 1960’s code, and 20,000 for the modern one (Livermore interviews).

核兵器のモデリングに必要な演算量について。この間に、モデルに利用されるメッシュも、2次元モデルから3次元モデルが必要になり、大幅に増加している。

Statically, up to 30% of the instructions in a Monte Carlo program may be branches (Michael interview).

メッシュ上で偏微分方程式の解を求める手法は、コードに分岐が少ないのに対して、核兵器モデリングでもうひとつ重要な計算手法であるモンテカルロ法は、本質的に分岐が多い。

One million IBM cards carried the requisite initial values, one card for each point in the computational mesh, and “the computations to be performed required the punching of intermediate output cards which were then resubmitted as input” (Stern 1981).

Los AlamosのFrankelとMetropolisが、最初にENIACで計算をしたモデルのこと。

Neither laboratory has bought one of the Japanese supercomputers, though here questions of nationalist protectionism come into play as well as the questions of the suitability of particular machines for the laboratories’ computational tasks.

Los AlamosとLivermoreは、スパコンの選択にうるさい。速いモデルだからという理由だけで計算機を選ぶことはなかった。当然、日本製のベクトルスパコンは選択されなかった。

The Control Data STAR-100 episode was pivotal, because it secured the commitment of the laboratories to mainstream supercomputing rather than the more massively parallel alternatives.

STAR-100というある意味画期的な計算機は、仕様を満たすことができなかったため、その後LAやLLNLはいわゆる普通の計算機(CDC Cyber205などのベクトル計算機)に強くコミットすることになった。実はこれらの核関連の研究は多くが機密であったため、本当に必要な「核」の部分のコードが公開されず、それをターゲットとしてのアーキテクチャやソフトの最適化が不可能という問題もあったようだ。そのため、LLNLではLivermore Loopsと呼ばれる、LLNLが必要としている計算の典型的なものを代表するベンチマークを用意するようになった。その後、色々な計算機はこの「ベンチマーク」の性能を競う風潮ができた。

The project was formulated in an “almost pathological atmosphere of optimism and its corollary, fear of being left behind” (Bashe et al. 1986); as outlined above, the designers were also trying to satisfy the needs of quite different kinds of users.

IBMのStretchについて。結果としてStretchは735個の命令をもつ計算機になってしまった。この直前に「VAXのアーキテクチャはDEC社内の官僚主義的な構造を反映している」というTom Westの言葉を、"The Soul of a New Machine"から引用している。IBMも同様の問題を抱えていたのだろう。

Livermore’s “Sidney Fernbach was reportedly one of the few people in the world from whom Seymour Cray would accept suggestions” (Lundstrom 1987, p. 90), and Cray took, and takes, care to become aware of the laboratories’ computational needs.

Cray-1は、LAやLLNLの要求にもとづいて設計された計算機ではなかった。実際、Crayは彼らのような核関連の研究だけでなく、NSAが必要としている計算にも配慮して、ある種のバランスのとれたアーキテクチャを開発することができて、そして成功した(Cray-1では)。それでも、この引用にあるように、LLNLによるCrayに対する影響が全くなかったわけではないらしい。"one of the few people"ねぇ。

The social limits of speed: development and use of supercomputers" by Elzen & Mackenzie (1994)

​http://dx.doi.org/10.1109/85.251854

を読みながらのメモ。以下、論文からの引用とコメント：

As one of Cray's famous maxims has it, he likes to start the design of a new-generation machine with "a clean sheet of paper."

CDC6600の成功後、次の計算機の設計方針についてのCrayの見解。Crayは過去のアーキテクチャとの互換は最重要視していなかった。結果的に、CrayはCDCの方針に反対し、CDCを退社してCRIを設立した。

Cray used to say that he knew all his users by their first names.

Crayは常に利用者(この場合アメリカの国立研究所の核関連の研究者達)のことを第一に考えて計算機を設計したらしい。

This was recognized even by IBM, prompting a famous acerbic memo to his staff from chairman Thomas J. Watson, Jr., inquiring why Cray's team of "34 people -- including the janitor" had outperformed the computing industry's mightiest corporation.

CDC6600が発売されたあとの当時のIBMのトップのメモ。"janitor"は守衛。

In 1980, it was unclear that the 200 kilowatts that were generated in such a small and densely packed computer could be dissipated.

3D基板で実装されたCray-2をどう冷却するのかという問題。結局フッ化炭素(液体)を使って冷却することになった。

In a sense, CRI gradually became the "prisoner of its own success." New innovations were constrained by the very socio-technical network that made CRI so successful.

これがこの論文の核心だろう。元々、CRIは、CrayがCDCの開発方針を受け入れることができず独立して設立した会社で、その目的はCrayの好きなように世界最速の計算機を作ることだった。Cray-1はそのもくろみ通りに、当時の世界最高速を実現した。一方でCRIは売り上げを伸ばすために、伝統的な顧客であった国立研究所以外にも、民間の会社への営業をするようになった。このときに"Hero-problems"という、その業界で重要な問題だが時間がかりすぎる問題を、CRIなら自社の計算機を使って高速化できる、というのを売りとして営業をした。色々と曲折はあったが、この戦略は成功し、民間のいろいろな大企業がCRIの計算機を導入した。その結果、それらの新しい顧客は必ずしも世界最速の計算機を求めているのではなく、ソフトウエア資産が継続的に利用可能かどうかや計算機の信頼性に重きをおいたため、CRIは、新しい世代の計算機のたびにアーキテクチャを一新するという、創業者のCray的なやりかたをとることが不可能になった。つまり、CRIは顧客の要望に基づいてハードとソフトを開発するようになったために、顧客はCRIのシステムに依存し、CRIはその呪縛から逃れることが不可能になった。これが90年代前半の状況。

"socio-technical network"とは、Crayとそのスパコンの利用者たちの関係性という意味での人間関係のつながりの意味と、上に書いたようなCRIの計算機とその顧客たちのとの強い関係性、さらにその影響下にある計算機のアーキテクチャやそのためのソフトウエア(OSやコンパイラ、ライブラリ、アプリまで)などの技術的な関係性、以上を総称した概念。技術的に世界最高速を目指すという、古き良きCrayら開拓者たちによるラディカルな計算機開発の時代は、技術的ではない制約を受けるようになり、1980年代の終わり頃には終焉しつつあった。それらの制約を避けて、新たな計算機を開発するには、CrayやSteve Chen(Cray Y-MPのアーキテクト)のように、CRIを退職して再び自分の会社を立ち上げるしかない。しかし、それが成功するかは別問題で、実際その二人の新しい計算機は失敗している。

すでに半年も前のワークショップだが、スライドが公開されていたので、いまさらながら閲覧してみた。

2011 Workshop on Achitechtures I: Exasclae and Beyond ​http://www.orau.gov/archI2011/

一通りみてみたが、内輪ネタがひどい(Murphy)のもあるし、ほとんど会社の技術的な宣伝しか書いてない(Borkar)のもあるし、スライドだけでは意味不明(Sterling)なものもあるし、色々。

Geistのスライドは"Resilience"をキーワードとして、exascale時代には今よりすべての部分(ハードもソフトも)でエラー率が桁で悪くなるので、どうすればよいかという話。ハード的に対応するのはもう不可能な時代であるから、また、OSやランタイムでどうにかするのもほとんど無理であるから、そもそものアルゴリズムをエラー前提で作る必要があるのでは、という提案。それとともに、大規模な並列計算機のシミュレーションをするためのソフトを作っている。プログラムのバグは再現できるから修正できるわけで、そもそも色々な「エラー」は再現することが困難であるからその場での対応はできず、せいぜいそれを検知して間違える前に戻ってやり直すというのが安直な解決策(checkpointing)ではある。しかし、今より大規模になるとcheckpointingをすることは事実上不可能になるので、エラー前提アルゴリズムを検討するために、アーキテクチャや計算機のシミュレーションが必要だということ。

Resnick("Within and Above the Exascale Program"のほう)の8ページにはぐうの音も出ない正論が。

Co-Design Is Not:

• Having a meeting upfront, and then parallel working groups that then mash separate results together into a compromise implementation
• Having monthly meetings—and then mostly ignoring any changed direction (except for complaints)

『「Co-Design」ってそもそもなんなの？』というようなとぼけた疑問を今更言い出したり、議論しているようでは話にならないのは自明。9ページには"Co-Design Should"の説明がある。究極的には、十分資金的なバックアップのあるベンチャー企業のような組織を立ち上げて、一人のリーダーが責任を負い、複数のグループをまとめて全ての最終決断をする、みたいな構造にしないとまともな計算機を作ることはできない。これは、スーパーコンピューターの歴史を振り返ると結構明白である(例えばBrooks & Blaauwの"Computer Architecture: Concepts and Evolution" ( ​http://dl.acm.org/citation.cfm?id=548907 )を参照。またトレイシー・キダーの「超マシン誕生」を読む。いますぐ)。それが、日本で公的な資金のものとで実現できるのかというのは、大きな問題ではある。

Shalfのスライドは実際の"Co-Design"について、かなり詳細な議論をしている。ShalfはLBNLで推進されているGreen Flashという気象の問題のための計算機開発のリーダーであり、その言葉には経験に裏打ちされた重みがある。x86命令のなかで全く使われていない命令のテーブルとか、それだけでも面白い。以前Green Flashについて見たスライドでは、とにかくTensillicaのシステムでどうにかする、という印象しかなかったが、このスライドは変化が見られる。ポイントは、計算機のためのCHIPを作るの時に、一番コストがかかるのはアーキテクチャや設計自体ではなく、その色々なレベルでの検証であるということ。その解決策のひとつは、すでに検証済みのIPを組み合わせるということで、この場合でも組み合わせたシステムの検証は必要だから、そこはFPGAでエミュレーションして、設計と検証のサイクルをできるだけ短い時間で行うという方法論が示されている。そのほか、COTSは今は"commodity”ではなくて、現在は検証済みIPが"commodity"なんだというのは鋭い(42ページ)。COTSではpetascaleでもまともなシステムを構築はできないのというのは全く正しくて、それをできると言い張るのは幻想と妄想の組み合わせにすぎない。

Shalfのインタビュー記事も興味深い。​http://insidehpc.com/2010/08/12/rock-stars-of-hpc-john-shalf/

自身が計算機工学や計算機科学の勉強や研究をしてきて、その途中で天体物理学や相対論のための並列シミュレーションフレームワークの開発に関わってきた経緯を話した上で、

I remember the perspective I had when I was in each of those different roles
(when in EE, I thought the scientists were all just bad programmers, 
and when working for the apps group, I thought the hardware architects 
were just idiots who would not listen to the needs of the application developers)

異なる立場というのは「計算機」側と「アプリ」側。

エントリのタイトルの意味は、この分野で研究者をやるのなら、コンピューターのことは所詮1年も勉強すれば超最先端までいけるのだから後回しにして、まずはもっと時間がかかる科学や工学の考え方や現在の問題点を理解することを最初にやりましょう、という大学院教育への個人的な意見。Co-Designと無関係ですか？

Holmらによる一連の論文をリストする。完璧ではないかもしれない。最初の論文をもとに、宇宙論的シミュレーションのため、短距離成分の密度分布としてS2 shapeを採用したP3Mの定式化についてはこちら(このpdfは個人的なメモであり随時アップデートされる。無保証)。

How to mesh up Ewald sums. I. A theoretical and numerical comparison of various particle mesh routines

Deserno and Holm 1998

​http://dx.doi.org/10.1063/1.477414

Ewald和のP3M法による定式化をあらためてやり直し、P3Mやその類の方法(PME, SPME)の中では、P3Mが最もよいということを示した。根本的には、新しいことをやっているわけではないが、Hockney&Eastwood(1988)と比べると、非常にわかりやすい。

How to mesh up Ewald sums. II. An accurate error estimate for the particle–particle–particle-mesh algorithm

Deserno and Holm 1998

​http://dx.doi.org/10.1063/1.477415

上の論文の続きで、そこでのP3Mの誤差評価の式は、Green関数と、meshの数、boxの大きさなどを決めると、mesh点での総和として数値的に評価できるが、その演算量は無視できるほど小さくはない。そこで、P3Mの誤差評価の式を近似的に計算する方法を提案。近似手法でも十分正確な誤差評価ができるので、必要な精度が与えられた時に、P3Mの最適なパラメーター(mesh数や短距離力の範囲)を楽に決める方法論を示した。

The optimal P3M algorithm for computing electrostatic energies in periodic systems

Ballenegger, Cerdà, Lenz and Holm, 2008

​http://dx.doi.org/10.1063/1.2816570 airxiv:0708028

P3Mでの最適化されたグリーン関数は、加速度の誤差を最小化するように選ばれているが、この論文では、分子動力学などで必要なエネルギーを計算するための、エネルギーの誤差を最小化するようなグリーン関数の定式化している。実際上は、結果はそれほど大きく違うわけではないようではある。

How to convert SPME to P3M: influence functions and error estimates

Vincent Ballenegger, Juan Jose Cerdà, and Christian Holm, 2012

​http://dx.doi.org/10.1021/ct2001792

SPME(Smoothed Particle Mesh Ewald)をP3Mのフレームワークで捉え直す。まだ読みかけ。

総論

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

凡例

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キーを使っている。

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よりもよい。

OpenCLでSPH法のプログラムを実装し、様々なプロセッサで性能を比較している。

このSPHのプログラムは、​http://arxiv.org/abs/1112.4539 の論文で発表したGPUでの直接treewalkを拡張したもので、この論文のAppendixにある、OpenCL Cで書かれたコードがベースとなっている。我々の手法は、treewalkだけをGPUなりmulti core CPUで実行するというものであり、我々の意見ではこうするのが一番性能がでる。この論文は重力相互作用場合の性能評価をしているが、この方法の拡張で、SPHの近接相互作用の和を直接計算することができる。粒子をキャッシュにやさしい順番に並び替えるということも同じく効果的である。 さらに、現時点のコードでは、ホスト計算機でのtreeのconstructionを含めて、全体的に実装しなおして高速化した。これは、GPUによりtreewalkが高速になった反面、今度はtree constructionがボトルネックになってきたからである。上記論文のデータでは、80万粒子の場合、tree constructionには約30％の時間がかかっている。アムダールの法則を持ち出すまでもなく、より大きい粒子の計算では、GPUの高速さ(ベクトル計算機として)の意義が失われてしまう。この論文のコードでは、tree constructionには、シングルスレッドで八分木のポインタに基づいた方法を使っていた。これはBarnesのオリジナルのFortranでの実装とほぼ同様である。今のコードではMorton keyを並列ソートしてから、直接Next/Moreのポインタを生成する手法とした(Next,Moreの意味は論文参照)。

このNext/Moreを使う方法は、文献によると"Brother/child representation"と呼ばれるらしい(参照： ​http://onlinelibrary.wiley.com/doi/10.1111/j.1467-8659.2007.01104.x/abstract )。あるいは、CS業界ではfully threaded treeと呼ばれるらしい。treewalkは再帰関数として記述することができ、非常に美しいアルゴリズムの一つであるが、再帰関数はそのまま扱うとスタックの操作が頻繁になるため、ループに比べると性能がよくないかもしれない。末尾再帰はループに変換することができる。つまりNext/Moreのポインタを使うことで、treewalkを再帰を使わずに、また明示的なスタック処理も必要なく扱うことができる。これはGPUに向いたアルゴリズムで(最初に提案されたのはDubinskiの修士論文(1988) or Makino&Hut 1989 ？かもしれない)、Makino(1990)にて提案されたtree法のベクトル化手法の一つである。実はこの時、他に二人の著者(BarnesとHernquist)がそれぞれ別のベクトル化を提案して、この三つの論文はまとめてJournal of Computational Physicsに掲載されている。このMakinoによる末尾再帰のループへの変換の応用といい、そもそもBarnesによる八分木ツリーの提案も、彼らがIASに滞在していた頃におこなわれている（ような気がする）。しかも、これはどちらもConnection Machineに関係していて、そしてそれはMITのLISP/Schemeハッカー達の多大なる影響のもとで発明されたと勝手に思っているが、どうだろうか？機会があれば、本人に尋ねてみよう。リンク。と、ものすごい脱線はさておき。

この方法の利点は、treeのデータ構造として必ずしも八分木を想定しなくてもよいことである。そもそもオリジナルのtree法では、treeの各ノードの粒子数が1個になるまで、再帰的にtreeデータを構築していた。だが、粒子数が増えてきた時には、1個になるまでノードを追加していくのは時間がかかるし、実際問題必要のないことが90年(?)頃に発見されて、特にGRAPEとtree法を組み合わせる場合、通常はノード内の粒子が数十個を下回ったらもうノードを作らないほうが速いし、力の精度もよい。こうすることで、単純には元々のtree法の再帰的動的であるという美しさが失われてしまうが、実際上速いのに越したことはない。我々の方法ではソートされた粒子のリストから、基本的には八分木のデータを生成する。ただし八分木のポインタは使わない。また、下の階層のノードでは、ノード内の粒子数が指定した値を下回った場合には、それ以上の階層を作るのをやめて、代わりにNextポインタで粒子を直接リンクする。既にソートされているものをポインタを使ってリンクしなおすことで、実は効率が落ちてしまうが、こうするとtreewalkカーネルは修正をする必要がない。これは、実効的には、最下層ノードはN分木になっていて、その上は八分木になっているようなツリーになる。またポインタを使うことで、上の階層の八分木も保持する必然はなくなる。幾何学的に局所性が保たれるのであればN分木に変更してもよいということになる。これが意味を持つのは並列化をする時である。

さて、性能評価は以下の図。図は随時アップデートされるかもしれない。やっている計算は、いわゆる冷たいガス球の収縮で、hydrodyamicsとself gravityを解いている。状態方程式は理想気体で、もちろん断熱である。粒子数を横軸として1 stepあたりの平均計算時間をプロットした。計算時間はメインループの計算時間より算出しているので、当然tree constructionやGPUとのデータ転送等すべてを含む。演算は全て単精度。IntelのCPUではIntelのSDKを、AMDのGPUとCPUはAMDのSDKを、NVIDIAのGPUではNVIDIAのSDKを使っている。

以下コメント：

• OpenCLのSDKは使いものになる。この結果は、コードベースはひとつで各社のSDKで問題なく動作した。
• GPUではTahitiが一番速い。
• FermiのほうがGT200より速い。
• 24 core機がCPUの中では一番速い。
• Sandybridge-Eは24 coreのOpteronとあまり性能が変わらない。
• APUの結果はGPUを使った場合。Nの小さい所で大健闘している。次世代のAPUは面白いかも。

。。。

コメントはこちら ​https://plus.google.com/117333869988556483969/posts/Y144K5p53LM

AMDのGPU用のデバイスドライバ、Catalyst 12.1が公開された。4x4 DGEMMカーネルの性能がどれくらい変わったかをベンチしてみた。12.1とは「CAL 1.4.1664 (VM)」。

最大性能は11.12では605 GFLOPSだったのが、12.1では675 GFLOPSまで上昇した。一割以上向上したことになる。生成されたコードは微妙に違う。入力のILコードはどちらも同じ。

11.12 ​https://gist.github.com/1632529

12.1 ​https://gist.github.com/1700860

メインループの部分、12.1では一度に全部のデータをロードしているのに対して、11.12では2回に分けてロードしている。ILの記述は11.12に近い。12.1はロードをまとめるような最適化をするようだ。

いずれの場合にも、メインループでネックになっていのるのは、小行列の行や列をロードするアドレス生成の部分。ここが長くなる理由は、今の行列格納形式がリニアなためで、この場合、小行列のロードは行列のサイズ分のstrideアクセスになる。一回のロードでは倍精度変数2語しかロードできないため、4語ロードのためには、アドレスが2つ必要であり、strideアクセス分のアドレス生成が必要になる。今のコードは1回のループアンローリングをしているので、アドレスは4個必要である。これをブロック化された格納形式にすると、この部分がリニアメモリの読み出しになり、ロード自体が高速になるはずで、かつアドレス計算も単純なインクリメントですむ。というような最適化はいずれ。

参考までに、SGEMMの性能は以下のような感じ。8x8のブロッキング。こちらはピーク性能比でかなり低いのでまだ最適化の余地が大きい。色々とテスト中。

コメントはこちらへ：​https://plus.google.com/117333869988556483969/posts/BpWbwwzggXR

先日発売されたAMD Radeon 7970の命令セット(ISA)について、実際のコードをみながらコメントしてみよう。例としてILでpixel shaderとして記述されたDGEMMカーネルから生成されたコードをgistに貼りつけた。利用したSDKのバージョンは"OpenCL 1.1 AMD-APP (831.4)"。

​https://gist.github.com/1632529

前説

まず、前提として、このカーネルでは4x4のレジスタブロッキングをしている。そのため、ループの中ではそれぞれ4ワードを行列A, Bから読み込み、外積形式出部分行列の積を計算し、和を保持する4x4 = 16ワードの部分行列Cをアップデートする構造となっている。詳しくは、

http://galaxy.u-aizu.ac.jp/note/wiki/Fast_GEMM_Implementation_On_Cypress

にあるスライド参照。このカーネルはスライドの説明にあるところのTNカーネルに対応する。

また我々の最新のDGEMMカーネルのパフォーマンスについては、

Multi-level Optimization of Matrix Multiplication for GPU-equipped Systems, K.Matsumoto etal. ​http://dx.doi.org/10.1016/j.procs.2011.04.036

を参照のこと。この論文では巨大行列のためのDGEMMの性能最適化について報告している。

概観

コードをみてすぐに分かるのは"v_"と"s_"のprefixがついた命令があることである。これはprefixから予想されるように、ベクトル命令とスカラー命令なのだろう。これらの命令は、パッと見たところ完全に混在している。レジスタは"vn"の形式で32bitごとにアクセスできることがわかる。また"v[n:n+1]"の形式は二つの連続する32bitレジスタを表すのだろうから、これは64bitの1ワードに相当する。このカーネルの場合これは倍精度変数である。

コード右側の""というコメント(と思われる)以降にある数字の意味は不明。普通に考えるとbinary形式だろうけど。

個々の命令については現時点ではdocumentが公開されていないので、名前から動作を推測するしかない。とはいえ、一部を除いて、どの命令も明示的な名前が付けられており、特記するまでもない。"image_sample"命令は、texture unit経由のデータ読み込みを意味すると思われる。

レビュー

コードの流れは単純で、最初に部分行列Cを初期化して(20 - 51行:32bitで32ワードを0で初期化)、52行のラベルからループに突入し、123行の"s_branch"命令がループの終了部分。それ以降は"alpha C + beta AB"を計算している。この最後の部分(125 - 252行)が妙に長くなっているが、カーネルあたり一度しか実行されないので、計算時間上のインパクトはない。実際、上記K.Matsumoto論文によると、GEMMのカーネルとしては、このカーネルのように"C <- alpha C + beta AB"の全てを計算するのは高速ではなく、ABを計算するカーネルを利用し、Cとの和はホストで計算したほうが高速である。ここではあくまで一番シンプルな実装を紹介しているにすぎない。いずれにしろ、この部分では部分行列Cをロードして、レジスタに保持されているカーネルで計算した結果との和をとる処理をするため、"tbuffer_load"や"tbuffer_store"という命令が使われている。

ということでDGEMMカーネルとして計算時間がかかるループの本体部分は52 - 123行である。このループの大部分は"v_fma_f64"という命令が占めている。これはその名の通りの倍精度FMA命令だろう。FMA命令はは素直に4 operandsの命令になっている(古い世代のCypressやCaymanでもFMA命令は4 operandsだった)。なお、AMD社は最新のx86_64のCPUでも4 operandsのFMA命令を採用している。一方Intel社は、将来的には3 operandsのFMA命令を採用することになっている。両社の方針の違いは、将来的にはx86_64 CPUに統合されるであろうGPUのアーキテクチャにも反映されている/深く関係しているのかもしれない。

53 - 60行はループの終了条件を判定している部分になる。"s_waitcnt"はよくわからない。こここでthreadが切り替わっているのかもしれない。

その意味では61行目のラベルからがループ本体であり、62 - 67行はデータをロードするアドレスの準備をしているようにみえる。この部分は冗長にみえる。

69 - 72行で行列Aと行列Bをロードしている。"image_sample"では、一番左のレジスタがdestinationだろうから、ひとつの命令あたり32bitで4ワードをロードしている。よって、この4命令で合計倍精度変数を8ワード、ロードすることになる。

74 - 90行で、ロードした変数どうしの外積形式による行列乗算がおこなわれる。"v_fma_f64"の個数は16個であり、これは4x4のレジスタブロッキングをしているので当然。ところどころwait用？の命令が挟み込まれている。

91 - 100行でまたなんらかのアドレスの準備をして、101 - 104行で再び行列AとBをロードする。実はこのカーネルでは2段のループアンローリングをしているのでこうなっている。実際74 - 90行と処理自体は同等だ。

不可解なのはループを抜けたあとの126 - 147行で、なぜかレジスタ同士で入れ替えをしている。入れ替えをする理由は、その後で結果をメモリに書き込むためで、その時にレジスタ上で連続するようにという理由付けはできるが、実際にはv15 - v36に連続しているものをv12 - v33にコピーしているので、よくわからない。

追加コメント

全体を見ると、これくらいであれば手で書こうと思えば書けなくもないかもしれない。ただ、アドレス計算はやっぱり面倒。

ループ本体の中にデータロードとFMA以外の命令が43%も含まれている。このコードの性能は最大600 GFLOPSであるが、これはピーク性能の63%に相当する。これは偶然なのかどうか、というところが、今後の最適化の指針になるだろう。

126 - 147行は明らかに不必要。

随時アップデートするかもー。

​http://gitlabhq.com/ のインストールは罠が多い。失敗を何度か繰り返した後に、結局この ​https://github.com/gitlabhq/gitlabhq/wiki/V2.0-easy-setup-for-ubuntu とおりに一字一句やればよい。Ubuntu 10.04 LTSで確認。この文書には書いていないがzlib1g-devをインストールしておかないと、ruby1.9がzlibなしになってしまい、うまく動かない。ちなみに"git"と"gitlabhq"というuserを追加する。他のドキュメントはここを省略しているので、"git"というuserのみでやろうとするとハマる。

次にwebrickで動作しているこのGitLabアプリをapacheを通して外から見えるようにしたい。ここでまた色々な罠が。最終的には、現バージョンのRails3(3.1.3)でURLにprefixをつけるには

​http://stackoverflow.com/a/7367886

が正解。config/routes.rbに

Gitlab::Application.routes.draw do
scope "/git" do
...略...
end
end

"scope"と"end"の二行を、内側に追加する。その上で、apacheのreverse proxyの設定は

ProxyPass        /git http://localhost:3000/git
ProxyPassReverse /git http://localhost:3000/git
ProxyPass        /assets http://localhost:3000/assets
ProxyPassReverse /assets http://localhost:3000/assets

とする必要がある(example.com/gitでアクセスする場合)。後者の設定は/asstes配下の画像ファイルやJSファイルのために必要。routes.rbを上記のように書き換えても、/assets配下のURLは。。。うーん。

magic fileがないみたいなエラーがでる場合

sudo gem install charlock_holmes

で、明示的に"charlock_holmes"すればよい。

2011年12月22日にRadeon HD 7900(Tahiti)について発表があった。アーキテクチャの概要については、日本語ではPC Watchの以下の記事にまとまっている。

​http://pc.watch.impress.co.jp/docs/column/kaigai/20111222_501138.html

GPUのアーキテクチャ名はSouthern Islandsとなる。ちなみに前世代(Radeon HD 6900: Cayman)はNorthern Islandsと呼ばれている。実際、AMDが12月に公開した文書​AMD Accelerated Parallel Processing OpenCL™ Programming Guide (v1.3f) の4-79ページに、新旧アーキテクチャの比較がまとめられている。

プログラミングの観点からすると、大きな違いはVLIW4と呼ばれるベクトル長4の演算ユニットの塊から、スカラー演算ユニット(言い換えるとベクトル長1)の塊へなったこと。そのためOpenCL用語でCompute Unit(CU)と呼ばれるCPUにおける「コア」に相当するユニットのベクトル長は、北島では4x16=64であったのが、南島では16になった。以下に、上記の文書にある図をコピーして貼り付ける。

北島の図はVLIW4のユニットであり、北島のCUはこれが16個まとまったものになる(そのためこの図は正確ではないことに注意)。南島の図はCU全体である。どこかで見かけたAMDの人のプレゼンテーション(検索中)によると、演算ユニット自体はVLIW4を流用しているとのこと。そのため、南島はベクトル長16のはずなので、図の「VECTOR ALU」はベクトル長4のはずであり、この部分は北島の図の「xyzw」演算器の部分を流用しているようだ。これらの図における演算器の最小単位は32 bit演算器である。VLIW4を流用していることから、浮動小数点単精度演算の性能と倍精度演算の性能比は、新旧で変わらず加算では2:1で、乗算またはFMAでは4:1となる。実際には、これに詳細不明のスカラー演算器の性能分が加わるが、ベクトル演算器とスカラー演算器の個数比は16:1のため、ピーク性能へのインパクトは大きくない。

もう一つ、この文書にはメモリ要素を比べたテーブルがある。

32 bit演算器あたりのレジスタ数は新旧で同じ256ワードであるが、北島はベクトル長4のベクトルレジスタであるのに対し、南島はレジスタもスカラーであるから、完璧にVLIWでベクトル化されているコードを除いて、レジスタの利用効率は南島のほうがよいはずである。それと、CUあたりの1次キャッシュが倍増し、CUあたりのLDS量も倍増して、さらにメモリチャンネルが4から6になったことで、2次キャッシュも増えている。キャッシュは北島では読み出し専用だったのが、読み書き可能になったとのこと。

Tahitiのリファレンス動作周波数は925MHzであり、この時の演算性能は単精度FMAで3788 GFLOPSになる。一方で、海外サイトのソース不明記事によると、オーバークロックされたGPUボードがすでに準備されているらしく、その動作周波数は1.3GHzに達するようだ。この場合の性能は約5.3/1.3 TFLOPS(単精度/倍精度FMA)になる。TSMCの28nmプロセスを待ったかいは非常にあるのかもしれない。

2010年11月に発表されたtop500からGPUクラスターによる結果を抜粋。ハイライトによるとtop500にはNVIDIAのGPUを使ったクラスターが10システム、AMDのGPUを使ったクラスターが1システムあるとのこと。であるが、全部を見つけることはできなかったので、以下ではそのうちの上位のシステムのみを抽出した。Fermiを使ったシステムが4システム、Radeon(Cypress)を使ったシステムが1システムである。

CPUコア数とGPUコア数はいずれもノードあたり。総コア数はノード数×(CPUコア数＋GPUコア数)と等しいはずなので、ノード数が公表されていない場合は、総コア数を信じてノード数を計算した。メモリは公表されているNmaxの値から逆算した値で、システムの総メモリ量ではない。22位のシステムは、top500で公表されている情報がSC10の会場で配られていた資料と異なるので、配られていた資料の数字を載せた。ただしノード数はRpeakから逆算した数(少し変？)。

ノードあたりのCPU,GPUの性能とGPUの性能の割合、そしてノードあたりのメモリ量は以下のとおりになる。

おおまかに言えば、4位145位のシステムはGPUヘビーであるのに対して22位のシステムはCPUヘビーであり、1位3位のシステムはその中間程度になっている。定性的には、Linpackベンチマークは1ノードあたりのCPUの演算性能の割合が多いほど、そしてノードあたりのメモリが多いほど、Rmax/Rpeak値はよくなる。その意味では145位のシステムは相当バランスを欠いており、メモリが少ないために、ノード数が少ないにもかかわらずRmax/Rpeak値が低い。3位のシステムも同様である。数字上ほどよくバランスが取れているのは1位と22位のシステムである。4位のシステムはメモリがもっと多く使えるのならより高効率になっただろう。ただし、Linpackベンチマークの観点でバランスを欠いているからといって、そのシステムが使いものにならないということはない。GPUクラスターにはある種の特化した目的があるはずなので、その目的が達成できるように設計されているのなら、ノードあたりのメモリが少なくとも問題はない。

FermiのDGEMM性能はそもそも1 GPUあたりのRmax/Rpeak値が60%しかないため、システムがGPUヘビーであるならば、全体のRmax/Rpeak値が60%を超えることは困難である。それに対してCypressでは最大では87%の効率であり(我々の実績値)、22位のシステムはCPUヘビーでもあるので、Rmax/Rpeak値が高くなるのは当然であろう。とはいえFermiのシステムと比べて、Linpackでの効率がよいのは明らかであり、Cypressアーキテクチャの優位性を示す。

22位のシステムについては、会場で配られていた資料によると、1ノードでのLinpack性能は563.2 GFLOPS(75.5%)を達成していて、現時点では1 GPUのシステムで最速であろう。ただしLU分解ではCPUも活用するため、異なるCPUやコア数のシステムを単純に比較することはできないことに注意。CypressでのDGEMMルーチンは我々のものではなく、彼らが独自で実装したものであると、議論をしたヨハン・ヴォルフガング・ゲーテ大学フランクフルト・アム・マインの研究者ら(なんとPh.D students！)は言っていた。

単位電力あたりの性能のランキングであるgreen500の結果もほぼ同時に公表されているが、消費電力測定のプロなどというものがいるかどうかは別としても、素人考えでも異なるサイトにおける全く異なるシステムの消費電力の測定を、意味のある形で科学的に比較することは単純ではない。よって、現時点でgrenn500の結果に意味があるかは大きな疑問符がつく。できるだけ科学的に比較できるような仕組みを真剣に議論する必要がある。

Fox, Williams, Messinaによる"Parallel Computing Works"(ISBN 1-55860-253-4 Morgan Kaufmann Publishers, Inc. 1994)の20章の訳を仮公開。

​https://docs.google.com/document/edit?id=1EdAV8_bwZwjiWZi46PQztlP3H6u2rmdND7r41nuvyI0&hl=ja

20章だけでなく他の章も面白いので、大著で通読するのは大変だし古い内容もあるが、最初の数章は読むべき価値はあると思う。一年前に書いたCosmic Cubeも、実はC3Pの前にCaltechで実現されたものである。そこからはじまったであろうQCD用に計算機を「手作り」するという考え方は、現在でもQPACEに生きている。

以下、最近知った蛇足情報。天文のN体シミュレーションをやっている人々の中では、並列ツリーコードで有名なJohn Salmonは、この本の最初の著者のGeoffrey C. FoxのところでPhDを取っていて、この本にもSalmonによる並列ツリーコードの記事がある。で、Salmonは現在、Antonを研究開発しているD.E.Shaw Researchで働いている。

​http://www.deshawresearch.com/members_compsci_salmon.html

まさに20章で述べられている意味でのparallel computing works！の実例の一人と言えるんではないだろうか。

『スーパーコンピュータ』シドニー・ファーンバック(長島重男訳)を某新刊書店で発見したので捕獲する。消費税が３％時代の値札修正シールがついた1988年4月1日初版本である。なぜこんな古い本が未だに陳列されているのかよくわからないが、内容は80年代前半のスーパーコンピュータのアーキテクチャについて。一番ページが費やされているのは、各社のベクトル計算機についてである。登場順に並べる。

最後の二つ以外はすべてベクトル計算機である。

各社の計算機紹介についての章の前に、ベクトル処理に関する章があり、それを読むことで誤解が解けた。それは、ベクトル処理を実現するための手段として、アレイ計算機とパイプラインベクトル計算機の二つの手段があり、通常のベクトル計算機(上記システム全て)はパイプラインベクトル計算機である。実際には、上記のシステムの中にはパイプラインを複数持つものもあるが、たかだか4本程度である。実際、SX-2は4本のパイプラインがベクトルを4入力ごとに処理する機構を持つ。アレイ計算機とは、概念図を見ると現代のGPUと似た構成で、相互結合ネットワークを介して主記憶と演算器が接続されていて、複数の演算器でベクトルデータを並列に処理する計算機である。例として挙げられている計算機はILLIAC IV, ICL DAP, Burroughs BSP, Goodyear MPPである。MPP以外にはこの本での紹介はない。80年代前半の時点で、普通のベクトル計算機のほうが有利(作りやすい？)であることは明らかだったようだ。

「GPUはベクトル計算機」と言った時に、個人的に利用したことのあるベクトル計算機であるVPPシリーズやSXシリーズを思い浮かべると、どうにも整合性がとれない。この本で他のベクトル計算機も含めて比べてみると、GPUと構造上で似ている点は全くないと結論してよいことがわかった。ちなみに、Cray, CDC, 日立, 富士通, NECの各システムを比較すると、どれも似たような技術を使った似たような構造の計算機であり、各システムに明確な差はないようであった。論理は100 - 1000ゲートの素子を組み合わせて動作速度は10nsくらい、主記憶はアクセス時間は10 - 60 nsのSRAMで大きさは100MBクラス、外部記憶にはDRAMが利用されていて大きさは1000MBクラスである。

GPUがベクトル処理を実現する手段としてのアレイ計算機であると考えれば、「GPUはベクトル計算機」という言い方でもしっくりくる。とはいえ、アレイ計算機として唯一紹介されているMPPの章を読むと、これは1 bitプロセッサを格子状に接続した計算機でほとんどFPGAのようなものであり、動作がアレイ計算機的になるのは理解できるとしても、その構造を考えるとあまりしっくりこない。ついでながら、日立のS-810は名前に反して他のベクトル計算機と同様の構成であった。そういう意味で、GPU的なベクトル処理機構というのは、80年代の前半の時点で既に死に絶えていたと言えるのかもしれない。

上記の中ではCYBERPLUSだけが他と比べると少し異質で、これは現代の分散並列計算機に近い。演算プロセッサはリングネットワークで接続されていて、データ転送とプロセッサでの演算は独立に動作するようだ。ただしそのプロセッサは、15個の演算ユニット(整数加算、整数乗算、シフト論理演算、メモリ、レジスタ、制御等等)をクロスバで接続したものであり、その構造はかなり奇異に感じる。その制御のための命令は合計で200 bit幅にもなるとのこと。さらに浮動小数点拡張をいれると240 bitとなるらしい。とはいえ、単に普通はユーザーからは隠されているマイクロコードが露出しているだけで、現代のプロセッサでも、水平なマイクロコードで制御するとこの位の情報量になるのかもしれないが。

IEEE Cluster 2009でのJohn Gustafson(の法則のGustafson)による"Future Clusters"と題したキーノートより。

世界最初のclusterは"Cosmic Cube"で、1984年にChuck Seitz と Geoffrey Fox が作ったQCD用の計算機。予算は5万ドル。スライドを撮った写真の肝心な部分が写ってない！ので、不明だが、8086プロセッサに8087のコプロセッサを使ったノードが、64ノードからなる。写真をみると、名前のとおり、個々のノードは現在のcube型PCとほぼ同じような姿形をしている。そのため、全体が非常にコンパクトにまとまっている。

​http://www.netlib.org/utk/lsi/pcwLSI/text/node13.html

もうひとつ、CPUでの「演算」の様々な部分における消費電力のおおまかな値。スライドを撮った写真より再構築。

iPhoneのカメラでは、前の方に座っても限界があった。

T2K用の最適化説明であるが、普通のLinuxにも適用可能。RDTSCの使い方もあり。

​http://www.cc.u-tokyo.ac.jp/publication/news/VOL10/No5/200809tuning.pdf

自律的に動作するMD専用計算機。

GRAPE的な力のパイプラインを内蔵し、長距離力用の演算回路も別口で持っている。GRAPEのように全体全の完全なN²演算ではなく、対称性を考慮した力の計算をし、さらに距離計算を低精度でおこなうことで、必要な近接相互作用のみを計算する仕組みがはいっている。そのため超速い。SC2009のBest Paper Awardにノミネートされた。

• ​http://en.wikipedia.org/wiki/Anton_(computer)
• ​SC2009のTechnical Paper
• ​Millisecond-Scale Molecular Dynamics Simulation on Anton

SC2009では、Anton関連でもうひとつ、FFTを高速に実装した論文もacceptされている。

wikipediaからリンクされていた2本の論文を読んだ。

• ポイントは力の計算だけを速くできたとしても、アムダールの法則により他の部分が足を引っ張るため、高速化には限界がある。そこで、他の部分もASICで実行できるようにすることで、Nが小さくかつ長時間積分に特化した計算機を作るということ。
• ASIC 90nm (dieの大きさは不明)
• ASICの基本動作速度は400 MHzで、HTISのみ800 MHzの倍速動作らしい。
• HTIS(high-throughput interaction subsystem)というforce pipelineとflexible subsystemという部分からなる
• HTISは二体力とmeshへの密度あり当てと、meshからの力の補完を計算する
  • fixed point integer for a position
  • 8bit for distant used in NT method (threahold距離より遠い粒子をフィルタリングする)
  • 26bit for force pipeline
  • 32 force pipelines (PPIP) on a chip (cf. MD-GRAPE3 20 pipelines @ 350 MHz)
• flexible subsystemが残りの計算をすべて行う
  • FFTによる長距離力
  • 数値積分
  • bond force
  • 拘束条件の適用
  • カスタマイズされたTensilicaのcoreに4-vector SIMD演算器を2個(こちらは独自設計らしい)
  • correction pipelineはPPIPと等価
• 使い方は詳しく書いてないが、少なくともTensilicaのcoreとSIMD演算器のプログラムが必要だろう。
• 比較に使われているMD-GRAPE3の結果が遅すぎるような気もする。

• Antonと比べるとMD-GRAPE3のアプローチは、システム全体の総演算速度では圧倒的に大きいが、Nが小さい系(<10⁵くらい？よくわからない)をミリ秒くらい計算したいという用途には、有効ではないということになる。この様な計算は、1 stepあたり10¹⁰演算を10⁹ステップ計算する必要があるので、総演算量は10¹⁹演算となる。これはNが単純に大きい大規模計算より困難があって、それはたった10¹⁰演算ごとにglobalな同期が必要であるので、1 PFLOS(10¹⁵)の計算速度をもつような超並列計算機では効率があがらないため。

OMakeを使うとc-x c-s c-z make [enter]のループから開放される！！！！

ソース　​http://omake.metaprl.org/index.html

ドキュメント　​http://omake.metaprl.org/manual/omake.html

日本語のまとめ ​http://unicus.ddo.jp/omake-wiki/index.php

だけでなくlatexのコンパイルからdviファイルそしてPDFの生成まで自動的にしてくれるのは、本当に便利なので、これだけでも使う価値がある。build systemとしては、まだテスト中。

​http://d.hatena.ne.jp/hayamiz/20081208/1228727272

top levelのOMakefileで、変数を定義するときには以下のように、".SUBDIRES"の前で定義しないと、sub directoryのbuild時に反映されない。top levelのOMakefileの例:

.PHONY: clean

LIBDIR = /usr/local/lib/foobar
.SUBDIRS: subdir1 subdir2

PowerXCell 8iをカスタムネットワーク(FPGAで実装)で接続したQCD用の計算機。

​http://en.wikipedia.org/wiki/QPACE

​http://www.itwm.fhg.de/hpc/workshop/mic/Qpace_(Dirk_Pleiter_-_Desy).pdf

概要は後者のプレゼンファイルがわかりやすい。PowerXCell 8iとVirtex5 LX110Tを同一ボード上に実装して、3次元のトーラスネットワークを実現。Virtex5から6本のlocal通信ネットワークと、1本のethernetを引き出す。localネットワークは10GbE PHYのチップ(PMC Ceirra PM8348, XAUI, 1 GB/sec(10B/8Bのあと))を6個外付けしている。

2009年2月のミーティングでのPDF ​http://www.fz-juelich.de/jsc/juice/eQPACE_Meeting_PDFs

今日配布されたuniv2000(​http://www.univ2000.com/ )のちらしに、 FireStream9270/9250の広告があった(これらはmodel numberから4400を引いてHDをつけたRadeonと同じものだろう)。ACUBEという会社が元の代理店。​http://www.acube-corp.com/products/firestream/

値段は「アカデミックプライス」だと超高いというほどでもなく、4400引いたものの秋葉原価格x5くらい。結果的に、値段はTeslaのアカデミック割引価格とほとんど同じになっている。パフォーマンスあたり価格はこちらのほうが上。上のURLには2年保証とあるが、アカデミックバージョンが2年保証かどうかはちらしには書いてなかった。

HDDを16台マウントするための専用金具 ​http://hddbancho.co.jp/