Extra Horsepower

The graph shows the relative performance (for sorting an array of 10 million integers) for Quicksort, versus 2-way, 3-way, and 4-way mergesorts that parallelize the last (and only the last) merge operation at each level of merging.

Before parallelizing 3- and 4-way merges (yuu-uuk), I wondered... would it be worth parallelizing all the 2-way-merge operations, not just the last one at each level of merging (and using a StackDispatcher too)...? The answer turned out to be... no. It's very slightly slower (possibly because the extra synchronization overhead, plus the binary-search overhead for finding where to split up the inputs when dividing up the processing, outweighs the cache benefit, if any).

I also figured it might be a good idea to *check* that stability is preserved. Luckily that's pretty easy, with template sort routines. I declared a type, Horse, generated arrays of Horse (see below) (by setting the value members, randomly permuting, and then writing the starting index of each element into the satellite members), and checked that stability was preserved after sorting.

template<class V, class S> class ValueAndSatellite
{
public:
  V value;
  S satellite;
  ValueAndSatellite()         { }
  ValueAndSatellite(V v, S s) { value=v; satellite=s; }
  ~ValueAndSatellite()        { }
  bool operator<  (const ValueAndSatellite<V,S> &rhs) { return value<rhs.value; }
  bool operator<= (const ValueAndSatellite<V,S> &rhs) { return value<=rhs.value; }
  bool operator== (const ValueAndSatellite<V,S> &rhs) { return value==rhs.value; }
};
typedef ValueAndSatellite<int, int> Horse;

Stability was preserved, for the 2- and 3-way mergesorts, but not for the 4-way mergesort. I haven't yet figured out how the 4-way mergesort went wrong. I'll have to post a correction when I find out.

Incidentally, with an array of Horse (each Horse is twice the size of an int, but has approximately the same comparison cost), Quicksort beat two-way, three-way, and even four-way merging comfortably. Even if all of the value members of the Horses in the array were distinct, Quicksort was 11% faster than a four-way merge.






For 3- and 4-way merging I've opted not to worry so much about trying to "split evenly in the presence of large numbers of duplicates". It's too hard to think through all the cases. Instead, I... just hope the splits are about right ("on average, if there aren't too many duplicates, they'll be okay"). And, at least for sorting integers on 4 cores,

• my 3-way mergesort beats Quicksort comfortably. But not by the 8% I had expected. Only by about 5%.


• my 4-way mergesort beats Quicksort by about 9% (I had expected 8%). But then, it isn't stable (ouch!)

Here's my current code (less classes, functions, and so forth that have appeared in previous posts)...

template<class T> void mergeForecastBackwardsRadix2External(T *a, T* aStop, T* b, T* bStop, T* dest)
{
  if (aStop[1]<=bStop[1])
  {
    //the "b"s will run out before the "a"s.
    for (;b>bStop;*dest--=*b--)
    {
      for (;*b<*a;*dest--=*a--);
    }
    for (;a>aStop;*dest--=*a--);
  }
  else
  {
    for (;a>aStop;*dest--=*a--)
    {
      for (;*a<=*b;*dest--=*b--);
    }
    for (;b>bStop;*dest--=*b--);
  }
}


template <class T> inline int indexOfFirstElementInBlockGreaterOrEqualTo(T *block, int count, T &v)
{
  int loIndex = 0;
  int hiIndex = count;
  while (loIndex<hiIndex)
  {
    int tryIndex = ( loIndex + hiIndex ) / 2;
    if ( v <= block[tryIndex]  )
      hiIndex  = tryIndex;
    else
      loIndex = tryIndex + 1;
  }
  return hiIndex;
}

template <class T> inline int indexOfFirstElementInBlockGreaterThan(T *block, int count, T &v)
{
  int loIndex = 0;
  int hiIndex = count;
  while (loIndex<hiIndex)
  {
    int tryIndex = ( loIndex + hiIndex ) / 2;
    if ( v < block[tryIndex]  )
      hiIndex  = tryIndex;
    else
      loIndex = tryIndex + 1;
  }
  return hiIndex;
}

template <class T> void paranoidMergeRadix2External(T *a, int aCount, T* b, int bCount, T* dest, bool bLeftToRight, INotifiable *t)
{
  if (aCount==0)
  {
    for (;bCount>0;--bCount) *dest++=*b++;
  }
  else if (bCount==0)
  {
    for (;aCount>0;--aCount) *dest++=*a++;
  }
  else if (bLeftToRight)
  {
    mergeForecastRadix2External(a, a+aCount, b, b+bCount, dest);
  }
  else
  {
    mergeForecastBackwardsRadix2External(a+aCount-1, a-1, b+bCount-1, b-1, dest+aCount+bCount-1);
  }
  if (t!=NULL)
  {
    t->notify();
  }
}

