OpenMPでクイックソート

なんかクイックソートブームが自分の中で起きてました。
早くするにはどうするんだろうと考えてて、アルゴリズム的にばしっと分割するんだからマルチスレッド行けるんじゃねと思ってOpenMP化してみた。

qsort.cpp

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

template<typename T>
class Rand{
  T x = 1;
public:
  T next() {
    x ^= x >> 12;
    x ^= x << 25;
    x ^= x >> 27;
    return x * 2685821657736338717;
  }
};

template<typename T>
void part(T a[], int& i, int& j) {
  const T pivot = a[i];
  while (1) {
    while (a[i] < pivot) i++;
    while (pivot < a[j]) j--;
    if (i >= j) break;
    const T tmp = a[i]; a[i] = a[j]; a[j] = tmp;
    i++; j--;
  }
}

template<typename T>
void _qsort(T a[], const int left, const int right) {
  if (left < right) {
    int i = left, j = right;
    part(a, i, j);
    _qsort(a, left, i - 1);
    _qsort(a, j + 1, right);
  }
}

template<typename T>
void qsort(T a[], const int left1, const int right1) {
  if (left1 >= right1) return;

  int i1 = left1, j1 = right1;
  part(a, i1, j1);
  const int left2 = left1, right2 = i1-1;
  const int left3 = j1+1, right3 = right1;

  int i2 = left2, j2 = right2;
  int i3 = left3, j3 = right3;
  part(a, i2, j2);
  part(a, i3, j3);

#pragma omp parallel shared(a)
  {
#pragma omp sections
    {
#pragma omp section
      _qsort(a, left2, i2 - 1);
#pragma omp section
      _qsort(a, j2 + 1, right2);。
#pragma omp section
      _qsort(a, left3, i3 - 1);
#pragma omp section
      _qsort(a, j3 + 1, right3);
    }
  }
}

int main(int argc, char** argv) {
  const int size = 10000000;

  auto rand = new Rand<int>();

  auto a = new int[size];
  for (int i = 0; i < size; i++) {
    a[i] = rand->next();
  }

  double st = omp_get_wtime();
  qsort(a, 0, size-1);
  double ut = omp_get_wtime() - st;
  printf("Used time         = %2.6lf (s)\n", ut);

  delete a;
  return 0;
}

やり方は一発目に呼ばれるqsort関数で配列を4つに分割して、それをスレッドにして普通のクイックソートである_qsort関数を実行する感じ。普通のクイックソートにラッパーかますイメージ。

qsort関数がごちゃってるのは本来、再帰な処理を一段階、無理やりインライン化してるから。まず全体を分割(part関数)して、さらにそれぞれを分割(part関数)することで4分割にして、_qsortに渡す、っていうラッパーがこの関数。

肝心の速度は割と早くなりました。