C言語のOpenMPのreductionで配列要素をreduction変数に指定する

OpenMP4.5からC,C++言語のOpenMPのreductionで配列の要素をreduction変数に指定できるようになったみたいなのでそのメモ。ちなみにFortranはかなり前から配列をreduction変数に指定できていたらしいがC言語では最近まで使えなかった(最近と言ってもgcc5.xのOpenMP4.5からは使えてたようだが)。

使い方

以下の足し合わせを例にとる。

• 正しく足し合わせができればすべての要素に 500*2=1000 が入る。

  int *array;
  array = (int*) malloc (sizeof(int)*nsize);

  int a=2;

  for (int i=0; i<500; i++) {
    for (int ii=0; ii<nsize; ii++) {
      array[ii] += a;
    }
  }

• reduction 指示句に配列を指定するとコンパイルの際にエラーがでる。

# pragma omp parallel for schedule(auto) reduction(+:array)
  for (int i=0; i<500; i++) {
    for (int ii=0; ii<nsize; ii++) {
      array[ii] += a;
    }
  }

エラー: user defined reduction not found for ‘array’
   | #pragma omp parallel for reduction(+:array)

• reduction 指示句に配列だけでなくサイズを指定すればコンパイルが通り、想定通りの動作をする。
  • サイズ指定の方法がPythonなどの書き方っぽくてあまりC言語的ではないが
  • サイズは変数でいい
  • 下記の例はarray[0:nsize]でもいい。試していないが配列の一部を指定するためにこのような範囲指定をするようになっていると思う。

# pragma omp parallel for schedule(auto) reduction(+:array[:nsize])
  for (int i=0; i<500; i++) {
    for (int ii=0; ii<nsize; ii++) {
      array[ii] += a;
    }
  }

動作確認

check.c

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

//#define __CHECK__

int check_array(int *array, int nsize)
{
# ifdef __CHECK__
  for (int ii=0; ii<nsize; ii++) {
    printf("%d %d\n",ii,array[ii]);
  }
# endif

  int array0=array[0];
  int flag=0;
  for (int ii=0; ii<nsize; ii++) {
    if(array0 != array[ii]) flag++;
  }

  return flag;
}

void zero_set(int *array, int nsize)
{
# pragma omp parallel for
  for (int i=0; i<nsize; i++) {
    array[i] = 0.0;
  }
}

void reduction_array(int *array, int nsize)
{
  zero_set(array, nsize);

  int a = 2;

//#pragma omp parallel for  // <- incorrect
//#pragma omp parallel for reduction(+:array)      // <- compile error
# pragma omp parallel for schedule(auto) reduction(+:array[:nsize])  // <- correct
  for (int i=0; i<500; i++) {
    for (int ii=0; ii<nsize; ii++) {
      array[ii] += a;
    }
  }

  int flag;
  flag = check_array(array, nsize);

  printf("correct flag : %d\n",flag);
}

int main(int argc, char **argv)
{
  int nsize;
  nsize = 1024*1024;

  int *array;
  array = (int*) malloc (sizeof(int)*nsize);

  reduction_array(array, nsize);

  free(array);
}

• 上記のファイルをコンパイルし実行
  • flagの値が0になっていれば全ての配列要素に正しい値が入っている
  • 700% になっているが一応この例では8スレッドで動いている

$ gcc -fopenmp -O3 check.c
$ time ./a.out
correct flag : 0
./a.out  0.10s user 0.00s system 698% cpu 0.014 total

• コメントアウトを変えてreduction指示句をなしにすると結果は正しくない

# pragma omp parallel for 
//#pragma omp parallel for schedule(auto) reduction(+:array[:nsize])

$ time ./a.out
correct flag : 1048568
./a.out  0.17s user 0.00s system 693% cpu 0.024 total

スタックサイズ

reduction指示句の変数は各スレッドのスタック領域に確保されるため大きなサイズの配列ではスタック領域が足りなくなる。上記の例ではnsize = 1024*1024*2;にするとsizeof(int)*nsize=8MBとなり環境変数OMP_STACKSIZEのデフォルトの4MBや8MBより大きいかほぼ同じになり segmentation faultで実行できない。その場合は環境変数を変更し、実行すればよい。OMP_STACKSIZEはreduction配列サイズよりある程度余裕を持たせる必要がある。

$ export OMP_STACKSIZE=16M
or
$ OMP_STACKSIZE=16M ./a.out

これでもsegmentation faultが出る場合はシステムのlimitを確認してみる。

$ limit             
cputime         unlimited
filesize        unlimited
datasize        unlimited
stacksize       8MB
coredumpsize    unlimited
memoryuse       unlimited
maxproc         63772
...

ここのstacksizeが優先されるので最低でもOMP_STACKSIZEより大きくするかunlimitedにしておく。

$ limit stacksize 16M
$ limit stacksize unlimited

その他の方法

参考までにreduction変数に配列が指定できなかったときのやり方と比較。reduction_arrayAは上で紹介したもの、reduction_arrayBはヒープ領域にスレッド数分の配列を確保しあとからターゲットの配列に逐次的に足し合わせ、reduction_arrayCは各スレッド内でヒープ領域に配列を確保し最後にcritical構文で足し合わせ。この例ではarrayBとarrayCはスタック領域を食い潰すわけではないのでOMP_STACKSIZEはデフォルトから変更の必要はない。

