【C++】x86のSIMDを触る【西田の戯言。】

西田の戯言。

なんか最近SIMDがアツいです。といいつつ、SIMDはきっとフレームワーク経由で触ることになると思います。今回は逆に直で触ってみようって話です。

SIMDって？

SIMDは、簡単に言うと、1命令で複数のデータに同じ処理を適用するというものです。フリンの分類の1種ですね。

軽くイメージを。

CPUのSIMDは、レジスタを区切ってレジスタよりも低精度の複数の値に対して、1命令で同じ命令を適用するというものです。

初期のSSE命令がこんな感じ。

これは、128-bitのレジスタを32-bitごとに区切って、一つの命令で32-bitの値に対して同じ演算を実行したり、ベクトル演算を適用する事ができるというものです。

あくまでこれは例であり、最大512-bitのレジスタで64個の8-bitの数を演算することもできます。

今回は、主にCPUのベクトル拡張について扱いますが、同じくSIMDに分類されるマトリックス演算器あるいは拡張では、1024-bitなどの大きな演算器を有しており、それで64個の16-bitの値を1命令で処理できることになります。

x86_64におけるSIMD拡張命令は大きく「SSE」「AVX」「AVX-512」と分けることができます。SSEは128-bit、AVXは256-bit、AVX-512は512-bitのレジスタを必要としています。

SSEが必要とする128-bitレジスタは「XMM」と呼ばれます。

AVXが必要とする256-bitレジスタは「YMM」と呼びます。XMMとYMMは共存ではなくAVX対応CPUがSSE命令を実行するときに、YMMの下位128-bitをXMMとして利用します。

AVX-512も同様に512-bitレジスタ「ZMM」を搭載していますが、SSE命令実行時には下位128-bitをXMMとして、AVX命令実行時には下位256-bitをYMMとして利用します。

実際には、SSEは「SSE」「SSE2」「SSE3」...など、AVXは「AVX」「AVX2」などというマイナーアップデート（？）が存在しており、世代やマイクロアーキテクチャによって対応する拡張命令セットに差があります。

C++で触るぜ

さて、ここからはこれらのSIMD命令を触っていきましょう。

動作環境

  • Ryzen 7 5700X(MMX/SSE1-4.2/SSSE3/AVX/AVX2)
  • Windows 11 Pro
  • gcc-g++

本当はAVX-512環境を用意したかったけど担当日までに用意できそうにないので、用意できたら検証します。

ヘッダーファイル

基本的に、WindowsとLinuxではCPUの組み込み命令を搭載したヘッダーファイルが既存で用意されているはずです。何もインストールせず、そのまま使えると思います（Intel Macが手元にないので、Macは検証できませんでした）。

#include<x86intrin.h>

x86intrin.hには、すべての機能が含まれていますが、各拡張命令ごとにヘッダーファイルが別れて存在してはいます。

  • mmintrin.h:MMX
  • xmmintrin.h:SSE
  • emmintrin.h:SSE2
  • pmmintrin.h:SSE3
  • tmmintrin.h:SSSE3
  • 2mmintrin.h:SSE4.1
  • nmmintrin.h:SSE4.2
  • wmmintrin.h:AES
  • immintrin.h:AVX/AVX2/FMA
  • zmmintrin.h:AVX-512

今回は基本的にAVX2を中心に話は進めていきます。

これらのヘッダーファイルに含まれている関数は「組み込み関数」（intrinsic）とも呼ばれており、CPUの機能が含まれています。

余談ですが、Armの場合、arm-neon.hなどを用います。残念ながら今回紹介するコードたちは使えません。

専用レジスタにロードする

まず第一歩、AVXレジスタにfloat型のベクトルを格納してみましょう

float型関数は32-bitですので、256-bitのAVXレジスタに8個入ります。

alignas(32) float a[8] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};
alignas(32) float b[8] = {8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0};

