# Implementing Partition with Prefix Scan

In my last post, I sketched out the algorithm for the `Partition`

operation within a Quicksort algorithm using Prefix Scan to reduce the span of the parallel execution when executing on a manycore processor. In this post, I present one implementation of that algorithm using OpenMP for the first level of task parallelism and Intel Threading Building Blocks (TBB) to give me manycore parallel operations.

As you see my code, you will soon realize that this is an initial implementation. My goal was to write more of a "proof of concept" than a well-tuned effort. While it is convenient for having the parallel algorithms that I need, TBB may not be the most efficient means of putting into practice those parallel parts that can utilize many cores. Something more SIMD-like (e.g., CUDA, Intel Array Building Blocks) could be more effective. I'm sure you will see many different improvements to be made to my code. Please take what I've shown here and feel free to make improvements on it.

Overall, the `Quicksort()`

function is called from within an OpenMP parallel region using some subset of all possible threads. The distribution of partitions to be sorted by these threads is done via a shared TBB `concurrent_queue`

container, `Q`

. The indices of the subarray to be sorted are held in the `struct`

type `qSortIndex`

.

tbb::concurrent_queue<qSortIndex> Q;

Before getting to the `Quicksort()`

code, the first thing to examine is the TBB class for the `parallel_scan()`

function. The code here is really nothing special. The addition operation is almost the default operation, so I simply took the sample code from the TBB documentation. I did modify the pointer names for the input array (`idata`

) and output array (`odata`

), but the rest is just the same.

// Prefix Scan class class pScan { int sum; int* const odata; const int* const idata; public: pScan(int *odata_, const int *idata_) : sum(0), idata(idata_), odata(odata_) {} template<typename Tag> void operator()(const blocked_range<int> &r, Tag) { int temp = sum; for (int i = r.begin(); i < r.end(); ++i) { temp = temp + idata[i]; if(Tag::is_final_scan()) odata[i] = temp; } sum = temp; } pScan(pScan& b, split) : sum(0), idata(b.idata), odata(b.odata) {} void reverse_join(pScan& a) {sum = a.sum + sum;} void assign(pScan& b) {sum = b.sum;} };

The next six code blocks all come from the `Quicksort()`

function. I've broken them out into separate segments to discuss each one without making you flip back to a monolithic block of code, trying to figure out what lines I'm writing about. I've used an iterative version of Quicksort, so there are no recursive calls to be implemented as new tasks for threads. As subsets are partitioned, the new subset indices are put into the queue, `Q`

. Threads that are executing in the parallel region that is the `Quicksort()`

function will pull an index pair from the queue, partition that subset, and place the indices of the two new subpartitions into the queue until the count of sorted elements equals the total number of original items to be sorted.

First up is the function header and local declarations.

void QuickSort(int *A, int N) { int p, r, q; qSortIndex d, d1, d2; int *LT = new int[N]; int *GT = new int[N]; int *B = new int[N]; int *C = new int[N]; int *Ld = new int[N]; int *Gd = new int[N];

The function takes in an array (`A`

) and the number of items to be sorted (`N`

) within the array. This array will hold the sorted elements when the function completes. The first two lines of declarations are for the index bounds of the partitioned subsets found and the `struct`

s that will hold these index values in and out of the shared queue.

The six arrays, each the same size as the original array, will hold the mask vector for less-than (`LT`

) and greater-than (`GT`

) elements compared to the pivot value, the results of the less-than prefix scan (`B`

) and greater-than prefix scan (`C`

), and the subset of elements less than (`Ld`

) and greater than (`Gd`

) the pivot value after the data is moved from `A`

.

while (1) { while (!Q.try_pop(d)) { if (Done) return; // if Done, no more need to process anything } p = d.lo; r = d.hi;

The iterative Quicksort algorithm uses an infinite loop, via the outer `while`

-loop, to keep threads busy. The first thing to be done in an iteration of this outer loop is to attempt to pull something from the queue of index pairs. If the queue is empty, the shared `Done`

flag is checked to determine if the sort has been completed. If the queue holds at least one thing, the thread that retrieves an index pair will know the indices of the next subset within `A`

to be partitioned.

if (p < r) { int *At = &A[p]; int pivot = At[0]; int n = r-p+1; parallel_for(blocked_range<int>(1, n), [&](const blocked_range<int>& r){ for (int i = r.begin(); i < r.end(); ++i) { if (At[i] < pivot) { LT[i] = 1; GT[i] = 0; } else { LT[i] = 0; GT[i] = 1; } } } ); LT[0] = GT[0] = -1;