はじめに
現代の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 追記
以下の記事に，より詳細な内容を書いたので，参考になるかもしれない．
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
#endif
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
#endif
#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() はアラインされた静的配列を作成するための関数なので，やや使い勝手が悪いといえる．
#include <cstddef>
#include <iostream>
#include <memory>
#include <type_traits>
#if defined(_MSC_VER) || defined(__MINGW32__)
# include <malloc.h>
#else
# include <cstdlib>
#endif
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
}
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);
}
static inline void
alignedFree(void* ptr) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
::_aligned_free(ptr);
#else
std::free(ptr);
#endif
}
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
#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
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
}
static inline void
alignedFree(void* ptr)
{
#if defined(_MSC_VER) || defined(__MINGW32__)
_aligned_free(ptr);
#else
free(ptr);
#endif
}
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
#endif
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
#elif defined(ENABLE_NEON)
# include <arm_neon.h>
#endif
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
}
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);
}
static inline void
alignedFree(void* ptr) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
::_aligned_free(ptr);
#else
std::free(ptr);
#endif
}
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
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
}
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
}
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
}
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++) {
sum = std::fma(a[i], b[i], sum);
}
return sum;
#endif
}
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
最近傍法による画像の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
#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
#elif defined(ENABLE_NEON)
# include <arm_neon.h>
#endif
#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
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
}
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);
}
static inline void
alignedFree(void* ptr) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
::_aligned_free(ptr);
#else
std::free(ptr);
#endif
}
#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
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
for (int i = 0; i < dstHeight; i++) {
for (int j = 0; j < dstWidth; j += INTERVAL * X_RATIO) {
#if defined(ENABLE_AVX512)
__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);
_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)
__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);
__m256i v256l = _mm256_permute2f128_si256(v256l_, v256u_, 0x20);
__m256i v256u = _mm256_permute2f128_si256(v256l_, v256u_, 0x31);
_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)
__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);
_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)
uint8x16_t v128 = vld1q_u8(&srcImageData[i / Y_RATIO * srcWidth + j]);
uint8x16x2_t v128x2 = vzipq_u8(v128, v128);
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
}
}
}
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
何の脈略も無しに，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
mov $0,%ecx
cpuid
まず 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
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
}
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);
}
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);
}
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
}
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);
}
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
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
}
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
}
これで，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に置いてある．
#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
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
}
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);
}
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);
}
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
}
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);
}
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);
}
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;
}
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;
}
static inline bool
isMmxAvailable() noexcept
{
return cpuidBit(1, 3, 23);
}
static inline bool
isSseAvailable() noexcept
{
return cpuidBit(1, 3, 25);
}
static inline bool
isSse2Available() noexcept
{
return cpuidBit(1, 3, 26);
}
static inline bool
isSse3Available() noexcept
{
return cpuidBit(1, 2, 0);
}
static inline bool
isSsse3Available() noexcept
{
return cpuidBit(1, 2, 9);
}
static inline bool
isSse41Available() noexcept
{
return cpuidBit(1, 2, 19);
}
static inline bool
isSse42Available() noexcept
{
return cpuidBit(1, 2, 20);
}
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);
}
static inline bool
isAvxAvailable() noexcept
{
return cpuidBit(1, 2, 28);
}
static inline bool
isAvx2Available() noexcept
{
return cpuidBit(7, 1, 5);
}
static inline bool
isFmaAvailable() noexcept
{
return cpuidBit(1, 2, 12);
}
static inline bool
isAvx512FAvailable() noexcept
{
return cpuidBit(7, 1, 16);
}
static inline bool
isAvx512BwAvailable() noexcept
{
return cpuidBit(7, 1, 30);
}
static inline bool
isAvx512CdAvailable() noexcept
{
return cpuidBit(7, 1, 28);
}
static inline bool
isAvx512DqAvailable() noexcept
{
return cpuidBit(7, 1, 17);
}
static inline bool
isAvx512ErAvailable() noexcept
{
return cpuidBit(7, 1, 27);
}
static inline bool
isAvx512Ifma52Available() noexcept
{
return cpuidBit(7, 1, 21);
}
static inline bool
isAvx512PfAvailable() noexcept
{
return cpuidBit(7, 1, 26);
}
static inline bool
isAvx512VlAvailable() noexcept
{
return cpuidBit(7, 1, 31);
}
static inline bool
isAvx512_4fmapsAvailable() noexcept
{
return cpuidBit(7, 3, 2);
}
static inline bool
isAvx512_4vnniwAvailable() noexcept
{
return cpuidBit(7, 3, 3);
}
static inline bool
isAvx512BitalgAvailable() noexcept
{
return cpuidBit(7, 2, 12);
}
static inline bool
isAvx512VpopcntdqAvailable() noexcept
{
return cpuidBit(7, 2, 14);
}
static inline bool
isAvx512VbmiAvailable() noexcept
{
return cpuidBit(7, 2, 1);
}
static inline bool
isAvx512Vbmi2Available() noexcept
{
return cpuidBit(7, 2, 6);
}
static inline bool
isAvx512VnniAvailable() noexcept
{
return cpuidBit(7, 2, 6);
}
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';
}
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());
}
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());
}
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() };
}
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]);
}
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);
}
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());
}
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
上記の関数を用いると，例えば，AVX2が利用可能であるかどうかは
auto hasAvx2 = isAvx2Available();
のようにして調べられる．
MSDNのcpuidのサンプルコード
ちなみに，MSDNにも __cpuid() を利用して利用可能なSIMD命令を調べるサンプルコードがある．
このサンプルコードはMSVCではコンパイルできるが，gccではコンパイルできない．
両者共にコンパイルできるようにするなら，以下のように書き直すとよい．
Wandboxでの実行結果はこのようになる．
#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
}
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
}
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
{
class InstructionSet_Internal;
public:
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;
cpuid(cpui, 0);
nIds_ = cpui[0];
for (int i = 0; i <= nIds_; ++i) {
cpuidex(cpui, i, 0);
data_.push_back(cpui);
}
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;
}
if (nIds_ >= 1) {
f_1_ECX_ = data_[1][2];
f_1_EDX_ = data_[1][3];
}
if (nIds_ >= 7) {
f_7_EBX_ = data_[7][1];
f_7_ECX_ = data_[7][2];
f_7_EDX_ = data_[7][3];
}
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);
}
if (static_cast<unsigned int>(nExIds_) >= 0x80000001) {
f_81_ECX_ = extdata_[1][2];
f_81_EDX_ = extdata_[1][3];
}
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_;
};
};
const InstructionSet::InstructionSet_Internal InstructionSet::CPU_Rep;
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
# ifndef NOMINMAX
# define NOMINMAX
# define CPUID_NOMINMAX_IS_NOT_DEFINED
# endif
# 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
# ifdef CPUID_NOMINMAX_IS_NOT_DEFINED
# undef CPUID_NOMINMAX_IS_NOT_DEFINED
# undef NOMINMAX
# endif
#endif
static inline bool
isCpuidSupported() noexcept
{
#if defined(__x86_64__) || defined(_WIN64) || defined(__MINGW64__)
# 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)
std::uint8_t code[] = {
0x9c,
0x9c,
0x58,
0x48, 0x89, 0xc1,
0x48, 0x35, 0x00, 0x00, 0x20, 0x00,
0x50,
0x9d,
0x9c,
0x58,
0x48, 0x31, 0xc8,
0x48, 0xc1, 0xe8, 0x15,
0x9d,
0xc3
};
::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
#else
# 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
#endif
}
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をプログラムに組み込めるようになるかもしれない．
参考文献
