SIMD 幅は広く、レジスタ数が少なく、演算器レイテンシが大き
く、L1 はあまり大きくなく、バンド幅もL2はそこそこあるがレイテンシが
大きくてバンド幅も小さいマシンを考えてみます。もちろん、具体的な
対象は A64fx です。「京」も同様な問題があったのですが、「陽解法規則格
子流体でピーク性能の 10ー15%しかでない」ということについてです。
B/F=4 ある古典ベクトルなら50-60% はでるわけで、B/F が 1/8 で 性能が
1/4 なら偉い、ともいえますが、一方、陽解法流体計算は、「理論的には」
必要 B/F は極めて小さいです。というのは、格子点あたり5変数に
対して、最近の高精度スキームなら1ステップで3000演算程度必要だからで
す。
つまり、理想的にキャッシュを使えるならば、80バイトの read/write に対し
て3000演算ですから、必要な B/F は 0.03 程度です。実際、2016年の
ゴードンベル賞を獲得した Sunway Taihulight 上の陰的スキームによるGCM
のグループは、陽解法スキームの実装の性能もだしており、20%を実現しています。
それに比べて10倍以上の B/F がある「京」や富岳で、
もっと高い効率を実現できていないのは何故か？という問題です。
Sunway で高い効率を実現できた大きな理由は、キャッシュの代わりに
主記憶とDMAで通信するローカルメモリを使わざるを得ない、というハードウェ
アアーキテクチャにあります。まともに性能がでるプログラムにするためには、
必要なデータをローカルメモリにもってきて、そこでできる計算は全部して、
それから結果を書き戻す、というスキームにせざるをえません。また、
通信はキャッシュによる非明示的なものではなくDMAによる明示的なものなの
で、キャッシュラインを踏むと無駄なところまできてしまうとか、そのために
必要なものがパージされるとかいったことがなく、細い主記憶バンド幅を最大
限有効に使うことができる、というか使わざるを得ないわけです。
ローカルメモリをもつアーキテクチャのもうひとつの利点は、ローカルメモリの高い B/F にあ
ります。演算ユニットは1つで、クロック毎に load/store のどちらかが演算
と同時に発行できるのが普通ですから、 B/F は4あります。この高い B/F で
ローカルメモリ全体をアクセスできるわけです。 A64fx は公開情報によると
L2 キャッシュの B/F は「1.2」 であり、Sunway Taihuligtht の 1/3 しかあ
りません。 L1 も、2load or 1 store なので平均をとると B/F が3、最悪ケー
スでは 2 になります。アーキテクチャの性質から普通のアーキテクチャより
ストア命令が多くなりがちな A64fx でストアのスループットが小さいのは残
念なことです。
とはいえ、おそらくキャッシュかどうかよりもこのローカルメモリバンド幅が
性能には重要であり、そのことは同様に L1キャッシュやローカルメモリのB/F
は4を実現している(LLC は2)PEZY-SC2 でやは陽解法流体コードで20% 前後の
効率を実現できていることが示していると考えられます。もちろん、SC2 のこ
の効率は、1段で空間4次時間3次を実現する新しい差分スキーム SL4TH3 と、
さらにtemporal blocking まで導入することで実現したもので、Sunway に比
べて主記憶のバンド幅の利用効率は悪くなっています。
SL4TH3 は collocation grid の5点差分での差分法で、有限体積法ではなく
リーマンソルバも使わないもので、衝撃波捕獲よりも、大規模計算での非圧縮
スキームの代替を狙うものです。5点差分では足りないなら7点にしましょうか
とかそういう話です。flux splitting をしないで、高次の cross derivative
も全部計算し、それから高次の時間導関数を直接計算することで、ルンゲクッ
タのような多段法を使わないで高次精度を実現します。
cross derivative は、例えば u_xy であれば、各点で u_x を計算して、
それをy方向5点で差分とるのと、逆の順番で計算するので、結果は同じになる
ので、座標について任意の順序で計算できます。
なお、時間については2階導関数までしか計算してなくて、予測子・修正子公
式としてエルミート公式を使うことで時間3次にしています。
以下、同じスキームを A64fx に実装することを考えてみます。
temporal blocking はしないとして、8MB の L2キャッシュを可能な限り
有効に使うためには、例えば x, y の2次面内でブロッキングして、
z 方向にスキャンしていくスキームが考えられます。空間微分の計算量を
最小化するためには、z方向の5面について、xy平面内の4次までの導関数
15個をもつことが望ましいです。但し、例えば x 方向だけもっておいてy方向
は再計算しても、演算量は極端には増えません。これは、空間微分から時間微
分に変換するところの計算量も十分大きいためです。以下では、15個の空間微
分を倍精度でもつ場合についてまず検討します。
3次元流体は5変数あるので、格子点あたり 75個の、さらにz方向の幅5を保持
すると375個の変数となり、ほぼ3KBです。従って、2500点、富岳が
12コアグループであることを考慮すると例えば 32x48 または 16x96 で
ブロッキングするのが妥当でしょう。y方向の幅が4または8、x方向は
SIMD 化する、ということです。
演算数を3000程度とすると、1面の演算数は 4.6e6 しかなく、理論ピークの
705Gflops では6マイクロ秒で終わります。いかに富岳のハードウェアバリア
機能が素晴らしいとはいえ、6マイクロ秒の中で複数回同期をとるようでは大
きな性能低下が避け難いので、複数面、あるいは z 方向のスキャン自体を
同期なしでできる必要があります。これは、x方向の差分計算を、コアの担当
部分からステンシルの幅だけはみだしてやって重複して結果をを保持しておけばよいこ
とになります。この分による必要メモリと計算量の増加は 16x96 のほうが
32x48 よりも理論的には小さくなります。
差分計算は、5語読んで、25演算しますが、そのうちFMA になるのが6演算しか
ないので最低19サイクル必要であり、理論限界性能は 66% になります。
空間導関数から時間導関数への変換は、基本的にはそれぞれの変数を掛けて足
すとかなので割合FMAが多いですが、逆に load limit になるので理論限界は
50% 程度かもしれません。また、問題になりえるのは L2 のリードレイテンシ
で、プリフェッチが必要かもしれません。
とはいえ、そういうわけで1段法ならば
20%よりもうちょっとよい効率を実現できる可能性が理論的にはある、と
いうことでしょう。
多段法でも、1段の中で上と同様なブロッキングを行なうことで、原理的には
現状より高い性能の実現が可能「かもしれません」
とはいえ、こんなコードを人間が書くのは正気の沙汰ではなく、コード生成系
が必須です。