# 誤り訂正シミュレーションのためのスレッド並列実践

大木です。質問が多くなったのでやりかたまとめます。
言語はC/C++を想定。サンプルコードはmain関数内に書いているけど外部関数でも大丈夫です。

## スレッド並列とは??
CPUの各スレッドにプログラムの作業を割り当てて、計算時間の効率化をする手法です。普通にプログラムを書くと、１スレッド内で計算されますが、現代のCPUはスレッド数が多いので計算資源を損しているのです。
CPUのスレッド数を比較してみましょう。
- Core i7-1068NG7(MacBook Pro 13インチのCPU上位モデル) : 4コア8スレッド
- AMD CPU Ryzen Threadripper 2990WX (研究室にある計算機) : 32コア64スレッド
- AMD CPU Ryzen Threadripper 3990WX (藤井二冠が使ってるという噂) : 64コア128スレッド

八嶋研には個人では所有しないようなPCが複数台あり、上手く使うことによって計算効率を数十倍にまで高めることができます。スリッパを個人所有している藤井二冠はぶっ飛んでますが...

## スレッド並列の手法
モンテカルロ・シミュレーションで用いそうなもののみ解説します。詳細はggってください。

### 基本構文
`omp.h`をインクルージョンして使いたいforループの前の行に`#pragma omp parallel for` と記述するだけです。簡単ですね。

#include <omp.h>

int main()
{

	#pragma omp parallel for
	for(int i=0; i<1000: i++){
		//　この内部のループが並列実行される
	}
}

使用するスレッド数も簡単に指定することができます。K個のスレッドで計算したいときは`#pragma omp parallel for num_threads(K)` としてあげればよいです。指定しない場合はCPUのスレッド全てを使って計算します。研究室共同の計算機ですので、スレッド数は適切な数に指定するようにしましょう。

### コンパイル方法
`gcc`の場合、`-fopenmp`コマンドをつけます。最適化オプション`-O3`などもつけとくと良いでしょう。

$ gcc -fopenmp -O3 hoge.c -o hoge
$ g++ -fopenmp -O3 hoge.cpp -o hoge

### シミュレーションする際の問題と解決手法

#### 問題:並列ループ内のカウント処理
たとえば以下のようなソースコードがあったとします。

#include <omp.h>
#include <iostream>

int main()
{
	int a = 0;
	#pragma omp parallel for
	for(int i=1; i<=100: i++){
		a += i+1;
	}
	cout << a << endl;
}

`1〜100`の数字をたすだけのプログラムです。`5050` と出るはずですが、並列化するとうまく行かないです。この原因は、各スレッドで`a`が非同期的に足されるからです。
並列の場合、変数`a`のメモリへの書き込みがスレッドごとで競合するため注意しなければなりません（読み込みだけなら大丈夫です）。
これと同じことが誤り率の計測でも起こります。

#### 解決方法① : 各スレッドで誤り率を配列にして最後に足し合わせる
各スレッドごとに変数を用意し、それをループ後(並列ループ外)に足し合わせれば良いという発想です。僕はこれでやってます。`omp_get_thread_num()` を用いて現在いるスレッド番号を取得するのがポイントです（スレッド１なら`1`）。

ブロック誤り率のシミュレーションを想定します。ビット誤り率計測も同様の処理で行います。

#include <omp.h>
#include <iostream>
#include <vector>
#define THREAD_NUM 10 // 使うスレッド数
#define NUMBER_OF_SIMULATE 1000 // シミュレーション回数

int main()
{

	vector<double> bler_array(THREAD_NUM); // ブロック誤り率をスレッドごとに格納
	double bler_avg = 0;
	
	#pragma omp parallel for num_threads(THREAD_NUM)
	for(int sim=0; sim<=NUMBER_OF_SIMULATE: sim++){ // simulate loop
		// Encode
		// Plus Noise
		// Decode
		if( estimated_info_bits != info_bits ){ // ブロック誤り発生
			bler_array[ omp_get_thread_num() ] += 1.0;
		}
	} // simulate loop end
	
	for(int i=0; i<K; i++){
		bler_avg += bler_array[i];
	}
	bler_avg /= NUMBER_OF_SIMULATE;
	
}

#### 解決方法② : reduction演算子を使う
`#pragma omp parallel for reduction(+:x)` で変数`x`を各スレッドで定義し、自動的に足し合わせてくれます。解決方法①を自動化したような形です。僕は怖いので使ってませんんが、ちゃんとやってくれるはずです。

#include <omp.h>
#include <iostream>
#include <vector>
#define THREAD_NUM 10 // 使うスレッド数
#define NUMBER_OF_SIMULATE 1000 // シミュレーション回数

int main()
{

	vector<double> bler_array(THREAD_NUM); // ブロック誤り率をスレッドごとに格納
	double bler = 0;
	
	#pragma omp parallel for reduction(+:bler) num_threads(THREAD_NUM)
	for(int sim=0; sim<=NUMBER_OF_SIMULATE: sim++){ // simulate loop
		// Encode
		// Plus Noise
		// Decode
		if( estimated_info_bits != info_bits ){ // ブロック誤り発生
			bler += 1.0;
		}
	} // simulate loop end
	
	bler /= NUMBER_OF_SIMULATE;
}

## 参考文献
OpenMPによるスレッド並列計算
http://www.eccse.kobe-u.ac.jp/assets/images/simulation_school/kobe-hpc-summer-basic-2018/KHPC2018-openmp.pdf