LoginSignup
6
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

OECUAdvent Calendar 2024

Day 10
@nkiib

【C++/x86】SIMDを触ってみようぜ

Last updated at Posted at 2024-12-09

1日目は、精度を高める方向のお話をしました。今回は、逆...とは言いませんが、精度を落としたうえで計算効率を高める内容です。

SIMDって？

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

軽くイメージを。

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

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

image.png

これは、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のこれを見ると良いです。

自動でベクトル・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化されていますね。

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

この節の記述において、@fujitanozomu様のご指摘を参考にさせていただきました。この場を借りてお礼申し上げます。

締め

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

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

アドベントカレンダーの担当日までにAVX-512の環境を用意できなかったので、用意できたら内容が充実するので、また数カ月後に見に来てくれると嬉しいです。

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

掲載当初、Int2を2-bitとしておりましたが、間違いです。Int2は2-Byteでした。

参考

6
1
8

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
6
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?