トップページへ戻る
過去の日記へ


予定


2月21日(水)

 光成さんによるx86/x64最適化勉強会8 レポートが出ています。 光成さんもそうですが、「x86/x64最適化勉強会」というと人外ガチ勢の集会というイメージだったので、そこで 僕が発表できたということは、少しでも人外ガチ勢に近づけた気がして嬉しいです。

 僕がまともにx86のSIMD化を始めたのは2016年9月のこの試行が初めてだと思うので、まだ2年たってない。2011年くらいにxmmを使ったSIMD化を試みているが、メモリアクセスがボトルネックでどうにもならなかった感じ。

 そういえば勉強会で「これ、業績になるんですか?」ってストレートに聞かれたけど、「ならないんじゃないですかね・・・」と答えておいた。 パディング付きAoS4要素パックとかは既に報告があるらしいし、新規性があるのは8要素パックくらい? でもまぁ新規性はともかく、実装の詳細を紹介した論文はどこかに通しておきたい。SIMD化の歴史の一つとして。 でもいま抱えてる論文がちっとも進まず、次の論文に取り掛かれないのであった。

 というわけで全く筆が進まない論文に取り組む。 ちなみに文献管理にJabRef使ってて、doi一発で取ってこられるから便利。さらに文献の省略ー展開が一発でできるのもいい感じなんだけど、 なぜかdoiで取ってきた文献の「The Journal of Chemical Physics」が省略できない。他の奴はできるのに。 色々調べたら「Journal of Chemical Physics」とTheを取ったらできるようになった。ふむ。


2月20日(火)

 少しずつでも頑張ろう。

 ローカルとリモートの区別がなく、ディレクトリ単位でコミットできて、全ての履歴が完全に一本道というSubversionの性質、わりと好きなんだけど、 .svnignoreの扱いだけは我慢ならない。gitみたいにデフォルト.svnignoreがあったらそれを参照、でいいじゃん。 もしくはそういう設定ができてしかるべきでしょ。 いちいち.svnignoreいじるたびに「svn propset svn:ignore -R -F .svnignore . 」みたいなコマンド打ってられないし、 そもそも覚えてなくて、毎回「どうやって書くんだっけ」って調べるハメに。


2月19日(月)

 土曜日の勉強会の配信動画、本当にきれいだな。良い機材を使ったとのこと。 後、スライドと講演者を両方同時に写しているのもすごい。他の人がほとんど動かなかったのに、 僕だけ歩き回ってて、撮影の人に申し訳ないことをしたな。あと、スライドの文字が細かいな、と思われる時には スライドだけを大写しにしたりとかをリアルタイムにやっててすごい。

 あと、勉強会でMeltdownとSpectreの講義があったのがすごく勉強になった。スパコン運用スタッフとしても非常にタイムリーだった。

 言い訳は見苦しいだけだとわかっているんだけど、最近は本当にいろいろ苦しい。これが僕の実力じゃないんですよ、と声を大にして言いたい(←今出ている結果を実力と呼びます)。なんとか立て直さないと。

 しかしサイボウズの東京オフィス、すごかったな。Googleにも感じるけど、なんというか、「優秀な奴を集めてうまいことマネタイズできていれば、後は好き勝手していい」みたいな感じで、その自由な雰囲気でまた優秀な奴が集まって、みたいなサイクルができてるんだろうな・・・。なんというか、ある意味「力で殴る」みたいな雰囲気。上手く表現できないが。

 翻って自分はどうなんだろうな。研究者としてにせよ、エンジニアとしてにせよ、「自分にどういう市場価値があるか」もっと言えば「どういう市場価値をどうやってつけるか」という問題か・・・。π型戦略といえば聞こえはいいけど、どっちつかずになる危険性と隣合わせなんだよな。

 あと、タスク管理が結局うまく言っていない。予め「この日の午前中はこれ、午後はこれ」くらいの粒度でタスクを割り当てるのが理想なんだけど、 そうなっておらず、結局朝に「その日できそうなタスク」を並べて、それらを消化しているだけになっている。 Todoistではなく、カレンダーベースのタスク管理の方がいいのかも。

 ・・・と思って調べてみたら、TodoistにGoogle Calendarとの連携機能があったのでそれを設定してみた。 これで一週間単位くらいでタスク管理してみて様子見だな・・・。

 「キラキラツールで生産性向上」とかそういうライフハック系、あまり好きじゃないんだけど、現状タスクがうまく回ってないのでいろいろ試すしかない。 ちなみに、仮にこれでうまく行ったとしても、「公私共に忙しいエンジニアが今すぐ導入すべきツール10選」みたいな記事は書きませんよ。

 そうそう、x86勉強会でもう一つ驚いたのはtanakmuraさんのこれ(4.6 KNLのFPUスループット)。KNLで二つあるFPU(VP0, VP1)のうち、VP1が5cycleに1cycle停止しており、かつスケジューラがそれを考慮に入れず命令を突っ込むから 結果として性能が落ちたように見えるという仮説。これが本当だとするとすごい。

 っていうかこういうの個人で調べちゃうの、すごいよな・・・。Ryzen SEGVバトルでも 命令実行アドレスがずれていることを突き止めた人とか。 なんかこの界隈の(に限らないけど)ガチ勢ってのは本当に怖い。

 ダメだ。絶不調だ。


2月18日(日)

 NOP


2月17日(土)

 x86/x64最適化勉強会に参加。

自分の発表スライドはこちらです→「AVX2/AVX-512を用いたLennard-Jones系ポテンシャルの力計算のSIMD化

 動画も配信されてたみたい。いまでもYouTubeで勉強会の様子が見られます。僕の発表はこのあたりからです。

 なんというか「こんにちは、物性研の渡辺です。こちらは物性研の公式マスコットキャラクター、その名も物性犬です」のくだりでだいたいひと笑い取れるから便利だ>物性犬

 某WGのためのマルチスレッド環境下のmalloc調査と、この勉強会ためのAVX2/AVX-512のSIMD化に相当時間をかけてきたのだが、どちらもまぁまぁ発表がうまくいったようでよかった。切り替えて次だ。


