用 OpenMP 實作平行化的 Quicksort

OpenMP 和 Pthread 不同，我們沒有辦法很精確地控制有多少的 Thread 正在執行。所以如果要實作 Quicksort，我們不能直接使用 Recursion。相反地，我們必須引入 Threading Pool 的概念，讓各個 Thread 主動去某個工作清單拿工作。概念程式碼如下：

void quicksort_impl(sort_value_t array[],
                    stack_t *s,
                    size_t *busy_threads) {
  int idle = 1;

  size_t begin = 0;
  size_t end = 0;

  /* We should loop forever until everyone is idle */
  while (1) {
    while (idle) {
      /* This thread is idling.
         Try to get a new task from the stack. */

#pragma omp critical
      if (!stack_empty(s)) {
        /* Get the pending sort range, start work now! */
        sort_range_t *range = stack_top(s);
        begin = range->begin;
        end = range->end;
        stack_pop(s);

        ++*busy_threads;
        idle = 0;
      }

      /* If everyone finished his work, then leave now ... */
      if (*busy_threads == 0) {
        return;
      }
    }

    size_t size = end - begin;

    if (size <= 50) {
      sort_value_t *seg = array + begin;

      /* Use backward insertion sort when the sort range is small */
      size_t i, j;
      for (i = 1; i < size; ++i) {
        sort_value_t cur = seg[i];
        /* Move the number backward if they are greater than cur */
        for (j = i; j > 0 && seg[j - 1] > cur; --j) {
          seg[j] = seg[j - 1];
        }
        seg[j] = cur;
      }

      /* Mark begin as end so that we can acquire next sort range */
      begin = end;

#pragma omp critical
      --*busy_threads;

      idle = 1;

      /* Finished our range continue to acquire next range */
      continue;
    }


    /* Partition our range */
    sort_value_t *seg = array + begin;

    size_t pivot_index = rand() % size;
    SWAP(seg[0], seg[pivot_index]);
    sort_value_t pivot = seg[0];

    size_t left = 1;
    size_t right = size - 1;

    while (left < right) {
      while (seg[left] <= pivot && left < right) { left++; }
      while (seg[right] >= pivot && left < right) { right--; }
      SWAP(seg[left], seg[right]);
    }

    if (seg[left] > pivot) {
      left--;
    }
    SWAP(seg[0], seg[left]);

    /* Push left subrange to stack */
#pragma omp critical
    {
      sort_range_t range;
      range.begin = begin;
      range.end = begin + left;
      stack_push(s, &range);
    }

    /* Continue to sort right subrange */
    begin += left + 1;
  }
}

void quicksort(sort_value_t array[], size_t size) {
  size_t busy_threads = 0;
  stack_t *s = stack_new();
  sort_range_t all = { 0 , size };
  stack_push(s, &all);

  /* Run quicksort_impl in parallel */
#pragma omp parallel
  quicksort_impl(array, s, &busy_threads);

  stack_delete(s);
}