template <class T> void parallelMergeRadix2External(T* a, int aCount, T* b, int bCount, T* dest, bool bLeftToRight, int processorCount, INotifiable *t)
{
  while (processorCount>1)
  {
    //Ick!  We need to split the merge into two parts.
    //We'll be kinda lazy.  We'll figure out where the median of the LHS would go in the RHS, and split that way.
    int middleLeft  = aCount/processorCount*(processorCount/2);
    int hiLeft = middleLeft+1+indexOfFirstElementInBlockGreaterThan(a+middleLeft+1, aCount-middleLeft-1, a[middleLeft]);  
      //index of first element >v, in a, or aCount if there is none
    int loRight   = indexOfFirstElementInBlockGreaterOrEqualTo(b, bCount, a[middleLeft]);
      //index of first element >=v, in b, or bCount if there is none

    //Normally: merge a[0..hiLeft-1] with b[0..loRight-1], so set... takeLeft=hiLeft.
    //But:      We want takeLeft + loRight to be roughly (aCount+bCount)/processorCount*(proessorCount/2), and takeLeft
    //      can be any value between loLeft and hiLeft.  

    int takeLeft     = hiLeft;
    int desiredSplit = (aCount+bCount)/processorCount*(processorCount/2);
    int loLeft       = -1;
    if (takeLeft + loRight > desiredSplit) //too many records will be going to the "left" merge operation
    {
      if (desiredSplit - loRight >= middleLeft)
        takeLeft = desiredSplit - loRight;    //since any value between middleLeft and hiLeft is fine for # records to take from a
      else
      {
        //otherwise: need to know the minimum number of records we could *safely* take from block "a".
        loLeft = indexOfFirstElementInBlockGreaterOrEqualTo(a, middleLeft, a[middleLeft]);
        if (desiredSplit - loRight > loLeft)
        {
          takeLeft = desiredSplit - loRight;
        }
        else
        {
          takeLeft = loLeft;
        }
      }
    }    

    parallelMergeRadix2External(a, takeLeft, b, loRight, dest, bLeftToRight, processorCount/2, t);
    a              += takeLeft;
    b              += loRight;
    aCount         -= takeLeft;
    bCount         -= loRight;
    dest           += takeLeft+loRight;
    processorCount -= processorCount/2;
  }
  gpDispatcher->queueOrDo(CallAction7(paranoidMergeRadix2External<T>, a, aCount, b, bCount, dest, bLeftToRight, t));
}

template <class T> void parallelMergesortSubArrayRadix2(T *src, T *dest, int elementCount
                            , int processorCount, int multiThreadingThreshold
                            , bool srcIsInput, bool isLastMerge, INotifiable *t)
{
  if (elementCount>multiThreadingThreshold)
  {  
    int halfCount = elementCount >> 1;
    IAction*     mergeAction;       //the merge itself
    INotifiable* pendingMerge;    //for notifying the count-down latch that prevents the merge running                

    if (isLastMerge)
    {
      INotifiable* notificationDemultiplexer = NewPendingNotification(processorCount, t);
      mergeAction = CallAction8(parallelMergeRadix2External<T>, src, halfCount, src+halfCount, elementCount-halfCount
              , dest, srcIsInput, processorCount, notificationDemultiplexer);
      pendingMerge = NewPendingAction(2, mergeAction);      

    }
    else 
    {
      if (srcIsInput)
      {
        mergeAction = CallAction5(mergeForecastRadix2External<T>,  src, src+halfCount, src+halfCount, src+elementCount, dest);
      }
      else
      {
        mergeAction = CallAction5(mergeForecastBackwardsRadix2External<T>, src+halfCount-1, src-1, src+elementCount-1, src+halfCount-1, dest+elementCount-1);
      }
      IAction* notifyingAction =  NewNotifyingAction(mergeAction,t);
      pendingMerge = NewPendingAction(2, notifyingAction);
    }        

    parallelMergesortSubArrayRadix2(dest, src, halfCount
      , processorCount, multiThreadingThreshold, !srcIsInput, false, pendingMerge);
    parallelMergesortSubArrayRadix2(dest+halfCount, src+halfCount, elementCount-halfCount
      , processorCount, multiThreadingThreshold, !srcIsInput, isLastMerge, pendingMerge);

  }
  else
  {
    gpDispatcher->queueOrDo( NewNotifyingAction( CallAction5(mergesortForecastRadix2External<T>, src, elementCount, dest, srcIsInput, 32), t));
  }
}

