OpenMPの基本的な使い方

サーバサイドの大規模計算・組み込み系・ゲームプログラミングなどの分野ではしばしば高速化を求められることがある。
その中でも、スレッド並列に関する高速化はマルチコア化の波で近年特に求められている。
この記事ではスレッド並列にOpenMPを利用した高速化手法についての基本的な部分を紹介していく。
*言葉の解釈によっては適切でない部分や、OpenMPの正しい使い方から外れた部分があるかもしれません。
*2017.02.16 記事の一部に書き途中・適切ではない記載があったので更新しました。

# スレッド並列化とは？
スレッド並列化とは、近年のPCによく使われているマルチコアのCPU上で行われる並列化のことである。CPUによっても変わるが、主流のIntel社製Core iシリーズのCPUなどでは1コアで2スレッドまでのスレッドを持つことができる。

例に、2016年の新モデルMacBookProの13インチモデルの標準CPU (Intel® Core™ i5-6360U Processor)を例にする。
このCPUは**クロック周波数が2.0GHz/1コア**で、2コア・4スレッドの環境である。
理論的には、逐次プログラム(並列化されていない)よりも最大で4倍の性能向上を期待できる。 [^1] 理論値上は2倍の性能が期待できる。

[^1]: ここでいう4倍の性能向上というのは厳密には間違えており、パイプラインが十分に埋まっていなかったことを期待してスレッドを1コアに複数立てることでパイプラインを埋めて、2倍の性能向上が最大でできるということである。

# スレッド並列のプログラミング
スレッド並列なプログラムを記載する場合はpthreadを使う手法が昔から存在している。
しかし、排他ロックやスレッドの指定など、並列化の手続きが非常に複雑でなかなか手を出したくない領域であった。
また、他のスレッド並列化の手法もあったが、仕様もバラバラでプログラマーの立ち位置からは使いにくいものであった。

# OpenMPによる並列化
OpenMPはスレッド並列なプログラムをディレクティブを挿入するだけで可能にしてくれる。
スレッド間のアクセス競合や依存関係などに気をつけて利用する必要があるものの、簡単にスレッド並列なプログラムが実装できる。
今では、多くのOSやプロセッサで利用することができる。
最近では、タスク並列を利用したプログラムにも対応してきている。

# OpenMPのコンパイルオプション (C言語の場合)
iccコンパイラ　　-openmp
gccコンパイラ　　-fopenmp

# OpenMPの使い方
OpenMPがよく利用されるのは、配列の加算や乗算である。
以下に簡単な例を示す。

```omp_samp.c
#include<stdio.h>
#include<omp.h>

int main(){
  int a[1000];
  int b[1000];
  int c[1000];
  int i;

  #pragma omp parallel for
  for(i=0;i<1000;i++){
    a[i] = i;
    b[i] = 1;
    c[i] = a[i] + b[i];
  }

  for(i=0;i<1000;i++){
    printf("%d\n",c[i]);
  }

  return 0;
}

```
このように記載する。
実に簡単で、上記の例は```c[i] = a[i] + b[i];```という演算を1000回行っている。
for文の前に```#pragma omp parallel for```というディレクティブをつけると``Core™ i5-6360U``ならば、4スレッド並列になるので、**演算250回/1スレッド**で演算を行ってくれる。各スレッドが並列に処理を行うため、1スレッドの担当する演算量は全体の1/4になる。

# ワークロードを手動で分散する方法
多くの並列化では、for文の演算を並列に計算したいという場合であるが、稀にfor文単位で並列化を行うことがとても厄介なパターンがある。

例えば、スレッド数がとても多い時などである。
今年から東京大学が運用しているReedBush-Uでは、Intel Xeon E5-2695v4が内部のプロセッサとして使われているが、18コア36スレッドである。
また、インテルのコプロセッサ(Xeon Phi Processor Knights Corner)では約60コア240スレッドもの高並列環境になってしまう。

このような場合にスレッド番号を指定してそれぞれのスレッドで別の配列を処理してあげることが可能である。
例のプログラムを以下に示す。

