std::variantと不動点演算子

実は、std::vectorC++17から不完全型をサポートしたので、std::variantと組み合わせると以下のようなコードが書ける。

struct config_t
{
    std::variant<bool, int, double, std::string, std::vector<config_t>> data;
};

これを上手く使えばJSONなどの木構造をかなり簡単に(ポインタを陽に触らずに)扱える。std::mapが不完全型をサポートしてくれていないのがつらいが、std::vector<std::pair<std::string, config_t>>とすることで強引に回避できる。

これをしたかったらしき人の質問を見つけた。

stackoverflow.com

これが通らないとのことだ。

class ScriptParameter;
using ScriptParameter = std::variant<bool, int, double, std::string, std::vector<ScriptParameter> >;

Best answerの人のいう通り、これは最初にclassとしてScriptParameterを設定しているから2重定義になるのであって、クラスの中に入れておけばよい。この回答のコードに関しては、継承することを想定していないクラスを継承するのはいつも少し不安になるので私はいつも所有関係にするが、まあ今回の主題ではない。

ここに異彩を放つ回答がついており、「不動点コンビネータを使え」と一括していた。TaPL輪読会真っ只中(最近再帰型を担当して発表した)なので不動点コンビネータとはそれなりに頻繁に触れ合っておりそこそこ身近に感じてはいたが、まさかこんなところでこれを使うコードが出てくると思っていなかったので見たときちょっと笑ってしまった。

この回答だ。https://stackoverflow.com/a/53504373

// non-recursive definition 
template<class T>
using Var = std::variant<int, bool, double, std::string, std::vector<T>>;

// tie the knot
template <template<class> class K>
struct Fix : K<Fix<K>>
{
   using K<Fix>::K;
};

using ScriptParameter = Fix<Var>;

しかも動くし。

wandbox.org

しかもこれは継承を使っているので、std::getなどの参照をとる関数にScriptParameterを直接渡せるし、メンバ関数はそのまま呼び出せる。面白いことを考える人がいるものだ。

アセンブリ解読:Fortran配列の中身

前回の続きです。

以下のコードを-O3コンパイルすると、

pure function nth(v, n) result(out)
    implicit none
    integer, intent(in) :: v(:)
    integer, intent(in) :: n
    integer :: out

    out = v(n)
    return
end function nth

こんなアセンブリになって驚いた。

nth_:
        mov     rdx, QWORD PTR [rdi+40]
        mov     eax, 1
        mov     rcx, QWORD PTR [rdi]
        test    rdx, rdx
        cmove   rdx, rax
        movsx   rax, DWORD PTR [rsi]
        imul    rax, rdx
        sub     rax, rdx
        mov     eax, DWORD PTR [rcx+rax*4]
        ret

参考: godbolt.org

これを解読していきたいと思う。

ここで、いくつか思い出さないといけないことがある。Fortranは多次元配列を組み込みで持っており、組み込み関数としてreshapeなんかも持っているということ。それと、Fortranでは特に指定されない限り配列の添字が1から始まることだ。これらを頭に入れて、このアセンブリを読んでいこう。

長いが、とりあえず何をしているかそのまま言葉にしてみる。

nth_:
        ; 第一引数から40バイト離れたところから64bitロードしてrdxに代入
        mov     rdx, QWORD PTR [rdi+40]
        ; raxの下位32bitに1を代入
        mov     eax, 1
        ; rcxに第一引数の先から64bitロード
        mov     rcx, QWORD PTR [rdi]
        ; rdxがゼロならZF=1に
        test    rdx, rdx
        ; ZF=1ならrdxにrax(今1が入っている)を代入
        cmove   rdx, rax
        ; raxに第二引数の先から32bitロードし、符号拡張して代入
        movsx   rax, DWORD PTR [rsi]
        ; rax *= rdx
        imul    rax, rdx
        ; rax -= rdx
        sub     rax, rdx
        ; rcxからrax*4バイト先の32bitをロードしてeaxに代入
        mov     eax, DWORD PTR [rcx+rax*4]
        ; return eax;
        ret

最後、戻り値レジスタeaxにアドレスの先からDWORD(32bit)ロードしてreturnしていることから、このときrcxが配列の先頭を指していて、raxが添字になっていることがわかる(4はintegerのサイズ)。raxには色々な計算がされているが、rcxは単にrdiの先からポインタがロードされているだけだ。ここから、第一引数は先頭ポインタを指すポインタであることがわかる。