template <class T>   void parallelMergesort(T *base, int elementCount, int processorCount, int multiThreadingThreshold)
{
  T *workArea = new T[elementCount];
  CountDown t(1);
  gpDispatcher->queueOrDo(CallAction8(parallelMergesortSubArrayRadix2<T>,workArea, base, elementCount, processorCount, multiThreadingThreshold, false, 
    true /*set to false to turn off parallization of the last merge on each level*//, (INotifiable*)&t));
  t.wait();
  delete [] workArea;
}

template <class T> void paranoidMergeRadix3External(T *a, int aCount, T* b, int bCount, T* c, int cCount, T* dest, bool bLeftToRight, INotifiable *t)
{
  if (aCount==0) 
    paranoidMergeRadix2External(b, bCount, c, cCount, dest, bLeftToRight, NULL);
  else if (bCount==0)
    paranoidMergeRadix2External(a, aCount, c, cCount, dest, bLeftToRight, NULL);
  else if (cCount==0)
    paranoidMergeRadix2External(a, aCount, b, bCount, dest, bLeftToRight, NULL);
  else if (bLeftToRight)
    mergeRadix3FastForward(a, a+aCount, b, b+bCount, c, c+cCount, dest);
  else
    mergeRadix3FastBackward(a+aCount-1, a-1, b+bCount-1, b-1, c+cCount-1, c-1, dest+aCount+bCount+cCount-1);

  if (t!=NULL)
  {
    t->notify();
  }
}

template <class T> void parallelMergeRadix3External(T* a, int aCount, T* b, int bCount, T* c, int cCount, T* dest, bool bLeftToRight, int processorCount, INotifiable *t)
{
  while (processorCount > 1)
  {
    int aMiddle  = aCount/processorCount*(processorCount/2);
    int bMiddle  = indexOfFirstElementInBlockGreaterThan(b, bCount, a[aMiddle]);
    int cMiddle  = indexOfFirstElementInBlockGreaterThan(c, cCount, a[aMiddle]);
        aMiddle += indexOfFirstElementInBlockGreaterThan(a+aMiddle, aCount-aMiddle, a[aMiddle]);
    parallelMergeRadix3External(a, aMiddle, b, bMiddle, c, cMiddle, dest, bLeftToRight, processorCount/2, t);
    a += aMiddle; aCount -= aMiddle;
    b += bMiddle; bCount -= bMiddle;
    c += cMiddle; cCount -= cMiddle;
    dest += aMiddle + bMiddle + cMiddle;
    processorCount -= processorCount/2;
  }
  gpDispatcher->queueOrDo(CallAction9(paranoidMergeRadix3External<T>, a, aCount, b, bCount, c, cCount, dest, bLeftToRight, t));
}


template <class T> void parallelMergesortSubArrayRadix3(T *src, int elementCount, T *dest, bool isSourceInput,
                          bool isLastMerge, int processorCount, int multiThreadingThreshold,
                          INotifiable* t)
{
  if (elementCount>multiThreadingThreshold)
  {  
    int stepSize = elementCount / 3 ;
    int twoStep  = stepSize + stepSize;
    IAction *mergeAction;
    INotifiable *pendingMerge;

    if (isLastMerge)
    {
      INotifiable *notificationDemultiplexer = NewPendingNotification(processorCount, t);
      mergeAction = CallAction10(parallelMergeRadix3External<T>, src, stepSize, src+stepSize, stepSize, src+twoStep, elementCount-twoStep
              , dest, isSourceInput, processorCount, notificationDemultiplexer);
      pendingMerge = NewPendingAction(3, mergeAction);      
    }
    else 
    {
      if (isSourceInput)
      {
        mergeAction = CallAction7(voidMergeRadix3FastForward<T>,  src, src+stepSize, src+stepSize, src+stepSize+stepSize, src+stepSize+stepSize,
                  src+elementCount, dest);
      }
      else
      {
        mergeAction = CallAction7(voidMergeRadix3FastBackward<T>, src+stepSize-1, src-1, src+stepSize+stepSize-1, src+stepSize-1,
                      src+elementCount-1, src+stepSize+stepSize-1, dest+elementCount-1);
      }
      IAction* notifyingAction =  NewNotifyingAction(mergeAction,t);
      pendingMerge = NewPendingAction(3, notifyingAction);
    }

    parallelMergesortSubArrayRadix3(dest, stepSize, src
      , !isSourceInput , false,       processorCount, multiThreadingThreshold, pendingMerge);
    parallelMergesortSubArrayRadix3(dest+stepSize, stepSize, src+stepSize
      , !isSourceInput , false,       processorCount, multiThreadingThreshold, pendingMerge);
    parallelMergesortSubArrayRadix3(dest+stepSize+stepSize, elementCount-stepSize-stepSize, src+stepSize+stepSize
      , !isSourceInput , isLastMerge, processorCount, multiThreadingThreshold, pendingMerge);
  }
  else
  {
    gpDispatcher->queueOrDo( NewNotifyingAction( CallAction4(mergeSortExternalRadix3Fast<T>, src, elementCount, dest, isSourceInput), t));
  }
}


