No Such Blog or Diary

メモリバンド幅って完全に使い切ることができるのかなぁとか思ったのでプログラム書いて試してみた．でもバンド幅を使い切るのはなかなか難しい．

並列プログラムを共有メモリ環境で作っているとメモリバンド幅が足りなくて並列化の効果が頭打ちになってるっぽいことがある．特に，処理するデータは沢山あるけれど個々の要素に適用する操作が単純だという場合に顕著に見られる．でも，並列プログラムの台数効果の頭打ちがメモリバンド幅が狭いせいだと主張するにはメモリバンド幅をちゃんと把握しておかねばならない．別にカタログスペック以上の幅があるわけではないのだから無駄だけど．

実験に使ったPCは次の２台．

PC1:
 Xeon X5550 (8MB L3 キャッシュ、2.66GHz、トリプルチャネル 32GB/s MC、6.4GT/s QPI) x2
 X5520 チップセット
 12GB (2GBx6) DDR3 RDIMM メモリ(1333MHz、ECC)
 Maximum bandwidth: 64GB/s

PC2:
 Xeon E5430 (2x6MB L2 キャッシュ、2.66GHz, 1333MHz FSB)x2
 5400 チップセット
 8GB (1GBx8) クワッドチャネル DDR2-SDRAM (667MHz、ECC)
    Maximum bandwidth: 10.6GB/s

PC1のメモリバンド幅は，単体CPUがトリプルチャネルで32GB/sが可能なところに３枚のDDR3 PC3-10600が付いてるので，単体で32GB/sになって２つで64GB/sになるはず．PC２の方は，FSBのバンド幅が 10.6GB/s で，クアッドチャネルに DDR2 PC2-5300 が８枚なので 21.2GB/sで，結局FSB側で押さえられてバンド幅が10.6GB/sになる．

で，下にあるプログラムで実験した結果：

PC1は各CPUにメモリがくっついているので，使うCPU数を倍にするとバンド幅も倍になる．一方，PC2は共通のFSBの上にメモリがあるのでCPUを増やしてもバンドは増えない．

で，今のところハードウェア仕様の最大バンドの60%位しか出せていない．どやったら良くなるのかなぁと．

以下，実験に使ったプログラム．メモリを読むのみ（runreadA）＆書くのみ（runwriteNT）．SSE命令で16Bごとの読み書き．また，書き込み時にはmovntpdでキャッシュを汚染しないようにした．さらに，スレッドごとにaffinityを指定したかったので，OpenMPでなくpthreadを使った（OpenMPでもできるのかもしれないけど）．また，各スレッドでメモリを割り当てているので，スレッドの乗っかっているCPUにつながっている物理メモリにデータが置かれると期待．

#ifndef N
#define N "536870912"        // バッファサイズ
#endif
 
#ifndef I
#define I "16"            // 繰り返し回数
#endif
 
#include <iostream>
#include <vector>
#include <string>
#include <pthread.h>
#include <stdio.h>
#include <sched.h>
#include <stdlib.h>
#include <time.h>
#include <sys/time.h>
double gettimeofday_sec()
{
  struct timeval tv;
  gettimeofday(&tv, NULL);
  return tv.tv_sec + (double)tv.tv_usec*1e-6;
}
 
 
extern "C" void runreadA(unsigned char *p);
  asm volatile("#here-readA");
  asm volatile("\n"
"runreadA:                   \n"
"    pushq    %rbx                 \n"
"    pushq    %rcx                 \n"
"    pushq    %rdx                 \n"
"                               \n"
"    leaq    "N"(%rdi), %rdx      \n"
"    xorq    %rbx, %rbx           \n"
".L10runrA:                  \n"
"    movq    %rdi, %rcx           \n"
".L9runrA:                   \n"
"    movapd    (%rcx), %xmm0    \n"
"    movapd  16(%rcx), %xmm1    \n"
"    movapd  32(%rcx), %xmm2    \n"
"    movapd  48(%rcx), %xmm3    \n"
"    movapd  64(%rcx), %xmm4    \n"
"    movapd  80(%rcx), %xmm5    \n"
"    movapd  96(%rcx), %xmm6    \n"
"    movapd 112(%rcx), %xmm7    \n"
"    addq    $128, %rcx           \n"
"    cmpq    %rdx, %rcx           \n"
"    jne    .L9runrA               \n"
"    addq    $1, %rbx             \n"
"    cmpq    $"      I ", %rbx    \n"
"    jne    .L10runrA              \n"
"                            \n"
"    popq    %rdx                 \n"
"    popq    %rcx                 \n"
"    popq    %rbx                 \n"
"    ret                        \n"
"");
 
