2021-01-03

CVPシミュレータを解読したメモ

新年になりました。今年もよろしくお願いします。

Value prediction（値予測）の世界的な大会、Second Championship Value Prediction @ HPCA'21というのが開催されているそうです。

値予測というのは投機的実行の一種で、レジスタに書き込まれる値を予測するものです。投機的実行の代表的な例として分岐予測がありますが、これは分岐命令のとび先を予測するものです。

どちらも、あくまで予測による投機的実行なので、実際に実行して予測が正しかったことを確認する必要があります。予測が間違っていた場合、その予測結果を（直接的にでも間接的にでも）使った命令を全てやり直す必要があります。特に、値予測の場合は分岐予測と違い、「なんでもいいから適当に予測しておく」のではなく「予測しない」という消極的選択肢を選んだ方が性能を高めることが多いです*1。

そういうわけで、値予測器の評価では（正答率やカバー率といった単純な指標ではなく、）実際にプログラムが高速に動作するかで評価を行う必要があります。そのために作られたシミュレータがCVP simulator（https://github.com/eric-rotenberg/CVP）です。

CVPシミュレータの特徴

CVPシミュレータの特徴は、以下の通りです。

Qualcomm社が提供するARM64のトレースを読み込んでシミュレーションする
- ARM64の命令は最大5つ*2のマイクロ命令に分解されてシミュレーションされる
（サイクル単位ではなく）命令単位でどのサイクルに何が起こるかをシミュレーションする
L1Iキャッシュをミスしない限り命令フェッチに時間がかからないとしている
分岐予測器にTAGE-SC-LとITTAGEを採用している
- BTBやRASは実装されていません

CVPシミュレータは、その構造上、以下の本質的な問題が存在します。

シミュレーションが正確ではない

サイクル順ではなくて命令順にシミュレーションしているため、（少なくとも）キャッシュの状態が正しくシミュレーションされません。具体的には、以下のようなパターンが問題になります。

あるロード命令由来のアクセスが、L1Dキャッシュミスを起こした後L2キャッシュでヒットする。
その次の命令が、L1Iキャッシュミスを起こした後L2キャッシュもミス。L3キャッシュのデータをL2キャッシュに挿入する。

このような場合、「L2キャッシュでヒットする」は「L2キャッシュに挿入する」の後に起こるかもしれず、その場合にはL2キャッシュの状態が変化しているので本当にヒットとは限りません。

命令キャッシュミスはそんなに発生しないので問題ないということなのでしょうか。しかし、配布されているトレースのうち一部のトレースでは23%も命令キャッシュミスが起きるらしい*3ので、「命令キャッシュミスはそんなに発生しない」という仮定は成り立っていないと思います。

とはいえ、サイクル順のシミュレーションをする場合はたくさんのコードを書かないといけないので、シミュレーションが多少不正確になっても問題ないとすれば、適当な妥協点なのかもしれません。

CVPシミュレータのおかしな点

CVPシミュレータの微妙な点（おそらく簡単に修正可能）は以下の通りです。

マイクロ命令への分解がナイーブ

ベースアドレス書き替え付きのロード命令がマイクロ命令に分解される場合、「ベースアドレス書き替えを実現する整数命令」と「ロード命令」に分けるべきですが、全てがロード命令に分解されます。意図していないアドレスの読み出しが行われるという問題もあります（aperaisさんも指摘しています→Aperais/cvp2v2 by aperais · Pull Request #3 · eric-rotenberg/CVP · GitHub）。

ベースアドレス更新のレイテンシに支配されているようにシミュレーションされているトレースでストライド予測を行ったらスコアが上がるのかもしれません。実際にはシミュレータの都合で不当に遅くなっているのを高速化しただけで、意味はないのですが……。

その他、フラグレジスタを介した依存関係を正しくシミュレーションできていません（トレースに記録されていないデータが必要なため、シミュレータで推測を行っているようですが、それが使われていません）。

ロード命令のアドレス変換には1cycleかかるのにストア命令のアドレス変換は0cycleで終わる

ロード命令の方はアドレス変換にかかる時間をシミュレーションしていますが、ストア命令の方はやっていません。

典型的にはストアアドレスの計算の方がデータ到着より早く終わるため、アドレス変換にかかる時間をカウントしていないのかもしれません。単に忘れていただけかもしれません。

ストアフォワーディングが0cycleで終わる

ストアを行ったサイクルにストアフォワーディングを待っているロード命令も終わります。

レジスタファイルの書き込みポート数をシミュレーションしていない

一サイクルに実行を開始できる命令の数（演算器の数）は正しくシミュレーションしていますが、一サイクルに実行を終了できる命令の数（レジスタファイルの書き込みポート数）はシミュレーションされていません。

固定レイテンシの命令だけであれば問題なさそうですが、キャッシュミスしたロード命令やストアフォワーディング待ちのロード命令は可変レイテンシであり、一サイクルに多数の命令が同時に終了します。

これは他のものと比べて実装難易度がやや高いですが、リプレイすべきサイクルを計算してリソース占有の予約を入れればシミュレーションできるはずです。

ここまで変な点をいくつか書きましたが、ChampSim（特にdevelopブランチは完全に壊れてしまっています）に比べたら全然ましです。