template <class T>   void parallelMergesortRadix3(T *base, int elementCount, int processorCount, int multiThreadingThreshold=50000)
{
  T *workArea = new T[elementCount];
  CountDown t(1);
  gpDispatcher->queueOrDo(CallAction8(parallelMergesortSubArrayRadix3<T>, workArea, elementCount, base, false, true, processorCount, multiThreadingThreshold, (INotifiable*)&t));
  t.wait();
  delete [] workArea;
}

template <class T> void paranoidMergeRadix4External(T *a, int aCount, T* b, int bCount, T* c, int cCount, T* d, int dCount, T* dest, bool bLeftToRight, INotifiable *t)
{
  if (aCount==0) 
    paranoidMergeRadix3External(b, bCount, c, cCount, d, dCount, dest, bLeftToRight, NULL);
  else if (bCount==0)
    paranoidMergeRadix3External(a, aCount, c, cCount, d, dCount, dest, bLeftToRight, NULL);
  else if (cCount==0)
    paranoidMergeRadix3External(a, aCount, b, bCount, d, dCount, dest, bLeftToRight, NULL);
  else if (dCount==0)
    paranoidMergeRadix3External(a, aCount, b, bCount, c, cCount, dest, bLeftToRight, NULL);
  else if (bLeftToRight)
    mergeRadix4FastForward(a, a+aCount, b, b+bCount, c, c+cCount, d, d+dCount, dest);
  else
    mergeRadix4FastBackward(a+aCount-1, a-1, b+bCount-1, b-1, c+cCount-1, c-1, d+dCount-1, d-1, dest+aCount+bCount+cCount+dCount-1);

  if (t!=NULL)
  {
    t->notify();
  }
}

template <class T> void parallelMergeRadix4External(
  T* a, int aCount, T* b, int bCount, T* c, int cCount, T* d, int dCount,
  T* dest, bool bLeftToRight, int processorCount, INotifiable *t)
{
  while (processorCount > 1)
  {
    int aMiddle  = aCount/processorCount*(processorCount/2);
    int bMiddle  = indexOfFirstElementInBlockGreaterThan(b, bCount, a[aMiddle]);
    int cMiddle  = indexOfFirstElementInBlockGreaterThan(c, cCount, a[aMiddle]);
    int dMiddle  = indexOfFirstElementInBlockGreaterThan(d, dCount, a[aMiddle]);
        aMiddle += indexOfFirstElementInBlockGreaterThan(a+aMiddle, aCount-aMiddle, a[aMiddle]);
    parallelMergeRadix4External(a, aMiddle, b, bMiddle, c, cMiddle, d, dMiddle, dest, bLeftToRight, processorCount/2, t);
    a += aMiddle; aCount -= aMiddle;
    b += bMiddle; bCount -= bMiddle;
    c += cMiddle; cCount -= cMiddle;
    d += dMiddle; dCount -= dMiddle;
    dest += aMiddle + bMiddle + cMiddle + dMiddle;
    processorCount -= processorCount/2;
  }
  gpDispatcher->queueOrDo(CallAction11(paranoidMergeRadix4External<T>, a, aCount, b, bCount, c, cCount, d, dCount, dest, bLeftToRight, t));
}

