SIMDの基礎と関数の実装例

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

生成コード
n=1
@L1:
vfmadd231ps zmm0, zmm0, zmm0
sub rcx, 1
jnz @L1
n=4
@L4:
sub rcx, 1
jnz @L4
21 / 58

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

s_xbyakによるループアンロール
Pythonによるマクロ的な使い方
例
は
v = [zmm0, zmm1, zmm2, zmm3]
unroll(vfmadd231ps)(v, v, zmm5)
と（たとえば）書けると便利
34 / 58

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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<vpermi2ps
512ビットZMMレジスタ2個をテーブルとして利用（レイテンシ3）
今回は（遅くなったので）使わない
44 / 58

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

ループ処理の概略(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<kmovd(k1, eax) でマスクレジスタを設定
zmm0|k1|T_z を使ってメモリの読み込みをする
57 / 58

View Slide

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

View Slide

SIMDの基礎と関数の実装例

SIMDの基礎と関数の実装例

herumi

More Decks by herumi

Other Decks in Programming

Featured

Transcript