// AVXレジスタにロード
__m256 vec_a = _mm256_load_ps(a);
__m256 vec_b = _mm256_load_ps(b);

__mm256型は、AVX用に用意されている新しい型で、32-bit浮動小数点数を8個格納します。

alignas(n)は、アラインメントを指定しているものです。アラインメントについては後述します。

alignas(n)の引数には、nバイト境界にアラインメントされるという値が入ります。具体的な値は以下のように設定してあげるといいです。

値
SSE4.2以前 16
AVX/AVX2 32
AVX-512 64

同様に以下の型が定義されています。

対応命令 型名 格納できる数値
SSE以降 __m128 4x FP32
SSE2以降 __m128d 2x FP64
SSE2以降 __m128i 16x 8-bit 整数
8x 16-bit 整数
4x 32-bit 整数
2x 64-bit 整数
AVX以降 __m256 8x FP32
AVX以降 __m256d 4x FP64
AVX以降 __m256i 32x 8-bit 整数
16x 16-bit 整数
8x 32-bit 整数
4x 64-bit 整数
AVX-512 __m512 16x FP32
AVX-512 __m512d 8x FP64
AVX-512 __m512i 64x 8-bit 整数
32x 16-bit 整数
16x 32-bit 整数
8x 64-bit 整数

ロード用の関数は、基本的に_mm{レジスタサイズ}_load_{格納する値の型}の形となります。

以下のような形です。

// YMMにdouble型を格納する(AVX)
__m256d = _mm_load_pd(a);

// ZMMに整数を格納する(AVX-512)
__m512i = _mm_load_si512(a);

格納する型の部分は様々ですが、基本float型ならps、double型ならpd、整数型ならsi{レジスタ幅}になります。

また、別の方法として、直で値を格納する方法もあります。

__m256 a = _mm256_set_ps(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0);

四則演算する

// ベクトル加算
__m256 vec_add_result = _mm256_add_ps(vec_a, vec_b);

// ベクトル減算
__m256 vec_sub_result = _mm256_sub_ps(vec_a, vec_b);

// ベクトル乗算
__m256 vec_mul_result = _mm256_mul_ps(vec_a, vec_b);

// ベクトル除算
__m256 vec_div_result = _mm256_div_ps(vec_a, vec_b);

基本的に256の数字を変えれば、ベクトル幅を変える事ができます。psの部分は、格納する型によって変えてください。

以下、足し算のみにします。

メモリにストア

__m256型などではそのまま出力することができません。なので、一度float型の配列に値を戻します。

// 結果を格納するベクトル
float add_result[8];

// メモリにストアする
__mm256_store_ps(dd_result, vec_add_result);

出力は普通にfloat型配列の出力なのでfor文を使うなどして出力してください。

結果は以下のとおりです。

add:9 9 9 9 9 9 9 9 
sub:-7 -5 -3 -1 1 3 5 7 
mul:8 14 18 20 20 18 14 8 
div:0.125 0.285714 0.5 0.8 1.25 2 3.5 8

コードの全体像

#include <x86intrin.h>
#include <iostream>

void output(float output[8]){
    for (int i = 0; i < 8; ++i) {
        std::cout << output[i] << " ";
    }
    std::cout << std::endl;
}

int main() {

    // YMMに格納する値の配列
    alignas(32)float a[8] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};
    alignas(32)float b[8] = {8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0};

    // 出力用の配列
    float add_result[8];
    float sub_result[8];
    float mul_result[8];
    float div_result[8];


    // AVXレジスタにロード
    __m256 vec_a = _mm256_load_ps(a);
    __m256 vec_b = _mm256_load_ps(b);

    // ベクトル加算
    __m256 vec_add_result = _mm256_add_ps(vec_a, vec_b);

    // ベクトル減算
    __m256 vec_sub_result = _mm256_sub_ps(vec_a, vec_b);

    // ベクトル乗算
    __m256 vec_mul_result = _mm256_mul_ps(vec_a, vec_b);

    // ベクトル除算
    __m256 vec_div_result = _mm256_div_ps(vec_a, vec_b);

    // メモリにストア
    _mm256_store_ps(add_result, vec_add_result);
    _mm256_store_ps(sub_result, vec_sub_result);
    _mm256_store_ps(mul_result, vec_mul_result);
    _mm256_store_ps(div_result, vec_div_result);
    
    // 出力
    std::cout << "add:";
    output(add_result);

    std::cout << "sub:";
    output(sub_result);

    std::cout << "mul:";
    output(mul_result);

    std::cout << "div:";
    output(div_result);

    return 0;
}