template <class T> void parallelMergesortSubArrayRadix4(T *src, int elementCount, T *dest, bool isSourceInput,
                          bool isLastMerge, int processorCount, int multiThreadingThreshold,
                          INotifiable* t)
{
  if (elementCount>multiThreadingThreshold)
  {  
    int stepSize  = elementCount >> 2 ;
    int twoStep   = stepSize + stepSize;
    int threeStep = twoStep  + stepSize;
    IAction *mergeAction;
    INotifiable *pendingMerge;

    if (isLastMerge)
    {
      INotifiable *notificationDemultiplexer = NewPendingNotification(processorCount, t);
      mergeAction = CallAction12(parallelMergeRadix4External<T>, src, stepSize, src+stepSize, stepSize
              , src+twoStep, stepSize, src+threeStep, elementCount-threeStep
              , dest, isSourceInput, processorCount, notificationDemultiplexer);
      pendingMerge = NewPendingAction(4, mergeAction);      
    }
    else 
    {
      if (isSourceInput)
      {
        mergeAction = CallAction9(voidMergeRadix4FastForward<T>,  src, src+stepSize, src+stepSize, src+twoStep, src+twoStep,
                  src+threeStep, src+threeStep, src+elementCount, dest);
      }
      else
      {
        mergeAction = CallAction9(voidMergeRadix4FastBackward<T>, src+stepSize-1, src-1, src+twoStep-1, src+stepSize-1,
                  src+threeStep-1, src+twoStep-1, src+elementCount-1, src+threeStep-1, dest+elementCount-1);
      }
      IAction* notifyingAction =  NewNotifyingAction(mergeAction,t);
      pendingMerge = NewPendingAction(4, notifyingAction);
    }

    parallelMergesortSubArrayRadix4(dest, stepSize, src
      , !isSourceInput , false,       processorCount, multiThreadingThreshold, pendingMerge);
    parallelMergesortSubArrayRadix4(dest+stepSize, stepSize, src+stepSize
      , !isSourceInput , false,       processorCount, multiThreadingThreshold, pendingMerge);
    parallelMergesortSubArrayRadix4(dest+twoStep, stepSize, src+twoStep
      , !isSourceInput , false,       processorCount, multiThreadingThreshold, pendingMerge);
    parallelMergesortSubArrayRadix4(dest+threeStep, elementCount-threeStep, src+threeStep
      , !isSourceInput , isLastMerge, processorCount, multiThreadingThreshold, pendingMerge);
  }
  else
  {
    gpDispatcher->queueOrDo( NewNotifyingAction( CallAction4(mergeSortExternalRadix4Fast<T>, src, elementCount, dest, isSourceInput), t));
  }
}

template <class T>   void parallelMergesortRadix4(T *base, int elementCount, int processorCount, int multiThreadingThreshold=50000)
{
  T *workArea = new T[elementCount];
  CountDown t(1);
  gpDispatcher->queueOrDo(CallAction8(parallelMergesortSubArrayRadix4<T>, workArea, elementCount, base, false, true, 
    processorCount, multiThreadingThreshold, (INotifiable*)&t));
  t.wait();
  delete [] workArea;
}

Tuning Array Mergesorts for speed - Parallel Merging

The multithreaded mergesort routines that I posted here a few days ago only multithreaded the sorting of partitions. They did not multithread the actual merge operations. I was lazy. How much did that omission cost? Well, I added some code to generate a gantt chart for each of the five threads involved in the four-core two-way mergesort. I'm not going to post the code that generates the gantt chart, or even show the chart, but the important thing that it showed was that the last action queued on the last worker thread to complete took 25% of the elapsed time (to sort 10 million records, with a multithreading threshold of 50,000, using 4 cores), because it was merging on one thread (and it wasn't just merging one level on one thread! It was merging at lots of levels).

Rewriting the algorithm so that the last merge operation at each level of merging (and *only* the last merge) can be split across multiple cores does indeed reduce the running time siginifantly (and if the merging is multithreaded, the speed-up for the overall sort, at least for sorting distinct integers, is much better. For 10 million integers, two, three and four core versions are roughly 1.98, 2.88, and 3.68 times faster, respectively, than the one-core version).

Way back in 2000, when I looked at multithreading the merge operations for two-way mergesort on a Pentium IV (which only had "a core and a half", given the limitations of the Pentium IV's hyperthreading), the performance benefit was marginal. But I should have done the calculations to predict the saving. Assuming all comparisons take equal time and that the comparisons required for the last few merge operations can be spread evenly across all the processors... the speed-up (on 4-cores) for 3- and 4-way mergesorts should be enough to make them both about 8% faster than Quicksort. Well, I don't think that can be right (Quicksort always wins!). But I'll have to code them to find out for sure.

I'm not certain the following code is correct (I haven't tested all the boundary cases yet. I think it will go wrong when hiRight==bCount), nor have I properly checked whether it even preserves stability, and I haven't written the code to split 3- and 4-way merges across multiple cores either.



