SIMDの基礎と関数の実装例

Slide 1

Slide 1 text

SIMDの基礎と関数の実装例サイボウズ・ラボ光成滋生計算科学技術特論A（2023） 2023/6/29 1 / 58

Slide 2

Slide 2 text

概要目次 1. x64アセンブラとSIMDの基本 2. AVX-512 3. exp/logの実装例 2 / 58

Slide 3

Slide 3 text

背景優れたライブラリやツールは沢山あるわざわざ自分で開発する必要はない行列計算における高速アルゴリズム１より引用が、それらのコードがアセンブリ言語（以下ASM）で記述されることが多いのも事実将来、もしかしたら自分が開発に関わるかもしれない ASM/SIMD（特にx64のAVX-512）を紹介する 3 / 58

Slide 4

Slide 4 text

自作アセンブラ 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

Slide 5

Slide 5 text

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

Slide 6

Slide 6 text

float型の加算 vaddps(z, x, y) z, x, yはZMMレジスタ xのi番目のfloat型の要素をとする for psはpacked single precision(float)の意味 x ... + + + + ... + y ... = = = = ... = z ... 6 / 58

Slide 7

Slide 7 text

添え字の種類浮動小数点数型 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

Slide 8

Slide 8 text

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

Slide 9

Slide 9 text

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

Slide 10

Slide 10 text

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

Slide 11

Slide 11 text

基本命令 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

Slide 12

Slide 12 text

条件分岐比較 : 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

Slide 13

Slide 13 text

フラグに基づく処理無条件ジャンプ 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

Slide 14

Slide 14 text

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

Slide 15

Slide 15 text

呼び出し規約 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

Slide 16

Slide 16 text

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

Slide 17

Slide 17 text

レイテンシとスループット CPUにおけるレイテンシ命令が発行されてからその実行が完了するまでの時間（クロックサイクル : 以下clkと記す）レイテンシが4なら実行されてから4clk後に結果を利用可能 IntelによるCPUのスループット同じ命令を続けて発行するときに待つclk 本来のスループットは単位時間あたりに処理できる量だが Intelのスループットは逆数(Reciprocal)スループットと言われることがある（文脈判断）17 / 58

Slide 18

Slide 18 text

スループットの例 FMA(Fused Multiply -Add) を計算する命令レイテンシ4 スループット0.5で命令同士が依存関係がない場合スループット0.5なので1clkごとに2命令発行できる依存関係が無いなら次の命令は1clk後に発行できる結果を利用できるのは4clk後になる cf. 第1回計算科学技術特論A（2023） 18 / 58

Slide 19

Slide 19 text

x64のFMA 命名規則 FMAは入力が3個なので組み合わせが名前に入っている vfmadd132ps( , , ) vfmadd213ps( , , ) vfmadd231ps( , , ) ±や加算・減算によるバリエーション vfmsub(xy-z), vfnmadd(-xy+z), vfnmsub(-xy-z)などがある 19 / 58

Slide 20

Slide 20 text

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

Slide 21

Slide 21 text

生成コード 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

Slide 22

Slide 22 text

実験結果 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

Slide 23

Slide 23 text

FMAの利用例多項式の評価次多項式のにおける値を計算するのときの求め方と変形するとして後ろから求める 1. 2. 3. 4. それぞれの計算にFMAを利用する同じ多項式に対する複数のによる評価はループアンロールしてレイテンシを隠す 23 / 58

Slide 24

Slide 24 text

係数の保持方法係数を全てレジスタに保持する余裕があるならレジスタに置いておくのがよいループアンロールするとレジスタが足りなくなる可能性がある係数をメモリに保持するメモリにSIMD分並べておいて読むメモリに1個分だけおいてSIMDレジスタに個数分だけ展開する（ブロードキャスト） 24 / 58

Slide 25

Slide 25 text

ブロードキャストの方法ブロードキャスト無しの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

Slide 26

Slide 26 text

実験結果ループ展開したときのそれぞれの処理時間 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

Slide 27

Slide 27 text

exp(x)の近似計算アルゴリズムを2ベキの形にする（）を整数と小数部分 ( )に分割するとするとの計算が整数なのでビット演算で処理できる ( )の計算 floatの精度 1b-23 = に近くなるまでローラン展開する 27 / 58