*1:実際、 https://www.microarch.org/cvp1/papers/Seznec_Unlimited_2020.pdfによれば、チューニングされた値予測器の正答率は99%を超えます。これは、プログラムの挙動が予測しやすいということではありません。値予測器には「予測しない」という選択肢が存在し、それを選ぶことが性能を高めるため、そうチューニングされているからです。

*2:SIMD レジスタを二つ読み込むload pair命令で、ベースレジスタ書き換え付きの場合に発生します。

*3:The 1st instruction prefetching championshipで使ったトレースはCVPで配布されているトレースを流用しているらしく、それを用いた評価（https://research.ece.ncsu.edu/ipc/wp-content/uploads/2020/05/eip.pdf）ではプリフェッチなしの場合23%も命令キャッシュミスが起こるとしています。

2020-12-20

CUDA libraryのexp(x)を読んだ

なんだか指数関数の実装ばかり見ているような気がしますが、GPGPU向けの実装はどうなっているんだろうなと思ったので見てみました。

Programming Guide :: CUDA Toolkit Documentationによれば、最大誤差1ULPを保証する実装になっているようです。 ULPの定義について詳しく語った論文を引用しているあたり、かなり丁寧なつくりになっていそうです。

中身を見てみたところ、11次の多項式近似をホーナー法で求め指数部をビット演算で操作する、というよくある感じの実装でした。

この方法でうまくいかない入力だった場合の例外処理が単精度浮動小数点数を使って実装されており面白かったです。これは、倍精度浮動小数点数演算器をなるべく使わないような実装をしているためだと考えられます。倍精度指数関数を計算する場合、倍精度浮動小数点数演算が大量に出てきますが、GPU上の倍精度浮動小数点数演算器は貴重なので、例外処理なんかに使いたくないということでしょう。

指数関数の入力は、大きく分けて以下の5つに分類されます*1。

小さすぎてexp(x)が0.0になるような入力（-Infを含む）
exp(x)が非正規化数になるような入力
exp(x)が正規化数になるような入力
大きすぎてexp(x)が+Infになるような入力（+Infを含む）
NaN（exp(NaN)=NaN）

この場合分けを、入力倍精度浮動小数点数の上位32bitを単精度浮動小数点数とみなして絶対値をとってから特定の定数と比較することで行っています。

「倍精度浮動小数点数の上位32bitを単精度浮動小数点数とみなして」というのがいきなり意味不明ですが、これはGPUで倍精度浮動小数点数を持つときのフォーマットを利用した技です。 GPU上で倍精度浮動小数点数を持つときには、隣り合う*232bitレジスタ二つに入れているようなので、これらのレジスタのうち片方を命令オペランドに指定することで、簡単に上位32bitを単精度浮動小数点数として扱えるということです。

倍精度浮動小数点数と単精度浮動小数点数だと指数部のbit数が違う点が気になるかもしれませんが、比較の両辺が正（+0.0含む、NaNを含まない）であればただの辞書順比較と同じなので大丈夫です。

そんなわけで、「元の倍精度浮動小数点数入力の上位32bitを単精度浮動小数点数とみなした時、その絶対値が特定の定数より小さければ、必ずexp(x)が正規化数になる」ように、「特定の定数」が選ばれているようです。実用上は大半の入力でこの条件が満たされるはずで、その場合にearly returnすれば追加の計算コストをほとんど払わずに済みます（この場合、単精度浮動小数点数の絶対値演算と比較演算がそれぞれ一回ずつで済みました）。

また、『元の倍精度浮動小数点数入力の上位32bitを単精度浮動小数点数とみなした時、その絶対値が別の特定の定数より大きい、またはその単精度浮動小数点数がNaNである*3場合、「小さすぎてexp(x)が0.0」「大きすぎてexp(x)が+Inf」「入力がNaNなのでexp(x)もNaN」のどれか』になるように「別の特定の定数」が選ばれているようです。このようなケースは例外ケースの中では頻出なので、ここも最適化テクニックで高速化されていました。具体的には、「exp(x)=x+(+Inf)を返す」という方法で「大きすぎてexp(x)が+Inf」と「入力がNaNなのでexp(x)がNaN」のケースをまとめて処理していました。どんな有限のxに対してもx+(+Inf)は+Infですし、xがNaNであればx+(+Inf)はNaNです*4。

それ以外のパターン、つまり「exp(x)が非正規化数になるような入力x」に加えて、単精度浮動小数点数を使ったがために厳密な場合分けができずに漏れてしまった「exp(x)がかなり大きな正規化数になるような入力x」「exp(x)がかなり小さな正規化数になるような入力x」「exp(x)が最も小さな非正規化数よりもやや小さくて0.0になるような入力x」の場合は、指数部のビット演算を工夫した後、浮動小数点数演算で最終結果を算出していました。これらのパターンが使われることはほとんどないですが、演算の量はそれなりに増えてしまうようです。

*1:実装アルゴリズムによっては、途中の計算結果がオーバーフローするが答えは正規化数になる、といったケースの例外処理が必要なこともあり得ます。

