OpenMPプログラムは，簡単にマルチスレッドを実現する仕組みです．OpenMPを有効にしてコンパイルしたプログラムは，起動するときに，環境変数OMP_NUM_THREADSを読み取り，その数値まではマルチスレッドが可能であるという条件でプログラムを実行します．具体的には, OpenMPに対応したプログラムを，

OSX:

まだサポートしやがらねえ．いい加減にしろApple　たぶんこいつら, C++なんかやめてSwift使えと誘導したいんだろうけど

Linux:

g++ -std=c++2a -fopenmp my_program.cpp -o ./my_program

CMake: CMakeLists.txtに以下を追加

find_package(OpenMP)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")

のようにOpenMP有効にしてコンパイルし，

OMP_NUM_THREAD=4 ./my_program

と起動すると, my_programが4スレッド利用可能となります．

以下では, 「OpenMPに対応したプログラム」の作成方法を示します．

絶対必要なもの

#include <omp.h>

OpenMP最初の一歩-FOR並列

OpenMPの利用は極楽です．並列化したいfor文があったとします．これをまず，整数で数えるタイプのfor文に書き換えます：

for(int i=0; i<data.size();i++){.....}

range-based forループは，2025年ごろに並列化できると聞きました．

さて，このfor文が, iを独立に計算してもOKであったとしましょう（ここの判断は，プログラマーの責任です．OpenMPは．「ここを並列に実行するとマズイ」とかの判断を行いません．計算結果がどうなろうと，並列に行えと指示されたら，並列に実行する，だけです）

書き換え後：

#pragma omp parallel for
for(int i=0; i<data.size();i++){.....}

pragma omp parallel for直後のfor文が，並列に実行されます．

実行の際には，OMP_NUM_THREAD環境変数で並列数を指定します：

sugimoto@sun3[22] $ OMP_NUM_THREADS=12 ./a.out

これで並列実行できます．

注意: 並列実行前に定義されていた変数は, すべてのスレッドで共有になる. 従って

int j;
.....
#pragma omp parallel for
for(i=0; i<40; i++){
    j=......  ループ変数i の値に依存する数式
    f[i]=....jの値に依存する数式.....
}

変数 j は全スレッドで共有されるので, 第二の式で, どの i に対応する j の値が用いられるのかは不定になる. 上の例では, 変数 j を，意味もなくプログラムの先頭で宣言したのが原因である. C++では, 変数をプログラムの先頭で定義すると，悪運しか招かない．変数は必要なところで「のみ」宣言するのが正しい:

.....
#pragma omp parallel for
for(i=0; i<40; i++){
    int j=......  i の値に依存する数式
    f[i]=....j.....
}

それでも，なおかつ，グローバル変数に値を代入しなければならないのであれば，private宣言を行う：

int j;
.....
#pragma omp parallel for private(j)
for(i=0; i<40; i++){
    j=......  i の値に依存する数式
    f[i]=....j.....
}

OpenMPでは, ループの並列化に必要な条件は次のものです：

ループを支配する繰り返し指定子 [ for(int i=0;....) とかの i ] をランダムにアクセスできる
- そりゃ, 1スレッド目が i=0 やってるから，隣のCoreが i=10 をやりたいので，i=10 に「サクッと」アクセス可能でないとだめ
- このため，イテレータを利用する奴はダメ. for(auto i=set.begin(); i<set.end(); i++) みたいなのだと，10番目にアクセスするために i++ を10回やらないとダメだからね．
サクッとアクセスできるなら，OpenMP可能である．といっても，それはやはり for(int i=0; i<1000; ...) みたいなのに限定されるわけだよ．
スレッド数の選択は重要です．例えば25セルの計算を12threadでやったら, 25=3x12-11なので16時間かかってしまい，13threadでやったら25=2x13-1であるので10時間で経産完了しまぴた．

OpenMPつぎの一歩 -SECTION並列

forループに乏しい場合でも, もうすこし利用できる場合はあります. たとえば

void my_function(){
    ....
    {
       some work
    }
    ....
    small work
    ....
    {
       other work
    }
    {
       yet another work
    }
    ....
}

some-work, other-work, yet-another-workは依存関係がないとします．これを並列実行するには

