4.4.3/parallel/balanced_quicksort.h

// -*- C++ -*-

// Copyright (C) 2007, 2008, 2009 Free Software Foundation, Inc.
//
// This file is part of the GNU ISO C++ Library.  This library is free
// software; you can redistribute it and/or modify it under the terms
// of the GNU General Public License as published by the Free Software
// Foundation; either version 3, or (at your option) any later
// version.

// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.

// Under Section 7 of GPL version 3, you are granted additional
// permissions described in the GCC Runtime Library Exception, version
// 3.1, as published by the Free Software Foundation.

// You should have received a copy of the GNU General Public License and
// a copy of the GCC Runtime Library Exception along with this program;
// see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
// <http://www.gnu.org/licenses/>.

/** @file parallel/balanced_quicksort.h
 *  @brief Implementation of a dynamically load-balanced parallel quicksort.
 *
 *  It works in-place and needs only logarithmic extra memory.
 *  The algorithm is similar to the one proposed in
 *
 *  P. Tsigas and Y. Zhang.
 *  A simple, fast parallel implementation of quicksort and
 *  its performance evaluation on SUN enterprise 10000.
 *  In 11th Euromicro Conference on Parallel, Distributed and
 *  Network-Based Processing, page 372, 2003.
 *
 *  This file is a GNU parallel extension to the Standard C++ Library.
 */

// Written by Johannes Singler.

#ifndef _GLIBCXX_PARALLEL_BALANCED_QUICKSORT_H
#define _GLIBCXX_PARALLEL_BALANCED_QUICKSORT_H 1

#include <parallel/basic_iterator.h>
#include <bits/stl_algo.h>

#include <parallel/settings.h>
#include <parallel/partition.h>
#include <parallel/random_number.h>
#include <parallel/queue.h>
#include <functional>

#if _GLIBCXX_ASSERTIONS
#include <parallel/checkers.h>
#endif

namespace __gnu_parallel
{
/** @brief Information local to one thread in the parallel quicksort run. */
template<typename RandomAccessIterator>
  struct QSBThreadLocal
  {
    typedef std::iterator_traits<RandomAccessIterator> traits_type;
    typedef typename traits_type::difference_type difference_type;

    /** @brief Continuous part of the sequence, described by an
    iterator pair. */
    typedef std::pair<RandomAccessIterator, RandomAccessIterator> Piece;

    /** @brief Initial piece to work on. */
    Piece initial;

    /** @brief Work-stealing queue. */
    RestrictedBoundedConcurrentQueue<Piece> leftover_parts;

    /** @brief Number of threads involved in this algorithm. */
    thread_index_t num_threads;

    /** @brief Pointer to a counter of elements left over to sort. */
    volatile difference_type* elements_leftover;

    /** @brief The complete sequence to sort. */
    Piece global;

