C++並列処理: openmp, thread, Eigen

# C++並列処理: openmp, thread, Eigen
毎回忘れてしまうC++の並列処理のメモ. 
よく使うやつだけ書きます. 

## Openmp
おそらく最も単純で誰でもすぐに使える並列化の方法. 
コンパイル時に, `-fopenmp`オプションをつける必要がある. 

### forを並列化する

#### 書き込み変数がスカラーの場合
***ベクトルの内積***を求める例で示す. 

```c++
double dot(const vector<double> &v1, const vector<double> &v2) {
	double ans = 0;
	#pragma omp parallel 
	{
		#pragma omp for reduction(+:ans)
		for (size_t i = 0; i < v1.size(); ++i){
			ans += v1[i]*v2[i];
		}
	}
	return ans;
}
```

- `#pragma omp parallel`で並列化対象を囲む
- `#pragma omp for`でforを並列化する
- CPUが同時に書き換える恐れのある変数`ans`は, `reduction`としてスレッドごとに別々に扱うようにする
- 最終的に, スレッドごとの`ans`を`+`でまとめればよいため, `reduction(+:ans)`と書く


#### 書き込み変数がVectorの場合
同時に書き換える恐れのある変数がvectorの場合は少し面倒になる. 
***ベクトル同士の要素ごとの積***を例に示す. 

```c++
#pragma omp declare reduction(vec_double_plus : std::vector<double> : \
                              std::transform(omp_out.begin(), omp_out.end(), omp_in.begin(), omp_out.begin(), std::plus<double>())) \
                    initializer(omp_priv = omp_orig)
vector<double> cwiseProduct(const vector<double> &v1, const vector<double> &v2) {
	vector<double> ans(v1.size(),0);
	#pragma omp parallel 
	{
		#pragma omp for reduction(vec_double_plus :ans)
		for (size_t i=0;i<v1.size();++i){
			ans[i]=v1[i]*v2[i];
		}
	}
    return ans;
}
```

- `#pragma omp declare reduction`でreductionで使う演算について定義する
- `vec_double_plus : std::vector<double>`で, `vector<double>`を引数にとる演算`vec_double_plus`を定義する(2つのベクトルの要素ごとの和)
   - `std::transform`には`#include <algorithm>`が必要
- `#pragma omp for reduction(vec_double_plus :ans)`で定義した演算を使う

**追記**: 
- 要素`ans[i]`自体に重複書き込みする可能性がなければ、reductionの記述は不要
    - この例では不要だった
- `#pragma omp parallel for`と1行で書くこともできる

#### 書き込み変数がMatrixの場合
EigenのMatrixXdを例に書く. 

```c++
MatrixXd A=MatrixXd::Zero(100,100);
#pragma omp declare reduction(+ : Eigen::MatrixXd : omp_out=omp_out+omp_in) initializer(omp_priv = omp_orig)
#pragma omp parallel for reduction(+:A)
for (int i=0;i<10;++i){
	A+=MatrixXd::Random(100,100);
}
```

- Eigenはoperator`+`が定義されているが, `#pragma omp declare`しないといけない
- `std::vector`よりもこちらのほうがわかりやすいのは, operator`+`がEigenで定義されているため
    - そのため, `std::vector`の方も演算子`+`を定義してからの方がわかりやすいかもしれない
- `#pragma omp parallel for`は`#pragma omp parallel`と`#pragma omp for`の省略表現


### タスクの並列化
複数のタスクを同時に実行したいときに有効. 

```c++
#pragma omp parallel
{
		#pragma omp sections
		{
				#pragma omp section
				{
					dot(v1,v2);
				}
				#pragma omp section
				{
					cwiseProduct(v1,v2);
				}
		}
}
```

- `#pragma omp parallel`で並列化範囲を宣言
- `#pragma omp sections`でセクションごとで並列化することを宣言
- `#pragma omp section`で1つのセクションを宣言する

## Thread
コンパイル時に`-pthread`が必要. 
`#include<thread>`が必要. 

### 戻り値なしの場合
以下の例では, 関数`addRef`の実行N回を別プロセスで行う. 
結果を受け取るには参照渡しをするしかない. 

```c++
void addRef(int x, int y,int &ans){
	ans = x + y;
}
int main() {
	const int N = 5;
	const int x=1, y=2;
	
	vector<thread> threads;
	vector<int> result(N);
	
	for (int i=0;i<N;++i)
		threads.push_back(thread(addRef, x, y, ref(result[i])));
	for (thread &t: threads)
		t.join();
	for (int r : result)
		cout<<r<<endl;
    return 0;
}
```

- 最初のforでスレッド作成: `thread(関数名, 引数1, 引数2, ...)`
    - `ref`を使わないと参照渡しできない
- 次のforでスレッド終了待ち: `join()`
- 最後のforで結果出力


### 戻り値ありの場合
ラムダ式で書けば受け取れる. 

```c++
int addReturn(int x, int y){
	return x + y;
}
int main() {
	const int N = 5;
	const int x=1, y=2;
	
	vector<thread> threads;
	vector<int> result(N);
	
	for (int i=0;i<N;++i)
		threads.push_back(thread([i, &result]{ result[i] = addReturn(x, y); }));
	for (thread &t: threads)
		t.join();
	for (int r:result)
		cout<<r<<endl;
    return 0;
}
```

- 最初のforをlambda式で覆っただけ
- [async](https://qiita.com/termoshtt/items/d3cb7fe226cdd498d2ef)を使ってもいい

## Eigen
[Eigen](http://eigen.tuxfamily.org/index.php?title=Main_Page)は, 行列を扱うライブラリ. 
ヘッダーをインクルードするだけで使える: (例)`g++ -I./Eigen test.cpp`と`#include <Eigen/Dense>`
クイックリファレンス: [AsciiQuickReference.txt](https://eigen.tuxfamily.org/dox/AsciiQuickReference.txt)

Eigenは`-fopenmp`オプションをつけるだけで自動的に並列化される. 
公式ドキュメント[Make Eigen run in parallel](https://eigen.tuxfamily.org/dox/TopicMultiThreading.html)を見ると, 以下の機能は並列化に対応しているようだ. 

- general dense matrix - matrix products
- PartialPivLU
- row-major-sparse * dense vector/matrix products
- ConjugateGradient with Lower|Upper as the UpLo template parameter.
- BiCGSTAB with a row-major sparse matrix format.
- LeastSquaresConjugateGradient


# 参考文献

- [C/C++でOpenMPを試してみる](http://onoono-life.blogspot.com/2014/04/ccopenmp.html)
- [Openmp and reduction on std::vector?](https://stackoverflow.com/questions/43168661/openmp-and-reduction-on-stdvector)
- [C++で簡単非同期処理](https://qiita.com/termoshtt/items/d3cb7fe226cdd498d2ef)
- [Make Eigen run in parallel
](https://eigen.tuxfamily.org/dox/TopicMultiThreading.html)