extern "C" void runwriteNT(unsigned char *p);
  asm volatile("#here-writeNT");
  asm volatile("\n"
"runwriteNT:                 \n"
"    pushq    %rbx                 \n"
"    pushq    %rcx                 \n"
"    pushq    %rdx                 \n"
"                               \n"
"    leaq    "N"(%rdi), %rdx      \n"
"    xorq    %rbx, %rbx           \n"
".L10runwNT:                 \n"
"    movq    %rdi, %rcx           \n"
".L9runwNT:                  \n"
"    movntpd %xmm0,    (%rcx)   \n"
"    movntpd %xmm1,  16(%rcx)   \n"
"    movntpd %xmm2,  32(%rcx)   \n"
"    movntpd %xmm3,  48(%rcx)   \n"
"    movntpd %xmm4,  64(%rcx)   \n"
"    movntpd %xmm5,  80(%rcx)   \n"
"    movntpd %xmm6,  96(%rcx)   \n"
"    movntpd %xmm7, 112(%rcx)   \n"
"    addq $128, %rcx            \n"
"    cmpq    %rdx, %rcx           \n"
"    jne    .L9runwNT              \n"
"    addq    $1, %rbx             \n"
"    cmpq    $"      I ", %rbx    \n"
"    jne    .L10runwNT             \n"
"                            \n"
"    popq    %rdx                 \n"
"    popq    %rcx                 \n"
"    popq    %rbx                 \n"
"    ret                        \n"
"");
 
struct info_t {
  int rank;
  int size;
  int n;
  int i;
  int cpu;
  double gb;
  void (*f)(unsigned char*);
  std::string str;
  pthread_barrier_t *st_barrier;
  pthread_barrier_t *et_barrier;
};
 
void *task(void *arg)
{
  info_t *info = ((info_t*)arg);
  int n = info->n;
  int i = info->i;
  int rank = info->rank;
  int cpu = info->cpu;
  unsigned char *p = (unsigned char*)malloc(n);
  //printf("%16lx\n", p); 
  for(int j = 0; j < n; j++) {
    p[j] = 0;
  }
 
  std::cout << "rank " << rank << " at cpu " << cpu << std::endl;
</div>
 
<div class="section">
  double st = 0;
 
  pthread_barrier_wait(info->st_barrier);
  if(rank==0) { st = gettimeofday_sec(); }
  info->f(p);
 
  pthread_barrier_wait(info->et_barrier);
 
  if(rank==0) {
    double et = gettimeofday_sec();
    std::cout << info->size << " " << (info->str) << " " << (et-st) << " " << (info->gb/(et-st))  << "GB/s" << std::endl;
  }
  free(p);
  return NULL;
}
 
