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

はじめに

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

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

項目	利用可能レジスタ
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 追記

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

• AVXとAVX-512のインタリーブ
• std::vectorなどでアラインされた領域を用いる

2019/02/03 追記

実行時にSSE/AVX等のx86/x64の命令が利用可能であるかをcpuidを用いて判断する方法について追記した． また，この記事中のユーティリティ関数をまとめたシングルヘッダファイルをkoturn/SimdUtilにて公開している．

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ではヘッダをインクルードするだけではSIMDの組み込み関数は利用できないため，以下のようにコンパイルオプションを指定する必要がある． 一方，MSVCはオプション指定をしなくてもSIMDの組み込み関数を利用できる．

なお，全てのx64プロセッサではSSE2までは利用できるため，gccであってもx64バイナリを生成するのであれば， -msse2 といったオプションの指定無しにSSE2までの組み込み関数が利用できるようだ．

gccの場合，コンパイラの自動ベクトル化でどの命令を利用するかの許可と利用可能な組み込み関数の許可を兼ねているのに対し，MSVCは自動ベクトル化でどの命令を利用するかの許可のみである． x86/x64においては，後述するcpuidによる実行時の利用可能なSIMD命令の判定が可能なため，MSVCの方が融通が利くように思われる．

命令セット	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のみ．確保サイズはアラインメントの整数倍に限る．C11/C++17の標準ライブラリ関数
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 アラインメントされたメモリを動的確保する関数
 * @tparam T  確保するメモリの要素型．この関数の返却値はT*
 * @param [in] nBytes     確保するメモリサイズ (単位はbyte)
 * @param [in] alignment  アラインメント (2のべき乗を指定すること)
 * @return  アラインメントし，動的確保されたメモリ領域へのポインタ
 */