namespace
{
  class PendingNotification: public INotifiable
  {
    protected:
      ThreadSafeCounter m_prerequisiteCount;
      INotifiable* m_target;
    public:
      PendingNotification(int prereqCount, INotifiable *target)
        : m_prerequisiteCount(prereqCount)
        , m_target(target)
      {
      }
      void notify()
      {
        if (--m_prerequisiteCount == 0 )
        {
          m_target->notify();
          delete this;
        }
      }
      bool notified()
      {
        return false;  //if it still exists, it hasn't been notified.
      }
      void wait()
      {
        return m_target->wait();  //it doesn't really make sense to wait on a PendingNotification
      }
  };
}

INotifiable *NewPendingNotification(int prerequisiteCount, INotifiable *target)
{
  return new PendingNotification(prerequisiteCount, target);
}

template <class T> void paranoidMergeRadix2External(T *a, int aCount, T* b, int bCount, T* dest, bool bLeftToRight, INotifiable *t)
{
  if (aCount==0)
  {
    for (;bCount>0;--bCount) *dest++=*b++;
  }
  else if (bCount==0)
  {
    for (;aCount>0;--aCount) *dest++=*a++;
  }
  else if (bLeftToRight)
  {
    mergeForecastRadix2External(a, a+aCount, b, b+bCount, dest);
  }
  else
  {
    mergeForecastBackwardsRadix2External(a+aCount-1, a-1, b+bCount-1, b-1, dest+aCount+bCount-1);
  }
  if (t!=NULL)
  {
    t->notify();
  }
}

template <class T> void parallelMergeRadix2External(T* a, int aCount, T* b, int bCount, T* dest, bool bLeftToRight, int processorCount, INotifiable *t)
{
  while (processorCount>1)
  {
    //Ick!  We need to split the merge into two parts.
    //We'll be kinda lazy.  We'll figure out where the median of the LHS would go in the RHS, and split that way.
    int middleLeft  = aCount/processorCount*(processorCount/2);
    int lowRight    = 0;
    int hiRight     = bCount;
    while (lowRight<hiRight)
    {
      int tryRight = ( lowRight + hiRight ) / 2;
      if ( a[middleLeft] <= b[tryRight]  )
        hiRight  = tryRight;
      else
        lowRight = tryRight + 1;
    }
    //At this point, we have either hiRight==bCount, or b[hiRight] < a[middleLeft] 
    while (middleLeft<aCount && a[middleLeft+1]<=b[hiRight]) middleLeft++;
    parallelMergeRadix2External(a, middleLeft+1, b, hiRight, dest, bLeftToRight, processorCount/2, t);
    a              += middleLeft+1;
    b              += hiRight;
    aCount         -= middleLeft+1;
    bCount         -= hiRight;
    dest           += middleLeft+hiRight+1;
    processorCount -= processorCount/2;
  }
  gpDispatcher->queueOrDo(CallAction7(paranoidMergeRadix2External<T>, a, aCount, b, bCount, dest, bLeftToRight, t));
}

template <class T> void parallelMergesortSubArrayRadix2(T *src, T *dest, int elementCount
                            , int processorCount, int multiThreadingThreshold
                            , bool srcIsInput, bool isLastMerge, INotifiable *t)
{
  if (elementCount>multiThreadingThreshold)
  {  
    int halfCount = elementCount >> 1;
    IAction*     mergeAction;       //the merge itself
    INotifiable* pendingMerge;    //for notifying the count-down latch that prevents the merge running                

    if (isLastMerge)
    {
      INotifiable* notificationDemultiplexer = NewPendingNotification(processorCount, t);
      mergeAction = CallAction8(parallelMergeRadix2External<T>, src, halfCount, src+halfCount, elementCount-halfCount
              , dest, srcIsInput, processorCount, notificationDemultiplexer);
      pendingMerge = NewPendingAction(2, mergeAction);      

    }
    else 
    {
      if (srcIsInput)
      {
        mergeAction = CallAction5(mergeForecastRadix2External<T>,  src, src+halfCount, src+halfCount, src+elementCount, dest);
      }
      else
      {
        mergeAction = CallAction5(mergeForecastBackwardsRadix2External<T>, src+halfCount-1, src-1, src+elementCount-1, src+halfCount-1, dest+elementCount-1);
      }
      IAction* notifyingAction =  NewNotifyingAction(mergeAction,t);
      pendingMerge = NewPendingAction(2, notifyingAction);
    }        

    parallelMergesortSubArrayRadix2(dest, src, halfCount
      , processorCount, multiThreadingThreshold, !srcIsInput, false, pendingMerge);
    parallelMergesortSubArrayRadix2(dest+halfCount, src+halfCount, elementCount-halfCount
      , processorCount, multiThreadingThreshold, !srcIsInput, isLastMerge, pendingMerge);

  }
  else
  {
    gpDispatcher->queueOrDo( NewNotifyingAction( CallAction5(mergesortForecastRadix2External<T>, src, elementCount, dest, srcIsInput, 32), t));
  }
}

