Tweet

以下,日記エントリからの include です. 書いてることが重複してたりするのは勘弁してね.

gcc-4.x での SIMD 命令サポート

「gcc 4 では MMX, 3DNow!, SSE などの SIMD 命令がサポートされた」「-msse オプションを付けてコンパイルすると SSE 命令が有効になる」なんて話を聞いたことがある人も多いと思います. が,これだけではマトモに SIMD 命令を生成してくれません

というわけで,そこらへんの話を調べた範囲で書いてみようと思います.

情報源

gcc の info がメインです. あと,各命令セットのドキュメントも用意しときましょう.

おおまかな枠組み

  1. SIMD のレジスタにあった型の変数を定義して
  2. ビルトイン関数で演算を記述
  3. -mXXX オプションを付けてコンパイル
    • -mmmx, -m3dnow, -msse, -msse2 など
  • 枠組みとしては汎用
    • x86 だけでなく,MIPS や XScale の wireless MMX もこの方法でサポートしてるそうです

ま,こんなところです.

3DNow! の例

何でいまどき 3DNow! かと言うと,今使ってる PC が Eden だったりするので. まぁ,SSE とかでも同じような要領でいけるはずです.

今回は,単精度浮動小数点の積和演算のプログラムを組んでみようと思います.

3DNow! 命令セット

マニュアルは AMD の web からダウンロードできます.

3DNow! 命令の作法

3DNow! では,MMX と同様に 80387 浮動小数点レジスタを利用します. が,好き勝手に 3DNow! 命令を使うと,80387 が混乱してしまいます.

というわけで,3DNow! 命令を使うセクションの前後には femms 命令でレジスタを初期化する必要があります. つまり

80387演算
femms
3DNow! 演算(80387 利用不可)
…
femms
80387演算

としなくてはなりません.

SSE 以降では専用レジスタがあるので,このような切り替えの作法は不要だと思います.

このような型を宣言します.

typedef float v2sf __attribute__ ((vector_size(8)));

float が 2 つくっついた型で,レジスタ 1 つ分に相当します.

宣言と代入

要するに typedef された型なので,宣言自体は

v2sf a;

でできます. 宣言と同時に初期値を代入する場合は

v2sf a = { 1.0, 2.0 };

と2つの実数を {} で囲んで渡します.

代入する時は

a = (v2sf){ 1.0, 2.0 };

とします.

値の取り出し

…どーやるんだ,これ? info を見ても書いてないようです.

float f;
f = a[1];

ではコンパイラにはねられます. が,gdb 上では有効だったりします.

しかたがないので,とりあえず

float f[2];
*(v2sf *)f = a;

とかやって,強引に取り出しています.

ちなみに

float *f = (float *)&a;

でも取り出すことはできます. この場合,「a のアドレス」を参照してしまっているので,a を register v2sf 型で宣言できなくなってしまいます. つまり,最適化の妨げになる可能性があるので,やらないほうがいいでしょう.

演算

gcc 上では __builtin_… 命令を使います.

今回使用したのは

v2sf __builtin_ia32_pfadd ( v2sf a, v2sf b )
a + b
v2sf __builtin_ia32_pfmul ( v2sf a, v2sf b )
a x b
v2sf __builtin_ia32_pfacc ( v2sf a, v2sf b )
{ a(上位) + a(下位), b(上位) + b(下位) }
void __builtin_ia32_femms()
80387 レジスタの初期化

です.

プログラム

ただ計算するだけでもつまんないので,80387 演算との速度を比較できるようにしてみました.

実行結果

コンパイルオプションを色々変えて試してみました. 環境は Eden 800MHz 上の Linux です. まずは,素で.

sh-3.2$ gcc -m3dnow benchmark.c 
sh-3.2$ ./a.out 
3DNow!:         11.700000       3998000.000000
387 double:     11.180000       3998000.000000
387 float:      11.160000       3998000.000000

2番目のカラムが実行時間,一番右が実行結果です. 3DNow!,全然速くないですね.それでは,最適化レベルをあげていってみましょう. まずは O1

