ナンクル力学系

Cでベクトル演算沢山やるような数値計算をするときに多次元配列を a[i][j] と a[i*Nj+j] の どちらで書くのが速いか気になったので試してみた．(Njは添え字jの数ね．) x 始めは，

• a[i*N+j]で書くよりa[i][j]で書くほうが速いかｗ ってか当たり前かな．毎回i*Nやるの無駄だし． — Twitter / tkf: a[i*N+j]で書くよりa[i][j]で書くほうが …

と思ってたけど，

• ポインタ参照(*)とアドレス取得(&)って四則演算くらい計算コストかかるもんなんだろうか — Twitter / tkf: ポインタ参照(*)とアドレス取得(&)って四則演算く …

とか良くわからないなと思ってたら

• @tkf CPU内部の四則演算より、メモリアクセスの方が、一般に時間はかかるでしょう。キャッシュ (先読み）メモリがあるのは、そのためでしょう。 — Twitter / hsugawa: @tkf CPU内部の四則演算より、メモリアクセスの …

と教えてもらって，一筋縄な問題じゃなさそうだと思ったので．

サンプルで作ったのはのRNN(recurrent neural network)． ソースは gist: 267098 – GitHub にある．結構キモいソースだと思うw．

1次元配列での実装(a[i*Nj+j])を rnn_ca1d.c に， 2次元配列での実装(a[i][j])を rnn_ca2d.c に書いてある．

rnn_ca1d.c は配列にアクセスするためのマクロを:

#define Wcc(i,j) self->wcc[ self->num_c*(i) + (j) ]
#define Bc(i)    self->bc[(i)]
#define Ec(i)    self->ec[(i)]
#define Uc(i,j)  self->uc[ self->num_c*(i) + (j) ]
#define Xc(i,j)  self->xc[ self->num_c*(i) + (j) ]

で書いていて， rnn_ca2d.c は:
#define Wcc(i,j) self->wcc[i][j]
#define Bc(i)    self->bc[i]
#define Ec(i)    self->ec[i]
#define Uc(i,j)  self->uc[i][j]
#define Xc(i,j)  self->xc[i][j]

で書いてある．あとの内容はほとんど同じ．
rnn_ca1d.c と rnn_ca2d.c を gcc と icc でそれぞれ -O2(デフォルトなはず) と -O3 オプションをつけてコンパイル．そして，実行時間を計ってみたのが以下:
tkf% make runtest 2> runtest.txt
for i in  rnn_ca1d-gcc-O2 rnn_ca1d-gcc-O3 rnn_ca1d-icc-O2 rnn_ca1d-icc-O3 rnn_ca
2d-gcc-O2 rnn_ca2d-gcc-O3 rnn_ca2d-icc-O2 rnn_ca2d-icc-O3; \
        do \
        printf "%s %d %d %d " \
        ./$i 30 1000 300 1>&2; \
        time ./$i 30 1000 300; \
        done
tkf% cat runtest.txt
./rnn_ca1d-gcc-O2 30 1000 300 1.65user 0.00system 0:01.65elapsed 99%CPU (0avgtex
t+0avgdata 0maxresident)k
0inputs+0outputs (0major+272minor)pagefaults 0swaps
./rnn_ca1d-gcc-O3 30 1000 300 1.56user 0.00system 0:01.57elapsed 99%CPU (0avgtext+0avgdata 0maxresident)k
(略)

結果が分かりにくいので整形してみる:
tkf% sed -e "N; s/\n//" runtest.txt | cut -f1,5 -d' '
./rnn_ca1d-gcc-O2 1.64user
./rnn_ca1d-gcc-O3 1.56user
./rnn_ca1d-icc-O2 0.80user
./rnn_ca1d-icc-O3 0.80user
./rnn_ca2d-gcc-O2 1.48user
./rnn_ca2d-gcc-O3 2.62user
./rnn_ca2d-icc-O2 1.07user
./rnn_ca2d-icc-O3 2.23user

分かることは，

icc で1次元配列(a[i*Nj+j])で実装したバージョンが一番速い．
2次元配列(a[i][j])での実装は，gcc/iccのどちらでも -O3 の パフォーマンスが落ちる．不思議．
gcc -O2 だと，2次元配列(a[i][j])のほうが若干速い．
※ 3回試してほとんど同じ結果だった．

ちなみに， make の時に出たメッセージを見ると2次元配列(a[i][j])での実装を-O3で コンパイルするとベクトル化されてるループが少ないことが分かる:
tkf% make all
gcc -lm -O2 rnn_ca1d.c -o rnn_ca1d-gcc-O2
gcc -lm -O3 rnn_ca1d.c -o rnn_ca1d-gcc-O3
icc -vec-report1 -O2 rnn_ca1d.c -o rnn_ca1d-icc-O2
rnn_ca1d.c(65): (col. 3) remark: ループがベクトル化されました。.
rnn_ca1d.c(67): (col. 5) remark: ループがベクトル化されました。.
rnn_ca1d.c(67): (col. 5) remark: ループがベクトル化されました。.
rnn_ca1d.c(20): (col. 3) remark: ループがベクトル化されました。.
rnn_ca1d.c(39): (col. 5) remark: ループがベクトル化されました。.
icc -vec-report1 -O3 rnn_ca1d.c -o rnn_ca1d-icc-O3
rnn_ca1d.c(65): (col. 3) remark: ループがベクトル化されました。.
rnn_ca1d.c(67): (col. 5) remark: ループがベクトル化されました。.
rnn_ca1d.c(67): (col. 5) remark: ループがベクトル化されました。.
rnn_ca1d.c(20): (col. 3) remark: ループがベクトル化されました。.
rnn_ca1d.c(39): (col. 5) remark: ループがベクトル化されました。.
gcc -lm -O2 rnn_ca2d.c -o rnn_ca2d-gcc-O2
gcc -lm -O3 rnn_ca2d.c -o rnn_ca2d-gcc-O3
icc -vec-report1 -O2 rnn_ca2d.c -o rnn_ca2d-icc-O2
rnn_ca2d.c(69): (col. 3) remark: ループがベクトル化されました。.
rnn_ca2d.c(71): (col. 5) remark: ループがベクトル化されました。.
rnn_ca2d.c(20): (col. 3) remark: ループがベクトル化されました。.
rnn_ca2d.c(39): (col. 5) remark: ループがベクトル化されました。.
icc -vec-report1 -O3 rnn_ca2d.c -o rnn_ca2d-icc-O3
rnn_ca2d.c(20): (col. 3) remark: ループがベクトル化されました。.
rnn_ca2d.c(39): (col. 5) remark: ループがベクトル化されました。.

iccの-O3オプションって-O2プラスアルファだと思ってたけど違うんだろうか． まあ，一次元配列で実装してれば問題ないことが分かったので良しとしよう!
あと，どの配列実装が速いかはたぶん計算に依存してるだろうから，この結果は 他の数値計算に適用できないと思う．これを書いてるときにちょっと添え字の 書き方間違えてて，その時の結果はかなり違ったし．だから，こういうテスト はなるべく実際の計算に近いソースで試すべきなんだろうな． 普通に行列xベクトルの計算プログラムじゃなくてRNNで試して良かった．