No Such Blog or Diary

EuroPar2007で発表のあった The Multi-Core Standard Template Library (MCSTL) が GCC に統合されている（されつつある？）らしいので，セミナーのネタにと試用してみた．GCCのオフィシャルページを見ると libstdc++ parallel mode という名前で統合中らしい．こいつは既存の STL を OpenMP による並列実装の STL に置換しましょうというもので，既存の STL で書かれたコードがそのままで並列プログラムになってくれる．

んで，試そうと思ったら最新のGCCリリースであるところの 4.2.2 にも入っていないということに GCC4.2.2 のビルドが終わった時点で気づき…　しょうがないので svn のリポジトリから最新のソースを落としてきて GCC をビルドした（８コアくらいあるマシンで make -j8 とかやるとビルドがとても速い）．とりあえずソースの状態で libstdc++-v3/include/parallel というディレクトリがあれば parallel mode が使える．

準備できたので次の std::sort を使ったプログラムを試してみた．

#include <iostream>
#include <algorithm>
#include <vector>
const int count = 1000000, num_tests = 10;
int main(int argc, char **argv)
{
    std::vector<int> v(count);
    for (int i = 0; i < num_tests; i++)
    {
        std::generate(v.begin(), v.end(), rand);
        std::sort(v.begin(), v.end());
        std::cout << "." << std::flush;
    }
}

コンパイルは parallel mode 用のフラグ（-D_GLIBCXX_PARALLEL -fopenmp）を立てればよいらしい．

 g++-XXX test.cpp -o testS -march=native  # sequential mode
 g++-XXX test.cpp -o testP -march=native -D_GLIBCXX_PARALLEL -fopenmp  # parallel mode

Quad の Xeon x2 という８コアで動かした結果は以下のとおり．

> OMP_NUM_THREADS=8 time ./testS   # sequential mode
..........        6.26 real         6.24 user         0.01 sys
> OMP_NUM_THREADS=8 time ./testP   # parallel mode
..........        1.39 real         5.66 user         0.75 sys

普通に４～５倍速くなってるのを確認．すげー．並列プログラミング簡単．とりあえず sort は簡単に並列に動いてくれた．

ついでに， parallel mode は他にも並列実装を持っているらしいので試してみた．

まず std::partial_sum．4.51秒が1.65秒になった．８コアで約四倍．

#include <iostream>
#include <algorithm>
#include <vector>
#include <numeric>
struct one {
  int operator()() const {
    return 1;
  }
};
const int count = 10000000, num_tests = 10;
int main(int argc, char **argv)
{
    std::vector<int> v(count);
    std::vector<int> r(count);
    std::generate(v.begin(), v.end(), one());
    for (int i = 0; i < num_tests; i++)
    {
        std::partial_sum(v.begin(), v.end(), r.begin());
        std::cout << "." << std::flush;
    }
}

次は std::accumulate．3.20秒が3.79秒に増えた．全要素の和をとるだけなので普通に考えると８コアで８倍いきそうだけど…　速くならない理由は不明．

#include <iostream>
#include <algorithm>
#include <vector>
#include <numeric>
struct one {
  int operator()() const {
    return 1;
  }
};
const int count = 10000000, num_tests = 10;
int main(int argc, char **argv)
{
    std::vector<int> v(count);
    std::generate(v.begin(), v.end(), one());
    for (int i = 0; i < num_tests; i++)
    {
        int sum = std::accumulate(v.begin(), v.end(), 0);
        std::cout << "." << std::flush;
    }
}

最後に試したのは std::for_each．3.81秒が3.75秒．全要素に独立に関数適用するだけなので８コアで８倍いきそうだけど…　速くならない理由は不明．

#include <iostream>
#include <algorithm>
#include <vector>
#include <numeric>
struct f {
  void operator()(int &x) const {
    x = x + 1;
  }
};
struct one {
  int operator()() const {
    return 1;
  }
};
const int count = 10000000, num_tests = 10;
int main(int argc, char **argv)
{
    std::vector<int> v(count);
    std::generate(v.begin(), v.end(), one());
    for (int i = 0; i < num_tests; i++)
    {
        std::for_each(v.begin(), v.end(), f());
        std::cout << "." << std::flush;
    }
}

ということで，簡単に並列化できそうな部分でいろいろと失敗した模様．使い方が悪いのか実装途中なのか？　よくわからんけどうまく動いた部分に関しては効果が確認できたのでよし．

Fortress のプログラムでインデントに TAB 使ってたら Parse エラーが出まくった．そして原因が TAB であることに気づくまでプログラムを何度も書き直す羽目に．Spec のどこかに TAB に意味があるって書いてあったっけかな？

BSP実装がなんか変な挙動を示すなぁ．AlltoAllで書きなおした方がいいのだろうか？とりあえずもう一セットやってみっか．

寝付けないので自前で考えた解き方を試しに実装してみた．解説聞いてないので合ってるかどうかは知らんけど．

