CV_8UC3 to rawをチューニングしたら、実行命令数が1/22になった(後編)

TL;DR
=====
- AVX2+SIMD256対応したら、実行命令数が元から見て1/22に削減できた

| Code                    | Ir        |
| -----------------------:|:----------|
| OpenCV 4.9.0            | 5,911,406 |
| SIMD無チューニング       | 2,453,868 |
| SSE3+SIMD128            | 1,150,789 |
| SSE4.1+SIMD128          |   325,695 |
| AVX2+SIMD256            |   260,544 |

これまでのあらすじ
=================
Wayland実装の中の、RAW変換部分をチューニングしていこう！！

前回は、SIMD無しでチューニングしていた。 https://qiita.com/hon_no_mushi/items/28be7acb2f8466291414

今回は、SIMD対応の話！

SIMDを使うレベルでのチューニング(SIMD128)
========================================
さて、ここまでは「普通」の実装でのチューニングをしてきた。

次は、SIMD実装をしていこう、使うのは「Universal Intrinsics」だ！ https://docs.opencv.org/4.x/d6/dd1/tutorial_univ_intrin.html

今度はソースコードから先に見ていこう。

```cpp
    // Convert from [b8:g8:r8] to [b8:g8:r8:x8]
    for (int y = 0; y < img_rows; y++)
    {
        const uint8_t* src = (uint8_t*)img.ptr(y);
        int x = 0;
#if CV_SIMD
#if CV_SIMD128
        for (; x < img_cols - 16; x+=16, src+=16*3, dst+=16*4)
        {
            cv::v_uint8x16 vB, vG, vR;
            cv::v_load_deinterleave(src, vB, vG, vR);     // BGR
            cv::v_store_interleave (dst, vB, vG, vR, vR); // BGRx (x is any).
        }
#endif // CV_SIMD128
#endif // CV_SIMD

        // tail
        for (; x < img_cols; x++, src+=3, dst+=4)
        {
            dst[0] = src[0];
            dst[1] = src[1];
            dst[2] = src[2];
        }
    }
```

- `CV_SIMD`は、SIMD命令を使える環境をターゲットにしますよ、という意味。
- `CV_SIMD128`は、レジスタ長さ128bitを使える環境をターゲットにしますよ、という意味。

ということで、uint8_tが16個ならんでいる、`v_uint8x16`型を使って、`vB`,`vG`,`vR`を定義する。

srcアドレスからの読み込み
------------------------
`src`のアドレスから内容をloadするときに"deinterleave"指定すると、B成分、G成分、R成分を抽出できる。

```plaintext
src = {
  [B0,G0,R0][B1,G1,R1][B2,G2,R2][B3,G3,R3][B4,G4,R4][B5,G5,R5][B6,G6,R6][B7,G7,R7]
  [B8,G8,R8][B9,G9,R9][Ba,Ga,Ra][Bb,Gb,Rb][Bc,Gc,Rc][Bd,Gd,Rd][Be,Ge,Re][Bf,Gf,Rf] }

> cv::v_load_deinterleave(src, vB, vG, vR);     // BGR

vB = { B0, B1, B2, B3, B4, B5, B6, B7, B8, B9, Ba, Bb, Bc, Bd, Be, Bf }
vG = { G0, G1, G2, G3, G4, G5, G6, G7, G8, G9, Ga, Gb, Gc, Gd, Ge, Gf }
vR = { R0, R1, R2, R3, R4, R5, R6, R7, R8, R9, Ra, Rb, Rc, Rd, Re, Rf }

```

dstアドレスへの書き込み
----------------------
`dst`のアドレスへ内容をstoreするときに"interleave"指定すると、B成分、G成分、R成分をまとめる事ができる。この時、4ch目にダミーとしてR成分を付けておけば、XRGBとして4chにすることができる。

```plaintext
vB = { B0, B1, B2, B3, B4, B5, B6, B7, B8, B9, Ba, Bb, Bc, Bd, Be, Bf }
vG = { G0, G1, G2, G3, G4, G5, G6, G7, G8, G9, Ga, Gb, Gc, Gd, Ge, Gf }
vR = { R0, R1, R2, R3, R4, R5, R6, R7, R8, R9, Ra, Rb, Rc, Rd, Re, Rf }

> cv::v_store_interleave(dst, vB, vG, vR, vR);     // BGRx (x is any)

dst = {
  [B0,G0,R0,R0][B1,G1,R1,R1][B2,G2,R2,R2][B3,G3,R3,B3]
  [B4,G4,R4,R4][B5,G5,R5,R5][B6,G6,R6,R6][B7,G7,R7,R7]
  [B8,G8,R8,R8][B9,G9,R9,R9][Ba,Ga,Ra,Ra][Bb,Gb,Rb,Rb]
  [Bc,Gc,Rc,Rc][Bd,Gd,Rd,Rd][Be,Ge,Re,Re][Bf,Gf,Rf,Rf] }

```