*2:論理レジスタ番号が2nと2n+1のレジスタという意味です。

*3:このような比較はunordered比較と言って一命令で実現可能です。

*4:xが-Infだとx+(+Inf)はNaNになりますが、まずxが負であれば0.0を返すというアルゴリズムなので、この計算結果が選ばれることはありません。

2020-11-15

分岐予測ミスペナルティを測ってみた

近年のほとんどのプロセッサには分岐予測機構が搭載されています。分岐予測機構は、この命令は分岐命令か否か、分岐命令だとしたらとび先はどこか、条件分岐ならそれが成立するかしないか、などを予測します。予測が外れていた場合、正しい命令のフェッチからやり直しになり、パイプラインが埋まるまでには時間がかかります。このパイプラインが埋まるまでにかかる時間を分岐予測ミスペナルティとします。

実験環境

Intel Core i5 7200U（Turbo Boost時3.1GHz）

実験1: リターンアドレスを破壊してみる

以下に示した実験コードでは、スタックに保存されているリターンアドレスを勝手に書き替えています。リターン命令で戻る先を関数コール命令の次ではないところにすることで、分岐予測を外します。

main:
        push rax
        mov  rcx, 1000000000
        call f
        pop  rax
        ret

f:
        sub rcx, 1
        jz L1
        sub QWORD PTR[rsp], 5
L1:
        ret

これを実行するのにかかった時間は10.8秒ほどでした。 3.1GHzで動いていることから、33G cycleほどかかったことがわかります。 1G回ループを実行していることから、分岐予測ミスペナルティは最大で33 cycle程度ということになります。

実際の分岐予測ミスペナルティは、ループの中にcall, sub, retといったメモリに触る命令があるため解析が難しいです。メモリアクセスには4cycleかかるとし、メモリを介した依存が正しく予測されているとすれば、callのストアに4cycle、subのロードに4cycle、subの演算に1cycle、subのストアに4cycle、retのロードに4cycle、で分岐予測ミスペナルティ以外に17cycleかかることになります。よって、分岐予測ミスペナルティ自体は16cycleということになります。

実験2: 疑似乱数を用いて方向予測器を外してみる

以下に示した実験コードでは、疑似乱数（64bitのxorshift）の下1bitを使った条件分岐があります。成立しても成立しなくても結局同じとび先なのですが、正しく予測できなければ分岐予測ミス扱いとなり、ペナルティが発生します。疑似乱数の下1bitを予測することは難しいので、50%の確率で分岐予測ミスが発生します。

main:
        mov     rcx, 2000000000
        mov     rax, 1

L0:
        test    eax, 1
        je      L1
L1:

        mov     rdx, rax
        shl     rdx, 7
        xor     rdx, rax
        mov     rax, rdx
        shr     rax, 9
        xor     rax, rdx

        sub     rcx, 1
        jne     L0
        ret

このコードを実行すると、8.45秒ほどかかりました。 3.1GHzで動いていることから、プログラムの実行に26G cycleかかったことになります。

mov rax, 1の部分をmov rax, 0に変更すると、xorshiftの計算自体は行われるものの常に0が出力され、分岐予測が当たるようになります。この場合には2.75秒かかりました。 3.1GHzで動いていることから、プログラムの実行に8G cycleかかったことになります。ループ一周当たりのレイテンシが4cycleであり2G回のループであるため、予想通りの結果です。

これらの差分から、分岐予測ミスペナルティは合計で18G cycle程度であることがわかります。分岐予測ミスが50%の確率で発生するので、分岐予測ミス回数は1G回です。よって、分岐予測ミスペナルティは18 cycleということになります。

まとめ

Intel Core i5 7200Uの分岐予測ミスペナルティは16～18cycle程度のようです。

2020-11-01

量子コンピュータの勉強（その1）

[1204.4221] Magic-state distillation with the four-qubit codeという論文のFig. 1(c)でやっていることを手計算で試して理解しました。

マジック状態とは

量子コンピュータの勉強をやると、Xゲート・Zゲート・Hゲート・CNOTゲート、などが出てきますが、これだけでは万能な量子計算はできません。

万能な量子計算をするためには、例えばブロッホ球上で45°回すみたいなゲートが追加であれば十分です。

しかし、そのようなゲートは、量子誤り訂正符号上で実現することが困難です。

この問題は、量子テレポーテーションの時と同様のうまい手順で解決することができます。

具体的には、「適切な状態」を事前に準備しておき、（簡単なゲート操作を行った後）適切な観測を行い、その観測結果によってゲートを通したり通さなかったりする、という手順です。

この時に必要となる「適切な状態」のことをマジック状態と呼びます。

以下の図では、マジック状態である $\left|H\right\rangle = \cos\left(\frac\pi8\right)\left|0\right\rangle+\sin\left(\frac\pi8\right)\left|1\right\rangle$ を用いたY軸周り45°回転を行う計算（以下図）をやってみます。

f:id:lpha_z:20201101224330p:plain — マジック状態|H>を用いたY軸周り45°回転手順（arXiv:1204.4221 の図1より一部抜粋）

落とし穴

論文中に書いてありますが、単独のYゲート操作と、角度パラメータ付きの $Y(\theta)$ ゲート操作は別物です。

