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/random_shuffle.h
26 *  @brief Parallel implementation of std::random_shuffle().
27 *  This file is a GNU parallel extension to the Standard C++ Library.
28 */
29
30// Written by Johannes Singler.
31
32#ifndef _GLIBCXX_PARALLEL_RANDOM_SHUFFLE_H
33#define _GLIBCXX_PARALLEL_RANDOM_SHUFFLE_H 1
34
35#include <limits>
36#include <bits/stl_numeric.h>
37#include <parallel/parallel.h>
38#include <parallel/random_number.h>
39
40namespace __gnu_parallel
41{
42  /** @brief Type to hold the index of a bin.
43    *
44    *  Since many variables of this type are allocated, it should be
45    *  chosen as small as possible.
46    */
47  typedef unsigned short _BinIndex;
48
49  /** @brief Data known to every thread participating in
50      __gnu_parallel::__parallel_random_shuffle(). */
51  template<typename _RAIter>
52    struct _DRandomShufflingGlobalData
53    {
54      typedef std::iterator_traits<_RAIter> _TraitsType;
55      typedef typename _TraitsType::value_type _ValueType;
56      typedef typename _TraitsType::difference_type _DifferenceType;
57
58      /** @brief Begin iterator of the __source. */
59      _RAIter& _M_source;
60
61      /** @brief Temporary arrays for each thread. */
62      _ValueType** _M_temporaries;
63
64      /** @brief Two-dimensional array to hold the thread-bin distribution.
65       *
66       *  Dimensions (_M_num_threads + 1) __x (_M_num_bins + 1). */
67      _DifferenceType** _M_dist;
68
69      /** @brief Start indexes of the threads' __chunks. */
70      _DifferenceType* _M_starts;
71
72      /** @brief Number of the thread that will further process the
73          corresponding bin. */
74      _ThreadIndex* _M_bin_proc;
75
76      /** @brief Number of bins to distribute to. */
77      int _M_num_bins;
78
79      /** @brief Number of bits needed to address the bins. */
80      int _M_num_bits;
81
82      /** @brief Constructor. */
83      _DRandomShufflingGlobalData(_RAIter& __source)
84      : _M_source(__source) { }
85    };
86
87  /** @brief Local data for a thread participating in
88      __gnu_parallel::__parallel_random_shuffle().
89    */
90  template<typename _RAIter, typename _RandomNumberGenerator>
91    struct _DRSSorterPU
92    {
93      /** @brief Number of threads participating in total. */
94      int _M_num_threads;
95
96      /** @brief Begin index for bins taken care of by this thread. */
97      _BinIndex _M_bins_begin;
98
99      /** @brief End index for bins taken care of by this thread. */
100      _BinIndex __bins_end;
101
102      /** @brief Random _M_seed for this thread. */
103      uint32_t _M_seed;
104
105      /** @brief Pointer to global data. */
106      _DRandomShufflingGlobalData<_RAIter>* _M_sd;
107    };
108
109  /** @brief Generate a random number in @c [0,2^__logp).
110    *  @param __logp Logarithm (basis 2) of the upper range __bound.
111    *  @param __rng Random number generator to use.
112    */
113  template<typename _RandomNumberGenerator>
114    inline int
115    __random_number_pow2(int __logp, _RandomNumberGenerator& __rng)
116    { return __rng.__genrand_bits(__logp); }
117
118  /** @brief Random shuffle code executed by each thread.
119    *  @param __pus Array of thread-local data records. */
120  template<typename _RAIter, typename _RandomNumberGenerator>
121    void
122    __parallel_random_shuffle_drs_pu(_DRSSorterPU<_RAIter,
123				     _RandomNumberGenerator>* __pus)
124    {
125      typedef std::iterator_traits<_RAIter> _TraitsType;
126      typedef typename _TraitsType::value_type _ValueType;
127      typedef typename _TraitsType::difference_type _DifferenceType;
128
129      _ThreadIndex __iam = omp_get_thread_num();
130      _DRSSorterPU<_RAIter, _RandomNumberGenerator>* __d = &__pus[__iam];
131      _DRandomShufflingGlobalData<_RAIter>* __sd = __d->_M_sd;
132
133      // Indexing: _M_dist[bin][processor]
134      _DifferenceType __length = (__sd->_M_starts[__iam + 1]
135				  - __sd->_M_starts[__iam]);
136      _BinIndex* __oracles = new _BinIndex[__length];
137      _DifferenceType* __dist = new _DifferenceType[__sd->_M_num_bins + 1];
138      _BinIndex* __bin_proc = new _BinIndex[__sd->_M_num_bins];
139      _ValueType** __temporaries = new _ValueType*[__d->_M_num_threads];
140
141      // Compute oracles and count appearances.
142      for (_BinIndex __b = 0; __b < __sd->_M_num_bins + 1; ++__b)
143	__dist[__b] = 0;
144      int __num_bits = __sd->_M_num_bits;
145
146      _RandomNumber __rng(__d->_M_seed);
147
148      // First main loop.
149      for (_DifferenceType __i = 0; __i < __length; ++__i)
150	{
151          _BinIndex __oracle = __random_number_pow2(__num_bits, __rng);
152          __oracles[__i] = __oracle;
153
154          // To allow prefix (partial) sum.
155          ++(__dist[__oracle + 1]);
156	}
157
158      for (_BinIndex __b = 0; __b < __sd->_M_num_bins + 1; ++__b)
159	__sd->_M_dist[__b][__iam + 1] = __dist[__b];
160
161#     pragma omp barrier
162
163#     pragma omp single
164      {
165	// Sum up bins, __sd->_M_dist[__s + 1][__d->_M_num_threads] now
166	// contains the total number of items in bin __s
167	for (_BinIndex __s = 0; __s < __sd->_M_num_bins; ++__s)
168          __gnu_sequential::partial_sum(__sd->_M_dist[__s + 1],
169					__sd->_M_dist[__s + 1]
170					+ __d->_M_num_threads + 1,
171					__sd->_M_dist[__s + 1]);
172      }
173
174#     pragma omp barrier
175
176      _SequenceIndex __offset = 0, __global_offset = 0;
177      for (_BinIndex __s = 0; __s < __d->_M_bins_begin; ++__s)
178	__global_offset += __sd->_M_dist[__s + 1][__d->_M_num_threads];
179
180#     pragma omp barrier
181
182      for (_BinIndex __s = __d->_M_bins_begin; __s < __d->__bins_end; ++__s)
183	{
184          for (int __t = 0; __t < __d->_M_num_threads + 1; ++__t)
185            __sd->_M_dist[__s + 1][__t] += __offset;
186          __offset = __sd->_M_dist[__s + 1][__d->_M_num_threads];
187	}
188
189      __sd->_M_temporaries[__iam] = static_cast<_ValueType*>
190	(::operator new(sizeof(_ValueType) * __offset));
191
192#     pragma omp barrier
193
194      // Draw local copies to avoid false sharing.
195      for (_BinIndex __b = 0; __b < __sd->_M_num_bins + 1; ++__b)
196	__dist[__b] = __sd->_M_dist[__b][__iam];
197      for (_BinIndex __b = 0; __b < __sd->_M_num_bins; ++__b)
198	__bin_proc[__b] = __sd->_M_bin_proc[__b];
199      for (_ThreadIndex __t = 0; __t < __d->_M_num_threads; ++__t)
200	__temporaries[__t] = __sd->_M_temporaries[__t];
201
202      _RAIter __source = __sd->_M_source;
203      _DifferenceType __start = __sd->_M_starts[__iam];
204
205      // Distribute according to oracles, second main loop.
206      for (_DifferenceType __i = 0; __i < __length; ++__i)
207	{
208          _BinIndex __target_bin = __oracles[__i];
209          _ThreadIndex __target_p = __bin_proc[__target_bin];
210
211          // Last column [__d->_M_num_threads] stays unchanged.
212	  ::new(&(__temporaries[__target_p][__dist[__target_bin + 1]++]))
213              _ValueType(*(__source + __i + __start));
214	}
215
216      delete[] __oracles;
217      delete[] __dist;
218      delete[] __bin_proc;
219      delete[] __temporaries;
220
221#     pragma omp barrier
222
223      // Shuffle bins internally.
224      for (_BinIndex __b = __d->_M_bins_begin; __b < __d->__bins_end; ++__b)
225	{
226          _ValueType* __begin =
227	    (__sd->_M_temporaries[__iam]
228	     + (__b == __d->_M_bins_begin
229		? 0 : __sd->_M_dist[__b][__d->_M_num_threads])),
230	    *__end = (__sd->_M_temporaries[__iam]
231		      + __sd->_M_dist[__b + 1][__d->_M_num_threads]);
232
233          __sequential_random_shuffle(__begin, __end, __rng);
234          std::copy(__begin, __end, __sd->_M_source + __global_offset
235		    + (__b == __d->_M_bins_begin
236		       ? 0 : __sd->_M_dist[__b][__d->_M_num_threads]));
237	}
238
239      for (_SequenceIndex __i = 0; __i < __offset; ++__i)
240	__sd->_M_temporaries[__iam][__i].~_ValueType();
241      ::operator delete(__sd->_M_temporaries[__iam]);
242    }
243
244  /** @brief Round up to the next greater power of 2.
245    *  @param __x _Integer to round up */
246  template<typename _Tp>
247    _Tp
248    __round_up_to_pow2(_Tp __x)
249    {
250      if (__x <= 1)
251	return 1;
252      else
253	return (_Tp)1 << (__rd_log2(__x - 1) + 1);
254    }
255
256  /** @brief Main parallel random shuffle step.
257    *  @param __begin Begin iterator of sequence.
258    *  @param __end End iterator of sequence.
259    *  @param __n Length of sequence.
260    *  @param __num_threads Number of threads to use.
261    *  @param __rng Random number generator to use.
262    */
263  template<typename _RAIter, typename _RandomNumberGenerator>
264    void
265    __parallel_random_shuffle_drs(_RAIter __begin, _RAIter __end,
266				  typename std::iterator_traits
267				  <_RAIter>::difference_type __n,
268				  _ThreadIndex __num_threads,
269				  _RandomNumberGenerator& __rng)
270    {
271      typedef std::iterator_traits<_RAIter> _TraitsType;
272      typedef typename _TraitsType::value_type _ValueType;
273      typedef typename _TraitsType::difference_type _DifferenceType;
274
275      _GLIBCXX_CALL(__n)
276
277      const _Settings& __s = _Settings::get();
278
279      if (__num_threads > __n)
280	__num_threads = static_cast<_ThreadIndex>(__n);
281
282      _BinIndex __num_bins, __num_bins_cache;
283
284#if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
285      // Try the L1 cache first.
286
287      // Must fit into L1.
288      __num_bins_cache =
289	std::max<_DifferenceType>(1, __n / (__s.L1_cache_size_lb
290					    / sizeof(_ValueType)));
291      __num_bins_cache = __round_up_to_pow2(__num_bins_cache);
292
293      // No more buckets than TLB entries, power of 2
294      // Power of 2 and at least one element per bin, at most the TLB size.
295      __num_bins = std::min<_DifferenceType>(__n, __num_bins_cache);
296
297#if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
298      // 2 TLB entries needed per bin.
299      __num_bins = std::min<_DifferenceType>(__s.TLB_size / 2, __num_bins);
300#endif
301      __num_bins = __round_up_to_pow2(__num_bins);
302
303      if (__num_bins < __num_bins_cache)
304	{
305#endif
306          // Now try the L2 cache
307          // Must fit into L2
308          __num_bins_cache = static_cast<_BinIndex>
309	    (std::max<_DifferenceType>(1, __n / (__s.L2_cache_size
310						 / sizeof(_ValueType))));
311          __num_bins_cache = __round_up_to_pow2(__num_bins_cache);
312
313          // No more buckets than TLB entries, power of 2.
314          __num_bins = static_cast<_BinIndex>
315	    (std::min(__n, static_cast<_DifferenceType>(__num_bins_cache)));
316          // Power of 2 and at least one element per bin, at most the TLB size.
317#if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
318          // 2 TLB entries needed per bin.
319          __num_bins = std::min(static_cast<_DifferenceType>(__s.TLB_size / 2),
320				__num_bins);
321#endif
322            __num_bins = __round_up_to_pow2(__num_bins);
323#if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
324	}
325#endif
326
327      __num_bins = __round_up_to_pow2(
328                        std::max<_BinIndex>(__num_threads, __num_bins));
329
330      if (__num_threads <= 1)
331      {
332        _RandomNumber __derived_rng(
333                            __rng(std::numeric_limits<uint32_t>::max()));
334	__sequential_random_shuffle(__begin, __end, __derived_rng);
335        return;
336      }
337
338      _DRandomShufflingGlobalData<_RAIter> __sd(__begin);
339      _DRSSorterPU<_RAIter, _RandomNumber >* __pus;
340      _DifferenceType* __starts;
341
342#     pragma omp parallel num_threads(__num_threads)
343      {
344	_ThreadIndex __num_threads = omp_get_num_threads();
345#       pragma omp single
346	{
347	  __pus = new _DRSSorterPU<_RAIter, _RandomNumber>[__num_threads];
348
349	  __sd._M_temporaries = new _ValueType*[__num_threads];
350	  __sd._M_dist = new _DifferenceType*[__num_bins + 1];
351	  __sd._M_bin_proc = new _ThreadIndex[__num_bins];
352	  for (_BinIndex __b = 0; __b < __num_bins + 1; ++__b)
353	    __sd._M_dist[__b] = new _DifferenceType[__num_threads + 1];
354	  for (_BinIndex __b = 0; __b < (__num_bins + 1); ++__b)
355	    {
356	      __sd._M_dist[0][0] = 0;
357	      __sd._M_dist[__b][0] = 0;
358	    }
359	  __starts = __sd._M_starts = new _DifferenceType[__num_threads + 1];
360	  int __bin_cursor = 0;
361	  __sd._M_num_bins = __num_bins;
362	  __sd._M_num_bits = __rd_log2(__num_bins);
363
364	  _DifferenceType __chunk_length = __n / __num_threads,
365	                         __split = __n % __num_threads,
366	                         __start = 0;
367	  _DifferenceType __bin_chunk_length = __num_bins / __num_threads,
368	                         __bin_split = __num_bins % __num_threads;
369	  for (_ThreadIndex __i = 0; __i < __num_threads; ++__i)
370	    {
371	      __starts[__i] = __start;
372	      __start += (__i < __split
373			  ? (__chunk_length + 1) : __chunk_length);
374	      int __j = __pus[__i]._M_bins_begin = __bin_cursor;
375
376	      // Range of bins for this processor.
377	      __bin_cursor += (__i < __bin_split
378			       ? (__bin_chunk_length + 1)
379			       : __bin_chunk_length);
380	      __pus[__i].__bins_end = __bin_cursor;
381	      for (; __j < __bin_cursor; ++__j)
382		__sd._M_bin_proc[__j] = __i;
383	      __pus[__i]._M_num_threads = __num_threads;
384	      __pus[__i]._M_seed = __rng(std::numeric_limits<uint32_t>::max());
385	      __pus[__i]._M_sd = &__sd;
386	    }
387	  __starts[__num_threads] = __start;
388	} //single
389          // Now shuffle in parallel.
390	__parallel_random_shuffle_drs_pu(__pus);
391      }  // parallel
392
393      delete[] __starts;
394      delete[] __sd._M_bin_proc;
395      for (int __s = 0; __s < (__num_bins + 1); ++__s)
396	delete[] __sd._M_dist[__s];
397      delete[] __sd._M_dist;
398      delete[] __sd._M_temporaries;
399
400      delete[] __pus;
401    }
402
403  /** @brief Sequential cache-efficient random shuffle.
404   *  @param __begin Begin iterator of sequence.
405   *  @param __end End iterator of sequence.
406   *  @param __rng Random number generator to use.
407   */
408  template<typename _RAIter, typename _RandomNumberGenerator>
409    void
410    __sequential_random_shuffle(_RAIter __begin, _RAIter __end,
411				_RandomNumberGenerator& __rng)
412    {
413      typedef std::iterator_traits<_RAIter> _TraitsType;
414      typedef typename _TraitsType::value_type _ValueType;
415      typedef typename _TraitsType::difference_type _DifferenceType;
416
417      _DifferenceType __n = __end - __begin;
418      const _Settings& __s = _Settings::get();
419
420      _BinIndex __num_bins, __num_bins_cache;
421
422#if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
423      // Try the L1 cache first, must fit into L1.
424      __num_bins_cache = std::max<_DifferenceType>
425	(1, __n / (__s.L1_cache_size_lb / sizeof(_ValueType)));
426      __num_bins_cache = __round_up_to_pow2(__num_bins_cache);
427
428      // No more buckets than TLB entries, power of 2
429      // Power of 2 and at least one element per bin, at most the TLB size
430      __num_bins = std::min(__n, (_DifferenceType)__num_bins_cache);
431#if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
432      // 2 TLB entries needed per bin
433      __num_bins = std::min((_DifferenceType)__s.TLB_size / 2, __num_bins);
434#endif
435      __num_bins = __round_up_to_pow2(__num_bins);
436
437      if (__num_bins < __num_bins_cache)
438	{
439#endif
440          // Now try the L2 cache, must fit into L2.
441          __num_bins_cache = static_cast<_BinIndex>
442	    (std::max<_DifferenceType>(1, __n / (__s.L2_cache_size
443						 / sizeof(_ValueType))));
444          __num_bins_cache = __round_up_to_pow2(__num_bins_cache);
445
446          // No more buckets than TLB entries, power of 2
447          // Power of 2 and at least one element per bin, at most the TLB size.
448          __num_bins = static_cast<_BinIndex>
449	    (std::min(__n, static_cast<_DifferenceType>(__num_bins_cache)));
450
451#if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
452          // 2 TLB entries needed per bin
453          __num_bins = std::min<_DifferenceType>(__s.TLB_size / 2, __num_bins);
454#endif
455          __num_bins = __round_up_to_pow2(__num_bins);
456#if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
457	}
458#endif
459
460      int __num_bits = __rd_log2(__num_bins);
461
462      if (__num_bins > 1)
463	{
464          _ValueType* __target =
465	    static_cast<_ValueType*>(::operator new(sizeof(_ValueType) * __n));
466          _BinIndex* __oracles = new _BinIndex[__n];
467          _DifferenceType* __dist0 = new _DifferenceType[__num_bins + 1],
468                	 * __dist1 = new _DifferenceType[__num_bins + 1];
469
470          for (int __b = 0; __b < __num_bins + 1; ++__b)
471            __dist0[__b] = 0;
472
473          _RandomNumber __bitrng(__rng(0xFFFFFFFF));
474
475          for (_DifferenceType __i = 0; __i < __n; ++__i)
476            {
477              _BinIndex __oracle = __random_number_pow2(__num_bits, __bitrng);
478              __oracles[__i] = __oracle;
479
480              // To allow prefix (partial) sum.
481              ++(__dist0[__oracle + 1]);
482            }
483
484          // Sum up bins.
485          __gnu_sequential::partial_sum(__dist0, __dist0 + __num_bins + 1,
486					__dist0);
487
488          for (int __b = 0; __b < __num_bins + 1; ++__b)
489            __dist1[__b] = __dist0[__b];
490
491          // Distribute according to oracles.
492          for (_DifferenceType __i = 0; __i < __n; ++__i)
493            ::new(&(__target[(__dist0[__oracles[__i]])++]))
494		_ValueType(*(__begin + __i));
495
496          for (int __b = 0; __b < __num_bins; ++__b)
497	    __sequential_random_shuffle(__target + __dist1[__b],
498					__target + __dist1[__b + 1], __rng);
499
500          // Copy elements back.
501          std::copy(__target, __target + __n, __begin);
502
503          delete[] __dist0;
504          delete[] __dist1;
505          delete[] __oracles;
506
507	  for (_DifferenceType __i = 0; __i < __n; ++__i)
508	    __target[__i].~_ValueType();
509          ::operator delete(__target);
510	}
511      else
512	__gnu_sequential::random_shuffle(__begin, __end, __rng);
513    }
514
515  /** @brief Parallel random public call.
516   *  @param __begin Begin iterator of sequence.
517   *  @param __end End iterator of sequence.
518   *  @param __rng Random number generator to use.
519   */
520  template<typename _RAIter, typename _RandomNumberGenerator>
521    inline void
522    __parallel_random_shuffle(_RAIter __begin, _RAIter __end,
523			      _RandomNumberGenerator __rng = _RandomNumber())
524    {
525      typedef std::iterator_traits<_RAIter> _TraitsType;
526      typedef typename _TraitsType::difference_type _DifferenceType;
527      _DifferenceType __n = __end - __begin;
528      __parallel_random_shuffle_drs(__begin, __end, __n,
529				    __get_max_threads(), __rng);
530    }
531}
532
533#endif /* _GLIBCXX_PARALLEL_RANDOM_SHUFFLE_H */
534