2月16日(金)

 体調がいまいち。


2月15日(木)

 C++でunordered_map使うたびに「contains」的なメソッドを書いている気がするんだけど、なにこれ? containsくらいデフォルトで実装してもバチは当たらないだろ、ともう5回くらい曰記に書いている気がしたが、 調べたらまだ一回しか書いてない?しかもこれはvectorの話か。

 あと、なんでstd::map::findはあってstd::vector::findは無いのよ?なんでよ? C++の言語設計者の美意識がわからん。

 あー、std::vectorでbool配列作った時にハマったことも思い出したぞ。 STLにはだいぶヘイトが溜まっている。僕はC++は好きじゃないんですよ。正直言って仕方なく使ってる。 Juliaでもなんでもいいんだけど、とりあえずそれなりの表現力があり、必要に応じてマシン語レベルのチューニングができるような 言語が出てきてくれれば・・・(他力本願)。

 そういえば最近のFortranどうなんだろ?少なくとも僕が知ってる頃はクラスのインスタンスが作れなかった(static扱い?)んだけど、 今調べたら2003からクラスのインスタンスが作れるようになった?文字列の扱いとかはもう少し楽になったのかしら。

 今のところ使ってて違和感が無いのはやっぱりRubyかなぁ。Pythonはstrやlen、joinとかの扱いが気持ち悪い(とか言うとPythonistaから袋叩きにあうらしいが)。


2月14日(水)

 Toggl使ってたら「in a jiffy」という言葉が出てきて、なんだそりゃ?と思って調べて見た。 とりあえず「すぐに」という意味なのだが、「jiffy」という言葉がわからない。 Online Etymology Dictionaryによると、 「origin unknown」とのことだが、もともと泥棒用語で、「lightning」を指す隠語だったという説が紹介されている。

 StackExchangeでは、 1780年の文献に「in a jiffy」の用法を見つけた、という報告がある。ここには「jeffy」「jiffey」といった他のバリアントも報告してあって面白い。 いずれにせよ、語源がはっきりしない単語のようだ。

 findの使い方をまったく覚えられず、使うたびにググっている。 で、今日もfindでググって出てきた「Linuxコマンド集 - 【 find 】 ファイルやディレクトリを検索する:ITpro」をクリックしたら タイトルだけでてきて「ここからは会員の登録が必要です。」と言われ、ひとつため息をついて「man find」した。


2月13日(火)

 一回休み。


2月12日(月)

 NOP


2月11日(日)

 やれやれ。


2月10日(土)

 NOP


2月9日(金)

 ようやく時間が取れて、別プロジェクトが少し進んだ。

 HaswellとSkylake、1ペアの力計算に何サイクルかかってるか計算してみるかな。 実際には周波数が揺らぐからあまり意味ないかもしれないけど。

 今、ペアの数は7839886。100回計算しているから、力の計算は783988600回している。 Haswellの最速実装ではそれに3309かかっているから、 ベースクロックが2.5GHzなんで8315000000.0サイクルくらいですか。 ってことは、10サイクル/ペア?あれ?思ったより早い。計算間違えたかな? LJの計算って加減乗算27回+除算1回で、1ペアあたり24バイトのloadと24バイトのストアがあるわりには結構早いな。 SIMD化がそれなりにうまくいってるってことか。

 同様に計算するとSkylake(Intel(R) Xeon(R) Gold 6148 CPU @ 2.40GHz)で7.8サイクル/ペア。 うーむ。理論演算性能がほぼ倍になったわりには・・・というところか。でもまぁ、AVX512でもっと動作周波数ががっつり 下がってると考えると、ざっくり5サイクル/ペアになって、倍くらいになる勘定。 でもターボブーストがかかってると考えるともっとサイクル数は増える。 実際、ログインノードだとベース2.1GHzなのに3.456GHzとか出る(perfによれば、だけど)。

 ちなみにXeon Phiは12.33サイクル/ペア。Xeonに比べて倍悪くなってないのは結構がんばってる気がする?

 昔、POWER6を使ってた頃にペアあたり50〜100サイクルと見積もったんだけど、だいぶ進化したなぁ。 それともやっぱりPOWER6が短距離古典MDに向いてなかったのか・・・。


2月8日(木)

 時間が取れない。


2月7日(水)

 いろいろ無駄に忙しい。lj_simdにavx2+swpがなかったので追加した。 Haswellではそれなりに効果があるのだが、Skylakeではあまり効果が無いな。


2月6日(火)

 数学の証明の最後に出てくる四角い記号、墓石記号とかハルモス記号とか呼ぶらしい。ハルモスはアメリカの数学者で、彼によれば「ハルモス記号そのものは私の発明ではない」そうだが、 もともと記事の最後を表す記号として使われていた□を「証明完了」の意味で使い始めたのは彼らしい。 彼の数学の教科書は広く読まれたらしいので、彼が広めたのは間違いなさそう。 あと、「if and only if」を「iff」と書くのも彼の発明らしい。

 なんか組み込み関数でAVX2やAVX-512を使っている場合、ループ内の独立なはずの行を手で入れ替えるとかなり性能に影響がある。 コンパイラは自分が吐いていないAVX2/AVX-512命令の順序の入れ替えができない or やりづらい?

 人外集会x86/x64最適化勉強会のスライドできた。わりとがんばって準備した。これ、論文になるかしら?


