読者です 読者をやめる 読者になる 読者になる

koturnの日記

転職したい社会人1年生の技術系日記

SIMDの組み込み関数のことはじめ

はじめに

現代のCPUではSIMD(Single Instruction Multiple Data)命令を利用することができる. SIMD命令とはその名の通り,ひとつの命令で複数のデータを処理するものである.

Intel系のCPUでは,MMX/SSE/AVX/AVX-512といったSIMD命令が利用可能であり,ARM CPUではNEONというSIMD命令が用意されている. 各SIMDSIMD用のレジスタの対応関係は以下のようになる.

項目 利用可能レジスタ
MMX 64bit のMMレジスタ
SSE 128bit のXMMレジスタ
AVX 256bit のYMMレジシタ
AVX-512 512bit のZMMレジシタ
ARM NEON 64bitのD(Double-Word)レジスタおよび128bitのQ(Quad-Word)レジスタ

これらのレジスタを用いて,例えば4つのint型を一気に処理するといったことを行うのがCPUにおけるSIMDである.

この記事では,このSIMD命令をC/C++から利用することについて記述する.

2017/02/20 追記

以下の記事に,より詳細な内容を書いたので,参考になるかもしれない.

SIMDをプログラム利用するには

SIMD命令というと小難しそうで,インラインアセンブラを利用しなければならないかというと,そうではない. C/C++から関数の形で利用できるように,各コンパイラで共通のAPIである組み込み関数が提供されている. 組み込み関数とはいえ関数なので,関数呼び出しの形で記述することになるが,実際に関数呼び出しが発生するわけではなく,インライン展開され,対応するアセンブラの命令へとコード生成される.

なお,SIMDレジスタに対して,メモリのロードやストアを行う場合,後述するように利用幅と同じ境界に配置されている位置に対して行う必要がある. 特に,MMX/SSE/AVX/AVX-512の場合,アラインメント条件を満たさなければ,SEGVで落ちる関数がある. 落ちない版の関数もあるが,そういった関数は落ちる関数より動作としては遅い.

ARM NEONは落ちる関数は無いが,アラインメント条件を満たしておいた方が高速に動作すると思われる.

インクルード

何はともあれ,まず組み込み関数が宣言されているヘッダをインクルードしなければ始まらない. 各SIMD命令セットとヘッダの対応関係は以下のようになる.

命令セット ヘッダファイル
MMX <mmintrin.h>
SSE <xmmintrin.h>
SSE2 <emmintrin.h>
SSE3 <pmmintrin.h>
SSSE3 <tmmintrin.h>
SSE4.1 <smmintrin.h>
SSE4.2 <nmmintrin.h>
AES <wmmintrin.h>
AVX, AVX2, FMA <immintrin.h>
AVX-512 <zmmintrin.h>
ARM NEON <arm_neon.h>

MMX/SSE/AVX/AVX-512関連のヘッダは多く,これらをいちいちインクルードするのは面倒である. 現実的にはまとめてインクルードすることが可能なヘッダを利用するのがよい. ただし,MSVCとgcc/clangでヘッダが異なるため,注意しなければならない.

環境 ヘッダファイル
MSVC <intrin.h>
gcc/clang <x86intrin.h>

具体的なインクルード部分のコードを書くと以下のようになる.

#ifdef _MSC_VER
#  include <intrin.h>
#else
#  include <x86intrin.h>
#endif

なお,gcc/clangでも,x64環境ならば <intrin.h> が存在するが,x86環境でも利用可能な方に合わせておく方が何かと都合が良いだろう.

コンパイルオプション

実はヘッダをインクルードするだけではSIMDの組み込み関数は利用できない. 以下のようにコンパイルオプションを指定する必要がある.

命令セット gccのオプション MSVCのオプション 定義されるマクロ
MMX -mmmx /arch:MMX __MMX__
SSE -msse /arch:SSE __SSE__
SSE2 -msse2 /arch:SSE2 __SSE2__
SSE3 -msse3 __SSE3__
SSSE3 -mssse3 __SSSE3__
SSE4.1 -msse4.1 __SSE4_1__
SSE4.2 -msse4.2 __SSE4_2__
AES -maes __AES__
AVX -mavx /arch:AVX __AVX__
AVX2 -mavx2 /arch:AVX2 __AVX2__
FMA -mfma __FMA__
AVX-512 -mavx512* ( *bw, cq, ed など) __AVX512*__
ARM NEON -mfpu=neon など __ARM_NEON または __ARM_NEON__

MMX/SSE/AVX/AVX-512関連のオプションは, -march=native-mtune=native などを指定することで,一括で上記のオプションのうち,利用可能なものを指定できる. ARM CPU環境のgccでは, -march=native-mtune=native と指定することができない場合があり,そのときは利用しているARM CPUに合わせて, -fpu=neon-fp-armv8 などと指定する必要がある(これはRaspberry Pi 3の例).