Slide 28

Slide 28 text

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

Slide 29

Slide 29 text

アルゴリズムのまとめ exp(x)の近似計算入力 : x 出力 : exp(x) 事前準備 : c[i] : 多項式の係数 1. 2. ここでは四捨五入関数。 3. 4. 5. 6. return . 29 / 58

Slide 30

Slide 30 text

の計算通常の方法 floatのビット表現 floatのビット表現符号s 指数部e 仮数部f ビット長 1 8 23 (n+127)<<23 がに対応する AVX-512の場合 n : float型の整数 w : float値 vscalefps(n, w) でとなる 30 / 58

Slide 31

Slide 31 text

四捨五入 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

Slide 32

Slide 32 text

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

Slide 33

Slide 33 text

実装例コア部分入力 : 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

Slide 34

Slide 34 text

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

Slide 35

Slide 35 text

メモリ参照のループアンロールマクロの続き 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

Slide 36

Slide 36 text

アンロール機能の例 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

Slide 37

Slide 37 text

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

Slide 38

Slide 38 text

アンロールされたコード複数のレジスタ 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

Slide 39

Slide 39 text

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

Slide 40

Slide 40 text

log(x)の近似計算アルゴリズムの方針を ( は整数 )の形に分解するに対するを求めることに専念するローラン展開では収束が遅い 40 / 58

Slide 41

Slide 41 text

収束を速くするために範囲を狭める天下りだがをの近似値とすると置くとなのでは0に近いと変形し、対数を取るともしを計算できるならよりもより0に近いについてを求めればよいの計算は[1, 2)の範囲なので上位9ビット(符号+指数部)は固定される仮数部の上位ビットでテーブルルックアップをすればが求まる 41 / 58

Slide 42

Slide 42 text

テーブルルックアップの仮数部の上位Lビットを利用するの下位ビットをマスクして近似値を得るをインデックスとする tbl1には , tbl2にはの値を入れておく 42 / 58

Slide 43

Slide 43 text

テーブルの作成 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

Slide 44

Slide 44 text

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<

Slide 45

Slide 45 text

仮数部と指数部の処理従来の方法 floatのビット表現を用いて仮数部と指数部を取り出す e = (x >> 23) - 127, m = x & mask(23) AVX-512の方法 vgetexpss # floatの指数部を取り出す vgetmantss # floatの仮数部を取り出すビット演算や余計な定数を保持する必要がなく高速 45 / 58

Slide 46

Slide 46 text

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

Slide 47

Slide 47 text

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

Slide 48

Slide 48 text

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

Slide 49

Slide 49 text

分岐が必要概念的なコード 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

Slide 50

Slide 50 text

AVX-512における分岐 64ビットのマスクレジスタk1, ..., k7 iビット目が 1ならZMMレジスタのi番目の操作を有効にする 0なら値を変更しない faddps(x|k, y, z) y ... z ... k ... x ... 50 / 58

Slide 51

Slide 51 text

分岐コードの挿入位置計算の流れ 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

Slide 52

Slide 52 text

挿入するコード |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

Slide 53

Slide 53 text

挿入結果挿入前と挿入後で相対誤差を確認する区間 [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

Slide 54

Slide 54 text

メモリアクセスにおけるマスクレジスタループ回数が32の倍数でない場合マスクレジスタを使って余分なアクセスをさせないマスクレジスタを使わない場合は境界をまたがないように注意 54 / 58

Slide 55

Slide 55 text

ゼロ化フラグに対応するレジスタを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

Slide 56

Slide 56 text

ループ処理の概略(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

Slide 57

Slide 57 text

ループ処理の概略(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<

Slide 58

Slide 58 text

まとめ AVX-512 レジスタ数が多いのでループアンロールが有効スループットとレイテンシを考慮するブロードキャスト命令の代わりにブロードキャストフラグでメモリアクセスを減らすマスクレジスタを使って分岐を避ける AVX-512専用命令を活用する整数部や小数部を取り出す : vrndscaleps, vreduceps 指数部や仮数部を取り出す : vgetexpss, vgetmantss 小さいテーブルルックアップ : vpermps, vpermi2psなど 58 / 58