2月6日(火)

 昨日の夜、AVX-512のレジスタ8本使って上位/下位の転置をしてバラして計算するコードをがんばって書いたが、さほど早くならなくて悲しい。 だいぶSIMD化に慣れてきた気がするが、キリが無いので現在KNLにおける最速実装と思われるgather/scatter+swpと、 gatherだけしてscatterしない版を試しておしまいにしよう。

 しかしインテルコンパイラの「v4dfとv8dfの混ざった含む関数オーバーロードをするとコンパイルできなくなる」バグが地味に不便。 これはGCCとの互換性確保のためなので、新しめのGCCにパスが通っていれば回避可能なのだが、いちいちmodule loadしないといけないのが面倒。 もしくは-no-gccでも回避可能だが、そうするとv4dfやv8dfが使えなくなって四則演算書くのが面倒に。

 gatherはするけど、scatterしないで転置してAoSで書き戻すルーチン書いてみたけど遅くて悲しい。 あと、書き戻しにマスクが必要だったので_mm256_mask_store_pd使ったらKNLでIllegal instructionと怒られた。そっか、そっちにはVLが無いのか。 せっかくのAVX-512なのに、SkylakeではAVX2実装が最速、ってのもちょっと悲しいなぁ。 要するにAVX-512を使うメリットと、その際に下がる動作周波数分の損益分岐点がどこか、という話だと思うんだけど、 「-qopt-zmm-usage=high」みたいなオプションがあって、デフォルトがlowな時点で、インテル側もあまりAVX-512で早くなるか自信が無いよ、ってことなんだろうけど。

 PCの見積もり依頼出した。

 よし、とりあえず現在知られているKNLでの最速実装を、ちゃんと自分の手で再実装したぞ。 ついでにSoAにもループの端数処理最適化+SWPを実装した。速度も中川さんの(逆数近似を使わない版の)最速実装とほぼ同じになった(本質的に同じコードだから当たり前だが)。 最終版(881c7a9495a024b1d96f3b28be763c9da46e9571)の結果を載せておこう。

# Oakforest-PACS 
# icpc (ICC) 18.0.1 20171018
# Intel(R) Xeon Phi(TM) CPU 7250 @ 1.40GHz
$ ./aos.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 37117 [ms]
N=119164, sorted 10431 [ms]
N=119164, sorted_swp 28309 [ms]
N=119164, avx2 13436 [ms]
N=119164, sorted_z 10420 [ms]
N=119164, sorted_z_avx2 14540 [ms]
N=119164, avx512 9527 [ms]
N=119164, avx512_loopopt 8430 [ms]
N=119164, avx512_loopopt_swp 6839 [ms]
N=119164, avx512_transpose 9759 [ms]

$ ./soa.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 19619 [ms]
N=119164, sorted 10259 [ms]
N=119164, sorted_swp 27508 [ms]
N=119164, avx2 21029 [ms]
N=119164, avx512 9628 [ms]
N=119164, avx512_loopopt 7916 [ms]
N=119164, avx512_loopopt_swp 7270 [ms]
# システムC
# icpc (ICC) 18.0.1 20171018
# Intel(R) Xeon(R) Gold 6130 CPU @ 2.10GHz
$ ./aos.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 5256 [ms]
N=119164, sorted 2919 [ms]
N=119164, sorted_swp 3819 [ms]
N=119164, avx2 1890 [ms]
N=119164, sorted_z 2852 [ms]
N=119164, sorted_z_avx2 1954 [ms]
N=119164, avx512 2698 [ms]
N=119164, avx512_loopopt 2301 [ms]
N=119164, avx512_loopopt_swp 1949 [ms]
N=119164, avx512_transpose 2459 [ms]

$ ./soa.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 4235 [ms]
N=119164, sorted 3014 [ms]
N=119164, sorted_swp 3676 [ms]
N=119164, avx2 3124 [ms]
N=119164, avx512 2424 [ms]
N=119164, avx512_loopopt 2231 [ms]
N=119164, avx512_loopopt_swp 2048 [ms]

 KNLでは逆数近似を使うともう少し(5%)程度早くなるらしいのだが、もう疲れた。 これでx86/x64最適化勉強会の準備はできたことにする。

 無駄なコードが入ったままorigin/masterにpushしてしまった。

      // --8<--

 SWP(ソフトウェアパイプライニング)をする時に、ループの半個ずらしをするのだが、そのキリトリセンを残したままだった。 面白いからそのままにしておくか。

 AVX2を使ったSIMD化を真面目にやり始めたのが2016年の9月らしい。そこからだいぶSIMD化したので、 ループの4倍/8倍展開や、SWPにだいぶ慣れた。特にSWPは、ステップごとにデバッグしながら進める手順を確立したのだが、 これ、忘れる前にどこかに残しておきたいなぁ。後、gatherで集めた8要素8個のzmmベクトルを4x4で上下で転置して _mm512_extractf64x4_pdでバラして_mm256_store_pdで格納する奴は、結局遅かったけどそれなりに頑張ったので、 いずれQiitaかどこかに葬っておきたい。

 よし、次だ次!論文と数独やるぞ。


2月5日(月)

 if文によるcontinueをゼロクリアに変えたら、インテルコンパイラの最新版でもgather/scatterが出るようになった。

void
force_sorted(void) {
  const int pn = particle_number;
  const int * __restrict sorted_list2 = sorted_list;
  for (int i = 0; i < pn; i++) {
    const double qx_key = q[i][X];
    const double qy_key = q[i][Y];
    const double qz_key = q[i][Z];
    const int np = number_of_partners[i];
    double pfx = 0;
    double pfy = 0;
    double pfz = 0;
    const int kp = pointer[i];
    for (int k = 0; k < np; k++) {
      const int j = sorted_list2[kp + k];
      double dx = q[j][X] - qx_key;
      double dy = q[j][Y] - qy_key;
      double dz = q[j][Z] - qz_key;
      double r2 = (dx * dx + dy * dy + dz * dz);
      //if (r2 > CL2) continue;  // ←これをやめて
      double r6 = r2 * r2 * r2;
      double df = ((24.0 * r6 - 48.0) / (r6 * r6 * r2)) * dt;
      if (r2 > CL2) df=0.0;  // ←ここでゼロクリア
      pfx += df * dx;
      pfy += df * dy;
      pfz += df * dz;
      p[j][X] -= df * dx;
      p[j][Y] -= df * dy;
      p[j][Z] -= df * dz;
    }
    p[i][X] += pfx;
    p[i][Y] += pfy;
    p[i][Z] += pfz;
  }
}

 NKLではgather/scatter両方使っており、Skylakeではgatherしか使っていない。 インテルコンパイラの-xと-axがよくわからない。