同時に、1行めからrdi+40の先を読み込んでいることに注目すると、第一引数は単なるポインタではなく、情報がつめ込まれた何らかの構造体の先頭を指しており、その1つめの要素が配列の先頭ポインタであることがわかる。

struct FortranIntArray {
    Int* array; // rdiの指す先
    // ここに32バイト分何かが入っている
    // が、今回は使われていない
    QWORD unknown; // rdi+40の指す先
};

このrdi+40が指しているものは添字の計算で使われているので、この辺で配列の形状ではないかと予想がつく。これがどのように使われているか追いかけてみよう。一部抜粋する。

        ; rdx = FortranIntArray.unknown;
        mov     rdx, QWORD PTR [rdi+40]
        ; rax = 1
        mov     eax, 1
        ; if (rdx == 0) {ZF = 1;} else {ZF = 0}
        test    rdx, rdx
        ; if (ZF == 1) {rdx = rax;}
        cmove   rdx, rax

この時点で、もしunknown == 0ならrdx = 1になっている。この仮定(unknown == 0)を置いたまま、最後まで読んでみよう。

nth_:
        ; rdx = 1になったとする。
        cmove   rdx, rax
        ; rax に第二引数の値(添字)を符号拡張して代入
        movsx   rax, DWORD PTR [rsi]
        ; rax *= 1
        imul    rax, rdx
        ; rax -= 1
        sub     rax, rdx
        ; rcxからrax*4バイト先の32bitをロードしてeaxに代入
        mov     eax, DWORD PTR [rcx+rax*4]
        ; return eax;
        ret

この場合、渡した添字の一つ前を見ていることになる。最初に思い出したとおり、Fortranでは配列の添字が1から始まる。ポインタは先頭を指しているので0から始まることになる。この不一致を解決するために、1を引く必要がある。

では、unknown != 0だったケースを考えてみよう。1だったらどうなるだろうか。rdxは非ゼロなので1のままになる。っと、これは0の場合と同じだ。0のときはraxが代入されて1になるので。

じゃあ例えば、unknown == 2だったら?

        ; rdxが非ゼロなので、rdx == 2のまま
        cmove   rdx, rax
        ; rax に第二引数の値(添字)を符号拡張して代入
        movsx   rax, DWORD PTR [rsi]
        ; rax *= 2
        imul    rax, rdx
        ; rax -= 2
        sub     rax, rdx
        ; rcxからrax*4バイト先の32bitをロードしてeaxに代入
        mov     eax, DWORD PTR [rcx+rax*4]
        ; return eax;
        ret

添字が2倍されてから2を引かれた。1番めの要素の場合、オフセットは0。2番めの要素の場合、オフセットは2。3番めの要素の場合、オフセットは4。これは、第一カラムに何要素あるかがunknownに入っていて、その分をスキップしているように見える。Fortranはcolumn-majorなので、配列のメモリ上での配置は以下のようになるからだ(わかりやすさのため、添字を0スタートにしている)。

0, 2, 4
1, 3, 5

と、ここまでで大分このコードが何をしているのかわかったわけだが、それはそれとしてFortranは配列の次元が違っていたらコンパイルエラーにしてくれる。なので実際にはa(2,3)みたいな配列をnthに渡すことはできない。多次元配列を渡すときは、多次元配列であるという型情報が関数のインターフェースに載るのだ。だからこの関数にはかならず1次元配列が来る。そう考えるとこの処理には何の意味もないのではないかという気もする。

