ナンクル力学系

学んだ事を書き連ねていこう。

Posts Tagged ‘C++

RailGun: C+Pythonでお手軽に数値計算プログラミングを加速させるライブラリ

leave a comment »

シミュレーションのプログラム書く時に,毎回同じようなコード書くのが嫌なので作ってみました.とりあえずドキュメント書いてみたから使ってみてよ!

Written by tkf

October 15, 2010 at 4:08 am

Posted in 研究日誌

Tagged with , ,

Cの配列はa[i][j]とa[i*Nj+j]のどちらが速いか試してみた

with 3 comments

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

と思ってたけど,

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

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

サンプルで作ったのはの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で試して良かった.

Written by tkf

January 1, 2010 at 11:00 pm

Posted in PC

Tagged with

Rで離散時間の力学系を計算するためにC拡張を書いた

leave a comment »

x_0 を与えたときに、 x_{t+1} = f(x_t) を使って 時系列 (x_0, x_1, x_2, \dots ) を計算するって単純なことが したいだけだったのに、上手くRでやる方法が思いつかなかった。

ので、Cで拡張作った。

とりあえず動いたver

Cのコード
 #include <R.h>
 #include <Rinternals.h>

SEXP dtds(SEXP _time, SEXP x0, SEXP fn, SEXP rho)
{
  int i, dim = length(x0);
  int t, time = INTEGER(_time)[0];
  SEXP xt, x1, x2;

  PROTECT(xt = allocMatrix(REALSXP, time, dim));
  for(i = 0; i < dim; i++){
    REAL(xt)&#91;time*i&#93; = REAL(x0)&#91;i&#93;;
  }
  for(t = 1; t < time; t++) {
    PROTECT( x1 = allocVector(REALSXP, dim) );
    for (i = 0; i < dim; ++i){
      REAL(x1)&#91;i&#93; = REAL(xt)&#91;t-1 + time*i&#93;;
    }
    defineVar(install("x"), x1, rho);
    UNPROTECT(1);

    PROTECT( x2 = coerceVector(eval(fn, rho), REALSXP) );
    for(i = 0; i < dim; i++){
      REAL(xt)&#91;t + time*i&#93; = REAL(x2)&#91;i&#93;;
    }
    UNPROTECT(1);
  }
  UNPROTECT(1);
  return(xt);
}
&#91;/sourcecode&#93;</div>
<div class="outline-5">
<h5>Rの関数(ラッパ?):</h5>

dyn.load("dtds.so")

dtds <- function(N, x0, fmap){
  .Call("dtds", as.integer(N), as.double(x0), quote(fmap(x)), new.env())
}

fmap <- function(x){
  x * (x - 1)
}

rms <- function(x,y){
  d <- c(x)-c(y)
  sqrt( d %*% d/length(d) )
}
&#91;/sourcecode&#93;</div>
<div class="outline-5">
<h5>実行例:</h5>