-axcode Tells the compiler to generate multiple, feature-specific auto-dispatch code paths for Intel(R) processors if there is a performance benefit.
-xcode Tells the compiler which processor features it may target, including which instruction sets and optimizations it may generate.

 同じバージョンのコンパイラに、同じコードを食わせているのに、Skylakeでは-xMIC-AVX512ではgatherしか使わないが、 -axMIC-AVX512ではgather/scatter/vpconflictdを出す。KNLでは-xMIC-AVX512だけで全部使う。 また、Skylakeでは-axCORE-AVX512ではgatherを使わないが、-xCORE-AVX512では使う。

 それはともかく、if文によるcontinueをゼロクリアに変えた場合、AoSでもSoAでもgather/scatterを使うようになった。

# Oakforest-PACS 
# if分岐をゼロクリアに修正後
# icpc -O3 -std=c++11 -w2 -w3 -diag-disable:remark -restrict -axMIC-AVX512
$ ./aos.out
Number of pairs: 7839886
N=119164, pair 37102 [ms]
N=119164, sorted 10647 [ms]
N=119164, sorted_swp 28328 [ms]
N=119164, avx2 13783 [ms]
N=119164, sorted_z 29226 [ms]
N=119164, sorted_z_avx2 18557 [ms]

$ ./soa.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 36015 [ms]
N=119164, sorted 10506 [ms]
N=119164, sorted_swp 27637 [ms]
N=119164, avx2 21222 [ms]
# システムC
# icpc (ICC) 18.0.1 20171018
# icpc -O3 -std=c++11 -w2 -w3 -diag-disable:remark -xCORE-AVX512 -qopt-zmm-usage=high 
$ ./aos.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 5240 [ms]
N=119164, sorted 2914 [ms]
N=119164, sorted_swp 3821 [ms]
N=119164, avx2 1866 [ms]
N=119164, sorted_z 4013 [ms]
N=119164, sorted_z_avx2 1914 [ms]

$ ./soa.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 4214 [ms]
N=119164, sorted 3009 [ms]
N=119164, sorted_swp 3668 [ms]
N=119164, avx2 3147 [ms]

 ちなみにSoA版に-axMIC-AVX512つけて無理やりgather/scatter/conflictd使わせてみたけど、 gatherしか使わない版に比べて
AoS: 2914 [ms] → 5816 [ms]
SoA: 3009 [ms] → 5794 [ms]

と遅くなる。なんでだ?

 もともとコンパイラがvpconflictdを使うのは、i粒子と相互作用するj粒子に重複が無いことをコンパイラが知らないからで、 vpconflictdで重複を確認したらシリアルコードに飛び、そうでない場合はscatterしてしまう、というコードを吐いているから。 しかし、Skylakeではgatherしか使わないので衝突検出は不要になる。KNLではSkylakeに比べてメモリのランダムな書き戻しが 重いため、vpconflictdを使ってなるべくscatterを使った方が得だけど、Skylakeではgatherは使うけどvpconflictd使うくらいなら そのまま書き戻した方が早い、とか?

 AVX-512でSoAはSIMD化できてAoSはSIMD化できなかった記憶があるんだけど、再現できない。 SoAだけif分岐のゼロクリアをしてて、それを見逃してて勘違いしたとかかな・・・。 うーん、でもそこはちゃんと調べたつもりなんだけど。もっと古いコンパイラを使ってたとか?

 いずれにせよコンパイラ依存性は強い。intel/2016.4.258 の方がintel/2018.1.163より有意に早いコードを吐く。

