OpenMPをうずまきナルトにたとえて話してみる...

はじめに

OpenMP は非常にカンタンに手が出せてしまいます!!
それは GCC とかに抱き合わせでついているので、環境構築コストが低いです。逆に MPI とかは、自分でちゃんと環境を導入する必要があり、動かすまでのコストがかなり掛かります!!

ただし、OpenMP も楽か?といわれるとこれはこれで結構めんどくさい。カンタンにかけてしまうから、ド素人が何も考えずに書いてしまうと遅くなる&結果がおかしくなります...

それを防ぐためにも、ある程度(atomicとか以外)はその機能について知っておかないとまずい...

概説

実は、というかふつーに OpenMP はスレッド並列です!
スレッドを立ててスレッドに仕事を割り振ることで計算を並列化しています。逆に、MPI とかはプロセス並列、といって、プロセスを立てることで並列化を実現しています。

この記事内での対応について

この記事内ではスレッドをうずまきナルトに例えて説明します!
そして、OpenMP でスレッドを立てることは影分身の術に対応します!!

初歩編

OpenMP を使うには...

基本的にこちらでは GCC を基本とします!!
Intel さんの立派なコンパイラとか世の中にはいろいろあるけど、このほうがベタでしょう。

その壱 : 単純に影分身の術!!

#include <stdio.h>

int main() {
  #pragma omp parallel
  {
    /**
     ** このブロック内がすべてのスレッドで処理される!
    **/
    printf("影分身の術!\n");
  }
  return 0;
}

これだけで並列化できてしまいます!!
詳しくいうと#pragma omp parallelによってスレッドが生成され、ブロック内の処理を各スレッドが実行します!

ヘッダに OpenMP について何もないけど大丈夫なの?

という人がいると思いますが大丈夫です! 次節で説明します!

実行してみよう!!

さぁ、術の印?はわかったので、実際に術を発動してみましょう!!

gcc <コード> -fopenmp

でコンパイルできます!!

そして,

./kagebunshin

で実行すると...

俺はうずまきナルトだってばよ!
俺はうずまきナルトだってばよ!
俺はうずまきナルトだってばよ!
俺はうずまきナルトだってばよ!
俺はうずまきナルトだってばよ!
俺はうずまきナルトだってばよ!
俺はうずまきナルトだってばよ!
俺はうずまきナルトだってばよ!

と表示されるはずです!!
これは、持っている CPU の論理スレッド数(例えば MacBook Pro 15 inch late2016 だと、4 コア 8 スレッド)だけスレッドが生成されます!

その弐 : スレッド数を表示しよう!

これだけだとどの影分身が言ってるかわかりません!!
ってことで、影分身の番号をつけてあげましょう!!

#include <stdio.h>
#include <omp.h>

int main() {
  #pragma omp parallel
  {
    /**
     ** このブロック内がすべてのスレッドで処理される!
    **/
    printf("俺は%d番目のうずまきナルトだってばよ!\n",omp_get_thread_num());
  }
  return 0;
}

これを実行すると

俺は1番目のうずまきナルトだってばよ!
俺は2番目のうずまきナルトだってばよ!
俺は5番目のうずまきナルトだってばよ!
俺は3番目のうずまきナルトだってばよ!
俺は0番目のうずまきナルトだってばよ!
俺は6番目のうずまきナルトだってばよ!
俺は4番目のうずまきナルトだってばよ!
俺は7番目のうずまきナルトだってばよ!

ほら、簡単でしょ?
ここで詳しく説明すると、omp ヘッダをインクルードして実行します!

その参 : うずまきナルト連打?

ここまでやってきたことは、あんま意味がありません!
ただ影分身をだして、そいつらに点呼させただけ。これじゃ並列化する意味がありませんよね?
ということで、並列化の意味があることをさせてみましょう!