A: LinkedList で普通に要素を抜いてくだけ．最大サイズでも一瞬で終わるのでこれで問題なさそう．

B: エラトステネスのふるいで素数求めたら終わり．与えられた数から前後の素数を探索するのは前後を順に見てくだけでいい．

どうせすぐ近くに素数があるだろうから．

C: 状態遷移行列を T乗すると考えるかT回の状態遷移を計算すると考えるか…　状態遷移行列が疎なのでグラフ上で状態遷移をシミュレートするイメージの実装．スタートに戻るマスへの遷移はスタートへ足しこむ．一回休みは一回休み用のノードを作って一度そのノードへ遷移するようにする．計算量は O((N+L)T) なので時間も問題ないでしょう．

D: 面倒なので考えてない．一つの三角形を決めたら他も決まるから全部の点を試せばいいのかな？

E: これもグラフを作るまでが面倒．グラフさえできれば最短路なので楽．

F: ソートした辺のリスト上で最大重み用のポインタと最小重み用のポインタをスライドさせてく．二つのポインタの間にある辺集合で全域木が作れるなら最小のポインタを進める．作れないときは最大側のポインタをするめる．全域木ができるかの判定は連結しているかの判定でできるので適当に Unon/Find でやった．これで全体としては O(E^2) だけど…　手元では最大サイズが一秒以内なのでギリギリオッケー．

G: 状態空間のサイズが 2^24 程度なので普通に幅優先でいい気がする．ただ最大サイズのサンプル入力で４秒かかるのでもう少しコードを特化した方がいいかもしれない．

H: 再帰下降で Exception なげて帰ってくる．Hashtable とか使わないとインデックスで困るけどそれだけだと思う．

I: 面倒なのでパス．多角形を周りから削っていって最後に残ったところか？何となく島が海に沈んでくイメージ…

J: 連立方程式立てて解けば終わりなのだけど…　有理数で一次方程式を解かねばならないし（doubleとかでもいいのかな？），そもそも方程式の係数を求めるのも面倒．コード書きたくないから二乗根の場合に関して手で解いてみただけ．久々に 4x4 の一次式を解いた．

とりあえず，ABCFGH の実装に（テレビ見をながらではあるけど）５時間ほど．やっぱ９問解いたチームってのはすごいなぁ．

とりあえずやってみることその一．演算子のオーバロードがダメと言われても試したくなるのが人間．

struct DOUBLE {
    double x;
    double y;
    DOUBLE() : x(0), y(0) { }
    DOUBLE(int d) {
        x = d;
        y = d;
    }
    DOUBLE operator+(const DOUBLE& a) const {
        DOUBLE z;
        z.x = x + a.x;
        z.y = y + a.y;
        return z;
    }
};
DOUBLE a1(int n, DOUBLE *a)
{
    int i;
    DOUBLE r = 0;
#pragma omp parallel for reduction(+:r)
    for(i = 0; i < n; i++) {
        r = r + a[i];
    }
    return r;
}

コンパイルしたら

 error: 'r' has invalid type for 'reduction' 

と文句を言われた．だめらしい．

次，中途半端にオーバーロードしてみる．

struct DOUBLE {
        double x;
        double y;
        DOUBLE() : x(0), y(0) { }
        DOUBLE(int d) {
                x = d;
                y = d;
        }
};
 
double operator+(double a, DOUBLE b) 
{
    return b.x - a;
}
 
double a1(int n, DOUBLE *a)
{
    int i;
    double r = 0;
#pragma omp parallel for reduction(+:r)
    for(i = 0; i < n; i++) {
        r = r + a[i];
    }
    return r;
}

コンパイルしても文句言われない．でも演算子の結合性がないので当然ながら結果はおかしい．各スレッドでの結果をマージする部分では通常の + になってしまうので当り前だけど．

ついでに，次の無意味なオーバーロードは正しい計算がされる．

double operator+(double a, DOUBLE b) 
{
    return a - b.x;
}

結局，全要素に - を map して + で reducction するだけだから．

結論：reduction に演算子オーバーロードはやっぱり使えなかった．でもこうなると行列を for 文で掛けまくるとかいう操作は並列化してくれないのかぁ．行列積自体は結合的だけど展開した式での各要素は reduction の式にならないし．使えん．

体力が尽きているので一番簡単な問題だけ解いて寝なおした．

250点問題：フォントを小さい方から試せばいいだけなのですぐ終わる．

500点問題：チェックしなければならない時刻は，任意の歩行者の組の間隔が車幅になるときと歩行者と横断歩道の開始点との差が車幅になるときなので，これらをチェックするだけ．でも書くのが面倒なので寝た．

1000点問題：螺旋かどうかの判定に線分の交差判定とか考えないといけないから…　パス．まあ，判定ができたとしても総当たりでいいのかどうかよくわからん．