sorted AoS (2016.4.258)  9621 +- 28[ms]
sorted AoS (2018.1.163) 10626 +-  2[ms]
sorted SoA (2016.4.258)  8867 +- 34[ms]
sorted SoA (2018.1.163) 10306 +- 49[ms]

 ただし、この頃のインテルコンパイラはバグが多かったと聞くので、バグを修正していったら性能も少し控えめになった可能性がある。

 というわけでとりあえずAVX-512のhand-simdizationのためにループを8倍展開して一個ずつデバッグして、とかやってるんだけど、 途中経過がこんな感じになるんすよ。

      const int j_1 = sorted_list2[kp + k];
      const int j_2 = sorted_list2[kp + k + 1];
      const int j_3 = sorted_list2[kp + k + 2];
      const int j_4 = sorted_list2[kp + k + 3];
      const int j_5 = sorted_list2[kp + k + 4];
      const int j_6 = sorted_list2[kp + k + 5];
      const int j_7 = sorted_list2[kp + k + 6];
      const int j_8 = sorted_list2[kp + k + 7];
      const auto vindex = _mm256_lddqu_si256((const __m256i*)(&sorted_list[kp + k]));
      const v8df vqxj = _mm512_i32gather_pd(vindex, &(q[X][0]), 8);
      const v8df vqyj = _mm512_i32gather_pd(vindex, &(q[Y][0]), 8);
      const v8df vqzj = _mm512_i32gather_pd(vindex, &(q[Z][0]), 8);
      const v8df vdx = vqxj - vqxi;
      const v8df vdy = vqyj - vqyi;
      const v8df vdz = vqzj - vqzi;
      const v8df vr2 = vdx*vdx + vdy*vdy + vdz*vdz;
      double dx_1 = q[X][j_1] - qx_key;
      double dy_1 = q[Y][j_1] - qy_key;
      double dz_1 = q[Z][j_1] - qz_key;
      double dx_2 = q[X][j_2] - qx_key;
      double dy_2 = q[Y][j_2] - qy_key;
      double dz_2 = q[Z][j_2] - qz_key;
      double dx_3 = q[X][j_3] - qx_key;
      double dy_3 = q[Y][j_3] - qy_key;
      double dz_3 = q[Z][j_3] - qz_key;
      double dx_4 = q[X][j_4] - qx_key;
      double dy_4 = q[Y][j_4] - qy_key;
      double dz_4 = q[Z][j_4] - qz_key;
      double dx_5 = q[X][j_5] - qx_key;
      double dy_5 = q[Y][j_5] - qy_key;
      double dz_5 = q[Z][j_5] - qz_key;
      double dx_6 = q[X][j_6] - qx_key;
      double dy_6 = q[Y][j_6] - qy_key;
      double dz_6 = q[Z][j_6] - qz_key;
      double dx_7 = q[X][j_7] - qx_key;
      double dy_7 = q[Y][j_7] - qy_key;
      double dz_7 = q[Z][j_7] - qz_key;
      double dx_8 = q[X][j_8] - qx_key;
      double dy_8 = q[Y][j_8] - qy_key;
      double dz_8 = q[Z][j_8] - qz_key;
      double r2_1 = (dx_1 * dx_1 + dy_1 * dy_1 + dz_1 * dz_1);
      double r2_2 = (dx_2 * dx_2 + dy_2 * dy_2 + dz_2 * dz_2);
      double r2_3 = (dx_3 * dx_3 + dy_3 * dy_3 + dz_3 * dz_3);
      double r2_4 = (dx_4 * dx_4 + dy_4 * dy_4 + dz_4 * dz_4);
      double r2_5 = (dx_5 * dx_5 + dy_5 * dy_5 + dz_5 * dz_5);
      double r2_6 = (dx_6 * dx_6 + dy_6 * dy_6 + dz_6 * dz_6);
      double r2_7 = (dx_7 * dx_7 + dy_7 * dy_7 + dz_7 * dz_7);
      double r2_8 = (dx_8 * dx_8 + dy_8 * dy_8 + dz_8 * dz_8);
      if(i==5){
        putsval(r2);
        puts(vr2);
        exit(1);
      }

 これでループを8倍展開して、8個のペアの距離の平方までできたところ。 最終的にデバッグ用のシリアル計算の部分は消えて、だいぶすっきりする予定なんだけど、作業途中はソースが10倍近く膨れ上がる。 単純作業なんだし、こういうのもうちょっとなんとかならんかな?という気はする。

 どうでもいいけど、デバッグ用にいろいろ出力する関数の名前、putsにしたのはあまりよくなかったかな。

 まずはSoAで単純8倍馬鹿SIMDルーチン作った。まずKNL。

$ ./soa.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 19609 [ms]
N=119164, sorted 10259 [ms]
N=119164, sorted_swp 27504 [ms]
N=119164, avx2 21025 [ms]
N=119164, avx512 9606 [ms]

 vpconflictdを使わなくて良い分、ちょっとだけ早くなる。

 しかし、SkylakeではAoS+AVX2の方が早い。

$ ./soa.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 4207 [ms]
N=119164, sorted 3011 [ms]
N=119164, sorted_swp 3654 [ms]
N=119164, avx2 3128 [ms]
N=119164, avx512 2395 [ms]

 SoAではAVX2実装(3128ms)よりAVX512実装(2395ms)の方が早いが、AoSのAVX2実装(1866)の方が早い。 あと試すのはループ端数処理最適化と8要素パック実装か。今日中にループ端数処理くらいまでいけるかな。

 ループ端数処理できた。基本的には中川さんのコードそのままだが、オリジナルがカットオフ処理と端数処理両方でblendしているのに対して、二つのマスクを_mm512_kandでANDとってから一度だけblendしたら、ほんのすこしだけ(2287ms→2248ms)早くなった。 (追記)オリジナルでも最速実装では_mm512_kand使ってた。いやん。

# システムC
$ ./soa.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 4211 [ms]
N=119164, sorted 3013 [ms]
N=119164, sorted_swp 3677 [ms]
N=119164, avx2 3111 [ms]
N=119164, avx512 2413 [ms]
N=119164, avx512_loopopt 2243 [ms]
# # Oakforest-PACS 
$ ./soa.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 19632 [ms]
N=119164, sorted 10260 [ms]
N=119164, sorted_swp 27549 [ms]
N=119164, avx2 21047 [ms]
N=119164, avx512 9635 [ms]
N=119164, avx512_loopopt 7958 [ms]

 ループ端数処理最適化の効果はKNLの方が大きいっぽい。Skylakeでは性能向上が7%程度だが、KNLでは20%以上早くなる。

 科研費の支払請求書(F-2)提出した。

 よし、なんとか今日中にAoSのAVX-512化+ループ端数最適化までやった。これでもうコードはいじらなくて良いはず。

# # Oakforest-PACS 
$ ./aos.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 37123 [ms]
N=119164, sorted 10525 [ms]
N=119164, sorted_swp 28478 [ms]
N=119164, avx2 13432 [ms]
N=119164, sorted_z 10418 [ms]
N=119164, sorted_z_avx2 14534 [ms]
N=119164, avx512 9523 [ms]
N=119164, avx512_loopopt 8432 [ms]

 ありゃ、SoAに負けてるぞ。swpしないとダメかな・・・。


2月4日(日)

 金曜日の「コンパイラのバージョン上げたら単純なループが早くなった」件、最初はAVX-512のgather/scatterを 使ったせいかとおもったが、よく考えたらその前にメモリバンド幅ネックになるはずなのでおかしい。 で、よく調べたらコンパイラが関数ポインタのインライン展開をして、ループの外側を交換していたからだった。

 時間計測ルーチンとしてこんなのを書いた。

