OpenMP超入門 その1


## OpenMPとは?
OpenMPのことをWikipediaさんで調べると...
> OpenMPは、並列コンピューティング環境を利用するために用いられる標準化された基盤。OpenMPは主に共有メモリ型並列計算機で用いられる。[OpenMP:Wikipedia](https://ja.wikipedia.org/wiki/OpenMP)

とかかれてますね... 

なんのこっちゃい, というと, 
**並列計算するともっとはやくなるお!! 拙者を使えばメモリも共有して並列計算できるから便利だお**
というわけです!!

対応言語としてはC/C++, fortranが対応しています!!

## OpenMPを使う羽目になったきっかけ
実はPythonで並列化できるだろうと思っていたら, メモリ共有に関して自分で適宜したClassではむりそうだったからです...

ホントは全部Pythonで書いてJupyterにまとめようと思ったのですが僕の見通しが悪かったみたいです.

## 環境
1. Windows
    * Cygwinで確認できてます!!
1. Mac
    * 悪いこといわないからbrew install gccとかしてgccに切り替えよう(なぜかMacではgccがclangのエイリアスになってる...)
1. Linux系統
    * 確認中(たぶんいける)

***

### すべての環境に関していえること

` gccは4.2以降 `

というか, 新しいバージョンであれば大丈夫, という認識でOKです...

## Step1 : OpenMPが動くか確認
今回, 実行についてC++, g++ 5,4,0でやっております...
### Hello World with OpenMP
比較的どこのサイトでもやっているやつ

```cpp
#include <iostream>
#include <omp.h> // OpenMPのヘッダをinclude

int main(){
//ここが並列処理される
  #pragma omp parallel
  {
    std::cout << "Hello World!\n";
  }
}
```

これを実行すると

```
Hello World!
Hello World!
Hello World!
Hello World!
```

と表示される. このままだと, 異なるスレッドに全く同じ処理を任せていることになる. 
### 配列を操作してみる

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

int main(){
  int A[10];
  #pragma omp parallel for
  for(int i = 0; i < 10; i++){
    A[i] = i;
  }

  for(int i = 0; i < 10; i++){
    std::cout << A[i] << "\n";
  }
}
```

上のようなプログラムがあったとすると, 

```
0
1
2
3
4
5
6
7
8
9
```

と表示してくれる. もちろん, 並列化されている. 
### 単純に並列化してみる
今から, 1,...1000を足してみる. 500500が答えだ. 
普通に書くと

```cpp
#include <iostream>

int main(){
  int sum = 0;
  for(int i = 0; i < 1001; i++){
    sum += i;
  }

  std::cout << "sum = " << sum <<"\n";

}
```

とかくはずだ. もちろん, これは正しい. しかし, これを並列化する際に

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

int main(){
  int sum = 0;
  #pragma omp parallel for
  for(int i = 0; i < 1001; i++){
    sum += i;
  }

  std::cout << "sum = " << sum <<"\n";

}
```

と書いてしまうと, 安定しない. 500500にならないときもある. 
正しくは, 

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

int main(){
  int sum = 0;
  #pragma omp parallel for reduction(+:sum)
  for(int i = 0; i < 1001; i++){
    sum += i;
  }

  std::cout << "sum = " << sum <<"\n";

}
```

ひとつの変数に足しこむときはreductionでその変数を設定する必要がある.

### もっと実用的っぽいコードを書いてみる...
ということで数値積分をやってみます...
*e*を求めてみましょう!!

```math
e - 1 = \int_{0}^1 e^x dx
```

なので, これを元に求めてみます. 
#### 逐次計算だと...

```cpp
#include <iostream>
#include <cmath>

#define N 10000

int main(){
  double step = 1.0 / N;
  double temp;
  double ans;
  for(int i = 0; i < N; i++){
    temp = (i + 0.5) * step;
    ans += exp(temp)/N;
  }
  std::cout << "ans = " << ans + 1 << "\n";
}
```

とかけます. 結果は

```
ans = 2.71828
```

です. 
#### 並列計算でやると...
***
いちばんやっちゃいけないやつ

```cpp
#include <iostream>
#include <cmath>
#include <omp.h>
#define N 10000

int main(){
  double step = 1.0 / N;
  double temp;
  double ans;
  #pragma omp parallel for
  for(int i = 0; i < N; i++){
    temp = (i + 0.5) * step;
    ans += exp(temp)/N;
  }
  std::cout << "ans = " << ans + 1 << "\n";
}
```

動かすと

```
ans = 2.0528
```

もちろん, 間違ってます.
***
惜しい!!って感じの間違え

```cpp
#include <iostream>
#include <cmath>
#include <omp.h>
#define N 10000

int main(){
  double step = 1.0 / N;
  double temp;
  double ans;
  #pragma omp parallel for reduction(+:ans)
  for(int i = 0; i < N; i++){
    temp = (i + 0.5) * step;
    ans += exp(temp)/N;
  }
  std::cout << "ans = " << ans + 1 << "\n";
}
```

~~まぁ, 動くには動くけど, お作法的にはお行儀がよくない. ちゃんと, Loop内で書き換えられる変数(用は使い捨ての変数とか)はちゃんとprivateしてあげたほうがいい.~~ (@yohhoy 様の指摘ありましたので変更します. )

実は, このままだと変数`temp`はshared変数ではないため, 実行に際して動作の保証ができないようです. そのため, 以下のようにしましょう...
***
ベストアンサー

```cpp
#include <iostream>
#include <cmath>
#include <omp.h>
#define N 10000

int main(){
  double step = 1.0 / N;
  double temp;
  double ans;
  #pragma omp parallel for reduction(+:ans) private(temp)
  for(int i = 0; i < N; i++){
    temp = (i + 0.5) * step;
    ans += exp(temp)/N;
  }
  std::cout << "ans = " << ans + 1 << "\n";
}
```

こう書いてあげたほうがいい. (パフォーマンスに響くんだとか)

## 最後に
ほんまにさわりだけやりました. たぶんこれより先の内容は別の記事にします!!