```no_Workload.c
#include<stdio.h>
#include<omp.h>
#define N 1000

void add(int *a, int *b, int *c);
void mul(int *a, int *b, int *d);
void add_add(int *a, int *b, int *e);
void mul_mul(int *a, int *b, int *f);

int main(){
  int a[N];
  int b[N];
  int c[N];
  int d[N];
  int e[N];
  int f[N];
  int i;

  #pragma omp parallel for
  for(i=0; i < N; i++){
    a[i] = i;
    b[i] = 1;
  }

  add(a,b,c);
  mul(a,b,d);
  add_add(a,b,e);
  mul_mul(a,b,f);


  for(i=0; i < N; i++){
    printf("c : %d\n",c[i]);
    printf("d : %d\n",d[i]);
    printf("e : %d\n",e[i]);
    printf("f : %d\n",f[i]);
  }

  return 0;
}

void add(int *a, int *b, int *c){
  int i;
  #pragma omp parallel for
  for(i=0; i < N; i++){
    c[i] = a[i] + b[i];
  }
}

void mul(int *a, int *b, int *d){
  int i;
  #pragma omp parallel for
  for(i=0; i < N; i++){
    d[i] = a[i] * b[i];
  }
}

void add_add(int *a, int *b, int *e){
  int i;
  #pragma omp parallel for
  for(i=0; i < N; i++){
    e[i] = a[i] + b[i] + b[i];
  }
}

void mul_mul(int *a, int *b, int *f){
  int i;
  #pragma omp parallel for
  for(i=0; i < N; i++){
    f[i] = a[i] * b[i] * b[i];
  }
}

```
上記のような各関数で複雑な処理のコードを作成した場合[^2]に、スレッド数が多いと十分にスレッド並列にできない状態で並列実行が終わってしまうことが多い。
なので、できるだけスレッドを細かく指定して割りあてる方が、実行効率が上がる可能性が高い。
具体的には以下のようなプログラムにおいてワークロードの手動分散が有効である。

[^2]: 今回の``no_Workload.c``はそんなに複雑ではないので関数をまとめればすぐにfor文でワークシェア構文を利用した並列化ができるが、もっと関数内で特殊な処理をしていた場合などは、依存関係などもありなかなか簡単に書き換えることができない。そのような場合にワークロードを手動で分散するのは有効的である。

```Workload.c
#include<stdio.h>
#include<omp.h>
#define N 1000

void add(int *a, int *b, int *c, int min, int max);
void mul(int *a, int *b, int *d, int min, int max);
void add_add(int *a, int *b, int *e, int min, int max);
void mul_mul(int *a, int *b, int *f, int min, int max);

int main(){
  int a[N];
  int b[N];
  int c[N];
  int d[N];
  int e[N];
  int f[N];
  int i;

  #pragma omp parallel for
  for(i=0; i < N; i++){
    a[i] = i;
    b[i] = 1;
  }

  #pragma omp parallel
  {
    if(omp_get_thread_num() == 0){
      add(a,b,c,0,250);
      mul(a,b,d,0,250);
      add_add(a,b,e,0,250);
      mul_mul(a,b,f,0,250);
    } else if(omp_get_thread_num() == 1){
      add(a,b,c,250,500);
      mul(a,b,d,250,500);
      add_add(a,b,e,250,500);
      mul_mul(a,b,f,250,500);
    } else if(omp_get_thread_num() == 2){
      add(a,b,c,500,750);
      mul(a,b,d,500,750);
      add_add(a,b,e,500,750);
      mul_mul(a,b,f,500,750);
    } else if(omp_get_thread_num() == 3){
      add(a,b,c,750,1000);
      mul(a,b,d,750,1000);
      add_add(a,b,e,750,1000);
      mul_mul(a,b,f,750,1000);
    }
  }


  for(i=0; i < N; i++){
    printf("c : %d\n",c[i]);
    printf("d : %d\n",d[i]);
    printf("e : %d\n",e[i]);
    printf("f : %d\n",f[i]);
  }

  return 0;
}

void add(int *a, int *b, int *c, int min, int max){
  int i;
  for(i=min; i < max; i++){
    c[i] = a[i] + b[i];
  }
}

void mul(int *a, int *b, int *d, int min, int max){
  int i;
  for(i=min; i < max; i++){
    d[i] = a[i] * b[i];
  }
}

void add_add(int *a, int *b, int *e, int min,int max){
  int i;
  for(i=min; i < max; i++){
    e[i] = a[i] + b[i] + b[i];
  }
}

void mul_mul(int *a, int *b, int *f, int min, int max){
  int i;
  for(i=min; i < max; i++){
    f[i] = a[i] * b[i] * b[i];
  }
}

```
このように分岐を行うことでそれぞれの命令を各スレッドに割り当てて計算させる。
しかし、このコードはプログラムの挙動をif文で分岐させているため、あまりスマートではない。
条件の分岐をそれぞれのスレッドIDで制御すた方がより簡潔にプログラムを記述できる。
以下に条件分岐部分を改善したコードを示す。関数部分は冗長になるので省略する。