    /** @brief Constructor.
     *  @param queue_size Size of the work-stealing queue. */
    QSBThreadLocal(int queue_size) : leftover_parts(queue_size) { }
  };

/** @brief Balanced quicksort divide step.
  *  @param begin Begin iterator of subsequence.
  *  @param end End iterator of subsequence.
  *  @param comp Comparator.
  *  @param num_threads Number of threads that are allowed to work on
  *  this part.
  *  @pre @c (end-begin)>=1 */
template<typename RandomAccessIterator, typename Comparator>
  typename std::iterator_traits<RandomAccessIterator>::difference_type
  qsb_divide(RandomAccessIterator begin, RandomAccessIterator end,
             Comparator comp, thread_index_t num_threads)
  {
    _GLIBCXX_PARALLEL_ASSERT(num_threads > 0);

    typedef std::iterator_traits<RandomAccessIterator> traits_type;
    typedef typename traits_type::value_type value_type;
    typedef typename traits_type::difference_type difference_type;

    RandomAccessIterator pivot_pos =
      median_of_three_iterators(begin, begin + (end - begin) / 2,
				end  - 1, comp);

#if defined(_GLIBCXX_ASSERTIONS)
    // Must be in between somewhere.
    difference_type n = end - begin;

    _GLIBCXX_PARALLEL_ASSERT(
           (!comp(*pivot_pos, *begin) && !comp(*(begin + n / 2), *pivot_pos))
        || (!comp(*pivot_pos, *begin) && !comp(*(end - 1), *pivot_pos))
        || (!comp(*pivot_pos, *(begin + n / 2)) && !comp(*begin, *pivot_pos))
        || (!comp(*pivot_pos, *(begin + n / 2)) && !comp(*(end - 1), *pivot_pos))
        || (!comp(*pivot_pos, *(end - 1)) && !comp(*begin, *pivot_pos))
        || (!comp(*pivot_pos, *(end - 1)) && !comp(*(begin + n / 2), *pivot_pos)));
#endif

    // Swap pivot value to end.
    if (pivot_pos != (end - 1))
      std::swap(*pivot_pos, *(end - 1));
    pivot_pos = end - 1;

    __gnu_parallel::binder2nd<Comparator, value_type, value_type, bool>
        pred(comp, *pivot_pos);

    // Divide, returning end - begin - 1 in the worst case.
    difference_type split_pos = parallel_partition(
        begin, end - 1, pred, num_threads);

    // Swap back pivot to middle.
    std::swap(*(begin + split_pos), *pivot_pos);
    pivot_pos = begin + split_pos;

#if _GLIBCXX_ASSERTIONS
    RandomAccessIterator r;
    for (r = begin; r != pivot_pos; ++r)
      _GLIBCXX_PARALLEL_ASSERT(comp(*r, *pivot_pos));
    for (; r != end; ++r)
      _GLIBCXX_PARALLEL_ASSERT(!comp(*r, *pivot_pos));
#endif

    return split_pos;
  }

/** @brief Quicksort conquer step.
  *  @param tls Array of thread-local storages.
  *  @param begin Begin iterator of subsequence.
  *  @param end End iterator of subsequence.
  *  @param comp Comparator.
  *  @param iam Number of the thread processing this function.
  *  @param num_threads
  *          Number of threads that are allowed to work on this part. */
template<typename RandomAccessIterator, typename Comparator>
  void
  qsb_conquer(QSBThreadLocal<RandomAccessIterator>** tls,
              RandomAccessIterator begin, RandomAccessIterator end,
              Comparator comp,
              thread_index_t iam, thread_index_t num_threads,
              bool parent_wait)
  {
    typedef std::iterator_traits<RandomAccessIterator> traits_type;
    typedef typename traits_type::value_type value_type;
    typedef typename traits_type::difference_type difference_type;

    difference_type n = end - begin;

    if (num_threads <= 1 || n <= 1)
      {
        tls[iam]->initial.first  = begin;
        tls[iam]->initial.second = end;

        qsb_local_sort_with_helping(tls, comp, iam, parent_wait);

        return;
      }

    // Divide step.
    difference_type split_pos = qsb_divide(begin, end, comp, num_threads);

#if _GLIBCXX_ASSERTIONS
    _GLIBCXX_PARALLEL_ASSERT(0 <= split_pos && split_pos < (end - begin));
#endif

    thread_index_t num_threads_leftside =
        std::max<thread_index_t>(1, std::min<thread_index_t>(
                          num_threads - 1, split_pos * num_threads / n));

#   pragma omp atomic
    *tls[iam]->elements_leftover -= (difference_type)1;

    // Conquer step.
#   pragma omp parallel num_threads(2)
    {
      bool wait;
      if(omp_get_num_threads() < 2)
        wait = false;
      else
        wait = parent_wait;

#     pragma omp sections
        {
#         pragma omp section
            {
              qsb_conquer(tls, begin, begin + split_pos, comp,
                          iam,
                          num_threads_leftside,
                          wait);
              wait = parent_wait;
            }
          // The pivot_pos is left in place, to ensure termination.
#         pragma omp section
            {
              qsb_conquer(tls, begin + split_pos + 1, end, comp,
                          iam + num_threads_leftside,
                          num_threads - num_threads_leftside,
                          wait);
              wait = parent_wait;
            }
        }
    }
  }

/**
  *  @brief Quicksort step doing load-balanced local sort.
  *  @param tls Array of thread-local storages.
  *  @param comp Comparator.
  *  @param iam Number of the thread processing this function.
  */
template<typename RandomAccessIterator, typename Comparator>
  void
  qsb_local_sort_with_helping(QSBThreadLocal<RandomAccessIterator>** tls,
                              Comparator& comp, int iam, bool wait)
  {
    typedef std::iterator_traits<RandomAccessIterator> traits_type;
    typedef typename traits_type::value_type value_type;
    typedef typename traits_type::difference_type difference_type;
    typedef std::pair<RandomAccessIterator, RandomAccessIterator> Piece;

    QSBThreadLocal<RandomAccessIterator>& tl = *tls[iam];

    difference_type base_case_n =
        _Settings::get().sort_qsb_base_case_maximal_n;
    if (base_case_n < 2)
      base_case_n = 2;
    thread_index_t num_threads = tl.num_threads;

    // Every thread has its own random number generator.
    random_number rng(iam + 1);

    Piece current = tl.initial;

    difference_type elements_done = 0;
#if _GLIBCXX_ASSERTIONS
    difference_type total_elements_done = 0;
#endif

    for (;;)
      {
        // Invariant: current must be a valid (maybe empty) range.
        RandomAccessIterator begin = current.first, end = current.second;
        difference_type n = end - begin;

        if (n > base_case_n)
          {
            // Divide.
            RandomAccessIterator pivot_pos = begin +  rng(n);

            // Swap pivot_pos value to end.
            if (pivot_pos != (end - 1))
              std::swap(*pivot_pos, *(end - 1));
            pivot_pos = end - 1;

            __gnu_parallel::binder2nd
                <Comparator, value_type, value_type, bool>
                pred(comp, *pivot_pos);

            // Divide, leave pivot unchanged in last place.
            RandomAccessIterator split_pos1, split_pos2;
            split_pos1 = __gnu_sequential::partition(begin, end - 1, pred);

            // Left side: < pivot_pos; right side: >= pivot_pos.
#if _GLIBCXX_ASSERTIONS
            _GLIBCXX_PARALLEL_ASSERT(begin <= split_pos1 && split_pos1 < end);
#endif
            // Swap pivot back to middle.
            if (split_pos1 != pivot_pos)
              std::swap(*split_pos1, *pivot_pos);
            pivot_pos = split_pos1;

            // In case all elements are equal, split_pos1 == 0.
            if ((split_pos1 + 1 - begin) < (n >> 7)
            || (end - split_pos1) < (n >> 7))
              {
                // Very unequal split, one part smaller than one 128th
                // elements not strictly larger than the pivot.
                __gnu_parallel::unary_negate<__gnu_parallel::binder1st
		  <Comparator, value_type, value_type, bool>, value_type>
		  pred(__gnu_parallel::binder1st
		       <Comparator, value_type, value_type, bool>(comp,
								  *pivot_pos));

                // Find other end of pivot-equal range.
                split_pos2 = __gnu_sequential::partition(split_pos1 + 1,
							 end, pred);
              }
            else
              // Only skip the pivot.
              split_pos2 = split_pos1 + 1;

            // Elements equal to pivot are done.
            elements_done += (split_pos2 - split_pos1);
#if _GLIBCXX_ASSERTIONS
            total_elements_done += (split_pos2 - split_pos1);
#endif
            // Always push larger part onto stack.
            if (((split_pos1 + 1) - begin) < (end - (split_pos2)))
              {
                // Right side larger.
                if ((split_pos2) != end)
                  tl.leftover_parts.push_front(std::make_pair(split_pos2,
							      end));

                //current.first = begin;	//already set anyway
                current.second = split_pos1;
                continue;
              }
            else
              {
                // Left side larger.
                if (begin != split_pos1)
                  tl.leftover_parts.push_front(std::make_pair(begin,
							      split_pos1));

                current.first = split_pos2;
                //current.second = end;	//already set anyway
                continue;
              }
          }
        else
          {
            __gnu_sequential::sort(begin, end, comp);
            elements_done += n;
#if _GLIBCXX_ASSERTIONS
            total_elements_done += n;
#endif

            // Prefer own stack, small pieces.
            if (tl.leftover_parts.pop_front(current))
              continue;

#           pragma omp atomic
            *tl.elements_leftover -= elements_done;

            elements_done = 0;

#if _GLIBCXX_ASSERTIONS
            double search_start = omp_get_wtime();
#endif

            // Look for new work.
            bool successfully_stolen = false;
            while (wait && *tl.elements_leftover > 0 && !successfully_stolen
#if _GLIBCXX_ASSERTIONS
              // Possible dead-lock.
              && (omp_get_wtime() < (search_start + 1.0))
#endif
              )
              {
                thread_index_t victim;
                victim = rng(num_threads);

                // Large pieces.
                successfully_stolen = (victim != iam)
                    && tls[victim]->leftover_parts.pop_back(current);
                if (!successfully_stolen)
                  yield();
#if !defined(__ICC) && !defined(__ECC)
#               pragma omp flush
#endif
              }

#if _GLIBCXX_ASSERTIONS
            if (omp_get_wtime() >= (search_start + 1.0))
              {
                sleep(1);
                _GLIBCXX_PARALLEL_ASSERT(omp_get_wtime()
					 < (search_start + 1.0));
              }
#endif
            if (!successfully_stolen)
              {
#if _GLIBCXX_ASSERTIONS
                _GLIBCXX_PARALLEL_ASSERT(*tl.elements_leftover == 0);
#endif
                return;
              }
          }
      }
  }

/** @brief Top-level quicksort routine.
  *  @param begin Begin iterator of sequence.
  *  @param end End iterator of sequence.
  *  @param comp Comparator.
  *  @param num_threads Number of threads that are allowed to work on
  *  this part.
  */
template<typename RandomAccessIterator, typename Comparator>
  void
  parallel_sort_qsb(RandomAccessIterator begin, RandomAccessIterator end,
                    Comparator comp,
                    thread_index_t num_threads)
  {
    _GLIBCXX_CALL(end - begin)

    typedef std::iterator_traits<RandomAccessIterator> traits_type;
    typedef typename traits_type::value_type value_type;
    typedef typename traits_type::difference_type difference_type;
    typedef std::pair<RandomAccessIterator, RandomAccessIterator> Piece;

    typedef QSBThreadLocal<RandomAccessIterator> tls_type;

    difference_type n = end - begin;

    if (n <= 1)
      return;

    // At least one element per processor.
    if (num_threads > n)
      num_threads = static_cast<thread_index_t>(n);

    // Initialize thread local storage
    tls_type** tls = new tls_type*[num_threads];
    difference_type queue_size = num_threads * (thread_index_t)(log2(n) + 1);
    for (thread_index_t t = 0; t < num_threads; ++t)
      tls[t] = new QSBThreadLocal<RandomAccessIterator>(queue_size);

    // There can never be more than ceil(log2(n)) ranges on the stack, because
    // 1. Only one processor pushes onto the stack
    // 2. The largest range has at most length n
    // 3. Each range is larger than half of the range remaining
    volatile difference_type elements_leftover = n;
    for (int i = 0; i < num_threads; ++i)
      {
        tls[i]->elements_leftover = &elements_leftover;
        tls[i]->num_threads = num_threads;
        tls[i]->global = std::make_pair(begin, end);

        // Just in case nothing is left to assign.
        tls[i]->initial = std::make_pair(end, end);
      }

    // Main recursion call.
    qsb_conquer(tls, begin, begin + n, comp, 0, num_threads, true);

#if _GLIBCXX_ASSERTIONS
    // All stack must be empty.
    Piece dummy;
    for (int i = 1; i < num_threads; ++i)
      _GLIBCXX_PARALLEL_ASSERT(!tls[i]->leftover_parts.pop_back(dummy));
#endif

    for (int i = 0; i < num_threads; ++i)
      delete tls[i];
    delete[] tls;
  }
} // namespace __gnu_parallel

#endif /* _GLIBCXX_PARALLEL_BALANCED_QUICKSORT_H */