上記の表では簡略に紹介したが,gccのAVX-512に関するオプションは以下のように多数ある.

  • -mavx512f
  • -mavx512er
  • -mavx512cd
  • -mavx512pf
  • -mavx512dq
  • -mavx512bw
  • -mavx512vl
  • -mavx512ifma
  • -mavx512vbmi

なお,現在のところAVX-512が利用できるCPUは限られている. -march=native を指定したとしても,AVX-512が有効にならない場合の方が多いので,上記のオプションを別途指定すると,コンパイルだけは通るだろう. しかし,非対応のCPUでAVX-512命令を実行したとても,以下のようなエラーメッセージが出力されるだろう. (これはMSYS2でzsh上で実行した結果である)

$ ./main.exe
zsh: illegal hardware instruction  ./main.exe

AVX-512の動作を確認するだけならば,Intel公式のエミュレータを利用するとよい. 予め,AVX-512命令が含まれる実行バイナリを生成し,以下のように実行する.

$ sde -- ./main.exe

変数のアラインメントを指定する

C++11,C11から言語の標準機能として,変数のアラインメントを指定することができるようになったが,それ以前は変数のアラインメントはコンパイラ独自の機能を利用しなければ,指定することができない. 古いコンパイラコンパイルすることを考慮すると,以下のように差を吸収するマクロを定義するとよい.

#include <cstddef>
#include <iostream>

#if defined(__cplusplus) && __cplusplus < 201103L
#  ifdef _MSC_VER
#    define alignas(n)  __declspec(align(n))
#  else
#    define alignas(n)  __attribute__((aligned(n)))
#  endif  // _MSC_VER
#endif  // defined(__cplusplus) && __cplusplus < 201103L


// 以下,利用コード


int
main()
{
  static const int ALIGN = 32;
  alignas(ALIGN) unsigned char array[10] = {0};
  if ((reinterpret_cast<std::ptrdiff_t>(array)) % ALIGN == 0) {
    std::cout << "Static array is " << ALIGN << " byte aligned.\n";
  } else {
    std::cout << "Static array is not " << ALIGN << " byte aligned.\n";
  }

  return 0;
}

上記はC++用だが,C言語なら以下のように定義するとよい.

#include <stddef.h>
#include <stdio.h>

#if defined(__STDC_VERSION__) && __STDC_VERSION__ < 201102L
#  ifdef _MSC_VER
#    define _Alignas(n)  __declspec(align(n))
#  else
#    define _Alignas(n)  __attribute__((aligned(n)))
#  endif  // _MSC_VER
#endif  // defined(__cplusplus) && __cplusplus < 201103L


/* 以下,利用コード */


#define ALIGN  32

int
main(void)
{
  _Alignas(ALIGN) unsigned char array[10] = {0};
  if ((ptrdiff_t) array % ALIGN == 0) {
    printf("Static array is %d byte aligned.\n", ALIGN);
  } else {
    printf("Static array is not %d byte aligned.\n", ALIGN);
  }

  return 0;
}

アラインされたメモリを動的確保する

通常のC/C++における std::malloc()std::calloc()new 等では16byteや32byte境界にアラインメントされたメモリを動的確保することはできない. 以下に示す専用のメモリ確保関数が必要となる. (C言語の場合, <cstdlib><stdlib.h> に読み換えること)

メモリ確保関数 メモリ解放関数 ヘッダ 特徴
_aligned_malloc() _aligned_free() <malloc.h> MSVCのみ.
posix_memalign() std::free() <cstdlib> gcc/clangのみ.
aligned_alloc() std::free() <cstdlib> gcc/clangのみ.確保サイズはアラインメントの整数倍に限る
memalign() std::free() <malloc.h> gcc/clangのみ.廃止されているとのこと.
_mm_malloc() _mm_free() <malloc.h> Intel CPUのみ.

種々のアラインされたメモリ確保関数があり,どれを利用すればいいか判断に困るかもしれない. しかし,おおまかには,以下のように利用する関数を判断すればよい.

  • MSVCなら _aligned_malloc()_aligned_free()
  • gcc/clangなら posix_memalign()std::free()

これを考慮し,条件コンパイルで利用する関数を分岐するラッパー関数を作るとよい. 簡単なコードは以下のようになる.

なお,C++11以降, std::align()std::aligned_storage() といった関数が利用できるが, std::align() は既に確保されたバッファの指定されたアドレスからポインタを進め,アラインメント条件を満たす位置のアドレスを返却するだけの関数であり, std::aligned_storage() はアラインされた静的配列を作成するための関数なので,やや使い勝手が悪いといえる.

// <type_traits> はC++11以降のものなので,それ以前でコンパイルしたい場合は関連部分を削除すること
#include <cstddef>
#include <iostream>
#include <memory>
#include <type_traits>
#if defined(_MSC_VER) || defined(__MINGW32__)
#  include <malloc.h>
#else
#  include <cstdlib>
#endif  // defined(_MSC_VER) || defined(__MINGW32__)


/*!
 * @brief アラインメントされたメモリを動的確保する関数
 * @param [in] size       確保するメモリサイズ (単位はbyte)
 * @param [in] alignment  アラインメント (2のべき乗を指定すること)
 * @return  アラインメントし,動的確保されたメモリ領域へのポインタ
 */
