LoginSignup
4
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@stanaka2

クイックソートをOpenMPでスレッド並列化

Last updated at Posted at 2019-09-07

C言語のクイックソートをOpenMPを用いたマルチスレッドで計算できるようにした。OpenMPの入れ子を使ったバージョンとタスクを使ったバージョンの二種類作った。結論としては両方共それなりにマルチスレッドの恩恵を得られるが、スレッドの使い方の性質上タスクバージョンのほうがいい。

下記のコードはOpenMP4.5で動作を確認しているがOpenMP3以上であれば使えるはず。

[2020.8.28] 修正 : Intel compilerでバグることがあったので関数宣言時の不要なconstを削除

説明

配列keyをソートする。

シングルスレッドクイックソート

いわゆる通常のクイックソート。マルチスレッドクイックソートで要素の長さがあるしきい値まで短くなったときに使用する。pivotの取り方は適当に選んだ3要素の中央値をmedian関数でセットしている。ここに出てきていない関数は下の方の性能測定コードの中を参照。

single_thread_qsort.c
static inline void qsort_partitioning(int *key, int *left, int *right,
				      int pivot)
{
  int ii, jj;
  ii = *left;
  jj = *right;

  // pivotより値が小さいか大きいかで要素を左か右に集める

  /* partitioning */
  while (1) {
    while (key[ii] < pivot)   ii++;
    while (pivot   < key[jj]) jj--;
    if (ii >= jj) break;

    const int tmpk = key[ii];

    key[ii] = key[jj];
    key[jj] = tmpk;

    ii++;
    jj--;
  }

  *left  = ii;
  *right = jj;
}

void single_thread_qsort(int *key, int left, int right)
{
  if (left < right) {
    int ii = left, jj = right;

    int pivot = median(key[ii],key[jj],key[(ii+jj)/2]);

    // pivotより値が小さいか大きいかで要素を左か右に集める
    qsort_partitioning(key, &ii, &jj, pivot);

    // 左側,右側で再帰呼出し
    single_thread_qsort(key, left, ii-1);
    single_thread_qsort(key, jj+1, right);
  }
}

OpenMP 入れ子版クイックソート

2分木的にスレッドが増えていき、最深部に到達するか分割された配列サイズがしきい値より小さくなった時点でシングルスレッドクイックソートに移行する。再帰呼び出しの中でスレッド生成をしていくので生成コストはかかるがせいぜい数段分しかスレッド生成しないので無視できると思っていい。それよりも2分木的にスレッドを増やす際になるべく正確な中央値を選ばないと2スレッドに分配する要素数に偏りができてしまい、設定したスレッドを使い切れない恐れがある。例えば最初に2分木を作る際に片方の要素数がものすごく少ない場合には以降に片側で新たなスレッドが生成されず結果的に半分のスレッドしか使えないということが考えられる。なのでシングルスレッドクイックソートに移行するまでは81点からおおよその中央値を求めるようにし、なるべく偏りをなくしている。経験的に3点中央値だと偏りが大きく理想的な並列はできないのでコア数、スレッド数に比べて全然速くならないことが多い。

parallel_nested_qsort.c

void parallel_nested_qsort(int* key, const int left, const int right)
{
  // 現在のnested flagが0か1かを確認。
  // 0ならnested flagを1にする。
  // すでに1の場合はこの関数の最後で0にしないためにこのようなことをしている。 
  int omp_nest_flag = omp_get_nested();
  if(omp_nest_flag==0) omp_set_nested(1);

  // 使えるスレッド数から入れ子の深さを設定。
  omp_set_max_active_levels( log2i(omp_get_max_threads()));

  // 入れ子版マルチスレッドクイックソート本体
  parallel_nested_qsort_internal(key, left, right);

  // 元のフラグ状態に戻す。
  if(omp_nest_flag==0) omp_set_nested(0);
}