> N  <- 100000
> x0 <- c(-0.1,0.1)
>
> xt <- matrix(rep(x0, N),N,byrow=T)
> system.time(for (t in 2:N){ xt[t,] <- fmap(xt&#91;t-1,&#93;) })
   ユーザ   システム       経過
     1.936      0.040      1.977
> print(xt[N,])
[1]  0.002237959 -0.002232851
>
> system.time(xs <- dtds(N, x0, fmap))
   ユーザ   システム       経過
      0.42       0.00       0.42
 # for使って計算したやつの4倍くらい?
> print(xs[N,])
[1]  0.002237959 -0.002232851
>
> print(rms(xt,xs))
     [,1]
[1,]    0
 # for使って計算したやつと同じ結果!
あれ?

4倍しかはやくならない?

くやしかったので少し改良したver

  • 毎ステップ allocVector しなくて良いはず
  • defineVar もボトルネックか?
Cのコード
 #include <R.h>
 #include <Rinternals.h>

SEXP dtds(SEXP _time, SEXP x0, SEXP fn, SEXP rho)
{
  int i, dim = length(x0);
  int t, time = INTEGER(_time)[0];
  SEXP xt, x1, x2;

  PROTECT(xt = allocMatrix(REALSXP, time, dim));
  PROTECT( x1 = allocVector(REALSXP, dim) );
  defineVar(install("x"), x1, rho);
  for(i = 0; i < dim; i++){
    REAL(xt)&#91;time*i&#93; = REAL(x0)&#91;i&#93;;
  }
  for(t = 1; t < time; t++) {
    for (i = 0; i < dim; ++i){
      REAL(x1)&#91;i&#93; = REAL(xt)&#91;t-1 + time*i&#93;;
    }

    PROTECT( x2 = coerceVector(eval(fn, rho), REALSXP) );
    for(i = 0; i < dim; i++){
      REAL(xt)&#91;t + time*i&#93; = REAL(x2)&#91;i&#93;;
    }
    UNPROTECT(1);
  }
  UNPROTECT(2);
  return(xt);
}
&#91;/sourcecode&#93;</div>
<div class="outline-5">
<h5>実行例:</h5>

> N  <- 100000
> x0 <- c(-0.1,0.1)
> #x0 <- c(0.1)
>
> xt <- matrix(rep(x0, N),N,byrow=T)
> system.time(for (t in 2:N){ xt[t,] <- fmap(xt&#91;t-1,&#93;) })
   ユーザ   システム       経過
     1.944      0.040      1.981
> #print(length(xt[1,]))
> #print(length(xt[,1]))
> print(xt[N,])
[1]  0.002237959 -0.002232851
>
> system.time(xs <- dtds(N, x0, fmap))
   ユーザ   システム       経過
     0.388      0.004      0.392
> #print(length(xs[1,]))
> #print(length(xs[,1]))
> print(xs[N,])
[1]  0.002237959 -0.002232851
>
> print(rms(xt,xs))
     [,1]
[1,]    0
うーん
  • ほとんど変化なし
  • やはり、evalがだめな感じがする

感想

PythonのC拡張比べてみて。(昔書いたPythonのC拡張はこんな感じ。 この例と同じ離散時間の系の計算だけど、fmapにあたる写像は決め打ち。)

  • PythonでC拡張書くより、全然書く量が少ないのは良い
  • ドキュメントがなさすぎるw
    • Writing R Extensions だけ
    • リファレンスが無いので、関数の役目が謎なままつかうことに
  • 関数とかマクロとかフリーダム
    • PythonのCコードは必ず「Py」から始まっていたり、マクロならそれ以下が大文字 になっていたりと安心する
    • まー、Cだし、コンパイラが怒るので問題なし?
  • 何回でも dyn.load(“dtds.so”) と dyn.unload(“dtds.so”) できるのはいいね!
    • Pythonだと、一度読み込んだCのモジュールは reload しても変化なしなので,Python自体を再起動する必要があった
    • デバッグ楽
  • Cコードレベルのミス(?)に気づくRいいね
    • 何度、 「 エラー: VECTOR_ELT() can only be applied to a ‘list’, not a ‘double’」 と怒られたことか。。。
  • 慣れたらサクサクかけそう!
    • 慣れる = Rのソース読めるくらいか?w

Written by tkf

March 11, 2009 at 7:17 pm

Posted in PC

Tagged with , ,

PythonのC拡張でAPIを公開する例

leave a comment »

PythonのC拡張を普通に書くと、その関数などは他のC拡張モジュールから使えない。という訳で他のC拡張から読み込めるようにAPIを公開する必要がある訳だけど、その作り方が少しトリッキーだったので公開してみる。

参考:拡張モジュールに C API を提供すると、NumPyのソース。

(追記。そういえばPython/C拡張の練習で書いたコードがあったので晒す:http://snipplr.com/view/8215/example-of-numpyc-api/ C APIは使ってないけど。)

あんまりPythonのC拡張どころか、C言語もそんなに触って無いのでおかしい所があるかも。ツッコミ大歓迎です!

準備するのは以下のファイル。son.Sonとfather.Fatherというモジュールを作ってFatherのメンバにSonをもってきてみた。

Read the rest of this entry »

Written by tkf

September 18, 2008 at 10:42 pm

Posted in programming

Tagged with , ,

STLのシャッフルアルゴリズムstd::random_shuffleの乱数にBoost/random.hppを使う

with 2 comments

C++の標準ライブラリSTLに入ってるシャッフルを使いたいんだけど,ちゃんとした乱数を使いたいってことでboostを使ってみることにした.

コードはいつものようにsnipplrに:

http://snipplr.com/view/5907/stdrandomshuffle-and-boostrandomhpp/

C++はちゃんと勉強せずにプログラムを書いたり読んだりしてるだけなので,変な所があるかもしれない.

参考(というか↓を組み合わせたら↑になるんだけどねw):

Written by tkf

April 19, 2008 at 1:34 pm

Posted in programming

Tagged with , ,