CEANによる配列操作 - いぐにさんにっき

はじめに

これはC++ Advent Calendar 2012の21日目の記事です。
本当はスパコン環境でのC++とかODEintについて書こうかと思っていたのですが、時間がなくなってしまったので予定よりも軽めの内容にします。もしかしたらODEintについて追記するかもしれません。しないかもしれません。

...というか書いてる途中で時間切れになったのでまた後で修正します（申し訳ない)

さてこのCEANというのは何かというと、C/C++ Extensions for Array Notationsのことです。Array Notation。つまり配列の記述方法です。残念ながらC++標準では使えません。Intel Compilerでのみ利用できる拡張機能です。Intel Compiler 13で使えることは確認してます。
拡張機能、しかも文法上の拡張なので使うと他コンパイラとの互換性が無くなります。そのため、この記事で紹介する機能は積極的な使用を推奨するものではありません。

Introduction

C++は実行速度がそこそこあるので数値計算分野にも使えるわけですが、現状として私の環境、特にシミュレーションでは数値計算はほぼFortranです。私も今年C++ほとんど書いてません。というか書けませんでした。スパコンでの利用を考えるとC++はなかなか難しい気もします。速度面での問題もあるのですが、やはりFortranは数値計算コードを書きやすい感があります。Fortranといっても昔のじゃないですよ、Fortran95(以降)です。ちょっとこの辺の例も書いてみます。(時間なくてコードチェックしてないのでミスしてたらすみません)
たとえば、３次元での計算コードとかだとほとんどの場合３重ループになるかと思います。

for(int i=0;i<NX;++i) {
  for(int j=0;j<NY;++j) {
    for(int k=0;k<NZ;++k) {
         fx[i][j][k] = (p[i+1][j][k]-p[i-1][j][k])*idx
    }
  }
}

Fortranで書くとこうなります。

forall(i=1:NX-1, j=0:NY, k=0:NZ)
   fx(i,j,k)=(p(i+1,j,k)-p(i-1,j,k))*invdx
end forall

並列化も期待できます。実際はもっと複雑な式になるのでforallを使いましたが、本当はこれだけならさらに簡単に書けます。

fx(1:NX-1, :, :) = (p(2:NX, :, :)-p(0:NX-2, :, :))*invdx

はい。コロンを使った書き方が出てきました。別の例を見てみます。めんどくさいのでC++での例は省略。

a(4:6,4:6) = b(4:6,4:6)  !  配列の一部をコピーします。
S = sum(c(c>0)) !  配列cの内、正の値を持つものの合計
d[0:NX:2]  = 1 ! 偶数番目に1を代入 
e(:) = sin(e(:))  !  e(i) = sin(e(i))をeの全要素に対して適用

こんな感じで、とにかく配列操作が強力です。こういう書き方ができるのでC++よりFortranで書いた方がかなり楽なことが多いです。数値計算のコードはあまり複雑にならないですし(多分)、やはりそれ専用の言語は強いです。例に上げた中でいくつかは、C++でもstd::valarrayとかのライブラリを使えば近いことはできます、sliceとかもちゃんとありますし。でもFortranほどシンプルには書けません。ちなみに、こういうコロンを使ったアクセス方法はPythonのnumpyもできますね。MATLABとかもできた気がします。便利です。C++でも使いたいです。FortranじゃなくてC++を使いたいんです。CEANならできます！しかもCEANはvectorizationもしてくれます。つまりSIMDで並列計算されます。言語拡張だけどちょっと触ってみてもいいんでない？

CEANのつかいかた。

最初に書きましたがIntel Compilerを使います。最近のならデフォルトで有効になってるはずです。追加で必要なコンパイルオプションはありません。
基本的には他言語のNotationと同じです。ただし、指定するのは範囲ではなくて、下限(lower bound)と長さ(length)です。

syntax定義：

 [ :  [: ] ]
             [ :  [: ] ] …

Simple example:

A[:]   // Aの全要素
B[6:2] // B[6]からB[8]まで。指定するのは「長さ」
C[:][5] // 多次元配列も可
D[0:8:2] // 偶数番目
E[:]*F[:] // 対応する要素同士で掛けます。内積とか外積ではない。
G[3:2][3:2] + H[4:2][2][6:2] // 次元やサイズがあってればよい。
I[:][:] + 3  // scalarを足したり掛けたりはOK

Build-in function:

リダクションとかは関数が用意されてます。
__sec_reduce_add(J[:] * K[:]) // 内積
ほかにもいろいろ。

後で書きなおす編

(...ここまで書いて時間なくなってきた。後でまた直そう)

残念ながらはポインタか配列型しかとれないのでstd::vectorとかは無理っぽい

//  std::vector<int> A{1,2,3,4,5,6,7,8,9}; // 無理だった
//  std::vector<int> A{9}; // 通るけど A(9)になった。えー
    std::vector<int> A(9); // サイズを指定するときは丸カッコですよね
    int idx[9];
    idx[:] = __sec_implicit_index(0); // idx[i] = i になる。便利。
                                      // 0から始まるという意味じゃなくて0番目のランクという意味です
                                      // 多次元配列だとさらに便利
//  A[1:4] += 1; // NG
//  A[idx] = idx[:]; // NG
    A[idx[:]] = idx[:]; //OK
    for(auto v: A) std::cout << v << ","; //0,1,2,3,4,5,6,7,8

Fortranの例で出した

e(:) = sin(e(:))

はBuiltin関数としてsinが使えるので

e[:] = sin(e[:])

と書けます。Builtinではない任意の関数で同様のことをするには

__declspec(vector) double sincos(double a, double b,){
   return sin(a)+cos(b);
}
//...
A[:] = sincos(B[:], C[:]);

でできます。都合により私の時間がないので、後で書きなおすことにして次のif文で締めます。

CEAN with if statement

いわゆるmaskingができます。流体解析の風上差分とかでは速度の符号によって計算式が変わったりすることもあるので、そういうときに便利かな？

    // CEAN with if statement
    int mask[8] = { 0,1,1,1,0,1,1,0 };
    if(mask[:]==0) { mask[:]=4; }                                                          
    else { mask[:]=8; }
    for(auto v: mask) cout << v << ","; // 4,8,8,8,4,8,8,4,

今回説明したようなことは次のリンク(PDF)にほとんど書いてありますので、
ぶっちゃけこのPDF見たほうが早いですね。時間ができたらまた続き書きます。ODEintは...年内は厳しいかなあ...
Intel Documentation
http://software.intel.com/sites/default/files/m/0/e/c/3/d/28991-Effective_Data_Parallelism_Using_CEAN.pdf

次は@hgodaiさんの「ゲームプログラミングの都市伝説」です。