template<typename T = void*, typename std::enable_if<std::is_pointer<T>::value, std::nullptr_t>::type = nullptr>
static inline T
alignedMalloc(std::size_t size, std::size_t alignment) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  return reinterpret_cast<T>(_aligned_malloc(size, alignment));
#else
  void* p;
  return reinterpret_cast<T>(posix_memalign(&p, alignment, size) == 0 ? p : nullptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


/*!
 * @brief アラインメントされたメモリを解放する関数
 * @param [in] ptr  解放対象のメモリの先頭番地を指すポインタ
 */
static inline void
alignedFree(void* ptr) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  _aligned_free(ptr);
#else
  std::free(ptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


// 以下,利用コード


/*!
 * @brief std::unique_ptr で利用するアラインされたメモリ用のカスタムデリータ
 */
struct AlignedDeleter
{
  void
  operator()(void* p) const noexcept
  {
    alignedFree(p);
  }
};


int
main()
{
  static constexpr int ALIGN = 32;
  std::unique_ptr<unsigned char[], AlignedDeleter> array(alignedMalloc<unsigned char*>(10 * sizeof(unsigned char), ALIGN));
  if (array.get() == nullptr) {
    std::cerr << "Failed to allocate memory" << std::endl;
    return 1;
  }

  if ((reinterpret_cast<std::ptrdiff_t>(array.get())) % ALIGN == 0) {
    std::cout << "Dynamic allocated memory is " << ALIGN << " byte aligned.\n";
  } else {
    std::cout << "Dynamic allocated memory is not " << ALIGN << " byte aligned.\n";
  }

  return 0;
}

このコードはC++11の範疇のものであるが,C言語の範囲で書き直すと以下のようになる. C99以降は inline が利用可能であるが,古いコンパイラを使用することを考慮し,置き換えるマクロを記述する.

#include <stdio.h>
#include <stddef.h>
#if defined(_MSC_VER) || defined(__MINGW32__)
#  include <malloc.h>
#else
#  include <stdlib.h>
#endif  /* defined(_MSC_VER) || defined(__MINGW32__) */

#ifndef __cplusplus
#  if defined(_MSC_VER)
#    define inline      __inline
#    define __inline__  __inline
#  elif !defined(__GNUC__) && !defined(__STDC_VERSION__) || __STDC_VERSION__ < 199901L
#    define inline
#    define __inline
#  endif
#endif


/*!
 * @brief アラインメントされたメモリを動的確保する関数
 * @param [in] size       確保するメモリサイズ (単位はbyte)
 * @param [in] alignment  アラインメント (2のべき乗を指定すること)
 * @return  アラインメントし,動的確保されたメモリ領域へのポインタ
 */
static inline void*
alignedMalloc(size_t size, size_t alignment)
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  return _aligned_malloc(size, alignment);
#else
  void* p;
  return posix_memalign((void**) &p, alignment, size) == 0 ? p : NULL;
#endif  /* _MSC_VER */
}


/*!
 * @brief アラインメントされたメモリを解放する関数
 * @param [in] ptr  解放対象のメモリの先頭番地を指すポインタ
 */
static inline void
alignedFree(void* ptr)
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  _aligned_free(ptr);
#else
  free(ptr);
#endif  /* _MSC_VER */
}


/* 以下,利用コード */


int
main(void)
{
  static const int ALIGN = 32;
  unsigned char* array = (unsigned char*) alignedMalloc(10 * sizeof(unsigned char), ALIGN);
  if (array == NULL) {
    fprintf(stderr, "Failed to allocate memory\n");
    return 1;
  }

  if (((ptrdiff_t) array) % ALIGN == 0) {
    printf("Dynamic allocated memory is %d byte aligned.\n", ALIGN);
  } else {
    printf("Dynamic allocated memory is not %d byte aligned.\n", ALIGN);
  }
  alignedFree(array);

  return 0;
}

なお,今回は適当なアラインメントを指定したが,実際にSSE/AVX/AVX-512/NEONを用いるときは,SSE/AVX/AVX-512/NEONの変数型からアラインメントを取得するとよい. 型や変数からアラインメントを得る機能はC++11およびC11以降であれば, alignof 演算子で取得でき,それ以前の環境であれば,コンパイラ拡張機能を用いることで取得できる. この差を吸収するなら,以下のようなマクロを定義するとよい.

#if defined(__cplusplus) && __cplusplus < 201103L
#  ifdef _MSC_VER
#    define alignof(n)  __alignof(n)
#  else
#    define alignof(n)  __alignof__(n)
#  endif  // _MSC_VER
#endif  // defined(__cplusplus) || cplusplus < 201103L

// alingas(alignof(__m256i)) のような形で使用

SSE/AVX/NEON のサンプルコード

簡単なサンプルコードをSSE/AVX/NEONの例として提示する. このコードはMSVC/gc/clangのいずれのコンパイラでもコンパイルすることができるようにしている.