void my_function(){
    ....
#pragma omp parallel sections num_threads(3)  ←ここで 3 スレッドの並走が始まる
{
    #pragma omp section
    {
       some work    ←最初のお兄ちゃんが取り組む．残り2スレッドは何もせず通過 
    }
    small workは{}からはみ出しているので, 全員が取り組む
    ....
    small work
    ....
    #pragma omp section
    {
       other work   ←2番目のお兄ちゃんが取り組む．最初のお兄ちゃんは上の仕事が終わったらここは通過
    }
    #pragma omp section
    {
       yet another work  ←末っ子が取り組む．上のお兄ちゃんたちはここをスキップ
    }
}                    ←ここで1スレッドを残して死滅する
    ....
}

これで分担できます．

OpenMP-TASK並列

sectionsでは, コードの一部を振り分けていました．多数のループを振り分けるには, TASKを利用します．

void my_function(){
    ....
for(auto &g:Set) {
}
#pragma omp parallel sections num_threads(3)  ←ここで 3 スレッドの並走が始まる
{
    #pragma omp section
    {
       some work    ←最初のお兄ちゃんが取り組む．残り2スレッドは何もせず通過 
    }
    ここは{}からはみ出しているので, 全員が取り組む
    #pragma omp section
    {
       other work   ←2番目のお兄ちゃんが取り組む．最初のお兄ちゃんは上の仕事が終わったらここは通過
    }
    #pragma omp section
    {
       yet another work  ←末っ子が取り組む．上のお兄ちゃんたちはここをスキップ
    }
}                    ←ここで1スレッドを残して死滅する
    ....
}

で

OpenMPの便利な知識

ここでは，並列化の際に困る色々な事態に対処する方法を述べます．

CRITICAL構文

parallelしてる途中に，例えば画面にメッセージを書くと, 全スレッドがバラバラに画面に書きます. だからHelloと書くと，画面には HeHlelolwlow と出力されます．そこで，並列中なんだけど，ここは1スレッドづつやってねーと指示できます．

#pragma omp parallel for
for(int i=.....){
    絶賛並列計算中
    #pragma tmp critical
    {   次のセクションは同時には1スレッドしか行わない
       if (result!=0) std::cout << "Fucking situation occurred" << std::end;
    }
    絶賛並列計算中
}

REDUTION構文

並列計算して合計とか出したいですよね．それはreduction 構文で可能になります.

double sum0=0.0, sum1=0.0;
#pragma omp parallel for reduction(+:sum0,sum1)
for(int i=.....){
   sum0+=Value.at(i)*Weight_0[i];
   sum0+=Value.at(i)*Weight_1[i];
}

上のプログラムでは, for文で並列計算が始まると, reductionで＋が指定されている変数sum0, sum1が各スレッドに確保され0に初期化されます. そして,

   sum0+=Value.at(i)*Weight_0[i];
   sum0+=Value.at(i)*Weight_1[i];

は各スレッドの変数に行われます. 並列ループが終了すると, 各スレッドの結果sum0, sum1の合計が, ユーザーが確保した変数sum0, sum1に代入されます．これをreduction指示なしで実行すると,

   スレッド0: sum0の値を取得
   スレッド1: sum0の値を取得
   スレッド0: 計算結果をsum0に代入
   スレッド2: sum0の値を取得
   スレッド1: 計算結果をsum0に上書き　　← スレッド0の作業が台無し
   スレッド3: sum0の値を取得

このような作業になり，結果がデタラメになります．演算子と初期値のペアは次のように定義されています

+	-	*	&	\|	^	&&	\|\|
0	0	1	~0	0	0	1	0

また, reduction構文は演算子以外に手続に対しても有効なものがあります. 使いそうなものとしては,

min : 並列ループの最後に各スレッドの結果の最小値をとる.

max: 並列ループの最後に各スレッドの結果の最大値をとる.

があります. 初期値についてはminは変数の型の表現できる最大の数, maxは変数の型の表現できる最小の数となります.

ゆうざあでふぁいんどreduction

OpenMP4.0から，user-defined reductionが使えます．いやそれなんのためにあるのか？

例えばある配列の最大値と最小値を求めることにしましょう．すると

double Amin=DBL_MAX,Amax=-DBL_MAX;
size_t Pmin=0,Pmax=0;
#pragma omp parallel for reduction(min:Amin) reduction(max:Amax)
for(size_t i=0;i<ARRAY.size();i++)
{
if (ARRAY[i]<Amin) {Amin=ARRAY[i];Pmin=i;}
if (ARRAY[i]>Amax) {Amax=ARRAY[i];Pmax=i;}
}
std::cout << "(OpenMP)Min value " << Amin << " at i= " << Pmin << std::endl;
std::cout << "(OpenMP)Max value " << Amax << " at i= " << Pmax << std::endl;

こんな感じで実行しますと,

