いままでIntelのCPU(x86)でSSEやAVXといったSIMDについては知っていたけれどARM系のCPUでもNEONと呼ばれるSIMDが使えるらしいことを知りました。手元に手頃なARMマシンであるラズパイがあるので勉強もかねて試してみたいと思います。
SIMDって？NEONって？
まずは基本的な用語をまとめてみたいと思います。
まずSIMDとはSingle Instruction Multiple Dataの略でその名の通り1回の命令で複数のデータを処理します。1命令で4つの値を同時に足し算するようなイメージです。例えばx = a + b, y = c + dのような計算が1命令でできます。
NEONとはARMのWEBサイトにもあるようにARMのSIMD拡張です。ARM系のプロセッサでSIMDを使うときに使用するアーキテクチャのことだと思っていいと思います。
ちなみにSIMDをサポートするIntel系(x86, x86_64)のアーキテクチャがSSEやAVXと呼ばれます。
ラズパイ3とNEON
ARM系プロセッサを搭載したマシンとして手元にラズパイ3があるのでこれを使用してNEONの動きを見てみたいと思います。
まずラズパイ3について簡単にまとめるとCPU(SoC)はBroadcom BCM2837を搭載しており、これはCortex-A53というコアを使用しています。Cortex-A53はARMv8-A64ビットの命令セットを備えています。ARMv7ではNEONはオプションでしたがARMv8では全てのコアがNEONをサポートしています。よって ラズパイ3では(ハードウェア的には)NEONを使用することができます。
NEONを使うコードを書いてみる
それでは実際にNEONを使用するコードを書いてみたいと思います。直接アセンブリで書くのは大変すぎるのでC言語でコードを書き、生成されたアセンブリを読むことで動きを確認し、実際に動かしてパフォーマンスを見てみるといった流れで進めていきます。
今回使用する環境はRaspberry Pi 3 Model BにGentoo Linuxを入れてあり、64ビットのカーネルが動いています。使用するコンパイラはgcc 12.2.1です。
NEONを使用した並列処理を活用できる計算としてベクトルの内積計算を例に行ってみます。ここでどの型で値を表現するかが重要なのですが、今回は単精度の浮動小数点数(32bitのfloat型)で表されるベクトルの計算を扱ってみます。
簡単のためにベクトルの要素数は100で固定とし、内積計算をするdot()関数を呼ぶプログラムを作成します。100としたのは4で割り切れるとメモリのアドレス的に嬉しいだろうという思いから4の倍数で適当な値を選びました。下記のようなmain.cを作成しました。
include <stdio.h>
float dot(float *vec1, float *vec2);
int main(int argc, char *argv[])
{
float v1[100];
float v2[100];
for(int i = 0; i < 100; i++) {
v1[i] = (float) i;
v2[i] = (float) i;
}
float result = dot(v1, v2);
printf("%f\n", result);
return 0;
}このプログラムから呼び出される、内積計算をする関数をdot.cに作成します。
float dot(float *vec1, float *vec2)
{
float s = 0.0f;
for(int i = 0; i < 100; i++) {
s += vec1[i] * vec2[i];
}
return s;
}そしてこれらをビルドします。
gcc -mcpu=cortex-a53 -O0 -c -o dot.o dot.c
gcc -mcpu=cortex-a53 -O0 -c -o main.o main.c
gcc -mcpu=cortex-a53 -o main dot.o main.oビルドされた実行ファイルmainを実行してみます。
$ ./main
328350.000000特になんてこともない、普通のCプログラムを実行した感じですね。意図したとおり内積の結果が得られました。この時どのような命令で動いているかを確認するためにアセンブリを見てみたいと思います。
アセンブリを得るためにはgccに-Sオプションをつけます。今回興味があるのは内積計算を行うdot.cなのでこれをアセンブリにコンパイルしてみます。
gcc -mcpu=cortex-a53 -O0 -S -o dot.s dot.c生成されたdot.sの中身を見てみます。ファイルの先頭のメタデータなど処理と直接関係しない場所は省略しています。
dot:
.LFB4350:
.cfi_startproc
sub sp, sp, #32
.cfi_def_cfa_offset 32
str x0, [sp, 8]
str x1, [sp]
str wzr, [sp, 24]
str wzr, [sp, 28]
b .L2
.L3:
ldrsw x0, [sp, 28]
lsl x0, x0, 2
ldr x1, [sp, 8]
add x0, x1, x0
ldr s1, [x0]
ldrsw x0, [sp, 28]
lsl x0, x0, 2
ldr x1, [sp]
add x0, x1, x0
ldr s0, [x0]
fmul s0, s1, s0
ldr s1, [sp, 24]
fadd s0, s1, s0
str s0, [sp, 24]
ldr w0, [sp, 28]
add w0, w0, 1
str w0, [sp, 28]
.L2:
ldr w0, [sp, 28]
cmp w0, 99
ble .L3
ldr s0, [sp, 24]
add sp, sp, 32
.cfi_def_cfa_offset 0
ret
.cfi_endprocまずは太字の箇所で、ループ用の変数[sp, 28]にゼロをセットして、インクリメントしながら99を超えるまで処理するループがあります。
dot:
.LFB4350:
.cfi_startproc
sub sp, sp, #32
.cfi_def_cfa_offset 32
str x0, [sp, 8]
str x1, [sp]
str wzr, [sp, 24]
str wzr, [sp, 28]
b .L2
.L3:
ldrsw x0, [sp, 28]
lsl x0, x0, 2
ldr x1, [sp, 8]
add x0, x1, x0
ldr s1, [x0]
ldrsw x0, [sp, 28]
lsl x0, x0, 2
ldr x1, [sp]
add x0, x1, x0
ldr s0, [x0]
fmul s0, s1, s0
ldr s1, [sp, 24]
fadd s0, s1, s0
str s0, [sp, 24]
ldr w0, [sp, 28]
add w0, w0, 1
str w0, [sp, 28]
.L2:
ldr w0, [sp, 28]
cmp w0, 99
ble .L3
ldr s0, [sp, 24]
add sp, sp, 32
.cfi_def_cfa_offset 0
ret
.cfi_endproc次は、ループ用の変数[sp, 28]から値をとってきて4倍(2ビット分左シフト)した値を[sp, 8]に加算しています。ここで[sp, 8]は引数のベクトルの配列の先頭アドレスだと思われます。同じことを[sp]に対しても行います。[sp]は2つ目の引数のアドレスだと思います。そうして算出したアドレスからs0とs1レジスタに値を読み込んでいます。s0とs1レジスタは32bitの浮動小数点用のレジスタでちょうどfloat型の値1つ分です。
dot:
.LFB4350:
.cfi_startproc
sub sp, sp, #32
.cfi_def_cfa_offset 32
str x0, [sp, 8]
str x1, [sp]
str wzr, [sp, 24]
str wzr, [sp, 28]
b .L2
.L3:
ldrsw x0, [sp, 28]
lsl x0, x0, 2
ldr x1, [sp, 8]
add x0, x1, x0
ldr s1, [x0]
ldrsw x0, [sp, 28]
lsl x0, x0, 2
ldr x1, [sp]
add x0, x1, x0
ldr s0, [x0]
fmul s0, s1, s0
ldr s1, [sp, 24]
fadd s0, s1, s0
str s0, [sp, 24]
ldr w0, [sp, 28]
add w0, w0, 1
str w0, [sp, 28]
.L2:
ldr w0, [sp, 28]
cmp w0, 99
ble .L3
ldr s0, [sp, 24]
add sp, sp, 32
.cfi_def_cfa_offset 0
ret
.cfi_endprocそれからs0とs1の内容をfmulで掛け算し、結果を総和を保管する変数[sp, 24]に保存します。最後にループのカウンタ[sp, 28]の値をインクリメントしてループ1回分の処理が終わります。
コンパイラが生成したアセンブリを眺めてきました。これは最適化なし(-O0オプション)のものなのでC言語との対応づけがわかりやすく簡単に理解できたと思います。次は強い最適化(-O3オプション)を行った処理と比較してみたいと思います。
gccに最適化オプション-O3を指定してアセンブリを出力させます。
gcc -mpu=cortex-a53 -O3 -S -o dot.s dot.c生成されたアセンブリdot.sを見てみます。
dot:
.LFB4350:
.cfi_startproc
movi v0.2s, #0
mov x2, 0
.p2align 3,,7
.L2:
ldr q2, [x0, x2]
ldr q1, [x1, x2]
add x2, x2, 16
cmp x2, 400
fmul v1.4s, v1.4s, v2.4s
fadd s0, s0, s1
dup s3, v1.s[1]
dup s2, v1.s[2]
dup s1, v1.s[3]
fadd s3, s3, s0
fadd s0, s2, s3
fadd s0, s0, s1
bne .L2
ret
.cfi_endproc様子がガラッと変わりとても短くなりました。なんとNEONを使用した並列化まで行われています。詳しく見ていきます。
まず引数のベクトルのポインタ2つからレジスタに値を読み込んでいます。ここで指定しているレジスタはq1とq2です。このq1とq2というレジスタは128bitのレジスタで32bitのfloat型4つ分です。次はポインタを読むオフセットに16を加算しています。128bit読み込んだので128bit/8=16byteとなっています。そしてオフセットと400を比較しています。16ずつ加算して要素数100/4=25回のループが終了するかの条件判断です。ジャンプ命令bneはずっと下にありますが途中で比較したフラグが変化しないためここで比較しています。
次は”fmul v1.4s, v1.4s, v2.4s”となっています。これがNEONを使用した並列の計算です！さきほど値をセットしたq1, q2レジスタとv1, v2レジスタは同じレジスタを示しています。なのでv1にはq1の値が入っています。v1の128bitのレジスタを32bitずつ4つの値と解釈するように伝えているのがv1.4sという書き方です。なのでこれは”v1 = v1 * v2″という計算していることになります。
その下のfaddとdup命令が7つ続いている処理は計算したv1に入っている4つのfloat型の値を順に足し合わせて(総和を求めて)合計を保存しているレジスタs0に値を入れています。ちなみにs0, …, s3は32bitの浮動小数点を扱うレジスタで、128bitのv0レジスタと共有しています。v0レジスタはプログラム先頭で0で初期化されていますのでそのタイミングでs0=0となっています。
組み込み関数でNEONを使用する
ここまでコンパイラの最適化で自動的に生成されたNEONを使用する処理を見てきましたが、C言語から直接NEONの処理を扱える組み込み関数があるのでそれを使用してみたいと思います。
さきほどのdot.cを下記のように書き換えます。
#include <arm_neon.h>
float dot(float *vec1, float *vec2)
{
float s = 0.0f;
for(int i = 0; i < 100 / 4; i++) {
float32x4_t v1, v2, v3;
v1 = vld1q_f32(vec1 + i * 4);
v2 = vld1q_f32(vec2 + i * 4);
v3 = vmulq_f32(v1, v2);
float32_t sum = vaddvq_f32(v3);
s += sum;
}
return s;
}まずC言語からNEONの組み込み関数を扱うためにarm_neon.hをincludeします。
そして今回は128bitのレジスタを使用してfloat型を4つ同時に処理したいのでfloat32x4_tという型でデータを扱います。このfloat32x4_tで値を読み込む関数はvld1q_f32()です。ベクトル2つ分の値を読み込んだらfmulq_f32()関数で要素ごとの積を求めます。
ここまでは先ほど見たアセンブリのコードと同じ処理ですが次に水平加算命令を使用します。水平加算命令とは1命令でレジスタにある複数の値の和を求めるものです。vaddvq_f32()関数を使用してfloat32x4_t型の4つの総和を求めてfloat32_t型の値にします。これは標準のfloatと同じ32bitの浮動小数点なので合計を保存する変数sに足しあわせます。
このプログラムをアセンブリにコンパイルしてみます。最適化を弱めて-O1としてみました。
gcc -mcpu=cortex-a53 -O1 -S dot.c生成されたコードがこちらです。
dot:
.LFB4351:
.cfi_startproc
mov x2, 0
movi v0.2s, #0
.L5:
ldr q1, [x0, x2]
ldr q2, [x1, x2]
fmul v1.4s, v1.4s, v2.4s
faddp v1.4s, v1.4s, v1.4s
faddp v1.4s, v1.4s, v1.4s
fadd s0, s0, s1
add x2, x2, 16
cmp x2, 400
bne .L5
ret
.cfi_endproc最適化-O1のレベルでは自動的に生成されないNEON命令が組み込み関数を使用することでしっかりと使われています。また、先ほどあった大量のfaddとdup命令がすっきりとしました。
ただ、組み込み関数で指定したfaddv命令が使われず、faddp命令2回になっているところが腑に落ちません。faddvは4つのfloat型の値の総和を1回で求めますがfaddpは4つのうち隣り合う2つ和を求めるため4つの総和を出すためには2回の命令が必要になってしまいます。
faddpのイメージとしては128bitのfloat value[4]に対して、
1回目で、value[0] = value[0] + value[1], value[1] = value[2] + value[3]とし、
2回目で、value[0] = value[0] + value[1]となって総和が求まるイメージです。
なぜこのようなコンパイル結果になったかはよく調べてみる必要がありそうです。
計算速度比較
最後にNEONを使用して計算速度がどの程度高速化するのか時間を測って比較してみたいと思います。
単純なループの処理、コンパイラの最適化で自動的にNEONを使用する、組み込み関数を使用してNEONを使うのプログラムを比較します。
時間計測のためにプログラムに計測用コードを追加します。下記のようにmain.cを変えました。太字が変更点です。
#include <stdio.h>
#include <time.h>
float dot(float *vec1, float *vec2);
int main(int argc, char *argv[])
{
struct timespec ts_start, ts_end;
float v1[100];
float v2[100];
for(int i = 0; i < 100; i++) {
v1[i] = (float) i;
v2[i] = (float) i;
}
timespec_get(&ts_start, TIME_UTC);
float result = dot(v1, v2);
timespec_get(&ts_end, TIME_UTC);
printf("time: %lu\n", ts_end.tv_nsec - ts_start.tv_nsec);
printf("%f\n", result);
return 0;
}これで時間計測を行います。計測は100回実行して平均値を算出しました。小数点以下切り捨てです。
| 条件 | 最適化レベル | 処理時間[ナノ秒] | 備考 |
| 単純なループ実装 | -O0 | 3703 | 最適化なし |
| NEON組み込み関数使用 | -O0 | 2721 | |
| 単純なループ実装 | -O1 | 1327 | 不要な変数が削除された |
| NEON組み込み関数使用 | -O1 | 810 | |
| 単純なループ実装 | -O2 | 980 | NEON使用コードが生成された |
| NEON組み込み関数使用 | -O2 | 788 | |
| 単純なループ実装 | -O3 | 950 | |
| NEON組み込み関数使用 | -O3 | 764 | NEON使用コードが生成された |
ちなみに計測した時のコードはこんな感じです。
for n in $(seq 1 100); do ./main | head -1; done | awk '{a+=$2}END{print a/NR}'注意点としては、ナノ秒単位の計測ですが、ラズパイでどの程度の精度で計測できているかが不明なのと、時間経過の算出で引き算していますが繰り上がりを考慮していないので外れ値が混ざる可能性があることが挙げられます。
結果を見るとNEON使用によって大きく高速化していることがわかります。また、それ以上に最適化で無駄な処理を削る方が大きく時間に影響しています。これは今回は100次元のベクトルの計算なので計算量が少なく、他の処理の方が大きく影響しているためだと考えられます。計算の高速化を比較するのであればベクトルの次元数を100から10000などに100倍ほど増やしてみると良いかもしれません。
まとめ
ARMのSIMDであるNEONを使用して計算を並列化して高速化してみる遊びができました。整数値を扱ってみることや掛け算以外の命令を使ってみるなど今回は触れなかったことがたくさんあります。ですが基本的な感覚はつかめたので満足です。機会があれば深堀りしていきます。
コメント