void parallel_nested_qsort_internal(int* key, int left, int right)
{
  int length = right - left;
  int active_level = omp_get_active_level();
  int max_level    = omp_get_max_active_levels();

  // nestが最も深いか要素数がしきい値より小さい場合はそれ以上nestしないで1スレッドで3点中央値を使ったクイックソート。
  // 関数を抜ける。
  if( (active_level==max_level) || (length < PARALLEL_QSORT_THRESHOLD) ) {
    single_thread_qsort(key, left, right);
    return;
  }

  int ii = left;
  int jj = right;
  
  // ネストレベルを増やして行く段階のpivotは81点中央値を採用
  int pivot = median81(key, ii, jj);

  // 上記の分岐がelseなら要素を2分割していく
  qsort_partitioning(key,&ii,&jj, pivot);

  // 2分割した配列が左と右に分かれるのでそれぞれにスレッドを割り当てる。
  // 各スレッドでparallel_nested_qsort_internal()を再帰的に呼び出し。
  // ここに到達したスレッドは毎回2スレッド生成する。
# pragma omp parallel num_threads(2)
  {
# pragma omp sections nowait
    {

# pragma omp section
      parallel_nested_qsort_internal(key, left, jj);
# pragma omp section
      parallel_nested_qsort_internal(key, ii, right);

    } // sections
  } // parallel
}

OpenMP Task版クイックソート

各スレッドにタスクを割り当てていき、分割された配列サイズがしきい値より小さくなった時点でシングルスレッドクイックソートに移行する。スレッド生成は最初の一回。再帰呼出しの中でタスク生成をしていくので生成コストはかかるが大きくはない。一般的にはスレッド生成速度>タスク生成速度らしいのでタスク構文は再帰呼出しを行う計算ではよく使われている。タスク版は使えるスレッドをどんどん使ってくれるので入れ子版で起こるような要素数の偏りによって使えるスレッド数が発生することはない。そのため81点中央値を取る必要もなく3点中央値で十分(むしろ81点だとそのコストで遅くることが多い)。

parallel_task_qsort.c
void parallel_task_qsort(int* key, int left, int right)
{
// はじめに並列領域を生成
# pragma omp parallel
  {

// とりあえずsingle構文で1スレッドだけで関数を呼び出す
# pragma omp single nowait
    {
      parallel_task_qsort_internal(key, left, right);
    }
  }
}

void parallel_task_qsort_internal(int* key, int left, int right)
{

  int length = right - left;

  // 要素数がしきい値より小さい場合はそれ以上再帰呼出ししないで1スレッドでクイックソート。
  // 関数を抜ける。
  if ( length < PARALLEL_QSORT_THRESHOLD ) {
    single_thread_qsort(key, left, right);
    return;
  }

  int ii = left;
  int jj = right;

  int pivot = median(key[ii],key[jj],key[(ii+jj)/2]);

  // 上記の分岐がelseなら要素を2分割していく
  qsort_partitioning(key,&ii,&jj, pivot);

// 左側の要素でタスク生成
# pragma omp task
  parallel_task_qsort_internal(key, left, jj);

// 右側の要素でタスク生成
# pragma omp task
  parallel_task_qsort_internal(key, ii, right);
}

コード例

parallel_qsort.c
# include <stdio.h>
# include <stdlib.h>
# include <stdint.h>
# include <omp.h>

// シングルスレッドクイックソートに移行するしきい値
// しきい値の長さのソートよりスレッドやタスクの生成コストが上回ると速くならないので小さ過ぎは注意
// 当然大きすぎてもマルチスレッドの恩恵は得られないので1000くらいでいい
# define PARALLEL_QSORT_THRESHOLD (1000)

// sortした値が正しいか確認したい場合はコメントアウトを外す
//#define __CHECK_SORT_DATA__

int log2i(const int index)
{
  int tmp = index;
  int level = 0;
  while (tmp >>= 1) ++level;

  return level;
}


static inline int median(const int a, const int b, const int c)
{
  if(a>b) {
    if(b>c)      return b;
    else if(c>a) return a;
    else         return c;

  }else{
    if(a>c)      return a;
    else if(c>b) return b;
    else         return c;
  }
}

static inline int median9(int *a, const int left, const int right)
{
  int del = (int)((right-left)/9);
  return median(median( a[left],       a[left+del],   a[left+2*del] ),
		median( a[left+3*del], a[left+4*del], a[left+5*del] ),
		median( a[left+6*del], a[left+7*del], a[left+8*del] ));
}


static inline int median27(int *a, const int left, const int right)
{
  int del = (int)((right-left)/3);
  return median(median9( a, left, left+del ),
		median9( a, left+del+1, left+2*del ),
		median9( a, left+2*del+1, right ));
}


static inline int median81(int *a, const int left, const int right)
{
  int del = (int)((right-left)/3);
  return median(median27( a, left, left+del ),
		median27( a, left+del+1, left+2*del ),
		median27( a, left+2*del+1, right ));
}