template<typename T = void>
static inline T*
alignedMalloc(std::size_t nBytes, std::size_t alignment = alignof(T)) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  return reinterpret_cast<T*>(::_aligned_malloc(nBytes, alignment));
#else
  void* p;
  return reinterpret_cast<T*>(::posix_memalign(&p, alignment, nBytes) == 0 ? p : nullptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


/*!
 * @brief アラインメントされたメモリを動的確保する関数．配列向けにalignedMallocの引数指定が簡略化されている
 * @tparam T  確保する配列の要素型．この関数の返却値はT*
 * @param [in] size       確保する要素数．すなわち確保するサイズは size * sizeof(T)
 * @param [in] alignment  アラインメント (2のべき乗を指定すること)
 * @return  アラインメントし，動的確保されたメモリ領域へのポインタ
 */
template<typename T>
static inline T*
alignedAllocArray(std::size_t size, std::size_t alignment = alignof(T)) noexcept
{
  return alignedMalloc<T>(size * sizeof(T), alignment);
}


/*!
 * @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(alignedAllocArray<unsigned char>(10, 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)


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


/*!
 * @brief アラインメントされたメモリを動的確保する関数．配列向けにalignedMallocの引数指定が簡略化されている
 * @tparam T  確保する配列の要素型．この関数の返却値はT*
 * @param [in] size       確保する要素数．すなわち確保するサイズは size * sizeof(T)
 * @param [in] alignment  アラインメント (2のべき乗を指定すること)
 * @return  アラインメントし，動的確保されたメモリ領域へのポインタ
 */
template<typename T>
static inline T*
alignedAllocArray(std::size_t size, std::size_t alignment = alignof(T)) noexcept
{
  return alignedMalloc<T>(size * sizeof(T), alignment);
}


/*!
 * @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(alignedAllocArray<float>(N_ELEMENT, ALIGN));
  std::unique_ptr<float[], AlignedDeleter> b(alignedAllocArray<float>(N_ELEMENT, 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 を付加するとよい．

なお，OpenCVの cv::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 アラインメントされたメモリを動的確保する関数
 * @tparam T  確保するメモリの要素型．この関数の返却値はT*
 * @param [in] nBytes     確保するメモリサイズ (単位はbyte)
 * @param [in] alignment  アラインメント (2のべき乗を指定すること)
 * @return  アラインメントし，動的確保されたメモリ領域へのポインタ
 */
template<typename T = void>
static inline T*
alignedMalloc(std::size_t nBytes, std::size_t alignment = alignof(T)) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  return reinterpret_cast<T*>(::_aligned_malloc(nBytes, alignment));
#else
  void* p;
  return reinterpret_cast<T*>(::posix_memalign(&p, alignment, nBytes) == 0 ? p : nullptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


/*!
 * @brief アラインメントされたメモリを動的確保する関数．配列向けにalignedMallocの引数指定が簡略化されている
 * @tparam T  確保する配列の要素型．この関数の返却値はT*
 * @param [in] size       確保する要素数．すなわち確保するサイズは size * sizeof(T)
 * @param [in] alignment  アラインメント (2のべき乗を指定すること)
 * @return  アラインメントし，動的確保されたメモリ領域へのポインタ
 */
template<typename T>
static inline T*
alignedAllocArray(std::size_t size, std::size_t alignment = alignof(T)) noexcept
{
  return alignedMalloc<T>(size * sizeof(T), alignment);
}


/*!
 * @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に関してはIntelのIntrinsics Guideを，ARM NEONに関してはARM NEON Intrinsicsを参照するとよい．

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

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

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

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

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

該当部分	内容
[xxx]	load や store など，行いたい命令がここにくる
[u]	u が付いている関数はアラインメント条件を満たしていなくても，SEGVで落ちない
[yyy]	引数の型によって変化する． ps なら __m128 ， pd なら __m128d ， si128 なら __m128i

ps ， pd はそれぞれ Precision Single， Precision Double の略であるそうだ． （si は調べていない）

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

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

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

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

SSEやAVXが実行時に利用可能かどうかを調べる

利用可能かどうかを調べるモチベーション

ここまではコンパイル時にどの命令を使用するかを指定することを前提にしていた． しかし，実行時にSIMD命令が利用可能かどうかを調べたい場合がある．

Linuxであれば，基本的にプログラムはその環境でコンパイルし，実行することが多いため，実行時にSIMD命令が利用可能かどうかを調べなくてもよいが，Windowsにおいてはある環境でコンパイルしたプログラムを様々な環境で動作させることが多いため，利用可否を調べる必要がある．

ここでは，SSE/AVX等のx86/x64におけるSIMD命令が利用可能かどうかを調べる方法を示す． （ARMのNEONについては未調査）

cpuid命令とcpuidの組み込み関数

答えは簡単でcpuid命令を利用するとよい． この命令はアセンブラでは1命令として用意されている．

mov $1,%eax  ; cpuidの引数1
mov $0,%ecx  ; cpuidの引数2
cpuid  ; これでeax, ebx, ecx, edxに結果が格納される

まず eax および ecx に取得したいCPUの情報に関する値をセットし，その後cpuid命令を実行すると，eax, ebx, ecx, edxに情報が返却される命令となっている．

アセンブラ，およびインラインアセンブラでなければ利用できないのかというとそうではなく，gcc, clang, MSVCであれば，cpuidの組み込み関数が用意されている． しかし，gcc/clangとMSVCで引数等が異なるため，以下のように統一して利用できるインライン関数を用意すると楽である．

#include <array>
#include <type_traits>

#if defined(__GNUC__)
#  include <cpuid.h>
#elif defined(_MSC_VER)
#  include <intrin.h>
#endif

/*!
 * @brief cpuidの実行結果を第一引数に格納する
 *
 * 実行結果のeaxをcpuInfo[0]，ebxをcpuInfo[1]，ecxをcpuInfo[2]，edxをcpuInfo[3]にコピーする
 *
 * @tparam T  int*
 * @param [out] cpuInfo  cpuidの結果格納先．cpuInfo[0]からcpuInfo[3]に結果が格納される．
 * @param [in]  eax  cpuidの引数
 */
template<
  typename T,
  typename std::enable_if<std::is_same<T, int*>::value, std::nullptr_t>::type = nullptr
>
static inline void
cpuid(T cpuInfo, int eax) noexcept
{
#if defined(__GNUC__)
  ::__cpuid(eax, cpuInfo[0], cpuInfo[1], cpuInfo[2], cpuInfo[3]);
#elif defined(_MSC_VER)
  ::__cpuid(cpuInfo, eax);
#endif  // defined(__GNUC__)
}


/*!
 * @brief cpuidの実行結果を第一引数の配列に格納する
 *
 * 実行結果のeaxをcpuInfo[0]，ebxをcpuInfo[1]，ecxをcpuInfo[2]，edxをcpuInfo[3]にコピーする
 *
 * @tparam kSize  配列サイズ
 * @param [out] cpuInfo  cpuidの結果格納先配列．要素数が4以上でなければコンパイルエラーとなる
 * @param [in]  eax  cpuidの引数
 */
template<std::size_t kSize>
static inline void
cpuid(int (&cpuInfo)[kSize], int eax) noexcept
{
  static_assert(kSize >= 4, "CPU info array size must be four or more");

  cpuid(&cpuInfo[0], eax);
}


/*!
 * @brief cpuidの実行結果を第一引数のstd::arrayに格納する
 *
 * 実行結果のeaxをcpuInfo[0]，ebxをcpuInfo[1]，ecxをcpuInfo[2]，edxをcpuInfo[3]にコピーする
 *
 * @tparam kSize  配列サイズ
 * @param [out] cpuInfo  cpuidの結果格納先配列．要素数が4以上でなければコンパイルエラーとなる
 * @param [in]  eax  cpuidの引数
 */
template<std::size_t kSize>
static inline void
cpuid(std::array<int, kSize>& cpuInfo, int eax) noexcept
{
  static_assert(kSize >= 4, "CPU info array size must be four or more");

  cpuid(cpuInfo.data(), eax);
}


/*!
 * @brief cpuidの実行結果を第一引数に格納する
 *
 * 実行結果のeaxをcpuInfo[0]，ebxをcpuInfo[1]，ecxをcpuInfo[2]，edxをcpuInfo[3]にコピーする
 *
 * @tparam T  int*
 * @param [out] cpuInfo  cpuidの結果格納先．cpuInfo[0]からcpuInfo[3]に結果が格納される．
 * @param [in]  eax  cpuidの引数
 * @param [in]  ecx  cpuidの引数
 */
template<
  typename T,
  typename std::enable_if<std::is_same<T, int*>::value, std::nullptr_t>::type = nullptr
>
static inline void
cpuidex(T cpuInfo, int eax, int ecx) noexcept
{
#if defined(__GNUC__)
  ::__cpuid_count(eax, ecx, cpuInfo[0], cpuInfo[1], cpuInfo[2], cpuInfo[3]);
#elif defined(_MSC_VER)
  ::__cpuidex(cpuInfo, eax, ecx);
#endif  // defined(__GNUC__)
}


/*!
 * @brief cpuidの実行結果を第一引数の配列に格納する
 *
 * 実行結果のeaxをcpuInfo[0]，ebxをcpuInfo[1]，ecxをcpuInfo[2]，edxをcpuInfo[3]にコピーする
 *
 * @tparam kSize  配列サイズ
 * @param [out] cpuInfo  cpuidの結果格納先配列．要素数が4以上でなければコンパイルエラーとなる
 * @param [in]  eax  cpuidの引数
 * @param [in]  ecx  cpuidの引数
 */
template<std::size_t kSize>
static inline void
cpuidex(int (&cpuInfo)[kSize], int eax, int ecx) noexcept
{
  static_assert(kSize >= 4, "[util::cpuidex] CPU info array size must be four or more");

  cpuidex(&cpuInfo[0], eax, ecx);
}


/*!
 * @brief cpuidの実行結果を第一引数のstd::arrayに格納する
 *
 * 実行結果のeaxをcpuInfo[0]，ebxをcpuInfo[1]，ecxをcpuInfo[2]，edxをcpuInfo[3]にコピーする
 *
 * @tparam kSize  配列サイズ
 * @param [out] cpuInfo  cpuidの結果格納先配列．要素数が4以上でなければコンパイルエラーとなる
 * @param [in]  eax  cpuidの引数
 * @param [in]  ecx  cpuidの引数
 */
template<std::size_t kSize>
static inline void
cpuidex(std::array<int, kSize>& cpuInfo, int eax, int ecx) noexcept
{
  static_assert(kSize >= 4, "[util::cpuidex] CPU info array size must be four or more");

  cpuidex(cpuInfo.data(), eax, ecx);
}

cpuid() はeaxを指定し，ecxは0として，第一引数にeaxからedxの値を順に格納する関数， cpuidex() は cpuid() のecx指定版である．

上記は第一引数に配列や sd::array を放り込んだとき，サイズをコンパイル時に判定するようにしてある． C言語用に書き直すなら以下のような単純な形でよい．

#if defined(__GNUC__)
#  include <cpuid.h>
#elif defined(_MSC_VER)
#  include <intrin.h>
#endif

#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 cpuidの実行結果を第一引数に格納する
 *
 * 実行結果のeaxをcpuInfo[0]，ebxをcpuInfo[1]，ecxをcpuInfo[2]，edxをcpuInfo[3]にコピーする
 *
 * @param [out] cpuInfo  cpuidの結果格納先．cpuInfo[0]からcpuInfo[3]に結果が格納される．
 * @param [in]  eax  cpuidの引数
 */
static inline void
cpuid(int* cpuInfo, int eax)
{
#if defined(__GNUC__)
  __cpuid(eax, cpuInfo[0], cpuInfo[1], cpuInfo[2], cpuInfo[3]);
#elif defined(_MSC_VER)
  __cpuid(cpuInfo, eax);
#endif  // defined(__GNUC__)
}


/*!
 * @brief cpuidの実行結果を第一引数に格納する
 *
 * 実行結果のeaxをcpuInfo[0]，ebxをcpuInfo[1]，ecxをcpuInfo[2]，edxをcpuInfo[3]にコピーする
 *
 * @param [out] cpuInfo  cpuidの結果格納先．cpuInfo[0]からcpuInfo[3]に結果が格納される．
 * @param [in]  eax  cpuidの引数
 * @param [in]  ecx  cpuidの引数
 */
static inline void
cpuidex(int* cpuInfo, int eax, int ecx)
{
#if defined(__GNUC__)
  __cpuid_count(eax, ecx, cpuInfo[0], cpuInfo[1], cpuInfo[2], cpuInfo[3]);
#elif defined(_MSC_VER)
  __cpuidex(cpuInfo, eax, ecx);
#endif  // defined(__GNUC__)
}

これで，CPUの情報を取得する準備はできた．

cpuidから取得できる情報から，どうすればSIMD命令が利用できるか判定できるかはcpuidについてのドキュメント等を参照するとよいが，表にまとめると以下の通りである．

SIMD命令	引数eax	引数ecx	レジスタとフラグビット
MMX	1	0	edx [bit 23]
SSE	1	0	edx [bit 25]
SSE2	1	0	edx [bit 26]
SSE3	1	0	ecx [bit 0]
SSSE3	1	0	ecx [bit 9]
SSE4.1	1	0	ecx [bit 19]
SSE4.2	1	0	ecx [bit 20]
SSE4A	0x80000001	0	ecx [bit 6]
AVX	1	0	ecx [bit 28]
AVX2	7	0	ebx [bit 5]
FMA	1	0	ecx [bit 12]
AVX512F	7	0	ebx [bit 16]
AVX512BW	7	0	ebx [bit 30]
AVX512CD	7	0	ebx [bit 28]
AVX512DQ	7	0	ebx [bit 17]
AVX512ER	7	0	ebx [bit 27]
AVX512IFMA52	7	0	ebx [bit 21]
AVX512PF	7	0	ebx [bit 26]
AVX512VL	7	0	ebx [bit 31]
AVX512_4FMAPS	7	0	edx [bit 2]
AVX512_4VNNIW	7	0	edx [bit 3]
AVX512BITALG	7	0	ecx [bit 12]
AVX512VPOPCNTDQ	7	0	ecx [bit 14]
AVX512VBMI	7	0	ecx [bit 1]
AVX512VBMI2	7	0	ecx [bit 6]
AVX512VNNI	7	0	ecx [bit 11]

ちなみに，x64ではSSE，SSE2は利用可能であるとのことなので，わざわざ判定する必要はない．

以上を踏まえて，以下のようなインライン関数を定義したヘッダファイルを用意しておくと便利である． なお，名前空間を加える等，多少改良したものをGitHubに置いてある．

// cpuid.hpp
#ifndef CPUID_HPP
#define CPUID_HPP

#include <algorithm>
#include <array>
#include <string>
#include <type_traits>
#include <utility>

#if defined(__GNUC__)
#  include <cpuid.h>
#elif defined(_MSC_VER)
#  include <intrin.h>
#endif


/*!
 * @brief cpuidの実行結果を第一引数に格納する
 *
 * 実行結果のeaxをcpuInfo[0]，ebxをcpuInfo[1]，ecxをcpuInfo[2]，edxをcpuInfo[3]にコピーする
 *
 * @tparam T  int*
 * @param [out] cpuInfo  cpuidの結果格納先．cpuInfo[0]からcpuInfo[3]に結果が格納される．
 * @param [in]  eax  cpuidの引数
 */
template<
  typename T,
  typename std::enable_if<std::is_same<T, int*>::value, std::nullptr_t>::type = nullptr
>
static inline void
cpuid(T cpuInfo, int eax) noexcept
{
#if defined(__GNUC__)
  ::__cpuid(eax, cpuInfo[0], cpuInfo[1], cpuInfo[2], cpuInfo[3]);
#elif defined(_MSC_VER)
  ::__cpuid(cpuInfo, eax);
#endif  // defined(__GNUC__)
}


/*!
 * @brief cpuidの実行結果を第一引数の配列に格納する
 *
 * 実行結果のeaxをcpuInfo[0]，ebxをcpuInfo[1]，ecxをcpuInfo[2]，edxをcpuInfo[3]にコピーする
 *
 * @tparam kSize  配列サイズ
 * @param [out] cpuInfo  cpuidの結果格納先配列．要素数が4以上でなければコンパイルエラーとなる
 * @param [in]  eax  cpuidの引数
 */
template<std::size_t kSize>
static inline void
cpuid(int (&cpuInfo)[kSize], int eax) noexcept
{
  static_assert(kSize >= 4, "CPU info array size must be four or more");

  cpuid(&cpuInfo[0], eax);
}


/*!
 * @brief cpuidの実行結果を第一引数のstd::arrayに格納する
 *
 * 実行結果のeaxをcpuInfo[0]，ebxをcpuInfo[1]，ecxをcpuInfo[2]，edxをcpuInfo[3]にコピーする
 *
 * @tparam kSize  配列サイズ
 * @param [out] cpuInfo  cpuidの結果格納先配列．要素数が4以上でなければコンパイルエラーとなる
 * @param [in]  eax  cpuidの引数
 */
template<std::size_t kSize>
static inline void
cpuid(std::array<int, kSize>& cpuInfo, int eax) noexcept
{
  static_assert(kSize >= 4, "CPU info array size must be four or more");

  cpuid(cpuInfo.data(), eax);
}


/*!
 * @brief cpuidの実行結果を第一引数に格納する
 *
 * 実行結果のeaxをcpuInfo[0]，ebxをcpuInfo[1]，ecxをcpuInfo[2]，edxをcpuInfo[3]にコピーする
 *
 * @tparam T  int*
 * @param [out] cpuInfo  cpuidの結果格納先．cpuInfo[0]からcpuInfo[3]に結果が格納される．
 * @param [in]  eax  cpuidの引数
 * @param [in]  ecx  cpuidの引数
 */
template<
  typename T,
  typename std::enable_if<std::is_same<T, int*>::value, std::nullptr_t>::type = nullptr
>
static inline void
cpuidex(T cpuInfo, int eax, int ecx) noexcept
{
#if defined(__GNUC__)
  ::__cpuid_count(eax, ecx, cpuInfo[0], cpuInfo[1], cpuInfo[2], cpuInfo[3]);
#elif defined(_MSC_VER)
  ::__cpuidex(cpuInfo, eax, ecx);
#endif  // defined(__GNUC__)
}


/*!
 * @brief cpuidの実行結果を第一引数の配列に格納する
 *
 * 実行結果のeaxをcpuInfo[0]，ebxをcpuInfo[1]，ecxをcpuInfo[2]，edxをcpuInfo[3]にコピーする
 *
 * @tparam kSize  配列サイズ
 * @param [out] cpuInfo  cpuidの結果格納先配列．要素数が4以上でなければコンパイルエラーとなる
 * @param [in]  eax  cpuidの引数
 * @param [in]  ecx  cpuidの引数
 */
template<std::size_t kSize>
static inline void
cpuidex(int (&cpuInfo)[kSize], int eax, int ecx) noexcept
{
  static_assert(kSize >= 4, "[util::cpuidex] CPU info array size must be four or more");

  cpuidex(&cpuInfo[0], eax, ecx);
}


/*!
 * @brief cpuidの実行結果を第一引数のstd::arrayに格納する
 *
 * 実行結果のeaxをcpuInfo[0]，ebxをcpuInfo[1]，ecxをcpuInfo[2]，edxをcpuInfo[3]にコピーする
 *
 * @tparam kSize  配列サイズ
 * @param [out] cpuInfo  cpuidの結果格納先配列．要素数が4以上でなければコンパイルエラーとなる
 * @param [in]  eax  cpuidの引数
 * @param [in]  ecx  cpuidの引数
 */
template<std::size_t kSize>
static inline void
cpuidex(std::array<int, kSize>& cpuInfo, int eax, int ecx) noexcept
{
  static_assert(kSize >= 4, "[util::cpuidex] CPU info array size must be four or more");

  cpuidex(cpuInfo.data(), eax, ecx);
}


/*!
 * @brief cpuidの実行結果のうち，指定レジスタの指定ビットが立っているかどうか調べる
 * @param [in] eax    cpuidの引数
 * @param [in] index  cpuidの結果のインデックス．0ならeax，1ならebx，2ならecx，3ならedx
 * @param [in] nBit       立っているかどうか調べたいビット
 * @return 指定レジスタの指定ビットが立っているならtrue，そうでなければfalse
 */
static inline bool
cpuidBit(int eax, int index, int nBit) noexcept
{
    std::array<int, 4> cpuInfo;
    cpuid(cpuInfo, eax);
    return (cpuInfo[index] & (1 << nBit)) != 0;
}


/*!
 * @brief cpuidの実行結果のうち，指定レジスタの指定ビットが立っているかどうか調べる
 * @param [in] eax    cpuidの引数
 * @param [in] ecx    cpuidの引数
 * @param [in] index  cpuidの結果のインデックス．0ならeax，1ならebx，2ならecx，3ならedx
 * @param [in] nBit  立っているかどうか調べたいビット
 * @return 指定レジスタの指定ビットが立っているならtrue，そうでなければfalse
 */
static inline bool
cpuidexBit(int eax, int ecx, int index, int nBit) noexcept
{
  std::array<int, 4> cpuInfo;
  cpuidex(cpuInfo, eax, ecx);
  return (cpuInfo[index] & (1 << nBit)) != 0;
}


/*!
 * @brief MMX命令が利用可能かどうかを調べる．
 * @return MMX命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isMmxAvailable() noexcept
{
  return cpuidBit(1, 3, 23);
}

/*!
 * @brief SSE命令が利用可能かどうかを調べる．
 * @return SSE命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isSseAvailable() noexcept
{
  return cpuidBit(1, 3, 25);
}

/*!
 * @brief SSE2命令が利用可能かどうかを調べる．
 * @return SSE2命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isSse2Available() noexcept
{
  return cpuidBit(1, 3, 26);
}

/*!
 * @brief SSE3命令が利用可能かどうかを調べる．
 * @return SSE3命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isSse3Available() noexcept
{
  return cpuidBit(1, 2, 0);
}

/*!
 * @brief SSSE3命令が利用可能かどうかを調べる．
 * @return SSSE3命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isSsse3Available() noexcept
{
  return cpuidBit(1, 2, 9);
}


/*!
 * @brief SSE4.1命令が利用可能かどうかを調べる．
 * @return SSE4.1命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isSse41Available() noexcept
{
  return cpuidBit(1, 2, 19);
}

/*!
 * @brief SSE4.2命令が利用可能かどうかを調べる．
 * @return SSE4.2命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isSse42Available() noexcept
{
  return cpuidBit(1, 2, 20);
}

/*!
 * @brief SSE4A命令が利用可能かどうかを調べる．
 * @return SSE4A命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isSse4aAvailable() noexcept
{
  std::array<int, 4> cpuInfo;
  cpuid(cpuInfo, 0x80000000);
  if (static_cast<unsigned int>(cpuInfo[0]) < 0x80000001U) {
    return false;
  }
  return cpuidBit(0x80000001, 2, 6);
}

/*!
 * @brief AVX命令が利用可能かどうかを調べる．
 * @return AVX命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isAvxAvailable() noexcept
{
  return cpuidBit(1, 2, 28);
}

/*!
 * @brief AVX2命令が利用可能かどうかを調べる．
 * @return AVX2命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isAvx2Available() noexcept
{
  return cpuidBit(7, 1, 5);
}

/*!
 * @brief FMA命令が利用可能かどうかを調べる．
 * @return FMA命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isFmaAvailable() noexcept
{
  return cpuidBit(1, 2, 12);
}

/*!
 * @brief AVX512F命令が利用可能かどうかを調べる．
 * @return AVX512F命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isAvx512FAvailable() noexcept
{
  return cpuidBit(7, 1, 16);
}

/*!
 * @brief AVX512BW命令が利用可能かどうかを調べる．
 * @return AVX512BW命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isAvx512BwAvailable() noexcept
{
  return cpuidBit(7, 1, 30);
}

/*!
 * @brief AVX512CD命令が利用可能かどうかを調べる．
 * @return AVX512CD命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isAvx512CdAvailable() noexcept
{
  return cpuidBit(7, 1, 28);
}

/*!
 * @brief AVX512DQ命令が利用可能かどうかを調べる．
 * @return AVX512DQ命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isAvx512DqAvailable() noexcept
{
  return cpuidBit(7, 1, 17);
}

/*!
 * @brief AVX512ER命令が利用可能かどうかを調べる．
 * @return AVX512ER命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isAvx512ErAvailable() noexcept
{
  return cpuidBit(7, 1, 27);
}

/*!
 * @brief AVX512IFMA52命令が利用可能かどうかを調べる．
 * @return AVX512IFMA52命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isAvx512Ifma52Available() noexcept
{
  return cpuidBit(7, 1, 21);
}

/*!
 * @brief AVX512PF命令が利用可能かどうかを調べる．
 * @return AVX512PF命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isAvx512PfAvailable() noexcept
{
  return cpuidBit(7, 1, 26);
}

/*!
 * @brief AVX512VL命令が利用可能かどうかを調べる．
 * @return AVX512VL命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isAvx512VlAvailable() noexcept
{
  return cpuidBit(7, 1, 31);
}

/*!
 * @brief AVX512_4FMAPS命令が利用可能かどうかを調べる．
 * @return AVX512_4FMAPS命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isAvx512_4fmapsAvailable() noexcept
{
  return cpuidBit(7, 3, 2);
}

/*!
 * @brief AVX512_4VNNIW命令が利用可能かどうかを調べる．
 * @return AVX512_4VNNIW命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isAvx512_4vnniwAvailable() noexcept
{
  return cpuidBit(7, 3, 3);
}

/*!
 * @brief AVX512BITALG命令が利用可能かどうかを調べる．
 * @return AVX512BITALG命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isAvx512BitalgAvailable() noexcept
{
  return cpuidBit(7, 2, 12);
}

/*!
 * @brief AVX512VPOPCNTDQ命令が利用可能かどうかを調べる．
 * @return AVX512VPOPCNTDQ命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isAvx512VpopcntdqAvailable() noexcept
{
  return cpuidBit(7, 2, 14);
}

/*!
 * @brief AVX512VBMI命令が利用可能かどうかを調べる．
 * @return AVX512VBMI命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isAvx512VbmiAvailable() noexcept
{
  return cpuidBit(7, 2, 1);
}

/*!
 * @brief AVX512VBMI2命令が利用可能かどうかを調べる．
 * @return AVX512VBMI2命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isAvx512Vbmi2Available() noexcept
{
  return cpuidBit(7, 2, 6);
}

/*!
 * @brief AVX512VNNI命令が利用可能かどうかを調べる．
 * @return AVX512VNNI命令が利用可能ならばtrue，そうでなければfalse．
 */
static inline bool
isAvx512VnniAvailable() noexcept
{
  return cpuidBit(7, 2, 6);
}


//// 以下はおまけ


/*!
 * @brief CPUのベンダIDを第一引数のポインタの指すメモリ領域にコピーする
 *
 * 先頭から13byteの上書きを行う
 *
 * @tparam T  char*
 * @param [out] vendorId  CPUのベンダID
 */
template<
  typename T,
  typename std::enable_if<std::is_same<T, char*>::value, std::nullptr_t>::type = nullptr
>
static inline void
copyCpuVendorId(T vendorId) noexcept
{
  std::array<int, 4> cpuInfo;
  cpuid(cpuInfo, 0);

  const auto p = reinterpret_cast<int*>(vendorId);
  p[0] = cpuInfo[1];
  p[1] = cpuInfo[3];
  p[2] = cpuInfo[2];
  vendorId[12] = '\0';
}


/*!
 * @brief CPUのベンダIDを第一引数の配列にコピーする
 *
 * 配列の要素数は13個以上でなければならない
 *
 * @tparam kSize  配列のサイズ
 * @param [out] vendorId  CPUのベンダID
 */
template<std::size_t kSize>
static inline void
copyCpuVendorId(char (&vendorId)[kSize]) noexcept
{
  static_assert(kSize >= 12, "CPU vendor ID array size must be 12 or more");

  copyCpuVendorId(vendorId.data());
}


/*!
 * @brief CPUのベンダIDを第一引数のstd::arrayにコピーする
 *
 * std::arrayの要素数は13個以上でなければならない
 *
 * @tparam kSize  配列のサイズ
 * @param [out] vendorId  CPUのベンダID
 */
template<std::size_t kSize>
static inline void
copyCpuVendorId(std::array<char, kSize>& vendorId) noexcept
{
  static_assert(kSize >= 12, "CPU vendor ID array size must be 12 or more");

  copyCpuVendorId(vendorId.data());
}


/*!
 * @brief CPUのベンダIDをstd::stringとして得る
 * @return CPUのベンダID
 */
static inline std::string
getCpuVendorId() noexcept
{
  std::array<char, 32> vendorId;
  std::fill(std::begin(vendorId), std::end(vendorId), '\0');

  copyCpuVendorId(vendorId);

  return std::string{ vendorId.data() };
}


/*!
 * @brief CPUのブランド文字列を第一引数のポインタの指すメモリ領域にコピーする
 * @tparam T  char*
 * @param [out] brandString  ブランド文字列出力先配列
 */
template<
  typename T,
  typename std::enable_if<std::is_same<T, char*>::value, std::nullptr_t>::type = nullptr
>
static inline void
copyCpuBrandString(T brandString) noexcept
{
  std::array<int, 4> cpuInfo;

  cpuid(cpuInfo, 0x80000000);
  if (static_cast<unsigned int>(cpuInfo[0]) < 0x80000004) {
    brandString[0] = '\0';
    return;
  }

  const auto p = reinterpret_cast<int*>(brandString);

  cpuid(cpuInfo, 0x80000002);
  std::copy(std::begin(cpuInfo), std::end(cpuInfo), &p[0]);

  cpuid(cpuInfo, 0x80000003);
  std::copy(std::begin(cpuInfo), std::end(cpuInfo), &p[cpuInfo.size()]);

  cpuid(cpuInfo, 0x80000004);
  std::copy(std::begin(cpuInfo), std::end(cpuInfo), &p[cpuInfo.size() * 2]);
}


/*!
 * @brief CPUのブランド文字列を第一引数の配列にコピーする
 * @param [out] brandString  ブランド文字列出力先配列
 */
template<std::size_t kSize>
static inline void
copyCpuBrandString(char (&brandstring)[kSize]) noexcept
{
  static_assert(kSize >= 64, "CPU brand string array size must be 64 or more");

  copyCpuBrandString(brandstring);
}


/*!
 * @brief CPUのブランド文字列を第一引数のstd::arrayにコピーする
 * @param [out] brandString  ブランド文字列出力先配列
 */
template<std::size_t kSize>
static inline void
copyCpuBrandString(std::array<char, kSize>& brandstring) noexcept
{
  static_assert(kSize >= 64, "CPU brand string array size must be 64 or more");

  copyCpuBrandString(brandstring.data());
}


/*!
 * @brief CPUのブランド文字列をstd::stringとして得る
 * @return CPUのブランド文字列
 */
static inline std::string
getCpuBrandString() noexcept
{
  std::array<char, 64> brandStringArray;
  std::fill(std::begin(brandStringArray), std::end(brandStringArray), '\0');

  copyCpuVendorId(brandStringArray);

  return std::string{ brandStringArray.data() };
}


#endif  // CPUID_HPP

上記の関数を用いると，例えば，AVX2が利用可能であるかどうかは

auto hasAvx2 = isAvx2Available();

のようにして調べられる．

MSDNのcpuidのサンプルコード

ちなみに，MSDNにも __cpuid() を利用して利用可能なSIMD命令を調べるサンプルコードがある． このサンプルコードはMSVCではコンパイルできるが，gccではコンパイルできない． 両者共にコンパイルできるようにするなら，以下のように書き直すとよい．

Wandboxでの実行結果はこのようになる．

// InstructionSet.cpp Compile by using: cl /EHsc /W4 InstructionSet.cpp
// processor: x86, x64
// Uses the __cpuid intrinsic to get information about
// CPU extended instruction set support.

#include <algorithm>
#include <array>
#include <bitset>
#include <iostream>
#include <string>
#include <vector>


#if defined(__GNUC__)
#  include <cpuid.h>
#elif defined(_MSC_VER)
#  include <intrin.h>
#endif


static inline void
cpuid(int* cpuInfo, int eax) noexcept
{
#if defined(__GNUC__)
  __cpuid(eax, cpuInfo[0], cpuInfo[1], cpuInfo[2], cpuInfo[3]);
#elif defined(_MSC_VER)
  __cpuid(cpuInfo, eax);
#endif  // defined(__GNUC__)
}


template<std::size_t kSize>
static inline void
cpuid(std::array<int, kSize>& cpuInfo, int eax) noexcept
{
  static_assert(kSize >= 4, "CPU info array size must be four or more");
  cpuid(cpuInfo.data(), eax);
}


static inline void
cpuidex(int* cpuInfo, int eax, int ecx) noexcept
{
#if defined(__GNUC__)
  __cpuid_count(eax, ecx, cpuInfo[0], cpuInfo[1], cpuInfo[2], cpuInfo[3]);
#elif defined(_MSC_VER)
  __cpuidex(cpuInfo, eax, ecx);
#endif  // defined(__GNUC__)
}


template<std::size_t kSize>
static inline void
cpuidex(std::array<int, kSize>& cpuInfo, int eax, int ecx) noexcept
{
  static_assert(kSize >= 4, "CPU info array size must be four or more");
  cpuidex(cpuInfo.data(), eax, ecx);
}


class InstructionSet
{
  // forward declarations
  class InstructionSet_Internal;

public:
  // getters
  static std::string Vendor() noexcept { return CPU_Rep.vendor_; }
  static std::string Brand() noexcept { return CPU_Rep.brand_; }
  static bool SSE3() noexcept { return CPU_Rep.f_1_ECX_[0]; }
  static bool PCLMULQDQ() noexcept { return CPU_Rep.f_1_ECX_[1]; }
  static bool MONITOR() noexcept { return CPU_Rep.f_1_ECX_[3]; }
  static bool SSSE3() noexcept { return CPU_Rep.f_1_ECX_[9]; }
  static bool FMA() noexcept { return CPU_Rep.f_1_ECX_[12]; }
  static bool CMPXCHG16B() noexcept { return CPU_Rep.f_1_ECX_[13]; }
  static bool SSE41() noexcept { return CPU_Rep.f_1_ECX_[19]; }
  static bool SSE42() noexcept { return CPU_Rep.f_1_ECX_[20]; }
  static bool MOVBE() noexcept { return CPU_Rep.f_1_ECX_[22]; }
  static bool POPCNT() noexcept { return CPU_Rep.f_1_ECX_[23]; }
  static bool AES() noexcept { return CPU_Rep.f_1_ECX_[25]; }
  static bool XSAVE() noexcept { return CPU_Rep.f_1_ECX_[26]; }
  static bool OSXSAVE() noexcept { return CPU_Rep.f_1_ECX_[27]; }
  static bool AVX() noexcept { return CPU_Rep.f_1_ECX_[28]; }
  static bool F16C() noexcept { return CPU_Rep.f_1_ECX_[29]; }
  static bool RDRAND() noexcept { return CPU_Rep.f_1_ECX_[30]; }
  static bool MSR() noexcept { return CPU_Rep.f_1_EDX_[5]; }
  static bool CX8() noexcept { return CPU_Rep.f_1_EDX_[8]; }
  static bool SEP() noexcept { return CPU_Rep.f_1_EDX_[11]; }
  static bool CMOV() noexcept { return CPU_Rep.f_1_EDX_[15]; }
  static bool CLFSH() noexcept { return CPU_Rep.f_1_EDX_[19]; }
  static bool MMX() noexcept { return CPU_Rep.f_1_EDX_[23]; }
  static bool FXSR() noexcept { return CPU_Rep.f_1_EDX_[24]; }
  static bool SSE() noexcept { return CPU_Rep.f_1_EDX_[25]; }
  static bool SSE2() noexcept { return CPU_Rep.f_1_EDX_[26]; }
  static bool FSGSBASE() noexcept { return CPU_Rep.f_7_EBX_[0]; }
  static bool BMI1() noexcept { return CPU_Rep.f_7_EBX_[3]; }
  static bool HLE() noexcept { return CPU_Rep.isIntel_ && CPU_Rep.f_7_EBX_[4]; }
  static bool AVX2() noexcept { return CPU_Rep.f_7_EBX_[5]; }
  static bool BMI2() noexcept { return CPU_Rep.f_7_EBX_[8]; }
  static bool ERMS() noexcept { return CPU_Rep.f_7_EBX_[9]; }
  static bool INVPCID() noexcept { return CPU_Rep.f_7_EBX_[10]; }
  static bool RTM() noexcept { return CPU_Rep.isIntel_ && CPU_Rep.f_7_EBX_[11]; }
  static bool AVX512F() noexcept { return CPU_Rep.f_7_EBX_[16]; }
  static bool AVX512DQ() noexcept { return CPU_Rep.f_7_EBX_[17]; }
  static bool RDSEED() noexcept { return CPU_Rep.f_7_EBX_[18]; }
  static bool ADX() noexcept { return CPU_Rep.f_7_EBX_[19]; }
  static bool AVX512IFMA() noexcept { return CPU_Rep.f_7_EBX_[21]; }
  static bool AVX512PF() noexcept { return CPU_Rep.f_7_EBX_[26]; }
  static bool AVX512ER() noexcept { return CPU_Rep.f_7_EBX_[27]; }
  static bool AVX512CD() noexcept { return CPU_Rep.f_7_EBX_[28]; }
  static bool SHA() noexcept { return CPU_Rep.f_7_EBX_[29]; }
  static bool AVX512BW() noexcept { return CPU_Rep.f_7_EBX_[30]; }
  static bool AVX512VL() noexcept { return CPU_Rep.f_7_EBX_[31]; }
  static bool PREFETCHWT1() noexcept { return CPU_Rep.f_7_ECX_[0]; }
  static bool AVX512VBMI() noexcept { return CPU_Rep.f_7_ECX_[1]; }
  static bool AVX512VBMI2() noexcept { return CPU_Rep.f_7_ECX_[6]; }
  static bool AVX512VNNI() noexcept { return CPU_Rep.f_7_ECX_[11]; }
  static bool AVX512BITALG() noexcept { return CPU_Rep.f_7_ECX_[12]; }
  static bool AVX512VPOPCNTDQ() noexcept { return CPU_Rep.f_7_ECX_[14]; }
  static bool AVX512_4VNNIW() noexcept { return CPU_Rep.f_7_EDX_[2]; }
  static bool AVX512_4FMAPS() noexcept { return CPU_Rep.f_7_EDX_[3]; }
  static bool LAHF() noexcept { return CPU_Rep.f_81_ECX_[0]; }
  static bool LZCNT() noexcept { return CPU_Rep.isIntel_ && CPU_Rep.f_81_ECX_[5]; }
  static bool ABM() noexcept { return CPU_Rep.isAMD_ && CPU_Rep.f_81_ECX_[5]; }
  static bool SSE4a() noexcept { return CPU_Rep.isAMD_ && CPU_Rep.f_81_ECX_[6]; }
  static bool XOP() noexcept { return CPU_Rep.isAMD_ && CPU_Rep.f_81_ECX_[11]; }
  static bool TBM() noexcept { return CPU_Rep.isAMD_ && CPU_Rep.f_81_ECX_[21]; }
  static bool SYSCALL() noexcept { return CPU_Rep.isIntel_ && CPU_Rep.f_81_EDX_[11]; }
  static bool MMXEXT() noexcept { return CPU_Rep.isAMD_ && CPU_Rep.f_81_EDX_[22]; }
  static bool RDTSCP() noexcept { return CPU_Rep.isIntel_ && CPU_Rep.f_81_EDX_[27]; }
  static bool _3DNOWEXT() noexcept { return CPU_Rep.isAMD_ && CPU_Rep.f_81_EDX_[30]; }
  static bool _3DNOW() noexcept { return CPU_Rep.isAMD_ && CPU_Rep.f_81_EDX_[31]; }

private:
  static const InstructionSet_Internal CPU_Rep;

  class InstructionSet_Internal
  {
  public:
    InstructionSet_Internal()
      : nIds_{0}
      , nExIds_{0}
      , vendor_{}
      , brand_{}
      , isIntel_{false}
      , isAMD_{false}
      , f_1_ECX_{0}
      , f_1_EDX_{0}
      , f_7_EBX_{0}
      , f_7_ECX_{0}
      , f_7_EDX_{0}
      , f_81_ECX_{0}
      , f_81_EDX_{0}
      , data_{}
      , extdata_{}
    {
      std::array<int, 4> cpui;

      // Calling __cpuid with 0x0 as the function_id argument
      // gets the number of the highest valid function ID.
      cpuid(cpui, 0);
      nIds_ = cpui[0];

      for (int i = 0; i <= nIds_; ++i) {
        cpuidex(cpui, i, 0);
        data_.push_back(cpui);
      }

      // Capture vendor string
      std::array<char, 0x20> vendor;
      std::fill(std::begin(vendor), std::end(vendor), '\0');
      *reinterpret_cast<int*>(&vendor[0]) = data_[0][1];
      *reinterpret_cast<int*>(&vendor[4]) = data_[0][3];
      *reinterpret_cast<int*>(&vendor[8]) = data_[0][2];
      vendor_ = std::string(vendor.data());
      if (vendor_ == "GenuineIntel") {
        isIntel_ = true;
      } else if (vendor_ == "AuthenticAMD") {
        isAMD_ = true;
      }

      // load bitset with flags for function 0x00000001
      if (nIds_ >= 1) {
        f_1_ECX_ = data_[1][2];
        f_1_EDX_ = data_[1][3];
      }

      // load bitset with flags for function 0x00000007
      if (nIds_ >= 7) {
        f_7_EBX_ = data_[7][1];
        f_7_ECX_ = data_[7][2];
        f_7_EDX_ = data_[7][3];
      }

      // Calling __cpuid with 0x80000000 as the function_id argument
      // gets the number of the highest valid extended ID.
      cpuid(cpui, 0x80000000);
      nExIds_ = cpui[0];

      std::array<char, 0x40> brand;
      std::fill(std::begin(brand), std::end(brand), '\0');

      for (int i = 0x80000000; i <= nExIds_; ++i) {
        cpuidex(cpui, i, 0);
        extdata_.push_back(cpui);
      }

      // load bitset with flags for function 0x80000001
      if (static_cast<unsigned int>(nExIds_) >= 0x80000001) {
        f_81_ECX_ = extdata_[1][2];
        f_81_EDX_ = extdata_[1][3];
      }

      // Interpret CPU brand string if reported
      if (static_cast<unsigned int>(nExIds_) >= 0x80000004) {
        std::copy(std::cbegin(extdata_[2]), std::cend(extdata_[2]), reinterpret_cast<int*>(&brand[0]));
        std::copy(std::cbegin(extdata_[3]), std::cend(extdata_[3]), reinterpret_cast<int*>(&brand[0] + sizeof(extdata_[0])));
        std::copy(std::cbegin(extdata_[4]), std::cend(extdata_[4]), reinterpret_cast<int*>(&brand[0] + sizeof(extdata_[0]) * 2));
        brand_ = std::string(brand.data());
      }
    };

    int nIds_;
    int nExIds_;
    std::string vendor_;
    std::string brand_;
    bool isIntel_;
    bool isAMD_;
    std::bitset<32> f_1_ECX_;
    std::bitset<32> f_1_EDX_;
    std::bitset<32> f_7_EBX_;
    std::bitset<32> f_7_ECX_;
    std::bitset<32> f_7_EDX_;
    std::bitset<32> f_81_ECX_;
    std::bitset<32> f_81_EDX_;
    std::vector<std::array<int, 4>> data_;
    std::vector<std::array<int, 4>> extdata_;
  };  // class InstructionSet_Internal
};  // class InstructionSet


// Initialize static member data
const InstructionSet::InstructionSet_Internal InstructionSet::CPU_Rep;


// Print out supported instruction set extensions
int
main()
{
  auto &outstream = std::cout;
  auto support_message = [&outstream](std::string isa_feature, bool is_supported) {
    outstream << isa_feature
              << (is_supported ? " supported" : " not supported")
              << std::endl;
  };

  std::cout << InstructionSet::Vendor() << std::endl;
  std::cout << InstructionSet::Brand() << std::endl;

  support_message("3DNOW", InstructionSet::_3DNOW());
  support_message("3DNOWEXT", InstructionSet::_3DNOWEXT());
  support_message("ABM", InstructionSet::ABM());
  support_message("ADX", InstructionSet::ADX());
  support_message("AES", InstructionSet::AES());
  support_message("AVX", InstructionSet::AVX());
  support_message("AVX2", InstructionSet::AVX2());
  support_message("AVX512CD", InstructionSet::AVX512CD());
  support_message("AVX512ER", InstructionSet::AVX512ER());
  support_message("AVX512F", InstructionSet::AVX512F());
  support_message("AVX512DQ", InstructionSet::AVX512DQ());
  support_message("AVX512IFMA", InstructionSet::AVX512IFMA());
  support_message("AVX512PF", InstructionSet::AVX512PF());
  support_message("AVX512BW", InstructionSet::AVX512BW());
  support_message("AVX512VL", InstructionSet::AVX512VL());
  support_message("AVX512VBMI", InstructionSet::AVX512VBMI());
  support_message("AVX512VBMI2", InstructionSet::AVX512VBMI2());
  support_message("AVX512VNNI", InstructionSet::AVX512VNNI());
  support_message("AVX512BITALG", InstructionSet::AVX512BITALG());
  support_message("AVX512VPOPCNTDQ", InstructionSet::AVX512VPOPCNTDQ());
  support_message("AVX512_4VNNIW", InstructionSet::AVX512_4VNNIW());
  support_message("AVX512_4FMAPS", InstructionSet::AVX512_4FMAPS());
  support_message("BMI1", InstructionSet::BMI1());
  support_message("BMI2", InstructionSet::BMI2());
  support_message("CLFSH", InstructionSet::CLFSH());
  support_message("CMPXCHG16B", InstructionSet::CMPXCHG16B());
  support_message("CX8", InstructionSet::CX8());
  support_message("ERMS", InstructionSet::ERMS());
  support_message("F16C", InstructionSet::F16C());
  support_message("FMA", InstructionSet::FMA());
  support_message("FSGSBASE", InstructionSet::FSGSBASE());
  support_message("FXSR", InstructionSet::FXSR());
  support_message("HLE", InstructionSet::HLE());
  support_message("INVPCID", InstructionSet::INVPCID());
  support_message("LAHF", InstructionSet::LAHF());
  support_message("LZCNT", InstructionSet::LZCNT());
  support_message("MMX", InstructionSet::MMX());
  support_message("MMXEXT", InstructionSet::MMXEXT());
  support_message("MONITOR", InstructionSet::MONITOR());
  support_message("MOVBE", InstructionSet::MOVBE());
  support_message("MSR", InstructionSet::MSR());
  support_message("OSXSAVE", InstructionSet::OSXSAVE());
  support_message("PCLMULQDQ", InstructionSet::PCLMULQDQ());
  support_message("POPCNT", InstructionSet::POPCNT());
  support_message("PREFETCHWT1", InstructionSet::PREFETCHWT1());
  support_message("RDRAND", InstructionSet::RDRAND());
  support_message("RDSEED", InstructionSet::RDSEED());
  support_message("RDTSCP", InstructionSet::RDTSCP());
  support_message("RTM", InstructionSet::RTM());
  support_message("SEP", InstructionSet::SEP());
  support_message("SHA", InstructionSet::SHA());
  support_message("SSE", InstructionSet::SSE());
  support_message("SSE2", InstructionSet::SSE2());
  support_message("SSE3", InstructionSet::SSE3());
  support_message("SSE4.1", InstructionSet::SSE41());
  support_message("SSE4.2", InstructionSet::SSE42());
  support_message("SSE4a", InstructionSet::SSE4a());
  support_message("SSSE3", InstructionSet::SSSE3());
  support_message("SYSCALL", InstructionSet::SYSCALL());
  support_message("TBM", InstructionSet::TBM());
  support_message("XOP", InstructionSet::XOP());
  support_message("XSAVE", InstructionSet::XSAVE());
}

cpuid命令自体が利用可能かどうかを調べる

cpuid命令自体が利用可能かどうかも調べる必要があるのではないか？と疑問を持たれる人もいるかもしれない． 実はその通りで，かなり昔のCPUではcpuid命令がなかったらしい．

cpuid命令が利用可能かどうかは，インテルのドキュメントに記載してあるように，eflagsの21bit目が変更可能であるかどうかを調べるとよい．

ただし，これはC言語，C++で記述することはできないので，インラインアセンブラに頼る必要がある．

#if defined(_MSC_VER) && defined(_WIN64)
#  ifndef WIN32_LEAN_AND_MEAN
#    define WIN32_LEAN_AND_MEAN
#    define CPUID_WIN32_LEAN_AND_MEAN_IS_NOT_DEFINED
#  endif  // !WIN32_LEAN_AND_MEAN
#  ifndef NOMINMAX
#    define NOMINMAX
#    define CPUID_NOMINMAX_IS_NOT_DEFINED
#  endif  // !NOMINMAX
#  include <windows.h>
#  ifdef CPUID_WIN32_LEAN_AND_MEAN_IS_NOT_DEFINED
#    undef CPUID_WIN32_LEAN_AND_MEAN_IS_NOT_DEFINED
#    undef WIN32_LEAN_AND_MEAN
#  endif  // CPUID_WIN32_LEAN_AND_MEAN_IS_NOT_DEFINED
#  ifdef CPUID_NOMINMAX_IS_NOT_DEFINED
#    undef CPUID_NOMINMAX_IS_NOT_DEFINED
#    undef NOMINMAX
#  endif  // CPUID_NOMINMAX_IS_NOT_DEFINED
#endif  // defined(_MSC_VER) && defined(_WIN64)


static inline bool
isCpuidSupported() noexcept
{
#if defined(__x86_64__) || defined(_WIN64) || defined(__MINGW64__)
  // x64とき (全てのIntel x64プロセッサではcpuid命令は利用可能なため，このように真面目に調べる必要はない)
#  if defined(__GNUC__)
  bool result;
  __asm__ __volatile__ (
    "pushfq\n\t"
    "pushfq\n\t"
    "pop %%rax\n\t"
    "mov %%rax, %%rcx\n\t"
    "xor $0x200000, %%rax\n\t"
    "push %%rax\n\t"
    "popfq\n\t"
    "pushfq\n\t"
    "pop %%rax\n\t"
    "xor %%rcx, %%rax\n\t"
    "shr $21, %%rax\n\t"
    "popfq\n\t"
    : "=a" (result)
    :
    : "cc", "%rcx");
  return result;
#  elif defined(_MSC_VER)
  // MSVCのx64ではインラインアセンブラを利用できないので、
  // マシンコード配列を用意し、そのメモリ領域に実行権限を与えて、
  // eflagsの21bit目が変更可能かどうかを調べる

  // cdecl function code
  std::uint8_t code[] = {
    0x9c,                                // pushfq
    0x9c,                                // pushfq
    0x58,                                // pop     rax
    0x48, 0x89, 0xc1,                    // mov     rcx,rax
    0x48, 0x35, 0x00, 0x00, 0x20, 0x00,  // xor     rax,200000h
    0x50,                                // push    rax
    0x9d,                                // popfq
    0x9c,                                // pushfq
    0x58,                                // pop     rax
    0x48, 0x31, 0xc8,                    // xor     rax,rcx
    0x48, 0xc1, 0xe8, 0x15,              // shr     rax,21
    0x9d,                                // popfq
    0xc3                                 // ret
  };
  ::DWORD oldProtect;
  ::VirtualProtect(code, sizeof(code), PAGE_EXECUTE_READWRITE, &oldProtect);
  const auto result = reinterpret_cast<bool(__cdecl*)()>(reinterpret_cast<unsigned char*>(code))();
  ::VirtualProtect(code, sizeof(code), oldProtect, &oldProtect);
  return result;
#  endif  // defined(__GNUC__)
#else
  // x86のとき
#  if defined(__GNUC__)
  bool result;
  __asm__ __volatile__ (
    "pushfl\n\t"
    "pushfl\n\t"
    "pop %%eax\n\t"
    "mov %%eax, %%ecx\n\t"
    "xorl $0x200000, %%eax\n\t"
    "push %%eax\n\t"
    "popfl\n\t"
    "pushfl\n\t"
    "pop %%eax\n\t"
    "xorl %%ecx, %%eax\n\t"
    "shrl $21, %%eax\n\t"
    "popfl\n\t"
    : "=a" (result)
    :
    : "cc", "%ecx");
  return result;
#  elif defined(_MSC_VER)
  bool result;
  __asm {
    pushfd
    pushfd
    pop eax
    mov ecx, eax
    xor eax, 200000h
    push eax
    popfd
    pushfd
    pop eax
    xor eax, ecx
    shr eax, 21
    mov result, al
    popfd
  }
  return result;
#  endif  // defined(__GNUC__)
#endif  // defined(__x86_64__) || defined(_WIN64) || defined(__MINGW64__)
}

Intelによると，全てのx64プロセッサでcpuid命令が利用可能であるため，x64の方のコードは不要で，常に true を返すようにしてもよい．

まとめ

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

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

特に，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をプログラムに組み込めるようになるかもしれない．

参考文献

• Intel Intrinsics Guide
• Intel® Software Development Emulator | Intel® Software
• さかな前線 » SSEとAVXで高次元ベクトルの内積計算を高速化してみた
• SSE.浮動小数点演算手動最適化は本当に効果的なのか - デー
• SIMD演算 - MUGI COM
• 概要: ストリーミング SIMD 拡張命令
• x86/x64 SIMD命令一覧表　（SSE～AVX2）
• ARM NEON Intrinsics - Using the GNU Compiler Collection (GCC)
• ARM NEON Development
• SIMD Assembly Tutorial: ARM NEON
• NEON を使用して Zynq-7000 AP SoC でのソフトウェア性能を向上
• ARM gcc バッドノウハウ集
• 2012 Intel® Processor Identification and the CPUID Instruction
• AMD CPUID Specification
• __cpuid, __cpuidex | Microsoft Docs