C OpenMP parallel quickSort

Question

Once again I'm stuck when using openMP in C++. This time I'm trying to implement a parallel quicksort.

Code:

#include <iostream>
#include <vector>
#include <stack>
#include <utility>
#include <omp.h>
#include <stdio.h>

#define SWITCH_LIMIT 1000

using namespace std;

template <typename T>
void insertionSort(std::vector<T> &v, int q, int r)
{
    int key, i;
    for(int j = q + 1; j <= r; ++j)
    {
        key = v[j];
        i = j - 1;
        while( i >= q && v[i] > key )
        {
            v[i+1] = v[i];
            --i;
        }
        v[i+1] = key;
    }
}

stack<pair<int,int> > s;

template <typename T>
void qs(vector<T> &v, int q, int r)
{
    T pivot;
    int i = q - 1, j = r;
    //switch to insertion sort for small data
    if(r - q < SWITCH_LIMIT) 
    {
        insertionSort(v, q, r);
        return;
    }

    pivot = v[r];
    while(true)
    {
        while(v[++i] < pivot);
        while(v[--j] > pivot);
        if(i >= j) break;
        std::swap(v[i], v[j]); 
    }
    std::swap(v[i], v[r]);

    #pragma omp critical
    {
        s.push(make_pair(q, i - 1));
        s.push(make_pair(i + 1, r));        
    }
}

int main()
{
    int n, x;
    int numThreads = 4, numBusyThreads = 0;
    bool *idle = new bool[numThreads];
    for(int i = 0; i < numThreads; ++i)
        idle[i] = true;
    pair<int, int> p;
    vector<int> v;
    cin >> n;
    for(int i = 0; i < n; ++i)
    {
        cin >> x;
        v.push_back(x);
    }
    cout << v.size() << endl;
    s.push(make_pair(0, v.size()));

    #pragma omp parallel shared(s, v, idle, numThreads, numBusyThreads, p) 
    {
        bool done = false;
        while(!done) 
        {
            int id = omp_get_thread_num();
            #pragma omp critical
            {
                if(s.empty() == false && numBusyThreads < numThreads) 
                {
                    ++numBusyThreads;
                    //the current thread is not idle anymore
                    //it will get the interval [q, r] from stack
                    //and run qs on it
                    idle[id] = false;
                    p = s.top();                    
                    s.pop();
                }
                if(numBusyThreads == 0)
                {
                    done = true;
                }
            }
            if(idle[id] == false)
            {

                qs(v, p.first, p.second);
                idle[id] = true;
                #pragma omp critical 
                --numBusyThreads;
            }

        }
    }
    return 0;
}

Algorithm:

To use openMP for a recursive function I used a stack to keep track of the next intervals on which the qs function should run. I manually add the 1st interval [0, size] and then let the threads get to work when a new interval is added in the stack.

The problem:

The program ends too early, not sorting the array after creating the 1st set of intervals ([q, i - 1], [i+1, r] if you look on the code. My guess is that the threads which get the work, considers the local variables of the quicksort function(qs in the code) shared by default, so they mess them up and add no interval in the stack.

How I compile:

g++ -o qs qs.cc -Wall -fopenmp

How I run:

./qs < in_100000 > out_100000

where in_100000 is a file containing 100000 on the 1st line followed by 100k intergers on the next line separated by spaces.

I am using gcc 4.5.2 on linux

Thank you for your help,

Dan

Answer 1

I didn't actually run your code, but I see an immediate mistake on p, which should be private not shared. The parallel invocation of qs: qs(v, p.first, p.second); will have races on p, resulting in unpredictable behavior. The local variables at qs should be okay because all threads have their own stack. However, the overall approach is good. You're on the right track.

---

Here are my general comments for the implementation of parallel quicksort. Quicksort itself is embarrassingly parallel, which means no synchronization is needed. The recursive calls of qs on a partitioned array is embarrassingly parallel.

However, the parallelism is exposed in a recursive form. If you simply use the nested parallelism in OpenMP, you will end up having thousand threads in a second. No speedup will be gained. So, mostly you need to turn the recursive algorithm into an interative one. Then, you need to implement a sort of work-queue. This is your approach. And, it's not easy.

For your approach, there is a good benchmark: OmpSCR. You can download at http://sourceforge.net/projects/ompscr/

In the benchmark, there are several versions of OpenMP-based quicksort. Most of them are similar to yours. However, to increase parallelism, one must minimize the contention on a global queue (in your code, it's s). So, there could be a couple of optimizations such as having local queues. Although the algorithm itself is purely parallel, the implementation may require synchronization artifacts. And, most of all, it's very hard to gain speedups.

---

However, you still directly use recursive parallelism in OpenMP in two ways: (1) Throttling the total number of the threads, and (2) Using OpenMP 3.0's task.

Here is pseudo code for the first approach (This is only based on OmpSCR's benchmark):

void qsort_omp_recursive(int* begin, int* end)
{
  if (begin != end) {
    // Partition ...

    // Throttling
    if (...)  {
      qsort_omp_recursive(begin, middle);
      qsort_omp_recursive(++middle, ++end);
    } else {

#pragma omp parallel sections nowait
      {
#pragma omp section
        qsort_omp_recursive(begin, middle);
#pragma omp section
        qsort_omp_recursive(++middle, ++end);
      }
    }
  }
}

In order to run this code, you need to call omp_set_nested(1) and omp_set_num_threads(2). The code is really simple. We simply spawn two threads on the division of the work. However, we insert a simple throttling logic to prevent excessive threads. Note that my experimentation showed decent speedups for this approach.

---

Finally, you may use OpenMP 3.0's task, where a task is a logically concurrent work. In the above all OpenMP's approaches, each parallel construct spawns two physical threads. You may say there is a hard 1-to-1 mapping between a task to a work thread. However, task separates logical tasks and workers.

Because OpenMP 3.0 is not popular yet, I will use Cilk Plus, which is great to express this kind of nested and recursive parallelisms. In Cilk Plus, the parallelization is extremely easy:

void qsort(int* begin, int* end)
{
  if (begin != end) {
    --end;
    int* middle = std::partition(begin, end,
      std::bind2nd(std::less<int>(), *end));
    std::swap(*end, *middle);

    cilk_spawn qsort(begin, middle);
    qsort(++middle, ++end);
    // cilk_sync; Only necessay at the final stage.
  }
}

I copied this code from Cilk Plus' example code. You will see a single keyword cilk_spawn is everything to parallelize quicksort. I'm skipping the explanations of Cilk Plus and spawn keyword. However, it's easy to understand: the two recursive calls are declared as logically concurrent tasks. Whenever the recursion takes place, the logical tasks are created. But, the Cilk Plus runtime (which implements an efficient work-stealing scheduler) will handle all kinds of dirty job. It optimally queues the parallel tasks and maps to the work threads.

Note that OpenMP 3.0's task is essentially similar to the Cilk Plus's approach. My experimentation shows that pretty nice speedups were feasible. I got a 3~4x speedup on a 8-core machine. And, the speedup was scale. Cilk Plus' absolute speedups are greater than those of OpenMP 3.0's.

The approach of Cilk Plus (and OpenMP 3.0) and your approach are essentially the same: the separation of parallel task and workload assignment. However, it's very difficult to implement efficiently. For example, you must reduce the contention and use lock-free data structures.