追記

本来はアラインメントを指定することが望まれるみたいです。指摘されて私自身初めて気が付きました。非常に恥ずかしい（ご指摘ありがとうございます）。

具体的な対応としては、alignas()を使用して指定する方法。ただし、これはfloat配列を格納する際などにしか使えない気がします。

AVXの場合、32を指定してあげると良いみたいです。

// YMMに格納する値の配列
alignas(32)float a[8] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};
alignas(32)float b[8] = {8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0};

このような形です。

実際に、アラインメントを指定して上げることで、より高速のデータ処理が期待されます。

他方で、アラインメントを気にせずにロード・ストアする方法もあります。_mm256_loadu_ps(a)や_mm256_storeu_ps(a,b)を使う方法です。

// AVXレジスタにロード
__m256 vec_a = _mm256_loadu_ps(a);
__m256 vec_b = _mm256_loadu_ps(b);

/* 中略 */

// メモリにストア
_mm256_storeu_ps(add_result, vec_add_result);
_mm256_storeu_ps(sub_result, vec_sub_result);
_mm256_storeu_ps(mul_result, vec_mul_result);
_mm256_storeu_ps(div_result, vec_div_result);

ただし、この方法では、アラインメントをしないため動的メモリを扱っている事になってしまい、メモリアクセスに時間がかかることになり、速度が出ないという問題もあります。

コンパイル

g++等でコンパイルするときには、使用してほしい拡張命令セットを指定すると良いです。

使用できるオプションは以下の通り

-mmmx,-msse,-msse2,-msse3,-mssse3,-msse4,-msse4a,-msse4.1,-msse4.2,-mavx,-mavx2,-mavx512f,-mavx512cd,-mavx512vl,-mavx512bw,-mavx512dq,-mavx512ifma,-mavx512vbmi

ちなみに、これらのオプションをつける事によって、組み込み関数を使用しなくとも、拡張命令を使用してくれる場合があります。

また、上記のオプションの代わりに、CPUのマイクロアーキテクチャをターゲット（-march=）として、-ftree-vectorizeというオプション、あるいは-O3最適化を指定することで、ターゲットにしたCPUに対して最適なベクトル化を行ってくれます。詳細は「自動でベクトル・SIMD化させる」の節にて。

そのあたりの詳細は、GNUのこれを見ると良いです。

https://gcc.gnu.org/onlinedocs/gcc/x86-Options.html

自動でベクトル・SIMD化させる

次に、組み込み関数を使わずにSIMD化させる方法について紹介します。

このようなソースコードを用意しました。

void store_ary(float* __restrict a, float* __restrict b){
    for(int i = 0; i < 100; ++i){
        a[i] = i;
        b[i] = 100 - i;
    }
}

void add(const float* __restrict a, const float* __restrict b, float* __restrict c){
    for (int i = 0; i < 100; ++i) {
        c[i] = a[i] + b[i];
    }
}

int main() {
    constexpr int n = 100;
    float a[n], b[n], c[n];

    store_ary(a, b);
    add(a, b, c);

    return 0;
}

このソースコードは、2つのfloat配列の要素を足していくという至ってシンプルなコード。これを、いくつかの方法でコンパイルしてみます。

説明の便宜上、A~Eの記号をつけてます。