```bash
valgrind --tool=callgrind ./a.out
==24916== Callgrind, a call-graph generating cache profiler
==24916== Copyright (C) 2002-2017, and GNU GPL'd, by Josef Weidendorfer et al.
==24916== Using Valgrind-3.22.0 and LibVEX; rerun with -h for copyright info
==24916== Command: ./a.out
==24916==
==24916== For interactive control, run 'callgrind_control -h'.
cv::currentUIFramework() returns WAYLAND
==24916==
==24916== Events    : Ir
==24916== Collected : 123413086
==24916==
==24916== I   refs:      123,413,086
kmtr@kmtr-VMware-Virtual-Platform:~/work/build4-main/temp$ callgrind_annotate callgrind.out.24916 | grep 8888 | head -1
 1,150,789 ( 0.93%)  /usr/lib/gcc/x86_64-linux-gnu/13/include/emmintrin.h:write_mat_to_xrgb8888(cv::Mat const&, void*)
```

SIMDを使うレベルでのチューニング(SIMD256)
========================================
私の中のハンチョウが囁いている...

"フフ……へただなあ、熊太郎くん。へたっぴさ……! 欲望の解放のさせ方がへた……。
 熊太郎君が本当に欲しいのは… **SIMD256** ……こっち……
これを **AVX2** でチンして ホッカホッカにしてさ、core i7で回したい……!だろ……？ "

ということで、`ccmake`コマンドを使って、BASELINEを引き上げる。

```
 CPU_BASELINE                     AVX2
 CPU_DISPATCH                     SSE4_1;SSE4_2;AVX;FP16;AVX2;AVX512_SKX
 ```

コードもSIMD256に対応させておく。

```cpp
static void write_mat_to_xrgb8888(cv::Mat const &img, void *data) {
    CV_CheckTrue(data != nullptr, "data must not be nullptr.");
    CV_CheckType(img.type(), img.type() == CV_8UC3, "Only 8UC3 images are supported.");

    int img_rows = img.rows;
    int img_cols = img.cols;
    uint8_t *dst = (uint8_t*) data;

    // to reduce calling img.ptr()
    if(img.isContinuous())
    {
        img_cols *= img_rows;
        img_rows  = 1;
    }

    // Convert from [b8:g8:r8] to [b8:g8:r8:x8]
    for (int y = 0; y < img_rows; y++)
    {
        const uint8_t* src = (uint8_t*)img.ptr(y);
        int x = 0;
#if CV_SIMD
#if CV_SIMD256
        for (; x < img_cols - 32; x+=32, src+=32*3, dst+=32*4)
        {
            cv::v_uint8x32 vB, vG, vR;
            cv::v_load_deinterleave(src, vB, vG, vR);     // BGR
            cv::v_store_interleave (dst, vB, vG, vR, vR); // BGRx (x is any).
        }
#endif // CV_SIMD256
#if CV_SIMD128
        for (; x < img_cols - 16; x+=16, src+=16*3, dst+=16*4)
        {
            cv::v_uint8x16 vB, vG, vR;
            cv::v_load_deinterleave(src, vB, vG, vR);     // BGR
            cv::v_store_interleave (dst, vB, vG, vR, vR); // BGRx (x is any).
        }
#endif // CV_SIMD128
#endif // CV_SIMD

        // tail
        for (; x < img_cols; x++, src+=3, dst+=4)
        {
            dst[0] = src[0];
            dst[1] = src[1];
            dst[2] = src[2];
        }
    }
}
```

```
 valgrind --tool=callgrind ./a.out
==28315== Callgrind, a call-graph generating cache profiler
==28315== Copyright (C) 2002-2017, and GNU GPL'd, by Josef Weidendorfer et al.
==28315== Using Valgrind-3.22.0 and LibVEX; rerun with -h for copyright info
==28315== Command: ./a.out
==28315==
==28315== For interactive control, run 'callgrind_control -h'.
cv::currentUIFramework() returns WAYLAND
==28315==
==28315== Events    : Ir
==28315== Collected : 122295687
==28315==
==28315== I   refs:      122,295,687
kmtr@kmtr-VMware-Virtual-Platform:~/work/build4-main/temp$ callgrind_annotate callgrind.out.28315 | grep 8888 | head -1
   260,544 ( 0.21%)  /usr/lib/gcc/x86_64-linux-gnu/13/include/avx2intrin.h:write_mat_to_xrgb8888(cv::Mat const&, void*)
```