AVX-512については,利用可能なCPUを搭載したマシンが手元に無いため割愛するが,AVXと同様のコードで記述できると思う. コンパイル時に以下のマクロを定義すると,対応した命令を用いたコードが有効化される. 有効化しようとしても,コンパイラが対応していない場合は,冒頭の部分でエラーが発生するはずだ. また,以下のいずれのマクロも定義しなかった場合,SIMDを用いないコードとなる.

マクロ 有効化されるSIMD
ENABLE_AVX AVX
ENABLE_SSE SSE
ENABLE_NEON ARM NEON

直接的に __AVX____SSE2__ 等のマクロが定義されているかどうかで判断しないのは,コンパイラAPIとして提供していたとしても,CPUが対応しておらず,SIMD命令を利用できない場合もあるからだ. また,AVXやSSEの切り替えが容易になり,ベンチマークテストがしやすいという利点もあるだろう.

さて,具体的には以下のようにオプションを指定してコンパイルするとよい. gccの場合は,

有効化する機能 コマンド
AVX-512 $ g++ -std=gnu++11 -march=native -mavx512f -DENABLE_AVX main.cpp -o main.o
AVX $ g++ -std=gnu++11 -march=native -DENABLE_AVX main.cpp -o main.o
SSE $ g++ -std=gnu++11 -march=native -DENABLE_SSE main.cpp -o main.o
ARM NEON $ g++ -std=gnu++11 -mfpu=neon-fp-armv8 -DENABLE_NEON main.cpp -o main.o
SIMDを利用しない $ g++ -std=gnu++11 main.cpp -o main.o

であり,MSVCの場合は,

有効化する機能 コマンド
AVX > cl.exe /arch:AVX /DENABLE_AVX main.cpp
SSE > cl.exe /arch:SSE2 /DENABLE_SSE main.cpp
SIMDを利用しない > cl.exe main.cpp

といった具合である.

ベクトルの内積計算

定番のベクトルの内積を計算するコードを示す. FMA(積和演算)が利用可能な場合は,そちらを用いて,高速に処理できるようにしてある. また,SIMDを用いない場合であっても,C++11/C11以降で <cmath> から提供されている std::fma() を用いることで,内積計算の高速化が期待できるようにする.

#if defined(ENABLE_AVX512) && !defined(__AVX512F__)
#  error Macro: ENABLE_AVX512 is defined, but unable to use AVX512F intrinsic functions
#elif defined(ENABLE_AVX) && !defined(__AVX__)
#  error Macro: ENABLE_AVX is defined, but unable to use AVX intrinsic functions
#elif defined(ENABLE_SSE) && !defined(__SSE2__)
#  error Macro: ENABLE_SSE is defined, but unable to use SSE intrinsic functions
#elif defined(ENABLE_NEON) && !defined(__ARM_NEON) && !defined(__ARM_NEON__)
#  error Macro: ENABLE_NEON is defined, but unable to use NEON intrinsic functions
#else


#include <cmath>
#include <cstddef>
#include <algorithm>
#include <iostream>
#include <memory>
#include <type_traits>
#if defined(_MSC_VER) || defined(__MINGW32__)
#  include <malloc.h>
#else
#  include <cstdlib>
#endif
#if defined(ENABLE_AVX512) || defined(ENABLE_AVX) || defined(ENABLE_SSE)
#  ifdef _MSC_VER
#    include <intrin.h>
#  else
#    include <x86intrin.h>
#  endif  // _MSC_VER
#elif defined(ENABLE_NEON)
#  include <arm_neon.h>
#endif  // defined(ENABLE_AVX512) || defined(ENABLE_AVX) || defined(ENABLE_SSE)

#if defined(__cplusplus) && __cplusplus < 201103L
#  ifdef _MSC_VER
#    define alignof(n)  __alignof(n)
#  else
#    define alignof(n)  __alignof__(n)
#  endif  // _MSC_VER
#endif  // defined(__cplusplus) || cplusplus < 201103L


/*!
 * @brief アラインメントされたメモリを動的確保する関数
 * @param [in] size       確保するメモリサイズ (単位はbyte)
 * @param [in] alignment  アラインメント (2のべき乗を指定すること)
 * @return  アラインメントし,動的確保されたメモリ領域へのポインタ
 */