void
measure(void(*pfunc)(), const char *name, int particle_number) {
  const auto s = std::chrono::system_clock::now();
  const int LOOP = 100;
  for (int i = 0; i < LOOP; i++) {
    pfunc();
  }
  const auto e = std::chrono::system_clock::now();
  const auto elapsed = std::chrono::duration_cast<std::chrono::milliseconds>(e - s).count();
  fprintf(stderr, "N=%d, %s %lld [ms]\n", particle_number, name, elapsed);
}

 関数のポインタpfuncをもらって、それを100回よんで実行時間を計測するものだが、インテルコンパイラは18.0.0.1から 関数ポインタのインライン展開ができるようになった。観測対象となるforce_pairはこんな関数になっている。

void
force_pair(void) {
  for (int k = 0; k < number_of_pairs; k++) {
    const int i = i_particles[k];
    const int j = j_particles[k];
    double dx = q[j][X] - q[i][X];
    double dy = q[j][Y] - q[i][Y];
    double dz = q[j][Z] - q[i][Z];
    double r2 = (dx * dx + dy * dy + dz * dz);
    if (r2 > CL2) continue;
    double r6 = r2 * r2 * r2;
    double df = ((24.0 * r6 - 48.0) / (r6 * r6 * r2)) * dt;
    p[i][X] += df * dx;
    p[i][Y] += df * dy;
    p[i][Z] += df * dz;
    p[j][X] -= df * dx;
    p[j][Y] -= df * dy;
    p[j][Z] -= df * dz;
  }
}

 純粋にペアの数(number_of_pairs)だけループを回して力を計算する関数。これをmeasureから100回呼んだわけだが、 インテルコンパイラはmeasureの中にforce_pairをインライン展開し、100回まわる観測のループとnumber_of_pairs回まわる ループを交換してしまった。アセンブリの該当部分はこうなっている。

        movslq    number_of_pairs(%rip), %rcx 
        xorl      %eax, %eax   
        testq     %rcx, %rcx 
        jle       ..B1.117 

(計算ルーチン)

..B1.114: 
        incb      %r8b 
        cmpb      $100, %r8b                                    #100回まわる観測ループ
        jb        ..B1.112 
..B1.115:
        incq      %rax
        cmpq      %rcx, %rax                                    #number_of_pairs回まわる計算ループ
        jb        ..B1.111 

確かにループが入れ替わってしまっている。本質的に同じ計算だが、i粒子とj粒子が完全にキャッシュに乗るので メモリアクセスのレイテンシが全く消えてしまう。インテルコンパイラはわりとこういう大胆な「ループ交換」をやるのでびっくりする。 GCC/clangもやるのかもしれないが、少なくとも僕が試したコードでは見たことがない。

 というわけで、ダミーのグローバル変数を入れてループ交換を阻害する。

// インテルコンパイラのループ交換最適化阻害のためのダミー変数
int sum = 0;
void
measure(void(*pfunc)(), const char *name, int particle_number) {
  const auto s = std::chrono::system_clock::now();
  const int LOOP = 100;
  for (int i = 0; i < LOOP; i++) {
    sum++; // ループ交換阻害
    pfunc();
  }
  const auto e = std::chrono::system_clock::now();
  const auto elapsed = std::chrono::duration_cast<std::chrono::milliseconds>(e - s).count();
  fprintf(stderr, "N=%d, %s %lld [ms]\n", particle_number, name, elapsed);
}
..B1.109:
        movslq    number_of_pairs(%rip), %rcx
(計算ルーチン)
..B1.114:
        incq      %rdi
        cmpq      %rcx, %rdi #number_of_pairs回まわる計算ループ
        jb        ..B1.112
..B1.116:
        incb      %al
        cmpb      $100, %al  #100回まわる観測ループ
        jb        ..B1.110 

 ちゃんと所望のループ構造になった。やれやれだ。もっときれいな方法があるのかもしれないけど、 ディレクティブとか入れるのは好きじゃないので、とりあえずこれでいいことにする。

 ずっと「AoS→SoA変換しただけでコンパイラがgather/scatterを吐いてすこぶる早くなった」という事象を再現できずに苦しんでいたんだけど、インテルコンパイラのバージョンの問題と判明。2016.4.258でSIMD化できてたコードがintel/2018.1.163でできなくなっている。 あと、生配列のままより、一度 int * __restrictedにポインタをコピーした方が早くなるのはちょっとどうかなぁ。

int sorted_list[MAX_PAIRS];

void
force_sorted(void) {
  const int pn = particle_number;
  for (int i = 0; i < pn; i++) {
    const int np = number_of_partners[i];
    const int kp = pointer[i];
    for (int k = 0; k < np; k++) {
      const int j = sorted_list[kp + k];
      //なんか計算
    }
  }
}

 みたいなコードがあるんだけど、これを

int sorted_list[MAX_PAIRS];

void
force_sorted(void) {
  const int pn = particle_number;
  const int * __restrict sorted_list2 = sorted_list;  //ここでrestrict修飾されたポインタに受ける
  for (int i = 0; i < pn; i++) {
    const int np = number_of_partners[i];
    const int kp = pointer[i];
    for (int k = 0; k < np; k++) {
      const int j = sorted_list2[kp + k]; //restrict修飾されたポインタを使う
      //なんか計算
    }
  }
}

 とするだけで1割以上早くなる。コンパイルオプションに-restrictつけてもダメ。 生配列の方が早いかと思ってそっちを使ってたのだが、少なくともこのケースでは__restrict修飾されたポインタの方が早い。

 改めてちゃんとした計測結果と、コンパイラのバージョン依存性を残しておこう。

# Oakforest-PACS 
# c5f29e916d1eaa3f227dfa0566ad63ebf103d509
# icpc -O3 -std=c++11 -w2 -w3 -diag-disable:remark -restrict -axMIC-AVX512 
# module load intel/2016.4.258
$ ./aos.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 37268 [ms]
N=119164, sorted 10717 [ms]
N=119164, sorted_swp 28808 [ms]
N=119164, avx2 13619 [ms]
N=119164, sorted_z 15550 [ms]
N=119164, sorted_z_avx2 19346 [ms]