```Workload_Improvement.c
#include<stdio.h>
#include<omp.h>
#define N 1000

void add(int *a, int *b, int *c, int min, int max);
void mul(int *a, int *b, int *d, int min, int max);
void add_add(int *a, int *b, int *e, int min, int max);
void mul_mul(int *a, int *b, int *f, int min, int max);

int main(){
  int a[N];
  int b[N];
  int c[N];
  int d[N];
  int e[N];
  int f[N];
  int i;

  #pragma omp parallel for
  for(i=0; i < N; i++){
    a[i] = i;
    b[i] = 1;
  }

  #pragma omp parallel
  {
    int thread_id = omp_get_thread_num();
    int num = N/4;
    int thread_part = num * thread_id;

    add(a,b,c,thread_part,thread_part + num);
    mul(a,b,d,thread_part,thread_part + num);
    add_add(a,b,e,thread_part,thread_part + num);
    mul_mul(a,b,f,thread_part,thread_part + num);
  }


  for(i=0; i < N; i++){
    printf("c : %d\n",c[i]);
    printf("d : %d\n",d[i]);
    printf("e : %d\n",e[i]);
    printf("f : %d\n",f[i]);
  }

  return 0;
}

```



# task並列化
タスク並列なコードに関しては、まだ具体的にコードを書いたことがないので的確なことは記載できないが、再帰構造を利用してタスク並列にプログラムを並列化する手法である。
コードを実際に書いて理解したら追記する予定である。

# sectionについて
各処理を別々に並列化するプログラムを記載する時に利用する。
~~あまり利用する機会がないので省略・・・~~
一見すると”ワークロードを手動で分散する方法”とあまり変わらないと感じるかもしれない。
しかし、実際の関係で見ると「ワークロードを手動で分散 > section」である。
なぜなら、sectionを利用すると１セクションで１スレッドしか割り当てられないからである。
では、どのような時に利用できるのか？
sectionを利用する場面は独立した処理を少ないスレッドで並列化したい時に特に利用価値が高いのではないかと感じた・・・(もっと正しい使い方があれば、コメントいただけると嬉しいです。)
sectionを利用する場合は、```#pragma omp parallel sections```というディレクティブと```#pragma omp parallel section```というディレクティブをセットで利用する。以下にsection分割で「Workload.c」と同じ処理を行なった場合のプログラムを示す。関数部分は冗長になるので省略する。

```Section.c
#include<stdio.h>
#include<omp.h>
#define N 1000

void add(int *a, int *b, int *c, int min, int max);
void mul(int *a, int *b, int *d, int min, int max);
void add_add(int *a, int *b, int *e, int min, int max);
void mul_mul(int *a, int *b, int *f, int min, int max);

int main(){
  int a[N];
  int b[N];
  int c[N];
  int d[N];
  int e[N];
  int f[N];
  int i;

  #pragma omp parallel for
  for(i=0; i < N; i++){
    a[i] = i;
    b[i] = 1;
  }

  #pragma omp parallel sections
  {
    #pragma omp section
    {
      add(a,b,c,0,250);
      mul(a,b,d,0,250);
      add_add(a,b,e,0,250);
      mul_mul(a,b,f,0,250);
    }
    #pragma omp section
    {
      add(a,b,c,250,500);
      mul(a,b,d,250,500);
      add_add(a,b,e,250,500);
      mul_mul(a,b,f,250,500);
    }
    #pragma omp section
    {
      add(a,b,c,500,750);
      mul(a,b,d,500,750);
      add_add(a,b,e,500,750);
      mul_mul(a,b,f,500,750);
    }
    #pragma omp section
    {
      add(a,b,c,750,1000);
      mul(a,b,d,750,1000);
      add_add(a,b,e,750,1000);
      mul_mul(a,b,f,750,1000);
    }
  }


  for(i=0; i < N; i++){
    printf("c : %d\n",c[i]);
    printf("d : %d\n",d[i]);
    printf("e : %d\n",e[i]);
    printf("f : %d\n",f[i]);
  }

  return 0;
}

```
このように実行する。

# まとめ
あまりまとまりなく書いてしまったが、一般的にはfor文の前に``#pragma omp parallel for``とディレクティブを追加するだけで良い。

個々の処理に応じて並列化したい場合は、``#pragma omp parallel``でスレッド分割を行った後に条件によってスレッドIDを指定すれば良い。


# 参考
・http://kaworu.jpn.org/c/pthread
・[Intel® Core™ i5-6360U Processor](http://ark.intel.com/ja/products/91156/Intel-Core-i5-6360U-Processor-4M-Cache-up-to-3_10-GHz)
・[Wikipedia-マルチコア](https://ja.wikipedia.org/wiki/マルチコア)
・[Wikipedia-OpenMP](https://ja.wikipedia.org/wiki/OpenMP)
・[東大のスパコン-ReedBush-U](http://www.cc.u-tokyo.ac.jp/system/reedbush/reedbush_intro.html)