array_reduction.c

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

//#define __CHECK__

int check_array(int *array, int nsize)
{
# ifdef __CHECK__
  for (int ii=0; ii<nsize; ii++) {
    printf("%d %d\n",ii,array[ii]);
  }
# endif

  int array0=array[0];
  int flag=0;
  for (int ii=0; ii<nsize; ii++) {
    if(array0 != array[ii]) flag++;
  }

  return flag;
}

void zero_set(int *array, int nsize)
{
# pragma omp parallel for
  for (int i=0; i<nsize; i++) {
    array[i] = 0.0;
  }
}


void reduction_arrayA(int *array, int nsize)
{
  zero_set(array, nsize);

  int a = 2;

  double start_time, end_time;
  start_time = omp_get_wtime();

# pragma omp parallel for schedule(auto) reduction(+:array[:nsize])  // <- correct
  for (int i=0; i<500; i++) {
    for (int ii=0; ii<nsize; ii++) {
      array[ii] += a;
    }
  }

  end_time = omp_get_wtime();

  int flag;
  flag = check_array(array, nsize);

  printf("arrayA : correct flag : %d\n",flag);
  printf("elapsed time : %g [sec]\n",end_time-start_time);
}


void reduction_arrayB(int *array, int nsize)
{
  zero_set(array, nsize);

  int a = 2;

  double start_time, end_time;
  start_time = omp_get_wtime();

  int num_threads = omp_get_max_threads();
  int *array_thread;
  array_thread = (int*) malloc (sizeof(int)*nsize*num_threads);

# pragma omp parallel for
  for (int i=0; i<nsize*num_threads; i++) array_thread[i] = 0.0;

# pragma omp parallel
  {
    int my_thread = omp_get_thread_num();

# pragma omp for schedule(auto)
    for (int i=0; i<500; i++) {

      for (int ii = 0; ii<nsize; ii++) {
        int idx = ii + my_thread*nsize;
        array_thread[idx] += a;
      }
    }
  } // omp parallel


# pragma omp parallel for schedule(auto)
  for (int ii=0; ii<nsize; ii++) {

    int aa=0;
    for(int it=0; it<num_threads; it++) {
      int idx = ii+it*nsize;
      aa += array_thread[idx];
    }
    array[ii] = aa;
  }

  free(array_thread);

  end_time = omp_get_wtime();

  int flag;
  flag = check_array(array, nsize);

  printf("arrayB : correct flag : %d\n",flag);
  printf("elapsed time : %g [sec]\n",end_time-start_time);
}


void reduction_arrayC(int *array, int nsize)
{
  zero_set(array, nsize);

  int a = 2;

  double start_time, end_time;
  start_time = omp_get_wtime();

# pragma omp parallel
  {
    int *array_thread;
    array_thread = (int*) malloc (sizeof(int)*nsize);

    for (int i=0; i<nsize; i++) array_thread[i] = 0.0;

# pragma omp for schedule(auto)
    for (int i=0; i<500; i++) {

      for (int ii = 0; ii<nsize; ii++) {
        array_thread[ii] += a;
      }
    }

# pragma omp critical
    for (int ii = 0; ii<nsize; ii++) {
      array[ii] += array_thread[ii];
    }

    free(array_thread);

  } // omp parallel

  end_time = omp_get_wtime();

  int flag;
  flag = check_array(array, nsize);

  printf("arrayC : correct flag : %d\n",flag);
  printf("elapsed time : %g [sec]\n",end_time-start_time);
}

int main(int argc, char **argv)
{
  int nsize;
  if(argc==2) nsize = atoi(argv[1]);
  else        nsize = 1024*1024*8;

  int *array;
  array = (int*) malloc (sizeof(int)*nsize);

  int num_threads = omp_get_max_threads();
  printf("num_threads = %d\n",num_threads);

  reduction_arrayA(array, nsize);
  reduction_arrayB(array, nsize);
  reduction_arrayC(array, nsize);

  free(array);
}

計測時間

# nsize=1024*1024*8`の場合
num_threads = 8
arrayA : correct flag : 0
elapsed time : 1.18341 [sec]
arrayB : correct flag : 0
elapsed time : 1.16044 [sec]
arrayC : correct flag : 0
elapsed time : 1.18492 [sec]

# nsize=1024*1024*128`の場合
num_threads = 8
arrayA : correct flag : 0
elapsed time : 19.0628 [sec]
arrayB : correct flag : 0
elapsed time : 18.7186 [sec]
arrayC : correct flag : 0
elapsed time : 19.2374 [sec]

サイズを変えたりして何回か計測してもだいたい上記のような結果になった。最も速いというわけではないが可読性の観点からすれば無理にarrayBのような書き方をしないでもいいと思う。arrayAとarrayCはスタック領域かヒープ領域に配列を取るかの違いで内部的にはcritical構文が使われているのかもしれない。

まとめ

• C,C++言語でもOpenMPのreduction変数に配列要素を指定できるようになった
• ヒープ領域を使う書き方よりかなり可読性は向上する
• limit stacksizeとOMP_STACKSIZEを適切に変更する必要がある

参考

• https://stackoverflow.com/questions/3775147/is-it-possible-to-do-a-reduction-on-an-array-with-openmp/40470790
• http://dvalts.io/optimisation/code/2016/11/06/OpenMP-array_reduction.html