$ ./soa.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 36101 [ms]
N=119164, sorted 9934 [ms]
N=119164, sorted_swp 27776 [ms]
N=119164, avx2 21296 [ms]

# intel/2018.1.163
$ ./aos.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 37090 [ms]
N=119164, sorted 31544 [ms]
N=119164, sorted_swp 28339 [ms]
N=119164, avx2 13778 [ms]
N=119164, sorted_z 29145 [ms]
N=119164, sorted_z_avx2 18515 [ms]

$ ./soa.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 36065 [ms]
N=119164, sorted 30902 [ms]
N=119164, sorted_swp 27500 [ms]
N=119164, avx2 20959 [ms]

 強いバージョン依存性があるけど、よくわからんから可視化してみるか。

 以前の経験ではAoSは古いバージョンでもpairをSIMD化できなかったのだが、いまやってみたらできてるな・・・。 なんだろう、オプションのせいだろうか。もうちゃんと詰めて追う元気が無い・・・。

# システムC
# icpc version 18.0.1 (gcc version 4.8.5 compatibility)
$ ./aos.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 5237 [ms]
N=119164, sorted 4635 [ms]
N=119164, sorted_swp 3806 [ms]
N=119164, avx2 1870 [ms]
N=119164, sorted_z 4015 [ms]
N=119164, sorted_z_avx2 1910 [ms]

$ ./soa.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 5106 [ms]
N=119164, sorted 4462 [ms]
N=119164, sorted_swp 3657 [ms]
N=119164, avx2 3149 [ms]

# module load intel/17.0.6 
$ ./aos.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 5259 [ms]
N=119164, sorted 4590 [ms]
N=119164, sorted_swp 3842 [ms]
N=119164, avx2 1864 [ms]
N=119164, sorted_z 4011 [ms]
N=119164, sorted_z_avx2 1904 [ms]

$ ./soa.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 5134 [ms]
N=119164, sorted 4511 [ms]
N=119164, sorted_swp 3734 [ms]
N=119164, avx2 3142 [ms]

 Skylakeではバージョン依存性は無いな。

# システムB
# icpc (ICC) 16.0.4 20160811
./aos.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 8348 [ms]
N=119164, sorted 7257 [ms]
N=119164, sorted_swp 6124 [ms]
N=119164, avx2 3766 [ms]
N=119164, sorted_z 6488 [ms]
N=119164, sorted_z_avx2 3468 [ms]

$ ./soa.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 8175 [ms]
N=119164, sorted 7170 [ms]
N=119164, sorted_swp 5752 [ms]
N=119164, avx2 5758 [ms]
# iMac
# Intel Core i5 @ 3.3 GHz
$ ./aos.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 6233 [ms]
N=119164, sorted 4561 [ms]
N=119164, sorted_swp 4183 [ms]
N=119164, avx2 1803 [ms]
N=119164, sorted_z 3869 [ms]
N=119164, sorted_z_avx2 1782 [ms]

$ ./soa.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 5974 [ms]
N=119164, sorted 4432 [ms]
N=119164, sorted_swp 3706 [ms]
N=119164, avx2 3357 [ms]

 後の自分のためにまとめておこう。

 ・・・ふぅ。手でgather/scatterを出すコードを試すのは後回しにしよう。


2月3日(土)

 NOP


2月2日(金)

 今日も雪。

 今日のSPAM。

Dear Dr.Huwaida M E Hassaneen ,

 誰だよ(エジプトのカイロ大学の人らしい)。

 他にも「multiphysics@list.uva.nl」というメーリングリストに無断で入れられており、 そこに論文投稿のお誘いSPAMが来たのだが、それに「なんで俺のところにこんなメールが来てるんだ?」 とか「fucking spammers...」とかそういう返事が次から次へと投稿され、それが全部 メーリングリストに流れるという地獄絵図。なんだかな。

 lj_simdを大幅にリファクタリング。 AoSとSoAで似たコードがかなりあったのをだいぶまとめた。また、sys/time.hからchronoを使うように修正。 intrinと書いていたのをavx2に修正(AVX-512版を追加する準備)。 念のため、システムBとCの結果を貼っておくか。

# 物性研システムB
$ ./aos.out
Number of pairs: 7839886
N=119164, pair 8188 [ms]
N=119164, sorted 7046 [ms]
N=119164, sorted_swp 5848 [ms]
N=119164, avx2 3508 [ms]
N=119164, sorted_z 6064 [ms]
N=119164, sorted_z_avx2 3274 [ms]

$ ./soa.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 9513 [ms]
N=119164, sorted 7183 [ms]
N=119164, sorted_swp 5629 [ms]
N=119164, avx2 5701 [ms]
# 物性研システムC
# g++ -O3 -std=c++11 -march=native
# g++ (GCC) 7.2.0
$ ./aos.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 6593 [ms]
N=119164, sorted 4953 [ms]
N=119164, sorted_swp 4518 [ms]
N=119164, avx2 1951 [ms]
N=119164, sorted_z 4202 [ms]
N=119164, sorted_z_avx2 1929 [ms]

$ ./soa.out 
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 6727 [ms]
N=119164, sorted 4861 [ms]
N=119164, sorted_swp 4092 [ms]
N=119164, avx2 3408 [ms]
# Oakforest-PACS 
# icpc -O3 -std=c++11 -w2 -w3 -diag-disable:remark -axMIC-AVX512 
# (最適化により正しくない評価になっている、追記参照)
$ ./aos.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 9783 [ms]
N=119164, sorted 31525 [ms]
N=119164, sorted_swp 28347 [ms]
N=119164, avx2 13414 [ms]
N=119164, sorted_z 29086 [ms]
N=119164, sorted_z_avx2 18522 [ms]

$ ./soa.out
Make pairlist.
Number of pairs: 7839886
N=119164, pair 9752 [ms]
N=119164, sorted 31004 [ms]
N=119164, sorted_swp 27573 [ms]
N=119164, avx2 20903 [ms]

  なにこれ?KNLの結果が以前より4倍近く早くなってるんですけど? コンパイラを古くして試して見る。

