自動ベクトル化とC言語のrestrict
単純なループ
与えられた2つの float の配列の要素ごとの和を計算する関数を考えます。
#include <stddef.h>
void add(size_t n, float *result, const float *a, const float *b) {
for (size_t i = 0; i < n; ++i) {
result[i] = a[i] + b[i];
}
}
C99の可変長配列 (VLA) を使うと、同等の関数を次のように書くこともできます。
#include <stddef.h>
void add_vla(size_t n, float result[n], const float a[n], const float b[n]) {
for (size_t i = 0; i < n; ++i) {
result[i] = a[i] + b[i];
}
}
SIMD intrinsicsの使用
最近のCPUにはSIMD命令が載っています。これらのSIMD命令を使って先ほどのコードを書き直してみましょう。
Intel SSEの場合：
#include <stddef.h>
#include <immintrin.h>
void add_sse(size_t n, float *result, const float *a, const float *b) {
size_t i = 0;
for (; i + 4 <= n; i += 4) {
__m128 va = _mm_loadu_ps(&a[i]);
__m128 vb = _mm_loadu_ps(&b[i]);
__m128 vr = _mm_add_ps(va, vb);
_mm_storeu_ps(&result[i], vr);
}
#pragma clang loop vectorize(disable)
for (; i < n; ++i) {
result[i] = a[i] + b[i];
}
}
Arm NEON/ASIMDの場合：
#include <stddef.h>
#include <arm_neon.h>
void add_neon(size_t n, float *result, const float *a, const float *b) {
size_t i = 0;
for (; i + 4 <= n; i += 4) {
float32x4_t va = vld1q_f32(&a[i]);
float32x4_t vb = vld1q_f32(&b[i]);
float32x4_t vr = vaddq_f32(va, vb);
vst1q_f32(&result[i], vr);
}
#pragma clang loop vectorize(disable)
for (; i < n; ++i) {
result[i] = a[i] + b[i];
}
}
コンパイラーの自動ベクトル化
最近のコンパイラーには自動ベクトル化という機能が備わっていて、最初のコードなら余裕でベクトル化してくれます（自動ベクトル化）。実際に生成されたコードを見てみましょう（生成コードを簡単にするためにloop unrollingは無効にしています）。
同じようなコードが生成されるかと思いきや、自動ベクトル化させた方が生成コードが長くなっています。なぜでしょうか？
実は元々のコードと手動ベクトル化したものでは、関数の意味が微妙に変わっています。具体的には、一部がオーバーラップした配列を与えた際の挙動が異なります。
#include <stdio.h>
int main(int argc, char *argv[]) {
{
float a[4] = {0.0f, 2.0f, 4.0f, 6.0f};
float b[5] = {0.0f, -1.0f, -2.0f, -3.0f, -4.0f};
add_p(4, b + 1, a, b);
for (int i = 0; i < 5; ++i) {
printf("%g\n", b[i]);
}
}
puts("---");
{
float a[4] = {0.0f, 2.0f, 4.0f, 6.0f};
float b[5] = {0.0f, -1.0f, -2.0f, -3.0f, -4.0f};
add_neon(4, b + 1, a, b); // or add_sse
for (int i = 0; i < 5; ++i) {
printf("%g\n", b[i]);
}
}
}
出力：
0
0
2
6
12
---
0
0
1
2
3
コンパイラーは最適化の際にコードの意味を変えてしまうといけないため、自動ベクトル化させた方は配列のオーバーラップの検査が入っているようです。
restrictキーワードの使用
入出力の配列がオーバーラップすることはないとコンパイラーに伝えることができれば、自動ベクトル化の際により簡潔なコードを出力できるようになります。C言語にはそのためのキーワード、 restrict があります。ポインターの場合は restrict は * の後に書きます。
#include <stddef.h>
void add_r(size_t n, float * restrict result, const float *a, const float *b) {
for (size_t i = 0; i < n; ++i) {
result[i] = a[i] + b[i];
}
}
配列の場合は、やや奇妙な文法ですが、 [] の中、要素数の前に書きます。
#include <stddef.h>
void add_r_vla(size_t n, float result[restrict n], const float a[n], const float b[n]) {
for (size_t i = 0; i < n; ++i) {
result[i] = a[i] + b[i];
}
}
生成されたコードを見ると、 restrict によって生成コードが短くなっていることがわかります：
restrict キーワードにはC++にはありませんが、各種コンパイラーが独自拡張として似たようなキーワードを提供している場合があるようです。
余談：__builtin_assume_aligned の使用
IntelのSSE/AVX系には、ポインターが16バイト（SSEの場合）／32バイト（AVXの場合）／64バイト（AVX-512の場合）にアラインされている場合に使えるロード・ストア命令（movaps）があります。元々のコードはポインターのアラインメントに仮定を置いていないためコンパイラーはベクトル化の際にそれらを使うことができませんが、コンパイラーにアラインメントの仮定を教えてやるとベクトル化の際にmovapsが使われるようになります。
GCC/Clangの場合は __builtin_assume_aligned を使うことでポインターのアラインメントをコンパイラーに教えてやることができます。
#include <stddef.h>
void add_aligned(size_t n, float * restrict result, const float *a, const float *b) {
result = __builtin_assume_aligned(result, 16); // 16バイトアラインされていることを仮定
a = __builtin_assume_aligned(a, 16);
b = __builtin_assume_aligned(b, 16);
for (size_t i = 0; i < n; ++i) {
result[i] = a[i] + b[i];
}
}
配列に対しても使えます。
#include <stddef.h>
void add_aligned_vla(size_t n, float result[restrict n], const float a[n], const float b[n]) {
result = __builtin_assume_aligned(result, 16);
a = __builtin_assume_aligned(a, 16);
b = __builtin_assume_aligned(b, 16);
for (size_t i = 0; i < n; ++i) {
result[i] = a[i] + b[i];
}
}
C++20には標準化された機能として std::assume_aligned があるようです。
Discussion