追記)後で思い出したが、Fortranにはスライスの時にストライドを設定できる。a(2:6:2)のようにすると、a(2), a(4), a(6)だけが使用可能になる。こういう場合だと、一次元配列なのにオフセットが必要になることがあり得る。(追記終わり

本当にオフセットが入っているのかわからなくなってきたので、無理やりCのコードをリンクして配列の先を見てみることにした。まずは一次元配列を見てみる。比較しやすいように、デフォルトのa(6)、3から始まるa(3:6)、0から始まるa(0:6)を見てみる。ここで、コンパイラを手元のものに変えたため、詳細なレイアウトは変わっているかも知れないことに注意したい。

integer :: a(6)
[rdi+0]  7ffc9ef56ee0
[rdi+8]  ffffffffffffffff
[rdi+16] 109
[rdi+24] 1
[rdi+32] 1
[rdi+40] 6
integer :: a2(3:6)
[rdi+0]  7ffc9ef56ed0
[rdi+8]  fffffffffffffffd
[rdi+16] 109
[rdi+24] 1
[rdi+32] 3
[rdi+40] 6
integer :: j3(0:6)
[rdi+0]  7fff7273e770
[rdi+8]  0
[rdi+16] 109
[rdi+24] 1
[rdi+32] 0
[rdi+40] 6

最初の部分は先頭ポインタだろう。2つめの64bitは添字から引き算するべき値に見える。配列が1から始まるときは-1(2の補数)で、3から始まるときは-3、0から始まるときは0だからだ。3つめは謎。4つめもよくわからない。5つめと6つめは最初と最後の添字だろう。これはわかりやすい。

二次元配列はどうだろう。

integer :: a(5,6)
[rdi+0]  7fff7273e810
[rdi+8]  fffffffffffffffa
[rdi+16] 10a
[rdi+24] 1
[rdi+32] 1
[rdi+40] 5
[rdi+48] 5
[rdi+56] 1
[rdi+64] 6

心の目に1,1,55,1,6という塊が見える。これは1次元配列の場合、a(6)1,1,6a(3:6)1,3,6a(0:6)1,0,6だったものに相当するのだろう。2次元のときに2つめの次元で5,1,6となっていて、最初の数字が5になっている。多分これは前までの次元の積で、この次元の配列を見に行こうとした時にスキップするべき塊のサイズを先に計算しているのではないか。

三次元配列を覗いてみる。a(5,6,7)としておけば、2次元めまでは同じ1,1,55,1,6で、3次元めで30,1,7になるだろう。

integer :: a(5,6,7)
[rdi+0]  7fff7273e420
[rdi+8]  ffffffffffffffdc
[rdi+16] 10b
[rdi+24] 1
[rdi+32] 1
[rdi+40] 5
[rdi+48] 5
[rdi+56] 1
[rdi+64] 6
[rdi+72] 1e
[rdi+80] 1
[rdi+88] 7

30は16 + 14なので16進で1e。そうなっているようだ。ところでポインタの次にある値は-36になっている。これは6*5+6だ。5x6x7の3次元配列の(1,1,1)番目の要素は、もし配列が0始まりなら5x6x1 + 6x1 + 1になる。この前半部分、5x6x1+6x1を前もって計算してあるのだろう。

あとは、配列の次元をどうやって記録しているかが知りたいところだ。これまでの情報から計算はできるが、毎回計算していたのでは遅かろう。どこかに即値があると思いたい。 よく見ると、[rdi+16]に入っている値が次元を上げると共に1増えていることに気づく。これっぽい。なぜそんな大きな値になっているのかはまだわからないが。


長くなってしまった。ここまでのことを踏まえると、最初の関数には配列の長さ情報を与えておくと命令数が2命令くらいになるに違いない。

pure function nth(v, n) result(out)
    implicit none
    integer, intent(in) :: v(3) ! 例えば
    integer, intent(in) :: n
    integer :: out

    out = v(n)
    return
end function nth

godbolt.org

なった。

nth_:
        movsx   rax, DWORD PTR [rsi]
        mov     eax, DWORD PTR [rdi-4+rax*4]
        ret

Fortranで関数に配列を渡すときは、もしわかっているのなら、配列の長さを陽に渡していたほうが速度面からよさそうに思える。

C++版と比較して少し不思議に思う所は、配列の添字をポインタで渡しているところだ。添字はintegerなので、明らかにレジスタに収まる。実際、C++で似たようなものを書くとnesiに乗った状態で(メモリアクセスせずに)持ってこられている。だがFortran版では[rsi]でメモリアクセスをしている。

int nth_int(const int* v, const std::int32_t n) noexcept
{
    return v[n];
}
nth_int(int const*, int):
        movsx   rsi, esi
        mov     eax, DWORD PTR [rdi+rsi*4]
        ret

ちょっと調べたところ、Fortran2003でVALUE属性というものが入ったらしい。これがあれば明示的に値渡しにできる。それ以外の場合は-O3でも参照渡しになるのはちょっとビビった。

www.nag-j.co.jp

試してみよう。

godbolt.org

値渡しになったようだ。レジスタの値がそのまま使われている。

pure function nth(v, n) result(out)
    implicit none
    integer, intent(in) :: v(3)
    integer, value, intent(in) :: n
    integer :: out

    out = v(n)
    return
end function nth
nth_:
        movsx   rsi, esi
        mov     eax, DWORD PTR [rdi-4+rsi*4]
        ret

細かなところだが意外と複雑で面白かった。わざわざFortranを書く状況なら速さが正義なはずなので、value属性や配列の長さに気をつけたほうがいいでしょうね。

x64で配列の添字にintを使うと遅い

こういうツイートを見た。気になったので試してみたら本当にそうらしい。

EAXとかRAXとかいうのはレジスタの名前だ。特にRAXは64ビットのレジスタで、その下位32ビットにEAXとしてアクセスできる。

 .-- RAX ------------.
 |        .----EAX---.
 |        |    .--AX-.
 |        |    |AH|AL|
64       32   16  8 bit

EAXに値をmovすると、上の32bitはゼロになる。だが、ゼロになると困る場合がある。32bitの符号付き整数を64bitに拡張して使いたいとしよう。 EAXに負の符号付き整数を入れた場合、例えば-1は0xFFFFFFFFだが、これを64bitレジスタに入れて上位ビットが0にされると0x0000'0000'FFFF'FFFFになってしまう。 すると今度は負の数でなく、2^32-1という値になってしまい、意味が変わってしまう。 そのため、CPUには符号を考慮してキャストするための命令が用意されており、この機能は符号拡張と呼ばれている。

intは多くの環境で32ビットであり、64ビット環境だとポインタが64ビットなので、配列の添字としてintを使おうとすると符号拡張命令が入って遅くなるということらしい。

確認してみよう。

godbolt.org

これが、

#include <vector>
#include <cstdint>

int nth_int(const std::vector<int>& v, const std::int32_t n) noexcept
{
    return v[n];
}

int nth_uint(const std::vector<int>& v, const std::uint32_t n) noexcept
{
    return v[n];
}

int nth_size_t(const std::vector<int>& v, const std::size_t n) noexcept
{
    return v[n];
}

int nth_uint64(const std::vector<int>& v, const std::uint64_t n) noexcept
{
    return v[n];
}
int nth_int64(const std::vector<int>& v, const std::int64_t n) noexcept
{
    return v[n];
}

こうなった。

nth_int(std::vector<int, std::allocator<int> > const&, int):
        mov     rax, QWORD PTR [rdi]
        movsx   rsi, esi
        mov     eax, DWORD PTR [rax+rsi*4]
        ret
nth_uint(std::vector<int, std::allocator<int> > const&, unsigned int):
        mov     rax, QWORD PTR [rdi]
        mov     esi, esi
        mov     eax, DWORD PTR [rax+rsi*4]
        ret
nth_size_t(std::vector<int, std::allocator<int> > const&, unsigned long):
        mov     rax, QWORD PTR [rdi]
        mov     eax, DWORD PTR [rax+rsi*4]
        ret
nth_uint64(std::vector<int, std::allocator<int> > const&, unsigned long):
        mov     rax, QWORD PTR [rdi]
        mov     eax, DWORD PTR [rax+rsi*4]
        ret
nth_int64(std::vector<int, std::allocator<int> > const&, long):
        mov     rax, QWORD PTR [rdi]
        mov     eax, DWORD PTR [rax+rsi*4]
        ret

確かに、int32_tを渡した時はmovsx(move sign extension)が呼ばれているし、そもそもesiからrsimovしていて無駄に見える。64ビット整数版より命令が一つ多い。

少しよくわからないのがuint32_tを渡した時で、符号なしなのでゼロで拡張してくれて良いのでは?と思うのだが、mov esi esiという謎の命令が入っている。なにこれ?意味なくね?

size_tは期待通り、直接アドレス操作に使われている。64bitの整数も同じくだ。

私はもとよりsize_tを使っているのだが、添字をキャッシュしておくようなコードを書いている時に32bitで足りるんだからメモリアクセスの効率を考えると32bitにした方がいいのかなあと考えて試したことがあった。 だがその時あまり速くならず、理由がよくわからないままそのブランチは切り捨てたのだが、もしかしたらここで追加されている命令のせいかもしれない。 おそらくその時足されていたであろう、uint32_tで追加される同じものを同じ所に代入している命令の意味はわからないが……。

ところで、私の知っているいくつかのFortran製のプロダクトでは皆単にintegerを使っていて、64ビットにしていることはほぼ無い。なのでこれは結構インパクトがあるのではと思って適当なコードを書いてみたら、なんだかよくわからなくなっていた。 なぜかは全くわからんが、命令数が異様に多い。C++版が2,3命令のところで9命令もある。-Ofastにしても変わらなかった。読み解くのも面倒になったので貼るだけにしておこう。あとで気が向いたら読む。

godbolt.org

thrustのバグを見つけて直した話

thrustがasync対応していたのを先月見つけた。

in-neuro.hatenablog.com

ドキュメントがまるで書かれていないところをコードを読み進めて試していたのだが、そこで不可解な動作を発見した。

以下のようなコードを書いたとする。thrustの関数を非同期に読んでそのfutureを返す関数だ。

#include <thrust/future.h>
#include <iostream>

// the definition is in a .cu file and will be linked later
thrust::system::cuda::unique_eager_future<int> do_something();

int main()
{
    auto f = do_something();
    std::cout << f.get() << std::endl;
    return 0;
}

これはコンパイルエラーになる。

$ g++-8 -std=c++11 -O2 -c main.cpp -I/usr/local/cuda/include/
In file included from /usr/local/cuda/include/thrust/system/cuda/detail/future.inl:26,
                 from /usr/local/cuda/include/thrust/system/cuda/future.h:71,
                 from /usr/local/cuda/include/thrust/future.h:55,
                 from main.cpp:1:
/usr/local/cuda/include/thrust/detail/event_error.h:110:19: error: explicit specialization of ‘template<class T> struct thrust::system::is_error_code_enum’ outside its namespace must use a nested-name-specifier [-fpermissive]
 template<> struct is_error_code_enum<event_errc> : true_type {};
                   ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

これは単純な話で、is_error_code_enumというトレイト構造体をevent_errcという型に対して特殊化しているのだが、その時にnamespaceを間違えていたのだ。構造体が所属する名前空間の外でその構造体を特殊化する時は、nested-name-specifierを明示する必要があり、今のthrustのバージョンではそれをしていなかった。

不思議なことに、これは翻訳単位を分けない限り気がつかない。同じファイルに関数の実体まで書いて完結させると問題なく通ってしまう。コンパイラの最適化周りの複雑な何かが絡んでいるのか規格通りの挙動なのかはよくわからないが、とにかくそうなってしまうので発見しづらい。テストコードはおそらく全て1ファイルで完結してしまっているのだろう。

というわけで、名前空間を補うだけのPull Reqを送った。サンプルコードも症状もバージョンも全部書いて、かつ問題点と解決策が自明なので一瞬でマージされるだろうと思ったら実際には1ヶ月近くかかった。

github.com

とはいえマージされたので、thrustのコントリビュータ一覧に名前が載った。小さなことではあるがそれなりに嬉しい。

tomlの複数行文字列内でのクオートの取り扱いに関するエッジケース

何が起きたのか

TOML規格レポジトリにこの記事のタイトルの通りのアップデートがあった。

github.com

基本的に、toml11ではTOMLのmasterブランチに入った変更は TOML11_USE_UNRELEASED_FEATURESdefine しないと有効にならないようにしている。だがこの手の"Clarification"の類は、もしtoml11がそれに則らない振る舞いをする場合、まだTOML側でリリースされていなくてもtoml11側の修正案件としている。これは「以前からこの振る舞いを想定していたが、これまで陽に書かれていなかっただけ」という解釈による。実際のところは最初から全てのエッジケースを想定していたと言うことはないだろうが。

このアップデートでは、以下のようなケースでの振る舞いが明確化された。

str = """"This," she said, "is just a pointless statement.""""
# 文字列の内容は"This," she said, "is just a pointless statement."
# になる(最初と最後の`"`は文字列の内容に含まれる)

何が問題だったのか

TOMLでは"""'''で始めた文字列はmultiline stringになる。そしてその場合は、"が文字列中に(エスケープなしで)現れることが許される。 ではこの場合、"""の直前や直後に"が来たらどうなるだろうか? 外側の"""をとるか内側の"""を取るかで、そのような文字列が文法的に正しいかが決まる。 以下に示すのは、文字列の最初と最後に"が来て欲しい比較的ありそうなケースだ。よく見ると、最初の"""が3つではなく4つになっている。

str = """"This," she said, "is just a pointless statement.""""

この場合のこの文字列の内容は "This," she said, "is just a pointless statement." が正しく、もっとも外側の """ デリミタが正しいデリミタの位置である。 このことが今までのspecだと曖昧だった。

toml11側の対応

これに今日気づいて、閉じる方の"""は最短マッチしてしまうと言うtoml11側の"不具合"を直した。"不具合"を"で括っているのは、v0.5.0時点の規格を厳密に参照するならここの動作は"処理系定義"に近いと感じているからだ。実際ClarificationはまだTOML仕様としてリリースされていないので。

あと、これに関する挙動を調べている時に、文字列のシリアライズ時にmultiline stringになるかどうかの場合分けとその場合のフォーマットがあまり美しくないなと感じたので、そこもあとで直す。

例えば、'で括られたliteral stringは一切のエスケープシーケンスを持たないので、literal stringの中で'を使おうとすると'''を使ってmultiline literal stringにするしかない。そこは正しく処理できるのだが、シリアライザで「multilineにするということは長いのだろう」という仮定がされており、開始即改行されてしまう。

str = '''
'That's still pointless', she said.'''

これはちょっと微妙だ。下のようになっていた方がカッコいい。

str = ''''That's still pointless', she said.'''

いやもしかしたら読みづらいかもな……?

あと、"で括られたbasic stringの場合は、"エスケープシーケンスを使うか全体を"""で括ってmultiline basic stringにするという2通りの選択肢がある。ここで常にエスケープシーケンスを使っていたが、バックスラッシュが多用されているとちょっと読みづらいので、少し考えるべきだろう。

TravisでOpenMPは使える

使えるみたいです。ただしOMP_NUM_THREADSを小さめの値、4とか、にしておかないといけないらしい。というのも、VMを動かしているハードのコア数(めちゃくちゃ多いんだと思う)を取得するからだそうで。CIに投げてる全員がそんな数のスレッドを立ち上げたらオーバーサブスクリプションどころの騒ぎじゃなくなるんだろう、多分。使われてる規模を考えると4くらいなら平気というのもすごい話ではある。

docs.travis-ci.com

それだけです。

OpenMPの挙動を見る

少し前、ちょっと面白い挙動に出くわしたので、順を追って説明しておこうと思った。知ってる人には当たり前のことが書かれているので下の方まで飛ばしてくれても構わない。

OpenMPはプラグマを書くだけでメモリ共有並列化ができるというやつで、手をつけやすいために広く使われている。例えば以下のような感じのコードが書ける。

double calc_nth_element(std::size_t n); // thread safe

int main()
{
    std::vector<double> v(N);
#pragma omp parallel for
    for(std::size_t i=0; i<N; ++i)
    {
        v[i] = calc_nth_element(i);
    }
}

メモリは共有されているので、ここで例えばvの同じ要素に書き込んだりすると詰む。

多くの資料でこのような単純な例は紹介されているが、若干複雑になってきた場合、例えばparallel領域から関数を読んで、その関数にもpragma omp forなどが付いている時にどう振る舞うのか、というようなことは紹介されていない。

parallel forからparallel forが書かれた関数を読んだ場合は紹介されていて、OMP_NESTEDtrueなら再帰的にスレッドグループが作られてしまい、ネストが深すぎるとすさまじい数のスレッドが立ち上がりシステムがハングするらしい。falseなら、単にシングルスレッドになるらしい。

だが、今気になっているのはparallel領域からparallelのついていないomp forだけが書かれた関数を読んだ場合だ。 並列化可能で互いに依存しないタスクが複数個あって、どれを実行するかが実行時にしか決まらない場合、そしてそれが複数ある場合、こういうことをしたくなる。 そういうタスクに同じインターフェースを持たせておいて、入力時に切り替えて、それをparallel領域内で呼び出すのだ。 普通のforparallel領域はそのブロックの終わりでjoinするので、大きなparallel領域内でomp for nowaitを使ってブロックせずに仕事が終わったスレッドから順に次の関数を実行できて欲しい。

だが、この手の用法を紹介している本や記事が異様に少ない、というか少し見た限りでは存在しない。GitHubopenmpを使っていそうなプロジェクトを調べてそのレポジトリ内のpragmaで検索しても誰も使っていない。 仕様書のforの部分をざっと見た限りでは、関数呼び出しについては何も書いていない。とりあえず動いているので使ったりしているが、心配なことに変わりはない。 何か知っている人がいたら教えて欲しい。

ここまで書いていて、nowaitありとなしでベンチマークを取っていないことに気付いた。あとでとっとこう。こんなことをする必要がないならしなければよいだけなので。

というわけでどうなるのかいっちょやってみよう。GCC 9.2.1を用意した。

: 一般に(実装が仕様を兼ねている言語を除いて)、言語仕様上保証されている動作と処理系の動作は必ずしも一致しない。その場合は処理系のバグとして解釈され、修正対象となる。多くの人が使っている処理系は、その「多くの人」の中に言語仕様を理解している人がいる可能性が高くなるため、そのようなバグは修正済みである確率が高く、言語仕様を調べる上でよいヒューリスティックになる。だが、言語仕様上処理系がどのように扱っても構わない場合というのも存在するため(処理系定義・未定義動作など)、言語処理系の動作から言語仕様、特にエッジケースでの挙動を理解したと考えるのは(多くの処理系が規格によく従っており複数の処理系が確立されている現代においては)基本的に愚行である。以下ではこの愚行を行っている点に注意していただきたい。

まず、どのスレッドから何が行われているかを表示するためにプリント関数を用意する。出力の順序がランダムにならないように、ロックを置いて毎回それを取得する。C++11を使っていたので再帰で出力しているが、C++17ならfold expressionを使おう。

std::mutex iomtx;

template<typename T1>
void print_impl(const T1& arg)
{
    std::cout << arg;
    return;
}

template<typename T1, typename T2, typename ... Ts>
void print_impl(const T1& arg1, const T2& arg2, const Ts& ... args)
{
    std::cout << arg1 << arg2;
    print_impl(args...);
    return;
}

template<typename ... Ts>
void print(Ts&& ... args)
{
    std::lock_guard<std::mutex> iolock(iomtx);
    std::cout << "from thread " << omp_get_thread_num() << ": ";

    print_impl(std::forward<Ts>(args) ...);
    std::cout << std::endl;
    return;
}

ではこれを使って色々どうなるか見てみよう。まず、普通にpragma omp parallel forをやってみる。あまりコア数が多いと見づらいので4スレッドに設定している。

void foo(const std::size_t n)
{
#pragma omp parallel for
    for(std::size_t i=0; i<n; ++i)
    {
        print("executing foo ", i, "-th loop");
    }
}

int main(int argc, char **argv)
{
    if(argc != 2)
    {
        return 1;
    }

    const std::size_t n = std::atoi(argv[1]);
    foo(n);
    return 0;
}
$ g++-9 -fopenmp -std=c++11 main.cpp
$ ./a.out 12
from thread 3: executing foo 9-th loop
from thread 3: executing foo 10-th loop
from thread 3: executing foo 11-th loop
from thread 2: executing foo 6-th loop
from thread 2: executing foo 7-th loop
from thread 2: executing foo 8-th loop
from thread 0: executing foo 0-th loop
from thread 0: executing foo 1-th loop
from thread 0: executing foo 2-th loop
from thread 1: executing foo 3-th loop
from thread 1: executing foo 4-th loop
from thread 1: executing foo 5-th loop

わかる。特に驚くようなことはない。

では少し書き換えて、parallel領域をmainの中に作り、pragma omp forだけが書かれたfooを呼んで見よう。

void foo(const std::size_t n)
{
    print("entered into foo");
#pragma omp for
    for(std::size_t i=0; i<n; ++i)
    {
        print("executing foo ", i, "-th loop");
    }
}

int main(int argc, char **argv)
{
    if(argc != 2)
    {
        return 1;
    }

    const std::size_t n = std::atoi(argv[1]);

#pragma omp parallel
    {
        foo(n);
    }
    return 0;
}

parallel領域では注釈がなければ全ての行を全スレッドが実行するので、全スレッドがfooを呼ぶ。出力は以下のようになる。

from thread 0: entered into foo
from thread 0: executing foo 0-th loop
from thread 0: executing foo 1-th loop
from thread 1: entered into foo
from thread 1: executing foo 3-th loop
from thread 3: entered into foo
from thread 3: executing foo 9-th loop
from thread 3: executing foo 10-th loop
from thread 3: executing foo 11-th loop
from thread 1: executing foo 4-th loop
from thread 1: executing foo 5-th loop
from thread 0: executing foo 2-th loop
from thread 2: entered into foo
from thread 2: executing foo 6-th loop
from thread 2: executing foo 7-th loop
from thread 2: executing foo 8-th loop

少し非直感的かも知れないが、期待通りの挙動をした。全スレッドがfooの中に入った後、for文の自分が担当する領域だけを実行して返ってきている。ループ全体が4回実行されるわけではない。

これは、あたかもfooがインライン展開されたかのような挙動だ。一応念の為、これがインライン展開されてからOpenMP化されたアーティフアクトでないことを確認する目的でインライン展開を抑制してみる。

void __attribute__((noinline)) foo(const std::size_t n)
{
    print("entered into foo");
#pragma omp for
    for(std::size_t i=0; i<n; ++i)
    {
        print("executing foo ", i, "-th loop");
    }
}

結果は同じだった。実行順序は違うが、それは並列化しているので当然です。

from thread 1: entered into foo
from thread 1: executing foo 3-th loop
from thread 1: executing foo 4-th loop
from thread 1: executing foo 5-th loop
from thread 2: entered into foo
from thread 2: executing foo 6-th loop
from thread 2: executing foo 7-th loop
from thread 0: entered into foo
from thread 0: executing foo 0-th loop
from thread 0: executing foo 1-th loop
from thread 3: entered into foo
from thread 3: executing foo 9-th loop
from thread 3: executing foo 10-th loop
from thread 3: executing foo 11-th loop
from thread 2: executing foo 8-th loop
from thread 0: executing foo 2-th loop

では、この関数をmaster領域から呼ぶとどうなるだろうか? master領域は、マスタースレッド(スレッド番号が0のスレッド)だけが実行することができる。そこからこの関数を呼ぶと?

void foo(const std::size_t n)
{
    print("entered into foo");
#pragma omp for
    for(std::size_t i=0; i<n; ++i)
    {
        print("executing foo ", i, "-th loop");
    }
}

int main(int argc, char **argv)
{
    if(argc != 2)
    {
        return 1;
    }

    const std::size_t n = std::atoi(argv[1]);

#pragma omp parallel
    {
#pragma omp master
        {
            foo(n);
        }
    }
    return 0;
}
$ ./a.out 12
from thread 0: entered into foo
from thread 0: executing foo 0-th loop
from thread 0: executing foo 1-th loop
from thread 0: executing foo 2-th loop
^C

0番は全体の1/4だけを実行する。残りの3/4は他のスレッドに割り当てられたタスクになるわけだが、この関数はmaster領域から呼ばれたため、0番以外のスレッドは動くことができない。というわけでこのプログラムはこのままハングしてしまった。

考えてみれば当たり前に思えるが、少し面白い挙動だ。

ちなみに、これを関数に分けずにコンパイルしようとすると、以下のようなエラーが発生する。

$ g++-9 -fopenmp -std=c++11 foo.cpp
foo.cpp: In function ‘int main(int, char**)’:
foo.cpp:59:9: error: work-sharing region may not be closely nested inside of work-sharing, ‘critical’, ‘ordered’, ‘master’, explicit ‘task’ or ‘taskloop’ region
   59 | #pragma omp for
      |         ^~~

もうちょっとだけ静的解析を頑張ってくれていれば関数越しにでもエラーを出せそうなんだが。 この挙動はどうにもコンパイラ側が、pragma omp forが単独で(omp parallel領域から呼ばれることを想定した)関数内に登場することを初めから考えていないようにも見える。

一つの関数の中でparallel領域を閉じたくない場合、OpenMPtaskを明示的に書いたり、あるいはTBBやThrust、cpp-taskflowなどのライブラリを使ったほうがよさそうな気がしてきた。

とはいえ、これはいち処理系での挙動なので(intelコンパイラでも同じく固まったし、この挙動に筋も通るが)、良い子の皆さんはいち処理系の挙動を調べて満足することなく仕様書を精読しましょう。

細かい状況でOpenMPがどう振る舞うのか概説書なんかを呼んでもよくわからないので、どこかで時間を見つけて仕様書をちゃんと読みたい……。

OpenMPの仕様書はここで見れます: www.openmp.org