# A:何も指定せずにコンパイル
g++ -o tmp.s ./tmp.cpp -S -march=znver3 -fopt-info-vec-optimized

# B:O1最適化を有効にしてコンパイル
g++ -O1 -o tmp.s ./tmp.cpp -S -march=znver3 -fopt-info-vec-optimized

# C:O2最適化を有効にしてコンパイル
g++ -O2 -o tmp.s ./tmp.cpp -S -march=znver3 -fopt-info-vec-optimized

# D:O2最適化を有効にしたうえで、-ftree-vectorize最適化も有効にする
g++ -O2 -o tmp.s ./tmp.cpp -S -march=znver3 -ftree-vectorize -fopt-info-vec-optimized

# E:O3最適化を有効にする
g++ -O3 -o tmp.s ./tmp.cpp -S -march=znver3 -fopt-info-vec-optimized

-fopt-info-vec-optimizedは、コンパイラがどのようにベクトル化を行ったかと言うのを出力してくれるオプションです。-Sはアセンブラ出力を行うオプションです。この2つのオプションについては、検証のために使用します。

-ftree-vectorizeは、ターゲットで使用可能なベクトル化を行うというオプションになります。-Oオプションは最適化レベルを示すオプションとなっています。今回のDとEの場合、ターゲットをRyzen 7 5700Xに合わせるために意図的に-march=znver3としていますので、AVX2を使用してベクトル化を行ってくれます。

最適化レベルとベクトル化は関係しています。Gentoo LinuxのWikiを参考にすると、-O2最適化で128-bit以内でベクトル化を行い、-O3最適化では-ftree-vectorizeがデフォルトで有効になっており、使用可能な範囲でベクトル化を行います。

-fopt-info-vec-optimizedの出力をそれぞれ見てみます。

#A
出力なし

#B
出力なし

#C
.\simd2.cpp:4:22: optimized: loop vectorized using 16 byte vectors
.\simd2.cpp:11:23: optimized: loop vectorized using 16 byte vectors

#D
.\simd2.cpp:4:22: optimized: loop vectorized using 32 byte vectors
.\simd2.cpp:4:22: optimized: loop vectorized using 16 byte vectors
.\simd2.cpp:11:23: optimized: loop vectorized using 32 byte vectors
.\simd2.cpp:11:23: optimized: loop vectorized using 16 byte vectors

#E
.\simd2.cpp:4:22: optimized: loop vectorized using 32 byte vectors
.\simd2.cpp:4:22: optimized: loop vectorized using 16 byte vectors
.\simd2.cpp:11:23: optimized: loop vectorized using 32 byte vectors
.\simd2.cpp:11:23: optimized: loop vectorized using 16 byte vectors

この結果を見ますと、ターゲット以外を何も指定しない場合とO1最適化に置いてはベクトル化は行われず、O2最適化では、128-bit（16-byte）の範囲でベクトル化、-ftree-vectorizeを有効にするかO3最適化を指定することで、AVXの256-bit（32-byte）での最適化が行われているということがわかります。

-ftree-vectorizeについては、ターゲットで使用できるベクトル化を使用するという事になっていますが、この拡張命令セットは指定することができます。これは、前節で触れたとおりです。具体的には-mavx2とつければ、AVX2の範囲でベクトル化してくれるし、-mavx512fとつければAVX-512の範囲でベクトル化してくれます。一例として実行環境はありませんが、-ftree-vectorizeとともに-mavx512fと-O3最適化をつけてコンパイルすると、

.\simd2.cpp:4:22: optimized: loop vectorized using 64 byte vectors
.\simd2.cpp:6:14: optimized: basic block part vectorized using 64 byte vectors
.\simd2.cpp:6:14: optimized: basic block part vectorized using 64 byte vectors
.\simd2.cpp:11:23: optimized: loop vectorized using 64 byte vectors
.\simd2.cpp:12:14: optimized: basic block part vectorized using 64 byte vectors