計算

上の量子ビットである $\left|H\right\rangle$ を、一般的に $\left|H\right\rangle = c\left|0\right\rangle+s\left|1\right\rangle$ と置いてみます。

下の量子ビットを $\left|H\right\rangle = a\left|0\right\rangle+b\left|1\right\rangle$ と置きます。

ここで、基底を計算基底（ $\left|0\right\rangle$ や $\left|1\right\rangle$ を用いた表現）からYの固有状態（ $\left|+y\right\rangle = \frac{ 1}{ \sqrt 2} \left|0\right\rangle + \frac{ i}{ \sqrt 2 }\left|1\right\rangle$ と $\left|-y\right\rangle = \frac{ 1}{ \sqrt 2} \left|0\right\rangle - \frac{ i}{ \sqrt 2 }\left|1\right\rangle$ ）を使った基底に変換します。これは、Yで測定するとYの固有状態のどちらかがその係数の絶対値の二乗の確率で出てくるため、その計算をやりやすくするためです。

すると、 $ca\left|0\right\rangle\left|0\right\rangle-isb\left|1\right\rangle\left|0\right\rangle+cb\left|0\right\rangle\left|1\right\rangle+isa\left|1\right\rangle\left|1\right\rangle$ $= \frac{ca - sb}{\sqrt 2}\left|+y\right\rangle\left|0\right\rangle+\frac{ca + sb}{\sqrt 2}\left|-y\right\rangle\left|0\right\rangle+\frac{cb + sa}{\sqrt 2}\left|+y\right\rangle\left|1\right\rangle+\frac{cb - sa}{\sqrt 2}\left|-y\right\rangle\left|1\right\rangle$ です。

これをYで測定して固有状態 $\left|+y\right\rangle$ が出た場合（確率1/2）、下の量子ビットは $(ca - sb)\left|0\right\rangle+(cb + sa)\left|1\right\rangle$ になっています。これは、上の量子ビットが $\left|H\right\rangle$ だった場合、Y軸周りに45°回転したことになります。したがって、特に追加のゲート操作は不要です*1。

一方、固有状態 $\left|-y\right\rangle$ が出た場合（確率1/2）、下の量子ビットは $(ca + sb)\left|0\right\rangle+(cb - sa)\left|1\right\rangle$ になっています。これは、上の量子ビットが $\left|H\right\rangle$ だった場合、-Y軸周りに45°回転したことになります。つまり、Y軸周りに-45°回転したことになります。したがって、追加でY軸周りに90°のゲート操作を行えば、合計でY軸周りに45°回転したことになる、という具合です。

結局、どちらの固有状態が出た場合でも、その観測結果に応じて適切にゲートを通したり通さなかったりすることで、ブロッホ球上でのY軸周り45°回転を行うことができました。

ちなみに、二つの測定結果はどちらも半々の確率で出ることから、測定結果から何らかの情報を得ることはできず、下の量子ビットに関する情報は漏れていない（＝下の量子ビットを破壊していない）ということも確認できます。

*1:論文ではpositive mesurementの場合にY(π/2)ゲートをかけるとしていましたが、逆だと思います。

2020-10-25

Intel MKLのSIMD exp関数を読んだ

Intel CPUで高速に動く、数学関数のSIMD実装に Intel Math Kernel Library というものがあります。

それに含まれる指数関数実装を読んでみました。

AVX-512（512bitレジスタがあるので、doubleを最大で八並列で演算できる）用の実装は、以下のようになっていました*1。