#include <stdio.h>
#define N 20
int main() {
  int a[N] = {0}, b[N] = {0}, c[N] = {0};
  int i = 0;
  //ここから並列化だってばよ!
  #pragma omp parallel
  {
    #pragma omp for
    for(i = 0; i < N; ++i){
      a[i]=b[i] = i;
      c[i] = a[i]*b[i];
    }
  }

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

結果はこうなります!

c[0] = 0
c[1] = 1
c[2] = 4
c[3] = 9
c[4] = 16
c[5] = 25
c[6] = 36
c[7] = 49
c[8] = 64
c[9] = 81
c[10] = 100
c[11] = 121
c[12] = 144
c[13] = 169
c[14] = 196
c[15] = 225
c[16] = 256
c[17] = 289
c[18] = 324
c[19] = 361

いままでは#pragma omp parallelだけだったけど、そのなかで#paragma omp forっていう謎の忍術があります!
この#paragma omp forが大切です!
なにやっているか、というと、

まず#pragma omp parallelでスレッドを生成します。そして、#paragma omp forと書かれている部分を各スレッドに分割します。そうすることで、それぞれのスレッドが別々の処理をすることになります。

に尽きます。

たとえば、2 コア 4 スレッドのパソコンでやると、

for(i = 0; i < 5; i++){
  a[i]=b[i] = i;
  c[i] = a[i]*b[i];
}
for(i = 5; i < 10; i++){
  a[i]=b[i] = i;
  c[i] = a[i]*b[i];
}
for(i = 10; i < 15; i++){
  a[i]=b[i] = i;
  c[i] = a[i]*b[i];
}
for(i = 15; i < 20; i++){
  a[i]=b[i] = i;
  c[i] = a[i]*b[i];
}

といったように、for文が 4 分割され、それぞれのスレッドで処理されます!
これはもっと楽にかけて、

#include <stdio.h>
#define N 20
int main() {
  int a[N] = {0}, b[N] = {0}, c[N] = {0};
  int i = 0;
  //ここから並列化だってばよ!
  #pragma omp parallel for
  for(i = 0; i < N; ++i){
    a[i]=b[i] = i;
    c[i] = a[i] * b[i];
  }

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

その四 : 螺旋丸を投げてみよう!!

ここが最初の難関かもしれないです。まず、うずまきナルト連打を真似て書いてみましょう!
やることは簡単!
単に$1,2,\ldots,n$の総和を計算するだけ!

#include <stdio.h>
#include <stdlib.h>

int main(int argc, char *argv[]) {
  int i;
  if(argc != 2) {
    printf("Nの値を入れないと計算しようがないってばよ...\n");
    exit(-1);
  }
  int N = atoi(argv[1]);
  long long int sum = 0;
  #pragma omp parallel for
  for(i = 0; i < N; ++i){
    sum += i;
  }
  printf("螺旋丸!! (威力%lld)\n",sum);
}

さぁ、これを$N=100000000$で試してみよう!
まずは、逐次で(-fopenmpのオプションを付けずに)やってみよう!

螺旋丸!! (威力4999999950000000)

わお、ひとりでも螺旋丸打てるやん...

じゃぁ、つぎに影分身使って螺旋丸打ってみよう!

螺旋丸!! (威力705116177563340)

あれ、桁が一つもちいさい...

おかしい理由

OpenMPはメモリ共有での並列化を行っています。つまり、変数の当てられているメモリは基本的にみんなの財産です!
今回おかしくなったのは、単純に、sumに当てられている変数にみんなでアクセスしちゃったから!
単純にいうと、数人の影分身が一斉にメモリの書き換えを行っても、競合が生じますよね...
これを防ぐためにOpenMPではReductionという機能があります!

Reduction

具体的には

#pragma omp parallel for reduction(演算子:変数)

と書きます!
こうすることで、競合した際のエラーがなくせます!
これを使ってみましょう!!

#include <stdio.h>
#include <stdlib.h>

int main(int argc, char *argv[]) {
  int i;
  if(argc != 2) {
    printf("Nの値を入れないと計算しようがないってばよ...\n");
    exit(-1);
  }
  int N = atoi(argv[1]);
  long long int sum = 0;
  #pragma omp parallel for reduction(+:sum)
  for(i = 0; i < N; ++i){
    sum += i;
  }
  printf("螺旋丸!! (威力%lld)\n",sum);
}

-fopenmpをつけてコンパイル&実行してみましょう!

螺旋丸!! (威力4999999950000000)

うまくいきました!