SIMDの基礎と関数の実装例

配信講義　計算科学技術特論A（2023）（2023年4月13日～7月27日・オンライン）
第11回 6月29日 SIMDの基礎と関数の実装例
https://www.r-ccs.riken.jp/outreach/schools/20230413-1/

Transcript

1. SIMDの基礎と関数の実装例 サイボウズ・ラボ 光成滋生 計算科学技術特論A（2023） 2023/6/29 1 / 58
2. 概要 目次 1. x64アセンブラとSIMDの基本 2. AVX-512 3. exp/logの実装例 2 / 58
3. 背景 優れたライブラリやツールは沢山ある わざわざ自分で開発する必要はない 行列計算における高速アルゴリズム１より引用 が、それらのコードがアセンブリ言語（以下ASM）で記述されることが多いのも事実 将来、もしかしたら自分が開発に関わるかもしれない ASM/SIMD（特にx64のAVX-512）を紹介する 3 / 58
4. 自作アセンブラ xbyak : x64用の動的なコードを生成するC++ライブラリ Intel oneAPI(oneDNN)のCPUエンジンで利用されている oneDNNはPyTorchやTensorFlowなどのIntel版で使われている チップレットになった「第4世代Xeon SP」、性能向上の鍵はAMXと4つのアクセラレータ pc.watch.impress.co.jpより引用 スーパーコンピュータ富岳用のxbyak_aarch64もある s_xbyak : Pythonでxbyakライクな記述ができるx64用の静的なアセンブラ 今回はこちらがメイン。詳細はXbyakライクなx64用静的ASM生成ツールs_xbyak 4 / 58
5. SIMD 一つの命令で複数のデータを同時に処理する AVX-512は512ビットのZMMレジスタを32個持つ zmm0, zmm1, ..., zmm31 整数型 8ビットを64個, 16ビットを32個, 32ビットを16個, 64ビットを8個扱える 浮動小数点数型 32ビット(float)を16個, 64ビット(double)を8個扱える Sapphire RapidsやGranite Rapidsは一部16ビットを32個として扱う命令もサポート Granite Rapids/Sierra ForestのAMX/AVXの新命令 基本の形 op z, x, y # opは命令, z, x, yはレジスタやメモリを指す z ← op(x, y) # xとyのopの結果をzに代入する 5 / 58
6. float型の加算 vaddps(z, x, y) z, x, yはZMMレジスタ xのi番目のfloat型の要素を とする for psはpacked single precision(float)の意味 x ... + + + + ... + y ... = = = = ... = z ... 6 / 58
7. 添え字の種類 浮動小数点数型 v+命令+{pd, ps} ps : float pd : double (double precision) 整数型 vp + 命令 + {b, w, d, q} b : byte (8ビット整数) w : word (16ビット整数) d : dword (32ビット整数) q : qword (64ビット整数) 7 / 58
8. SIMDを使う 配列同士の加算関数 void addc(float *z, const float *x, const float *y, size_t n) { for (size_t i = 0; i < n; i++) { z[i] = x[i] + y[i]; } } 簡単にするための仮定 nは正で16の倍数とする x[0..n], y[0..n], z[0..n]は互いにオーバーラップしていない restrictとする x, y, zは64バイト境界にアラインされている（ことが望ましい） 8 / 58
9. SIMD(AVX-512)での実装例 s_xbyakによる実装例 with FuncProc('add_avx512'): # add_avx512という関数を作る with StackFrame(4, vNum=1, vType=T_ZMM) as sf: # 関数の引数は4個. ZMMレジスタを1個使う pz = sf.p[0] # 1番目の引数はzへのポインタ（を示すレジスタ） px = sf.p[1] # 2番目の引数はxへのポインタ py = sf.p[2] # 3番目の引数はyへのポインタ n = sf.p[3] # 4番目の引数はn lp = Label() # ループ用ラベル L(lp) # ループ用ラベルの定義 vmovups(zmm0, ptr(px)) # zmm0 = *px vaddps(zmm0, zmm0, ptr(py)) # zmm0 += *py as float vmovups(ptr(pz), zmm0) # *pz = zmm0 add(px, 64) # px += 64 add(py, 64) # py += 64 add(pz, 64) # pz += 64 sub(n, 16) # n -= 16 jnz(lp) # n != 0 ならlpに戻って繰り返す 9 / 58
10. x64 CPUの基礎 汎用レジスタ 64ビット整数 : rax, rbx, rcx, rdx, rsi, rdi, rbp, r8-r15 32ビット整数 : eax, ebx, ecx, edx, esi, edi, ebp, r8d-r15d 64ビットレジスタの下位32ビットを使う |63 32|31 16|15 8|7 0| | rax | | | eax | | | | ax | | | | ah | al | スタックレジスタ : rsp（ローカル変数など保持するスタック領域を指す） その他のレジスタ フラグレジスタEFLAGS（ZF, CFなどのフラグの集合 - 後述） SIMDでない浮動小数点数型もあるが略 10 / 58
11. 基本命令 Intel ASMの基本の形 op x, y # x ← op(x, y) xやyはレジスタやメモリを指す xとyのopの操作結果がxに代入される 殆どの命令はEFLAGSを実行ごとに更新する mov x, y # x ← y mov(x, ptr(y)) は x = *(uint64_t*)y; の意味（xが64ビットの場合） xがSIMDレジスタのときはvmovupsなどを使う 算術演算(add, sub, ...)や論理演算(and, xor, or, ...) add x, y # x ← x + y and x, y # x ← x & y 11 / 58
12. 条件分岐 比較 : cmp x, y x - yの結果をEFLAGSにセットするが演算結果はxに反映しない ZF(Zero Flag) : x == yなら1 CF(Carry Flag) : x < yなら1 (unsigned) cf. sub x, yは演算結果がxに反映される（x ← x - y） テスト : test x, y x & yの結果をEFLAGSにセットするが演算結果はxに反映されない ZF : x & y == 0なら1 and x, yは同じ値をEFLAGSにセットしてxに反映される（x ← x & y） その他の命令 多くの演算はEFLAGSを変更するが主に使われるのは上記2種類 12 / 58
13. フラグに基づく処理 無条件ジャンプ jmp label : いつでもlabelにジャンプ 条件ジャンプ j<なんとか> label : <なんとか>が成り立てばlabelにジャンプ je : ZF=1（結果が等しい, 結果がゼロ）ならジャンプ jg : 符号ありで大きければジャンプ ja : 符号無しで大きければジャンプ 他にjge, jne, jleなどなど 条件mov cmovz x, y : ZF=1ならmov x, y など条件ジャンプと同じくいろいろある 13 / 58
14. x64 Windows上のMASM用出力 -m masmオプションの結果(*.asmファイル) _text$x segment align(64) execute # 64バイト境界アライン可能な実行可能なセグメント指定 add_avx512 proc export # 関数の始まり @L1: vmovups zmm0, zmmword ptr [rdx] # 関数の2番目の引数xはrdx vaddps zmm0, zmm0, zmmword ptr [r8] # 3番目の引数yはr8 vmovups zmmword ptr [rcx], zmm0 # 1番目の引数zはrcx add rdx, 64 add r8, 64 add rcx, 64 sub r9, 16 # 4番目の引数nはr9 jnz @L1 ret add_avx512 endp # 関数の終わり _text$x ends # セグメントの終わり end # ファイルの終わり Visual Studio付属のml64.exeでアセンブルしてobjファイルを作る C（や他の言語）から呼び出して使う 14 / 58
15. 呼び出し規約 CからASMを呼び出す場合 引数に割り当てられるレジスタや保存すべきレジスタがある 関数の引数（整数型・ポインタ型4個までの例） Windows : rcx, rdx, r8, r9 Linux : rdi, rsi, rdx, rcx 呼び出された関数内で保存すべきレジスタ Windows : rbx, rbp, rdi, rsi, r12-r15, rsp Linux : rbx, rbp, r12-r15, rsp 整数型の返り値はraxに入る 詳細 Windows x64 ソフトウェアの呼び出し規約 x64 Linux System V Application Binary Interface 15 / 58
16. x64 Linux/macOS上のGAS用出力 -m gasオプションの結果(*.Sファイル) ... 略 #endif .global PRE(add_avx512) # PRE(add_avx512)のシンボルを外部から参照可能にする PRE(add_avx512): # PRE(add_avx512)のシンボルを定義する(PREはmacOSでアンダースコアを付与する) .L1: vmovups (%rsi), %zmm0 # 2番目の引数はrsi vaddps (%rdx), %zmm0, %zmm0 # 3番目の引数はrdx vmovups %zmm0, (%rdi) # 1番目の引数はrdi add $64, %rsi add $64, %rdx add $64, %rdi sub $16, %rcx jnz .L1 ret SIZE(add_avx512) # Linux用GASで関数のサイズを設定する LinuxのGNU assemblerは引数の順序が逆になる（細かな差異は構文差異参照） Windows MASMとLinux GASの違いを吸収している 16 / 58
17. レイテンシとスループット CPUにおけるレイテンシ 命令が発行されてからその実行が完了するまでの時間（クロックサイクル : 以下clkと記す） レイテンシが4なら実行されてから4clk後に結果を利用可能 IntelによるCPUのスループット 同じ命令を続けて発行するときに待つclk 本来のスループットは単位時間あたりに処理できる量だが Intelのスループットは逆数(Reciprocal)スループットと言われることがある（文脈判断）17 / 58
18. スループットの例 FMA(Fused Multiply -Add) を計算する命令 レイテンシ4 スループット0.5で命令同士が依存関係がない場合 スループット0.5なので1clkごとに2命令発行できる 依存関係が無いなら次の命令は1clk後に発行できる 結果を利用できるのは4clk後になる cf. 第1回 計算科学技術特論A（2023） 18 / 58
19. x64のFMA 命名規則 FMAは入力が3個なので組み合わせが名前に入っている vfmadd132ps( , , ) vfmadd213ps( , , ) vfmadd231ps( , , ) ±や加算・減算によるバリエーション vfmsub(xy-z), vfnmadd(-xy+z), vfnmsub(-xy-z)などがある 19 / 58
20. FMAのスループットを確認する実験 実験コード ループ内に互いに依存関係の無いFMAを複数個並べて速度評価する # func{n}(int c); // n個のFMAをc回ループする with FuncProc(f'func{n}'): with StackFrame(1, vNum=n+1, vType=T_ZMM) as sf: c = sf.p[0] lp = Label() for i in range(n): vxorps(Zmm(i), Zmm(i), Zmm(i)) # ZMMレジスタクリア align(32) L(lp) for i in range(n): vfmadd231ps(Zmm(i), Zmm(i), Zmm(i)) # n個のFMAを発行 sub(c, 1) jnz(lp) # c回ループする コード全体はfma Xeon Platinum 8280 Turbo Boost offで測定 20 / 58
21. 生成コード n=1 @L1: vfmadd231ps zmm0, zmm0, zmm0 sub rcx, 1 jnz @L1 n=4 @L4: vfmadd231ps zmm0, zmm0, zmm0 vfmadd231ps zmm1, zmm1, zmm1 vfmadd231ps zmm2, zmm2, zmm2 vfmadd231ps zmm3, zmm3, zmm3 sub rcx, 1 jnz @L4 21 / 58
22. 実験結果 nを増やしたときの1ループあたりの処理時間 n 1 2 3 4 5 6 7 8 9 10 clk 4.1 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.5 5.0 n=8まで互いに独立なFMAを発行しても時間が変化しない（レイテンシ4で処理できている） ループアンロールが重要 FMAの連続発行の流れ 22 / 58
23. FMAの利用例 多項式の評価 次多項式 の における値を計算する のときの求め方 と変形する として後ろから求める 1. 2. 3. 4. それぞれの計算にFMAを利用する 同じ多項式に対する複数の による評価はループアンロールしてレイテンシを隠す 23 / 58
24. 係数の保持方法 係数を全てレジスタに保持する 余裕があるならレジスタに置いておくのがよい ループアンロールするとレジスタが足りなくなる可能性がある 係数をメモリに保持する メモリにSIMD分並べておいて読む メモリに1個分だけおいてSIMDレジスタに個数分だけ展開する（ブロードキャスト） 24 / 58
25. ブロードキャストの方法 ブロードキャスト無しのFMA # coeff # |c0|c0|c0|...|c0|c1|c1|... vfmadd231ps(zmm0, zmm1, ptr(coeff)) # coeffから16個分の同じfloatを読み込みFMA ブロードキャスト命令vbroadcastss # coeff # |c0|c1|c2|... vbroadcastss(zmm2, ptr(coeff)) # coeffから1個分のfloatを読み込みzmm2に展開する vfmadd231ps(zmm0, zmm1, zmm2) FMA+ptr_b（ブロードキャストフラグ） vfmadd231ps(zmm0, zmm1, ptr_b(coeff)) # coeffから1個分のfloatを読み込み展開してFMA ptr_bはxbyakの記法. MASMはdword bcstでGASは{1to16}など 25 / 58
26. 実験結果 ループ展開したときのそれぞれの処理時間 N 1 2 3 4 5 6 7 8 9 全部レ ジスタ 4.13 4.02 4.00 4.00 4.00 4.00 4.00 4.00 4.57 broad cast 4.13 4.01 4.00 4.00 4.00 4.06 4.22 4.74 5.25 ptr_b 4.14 4.01 4.00 4.00 4.00 4.00 4.06 4.38 4.97 全部レジスタならN=8がよい メモリから読みながらならptr_bを使うのがよい 26 / 58
27. exp(x)の近似計算 アルゴリズム を2ベキの形にする （ ） を整数 と小数部分 ( )に分割する とすると の計算 が整数なのでビット演算で処理できる ( )の計算 floatの精度 1b-23 = に近くなるまでローラン展開する 27 / 58
28. Sollya 近似計算のためのソフトウェア https://www.sollya.org/ 数学の公式のローラン展開よりもよりよい値を提示してくれる 次数の見積もり >guessdegree(2^x,[-0.5,0.5],1b-23); [5;5] 多項式近似 >fpminimax(2^x,[|1,2,3,4,5|],[|D...|],[-0.5,0.5],1); 1 + x * (0.69314697759916432673321651236619800329208374023437 + x * (0.24022242085378028852993281816452508792281150817871 + x * (5.5507337432541360711102385039339424110949039459229e-2 + x * (9.6715126395259202324306002651610469911247491836548e-3 + x * 1.326472719636653634089906717008489067666232585907e-3)))) 28 / 58
29. アルゴリズムのまとめ exp(x)の近似計算 入力 : x 出力 : exp(x) 事前準備 : c[i] : 多項式の係数 1. 2. ここで は四捨五入関数。 3. 4. 5. 6. return . 29 / 58
30. の計算 通常の方法 floatのビット表現 floatのビット表現 符号s 指数部e 仮数部f ビット長 1 8 23 (n+127)<<23 が に対応する AVX-512の場合 n : float型の整数 w : float値 vscalefps(n, w) で となる 30 / 58
31. 四捨五入 cvtps2dq float→int型変換 : AVX-512で丸めモードを指定できる 丸めモード(ctrl) round-to-nearest-even : 一番近い整数(0.5は偶数に)に丸める round-toward-zero : 0に近い方向に丸める round-down : -∞方向に丸める round-up : +∞方向に丸める 結果はint型なのでfloat型への変換が必要 vrndscaleps 結果はfloat型 31 / 58
32. vreduceps 先に小数部を求める命令 レイテンシとスループット 命令 レイテンシ スループット vrndscaleps 8 1 vreduceps 4 0.5-1 先に整数部分 先に小数部分 n ← vrndscaleps(x) a ← vreduceps(x) a ← x - n n ← x - a 先に小数部分を求める方が速くて都合がよい（nよりaが先に必要） 32 / 58
33. 実装例 コア部分 入力 : v0 出力 : v0 = exp(v0) 事前準備 : log2_e, expCoeff[5]のレジスタに値を設定しておく # v0 = exp(v0) vmulps(v0, v0, self.log2_e) vreduceps(v1, v0, 0) # a = x - n vsubps(v0, v0, v1) # n = x - a = round(x) vmovaps(v2, self.expCoeff[5]) for i in range(4, -1, -1): vfmadd213ps(v2, v1, self.expCoeff[i]) vscalefps(v0, v2, v0) # v2 * 2^n 33 / 58
34. s_xbyakによるループアンロール Pythonによるマクロ的な使い方 例 vfmadd231ps zmm0, zmm0, zmm5 vfmadd231ps zmm1, zmm1, zmm5 vfmadd231ps zmm2, zmm2, zmm5 vfmadd231ps zmm3, zmm3, zmm5 は v = [zmm0, zmm1, zmm2, zmm3] unroll(vfmadd231ps)(v, v, zmm5) と（たとえば）書けると便利 34 / 58
35. メモリ参照のループアンロール マクロの続き vfaddps(zmm0, zmm0, ptr(rax)) vfaddps(zmm1, zmm1, ptr(rax+64)) vfaddps(zmm2, zmm2, ptr(rax+128)) なら v = [zmm0, zmm1, zmm2] unroll(vfaddps)(v, v, ptr(rax)) 自動的にオフセットを付与してほしい ptr_bを使っていたらオフセットはつけないでほしい? 35 / 58
36. アンロール機能の例 Pythonの自由度の高さを利用する def Unroll(n, op, *args, addrOffset=None): xs = list(args) for i in range(n): ys = [] for e in xs: if isinstance(e, list): # 引数が配列ならi番目を利用する ys.append(e[i]) elif isinstance(e, Address): # 引数がアドレスなら if addrOffset == None: if e.broadcast: addrOffset = 0 # broadcastモードならオフセット0（このあたりは好みで） else: addrOffset = SIMD_BYTE # そうでないときはSIMDのサイズずらす(addrOffsetで細かい制御はできる) ys.append(e + addrOffset*i) else: ys.append(e) op(*ys) 36 / 58
37. unrollの利用例 ヘルパー関数 アンロール回数Nを引数として、オペランドopを取り引数に対してアンロールする関数を返す def genUnrollFunc(n): """ return a function takes op and outputs a function that takes *args and outputs n unrolled op """ def fn(op, addrOffset=None): def gn(*args): Unroll(n, op, *args, addrOffset=addrOffset) return gn return fn # v0 = [zmm0, zmm1, ...], v1 = [zmm4, zmm5, ...], ... un = genUnrollFunc(n) # アンロール回数を指定する un(vmulps)(v0, v0, self.log2_e) un(vreduceps)(v1, v0, 0) # a = x - n un(vsubps)(v0, v0, v1) # n = x - a = round(x) 37 / 58
38. アンロールされたコード 複数のレジスタ v0=[zmm0, zmm1, zmm2] un(vmulps)(z0, z0, log2_e) → vmulps(zmm0, zmm0, log2_e) vmulps(zmm1, zmm1, log2_e) vmulps(zmm2, zmm2, log2_e) メモリ参照 v0=[zmm0, zmm1, zmm2] v1=[zmm3, zmm4, zmm5] un(vfmadd231ps)(v0, v1, ptr(rax)) → vfmadd231ps(zmm0, zmm3, ptr(rax)) vfmadd231ps(zmm1, zmm4, ptr(rax+64)) vfmadd231ps(zmm2, zmm5, ptr(rax+128)) ptrの代わりにptr_bを使うとオフセットは0のまま 38 / 58
39. exp(x)のベンチマーク アンロール回数Nを変えて測定 N 1 2 3 4 5 6 7 8 全てレ ジスタ 17.91 15.89 14.14 13.85 13.68 13.08 13.03 13.78 定数を ptr_bで 読む 18.06 16.21 14.82 14.37 14.54 14.61 14.66 16.19 全てレジスタに置く場合はN=7がよかった あまりアンロールすると命令キャッシュを圧迫してしまう 一部の定数をメモリに置いた場合はN=4がよかった 全体のコードはfmathのgen_fmath.py 39 / 58
40. log(x)の近似計算 アルゴリズムの方針 を ( は整数 )の形に分解する に対する を求めることに専念する ローラン展開 では収束が遅い 40 / 58
41. 収束を速くするために範囲を狭める 天下りだが を の近似値とする と置くと なので は0に近い と変形し、対数を取ると もし を計算できるなら よりもより0に近い について を求めればよい の計算 は[1, 2)の範囲なので上位9ビット(符号+指数部)は固定される 仮数部の上位ビットでテーブルルックアップをすれば が求まる 41 / 58
42. テーブルルックアップ の仮数部の上位Lビットを利用する の下位 ビットをマスクして近似値 を得る をインデックスとする tbl1には , tbl2には の値を入れておく 42 / 58
43. テーブルの作成 Pythonコード self.logTbl1 = [] self.logTbl2 = [] self.L = 4 LN = 1 << self.L for i in range(LN): u = (127 << 23) | ((i*2+1) << (23 - self.L - 1)) v = 1 / uint2float(u) v = uint2float(float2uint(v)) # (X) 32ビットのfloat値に強制変換 self.logTbl1.append(v) self.logTbl2.append(math.log(v)) uint2floatはuint32_tの値に対応するfloat値を得る関数 float2uintはfloat値に対応するuint32_tの値を得る関数 (X)の意味 PythonのfloatはCにおけるdoubleなので、そのまま使うとインデックスが（少し）ずれる 43 / 58
44. AVX-512におけるテーブルルックアップ vgather命令 vgatherdd(out|k, ptr(rax + idx*4)) # out[i] = *(rax + idx[i]) 複数のインデックスをSIMDレジスタに入れて1命令で複数のテーブル参照を行う 便利だが遅い（レイテンシ20以上+メモリアクセス） vpermps vpermps(x, y, z) # x[i] = y[z[i]] floatが16(=1<<L)個入るZMMレジスタをテーブルとして利用（レイテンシ3） vpermi2ps 512ビットZMMレジスタ2個をテーブルとして利用（レイテンシ3） 今回は（遅くなったので）使わない 44 / 58
45. 仮数部と指数部の処理 従来の方法 floatのビット表現を用いて仮数部と指数部を取り出す e = (x >> 23) - 127, m = x & mask(23) AVX-512の方法 vgetexpss # floatの指数部を取り出す vgetmantss # floatの仮数部を取り出す ビット演算や余計な定数を保持する必要がなく高速 45 / 58
46. log(x)の実装例 ローラン展開をする前処理 入力 : v0 出力 : v1 一時変数 : v2, v3 定数 : L = 4, tbl1, tbl2, one=1.0 un(vgetexpps)(v1, v0) # n un(vgetmantps)(v0, v0, 0) # a un(vpsrad)(v2, v0, 23 - self.L) # d = v0 >> (23-L) un(vpermps)(v3, v2, self.tbl1) # b = tbl1[d] = 1/a un(vfmsub213ps)(v0, v3, self.one) # c = a * b - 1 un(vpermps)(v3, v2, self.tbl2) # log_b = tbl2[d] = log(1/a) un(vfmsub132ps)(v1, v3, ptr_b(rip+self.LOG2)) # z = n * log2 - log_b ここからv0に対して log(v0) をローラン展開で求める 46 / 58
47. log(1+x)のローラン展開 4次で打ち切った数式の場合 4次で打ち切って係数を調整した場合 // 相対誤差 float relErr(float x, float y) { return x == 0 ? 0 : std::fabs(x - y) / std::fabs(x); } 区間[1-1/32,1+1/32]でxを1e-6ずつ変更したときのstd::logとの誤差 式 数式 補正 最大誤差 3.54e-7 2.44e-7 平均誤差 4.91e-8 3.97e-8 47 / 58
48. x~1のときの計算方法 精度劣化 なので のときは直接ローラン展開を利用した方がよい テーブル参照の処理が入ると精度が落ちる 区間内でxを1e-6ずつ増やしながらstd::log(x)との相対誤差を求める 区間 [0.99,1.01] [2,3] 最大誤差 2.80e-2 1.19e-7 平均誤差 2.61e-5 2.38e-8 区間[0.99,1.01]は区間[2,3]に比べて精度劣化が大きい 48 / 58
49. 分岐が必要 概念的なコード if math.abs(x-1) < B: # Bは適当な境界値 return log(1+x)のローラン展開 else: return テーブルルックアップの後ローラン展開 SIMDでの分岐 要素ごとに分岐するのは困難 両方の計算をしておき、結果をマスクで選択する v1 = compute_if_true(x) v2 = compute_if_false(x) cond = math.abs(x-1) < B # x~1の判定 mask = 0xffffffff if cond else 0 return (v1 & mask) | (v2 & ~mask) 49 / 58
50. AVX-512における分岐 64ビットのマスクレジスタk1, ..., k7 iビット目が 1ならZMMレジスタのi番目の操作を有効にする 0なら値を変更しない faddps(x|k, y, z) y ... z ... k ... x ... 50 / 58
51. 分岐コードの挿入位置 計算の流れ 1. 入力 : v0 = x 2. テーブルルックアップ処理後 3. v0 ← c = ab - 1 4. v1 ← n * log2 - log_b 5. （挿入）ここで なら v0 ← x-1, v1 ← 0 とする 6. v2 ← log(1+v0) # ローラン展開 7. v0 ← v2 + v1 # 出力 log(x) 51 / 58
52. 挿入するコード |x-1|< B を確認する keepX ← v0 un(vsubps)(v2, keepX, self.one) # x-1 un(vandps)(v3, v2, ptr_b(rip+self.C_0x7fffffff)) # |x-1| un(vcmpltps)(vk, v3, ptr_b(rip+self.BOUNDARY)) un(vmovaps)(zipOr(v0, vk), v2) # c = v0 = x-1 un(vxorps)(zipOr(v1, vk), v1, v1) # z = 0 floatの絶対値は最上位ビットを落とす（0x7fffffffとand） vcmpltpsで比較してマスクレジスタvkに結果を格納 zipOrはv0とvkの要素をorするヘルパー関数 def zipOr(v, k): r = [] for i in range(len(v)): r.append(v[i]|k[i]) return r 52 / 58
53. 挿入結果 挿入前と挿入後で相対誤差を確認する 区間 [0.99,1.01] [2,3] 最大誤差 2.80e-2 → 1.19e-7 1.19e-7 平均誤差 2.61e-5 → 3.02e-8 2.38e-8 [0.99,1.01]の精度が改善された 速度 - 挿入前 挿入後 clk 14.96 19.13 アンロール回数n = 4のとき約4clk遅くなった 相対誤差ではなく絶対誤差で十分なとき（logの和など）は補正しない選択枝もある 53 / 58
54. メモリアクセスにおけるマスクレジスタ ループ回数が32の倍数でない場合 マスクレジスタを使って余分なアクセスをさせない マスクレジスタを使わない場合は境界をまたがないように注意 54 / 58
55. ゼロ化フラグ に対応するレジスタを0クリアする vmovups(x|k|T_z,ptr(rax)) T_z はxbyakの表記, MASM/GASでは {z} マージ（デフォルト）はレジスタの依存関係が発生するので遅くなる可能性 ゼロ化フラグをつけると実行前のxの値に依存しないので高速に処理できる vmovups(x|k|T_z, ptr(m)) x ... m ... k ... x ... 55 / 58
56. ループ処理の概略(1/2) メインと端数処理に分割 メイン処理 : 16N（Nはアンロール回数）ずつ要素を処理する mov(rcx, n) # n は配列の要素数 jmp(check1L) align(32) L(lpUnrollL) un(vmovups)(v0, ptr(src)) # 16N個分のfloatを読み込む add(src, 64*unrollN) # 64Nバイト分srcアドレスを進める func(unrollN, args) # 関数本体を呼び出す un(vmovups)(ptr(dst), v0) # 16N個結果を書き込む add(dst, 64*unrollN) # 64Nバイト分dstアドレスを進める sub(n, 16*unrollN) # 16Nカウンタを減らす L(check1L) cmp(n, 16*unrollN) jae(lpUnrollL) # n >= 16*Nの間ループ 理想はsrcを64バイトアライメントする処理を入れる（vmovupsでも十分速いけど） 注意 : マスクレジスタを使うと遅くなるので不要なときは使わない 56 / 58
57. ループ処理の概略(2/2) 端数処理 0 < ecx < 16 and_(ecx, 15) jz(exitL) mov(eax, 1) # eax = 1 shl(eax, cl) # eax = 1 << n sub(eax, 1) # eax = (1 << n) - 1 kmovd(k1, eax) # k1 = eax vmovups(zm0|k1|T_z, ptr(src)) func(1, args) # 関数本体を呼び出す vmovups(ptr(dst)|k1, zm0) L(exitL) (1<<n)-1 でn個の1が立つビットマスクを作る kmovd(k1, eax) でマスクレジスタを設定 zmm0|k1|T_z を使ってメモリの読み込みをする 57 / 58
58. まとめ AVX-512 レジスタ数が多いのでループアンロールが有効 スループットとレイテンシを考慮する ブロードキャスト命令の代わりにブロードキャストフラグでメモリアクセスを減らす マスクレジスタを使って分岐を避ける AVX-512専用命令を活用する 整数部や小数部を取り出す : vrndscaleps, vreduceps 指数部や仮数部を取り出す : vgetexpss, vgetmantss 小さいテーブルルックアップ : vpermps, vpermi2psなど 58 / 58