template <class T>   void parallelMergesort(T *base, int elementCount, int processorCount, int multiThreadingThreshold)
{
  T *workArea = new T[elementCount];
  CountDown t(1);
  gpDispatcher->queueOrDo(CallAction8(parallelMergesortSubArrayRadix2<T>,workArea, base, elementCount, processorCount, multiThreadingThreshold, false, true, (INotifiable*)&t));
  t.wait();
  delete [] workArea;
}

Multithreading other sort routines

With all that infrastructure for writing multithreading sorts available (see my last few posts), it's not difficult to write multithreaded versions of other sorting algorithms. This post will cover multithreaded versions of combsort and shellsort. I can warn in advance that... they suck. I've graphed them running on one, two, three and four cores, against one- (and sometimes two-) core classic Quicksort. They don't do very well at all.

It's easy to see how to parallelize shellsort and combsort. Both algorithms use a series of "gaps" (or "intervals"). For each gap, g shellsort fully sorts each "slice" of elements with indexes separated by g, and combsort partially sorts each "slice". For each gap, g, each of the g separate slices can be processed independently. But they're slow. The examples I provide here could be tuned better.

You can see how they might be tuned better by looking at, for example, the fourth graph, showing a breakdown of the time spent for each combsort pass. There, it is clear that an insertion sort should be used much sooner, say when the gap reaches 17. But that wouldn't really improve things much (my own tests suggest by about 17%, for sorting 9 million integers on 4 cores). Shellsort and Combsort just do not parallelize well in practice. Shellsort in particular parallelizes poorly.

There are two versions of Combsort shown here, because the stock-standard Combsort has truly awful performance on average (and a very bad worst-case). A slightly different version that alternates between left-to-right and right-to-left scans has a markedly better average case (and a way better worst case).

I'll start with a parallel insertion sort (!), because that's how you finish off a parallel combsort or a parallel shell sort. This implementation is only worth calling if the array is "nearly" sorted, as is the case "at the tail end" of a shellSort or a combSort.

template <class T>   void parallelInsertionSort(T *base, int elementCount, int processorCount, int multiThreadingThreshold=50000
                         , INotifiable *notifyMe=NULL)
{
  CountDown task;
  if (elementCount<processorCount*64 || processorCount<2)
  {
    insertionSort(base, elementCount);
  }
  else
  {
    task += processorCount;
    for (int p=0;p<processorCount;++p)
    {
      int i = (elementCount * p     ) / processorCount;
      int j = (elementCount * (p+1) ) / processorCount;
      gpDispatcher->queueOrDo( NewNotifyingAction(CallAction2(insertionSort<T>,base+i,j-i), &task ) );
    }
    task.wait();
    for (int p=1;p<processorCount;++p)
    {
      int i = (elementCount * p     ) / processorCount;
      int k = (elementCount * (p+1) ) / processorCount;
      T v;
      //insert the elements in the range [i] through [k-1]
      //into the range [0] through [i-1], given that the
      //two ranges, [0..(i-1)] and [i..(k-1)] are 
      //independently sorted.
      for (int j=i;j<k;++j)
      {
        v = base[j];
        int h;
        for (h=j;h>0;--h)
          if (v<base[h-1])
            base[h] = base[h-1];
          else
            break;
        if (h==j)
        {
          //as soon as the element that *was* at [i-1]
          //is placed correctly, it is not necessary to
          //look at [(i+1)..(k-1)].  They're placed
          //where they should be already.
          break;
        }
        base[h] = v;
      }
    }    
  }
  if (notifyMe!=NULL)
    notifyMe->notify();
}

Now for the parallel shellsort:

template <class T> void shellSortSlices(T *base, int elementCount, int step, int firstChain, int lastChain)
{
  for (int i=firstChain+step;i<elementCount;i+=step)
  {
    int k=i + lastChain - firstChain;
    if (elementCount<=k) 
    {
      k=elementCount-1;
    }
    for (int j=i;j<=k;++j)
    {
      int h;
      T v = base[j];
      for (h=j;h>=step;h-=step)
        if (v<base[h-step]) 
          base[h]=base[h-step];
        else
          break;
      base[h] = v;
    }
  }
}