static inline void qsort_partitioning(int *key, int *left, int *right,
				      int pivot)
{
  int ii, jj;
  ii = *left;
  jj = *right;

  /* partitioning */
  while (1) {
    while (key[ii] < pivot)   ii++;
    while (pivot   < key[jj]) jj--;
    if (ii >= jj) break;

    const int tmpk = key[ii];

    key[ii] = key[jj];
    key[jj] = tmpk;

    ii++;
    jj--;
  }

  *left  = ii;
  *right = jj;
}

void single_thread_qsort(int *key, int left, int right)
{
  if (left < right) {
    int ii = left, jj = right;

    int pivot = median(key[ii],key[jj],key[(ii+jj)/2]);
    qsort_partitioning(key, &ii, &jj, pivot);

    single_thread_qsort(key, left, ii-1);
    single_thread_qsort(key, jj+1, right);
  }
}


void parallel_nested_qsort_internal(int* key, int left, int right)
{
  int length = right - left;
  int active_level = omp_get_active_level();
  int max_level    = omp_get_max_active_levels();

  if( (active_level==max_level) || (length < PARALLEL_QSORT_THRESHOLD) ) {
    single_thread_qsort(key, left, right);
    return;
  }

  int ii = left;
  int jj = right;

  int pivot = median81(key, ii, jj);

  qsort_partitioning(key,&ii,&jj, pivot);

# pragma omp parallel num_threads(2)
  {
# pragma omp sections nowait
    {

# pragma omp section
      parallel_nested_qsort_internal(key, left, jj);
# pragma omp section
      parallel_nested_qsort_internal(key, ii, right);

    } // sections
  } // parallel
}

void parallel_nested_qsort(int* key, int left, int right)
{
  int omp_nest_flag = omp_get_nested();
  if(omp_nest_flag==0) omp_set_nested(1);

  omp_set_max_active_levels( log2i(omp_get_max_threads()));

  parallel_nested_qsort_internal(key, left, right);

  if(omp_nest_flag==0) omp_set_nested(0);
}


void parallel_task_qsort_internal(int* key, int left, int right)
{

  int length = right - left;

  if ( length < PARALLEL_QSORT_THRESHOLD ) {
    single_thread_qsort(key, left, right);
    return;
  }

  int ii = left;
  int jj = right;

  int pivot = median(key[ii],key[jj],key[(ii+jj)/2]);

  qsort_partitioning(key,&ii,&jj, pivot);

# pragma omp task
  parallel_task_qsort_internal(key, left, jj);

# pragma omp task
  parallel_task_qsort_internal(key, ii, right);

}

void parallel_task_qsort(int* key, int left, int right)
{
# pragma omp parallel
  {
# pragma omp single nowait
    {
      parallel_task_qsort_internal(key, left, right);
    }
  }
}



int asc_int(const void *a, const void *b)
{
    return *(int*)a - *(int*)b;
}

void stdlib_qsort(int *key, int nsize)
{
  qsort(key, nsize, sizeof(int), asc_int);
}



void set_data(int *key, int n)
{
  srand(0);

  for(int i=0;i<n;i++) {
    key[i] = rand();
  }

# ifdef __CHECK_SORT_DATA__
  printf("key     : %d %d %d %d %d\n",key[0],key[n/10],key[n/2],key[(int)(0.9*n)],key[n-1]);
# endif
}


void compare_data(int *key, int *sol_key, int n)
{
  int diff_key=0;

# pragma omp parallel for reduction(+:diff_key)
  for(int ix=0; ix<n; ix++) {
    if(key[ix] != sol_key[ix]) diff_key++;
  }

  printf("key     : %d %d %d %d %d\n",key[0],key[n/10],key[n/2],key[(int)(0.9*n)],key[n-1]);
  printf("sol_key : %d %d %d %d %d\n",sol_key[0],sol_key[n/10],sol_key[n/2],sol_key[(int)(0.9*n)],sol_key[n-1]);

  printf("diff %d\n", diff_key);
}


