C++
AVX
AVX2

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

More than 1 year has passed since last update.

はじめに

そうとうがんばって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坂をよ・・・

こんどこそおしまい?