sh-3.2$ gcc -O1 -m3dnow benchmark.c 
sh-3.2$ ./a.out 
3DNow!:         1.270000        3998000.000000
387 double:     3.070000        3998000.000000
387 float:      2.540000        3998000.000000

さらに O2

sh-3.2$ gcc -O2 -m3dnow benchmark.c 
sh-3.2$ ./a.out 
3DNow!:         1.280000        3998000.000000
387 double:     3.050000        3998000.000000
387 float:      2.570000        3998000.000000

O1 と変化無いですね. ループを展開してみましょう.

sh-3.2$ gcc -funroll-loops -O2 -m3dnow benchmark.c 
sh-3.2$ ./a.out 
3DNow!:         0.890000        3998000.000000
387 double:     3.050000        3998000.000000
387 float:      2.540000        3998000.000000

こいつは効果があるようです. 80387 命令比で 2.8 倍の速度で計算できてます.

gcc 4.x の 3DNow! 命令サポートについての考察

日記/2008-09-03/gcc-4.x での SIMD 命令サポートでは,3DNow! 命令セットで簡単なベンチマークプログラムを書いて,x87 命令に対しパフォーマンスアップが確認できました.

今回は,もう少し突っ込んで見ていきます.

v2sf 型の代入

-m3dnow オプションを付けると,実は __builtin__ia32_… 関数を明示的に使わないところでも拡張命令が使用されてしまいます. v2sf 型変数への代入部分です.

例えば

v2sf a;
a = (v2sf){ 1.0, 2.0 };

というコードがあると,ここで movq 命令が使用されます(movq は本当は MMX 命令なのですが,gcc ではこれも 3DNow! 命令セットの一部とみなしてしまうようです). したがって

float f = 3.0;
…
a = (v2sf){ f + 1.0, 2.0 };

は,x87 命令と 3DNow! 命令の「混ぜたら危険」ルールに反することになります. 本来は gcc が良きに計らって欲しいところですが,残念ながらケアしてくれないようです.

ループ展開の効果

ベンチマークの最後で -funroll-loops オプションを付けた場合の実行結果を示しました. 3DNow! 命令ではパフォーマンスが上昇したのに対し,x87 命令ではパフォーマンスはほとんど差がありませんでした. なぜでしょう?

CPU の中を覗いてきたわけではないので断言はできないのですが,ある程度想像はつきます. x87 では,中の構造はスタックマシン,要するに逆ポーランド電卓です. スタックに値を詰め込み,スタックトップに対して演算を指示することにより計算を行います. 複数の演算を続けて行う場合,前の演算結果に依存することになります. このため,パイプラインによる並列処理の恩恵を受けにくくなっているのではないかと想像できます.

一方,3DNow! では普通のレジスタマシンです. 8 本の MMX レジスタを自在に使って計算を行います. 実際に逆アセンブルしてみたところ,ループ展開を行うと,複数の MMX レジスタを使って計算を行うコードが生成されていました. プログラムの並列性が上がったために複数の演算が同時進行することになり,パフォーマンスが上がったのではないかと思われます.

gcc-4.x の SIMD 命令サポート その2

新しいマシンを買い,Linux マシンが Celeron D にグレードアップしてしまいました. というわけで,今回は SSE/SSE2 について.

SSE レジスタ

SSE では 128 ビット長の専用レジスタを 8 本新設しています. つまり,前回の 3DNow! とは違い,x87 演算と分離する必要はありません.

で,レジスタの中身ですが,SSE レジスタ 1 本に

  • SSE では 32bit 単精度 (float) 型を 4 本
  • SSE2 では 64bit 倍精度 (double) 型を 2 本

を格納し,SIMD 演算を行うことができます.

参考文献

参考となる書籍など.

おおまかな枠組み

前回と一緒で

  1. SIMD のレジスタにあった型の変数を定義して
  2. ビルトイン関数で演算を記述
  3. -mXXX オプションを付けてコンパイル

です.

プログラム例

前回と同様にベンチマークを兼ねてプログラムを作ってみました.