えーっと、単純な命令数で比較なので、一概には言えないけど・・・ 半分以下？（大混乱）

ちょっともって、元々のbuild configurationに戻ると…
```
--   CPU/HW features:
--     Baseline:                    SSE SSE2 SSE3
--       requested:                 SSE3
--     Dispatched code generation:  SSE4_1 SSE4_2 FP16 AVX AVX2 AVX512_SKX
--       requested:                 SSE4_1 SSE4_2 AVX FP16 AVX2 AVX512_SKX
--       SSE4_1 (18 files):         + SSSE3 SSE4_1
--       SSE4_2 (2 files):          + SSSE3 SSE4_1 POPCNT SSE4_2
--       FP16 (1 files):            + SSSE3 SSE4_1 POPCNT SSE4_2 FP16 AVX
--       AVX (9 files):             + SSSE3 SSE4_1 POPCNT SSE4_2 AVX
--       AVX2 (38 files):           + SSSE3 SSE4_1 POPCNT SSE4_2 FP16 FMA3 AVX AVX2
--       AVX512_SKX (8 files):      + SSSE3 SSE4_1 POPCNT SSE4_2 FP16 FMA3 AVX AVX2 AVX_512F AVX512_COMMON AVX512_SKX
```

**SSE3** !?!? Pentium4時代の化石じゃないですかね！！（目血走り）

SSE3の場合
----------
https://github.com/opencv/opencv/blob/4.9.0/modules/core/include/opencv2/core/hal/intrin_sse.hpp

40 instruction / 16 elements

```cpp
inline void v_store_interleave( uchar* ptr, const v_uint8x16& a, const v_uint8x16& b,
                                const v_uint8x16& c, hal::StoreMode mode = hal::STORE_UNALIGNED)
{
#if CV_SSE4_1
...
#elif CV_SSSE3
...
#else
    __m128i z = _mm_setzero_si128();
    __m128i ab0 = _mm_unpacklo_epi8(a.val, b.val);
    __m128i ab1 = _mm_unpackhi_epi8(a.val, b.val);
    __m128i c0 = _mm_unpacklo_epi8(c.val, z);
    __m128i c1 = _mm_unpackhi_epi8(c.val, z);

    __m128i p00 = _mm_unpacklo_epi16(ab0, c0);
    __m128i p01 = _mm_unpackhi_epi16(ab0, c0);
    __m128i p02 = _mm_unpacklo_epi16(ab1, c1);
    __m128i p03 = _mm_unpackhi_epi16(ab1, c1);

    __m128i p10 = _mm_unpacklo_epi32(p00, p01);
    __m128i p11 = _mm_unpackhi_epi32(p00, p01);
    __m128i p12 = _mm_unpacklo_epi32(p02, p03);
    __m128i p13 = _mm_unpackhi_epi32(p02, p03);

    __m128i p20 = _mm_unpacklo_epi64(p10, p11);
    __m128i p21 = _mm_unpackhi_epi64(p10, p11);
    __m128i p22 = _mm_unpacklo_epi64(p12, p13);
    __m128i p23 = _mm_unpackhi_epi64(p12, p13);

    p20 = _mm_slli_si128(p20, 1);
    p22 = _mm_slli_si128(p22, 1);

    __m128i p30 = _mm_slli_epi64(_mm_unpacklo_epi32(p20, p21), 8);
    __m128i p31 = _mm_srli_epi64(_mm_unpackhi_epi32(p20, p21), 8);
    __m128i p32 = _mm_slli_epi64(_mm_unpacklo_epi32(p22, p23), 8);
    __m128i p33 = _mm_srli_epi64(_mm_unpackhi_epi32(p22, p23), 8);

    __m128i p40 = _mm_unpacklo_epi64(p30, p31);
    __m128i p41 = _mm_unpackhi_epi64(p30, p31);
    __m128i p42 = _mm_unpacklo_epi64(p32, p33);
    __m128i p43 = _mm_unpackhi_epi64(p32, p33);

    __m128i v0 = _mm_or_si128(_mm_srli_si128(p40, 2), _mm_slli_si128(p41, 10));
    __m128i v1 = _mm_or_si128(_mm_srli_si128(p41, 6), _mm_slli_si128(p42, 6));
    __m128i v2 = _mm_or_si128(_mm_srli_si128(p42, 10), _mm_slli_si128(p43, 2));
#endif

...
}

```