template <class T>   void parallelShellSort(T *base, int elementCount, int processorCount, int multiThreadingThreshold)
{      
  int gaps[] = { 5,13,37,97,251,631,1579,3947,9871,24677,61703,154267,385709,964283,2410711,
          6026827,15067067,37667713,94169297,235423249,588558169,1471395493 };
  int f;
  for (f=0; gaps[f]<elementCount; f++);
  f--;  
  CountDown task;
  for (;f>=0;f--)
  {
    int g = gaps[f];
    task += processorCount;
    for (int p=0; p<processorCount; p++)
    {
      int i =  p*g / processorCount;
      int j = (p+1)*g / processorCount - 1;
      gpDispatcher->queueOrDo ( NewNotifyingAction(CallAction5(shellSortSlices<T>, base, elementCount, g, i, j), &task));
    }        
    task.wait();
  }
  parallelInsertionSort(base, elementCount, processorCount, multiThreadingThreshold, &task);  
}

And the parallel combsort:


int alternatingGaps[] = {2,3,5,7,11,17,23,37,53,79,113,163,229,331,463,653,
  919,1289,1811,2539,3557,4987,6983,9781,13693,19181,26861,37607,52667,
  73751,103289,144611,202471,283463,396871,555637,777901,1089091,1524763,
  2134697,2988607,4184087,5857727,8200847,11481199,16073693,22503181,31504453,
  44106241,61748749,86448259,121027561,169438579,237214009,332099617,464939459,
  650915219,911281291,1275793807,1786111267};
  //These are primes near powers of P=1.4
  //Alternating combsort can make do with slightly fewer passes than
  //vanilla combsort, as it does better at moving elements that belong
  //low, low.  Vanilla combsort is nowhere near as good at that.
  //Vanilla combsort moves elements that belong high, high, very efficiently,
  //and "hopes the values that belong low take care of themselves".
  //But they don't.

int combSortIntervals[] = {
  2,3,5,7,11,17,23,29,37,53,71,97,127,167,223,293,383,499,653,853,
  1109,1447,1889,2459,3203,4177,5431,7069,9199,11959,15551,20219,  
  26293,34183,44449,57787,75133,97673,126989,165089,214631,279023,
  362741,471571,613049,796967,1036067,1346899,1750979,2276293,2959183,
  3846943,5001049,6501367,8451791,10987343,14283559,18568637,
  24139231,31381043,40795357,53033969,68944157,89627413,116515667,
  151470353,196911461,255984881,332780323,432614407,562398721,731118341,
  950453827,1235589917,1606266791,2088146713 };
  //These are the primes near powers of P=1.3
  //Combsort runs very slightly faster, on average, if all the 
  //intervals are prime.

template <class T> void combPassSlices(T *base, int elementCount, int step, int firstChain, int lastChain)
{
  T v;
  for (int i=firstChain+step;i<elementCount;i+=step)
  {
    int k=i + lastChain - firstChain;
    if (elementCount<=k) 
      k=elementCount-1;
    for (int j=i;j<=k;++j)
    {
      if (base[j]<base[j-step])
      {
        v = base[j];
        base[j] = base[j-step];
        base[j-step] = v;
      }
    }
  }
}

template <class T> void downwardCombPassSlices(T *base, int elementCount, int step, int firstChain, int lastChain)
{
  T v;
  int h,i,j;
  for (j=elementCount-firstChain-1;j>=step;j-=step)
  {
    h = j + firstChain - lastChain;
    if (h<step) 
    {
      h=step;
    }
    for (i=j;i>=h;--i)
    {
      if (base[i]<base[i-step])
      {
        v = base[i];
        base[i] = base[i-step];
        base[i-step] = v;
      }
    }
  }
}

template <class T>   void parallelCombSort(T *base, int elementCount, int processorCount, int multiThreadingThreshold, bool alternating)
{    
  int *gaps = combSortIntervals;
  if (alternating)
    gaps = alternatingGaps;

  int f;
  for (f=0; gaps[f]<elementCount; f++);
  f--;  
  CountDown task;
  for (;f>=0;f--)
  {
    int g = gaps[f];
    task += processorCount;
    for (int p=0; p<processorCount; p++)
    {
      int i =  p*g / processorCount;
      int j = (p+1)*g / processorCount - 1;
      if (alternating && (f&1)!=0)
        gpDispatcher->queueOrDo (NewNotifyingAction(CallAction5(downwardCombPassSlices<T>, base, elementCount, g, i, j), &task));
      else
        gpDispatcher->queueOrDo (NewNotifyingAction(CallAction5(combPassSlices<T>, base, elementCount, g, i, j), &task));
    }        
    task.wait();
  }
  parallelInsertionSort(base, elementCount, processorCount, multiThreadingThreshold);  
}