LJの力計算のSIMD化ステップ・バイ・ステップ その6

はじめに

そうとうがんばってSIMD化して、もう速くならないだろうと思ってたらkohnakagawaさんから「もう少し高速化しました」というプルリクが来て慌てた話。

コードは https://github.com/kaityo256/lj_simdstep においてある。

どこを直したか

一度に256ビット取ってくる命令を使うため、配列をパディングして、(x,y,z,0)という並びにしていた。それを4要素のベクトルとしてロードして、相対座標ベクトルを作るのだが、ループを4倍展開しているので、

\[(dx_1, dy_1, dz_1, 0)\\ (dx_2, dy_2, dz_2, 0)\\ (dx_3, dy_3, dz_3, 0)\\ (dx_4, dy_4, dz_4, 0)\]

という四本のベクトルができる。これを、それぞれ自乗してからシャッフルして足すことで、相対距離を作っていた。

\[(dr^2_1, dr^2_1, dr^2_1,0) \\ (dr^2_2, dr^2_2, dr^2_2,0) \\ (dr^2_3, dr^2_3, dr^2_3,0) \\ (dr^2_4, dr^2_4, dr^2_4,0) \\\]

で、これを転置して、

\[v_r = (dr^2_1, dr^2_2, dr^2_3, dr^4)\]

を作ってから力の計算を4並列でやっていた。しかし、どうせ最後に転置するなら、最初に転置した方が計算量が減る。つまり、

\[(dx_1, dy_1, dz_1, 0)\\ (dx_2, dy_2, dz_2, 0)\\ (dx_3, dy_3, dz_3, 0)\\ (dx_4, dy_4, dz_4, 0)\]

を転置して、

\[v_x = (dx_1, dx_2, dx_3, dx_4)\\ v_y = (dy_1, dy_2, dy_3, dy_4)\\ v_z = (dz_1, dz_2, dz_3, dy_4)\\ v_0 = (0, 0, 0, 0)\]

を得たら、あとは

\[v_r = v_x^2+v_y^2+v_z^2\]
でできるじゃん、というもの。言われてみればそれはそうだ。なんで気が付かなかったんだろ○ ￣ ＿

結果

というわけで上記を実装したものをこれまでと同条件で測定。

計算手法 実行時間[sec] 備考
pair 8.201647 ペアごとにループをまわしたもの
pair_swp 7.094673 上記をSWPしたもの
pair_swp_intrin 4.282396 上記を4倍展開してSIMD化したもの
sorted 7.040513 i粒子でソートしたもの
sorted_intrin 3.889211 i粒子でソートしたものをSIMD化したもの
sorted_swp 5.806741 i粒子でソートしたものをSWPしたもの
sorted_swp_intrin 3.272044 i粒子でソートしたものをSWPしたものをSIMD化したもの
sorted_swp_intrin_mat_transp 2.980598 上記の転置順序を変えたもの

10%弱高速化されましたな。

コードは以下に。

https://github.com/kaityo256/lj_simdstep/tree/master/step6

まとめ

もともと、4倍展開する前にSIMD化してたものを4倍展開してから馬鹿SIMD化したので、いろいろ計算した後、最後に転置する形しか考えつかなかった。でも、プルリクもらってよく考えたら、先に転置しちゃった方が計算量少ないから速い。うーん、まだまだ修行が足りませんな。

この様子じゃ、データ構造とか工夫したり、ループのほどきかた工夫したりしたら、まだまだ高速化の余地がありそうですね。

オレはようやくのぼりはじめたばかりだからな このはてしなく遠いSIMD坂をよ・・・

こんどこそおしまい？