AVX2の場合
----------
https://github.com/opencv/opencv/blob/4.9.0/modules/core/include/opencv2/core/hal/intrin_avx.hpp#L2749-L2778

17 instruction / 32 elements

```
inline void v_store_interleave( uchar* ptr, const v_uint8x32& a, const v_uint8x32& b, const v_uint8x32& c,
                                hal::StoreMode mode=hal::STORE_UNALIGNED )
{
    const __m256i sh_b = _mm256_setr_epi8(
            0, 11, 6, 1, 12, 7, 2, 13, 8, 3, 14, 9, 4, 15, 10, 5,
            0, 11, 6, 1, 12, 7, 2, 13, 8, 3, 14, 9, 4, 15, 10, 5);
    const __m256i sh_g = _mm256_setr_epi8(
            5, 0, 11, 6, 1, 12, 7, 2, 13, 8, 3, 14, 9, 4, 15, 10,
            5, 0, 11, 6, 1, 12, 7, 2, 13, 8, 3, 14, 9, 4, 15, 10);
    const __m256i sh_r = _mm256_setr_epi8(
            10, 5, 0, 11, 6, 1, 12, 7, 2, 13, 8, 3, 14, 9, 4, 15,
            10, 5, 0, 11, 6, 1, 12, 7, 2, 13, 8, 3, 14, 9, 4, 15);


    __m256i b0 = _mm256_shuffle_epi8(a.val, sh_b);
    __m256i g0 = _mm256_shuffle_epi8(b.val, sh_g);
    __m256i r0 = _mm256_shuffle_epi8(c.val, sh_r);


    const __m256i m0 = _mm256_setr_epi8(0, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0,
                                               0, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0);
    const __m256i m1 = _mm256_setr_epi8(0, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0,
                                               0, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0);


    __m256i p0 = _mm256_blendv_epi8(_mm256_blendv_epi8(b0, g0, m0), r0, m1);
    __m256i p1 = _mm256_blendv_epi8(_mm256_blendv_epi8(g0, r0, m0), b0, m1);
    __m256i p2 = _mm256_blendv_epi8(_mm256_blendv_epi8(r0, b0, m0), g0, m1);


    __m256i bgr0 = _mm256_permute2x128_si256(p0, p1, 0 + 2*16);
    __m256i bgr1 = _mm256_permute2x128_si256(p2, p0, 0 + 3*16);
    __m256i bgr2 = _mm256_permute2x128_si256(p1, p2, 1 + 3*16);
```

ついでに以下が、SSSE4.1 + SIMD128での検証結果。
```
valgrind --tool=callgrind ./a.out
==35400== Callgrind, a call-graph generating cache profiler
==35400== Copyright (C) 2002-2017, and GNU GPL'd, by Josef Weidendorfer et al.
==35400== Using Valgrind-3.22.0 and LibVEX; rerun with -h for copyright info
==35400== Command: ./a.out
==35400==
==35400== For interactive control, run 'callgrind_control -h'.
cv::currentUIFramework() returns WAYLAND
==35400==
==35400== Events    : Ir
==35400== Collected : 122546236
==35400==
==35400== I   refs:      122,546,236
kmtr@kmtr-VMware-Virtual-Platform:~/work/build4-main/temp$ callgrind_annotate callgrind.out.35400 | grep 8888 | head
-1
   325,695 ( 0.27%)  /usr/lib/gcc/x86_64-linux-gnu/13/include/emmintrin.h:write_mat_to_xrgb8888(cv::Mat const&, void*)
```

SIMD実装同士の比較
------------------

命令セット | Instruction数 | elements | instruction/32 elements| (参考) Ir |
----------:|---:|---:|---:|---:|
SSE3+SIMD128 | 40 | 16 | 80 | 1,150,789 |
SSE4.1+SIMD128 | 11 | 16 | 22 |   325,695 |
AVX2+SIMD256 | 17 | 32 | 17 |   260,544 |

全体のまとめ
============

| Code                    | Ir        |
| -----------------------:|:----------|
| OpenCV 4.9.0            | 5,911,406 |
| SIMD無チューニング       | 2,453,868 |
| SSE3+SIMD128            | 1,150,789 |
| SSE4.1+SIMD128          |   325,695 |
| AVX2+SIMD256            |   260,544 |

ということで、 5911406 / 260544 = 22.689 ということで、命令数が1/22に減りました。

（命令数が1/22になっても、Clock per instructionだのなんだので、性能が22倍になっているわけではない・・・）

以上です、ご精読ありがとうございました！！