# Oakforest-PACS 
$ module purge
$ module load intel/2016.4.258 
$ make
$ ./aos.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 37332 [ms]
N=119164, sorted 12550 [ms]
N=119164, sorted_swp 28808 [ms]
N=119164, avx2 13646 [ms]
N=119164, sorted_z 15637 [ms]
N=119164, sorted_z_avx2 19381 [ms]

$ ./soa.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 36110 [ms]
N=119164, sorted 12219 [ms]
N=119164, sorted_swp 27783 [ms]
N=119164, avx2 21296 [ms]

 最新版(intel/2018.1.163)で再ビルドして試したけど間違いない。 コンパイラの性格が劇的に変わっている。

 古→新について、「pair (ペアごとに力を計算)」が劇的に早くなり、「sorted (i粒子をキーにソート)」が劇的に遅くなっている。 ちゃんとasm見てないけど、速度的にgather/scatterの使用の有無と見て間違いないだろう。 っていうかコンパイラの進歩すげーな。あのアクセスパターンでgather/scatter吐けるのか。 でもなんでi粒子でソートした奴が遅くなったんだろ?

 AVX2を直接書いちゃってる奴の速度が変わらないのは、そりゃそうだよね、という感じ。 sorted_swpは手でソフトウェアパイプライニング(一個ずらし)をかけたもので、これは 旧版では遅くなる(ループボディが複雑になって最適化がかからないせいと思われる)。 新版ではsortedが遅いので、相対的にSWPが効いて少し早くなる。 しかし、最も単純な書き方が最も早いってきついな。これから試すAVX-512組み込み版でひっくり返せなかったら 「人類はコンパイラに負けました」ってことで、それはそれで喜ばしいことではあるけれど・・・。

 念のため物性研システムC(Skylake)でインテルコンパイラで再度試す。

# 物性研システムC
# icpc -O3 -std=c++11 -w2 -w3 -diag-disable:remark -xCORE-AVX512 -qopt-zmm-usage=high
# Intel(R) Xeon(R) Gold 6130 CPU @ 2.10GHz
# icpc (ICC) 18.0.1 20171018
# (最適化により正しくない評価になっている、追記参照)
$ ./aos.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 1908 [ms]
N=119164, sorted 4620 [ms]
N=119164, sorted_swp 3802 [ms]
N=119164, avx2 1861 [ms]
N=119164, sorted_z 4005 [ms]
N=119164, sorted_z_avx2 1908 [ms]

$ ./soa.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 1876 [ms]
N=119164, sorted 4586 [ms]
N=119164, sorted_swp 3660 [ms]
N=119164, avx2 3145 [ms]

 icpc+単純ループがGCC+AVX2の最速実装を越えてきますね・・・。ふむ。

 Haswellでも確認。

# 物性研システムB
# 
# icpc -O3 -std=c++11 -w2 -w3 -diag-disable:remark -xHOST
# icpc (ICC) 16.0.4 20160811
$ ./aos.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 8299 [ms]
N=119164, sorted 7144 [ms]
N=119164, sorted_swp 5917 [ms]
N=119164, avx2 3555 [ms]
N=119164, sorted_z 6186 [ms]
N=119164, sorted_z_avx2 3315 [ms]

$ ./soa.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 8056 [ms]
N=119164, sorted 7126 [ms]
N=119164, sorted_swp 5650 [ms]
N=119164, avx2 5716 [ms]

# icpc (ICC) 18.0.1 20171018
$ (最適化により正しくない評価になっている、追記参照)
$ ./aos.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 2538 [ms]
N=119164, sorted 7076 [ms]
N=119164, sorted_swp 5680 [ms]
N=119164, avx2 5834 [ms]

$ ./soa.out
A pair-file is found. I use it.
Number of pairs: 7839886
N=119164, pair 2533 [ms]
N=119164, sorted 7058 [ms]
N=119164, sorted_swp 5710 [ms]
N=119164, avx2 5859 [ms]

なんだ?インテルコンパイラのバージョンを上げたらpairが劇的に早くなったぞ。 こっちはHaswellだからgather/scatterが無いのに怪しい。コンパイラが何かを見抜いてインチキしているっぽいな・・・。 (追記)インテルコンパイラが関数ポインタをインライン展開することでおかしな最適化が入ったのが原因。


2月1日(木)

 一月が終わってしまった。

 なんだかさ、多分人生って100年ないわけじゃん。研究人生も、あと数十年。例えば65定年として25年ちょっと。 25*12=300ということで、あと300ヶ月。先月はその1/300だったわけだけど、「貴重な人生の1/300を費やして何をしたのか?」 と聞かれてしまうと「うーん」と思ってしまう。最近は特に時間が取れずに苦しい。

 生ゴミ処理機を買う時に柏市の助成金を利用したのだが、その状況調査が来たので答えて出す。 子供がいる家庭ではどうしても生ゴミが出やすいので非常に助かる。 これ以外にも柏市には子育て用の様々な助成金や制度があり、子育てがし易い環境と感じる。

 昨日のSPAM、議論しているサイトを見つけた。Flaky Academic Journals。掲載されているメールも一語一句同じだった。 まぁ要するに怪しい雑誌だったということ。

 「査読システムは本当に必要なのか?」と真正面から問われてしまい、答えに窮する。僕は必要だと思っているが、それをうまく言語化できていない。ほとんどの論文がpublish前にarXivで読まれている時代、「どの雑誌に通ったか」というのは研究/研究者の評価にしか使われていない気がする。また、こうして論文投稿しませんか/本を書きませんかSPAMが大量に来るということは、それが商売になるからであり、 裏を返せば形だけの査読をして、査読論文の数を稼ぎたい研究者が多い、ということになる。 その辺考え出すと「うーん」となってしまう。


トップページへ戻る