お題も前回と同じく積和演算です.

型宣言

まずは SSE での「単精度 x 4 本」の構成の変数の型宣言から.

typedef float v4sf __attribute__ ((vector_size(16))) __attribute__((aligned(16)));

前回と違っておしりに

__attribute__((aligned(16)))

なんて属性がついてますね. これは

変数を 16 で割り切れるアドレスに割り当てよ

という属性です. これは,実は SSE レジスタのロード・ストア時においてこの様な制限があるために,このような属性を追加しています. この制限を無視するとどうなるかというと,Segmentation fault が発生します. ARM や MIPS などの変数のアラインメントと似たような感じです.

ここで 1 つ注意.この aligned(16) 属性は

v4sf a;

のような変数の宣言では有効に機能しますが

v4sf *a = malloc ( sizeof(v4sf) );

のような動的メモリ割り当てではうまく働きません.また

float a[2];
*(v4sf *)&a[0] = b;

のようなキャストもだめです. つまり,既に割り当てられている領域を「16 で割り切れるアドレスに配置換え」する力はありません.

ついでに「倍精度 x 2」の型宣言.

typedef double v2df __attribute__ ((vector_size(16))) __attribute__((aligned(16)));

制限・注意は v4sf 型と一緒です.

SSE/SSE2 命令

まずは v4sf 型について. これらは SSE 命令セットに属します.

v4sf __builtin_ia32_addps (v4sf a, v4sf b)
v4sf 型に対する SIMD の足し算です.a + b の結果を返します.
v4sf __builtin_ia32_mulps (v4sf a, v4sf b)
v4sf 型に対する SIMD 掛け算です.a * b の結果を返します.

v2df 型について.これらは SSE2 命令セットです.

v2df __builtin_ia32_addpd (v2df a, v2df b)
v2df 型に対する SIMD 足し算.a + b の結果を返します.
v2df __builtin_ia32_mulpd (v2df a, v2df b)
v2df 型に対する SIMD 掛け算.a * b の結果を返します.

これらの明示的なビルトイン関数以外にも例えば

v2sf a = { 1.0, 2.0, 3.0, 4.0 };

のような代入文でも SSE 命令が使用されます. つまり,単なる代入においても右辺と左辺の変数領域は 16 で割り切れるアドレスでなければなりません. うっかりすると忘れてしまうので注意しましょう.

ベンチマーク環境

結果に先立って,ベンチマーク環境について書いときます.

$ cat /proc/cpuinfo
processor	: 0
vendor_id	: GenuineIntel
cpu family	: 15
model		: 4
model name	: Intel(R) Celeron(R) CPU 2.66GHz
stepping	: 9
cpu MHz		: 2667.096
cache size	: 256 KB
fdiv_bug	: no
hlt_bug		: no
f00f_bug	: no
coma_bug	: no
fpu		: yes
fpu_exception	: yes
cpuid level	: 5
wp		: yes
flags		: fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts \
acpi mmx fxsr sse sse2 ss ht tm pbe lm constant_tsc pebs bts pni monitor ds_cpl tm2 cid cx16 xtpr \
lahf_lm
bogomips	: 5340.17
clflush size	: 64
power management:

ベンチマーク結果

まずは最適化オプションなしでこんな感じ.

$ gcc -msse2 benchmark.c
$ ./a.out 
SSE double:	9.620000	3998000.000000
norm double:	20.890000	3998000.000000
SSE float:	4.910000	3998000.000000
norm float:	20.560000	3998000.000000

中央のカラムが演算にかかった CPU 時間で,単位は秒です. 最適化オプションを付けてみましょう.

$ gcc -msse2 -O2 benchmark.c
$ ./a.out 
SSE double:	2.680000	3998000.000000
norm double:	4.910000	3998000.000000
SSE float:	1.070000	3998000.000000
norm float:	4.760000	3998000.000000

かなり違いますね.

ループ展開もしてみましょう.