template<typename T = void*, typename std::enable_if<std::is_pointer<T>::value, std::nullptr_t>::type = nullptr>
static inline T
alignedMalloc(std::size_t size, std::size_t alignment) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  return reinterpret_cast<T>(_aligned_malloc(size, alignment));
#else
  void* p;
  return reinterpret_cast<T>(posix_memalign(&p, alignment, size) == 0 ? p : nullptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


/*!
 * @brief アラインメントされたメモリを解放する関数
 * @param [in] ptr  解放対象のメモリの先頭番地を指すポインタ
 */
static inline void
alignedFree(void* ptr) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  _aligned_free(ptr);
#else
  std::free(ptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


/*!
 * @brief std::unique_ptr で利用するアラインされたメモリ用のカスタムデリータ
 */
struct AlignedDeleter
{
  void
  operator()(void* p) const noexcept
  {
    alignedFree(p);
  }
};


#if defined(ENABLE_AVX512)
static constexpr int ALIGN = alignof(__m512);
#elif defined(ENABLE_AVX)
static constexpr int ALIGN = alignof(__m256);
#elif defined(ENABLE_SSE)
static constexpr int ALIGN = alignof(__m128);
#elif defined(ENABLE_NEON)
static constexpr int ALIGN = alignof(float32x4_t);
#else
static constexpr int ALIGN = 8;
#endif  // defined(ENABLE_AVX512)


/*!
 * @brief 内積計算を行う関数
 * @param [in] a  ベクトルその1
 * @param [in] b  ベクトルその2
 * @param [in] n  ベクトルのサイズ
 * @return  内積
 */
static inline float
innerProduct(const float* a, const float* b, std::size_t n)
{
#if defined(ENABLE_AVX512)
  static constexpr std::size_t INTERVAL = sizeof(__m512) / sizeof(float);
  __m512 sumx16 = {0};
  for (std::size_t i = 0; i < n; i += INTERVAL) {
    __m512 ax16 = _mm512_load_ps(&a[i]);
    __m512 bx16 = _mm512_load_ps(&b[i]);
#  ifdef __FMA__
    sumx16 = _mm512_fmadd_ps(ax16, bx16, sumx16);
#  else
    sumx16 = _mm512_add_ps(sumx16, _mm512_mul_ps(ax16, bx16));
#  endif  // __FMA__
  }

  alignas(ALIGN) float s[INTERVAL] = {0};
  _mm512_store_ps(s, sumx16);

  std::size_t offset = n - n % INTERVAL;
  return std::inner_product(
      a + offset,
      a + n,
      b + offset,
      std::accumulate(std::begin(s), std::end(s), 0.0f));
#elif defined(ENABLE_AVX)
  static constexpr std::size_t INTERVAL = sizeof(__m256) / sizeof(float);
  __m256 sumx8 = {0};
  for (std::size_t i = 0; i < n; i += INTERVAL) {
    __m256 ax8 = _mm256_load_ps(&a[i]);
    __m256 bx8 = _mm256_load_ps(&b[i]);
#  ifdef __FMA__
    sumx8 = _mm256_fmadd_ps(ax8, bx8, sumx8);
#  else
    sumx8 = _mm256_add_ps(sumx8, _mm256_mul_ps(ax8, bx8));
#  endif  // __FMA__
  }

  alignas(ALIGN) float s[INTERVAL] = {0};
  _mm256_store_ps(s, sumx8);

  std::size_t offset = n - n % INTERVAL;
  return std::inner_product(
      a + offset,
      a + n,
      b + offset,
      std::accumulate(std::begin(s), std::end(s), 0.0f));
#elif defined(ENABLE_SSE)
  static constexpr std::size_t INTERVAL = sizeof(__m128) / sizeof(float);
  __m128 sumx4 = {0};
  for (std::size_t i = 0; i < n; i += INTERVAL) {
    __m128 ax4 = _mm_load_ps(&a[i]);
    __m128 bx4 = _mm_load_ps(&b[i]);
#  ifdef __FMA__
    sumx4 = _mm_fmadd_ps(ax4, bx4, sumx4);
#  else
    sumx4 = _mm_add_ps(sumx4, _mm_mul_ps(ax4, bx4));
#  endif  // __FMA__
  }

  alignas(ALIGN) float s[INTERVAL] = {0};
  _mm_store_ps(s, sumx4);
  float sum = std::accumulate(std::begin(s), std::end(s), 0.0f);

  std::size_t offset = n - n % INTERVAL;
  return std::inner_product(
      a + offset,
      a + n,
      b + offset,
      std::accumulate(std::begin(s), std::end(s), 0.0f));
#elif defined(ENABLE_NEON)
  static constexpr std::size_t INTERVAL = sizeof(float32x4_t) / sizeof(float);
  float32x4_t sumx4 = {0};
  for (std::size_t i = 0; i < n; i += INTERVAL) {
    float32x4_t ax4 = vld1q_f32(&a[i]);
    float32x4_t bx4 = vld1q_f32(&b[i]);
    sumx4 = vmlaq_f32(sumx4, ax4, bx4);
  }

  std::size_t offset = n - n % INTERVAL;
  return std::inner_product(
      a + offset,
      a + n,
      b + offset,
      std::accumulate(std::begin(s), std::end(s), 0.0f));
#else
  float sum = 0.0f;
  for (std::size_t i = 0; i < n; i++) {
    // <cmath>のstd::fma関数を用いると,積和演算がハードウェアのサポートを受けることを期待できる
    // 処理としては, sum += a[i] * b[i]; と同じ
    sum = std::fma(a[i], b[i], sum);
  }
  return sum;
#endif  // defined(ENABLE_AVX512)
}


int
main()
{
  static constexpr int N_ELEMENT = 256;

  std::unique_ptr<float[], AlignedDeleter> a(alignedMalloc<float*>(N_ELEMENT * sizeof(float), ALIGN));
  std::unique_ptr<float[], AlignedDeleter> b(alignedMalloc<float*>(N_ELEMENT * sizeof(float), ALIGN));
  for (int i = 0; i < N_ELEMENT; i++) {
    a[i] = static_cast<float>(i);
    b[i] = static_cast<float>(i);
  }
  std::cout << innerProduct(a.get(), b.get(), N_ELEMENT) << std::endl;

  return 0;
}


#endif  // defined(ENABLE_AVX512) && !defined(__AVX512F__)

最近傍法による画像の2倍拡大

最近傍法,すなわち単純なピクセルコピーのみを行って,8bitグレースケール画像を2倍に拡大するコードを記述する. 2倍拡大という条件に限定すれば,出力先画像のインデックス値のとる値が単純になるので,SIMDで簡単に処理を記述できる.

読み込む画像ファイル名は test.jpg とし,読み込みにOpenCVを用いる. 画像ファイルの横幅は,16または32の倍数でなければならない.

コンパイルは以下のようにするとよい.

$ g++ -std=gnu++11 main.cpp -march=native -DENABLE_AVX -I/usr/include/opencv -I/usr/include/opencv2 -lopencv_core -lopencv_highgui -lopencv_imgcodecs -o main.o

AVX-512を利用する場合は, -mavx512vbmi -DENABLE_AVX512 を付加するとよい.

なお,OpenCVcv::Matカスタムアロケータを適用することができるらしいが,コードが煩雑になりそうなので,SSE/AVXにおいてはアラインメント条件を満たさなくてもよい関数を用いている.

#if defined(ENABLE_AVX512) && !defined(__AVX512F__)
#  error Macro: ENABLE_AVX512 is defined, but unable to use AVX512F intrinsic functions
#elif defined(ENABLE_AVX) && !defined(__AVX__)
#  error Macro: ENABLE_AVX is defined, but unable to use AVX intrinsic functions
#elif defined(ENABLE_SSE) && !defined(__SSE2__)
#  error Macro: ENABLE_SSE is defined, but unable to use SSE intrinsic functions
#elif defined(ENABLE_NEON) && !defined(__ARM_NEON) && !defined(__ARM_NEON__)
#  error Macro: ENABLE_NEON is defined, but unable to use NEON intrinsic functions
#else  // defined(ENABLE_AVX512) && !defined(__AVX512F__)


#include <cmath>
#include <cstddef>
#include <iostream>
#include <memory>
#include <type_traits>
#if defined(_MSC_VER) || defined(__MINGW32__)
#  include <malloc.h>
#else
#  include <cstdlib>
#endif
#if defined(ENABLE_AVX512) || defined(ENABLE_AVX) || defined(ENABLE_SSE)
#  ifdef _MSC_VER
#    include <intrin.h>
#  else
#    include <x86intrin.h>
#  endif  // _MSC_VER
#elif defined(ENABLE_NEON)
#  include <arm_neon.h>
#endif  // defined(ENABLE_AVX512) || defined(ENABLE_AVX) || defined(ENABLE_SSE)

#include <opencv2/opencv.hpp>


#if defined(_MSC_VER) && _MSC_VER >= 1400 || \
  defined(__GNUC__) && defined(__GNUC_MINOR__) && (__GNUC__ > 2 || __GNUC__ == 2 && __GNUC_MINOR__ >= 92)
#  define restrict  __restrict
#else
#  define restrict
#endif


/*!
 * @brief アラインメントされたメモリを動的確保する関数
 * @param [in] size       確保するメモリサイズ (単位はbyte)
 * @param [in] alignment  アラインメント (2のべき乗を指定すること)
 * @return  アラインメントし,動的確保されたメモリ領域へのポインタ
 */
template<typename T = void*, typename std::enable_if<std::is_pointer<T>::value, std::nullptr_t>::type = nullptr>
static inline T
alignedMalloc(std::size_t size, std::size_t alignment) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  return reinterpret_cast<T>(_aligned_malloc(size, alignment));
#else
  void* p;
  return reinterpret_cast<T>(posix_memalign(&p, alignment, size) == 0 ? p : nullptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


/*!
 * @brief アラインメントされたメモリを解放する関数
 * @param [in] ptr  解放対象のメモリの先頭番地を指すポインタ
 */
static inline void
alignedFree(void* ptr) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  _aligned_free(ptr);
#else
  std::free(ptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


#if defined(ENABLE_AVX512)
static constexpr int ALIGN = alignof(__m512i);
#elif defined(ENABLE_AVX)
static constexpr int ALIGN = alignof(__m256i);
#elif defined(ENABLE_SSE)
static constexpr int ALIGN = alignof(__m128i);
#elif defined(ENABLE_NEON)
static constexpr int ALIGN = alignof(uint8x16_t);
#else
static constexpr int ALIGN = 8;
#endif  // defined(ENABLE_AVX512)


/*!
 * @brief 入力画像データを最近傍法により,2倍のサイズに拡大する
 * @param [out] dstImageData  出力画像データ領域の先頭へのポインタ
 * @param [in]  dstWidth      出力画像データの横幅
 * @param [in]  dstHeight     出力画像データの縦幅
 * @param [in]  srcImageData  入力画像データ領域の先頭へのポインタ
 * @param [in]  srcWidth      入力画像データの横幅
 * @param [in]  srcHeight     入力画像データの縦幅
 * @return  アラインメントし,動的確保されたメモリ領域へのポインタ
 */
static inline void
scale2x(
    unsigned char* restrict dstImageData,
    int dstWidth,
    int dstHeight,
    const unsigned char* restrict srcImageData,
    int srcWidth,
    int srcHeight) noexcept
{
  static constexpr int X_RATIO = 2;
  static constexpr int Y_RATIO = 2;
#if defined(ENABLE_AVX512)
  static constexpr int INTERVAL = sizeof(__m512i) / sizeof(unsigned char);
  static const __m512i LOWIDX = _mm512_setr_epi64(
      0x4303420241014000,
      0x4707460645054404,
      0x4b0b4a0a49094808,
      0x4f0f4e0e4d0d4c0c,
      0x5313521251115010,
      0x5717561655155414,
      0x5b1b5a1a59195818,
      0x5f1f5e1e5d1d5c1c);
  static const __m512i HIGHIDX = _mm512_setr_epi64(
      0x6323622261216020,
      0x6727662665256424,
      0x6b2b6a2a69296828,
      0x6f2f6e2e6d2d6c2c,
      0x7333723271317030,
      0x7737763675357434,
      0x7b3b7a3a79397838,
      0x7f3f7e3e7d3d7c3c);
#elif defined(ENABLE_AVX)
  static constexpr int INTERVAL = sizeof(__m256i) / sizeof(unsigned char);
#elif defined(ENABLE_SSE)
  static constexpr int INTERVAL = sizeof(__m128i) / sizeof(unsigned char);
#elif defined(ENABLE_NEON)
  static constexpr int INTERVAL = sizeof(uint8x16_t) / sizeof(unsigned char);
#else
  static constexpr int INTERVAL = sizeof(unsigned char);
#endif  // defined(ENABLE_AVX512)

  for (int i = 0; i < dstHeight; i++) {
    for (int j = 0; j < dstWidth; j += INTERVAL * X_RATIO) {
#if defined(ENABLE_AVX512)
      // 64pixel分の画素データをロード
      __m512i v512 = _mm512_loadu_si512(reinterpret_cast<const __m512i*>(&srcImageData[i / Y_RATIO * srcWidth + j / X_RATIO]));
      // インタリーブ
      __m512i v512l = _mm512_permutex2var_epi8(v512, LOWIDX, v512);
      __m512i v512u = _mm512_permutex2var_epi8(v512, HIGHIDX, v512);
      // 64pixel x 2のデータを書き込み
      _mm512_storeu_si512(reinterpret_cast<__m512i*>(&dstImageData[i * dstWidth + j + sizeof(__m512i) * 0]), v512l);
      _mm512_storeu_si512(reinterpret_cast<__m512i*>(&dstImageData[i * dstWidth + j + sizeof(__m512i) * 1]), v512u);
#elif defined(ENABLE_AVX)
      // 32pixel分の画素データをロード
      __m256i v256 = _mm256_loadu_si256(reinterpret_cast<const __m256i*>(&srcImageData[i / Y_RATIO * srcWidth + j / X_RATIO]));
      // インタリーブ
      __m256i v256l_ = _mm256_unpacklo_epi8(v256, v256);
      __m256i v256u_ = _mm256_unpackhi_epi8(v256, v256);
      // 上下128bit交換
      __m256i v256l = _mm256_permute2f128_si256(v256l_, v256u_, 0x20);
      __m256i v256u = _mm256_permute2f128_si256(v256l_, v256u_, 0x31);
      // 32pixel x 2のデータを書き込み
      _mm256_storeu_si256(reinterpret_cast<__m256i*>(&dstImageData[i * dstWidth + j + sizeof(__m256i) * 0]), v256l);
      _mm256_storeu_si256(reinterpret_cast<__m256i*>(&dstImageData[i * dstWidth + j + sizeof(__m256i) * 1]), v256u);
#elif defined(ENABLE_SSE)
      // 16pixel分の画素データをロード
      __m128i v128 = _mm_loadu_si128(reinterpret_cast<const __m128i*>(&srcImageData[i / Y_RATIO * srcWidth + j / X_RATIO]));
      // インタリーブ
      __m128i v128l = _mm_unpacklo_epi8(v128, v128);
      __m128i v128u = _mm_unpackhi_epi8(v128, v128);
      // 16pixel x 2のデータを書き込み
      _mm_storeu_si128(reinterpret_cast<__m128i*>(&dstImageData[i * dstWidth + j + sizeof(__m128i) * 0]), v128l);
      _mm_storeu_si128(reinterpret_cast<__m128i*>(&dstImageData[i * dstWidth + j + sizeof(__m128i) * 1]), v128u);
#elif defined(ENABLE_NEON)
      // 16pixel分の画素データをロード
      uint8x16_t v128 = vld1q_u8(&srcImageData[i / Y_RATIO * srcWidth + j]);
      // インタリーブ
      uint8x16x2_t v128x2 = vzipq_u8(v128, v128);
      // 16pixel x 2のデータを書き込み
      vst1q_u8(dstImageData[i * dstWidth + j + sizeof(uint8x16_t) * 0], v128x2.val[0]);
      vst1q_u8(dstImageData[i * dstWidth + j + sizeof(uint8x16_t) * 1], v128x2.val[1]);
#else
      dstImageData[i * dstWidth + j] = srcImageData[i / Y_RATIO * srcWidth + j];
#endif  // defined(ENABLE_AVX512)
    }
  }
}


int
main()
{
  cv::Mat img = cv::imread("test.jpg", 0);
  if (img.data == nullptr) {
    std::cerr << "Cannot open image file: test.jpg" << std::endl;
    return 1;
  }
  cv::Mat scaledImg(cv::Size(img.cols * 2, img.rows * 2), CV_8UC1);
  scale2x(scaledImg.data, scaledImg.cols, scaledImg.rows, img.data, img.cols, img.rows);

  cv::namedWindow("src", CV_WINDOW_AUTOSIZE);
  cv::namedWindow("scaled", CV_WINDOW_AUTOSIZE);
  cv::imshow("src", img);
  cv::imshow("scaled", scaledImg);
  std::cout << "Please hit any key on the window to exit this program" << std::endl;
  cv::waitKey(0);

  return 0;
}


#endif  // defined(ENABLE_AVX512) && !defined(__AVX512F__)

何の脈略も無しに,SSE/AVXやARM NEONの組み込み関数や型を利用したが,SSE/AVXに関してはIntelIntrinsics Guideを,ARM NEONに関してはARM NEON Intrinsicsを参照するとよい.

SSE/AVXの変数型は以下の通り.

内容
__m128 float 型4個分
__m128d double 型2個分
__m128i 整数型 (intunsigned char などを格納できる)
__m256 float 型8個分
__m256d double 型4個分
__m256i 整数型 (intunsigned char などを格納できる)
__m512 float 型16個分
__m512d double 型8個分
__m512i 整数型 (intunsigned char などを格納できる)

SSE/AVXの組み込み関数は基本的に

  • SSEの場合, _mm_[xxx]{[u]}_[yyy]
  • AVXの場合, _mm256_[xxx]{[u]}_[yyy]
  • AVX-512の場合, _mm512_[xxx]{[u]}_[yyy]

の形式で命名されている. [xxx][{u}][yyy] の部分については以下の通り.

該当部分 内容
[xxx] loadstore など,行いたい命令がここにくる
[u] u が付いている関数はアラインメント条件を満たしていなくても,SEGVで落ちない
[yyy] 引数の型によって変化する. ps なら __m128pd なら __m128dsi128 なら __m128i

pspd はそれぞれ Precision Single, Precision Double の略であるそうだ. (si は調べていない)

ARM NEONの変数型は見た目通り, [xxx][size]x[NNN]{x[MMM]} の形式となっている.

該当部分 内容
[xxx] uintintfloat などのベクタの1要素の型がここにくる
[size] ベクタの要素型1つのサイズ (単位はbit)
[NNN] ベクタ要素の個数
[MMM] インタリーブ用にくっつけたNEONレジスタの個数.2から4までの値ろ取り,1つの場合は省略される

ARM NEONの組み込み関数も直感的に利用できる命名で, v[xxx]{[q]}_{yyy} となっている.

該当部分 内容
[xxx] addld など,行いたい命令がここにくる
[q] qが付いていればQレジスタ(128bit)を用いる命令,付いていないならばDレジスタ(64bit)を用いる命令
[yyy] 引数の型によって変化する. u8s16f32 など

まとめ

この記事では以下のことを紹介した.

  • SIMDの概要
  • SIMDの組み込み関数の利用方法
  • コンパイラの差を吸収するアラインメントの指定方法
  • ベクトルの内積を計算するサンプルコード

特に,SIMDの組み込み関数の利用方法を簡単にまとめると以下のようになる.

  • alignas(alignof(__m256i)) ... の形で,変数のアラインメント指定
    • 古いMSVCなら __declspec(align(32))
    • 古いgccなら __attribute__((aligned(32)))
  • gcc
    • #include <x86intrin>
    • $ g++ -march=native ...
    • pisix_memalign() でアラインされた動的メモリ確保, std::free() で解放
  • MSVC
    • #include <intrin>
    • > cl.exe /arch:AVX2 ...
    • _aligned_malloc() でアラインされた動的メモリ確保, _aligned_free() で解放
  • AVX-512非対応のCPUでAVX-512をテストする場合は,Intelのエミュレータを利用

この記事はあくまでSIMDの基礎に過ぎないが,あとは組み込み関数を調べ,うまく組み合わせることで,SIMDをプログラムに組み込めるようになるかもしれない.

参考文献

広告を非表示にする