int main(int argc, char *argv[])
{
  if(argc<=1) {
    std::cout << argv[0] << " p" << std::endl;
    return 0;
  }
  int n = atoi(N);
  int i = atoi(I);
  int p = atoi(argv[1]);
  double gb = (1./1024.0/1024.0/1024.0*i*n*p);
  std::cout << "# run (buffer) size: "<< (1./1024.0/1024.0*n) <<"MB." << std::endl;
  std::cout << "# iteration: "<< i <<" times" << std::endl;
  std::cout << "# "<< (1./1024.0/1024.0/1024.0*i*n)<<"GB/core will be transfered." << std::endl;
  std::cout << "# in total "<< gb <<"GB will be transfered." << std::endl;
 
  pthread_barrier_t st_barrier, et_barrier;
 
  pthread_t ps[p];
  info_t is[p];
  cpu_set_t cs[p];
  for(int k = 0; k < p; k++) {
    is[k].rank = k;
    is[k].size = k;
    is[k].n = n;
    is[k].i = i;
    is[k].gb = gb;
    is[k].st_barrier = &st_barrier;
    is[k].et_barrier = &et_barrier;
    is[k].cpu = argc > 2 ? k : ((k&1)<<2) | ((k&0x6)>>1);
 
    CPU_ZERO(&cs[k]);
    CPU_SET(is[k].cpu, &cs[k]);
  } 
  std::vector<std::string> ss;
  std::vector<void (*) (unsigned char*)> fs;
  fs.push_back(runreadA);  ss.push_back("ra");
  fs.push_back(runwriteNT); ss.push_back("wnt");
 
  for(int l = 0; l < fs.size(); l++) {
    pthread_barrier_init(&st_barrier, NULL, p);
    pthread_barrier_init(&et_barrier, NULL, p);
    for(int k = 0; k < p; k++) {
      is[k].f = fs[l];
      is[k].str = ss[l];
    }  
</div>
 
<div class="section">
    for(int k = 0; k < p; k++) {
      pthread_create(&ps[k], NULL, task, &is[k]);
      pthread_setaffinity_np(ps[k], sizeof(cpu_set_t), &cs[k]);
    }
    for(int j = 0; j < p; j++){
      int *tmp;
      pthread_join(ps[j], (void**)&tmp);
    }
    pthread_barrier_destroy(&st_barrier);
    pthread_barrier_destroy(&et_barrier);
  }
}

valgrind --tool=cachegrind でキャシュミスを調べようと思ってプログラムを走らせたのだが…　計算が終了する気配が無い．データサイズ小さくすりゃ速くなるだろうけど，そしたらデータがキャッシュ乗り切ってしまうので本末転倒．

まあ，キャッシュ効きまくり状態でも50秒かかるのだからシミュレーションでキャッシュの計測とかしてたらなかなか終わらんか．うーん，この実験を繰り返すのは苦痛だなぁ．

ということで，せっかくGoを勉強したのでゴルフ場で（頭と指の）運動をしてみた．

大人気なく問題を解きまくってGoで解いた人ランキング一位になった．

以下，ショートコーディングのための技術．

• Raw String Literal を使うべし． 特に，改行にエスケープがいらないため，

  s:="hoge\nhuga"

  が

  s:=`hoge
  huga`

  でいいので1Bの短縮．他にもバイナリデータを文字列として埋め込んでしまうという…