$ gcc -msse2 -O2 -funroll-loops benchmark.c
imai@JR0BAK:~/progs/sse$ ./a.out 
SSE double:	2.380000	3998000.000000
norm double:	4.850000	3998000.000000
SSE float:	1.060000	3998000.000000
norm float:	4.670000	3998000.000000

ちょっぴり速くなりましたが,3DNow! のときのように劇的な効果はないようです. ループ展開したことで SSE/SSE2 の演算では多くの SSE レジスタを使用するコードが生成されています. にもかかわらずあまり速度が上がってないのは,実は SSE/SSE2 でも演算ユニットは 1 セットしかないのかもしれません.

最後の例で SSE/SSE2 と x87 の速度を比較すると,SSE/SSE2 の SIMD 演算を用いたことにより

  • 単精度で 4.4 倍
  • 倍精度で 2.0 倍

のスピードアップが図られたことになります.

gcc-4.x での SSE/SSE2 浮動小数点演算

前回は,SSE/SSE2 の SIMD 演算機能を使ってみました. が,SIMD 機能を利用するにはプログラムを書き直さなくてはいけません. というわけで,今回はお手軽な方法を.

-mfpmath オプション

gcc-4.x のコンパイルオプションで,double 型や float 型の浮動小数点演算の方法を指定できます.

-mfpmath=387
浮動小数点演算を x87 命令に変換します.
-mfpmath=sse
浮動小数点演算を SSE/SSE2 命令に変換します.-msse や -msse2 オプションも必要です.
-mfpmath=sse,387
浮動小数点演算を x87 と SSE/SSE2 命令の両方に変換します.

例えば

double a,b,c;
c = a + b;

というソースコードを

$ gcc -mfpmath=sse -msse2 foo.c

とオプションでコンパイルしてやると,SSE2 命令を使った実行プログラムが生成されます. つまり,ソースコードの修正無しにリコンパイルするだけでパフォーマンスの向上が期待できることになります.

ベンチマーク

例によってベンチマーク.

前回のプログラムをそのまま使いまわします.

まずは x87 演算.

$ gcc -msse2 -O2 -funroll-loops benchmark.c
$ ./a.out 
SSE double:	2.280000	3998000.000000
norm double:	5.830000	3998000.000000
SSE float:	1.050000	3998000.000000
norm float:	4.700000	3998000.000000

SSE/SSE2 命令を使うようオプションを付けてみます. norm double, norm float の演算が SSE/SSE2 命令でコンパイルされます.

$ gcc -msse2 -O2 -funroll-loops -mfpmath=sse benchmark.c  
$ ./a.out 
SSE double:	2.320000	3998000.000000
norm double:	4.000000	3998000.000000
SSE float:	1.010000	3998000.000000
norm float:	3.790000	3998000.000000
  • 単精度で 1.2 倍
  • 倍精度で 1.5 倍

の速度向上が見られました. が,SSE double, SSE float の値には遠く及びません.

でも,まぁ,コンパイルし直すだけでこれだけ速度が上がれば御の字とも言えるでしょう.

SIMD ではない

というのも,実はこのオプションで生成されるのはスカラ演算命令で,SIMD 命令ではありません. つまり,SSE レジスタを使って x87 のような浮動小数点演算を行う命令を生成するだけで,データ構造を束ねて SIMD 演算してくれるような機能はありません. が,その代わりに例の「16 で割り切れるアドレス」の制限は生じません.

double 型の精度

double 型の場合,メモリ上では 64 bit の浮動小数点数ですが,x87 レジスタ上では 80 bit の浮動小数点で演算されます. このため,演算結果をメモリに格納する時点で有効桁数が落ちることになります. 最適化によってメモリに格納するタイミングが変わったりすると計算結果が変わってしまう,という現象は有名ですね.

が,SSE レジスタ上での倍精度浮動小数点数は 64 bit 長です. つまり,倍精度演算では厳密には x87 演算と結果が異なるかもしれません. でも,倍精度演算の精度は「最適化オプションでも結果が変わる」というレベルなので,大抵のプログラムは SSE2 命令を使ったからといって直ちに誤動作することは無いと思います.


トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2008-09-22 (月) 00:19:36 (3408d)