mla -積和命令-

# はじめに
 - この記事は[ひとりNEONアドベントカレンダー2020](https://qiita.com/advent-calendar/2020/neon-only) 11日目の記事です
 - [昨日](https://qiita.com/tomoaki_teshima/items/0c2e1690a1e6efb60f0b)は`store`命令を紹介した。
 - みんな大好き積和命令を紹介する

![mlaを含む命令一覧](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/12162/b7050604-f2dc-1801-9c60-af231197132e.png)

```
$ grep ^v.*mla /usr/lib/gcc/aarch64-linux-gnu/7.5.0/include/arm_neon.h | cut -f 1 -d ' ' | sed -e 's/[suf][0-9]\+//g' | sort | uniq -c
      8 vmla_
      5 vmla_lane_
      5 vmla_laneq_
      5 vmla_n_
      6 vmlal_
      6 vmlal_high_
      4 vmlal_high_n_
      4 vmlal_n_
      8 vmlaq_
      5 vmlaq_lane_
      5 vmlaq_laneq_
      5 vmlaq_n_
      2 vqdmlal_
      2 vqdmlal_high_
      2 vqdmlal_high_lane_
      2 vqdmlal_high_laneq_
      2 vqdmlal_high_n_
      2 vqdmlal_lane_
      2 vqdmlal_laneq_
      2 vqdmlal_n_
      1 vqdmlalh_
      1 vqdmlalh_lane_
      1 vqdmlalh_laneq_
      1 vqdmlals_
      1 vqdmlals_lane_
      1 vqdmlals_laneq_
      2 vqrdmlah_
      2 vqrdmlah_lane_
      2 vqrdmlah_laneq_
      1 vqrdmlahh_
      1 vqrdmlahh_lane_
      1 vqrdmlahh_laneq_
      2 vqrdmlahq_
      2 vqrdmlahq_lane_
      2 vqrdmlahq_laneq_
      1 vqrdmlahs_
      1 vqrdmlahs_lane_
      1 vqrdmlahs_laneq_
```

# `mla`命令

 - 積和命令は、`mla`と`mlaq`が基本形で、`a`、`b`、`c`の3パラメータを受け取り、`a + b * c`を返す
 - `mla_n`、`mla_lane`、`mla_laneq`は、`c`がベクトルでなく、1つの値(ベクトルの中のレーンで指定したり、単一の`int`で指定したりする)
 - `mlal`は、入力のレーン数を半減する代わりに、乗算結果がオーバーフローしないように引数`b``c`のbit拡張を行う
    - `mlal`は下位半分のレーンを利用し、`mlal_high`は上位半分のレーンを使う
    - `mlal_n`と`mlal_high_n`は`b`の半分のレーンを使い、`c`はベクトルでなく単一の値を使う
 - `vqdmla`命令は`a+b*c`でなく、`a+b*c*2`を計算する。
 - `vqrdmla`命令は`a+b*c`でなく、`a+b*c*2+(const)`を計算する。

# `vqdmlal`

```vqdmlal.cpp
#include <arm_neon.h>
#include <iostream>


#define KERNEL(suffix1, suffix2, nSize, type1, type2, vtype1, vtype2) \
void kernel_##suffix1()                                               \
{                                                                     \
        type1 src0[] = { 198, 207,  51, 984};                         \
        type1 src1[] = { 122, 920, 151, 129};                         \
        type2 acc [] = { 10000, 10000, 10000, 10000, };               \
        type2 res [] = {   0,   0,   0,   0};                         \
        int ## vtype2 ## _t a = vld1q_## suffix2(acc);                \
        int ## vtype1 ## _t b = vld1_## suffix1(src0);                \
        int ## vtype1 ## _t c = vld1_## suffix1(src1);                \
        int ## vtype2 ## _t d = vqdmlal_## suffix1(a, b, c);          \
        vst1q_ ## suffix2(res, d);                                    \
        for( auto i = 0;i < nSize;i++)                                \
        {                                                             \
            std::cout << i << ":" << acc[i] << '\t' << src0[i] << '\t' << src1[i] << '\t' << res[i] << std::endl;  \
        }                                                             \
}

KERNEL(s32, s64, 2, int,   int64_t, 32x2, 64x2);
KERNEL(s16, s32, 4, short, int,     16x4, 32x4);

int main(int argc, char** argv)
{
        kernel_s32();
        kernel_s16();
}
```

```
0:10000 198     122     58312
1:10000 207     920     390880
0:10000 198     122     58312
1:10000 207     920     390880
2:10000 51      151     25402
3:10000 984     129     263872
```

 - `a + b * c * 2`が実行される。`b*c*2`はオーバーフローに備えて乗算前に符号拡張される。`vq`が接頭辞なので、最後の加算は飽和演算。
接頭辞が純粋に`v`だけの非飽和演算は存在しない。

## `vqdmlalh` `vqdmlals`
`v` `q` `d` `mla` `l` なので、ベクトル演算、飽和演算、倍加算、積和演算、bit拡張。最後の`h``s`は、`half word`(16bit)`single word`(32bit)に対応する接尾辞で、`_s16`、`_s32`に相当する。`v`が接頭辞だが、単体の `a + b * c * 2`しか行わない。つまり、厳密にはSIMD命令では無い。

# `vqrdmlah`
 - v8.1拡張命令

```arm_neon.h
#pragma GCC target ("+nothing+rdma")
vqrdmlah_s16 (int16x4_t __a, int16x4_t __b, int16x4_t __c)
vqrdmlah_s32 (int32x2_t __a, int32x2_t __b, int32x2_t __c)
vqrdmlahq_s16 (int16x8_t __a, int16x8_t __b, int16x8_t __c)
vqrdmlahq_s32 (int32x4_t __a, int32x4_t __b, int32x4_t __c)
vqrdmlah_laneq_s16 (int16x4_t __a, int16x4_t __b, int16x8_t __c, const int __d)
vqrdmlah_laneq_s32 (int32x2_t __a, int32x2_t __b, int32x4_t __c, const int __d)
vqrdmlahq_laneq_s16 (int16x8_t __a, int16x8_t __b, int16x8_t __c, const int __d)
vqrdmlahq_laneq_s32 (int32x4_t __a, int32x4_t __b, int32x4_t __c, const int __d)
vqrdmlah_lane_s16 (int16x4_t __a, int16x4_t __b, int16x4_t __c, const int __d)
vqrdmlah_lane_s32 (int32x2_t __a, int32x2_t __b, int32x2_t __c, const int __d)
vqrdmlahq_lane_s16 (int16x8_t __a, int16x8_t __b, int16x4_t __c, const int __d)
vqrdmlahq_lane_s32 (int32x4_t __a, int32x4_t __b, int32x2_t __c, const int __d)
vqrdmlahh_s16 (int16_t __a, int16_t __b, int16_t __c)
vqrdmlahh_lane_s16 (int16_t __a, int16_t __b, int16x4_t __c, const int __d)
vqrdmlahh_laneq_s16 (int16_t __a, int16_t __b, int16x8_t __c, const int __d)
vqrdmlahs_s32 (int32_t __a, int32_t __b, int32_t __c)
vqrdmlahs_lane_s32 (int32_t __a, int32_t __b, int32x2_t __c, const int __d)
vqrdmlahs_laneq_s32 (int32_t __a, int32_t __b, int32x4_t __c, const int __d)
#pragma GCC pop_options
```

# `fma`命令
 - 浮動小数点型向けには、`fma`命令も提供される

```
$ grep ^vfma /usr/lib/gcc/aarch64-linux-gnu/7.5.0/include/arm_neon.h
vfma_f64 (float64x1_t __a, float64x1_t __b, float64x1_t __c)
vfma_f32 (float32x2_t __a, float32x2_t __b, float32x2_t __c)
vfmaq_f32 (float32x4_t __a, float32x4_t __b, float32x4_t __c)
vfmaq_f64 (float64x2_t __a, float64x2_t __b, float64x2_t __c)
vfma_n_f32 (float32x2_t __a, float32x2_t __b, float32_t __c)
vfma_n_f64 (float64x1_t __a, float64x1_t __b, float64_t __c)
vfmaq_n_f32 (float32x4_t __a, float32x4_t __b, float32_t __c)
vfmaq_n_f64 (float64x2_t __a, float64x2_t __b, float64_t __c)
vfma_lane_f32 (float32x2_t __a, float32x2_t __b,
vfma_lane_f64 (float64x1_t __a, float64x1_t __b,
vfmad_lane_f64 (float64_t __a, float64_t __b,
vfmas_lane_f32 (float32_t __a, float32_t __b,
vfma_laneq_f32 (float32x2_t __a, float32x2_t __b,
vfma_laneq_f64 (float64x1_t __a, float64x1_t __b,
vfmad_laneq_f64 (float64_t __a, float64_t __b,
vfmas_laneq_f32 (float32_t __a, float32_t __b,
vfmaq_lane_f32 (float32x4_t __a, float32x4_t __b,
vfmaq_lane_f64 (float64x2_t __a, float64x2_t __b,
vfmaq_laneq_f32 (float32x4_t __a, float32x4_t __b,
vfmaq_laneq_f64 (float64x2_t __a, float64x2_t __b,
vfma_f16 (float16x4_t __a, float16x4_t __b, float16x4_t __c)
vfmaq_f16 (float16x8_t __a, float16x8_t __b, float16x8_t __c)
vfmah_lane_f16 (float16_t __a, float16_t __b,
vfmah_laneq_f16 (float16_t __a, float16_t __b,
vfma_lane_f16 (float16x4_t __a, float16x4_t __b,
vfmaq_lane_f16 (float16x8_t __a, float16x8_t __b,
vfma_laneq_f16 (float16x4_t __a, float16x4_t __b,
vfmaq_laneq_f16 (float16x8_t __a, float16x8_t __b,
vfma_n_f16 (float16x4_t __a, float16x4_t __b, float16_t __c)
vfmaq_n_f16 (float16x8_t __a, float16x8_t __b, float16_t __c)
```
 - しかし、謎に`mla`命令も浮動小数点型向けに提供される

```
$ grep ^vmla.*_f /usr/lib/gcc/aarch64-linux-gnu/7.5.0/include/arm_neon.h
vmla_n_f32 (float32x2_t a, float32x2_t b, float32_t c)
vmlaq_n_f32 (float32x4_t a, float32x4_t b, float32_t c)
vmla_f32 (float32x2_t a, float32x2_t b, float32x2_t c)
vmla_f64 (float64x1_t __a, float64x1_t __b, float64x1_t __c)
vmlaq_f32 (float32x4_t a, float32x4_t b, float32x4_t c)
vmlaq_f64 (float64x2_t a, float64x2_t b, float64x2_t c)
vmla_lane_f32 (float32x2_t __a, float32x2_t __b,
vmla_laneq_f32 (float32x2_t __a, float32x2_t __b,
vmlaq_lane_f32 (float32x4_t __a, float32x4_t __b,
vmlaq_laneq_f32 (float32x4_t __a, float32x4_t __b,
```
 - ん？じゃあ`mlaq_f32`と`fmaq_f32`の違いは？
 - [`mlaq_f32`](https://developer.arm.com/architectures/instruction-sets/simd-isas/neon/intrinsics?search=vmlaq_f32)と[`fmaq_f32`](https://developer.arm.com/architectures/instruction-sets/simd-isas/neon/intrinsics?search=vfmaq_f32)の公式ドキュメントを見ると
 - [`mlaq_f32`](https://developer.arm.com/architectures/instruction-sets/simd-isas/neon/intrinsics?search=vmlaq_f32)
 
> Floating-point multiply-add to accumulator
 
  - [`fmaq_f32`](adadad)

> Floating-point fused multiply-add to accumulator

 - となっている。一応`arm_neon.h`も見てみよう

```arm_neon.h
vmlaq_f32 (float32x4_t a, float32x4_t b, float32x4_t c)
{
  return a + b * c;
}
:
vfmaq_f32 (float32x4_t __a, float32x4_t __b, float32x4_t __c)
{
  return __builtin_aarch64_fmav4sf (__b, __c, __a);
}
```
 - なるほど、解説の通り、`mla`命令は本当にただ掛けて足しているだけである
 - 一方で`fma`命令は`fmav4sf`というビルトイン関数を呼んでおり、こちらはたしかにfusedされてそうである。
 - 実際にサンプルコードで試してみよう

```mla.cpp
union suf32
{
        int32_t s;
        uint32_t u;
        float f;
};

void dumpHex(float32x4_t& v)
{
        float work[4];
        vst1q_f32(work, v);
        for(int i = 0;i < 4;i++)
        {
                suf32 conv;
                conv.f = work[i];
                std::cout << i << "\t0x" << std::hex << conv.u << '\t' << conv.f << std::endl;
        }
}

int main(int argc, char**argv)
{
    float src0[] = {-10703290.9296708,-10703290.9296708,-10703290.9296708,-10703290.9296708,};
    float src1[] = {3333.3333, 3333.3333, 3333.3333, 3333.3333,};
    float src2[] = {3210.9876, 3210.9876, 3210.9876, 3210.9876,};
    float32x4_t vsrc0 = vld1q_f32(src0);
    float32x4_t vsrc1 = vld1q_f32(src1);
    float32x4_t vsrc2 = vld1q_f32(src2);
    {
            float32x4_t v_dst = vmlaq_f32(vsrc0, vsrc1, vsrc2); // vsrc0 + vsrc1 * vsrc2
            dumpHex(v_dst);
    }
    {
            float32x4_t v_dst = vfmaq_f32(vsrc0, vsrc1, vsrc2); // vsrc0 + vsrc1 * vsrc2
            dumpHex(v_dst);
    }
    return 0;
}
```
 - このコードは、fusedされてるどうかを判定するために、両方のintrinsicを利用する
    - また、fusedされてる場合は桁落ちしないが、ただ掛けて足すだけだと桁落ちが発生する。
    - 桁落ちを再現するために、積は$3333.3333 \times 3210.9876 = 10703291.89296708$、積和の結果は$0.96329628$と、 **積に対して演算結果が非常に小さく**なるように設定した
    - 通常の`float`の演算では、桁落ちが発生し、結果が一致しない。
    - fusedされてる場合は途中の積が丸められないため、小さい桁まで保持される
 - では実際に実行してみよう
 - 結論から先に書くと、最適化レベルに応じて`mla`命令の挙動が変わったが、`fma`命令の挙動は変わらなかった。
 - その点は注意が必要である（なお、GCC 7.5.0で確認)
 - 最適化レベル`-O2`

```
0       0x3f117011      0.568116    // mlaの結果
1       0x3f117011      0.568116    // mlaの結果
2       0x3f117011      0.568116    // mlaの結果
3       0x3f117011      0.568116    // mlaの結果
0       0x3f117011      0.568116    // fmaの結果
1       0x3f117011      0.568116    // fmaの結果
2       0x3f117011      0.568116    // fmaの結果
3       0x3f117011      0.568116    // fmaの結果
```
 - 最適化レベル`-O1`

```
0       0x3f800000      1           // mlaの結果
1       0x3f800000      1           // mlaの結果
2       0x3f800000      1           // mlaの結果
3       0x3f800000      1           // mlaの結果
0       0x3f117011      0.568116    // fmaの結果
1       0x3f117011      0.568116    // fmaの結果
2       0x3f117011      0.568116    // fmaの結果
3       0x3f117011      0.568116    // fmaの結果
```
 - 10進数表記と16進数表記を併記した。
 - 最適化レベル`-O2`では、`mla`の結果と`fma`の結果は`bit-exact`である。
 - 一方で、`-O1`では`mla`と`fma`はピッタリは一致しない。`mla`の結果の、下位bitが多数0になっている
    - これが桁落ちと呼ばれる現象
 - という訳で
    - `mla`命令は掛けて足すだけ（ただしコンパイラによっては最適化によってfusedされることがある
    - `fma`命令はかならずfusedされる
    - ことが分かった

# おわりに
- 積和命令を紹介した
- 明日も手島の執筆の予定ですが、そろそろネタが尽きるので、明日は何を書くか未定です