• 外部パッケージをローカルの名前空間に展開すべし．
  インポート時の束縛名にドットを使うと，他パッケージの関数などをローカルの名前空間に展開出来る．なので，

  import"fmt"func m(){fmt.Print(“c")}

  が

  import."fmt"func m(){Print(“c")}

  に短縮できる．複数のパッケージをローカルに展開できるが，同一名があると衝突するので，解決に短い名前を付ける必要がある．ついでにインポートのファイル名の後には空白がいらないという…

• 関数変数を使うべし．
  関数も変数に入れられる．なので，

  x:=Sprintf("%c",c);y:=Sprintf("%c",d);

  が

  f=Sprintf;x:=f("%c",c);y:=f("%c",d);

  に短縮できる．でもprintlnとかの特別な関数は束縛出来ないという…

• 依存の無い（初期値）代入はまとめるべし．
  変数への値代入は，型に関係なくカンマでつないで一気に書ける．なので，

  x:=1;ｙ:="str";z:=2

  が

  x,y,z:=1,"str",2

  に短縮できる．でも，複数の値を返す式が右辺にあると怒られるという…　（複数値のコンテキストのネストが出来ないらしい） ついでに，複数値の代入文は右辺が全て評価されてから左辺へと代入されるので，

  t:=x;x=y;y=t

  変数のスワップなどは複数値の代入を使って

  x,y=y,x

  のように短く書ける．

• Sliceは上手に作るべし．
  Sliceの生成にはmakeを使う方法とcomposite literalを使う方法がある．なので，小さなsliceをつくるなら

  b:=make([]byte,1)

  でなく

  b:=[]byte{1}

  にする．また，大きなsliceも

  b:=make([]byte,10000)

  でなく

  b:=make([]byte,1e4)

  にする．

並列に動くっぽいけど本当にそうなのか気になったので，色々書けるようにと勉強してみた．

セミコロン少なめだし，カッコ少なめだし，型を書くのも少なめだし，タプルあるし，型システム単純だし，例外ないし，オーバーロードないし，generics 入ってないし，GCだし，メソッド定義にはレシーバ書くし，継承ないし，interface とかメソッドさえ揃っていれば継承したことになるし，embedded types とか楽すぎだし，宣言したのに使ってない変数があるとコンパイルしてくれないし，スコープ上書きしたときに何も言ってくれないし，速度も十分速いし．スクリプト言語の気分で速いプログラムが書けるような気がする．一方で，generic types がないので高レベルな汎用ルーチンが書きづらい．

で，そんなことはさておきとりあえずベンチマークということで，n-queen問題のプログラムを書いてみた：

package main
import ("fmt"; "flag")
 
func nqueen(n, i, u, l, r int) int {
  if i==n { return 1 } 
  nl := (l << 1) | 1;
  nr := ((r|(1<<uint(n))) >> 1);
  ni := i + 1;
  p := u & nl & nr;
  c := 0;
  for p!=0 {
    lb := (-p) & p;
    p ^= lb;
    c += nqueen(n, ni, u ^ lb, nl ^ lb, nr ^ lb)
  }
  return c
}
 
type State struct { n, i, u, l, r int }
func nqueenRes(ch chan int, n, i, u, l, r int) { 
  ch <- nqueen(n, i, u, l, r)
 }
 
func gentasks(cs *[] chan int, lim, n, i, u, l, r int) {
  if i==lim {
    lx := len(*cs);
    *cs = (*cs)[0:lx+1];
    (*cs)[lx] = make (chan int);
    go nqueenRes((*cs)[lx], n, i, u, l, r);    
    return;
  }
  nl := (l << 1) | 1;
  nr := ((r|(1<<uint(n))) >> 1);
  ni := i + 1;
  p := u & nl & nr;
  for p!=0 {
    lb := (-p) & p;
    p ^= lb;
    gentasks(cs, lim, n, ni, u ^ lb, nl ^ lb, nr ^ lb)
  }
}
 
func nqueensGo(n int) int {
  b := (1<<uint(n)) - 1;
  cs := make( [] chan int, 0, n * n * n);
  gentasks(&cs, 3, n, 0, b, b, b);
 
  sum := 0;
  for _, c := range cs { sum += <- c }
 
  return sum;
}
 
var n = flag.Int("n", 16, "n")
 
func main() {
  flag.Parse();
  fmt.Println(nqueensGo(*n));
}

ついでに比較対象のC++プログラム（OpenMP使用）：

#include<iostream>
#include<vector>
#include<omp.h>
#include<stdlib.h>
 
using namespace std;
 
int nqueen(int n, int i, int u, int l, int r) 
{
  if(i==n) return 1; 
  const int nl = (l << 1) | 1;
  const int nr = ((r|(1<<n)) >> 1);
  const int ni = i + 1;
  int p = u & nl & nr;
  int c = 0;
  while(p!=0) {
    const int lb = (-p)&p;
    p ^= lb;
    c += nqueen(n, ni, u ^ lb, nl ^ lb, nr ^ lb);
  }
  return c;
}
struct state {
  int n, i, u, l, r;
  state(int n, int i, int u, int l, int r) : n(n), i(i), u(u), l(l), r(r) {}
};
int nqueen_res(const state &st) { return nqueen(st.n, st.i, st.u, st.l, st.r); }
 
void gentasks(vector<state> &xs, int lim, int n, int i, int u, int l, int r) 
{
  if(i==lim) {
    xs.push_back(state(n, i, u, l, r));
    return;
  }
  const int nl = (l << 1) | 1;
  const int nr = ((r|(1<<n)) >> 1);
  const int ni = i + 1;
  int p = u & nl & nr;
  int c = 0;
  while(p!=0) {
    const int lb = (-p)&p;
    p ^= lb;
    gentasks(xs, lim, n, ni, u ^ lb, nl ^ lb, nr ^ lb);
  }
}
 
int nqueensOMP(int n)
{
  int b = (1<<n) - 1;
  vector<state> xs;
  gentasks(xs, 3, n, 0, b, b, b);
 
  int nx = xs.size();
  int sum = 0;
#pragma omp parallel for reduction(+:sum)
  for(int i = 0; i < nx; i++) {
    sum = sum + nqueen_res(xs[i]);
  }
  return sum;
 
}
 
int main(int argc, char *argv[])
{
  int n = argc > 1 ? atoi(argv[1]) : 16;
  int c = nqueensOMP(n);
  std::cout << n << " " << c << std::endl;
  return 0;
}

とりあえず実験結果．OpenMPのスケジューリングを変更して性能を改善したものの結果は下のほう．

8コアのマシン(Xeon E5430 @ 2.66GHz x2)での実験結果（n=1７）．C++ はg++4.5.0 (exeperimental) に -O3，Go は 6g/6l でコンパイル．OMP_NUM_THREADS と GOMAXPROCS でスレッド数を指定．speedup は C++ の１スレッドに対する速度比．

どうやらC++の1.3倍程度の遅さらしい．問題が問題だけに台数効果はほぼリニアに出ている（スレッド数増やした時のC++の伸び具合が今ひとつだが…　スケジューリング変えるべきか？）．Goは性能的に何の問題もなく優秀だねぇ．まあ，この手の極単純なループの並列化に関してはOpenMPの方が楽ではあるけれど．

スレッド数増やしたときにGoの方が速くなっているのはスケジューリングがうまいからだろうなぁ．OpenMPもラウンドロビンとかにすれば改善するかも？

ということで，サイズ１でラウンドロビンするように

#pragma omp parallel for reduction(+:sum) schedule(dynamic)

に変えてC++の実験を取り直した．

8コアマシンで．

24コアマシンで．

うん，常にC++の方が速い結果になったので問題なし．スケジューリングのオーバヘッドが無視出来る問題だから，スケジューリングの柔軟さがやっぱ効いてくるねと．

あー，New wording for C++0x Lambdas (rev. 2) (N2927) にλ式はdecltypeの引数にはなれないと書いてあった．残念．

In addition, this rewrite adds the restriction that lambda expressions cannot be used in the operand of a sizeof operator, alignof operator, or decltype specifier.That restriction—suggested by Doug Gregor and John Spicer—avoids severe implementation difficulties with template argument deduction.

ということは，λ式を返す関数は戻り値の型にautoと書けないのか．λ式の正確な型の書き方を知る必要があるな．できるかどうかすら分からないけど．

うーむ，decltypeに直接λ式を入れるだけでなく，部分式にすらλ引きが許されないのか…．しかもN2927とかN2550とか読む限りではλ式が変換されるクロージャオブジェクトにはユニークだけど名無しの型がつくそうで…．やっぱλ式を返すにはstd::functionに突っ込むしかないのかな（何の問題もないのだけど何か負けた気分）．

下のようにdecltypeの中にλ式を書いたらコンパイル通らなかった．うーん，expressionなら何でもよいと思ったのだが何か間違っただろうか？

#include <iostream>
 
template <typename T>
auto func(T i) -> decltype(i) { return i; }
 
template <typename T>
auto func1(T i) -> decltype([=](int j){ return i+j; })
 { return [=](int j){ return i+j; }; }
 
int main(int argc, char *argv[])
{
        func(1);
        func1(1);
        return 0;
}

上側のfuncの定義は通る．でも下のfunc1の定義が通らない．expected primary-expression before ')' token （λ式の後の)の手前に式が必要）とか言われる．

そしてλ式とテンプレートの同時使用の仕方が分からない．ポリモーフィックなλ式は作れるのか？