Home | History | Annotate | Download | only in parallel
      1 // -*- C++ -*-
      2 
      3 // Copyright (C) 2007-2013 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 
     40 namespace __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