int main(int argc, char **argv)
{
  double start_time, end_time;

  printf("omp_get_max_threads, max_nested_level : %d, %d\n",omp_get_max_threads(),log2i(omp_get_max_threads()));

  int *key, *sol_key;

  int n;
  if(argc==2) n = atoi(argv[1]);
  else        n = 1<<20;

  key = (int*)malloc(sizeof(int)*n);
  sol_key = (int*)malloc(sizeof(int)*n);

  /* random data set */
  set_data(key, n);

  start_time = omp_get_wtime();

  stdlib_qsort(key, n);

  end_time = omp_get_wtime();

  printf("# stdlib qsort :: %14.6e [sec]\n",  end_time-start_time);

# ifdef __CHECK_SORT_DATA__
  for(int ix=0; ix<n; ix++) {
    sol_key[ix] = key[ix];
  }
# endif

  /* random data reset with same seed */
  set_data(key, n);

  start_time = omp_get_wtime();

  single_thread_qsort(key, 0, n-1);

  end_time = omp_get_wtime();

  printf("# single thread qsort :: %14.6e [sec]\n",  end_time-start_time);

# ifdef __CHECK_SORT_DATA__
  compare_data(key, sol_key, n);
# endif

  /* random data reset with same seed */
  set_data(key, n);

  start_time = omp_get_wtime();

  parallel_nested_qsort(key, 0, n-1);

  end_time = omp_get_wtime();

  printf("# parallel_nested_qsort :: %14.6e [sec]\n",  end_time-start_time);

# ifdef __CHECK_SORT_DATA__
  compare_data(key, sol_key, n);
# endif

  /* random data reset with same seed */
  set_data(key, n);

  start_time = omp_get_wtime();

  parallel_task_qsort(key, 0, n-1);

  end_time = omp_get_wtime();

  printf("# parallel_task_qsort :: %14.6e [sec]\n", end_time-start_time);

# ifdef __CHECK_SORT_DATA__
  compare_data(key, sol_key, n);
# endif

  free(key);
  free(sol_key);

  return EXIT_SUCCESS;
}

$ gcc -std=c11 -O3 -fopenmp parallel_qsort.c
$ ./a.out # デフォルトのスレッド数で10^6要素くらい
$ OMP_NUM_THREADS=16 ./a.out $((10**8))  # スレッド数と要素数を指定して実行

性能評価

上記のコードを若干修正したもので性能評価。初期にそれぞれ同じ乱数seedを与えて計算。乱数seedを10個与えてその平均値を取っている。計算環境は4コア/8スレッドのハイパースレッディングを使ったマシン。nested qsortとtask qsortはOpenMP8スレッド数計算。

計測時間の比較

stdlibのqsort、シングルスレッドqsort、入れ子を使ったqsort(8スレッド)、taskを使ったqsort(8スレッド)の比較。単位は秒。

size stdlib single thread nested task
1,000 2.30e-4 3.57e-5 3.63e-5 3.36e-5
10,000 9.19e-4 4.47e-4 1.08e-3 1.76e-4
100,000 8.41e-3 5.12e-3 2.05e-3 1.93e-3
1,000,000 1.01e-1 6.00e-2 1.90e-2 1.32e-2
10,000,000 1.19e 6.89e-1 1.92e-1 1.50e-1
100,000,000 1.37e+1 7.73 2.08 1.59
  • 参考としてstdlibのものを載せたが様々な型に対応させるようにいちいち関数を読み込んでいるため非常に遅い
  • task版のnested版より基本的に1.3倍くらい高速
    • nested版はシングルスレッドに移行した後に処理が終了したスレッドは仕事をしないから
    • task版は常に全スレッドが働いている

スレッド数1を基準とした際の加速率

size single thread nested task
1,000 1 0.98 1.18
10,000 1 0.41 2.53
100,000 1 2.50 2.68
1,000,000 1 3.15 4.55
10,000,000 1 3.59 4.59
100,000,000 1 3.72 4.86
  • task版の方が高速
  • 物理コアが4個なので4倍が限界だと思ったがtask版は(何回か試しても)コア数以上に高速にソートできている
    • OSレベルで8スレッドで動いているので1スレッドクイックソートが若干遅めなのかもしれない
    • ハイパースレッディングをオフにしたら変わるかもしれない

まとめ

  • OpenMPを使いクイックソートをスレッド並列化した
  • OpenMPのnestedとtaskを使った二つの手法を紹介
  • nested版は常に配列の中央値を選べれば理論的にはスレッド数倍速くなるが現実的には使われないスレッドが足を引っ張る
  • スレッドの使い方の性質上task版を使ったほうが高速

参考

4
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
4
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?