(serial)Min value -0.99999999987699972337 at i= 2163921589 <-- 念の為並列化しないやつ
(serial)Max value 0.99999999987684240477 at i= 2251222562 <-- 念の為並列化しないやつ
(OpenMP)Min value -0.99999999987699972337 at i= 1026610216
(OpenMP)Max value 0.99999999987684240477 at i= 14200192

ふむ．値はあってます. が, 場所が計算ミスですね！と，申しますのも，こいつらreductionは, AminとAmaxの面倒しか見ませんので・・・どこで最大だったかとかのPmin, Pmaxは放置なんですね.

じゃ，どうすりゃいいんだ．Pmin, Pmax は最大値を取るわけではないので・・・ええっとえええっと．困ってしまいます．

そこで，おいらのいう通りのreductionを定義したくなるわけです．こいつの利点は

max/min演算の結果として，インスタンスの代入が使える

ので, 値と一緒にindexも記録するクラスでも作ってしまえばOK. 具体的には

class AssHole
{
public:
double value;
size_t pos;
AssHole(double X){value=X;}
};
#pragma omp declare reduction(my_min:class AssHole:omp_out=(omp_in.value<omp_out.value)?omp_in:omp_out) initializer(omp_priv=omp_orig)

と準備します．AssHoleクラスは説明不要ですよね! で, pragma呪文ですが, これがあんたが定義したreduction演算の名前. すると, 並列計算を指示するやつが

#pragma omp parallel for reduction(あんたが定義したreduction演算:計算してポチい変数)

omp_out は予約語で「計算結果を代入するとこ」を表します．で，計算してポチい変数として「ループでやってくる値」を表す予約語がomp_in. ここでは, valueを比較して小さいならomp_inをomp_outに代入する．そうでなければ, omp_outをomp_outに代入する.(インスタンスを丸ごとコピーすることになっている)

謎めいているinitializerですが，最小値を計算するので, 初期値を設定してから比較しないとダメなんで，omp_priv=omp_orig と書くと，計算してポチい変数が最初に持っている値 omp_orig が，各スレッドの一時変数omp_privにコピーされて計算されます．というわけで，

AssHole XXX(DBL_MAX); <--- でかい値が, XXX.value の初期値になる
#pragma omp parallel for reduction(my_min:XXX)
for(size_t i=0;i<ARRAY.size();i++)
{
if (ARRAY[i].value<XXX.value) {XXX.value=ARRAY[i]; XXX.pos=i;}
}
std::cout << "(OMP-2)Min value " << XXX.value << " at i= " << XXX.pos << std::endl;

となるわけです．すると，

(serial)Min value -0.99999999987699972337 at i= 2163921589
(OpenMP)Min value -0.99999999987699972337 at i= 790007850　　　＜ーー間違えた
(OMP-2)Min value -0.99999999987699972337 at i= 2163921589　　　＜ーーちゃんとあってる！

間違う時には，結果がCPUの気分というか，そのコアがその時やっている雑用の量に依存しますので，結果が毎回変わるのが味噌ですね．

OpenMPバージョン

いろいろな都合で, OpenMPのバージョンが知りたくなった場合, _OPENMP ディレクティブを使いましょう.

std::string openmp_version()
{
#ifdef _OPENMP
switch (_OPENMP)
{
case 201611:return "OpenMP5.0TR1";
case 201511:return "OpenMP4.5";
case 201307:return "OpenMP4.0";
case 201107:return "OpenMP3.1";
case 200805:return "OpenMP3.0";
case 200505:return "OpenMP2.5";
default: return "OpenMP(unknown)";
}
#else
return "NoOpenMP";
#endif
}

デバッグは-Wunknown-pragmas

OSやコンパイラーでサポートされていない機能は使えません．これを調べるには, -Wall をつけてビルドしてみましょう；

as2022-1 $ g++ -Ofast reduction.cpp -fopenmp -o red -Wall
reduction.cpp:23: 警告: #pragma omp rarallel を無視します [-Wunknown-pragmas]
#pragma omp rarallel for reduction(min:Amin,max:Amax)

なんと, reduction(min:...) がサポートされていない? CentOS8 Stream で gcc Version 8.5.0, libgomp-8.5.0 で，

OpenMP4.5までサポートされている(C/C++でも min, maxが使える，と書いてある)

はずなんですけど，ねえ．

これは実は, pragma omp rarallel for reduction(min:Amin,max:Amax) ってのがバグっているので, pragma omp rarallel for reduction(min:Amin) reduction(max:Amax) と書けば良いのですね．

これいいな

C++OpenMP