00000000004013d0 <__svml_exp8_z0>:
  4013d0:  f3 0f 1e fa                     repz nop     edx
  4013d4:  62 f1 7c 48 10 25 62 29 00 00   vmovups      zmm4 , [rip+0x2962]
  4013de:  62 f1 7c 48 10 0d 98 29 00 00   vmovups      zmm1 , [rip+0x2998]
  4013e8:  62 f1 7c 48 10 15 ce 29 00 00   vmovups      zmm2 , [rip+0x29ce]
  4013f2:  62 f1 7c 48 10 1d 04 2a 00 00   vmovups      zmm3 , [rip+0x2a04]
  4013fc:  62 f1 7c 48 10 2d 7a 2a 00 00   vmovups      zmm5 , [rip+0x2a7a]
  401406:  62 71 7c 48 10 2d b0 2a 00 00   vmovups      zmm13, [rip+0x2ab0]
  401410:  62 f2 fd 78 a8 e1               vfmadd213pd  zmm4 , zmm0 , zmm1 , {rz-sae}
  401416:  62 f1 7c 48 10 35 e0 2a 00 00   vmovups      zmm6 , [rip+0x2ae0]
  401420:  62 71 7c 48 10 0d 16 2b 00 00   vmovups      zmm9 , [rip+0x2b16]
  40142a:  62 71 7c 48 10 1d 4c 2b 00 00   vmovups      zmm11, [rip+0x2b4c]
  401434:  62 71 7c 48 10 25 02 28 00 00   vmovups      zmm12, [rip+0x2802]
  40143e:  62 71 dd 18 5c f1               vsubpd       zmm14, zmm4 , zmm1 , {rn-sae}
  401444:  62 72 dd 48 7f 25 32 28 00 00   vpermt2pd    zmm12, zmm4 , [rip+0x2832]
  40144e:  62 d2 ed 18 bc c6               vfnmadd231pd zmm0 , zmm2 , zmm14, {rn-sae}
  401454:  62 d2 e5 18 bc c6               vfnmadd231pd zmm0 , zmm3 , zmm14, {rn-sae}
  40145a:  62 f1 fd 48 54 3d dc 29 00 00   vandpd       zmm7 , zmm0 , [rip+0x29dc]
  401464:  62 71 c5 18 59 c7               vmulpd       zmm8 , zmm7 , zmm7 , {rn-sae}
  40146a:  62 72 d5 18 b8 ef               vfmadd231pd  zmm13, zmm5 , zmm7 , {rn-sae}
  401470:  62 72 cd 18 b8 cf               vfmadd231pd  zmm9 , zmm6 , zmm7 , {rn-sae}
  401476:  62 71 bd 18 59 d7               vmulpd       zmm10, zmm8 , zmm7 , {rn-sae}
  40147c:  62 72 bd 18 a8 df               vfmadd213pd  zmm11, zmm8 , zmm7 , {rn-sae}
  401482:  62 52 bd 18 a8 e9               vfmadd213pd  zmm13, zmm8 , zmm9 , {rn-sae}
  401488:  62 52 ad 18 a8 eb               vfmadd213pd  zmm13, zmm10, zmm11, {rn-sae}
  40148e:  62 52 9d 18 a8 ec               vfmadd213pd  zmm13, zmm12, zmm12, {rn-sae}
  401494:  62 d2 95 18 2c c6               vscalefpd    zmm0 , zmm13, zmm14, {rn-sae}
  40149a:  c3                              ret

この実装では、 $\exp(x)$ を求めるとき、まず $\frac{x}{\log 2}$ を1/16の位まで求めます。この値の整数部分は、最後のvscalefpdで使用されます。小数部分は、vpermt2pdを使った16エントリの表引きに使われます。

$x$ を $\frac{1}{\log 2}$ で割ったあまりを $t$ とすると、絶対値は $\frac{\log 2}{16}$ 未満のはずなので、 $\exp(t)$ を六次の多項式を使って近似します。

vandpdは、 $x$ の絶対値が大きすぎて $\frac{x}{\log 2}$ を1/16の位まで求める計算がおかしくなった時に $t$ がおかしな値となり、最終結果がおかしくなることを防ぐ役割があるようです。 vandpdで使われるマスク定数は0xbfffffffffffffffであり、絶対値が2以上の数値やNaN・±Infを絶対値が2未満の数に変更する効果があります。 $x$ の絶対値が大きすぎる場合はどうせvscalefpdで0か+Infになるので、適当に小さくすればよいということのようです。

こうして見てみると、vperfmt2命令やvscalef命令は指数関数を実装するためにあるという感じの命令ですね……。

*1:マスクレジスタを使わないことが確定している場合の実装です。

2020-10-18

Intel CPU の分岐予測器を調べてみた（ループの場合だけ）

なんだか分岐予測器で当てれそうなところを当ててくれなかったことがあったので、Intel CPUの分岐予測器はどんな時に分岐予測を当てることができるのかを調べてみました。

実際に調べてみると、ループのような単純な場合であっても非常に複雑怪奇な挙動を示すため、内部構造がどうなっているか全く推測できませんでした。

とりあえず、結果だけまとめておきます。

実験方法

単純なループを10億回くりかえすだけの以下のアセンブリコードを実行したときの分岐予測ミス回数を計測しました。

        .globl  main
main:
        xorl      %eax, %eax
L6:
        xorl      %edx, %edx
L7:
        incl      %edx
        cmpq      $__, %rdx
        jb        L7
        incl      %eax
        cmpl      $1000000000, %eax
        jb        L6
        ret

ループ回数は、以下のスクリプトで設定しました。

for i in $(seq 1 100); do
 sed s/__/$i/ main.s > m.s
 gcc m.s
 perf stat ./a.out
done

実験に使用したCPU

ここ10年くらいのXeonで実験することができました。

Xeon X5550（Neharem EP、2009年）
Xeon E5-2690（Sandy Bridge EP、2012年）
Xeon E5-2643 v3（Haswell EP、2014年）
Xeon E5-1620 v4（Broadwell EP、2016年）
Xeon Gold 5115（Skylake SP、2017年）
Xeon Gold 6126（Skylake SP、2017年）
Xeon platinum 8280（Canscadelake SP、2019年）

結果

Xeon X5550（Neharem EP、2009年）

ループが64周以下ならば、分岐予測ミスは発生しない。
ループが65周以上であれば、（おそらくループ脱出時に）分岐予測ミスが1回発生する。

分岐履歴長が63であることを示唆しているように思えます。

非常に素直な結果になりました。

Xeon E5-2690（Sandy Bridge EP、2012年）

ループが33周以下ならば、分岐予測ミスは発生しない。
ループが34周以上であれば、（おそらくループ脱出時に）分岐予測ミスが1回発生する。

