AVX-512はFMA（fused multiply-add、 融合積和）命令を実装しています。FMA命令を数式にしたものを以下に示します。
a × b + k
式が表すように、積と和が合体した演算です。性能が向上すだけでなく、丸め（浮動小数点数で正確に表せない値を近い値にする）が1回しか発生しないため、積と和で実装したときに比べ精度が向上します。
FMAに関しては、精度、速度、コンパイラなのサポート状況について多くの議論がなされており、FMA自体の議論については、そちらへ譲り、本書ではAVX-512でサポートしているベクトル対応したFMAについて紹介します。
64ビット浮動小数点
呼び出し側
以降に、呼び出し側のC++コードを示します。一次元配列から指定された値のインデックスと、その数を取得します。
common.hをincludeしていますがSIMDで記述したアセンブリ言語の関数と等価な機能を持つC++言語で記述した共通関数のヘッダです。最下部にソースコードを示します。条件は関数名に含まれます。
#include "..\common.h" // TEMPLATES #define T double extern "C" { void afmaaddpd (const T*, const T*, T*, const T, const size_t); void afmasubpd (const T*, const T*, T*, const T, const size_t); void afnmaaddpd (const T*, const T*, T*, const T, const size_t); void afnmasubpd (const T*, const T*, T*, const T, const size_t); void afmaddsubpd(const T*, const T*, T*, const T, const size_t); void afmsubaddpd(const T*, const T*, T*, const T, const size_t); } typedef void (*Dfunc)(const T*, const T*, T*, const T, const size_t); extern "C" Dfunc afunc[] = { afmaaddpd, afmasubpd, afnmaaddpd, afnmasubpd, afmaddsubpd, afmsubaddpd }; Dfunc cfunc[] = { cfmaadd, cfmasub, cfnmaadd, cfnmasub, cfmaddsub, cfmsubadd }; // main int main(void) { const size_t ArrLen = 4096; static_assert(ArrLen % 16 == 0, "number of elements must be an integral multiple of 16."); T a[ArrLen], b[ArrLen], c[ArrLen], v[ArrLen], k = 11.5f; init(a, b, ArrLen); for (int i = 0; i < sizeof(cfunc) / sizeof(*cfunc);i++) { cout << "---[" << i << "]--- "; cfunc[i](a, b, c, k, ArrLen); afunc[i](a, b, v, k, ArrLen); verify(c, v, ArrLen); } return 0; }
簡単な処理ですので説明は省きます。
呼び出され側
本関数へ渡される内容とレジスターの対応を表で示します。
| レジスタ | 内容 |
|---|---|
| rcx | 入力配列 a の先頭アドレス |
| rdx | 入力配列 b の先頭アドレス |
| r8 | 出力配列 c の先頭アドレス |
| zmm2 | 最下位要素に定数 k |
| r9 | 入力配列の要素数 |
64ビット浮動小数点へ対応ライブラリを、マクロを使って開発したものを示します。FMAに限らず、いくつかの操作を行います。このような関数をマクロで記述したものを示します。
mymacro macro MNAME, MINST public MNAME align 16 MNAME proc mov r9, qword ptr [rsp+40] ; r9 = length vbroadcastsd zmm2, xmm3 ; zmm2 = Scalar Value xor rax, rax ; clear index loop_f: vmovupd zmm0, zmmword ptr [rcx+rax*8] ; load a[] vmovupd zmm1, zmmword ptr [rdx+rax*8] ; load b[] MINST zmm0, zmm1, zmm2 ; zmm0 = ([-]a*b) <op> k, ... vmovupd zmmword ptr [r8+rax*8], zmm0 ; stote c[] add rax, 8 cmp rax, r9 jb short loop_f ret MNAME endp endm ;------------------------------------------------------------------- ; code _TEXT segment mymacro afmaaddpd, vfmadd213pd mymacro afmasubpd, vfmsub213pd mymacro afnmaaddpd, vfnmadd213pd mymacro afnmasubpd, vfnmsub213pd mymacro afmaddsubpd, vfmaddsub213pd mymacro afmsubaddpd, vfmsubadd213pd _TEXT ends end
mymacro マクロを定義し、これを使ってafmaaddpd、afmasubpd、afnmaaddpd、afnmasubpd、afmaddsubpd、afmsubaddpd関数を作ります。各関数の機能を表で示します。
| 関数名 | 機能 |
|---|---|
| afmaaddpd | a*b + k |
| afmasubpd | a*b - k |
| afnmaaddpd | -a*b + k |
| afnmasubpd | -a*b - k |
| afmaddsubpd | 偶数要素: a*b + c、奇数要素: a*b - k |
| afmsubaddpd | 偶数要素 a*b - c、奇数要素: a*b + k |
分かり易くするため、afmaaddpdが展開されたソースリストを示します。
public afmaaddpd align 16 afmaaddpd proc mov r9, qword ptr [rsp+40] ; r9 = length vbroadcastsd zmm2, xmm3 ; zmm2 = Scalar Value xor rax, rax ; clear index loop_f: vmovupd zmm0, zmmword ptr [rcx+rax*8] ; load a[] vmovupd zmm1, zmmword ptr [rdx+rax*8] ; load b[] vfmadd213pd zmm0, zmm1, zmm2 ; zmm0 = a*b + k vmovupd zmmword ptr [r8+rax*8], zmm0 ; stote c[] add rax, 8 cmp rax, r9 jb short loop_f ret afmaaddpd endp
r9レジスターへ要素数を読み込みます。zmm2レジスターへvbroadcastsd命令でzmm3レジスターの最下位要素に格納されている定数kの全要素へブロードキャストします。
一次元配列のインデックスを管理するrax レジスターをxor命令でゼロクリアします。
zmm0レジスターへ配列 a の8要素、zmm1レジスターへ配列 b の8要素を読み込みます。次に、vfmadd213pd命令で「a*b + k」を行い、zmm0レジスターへ結果を得ます。この値を、配列 c へ書き込みます。
インデックスを保持するraxレジスターへ8を加算し、これを配列の要素数を保持するr9レジスターと比較し、すべて処理したか調べます。処理完了なら関数を抜けます。そうでなければ、次の要素の処理へ移ります。
vfmadd213pd 命令
この命令は、三つのソースオペランドを使用してパックド倍精度浮動小数点値に対し、融合積和（FMA: fused multiply-add）を実行します。デスティネーションオペランドは、最初のソースオペランドでもあります。2 番目のオペランドはSIMD レジスターでなければなりません。3 番目のソースオペランドは、SIMD レジスターまたはメモリーロケーションです。
8 個（レジスターにより2 か4 のときもある）のパックド倍精度浮動小数点値を、
1 番目のオペランドと2番目のオペランドの8 個のパックド倍精度浮動小数点値乗算し、無限精度の中間結果を3 番目オペランドと加算します。
vfmadd213pd zmm0, zmm1, zmm2
と指定した場合、zmm0とzmm1レジスターのパックド倍精度浮動小数点値を乗算し、zmm2レジスターに加算し、結果をzmm0レジスターへ格納します。命令の213 の部分を、132 や231 へ書き換えるとオペランドの対応が変わります。この命令は、マスクレジスターやゼロマスクレジスターにも対応します。
common.h、共通関数
ヘッダーファイルの一部を示します。簡単な内容なので、説明は省きします。
#include <iostream> : // C++ template <typename T> void cfmaadd(const T* a, const T* b, T* c, const T k, const size_t length) { for (size_t i = 0; i < length; i++) { c[i] = a[i] * b[i] + k; } } // C++ template <typename T> void cfmasub(const T* a, const T* b, T* c, const T k, const size_t length) { for (size_t i = 0; i < length; i++) { c[i] = a[i] * b[i] - k; } } // C++ template <typename T> void cfnmaadd(const T* a, const T* b, T* c, const T k, const size_t length) { for (size_t i = 0; i < length; i++) { c[i] = -a[i] * b[i] + k; } } // C++ template <typename T> void cfnmasub(const T* a, const T* b, T* c, const T k, const size_t length) { for (size_t i = 0; i < length; i++) { c[i] = -a[i] * b[i] - k; } } // C++ template <typename T> void cfmaddsub(const T* a, const T* b, T* c, const T k, const size_t length) { for (size_t i = 0; i < length; i++) { if (i % 2 == 0) c[i] = a[i] * b[i] - k; else c[i] = a[i] * b[i] + k; } } // C++ template <typename T> void cfmsubadd(const T* a, const T* b, T* c, const T k, const size_t length) { for (size_t i = 0; i < length; i++) { if (i % 2 == 0) c[i] = a[i] * b[i] + k; else c[i] = a[i] * b[i] - k; } } :
すべての型に対応させるため、テンプレート関数とします。cfmaadd、cfmasub、cfnmaadd、cfnmasub、cfmaddsub、そしてcfmsubadd関数を呼び出して得た結果と、アセンブリ言語で開発した関数を呼び出した結果を比較します。単純比較すると誤差の関係でエラーになるときがあるため、誤差を考慮して比較します。
cfmaadd、cfmasub、cfnmaadd、cfnmasub、cfmaddsub、そしてcfmsubadd関数は、アセンブリ言語で開発した関数をC++言語で記述したものです。
このプログラムのビルドし、実行した様子を示します。
実行結果
C:\>ml64 /c fmaPDAsm.asm
…
C:\>cl /O2 /EHsc fmaPD.cpp fmaPDAsm.obj
…
C:\>fmaPD
---[0]--- Ok!
---[1]--- Ok!
---[2]--- Ok!
---[3]--- Ok!
---[4]--- Ok!
---[5]--- Ok!
正常に処理されたようです。
