D言語で多次元配列ライブラリを書いた #dlang

はじめに

D言語にはmirという多次元配列のパッケージを提供する数値計算用のライブラリがある。が、実は使ったことない。
D言語で数値計算 mir-algorithmを読むとなんだか難しそうで、使う前から諦めてしまった。mapとかreduceとか便利だが、バグが入り込みやすいところで多用したくない。
なんだかんだfor文（foreach文）で書くのは分かりやすい。

しかし多次元配列のライブラリは欲しい。
自分はいままでFortranで流体用の計算コードを書いたりしているのだけれど、D言語で流体計算やりたいなと思っていて、その為にも多次元配列を簡単に使うためのライブラリが欲しい。

そんなわけでggedというライブラリを書いたのでその紹介である。

多次元配列

もともとD言語は多次元配列を扱う機能がある、気がする。たとえば

auto array = new double[][](3,5,10); /// 3x5x10の配列

こう配列をつくることができる。ただこれを扱うには

foreach(i;0..3){
    foreach(j;0..5){
        foreach(k;0..10){
            // array[i][j][k]に関する処理
        }
    }
}

と、まあ面倒くさい。何が嫌って

foreachが入れ子になる。
3,5,10と要素数に注意しなくてはならない。

何故かというと、D言語の多次元配列は正確にはジャグ配列(jagged array)と呼ばれる配列の配列でしかないからだ。
今は欲しいのは「jaggedじゃない配列」である。

ggedライブラリ

基本的な使い方

ggedで多次元配列はこのように宣言する。

Gged!(double,3) ggarray = gged!double(3,5,10); /// 3x5x10の配列

ここで型に出てくる3は,配列の（3,5,10）の長さ、つまり次元数だ。

gged配列の要素を扱うには次のようにする。

foreach(i,j,k;ggarray)
{
    /// ggarray[i,j,k]に関する処理
}

すっきり！

arrayが要素数の情報を持っているので、要素数なんて気にする必要が無いのだ。
ちなみに処理内で要素数が欲しい場合はこのように書ける。

foreach(i,j,k;ggarray)
{
    ggarray[i,j,k] = 1.*i/i.max + 1.*j/j.max + 1.*k/k.max;
}

途中で処理を挟むとき

ところで、こんな風に書きたいときもある。

foreach(i;0..3){
    foreach(j;0..5){
        f();/// kのループが始まる前にやっておきたい処理
        foreach(k;0..10){
            // array[i][j][k]に関する処理
        }
        g();/// kのループが終わった後にやっておきたい処理
    }
}

これに対してはこう書ける。

foreach(i,j,k;ggarray.Serial)
{
    if(j.once)
    {
        /// jのループ時に一度だけ行う処理
    }
    /// ggarray[i,j,k]に関する処理
    if(k.last)
    {
        /// kのループが終わるときに行う処理
    }
}

ここでSerialと書いてるのは、ggedのループをforeachを入れ子構造で書いたとき同じ順序で実行するという宣言である。
gged内では一次元配列でデータを保持しているので、Serialを使わない時はstd.parallelismをつかって並列化している。

さらに簡便に書く

次元数を気にしたくない時もあるかもしれない。

auto A = gged!double(1,2,3);
auto B = gged!double(1,2,3,4);

foreach(i,j,k;A){
    /// A[i,j,k] に関する処理
}
foreach(i,j,k,l;B){
    /// B[i,j,k,l] に関する処理
}

配列の次元数に合わせて、i,j,k...と書かなきゃいけないのも面倒くさい時には、こう書ける。

foreach(ijk; A)
{
    /// A[ijk] に関する処理
}
foreach(ijk; B)
{
    /// B[ijk] に関する処理
}

このときijkはi,j,kの情報をまとめた構造体となっており、配列のように扱うことができる。

foreach(ijk; A)
{
    assert(A[ijk] == A[ijk[0],ijk[1],ijk[2]]);
}

まとめ

ジャグじゃない配列を扱うggedライブラリの紹介をした。
今後、流体計算コード書いていくときに欲しいと思った機能を盛り込んでいく予定だ。
実行速度のことは……とりあえず二の次で。