分岐履歴長が32であることを示唆しているように思えます。

非常に素直な結果ですが、なぜ3年前のプロセッサよりも劣化しているのでしょうか。ループ以外の分岐予測の精度が上がったと信じたいです。

Xeon E5-2643 v3（Haswell EP、2014年）

ループが35周、40周、42周、44周、の場合にはループ一周につき1回以上の分岐予測ミスが発生する。
ループが53周、54周、55周、60周～74周、91周～94周の場合はループ一周につき1回未満だが0.1回以上の分岐予測ミスが発生する。
それ以外の場合でループが59周以下であれば、分岐予測ミスはほとんど発生しない。
ループが75周～90周の場合、ループ一周につき1.2回程度の分岐予測ミスが発生する。
ループが95周以上であれば、（おそらくループ脱出時に）分岐予測ミスが1回発生する。

わけの分からない結果になりました。

Xeon E5-1620 v4（Broadwell EP、2016年）

やるたびに結果が変わってよくわかりませんが、

ループが35周、40周、42周、44周の場合はループ一周につき1回程度の分岐予測ミスが発生しやすい
ループが40周以上だと当てられなくなる時と、60周以上だと当てられなくなる時がある

ループ35周、40周、42周、44周の時に当てられないという挙動が、Xeon E5-2643 v3（Haswell EP、2014年）と酷似しています。

Xeon Gold 5115（Skylake SP、2017年）と Xeon Gold 6126（Skylake SP、2017年）

やはりやるたびに結果が変わってよくわかりませんが、

ループが35周未満であれば、分岐予測ミスはほとんど発生しない
ループが35周のとき、分岐予測ミスが発生しやすい

やはり35周のところに何かがあるようです。

Xeon platinum 8280（Canscadelake SP、2019年）

ループが23周以下ならば、分岐予測ミスは発生しない。
ループが24周～28周の場合はループ一周につき1回未満だが0.05回以上の分岐予測ミスが発生する。
ループが29周以上であれば、（おそらくループ脱出時に）分岐予測ミスが1回発生する。

素直な結果に戻りましたが、7年前のプロセッサよりもさらに劣化しています。 "improved branch predictor" とはいったい何だったのでしょう……

まとめ

Intel CPU にはループ予測器はついていなさそう
ループ脱出の予測に関しては、Intel CPUはここ10年で全く進化していない（どころか当てられるループ周回数が減少している）
30周程度のループでも当てられないことがある

35周の時に当てられない現象は、分岐命令のアドレスなどを変えて再検証する必要あり

2020-10-04

「単純和」プログラムの速度をJuliaとC言語で比較した

先週、浮動小数点数を足していくプログラムについて、Juliaで書いた場合とC言語で書いた場合とで速度を比較しました（浮動小数点数を足すだけのプログラムがC言語よりJuliaのほうが速い、について - よーる）。

この記事について、黒木玄さん（黒木玄 Gen Kuroki (@genkuroki) | Twitter）から「Juliaでは単純な和には普通にsimdを使うので、JuliaとCの両方でsimdを使って比較した方が良かったと思います。」というコメントをいただきました。先週の記事で比較に用いたプログラムには並列性がないため、SIMD (single instruction multiple data) 命令を使った高速化は不可能なはずですが、Juliaでsimdを使うと10倍高速化するそうです。一体どういうことなのだろうと思ってJuliaのマニュアル（Base.SimdLoop.@simd (macro)）を見てみると、@simdというアノテーションは浮動小数点数演算の順序を変える効果があるそうです。

浮動小数点数演算は一般に結合則が成り立ちません。したがって、浮動小数点演算の順序を変えてしまうと、プログラムの意味が変わってしまいます（一般に計算結果も変わってしまいます）。先週の記事で比較したプログラムはあくまで「浮動小数点数を前から順番に足していく」プログラムであって、「いい感じに浮動小数点数の和を求める」プログラムではありません。

この辺は、多くの読者が勘違いされているところだと感じました。

とはいえ、「浮動小数点数を前から順番に足していった値を正確に求めたい」というワークロードが常識的ではないのも事実です。そこで、「浮動小数点数を順番を気にせずに足した値を求める」というワークロードで、JuliaとC言語の速度比較をしてみます。

使ったソースコード

ちょうど黒木玄さんがJulia側のソースコードを示してくださいました（simd example · GitHub）ので、それを使うことにします*1。ただし、実行時間が短すぎて正確な評価ができなかったので、外側ループの回数を100倍にしました。

C言語側のソースコードは、私が素直に移植したものになります。

function test_inbounds_simd()
    a = fill(1e-10, 1000)
    s = 0.0
    for j = 1:1000000000
        @inbounds @simd for i = 1:1000
            s += a[i]
        end 
    end
    s
end
test_inbounds_simd()

double test() {
  unsigned long long i, j;
  double a[1000];
  double s = 0.0;
  for( i = 0; i < 1000; ++i ) { a[i] = 1e-10; }

  for( j = 0; j < 1000000000; ++j ) {
    for( i = 0; i < 1000; ++i ) {
      s += a[i];
    }
  }

  return s;
}

int main() {
  return test();
}

使用したハードウェア

