1// -*- C++ -*-
2
3// Copyright (C) 2007-2015 Free Software Foundation, Inc.
4//
5// This file is part of the GNU ISO C++ Library.  This library is free
6// software; you can redistribute it and/or modify it under the terms
7// of the GNU General Public License as published by the Free Software
8// Foundation; either version 3, or (at your option) any later
9// version.
10
11// This library is distributed in the hope that it will be useful, but
12// WITHOUT ANY WARRANTY; without even the implied warranty of
13// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14// General Public License for more details.
15
16// Under Section 7 of GPL version 3, you are granted additional
17// permissions described in the GCC Runtime Library Exception, version
18// 3.1, as published by the Free Software Foundation.
19
20// You should have received a copy of the GNU General Public License and
21// a copy of the GCC Runtime Library Exception along with this program;
22// see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
23// <http://www.gnu.org/licenses/>.
24
25/** @file parallel/partition.h
26 *  @brief Parallel implementation of std::partition(),
27 *  std::nth_element(), and std::partial_sort().
28 *  This file is a GNU parallel extension to the Standard C++ Library.
29 */
30
31// Written by Johannes Singler and Felix Putze.
32
33#ifndef _GLIBCXX_PARALLEL_PARTITION_H
34#define _GLIBCXX_PARALLEL_PARTITION_H 1
35
36#include <parallel/basic_iterator.h>
37#include <parallel/sort.h>
38#include <parallel/random_number.h>
39#include <bits/stl_algo.h>
40#include <parallel/parallel.h>
41
42/** @brief Decide whether to declare certain variables volatile. */
43#define _GLIBCXX_VOLATILE volatile
44
45namespace __gnu_parallel
46{
47  /** @brief Parallel implementation of std::partition.
48    *  @param __begin Begin iterator of input sequence to split.
49    *  @param __end End iterator of input sequence to split.
50    *  @param __pred Partition predicate, possibly including some kind
51    *         of pivot.
52    *  @param __num_threads Maximum number of threads to use for this task.
53    *  @return Number of elements not fulfilling the predicate. */
54  template<typename _RAIter, typename _Predicate>
55    typename std::iterator_traits<_RAIter>::difference_type
56    __parallel_partition(_RAIter __begin, _RAIter __end,
57			 _Predicate __pred, _ThreadIndex __num_threads)
58    {
59      typedef std::iterator_traits<_RAIter> _TraitsType;
60      typedef typename _TraitsType::value_type _ValueType;
61      typedef typename _TraitsType::difference_type _DifferenceType;
62
63      _DifferenceType __n = __end - __begin;
64
65      _GLIBCXX_CALL(__n)
66
67      const _Settings& __s = _Settings::get();
68
69      // shared
70      _GLIBCXX_VOLATILE _DifferenceType __left = 0, __right = __n - 1,
71                                        __dist = __n,
72                                        __leftover_left, __leftover_right,
73                                        __leftnew, __rightnew;
74
75      // just 0 or 1, but int to allow atomic operations
76      int* __reserved_left = 0, * __reserved_right = 0;
77
78      _DifferenceType __chunk_size = __s.partition_chunk_size;
79
80      //at least two chunks per thread
81      if (__dist >= 2 * __num_threads * __chunk_size)
82#       pragma omp parallel num_threads(__num_threads)
83	{
84#         pragma omp single
85	  {
86	    __num_threads = omp_get_num_threads();
87	    __reserved_left = new int[__num_threads];
88	    __reserved_right = new int[__num_threads];
89
90	    if (__s.partition_chunk_share > 0.0)
91	      __chunk_size = std::max<_DifferenceType>
92		(__s.partition_chunk_size, (double)__n
93		 * __s.partition_chunk_share / (double)__num_threads);
94	    else
95	      __chunk_size = __s.partition_chunk_size;
96	  }
97
98	  while (__dist >= 2 * __num_threads * __chunk_size)
99	    {
100#             pragma omp single
101	      {
102		_DifferenceType __num_chunks = __dist / __chunk_size;
103
104		for (_ThreadIndex __r = 0; __r < __num_threads; ++__r)
105		  {
106		    __reserved_left [__r] = 0; // false
107		    __reserved_right[__r] = 0; // false
108		  }
109		__leftover_left = 0;
110		__leftover_right = 0;
111	      } //implicit barrier
112
113	      // Private.
114	      _DifferenceType __thread_left, __thread_left_border,
115		              __thread_right, __thread_right_border;
116
117	      __thread_left = __left + 1;
118	      // Just to satisfy the condition below.
119	      __thread_left_border = __thread_left - 1;
120
121	      __thread_right = __n - 1;
122             // Just to satisfy the condition below.
123	      __thread_right_border = __thread_right + 1;
124
125	      bool __iam_finished = false;
126	      while (!__iam_finished)
127		{
128		  if (__thread_left > __thread_left_border)
129		    {
130                      _DifferenceType __former_dist =
131                              __fetch_and_add(&__dist, -__chunk_size);
132                      if (__former_dist < __chunk_size)
133                        {
134                          __fetch_and_add(&__dist, __chunk_size);
135                          __iam_finished = true;
136                          break;
137                        }
138                      else
139                        {
140                          __thread_left =
141                                  __fetch_and_add(&__left, __chunk_size);
142                          __thread_left_border =
143                                  __thread_left + (__chunk_size - 1);
144                        }
145		    }
146
147		  if (__thread_right < __thread_right_border)
148		    {
149                      _DifferenceType __former_dist =
150                              __fetch_and_add(&__dist, -__chunk_size);
151                      if (__former_dist < __chunk_size)
152                        {
153                          __fetch_and_add(&__dist, __chunk_size);
154                          __iam_finished = true;
155                          break;
156                        }
157                      else
158                        {
159                          __thread_right =
160                                  __fetch_and_add(&__right, -__chunk_size);
161                          __thread_right_border =
162                                  __thread_right - (__chunk_size - 1);
163                        }
164		    }
165
166		  // Swap as usual.
167		  while (__thread_left < __thread_right)
168		    {
169		      while (__pred(__begin[__thread_left])
170			     && __thread_left <= __thread_left_border)
171			++__thread_left;
172		      while (!__pred(__begin[__thread_right])
173			     && __thread_right >= __thread_right_border)
174			--__thread_right;
175
176		      if (__thread_left > __thread_left_border
177			  || __thread_right < __thread_right_border)
178			// Fetch new chunk(__s).
179			break;
180
181		      std::iter_swap(__begin + __thread_left,
182                             __begin + __thread_right);
183		      ++__thread_left;
184		      --__thread_right;
185		    }
186		}
187
188	      // Now swap the leftover chunks to the right places.
189	      if (__thread_left <= __thread_left_border)
190#               pragma omp atomic
191		++__leftover_left;
192	      if (__thread_right >= __thread_right_border)
193#               pragma omp atomic
194		++__leftover_right;
195
196#             pragma omp barrier
197
198              _DifferenceType
199                    __leftold = __left,
200                    __leftnew = __left - __leftover_left * __chunk_size,
201                    __rightold = __right,
202                    __rightnew = __right + __leftover_right * __chunk_size;
203
204	      // <=> __thread_left_border + (__chunk_size - 1) >= __leftnew
205	      if (__thread_left <= __thread_left_border
206		  && __thread_left_border >= __leftnew)
207		{
208		  // Chunk already in place, reserve spot.
209		__reserved_left[(__left - (__thread_left_border + 1))
210				/ __chunk_size] = 1;
211		}
212
213	      // <=> __thread_right_border - (__chunk_size - 1) <= __rightnew
214	      if (__thread_right >= __thread_right_border
215		  && __thread_right_border <= __rightnew)
216		{
217		  // Chunk already in place, reserve spot.
218		  __reserved_right[((__thread_right_border - 1) - __right)
219				   / __chunk_size] = 1;
220		}
221
222#             pragma omp barrier
223
224	      if (__thread_left <= __thread_left_border
225		  && __thread_left_border < __leftnew)
226		{
227		  // Find spot and swap.
228		  _DifferenceType __swapstart = -1;
229                  for (int __r = 0; __r < __leftover_left; ++__r)
230                    if (__reserved_left[__r] == 0
231                        && __compare_and_swap(&(__reserved_left[__r]), 0, 1))
232                      {
233                        __swapstart = __leftold - (__r + 1) * __chunk_size;
234                        break;
235                      }
236
237#if _GLIBCXX_ASSERTIONS
238		  _GLIBCXX_PARALLEL_ASSERT(__swapstart != -1);
239#endif
240
241		  std::swap_ranges(__begin + __thread_left_border
242				   - (__chunk_size - 1),
243				   __begin + __thread_left_border + 1,
244				   __begin + __swapstart);
245		}
246
247	      if (__thread_right >= __thread_right_border
248		  && __thread_right_border > __rightnew)
249		{
250		  // Find spot and swap
251		  _DifferenceType __swapstart = -1;
252                  for (int __r = 0; __r < __leftover_right; ++__r)
253                    if (__reserved_right[__r] == 0
254                        && __compare_and_swap(&(__reserved_right[__r]), 0, 1))
255                      {
256                        __swapstart = __rightold + __r * __chunk_size + 1;
257                        break;
258                      }
259
260#if _GLIBCXX_ASSERTIONS
261		  _GLIBCXX_PARALLEL_ASSERT(__swapstart != -1);
262#endif
263
264		  std::swap_ranges(__begin + __thread_right_border,
265				   __begin + __thread_right_border
266				   + __chunk_size, __begin + __swapstart);
267	      }
268#if _GLIBCXX_ASSERTIONS
269#             pragma omp barrier
270
271#             pragma omp single
272	      {
273		for (_DifferenceType __r = 0; __r < __leftover_left; ++__r)
274		  _GLIBCXX_PARALLEL_ASSERT(__reserved_left[__r] == 1);
275		for (_DifferenceType __r = 0; __r < __leftover_right; ++__r)
276		  _GLIBCXX_PARALLEL_ASSERT(__reserved_right[__r] == 1);
277	      }
278#endif
279
280	      __left = __leftnew;
281	      __right = __rightnew;
282              __dist = __right - __left + 1;
283	    }
284
285#           pragma omp flush(__left, __right)
286	} // end "recursion" //parallel
287
288        _DifferenceType __final_left = __left, __final_right = __right;
289
290	while (__final_left < __final_right)
291	  {
292	    // Go right until key is geq than pivot.
293	    while (__pred(__begin[__final_left])
294		   && __final_left < __final_right)
295	      ++__final_left;
296
297	    // Go left until key is less than pivot.
298	    while (!__pred(__begin[__final_right])
299		   && __final_left < __final_right)
300	      --__final_right;
301
302	    if (__final_left == __final_right)
303	      break;
304	    std::iter_swap(__begin + __final_left, __begin + __final_right);
305	    ++__final_left;
306	    --__final_right;
307	  }
308
309	// All elements on the left side are < piv, all elements on the
310	// right are >= piv
311	delete[] __reserved_left;
312	delete[] __reserved_right;
313
314	// Element "between" __final_left and __final_right might not have
315	// been regarded yet
316	if (__final_left < __n && !__pred(__begin[__final_left]))
317	  // Really swapped.
318	  return __final_left;
319	else
320	  return __final_left + 1;
321    }
322
323  /**
324    *  @brief Parallel implementation of std::nth_element().
325    *  @param __begin Begin iterator of input sequence.
326    *  @param __nth _Iterator of element that must be in position afterwards.
327    *  @param __end End iterator of input sequence.
328    *  @param __comp Comparator.
329    */
330  template<typename _RAIter, typename _Compare>
331    void
332    __parallel_nth_element(_RAIter __begin, _RAIter __nth,
333			   _RAIter __end, _Compare __comp)
334    {
335      typedef std::iterator_traits<_RAIter> _TraitsType;
336      typedef typename _TraitsType::value_type _ValueType;
337      typedef typename _TraitsType::difference_type _DifferenceType;
338
339      _GLIBCXX_CALL(__end - __begin)
340
341      _RAIter __split;
342      _RandomNumber __rng;
343
344      const _Settings& __s = _Settings::get();
345      _DifferenceType __minimum_length = std::max<_DifferenceType>(2,
346        std::max(__s.nth_element_minimal_n, __s.partition_minimal_n));
347
348      // Break if input range to small.
349      while (static_cast<_SequenceIndex>(__end - __begin) >= __minimum_length)
350	{
351          _DifferenceType __n = __end - __begin;
352
353          _RAIter __pivot_pos = __begin + __rng(__n);
354
355          // Swap __pivot_pos value to end.
356          if (__pivot_pos != (__end - 1))
357            std::iter_swap(__pivot_pos, __end - 1);
358          __pivot_pos = __end - 1;
359
360          // _Compare must have first_value_type, second_value_type,
361          // result_type
362          // _Compare ==
363          // __gnu_parallel::_Lexicographic<S, int,
364	  //                                __gnu_parallel::_Less<S, S> >
365          // __pivot_pos == std::pair<S, int>*
366          __gnu_parallel::__binder2nd<_Compare, _ValueType, _ValueType, bool>
367            __pred(__comp, *__pivot_pos);
368
369          // Divide, leave pivot unchanged in last place.
370          _RAIter __split_pos1, __split_pos2;
371          __split_pos1 = __begin + __parallel_partition(__begin, __end - 1,
372							__pred,
373							__get_max_threads());
374
375          // Left side: < __pivot_pos; __right side: >= __pivot_pos
376
377          // Swap pivot back to middle.
378          if (__split_pos1 != __pivot_pos)
379            std::iter_swap(__split_pos1, __pivot_pos);
380          __pivot_pos = __split_pos1;
381
382          // In case all elements are equal, __split_pos1 == 0
383          if ((__split_pos1 + 1 - __begin) < (__n >> 7)
384              || (__end - __split_pos1) < (__n >> 7))
385            {
386              // Very unequal split, one part smaller than one 128th
387              // elements not strictly larger than the pivot.
388              __gnu_parallel::__unary_negate<__gnu_parallel::
389        	__binder1st<_Compare, _ValueType,
390		            _ValueType, bool>, _ValueType>
391        	__pred(__gnu_parallel::__binder1st<_Compare, _ValueType,
392		       _ValueType, bool>(__comp, *__pivot_pos));
393
394              // Find other end of pivot-equal range.
395              __split_pos2 = __gnu_sequential::partition(__split_pos1 + 1,
396							 __end, __pred);
397            }
398          else
399            // Only skip the pivot.
400            __split_pos2 = __split_pos1 + 1;
401
402          // Compare iterators.
403          if (__split_pos2 <= __nth)
404            __begin = __split_pos2;
405          else if (__nth < __split_pos1)
406            __end = __split_pos1;
407          else
408            break;
409	}
410
411      // Only at most _Settings::partition_minimal_n __elements __left.
412      __gnu_sequential::nth_element(__begin, __nth, __end, __comp);
413    }
414
415  /** @brief Parallel implementation of std::partial_sort().
416  *  @param __begin Begin iterator of input sequence.
417  *  @param __middle Sort until this position.
418  *  @param __end End iterator of input sequence.
419  *  @param __comp Comparator. */
420  template<typename _RAIter, typename _Compare>
421    void
422    __parallel_partial_sort(_RAIter __begin,
423			    _RAIter __middle,
424			    _RAIter __end, _Compare __comp)
425    {
426      __parallel_nth_element(__begin, __middle, __end, __comp);
427      std::sort(__begin, __middle, __comp);
428    }
429
430} //namespace __gnu_parallel
431
432#undef _GLIBCXX_VOLATILE
433
434#endif /* _GLIBCXX_PARALLEL_PARTITION_H */
435