と表示されます。512-bit（64-byte）の範囲でベクトル化が行われていることがわかります。

では、Eのコマンドでコンパイルしたアセンブラ出力の一部（add関数）を見てみましょう。

_Z3addPKfS0_Pf:
.LFB1860:
    .seh_endprologue
    vmovups (%rcx), %ymm1
    vmovups 32(%rcx), %ymm2
    vaddps  (%rdx), %ymm1, %ymm0
    vmovups 64(%rcx), %ymm3
    vmovups 96(%rcx), %ymm4
    vmovups 128(%rcx), %ymm5
    vmovups 160(%rcx), %ymm1
    vmovups %ymm0, (%r8)
    vaddps  32(%rdx), %ymm2, %ymm0
    vmovups 192(%rcx), %ymm2
    vmovups %ymm0, 32(%r8)
    vaddps  64(%rdx), %ymm3, %ymm0
    vmovups 224(%rcx), %ymm3
    vmovups %ymm0, 64(%r8)
    vaddps  96(%rdx), %ymm4, %ymm0
    vmovups 256(%rcx), %ymm4
    vmovups %ymm0, 96(%r8)
    vaddps  128(%rdx), %ymm5, %ymm0
    vmovups 288(%rcx), %ymm5
    vmovups %ymm0, 128(%r8)
    vaddps  160(%rdx), %ymm1, %ymm0
    vmovups 320(%rcx), %ymm1
    vmovups %ymm0, 160(%r8)
    vaddps  192(%rdx), %ymm2, %ymm0
    vmovups %ymm0, 192(%r8)
    vaddps  224(%rdx), %ymm3, %ymm0
    vmovups %ymm0, 224(%r8)
    vaddps  256(%rdx), %ymm4, %ymm0
    vmovups %ymm0, 256(%r8)
    vaddps  288(%rdx), %ymm5, %ymm0
    vmovups %ymm0, 288(%r8)
    vaddps  320(%rdx), %ymm1, %ymm0
    vmovups %ymm0, 320(%r8)
    vmovups 352(%rcx), %ymm2
    vaddps  352(%rdx), %ymm2, %ymm0
    vmovups %ymm0, 352(%r8)
    vmovups 384(%rdx), %xmm0
    vaddps  384(%rcx), %xmm0, %xmm0
    vmovups %xmm0, 384(%r8)
    vzeroupper
    ret
    .seh_endproc
    .def    __main; .scl    2;  .type   32; .endef
    .section    .text.startup,"x"
    .p2align 4
    .globl  main
    .def    main;   .scl    2;  .type   32; .endef
    .seh_proc   main

これを見たところvmovupsとvaddpsという命令が使用されています。vmovups命令では、8つの32-bitの値（float）を256-bitのYMMレジスタにデータをロードし、vaddpsでは、その8つの値が格納されたYMMを足し合わせていくというものです。

YMMには8つの値を格納しているので、vaddps命令一つで、8つの演算を同時に処理する事ができ、SIMD化されていますね。

このような形でベクトル化させる場合、ソースコード側の配慮も必要なので、注意しましょう。

締め

結果ではなく手段の話でした。なので、実行時間の違いというのはありません。申し訳ないです。

しかし、使えるならSIMD演算を用いると、シングルコアでも並列度が上がって高速化することがあるかもしれませんのでおすすめです。

以下、参考を載せておきます。基本的に、Intel C Compiler向けの内容ですが、GCCやClang、なんならVCでも使えます。

参考

https://www.intel.com/content/www/us/en/docs/intrinsics-guide/index.html

https://jp.xlsoft.com/documents/intel/compiler/19/cpp_19_win_lin/index.htm#GUID-014B2DB6-B363-4CEB-97EC-74CD6A018106.html

https://learn.microsoft.com/ja-jp/cpp/intrinsics/x86-intrinsics-list?view=msvc-170