Intel(R) Xeon(R) Platinum 8280 CPU @ 2.70GHz
- 実際には3.6GHzくらいで動いているっぽい(?)
- AVX-512使用可能（一命令で倍精度浮動小数点数8個を同時に演算可能）
- FMA実行ユニット2基（1cycleに浮動小数点数加算命令が2回実行できる）
- 浮動小数点数加算命令のレイテンシは4cycle

使用したソフトウェア

julia version 1.5.2（内部で使用しているLLVMのバージョンはLLVM version 9.0.1jl）
clang version 10.0.0
icc (ICC) 19.0.5.281 20190815
perf
gdb

理論限界

ベンチマークプログラムは1000G個の浮動小数点数を加算するものです。実験に使用したプロセッサは1サイクルに16個の浮動小数点数を加算することができますので、理論限界は62.5Gサイクルです。 3.6GHzで動いているとすると17秒程度かかることになりますが、周波数が安定しないので実時間は参考程度ということにします。

Juliaの結果

199G cycles
278G insns
43G branches
55秒くらい

JIT コンパイルされた機械語をgdbで確認してみると、以下のようになっていました。

   0x00007fffc565f0d0:  cmp    rax,0x3b9aca00
   0x00007fffc565f0d6:  lea    rax,[rax+0x1]
   0x00007fffc565f0da:  je     0x7fffc565f187
   0x00007fffc565f0e0:  vmovq  xmm0,xmm0
   0x00007fffc565f0e4:  vxorpd xmm1,xmm1,xmm1
   0x00007fffc565f0e8:  mov    esi,0x18
   0x00007fffc565f0ed:  vxorpd xmm2,xmm2,xmm2
   0x00007fffc565f0f1:  vxorpd xmm3,xmm3,xmm3
   0x00007fffc565f0f5:  nop    WORD PTR cs:[rax+rax*1+0x0]
   0x00007fffc565f0ff:  nop
   0x00007fffc565f100:  vaddpd zmm0,zmm0,ZMMWORD PTR [rcx+rsi*8-0xc0]
   0x00007fffc565f108:  vaddpd zmm1,zmm1,ZMMWORD PTR [rcx+rsi*8-0x80]
   0x00007fffc565f110:  vaddpd zmm2,zmm2,ZMMWORD PTR [rcx+rsi*8-0x40]
   0x00007fffc565f118:  vaddpd zmm3,zmm3,ZMMWORD PTR [rcx+rsi*8]
   0x00007fffc565f11f:  add    rsi,0x20
   0x00007fffc565f123:  cmp    rsi,0x3f8
   0x00007fffc565f12a:  jne    0x7fffc565f100
   0x00007fffc565f12c:  vaddpd zmm0,zmm1,zmm0
   0x00007fffc565f132:  vaddpd zmm0,zmm2,zmm0
   0x00007fffc565f138:  vaddpd zmm0,zmm3,zmm0
   0x00007fffc565f13e:  vextractf64x4 ymm1,zmm0,0x1
   0x00007fffc565f145:  vaddpd zmm0,zmm0,zmm1
   0x00007fffc565f14b:  vextractf128 xmm1,ymm0,0x1
   0x00007fffc565f151:  vaddpd xmm0,xmm0,xmm1
   0x00007fffc565f155:  vpermilpd xmm1,xmm0,0x1
   0x00007fffc565f15b:  vaddsd xmm0,xmm0,xmm1
   0x00007fffc565f15f:  test   dl,dl
   0x00007fffc565f161:  jne    0x7fffc565f0d0
   0x00007fffc565f167:  xor    esi,esi
   0x00007fffc565f169:  nop    DWORD PTR [rax+0x0]
   0x00007fffc565f170:  vaddsd xmm0,xmm0,QWORD PTR [rcx+rsi*8+0x1f00]
   0x00007fffc565f179:  inc    rsi
   0x00007fffc565f17c:  cmp    rsi,0x8
   0x00007fffc565f180:  jne    0x7fffc565f170
   0x00007fffc565f182:  jmp    0x7fffc565f0d0

この機械語コードから予測すると、1000個の浮動小数点数を足すのに

192cycle（zmmレジスタを使った主ループが4cycle×31回、マージしていく部分が36cycle、端数処理ループが4cycle×8回）
273命令（準備に11命令、主ループが7命令×31回、マージしていく部分が13命令、端数処理ループが4命令×8回）
42分岐命令

が必要です。

これはperfコマンドの結果と整合します。

理論限界と大きく差がついている理由は、以下の通りです。

主ループのアンローリングが不足している
- 4cycleレイテンシを隠蔽しつつFMAユニット2基をフル稼働させるためには、依存関係にない命令が最低でも8つ必要
端数処理が直列化されている
- zmmレジスタを使って並列化可能

C言語の結果（clangでコンパイルした場合）

コンパイルオプションは、-O2 -mavx512f -funsafe-math-optimizationsです*2。 C言語の範囲では「浮動小数点数演算の順序を自由に変更してかまわない」という記述をすることができないので、-funsafe-math-optimizationsをつけます*3 *4。

194G cycles
213G insns
32G branches
55秒くらい

Juliaで書いた場合と同じくらいの速度になりました。それもそのはずで、JuliaのJITエンジンはLLVM、つまりclangのバックエンドと同じです。というわけでコード生成結果や実行速度が大きく変わるはずはないのです。

命令数が少なくなっているのは、JIT コンパイルより時間をかけてコンパイルすることできれいなコードを生成できたことによるものと思われます。

.LBB0_5:
        vaddpd  zmm0, zmm2, zmm0
        vaddpd  zmm0, zmm3, zmm0
        vaddpd  zmm0, zmm4, zmm0
        vextractf64x4   ymm2, zmm0, 1
        vaddpd  zmm0, zmm0, zmm2
        vextractf128    xmm2, ymm0, 1
        vaddpd  xmm0, xmm0, xmm2
        vpermilpd       xmm2, xmm0, 1
        vaddsd  xmm0, xmm0, xmm2
        vaddsd  xmm0, xmm0, xmm1
        vaddsd  xmm0, xmm0, xmm1
        vaddsd  xmm0, xmm0, xmm1
        vaddsd  xmm0, xmm0, xmm1
        vaddsd  xmm0, xmm0, xmm1
        vaddsd  xmm0, xmm0, xmm1
        vaddsd  xmm0, xmm0, xmm1
        vaddsd  xmm0, xmm0, xmm1
        add     rax, 1
        cmp     rax, 1000000000
        je      .LBB0_6
.LBB0_3:
        vmovq   xmm0, xmm0
        vxorpd  xmm2, xmm2, xmm2
        mov     ecx, 56
        vxorpd  xmm3, xmm3, xmm3
        vxorpd  xmm4, xmm4, xmm4
        .p2align        4, 0x90
.LBB0_4:
        vaddpd  zmm0, zmm0, zmmword ptr [rsp + 8*rcx - 576]
        vaddpd  zmm2, zmm2, zmmword ptr [rsp + 8*rcx - 512]
        vaddpd  zmm3, zmm3, zmmword ptr [rsp + 8*rcx - 448]
        vaddpd  zmm4, zmm4, zmmword ptr [rsp + 8*rcx - 384]
        cmp     rcx, 1016
        je      .LBB0_5
        vaddpd  zmm0, zmm0, zmmword ptr [rsp + 8*rcx - 320]
        vaddpd  zmm2, zmm2, zmmword ptr [rsp + 8*rcx - 256]
        vaddpd  zmm3, zmm3, zmmword ptr [rsp + 8*rcx - 192]
        vaddpd  zmm4, zmm4, zmmword ptr [rsp + 8*rcx - 128]
        add     rcx, 64
        jmp     .LBB0_4

zmmレジスタを使った主ループが二倍アンローリングされたほか、端数ループが八倍（完全）アンローリングされています。これらが命令数削減に役立っているようです。

この機械語コードから予測すると、1000個の浮動小数点数を足すのに

192cycle（zmmレジスタを使った主ループが8cycle×15.5回、マージしていく部分が36cycle、端数処理が4cycle×8回）
213命令（ループの外が27命令、ループの中が12命令×15.5回）
32分岐命令

が必要です。

これはperfコマンドの結果と整合します。

理論限界と大きく差がついている理由は、Juliaの場合と全く同様です。

どうでもいいですが、配列の最後のほうがループ内定数とみなされてxmmレジスタに乗っていますね……。これは微妙にずるな気がしますが計算時間や命令数に変わりはないのでよしとしましょう。

C言語の本気（iccでコンパイルした場合）

C言語の強みは何といっても ハードウェアを作っている会社が専用コンパイラを提供してくれる ことにあります。

Intel社が作ったコンパイラ、iccを使ってコンパイルしてみましょう。コンパイルオプションは、-axCORE-AVX512 -qopt-zmm-usage=highです。

結果は以下の通りです。

62.7G cycles
172G insns
15.6G branches
18秒くらい

理論限界の62.5Gサイクルにかなり近い速度が達成できていることがわかります*5。

最内ループを確認してみると、お互いに依存しない8個のvaddpd命令になっていました。さすがにこの辺の最適化はちゃんとしていますね。

まとめ

Juliaの速度はclangでコンパイルしたC言語と同程度
「単純和」プログラムをC言語で書いてiccでコンパイルした場合、clangでコンパイルした場合やJuliaで書いた場合より3倍くらい速い
コンパイラが自動でやってくれる程度のSIMD化をする手間は、JuliaもC言語も大して変わらないと思う

*1:この手の言語間の速度比較では、不適切なソースコードが使われることにより特定の言語の性能が不当に低く出てしまうということがありがちですが、そういうことにならない有意義な比較ができることを感謝いたします。

*2:-O3や-Ofastだと完全にアンローリングされてしまって見づらかったので-O2にしました（実行にかかるサイクル数は同じでした）。

*3:-funsafe-math-optimizationsは-ffast-mathに含まれるので、-ffast-mathでもよいです。

*4:-funsafe-math-optimizationsは複数のオプションを有効にしますが、実際に役立っているのは-fassociative-math -fno-signed-zerosのようです。

*5:ループ交換が行われているので本来とは全く異なる順序で計算されていますが、加算は1000G回行われています。