libstdc++
random_shuffle.h
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
00004 //
00005 // This file is part of the GNU ISO C++ Library.  This library is free
00006 // software; you can redistribute it and/or modify it under the terms
00007 // of the GNU General Public License as published by the Free Software
00008 // Foundation; either version 3, or (at your option) any later
00009 // version.
00010 
00011 // This library is distributed in the hope that it will be useful, but
00012 // WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014 // General Public License for more details.
00015 
00016 // Under Section 7 of GPL version 3, you are granted additional
00017 // permissions described in the GCC Runtime Library Exception, version
00018 // 3.1, as published by the Free Software Foundation.
00019 
00020 // You should have received a copy of the GNU General Public License and
00021 // a copy of the GCC Runtime Library Exception along with this program;
00022 // see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
00023 // <http://www.gnu.org/licenses/>.
00024 
00025 /** @file parallel/random_shuffle.h
00026  *  @brief Parallel implementation of std::random_shuffle().
00027  *  This file is a GNU parallel extension to the Standard C++ Library.
00028  */
00029 
00030 // Written by Johannes Singler.
00031 
00032 #ifndef _GLIBCXX_PARALLEL_RANDOM_SHUFFLE_H
00033 #define _GLIBCXX_PARALLEL_RANDOM_SHUFFLE_H 1
00034 
00035 #include <limits>
00036 #include <bits/stl_numeric.h>
00037 #include <parallel/parallel.h>
00038 #include <parallel/random_number.h>
00039 
00040 namespace __gnu_parallel
00041 {
00042   /** @brief Type to hold the index of a bin.
00043     *
00044     *  Since many variables of this type are allocated, it should be
00045     *  chosen as small as possible.
00046     */
00047   typedef unsigned short _BinIndex;
00048 
00049   /** @brief Data known to every thread participating in
00050       __gnu_parallel::__parallel_random_shuffle(). */
00051   template<typename _RAIter>
00052     struct _DRandomShufflingGlobalData
00053     {
00054       typedef std::iterator_traits<_RAIter> _TraitsType;
00055       typedef typename _TraitsType::value_type _ValueType;
00056       typedef typename _TraitsType::difference_type _DifferenceType;
00057 
00058       /** @brief Begin iterator of the __source. */
00059       _RAIter& _M_source;
00060 
00061       /** @brief Temporary arrays for each thread. */
00062       _ValueType** _M_temporaries;
00063 
00064       /** @brief Two-dimensional array to hold the thread-bin distribution.
00065        *
00066        *  Dimensions (_M_num_threads + 1) __x (_M_num_bins + 1). */
00067       _DifferenceType** _M_dist;
00068 
00069       /** @brief Start indexes of the threads' __chunks. */
00070       _DifferenceType* _M_starts;
00071 
00072       /** @brief Number of the thread that will further process the
00073           corresponding bin. */
00074       _ThreadIndex* _M_bin_proc;
00075 
00076       /** @brief Number of bins to distribute to. */
00077       int _M_num_bins;
00078 
00079       /** @brief Number of bits needed to address the bins. */
00080       int _M_num_bits;
00081 
00082       /** @brief Constructor. */
00083       _DRandomShufflingGlobalData(_RAIter& __source)
00084       : _M_source(__source) { }
00085     };
00086 
00087   /** @brief Local data for a thread participating in
00088       __gnu_parallel::__parallel_random_shuffle().
00089     */
00090   template<typename _RAIter, typename _RandomNumberGenerator>
00091     struct _DRSSorterPU
00092     {
00093       /** @brief Number of threads participating in total. */
00094       int _M_num_threads;
00095 
00096       /** @brief Begin index for bins taken care of by this thread. */
00097       _BinIndex _M_bins_begin;
00098 
00099       /** @brief End index for bins taken care of by this thread. */
00100       _BinIndex __bins_end;
00101 
00102       /** @brief Random _M_seed for this thread. */
00103       uint32_t _M_seed;
00104 
00105       /** @brief Pointer to global data. */
00106       _DRandomShufflingGlobalData<_RAIter>* _M_sd;
00107     };
00108 
00109   /** @brief Generate a random number in @c [0,2^__logp).
00110     *  @param __logp Logarithm (basis 2) of the upper range __bound.
00111     *  @param __rng Random number generator to use.
00112     */
00113   template<typename _RandomNumberGenerator>
00114     inline int
00115     __random_number_pow2(int __logp, _RandomNumberGenerator& __rng)
00116     { return __rng.__genrand_bits(__logp); }
00117 
00118   /** @brief Random shuffle code executed by each thread.
00119     *  @param __pus Array of thread-local data records. */
00120   template<typename _RAIter, typename _RandomNumberGenerator>
00121     void 
00122     __parallel_random_shuffle_drs_pu(_DRSSorterPU<_RAIter,
00123                      _RandomNumberGenerator>* __pus)
00124     {
00125       typedef std::iterator_traits<_RAIter> _TraitsType;
00126       typedef typename _TraitsType::value_type _ValueType;
00127       typedef typename _TraitsType::difference_type _DifferenceType;
00128 
00129       _ThreadIndex __iam = omp_get_thread_num();
00130       _DRSSorterPU<_RAIter, _RandomNumberGenerator>* __d = &__pus[__iam];
00131       _DRandomShufflingGlobalData<_RAIter>* __sd = __d->_M_sd;
00132 
00133       // Indexing: _M_dist[bin][processor]
00134       _DifferenceType __length = (__sd->_M_starts[__iam + 1]
00135                   - __sd->_M_starts[__iam]);
00136       _BinIndex* __oracles = new _BinIndex[__length];
00137       _DifferenceType* __dist = new _DifferenceType[__sd->_M_num_bins + 1];
00138       _BinIndex* __bin_proc = new _BinIndex[__sd->_M_num_bins];
00139       _ValueType** __temporaries = new _ValueType*[__d->_M_num_threads];
00140 
00141       // Compute oracles and count appearances.
00142       for (_BinIndex __b = 0; __b < __sd->_M_num_bins + 1; ++__b)
00143     __dist[__b] = 0;
00144       int __num_bits = __sd->_M_num_bits;
00145 
00146       _RandomNumber __rng(__d->_M_seed);
00147 
00148       // First main loop.
00149       for (_DifferenceType __i = 0; __i < __length; ++__i)
00150     {
00151           _BinIndex __oracle = __random_number_pow2(__num_bits, __rng);
00152           __oracles[__i] = __oracle;
00153 
00154           // To allow prefix (partial) sum.
00155           ++(__dist[__oracle + 1]);
00156     }
00157 
00158       for (_BinIndex __b = 0; __b < __sd->_M_num_bins + 1; ++__b)
00159     __sd->_M_dist[__b][__iam + 1] = __dist[__b];
00160 
00161 #     pragma omp barrier
00162 
00163 #     pragma omp single
00164       {
00165     // Sum up bins, __sd->_M_dist[__s + 1][__d->_M_num_threads] now
00166     // contains the total number of items in bin __s
00167     for (_BinIndex __s = 0; __s < __sd->_M_num_bins; ++__s)
00168           __gnu_sequential::partial_sum(__sd->_M_dist[__s + 1],
00169                     __sd->_M_dist[__s + 1]
00170                     + __d->_M_num_threads + 1,
00171                     __sd->_M_dist[__s + 1]);
00172       }
00173 
00174 #     pragma omp barrier
00175 
00176       _SequenceIndex __offset = 0, __global_offset = 0;
00177       for (_BinIndex __s = 0; __s < __d->_M_bins_begin; ++__s)
00178     __global_offset += __sd->_M_dist[__s + 1][__d->_M_num_threads];
00179 
00180 #     pragma omp barrier
00181 
00182       for (_BinIndex __s = __d->_M_bins_begin; __s < __d->__bins_end; ++__s)
00183     {
00184           for (int __t = 0; __t < __d->_M_num_threads + 1; ++__t)
00185             __sd->_M_dist[__s + 1][__t] += __offset;
00186           __offset = __sd->_M_dist[__s + 1][__d->_M_num_threads];
00187     }
00188 
00189       __sd->_M_temporaries[__iam] = static_cast<_ValueType*>
00190     (::operator new(sizeof(_ValueType) * __offset));
00191 
00192 #     pragma omp barrier
00193 
00194       // Draw local copies to avoid false sharing.
00195       for (_BinIndex __b = 0; __b < __sd->_M_num_bins + 1; ++__b)
00196     __dist[__b] = __sd->_M_dist[__b][__iam];
00197       for (_BinIndex __b = 0; __b < __sd->_M_num_bins; ++__b)
00198     __bin_proc[__b] = __sd->_M_bin_proc[__b];
00199       for (_ThreadIndex __t = 0; __t < __d->_M_num_threads; ++__t)
00200     __temporaries[__t] = __sd->_M_temporaries[__t];
00201 
00202       _RAIter __source = __sd->_M_source;
00203       _DifferenceType __start = __sd->_M_starts[__iam];
00204 
00205       // Distribute according to oracles, second main loop.
00206       for (_DifferenceType __i = 0; __i < __length; ++__i)
00207     {
00208           _BinIndex __target_bin = __oracles[__i];
00209           _ThreadIndex __target_p = __bin_proc[__target_bin];
00210 
00211           // Last column [__d->_M_num_threads] stays unchanged.
00212       ::new(&(__temporaries[__target_p][__dist[__target_bin + 1]++]))
00213               _ValueType(*(__source + __i + __start));
00214     }
00215 
00216       delete[] __oracles;
00217       delete[] __dist;
00218       delete[] __bin_proc;
00219       delete[] __temporaries;
00220 
00221 #     pragma omp barrier
00222 
00223       // Shuffle bins internally.
00224       for (_BinIndex __b = __d->_M_bins_begin; __b < __d->__bins_end; ++__b)
00225     {
00226           _ValueType* __begin =
00227         (__sd->_M_temporaries[__iam]
00228          + (__b == __d->_M_bins_begin
00229         ? 0 : __sd->_M_dist[__b][__d->_M_num_threads])),
00230         *__end = (__sd->_M_temporaries[__iam]
00231               + __sd->_M_dist[__b + 1][__d->_M_num_threads]);
00232 
00233           __sequential_random_shuffle(__begin, __end, __rng);
00234           std::copy(__begin, __end, __sd->_M_source + __global_offset
00235             + (__b == __d->_M_bins_begin
00236                ? 0 : __sd->_M_dist[__b][__d->_M_num_threads]));
00237     }
00238 
00239       for (_SequenceIndex __i = 0; __i < __offset; ++__i)
00240     __sd->_M_temporaries[__iam][__i].~_ValueType();
00241       ::operator delete(__sd->_M_temporaries[__iam]);
00242     }
00243 
00244   /** @brief Round up to the next greater power of 2.
00245     *  @param __x _Integer to round up */
00246   template<typename _Tp>
00247     _Tp 
00248     __round_up_to_pow2(_Tp __x)
00249     {
00250       if (__x <= 1)
00251     return 1;
00252       else
00253     return (_Tp)1 << (__rd_log2(__x - 1) + 1);
00254     }
00255 
00256   /** @brief Main parallel random shuffle step.
00257     *  @param __begin Begin iterator of sequence.
00258     *  @param __end End iterator of sequence.
00259     *  @param __n Length of sequence.
00260     *  @param __num_threads Number of threads to use.
00261     *  @param __rng Random number generator to use.
00262     */
00263   template<typename _RAIter, typename _RandomNumberGenerator>
00264     void
00265     __parallel_random_shuffle_drs(_RAIter __begin, _RAIter __end,
00266                   typename std::iterator_traits
00267                   <_RAIter>::difference_type __n,
00268                   _ThreadIndex __num_threads,
00269                   _RandomNumberGenerator& __rng)
00270     {
00271       typedef std::iterator_traits<_RAIter> _TraitsType;
00272       typedef typename _TraitsType::value_type _ValueType;
00273       typedef typename _TraitsType::difference_type _DifferenceType;
00274 
00275       _GLIBCXX_CALL(__n)
00276 
00277       const _Settings& __s = _Settings::get();
00278 
00279       if (__num_threads > __n)
00280     __num_threads = static_cast<_ThreadIndex>(__n);
00281 
00282       _BinIndex __num_bins, __num_bins_cache;
00283 
00284 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
00285       // Try the L1 cache first.
00286 
00287       // Must fit into L1.
00288       __num_bins_cache =
00289     std::max<_DifferenceType>(1, __n / (__s.L1_cache_size_lb
00290                         / sizeof(_ValueType)));
00291       __num_bins_cache = __round_up_to_pow2(__num_bins_cache);
00292 
00293       // No more buckets than TLB entries, power of 2
00294       // Power of 2 and at least one element per bin, at most the TLB size.
00295       __num_bins = std::min<_DifferenceType>(__n, __num_bins_cache);
00296 
00297 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
00298       // 2 TLB entries needed per bin.
00299       __num_bins = std::min<_DifferenceType>(__s.TLB_size / 2, __num_bins);
00300 #endif
00301       __num_bins = __round_up_to_pow2(__num_bins);
00302 
00303       if (__num_bins < __num_bins_cache)
00304     {
00305 #endif
00306           // Now try the L2 cache
00307           // Must fit into L2
00308           __num_bins_cache = static_cast<_BinIndex>
00309         (std::max<_DifferenceType>(1, __n / (__s.L2_cache_size
00310                          / sizeof(_ValueType))));
00311           __num_bins_cache = __round_up_to_pow2(__num_bins_cache);
00312 
00313           // No more buckets than TLB entries, power of 2.
00314           __num_bins = static_cast<_BinIndex>
00315         (std::min(__n, static_cast<_DifferenceType>(__num_bins_cache)));
00316           // Power of 2 and at least one element per bin, at most the TLB size.
00317 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
00318           // 2 TLB entries needed per bin.
00319           __num_bins = std::min(static_cast<_DifferenceType>(__s.TLB_size / 2),
00320                 __num_bins);
00321 #endif
00322             __num_bins = __round_up_to_pow2(__num_bins);
00323 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
00324     }
00325 #endif
00326 
00327       __num_bins = __round_up_to_pow2(
00328                         std::max<_BinIndex>(__num_threads, __num_bins));
00329 
00330       if (__num_threads <= 1)
00331       {
00332         _RandomNumber __derived_rng(
00333                             __rng(std::numeric_limits<uint32_t>::max()));
00334     __sequential_random_shuffle(__begin, __end, __derived_rng);
00335         return;
00336       }
00337 
00338       _DRandomShufflingGlobalData<_RAIter> __sd(__begin);
00339       _DRSSorterPU<_RAIter, _RandomNumber >* __pus;
00340       _DifferenceType* __starts;
00341 
00342 #     pragma omp parallel num_threads(__num_threads)
00343       {
00344     _ThreadIndex __num_threads = omp_get_num_threads();
00345 #       pragma omp single
00346     {
00347       __pus = new _DRSSorterPU<_RAIter, _RandomNumber>[__num_threads];
00348       
00349       __sd._M_temporaries = new _ValueType*[__num_threads];
00350       __sd._M_dist = new _DifferenceType*[__num_bins + 1];
00351       __sd._M_bin_proc = new _ThreadIndex[__num_bins];
00352       for (_BinIndex __b = 0; __b < __num_bins + 1; ++__b)
00353         __sd._M_dist[__b] = new _DifferenceType[__num_threads + 1];
00354       for (_BinIndex __b = 0; __b < (__num_bins + 1); ++__b)
00355         {
00356           __sd._M_dist[0][0] = 0;
00357           __sd._M_dist[__b][0] = 0;
00358         }
00359       __starts = __sd._M_starts = new _DifferenceType[__num_threads + 1];
00360       int __bin_cursor = 0;
00361       __sd._M_num_bins = __num_bins;
00362       __sd._M_num_bits = __rd_log2(__num_bins);
00363 
00364       _DifferenceType __chunk_length = __n / __num_threads,
00365                              __split = __n % __num_threads,
00366                              __start = 0;
00367       _DifferenceType __bin_chunk_length = __num_bins / __num_threads,
00368                              __bin_split = __num_bins % __num_threads;
00369       for (_ThreadIndex __i = 0; __i < __num_threads; ++__i)
00370         {
00371           __starts[__i] = __start;
00372           __start += (__i < __split
00373               ? (__chunk_length + 1) : __chunk_length);
00374           int __j = __pus[__i]._M_bins_begin = __bin_cursor;
00375 
00376           // Range of bins for this processor.
00377           __bin_cursor += (__i < __bin_split
00378                    ? (__bin_chunk_length + 1)
00379                    : __bin_chunk_length);
00380           __pus[__i].__bins_end = __bin_cursor;
00381           for (; __j < __bin_cursor; ++__j)
00382         __sd._M_bin_proc[__j] = __i;
00383           __pus[__i]._M_num_threads = __num_threads;
00384           __pus[__i]._M_seed = __rng(std::numeric_limits<uint32_t>::max());
00385           __pus[__i]._M_sd = &__sd;
00386         }
00387       __starts[__num_threads] = __start;
00388     } //single
00389           // Now shuffle in parallel.
00390     __parallel_random_shuffle_drs_pu(__pus);
00391       }  // parallel
00392 
00393       delete[] __starts;
00394       delete[] __sd._M_bin_proc;
00395       for (int __s = 0; __s < (__num_bins + 1); ++__s)
00396     delete[] __sd._M_dist[__s];
00397       delete[] __sd._M_dist;
00398       delete[] __sd._M_temporaries;
00399 
00400       delete[] __pus;
00401     }
00402 
00403   /** @brief Sequential cache-efficient random shuffle.
00404    *  @param __begin Begin iterator of sequence.
00405    *  @param __end End iterator of sequence.
00406    *  @param __rng Random number generator to use.
00407    */
00408   template<typename _RAIter, typename _RandomNumberGenerator>
00409     void
00410     __sequential_random_shuffle(_RAIter __begin, _RAIter __end,
00411                 _RandomNumberGenerator& __rng)
00412     {
00413       typedef std::iterator_traits<_RAIter> _TraitsType;
00414       typedef typename _TraitsType::value_type _ValueType;
00415       typedef typename _TraitsType::difference_type _DifferenceType;
00416 
00417       _DifferenceType __n = __end - __begin;
00418       const _Settings& __s = _Settings::get();
00419 
00420       _BinIndex __num_bins, __num_bins_cache;
00421 
00422 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
00423       // Try the L1 cache first, must fit into L1.
00424       __num_bins_cache = std::max<_DifferenceType>
00425     (1, __n / (__s.L1_cache_size_lb / sizeof(_ValueType)));
00426       __num_bins_cache = __round_up_to_pow2(__num_bins_cache);
00427 
00428       // No more buckets than TLB entries, power of 2
00429       // Power of 2 and at least one element per bin, at most the TLB size
00430       __num_bins = std::min(__n, (_DifferenceType)__num_bins_cache);
00431 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
00432       // 2 TLB entries needed per bin
00433       __num_bins = std::min((_DifferenceType)__s.TLB_size / 2, __num_bins);
00434 #endif
00435       __num_bins = __round_up_to_pow2(__num_bins);
00436 
00437       if (__num_bins < __num_bins_cache)
00438     {
00439 #endif
00440           // Now try the L2 cache, must fit into L2.
00441           __num_bins_cache = static_cast<_BinIndex>
00442         (std::max<_DifferenceType>(1, __n / (__s.L2_cache_size
00443                          / sizeof(_ValueType))));
00444           __num_bins_cache = __round_up_to_pow2(__num_bins_cache);
00445 
00446           // No more buckets than TLB entries, power of 2
00447           // Power of 2 and at least one element per bin, at most the TLB size.
00448           __num_bins = static_cast<_BinIndex>
00449         (std::min(__n, static_cast<_DifferenceType>(__num_bins_cache)));
00450 
00451 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
00452           // 2 TLB entries needed per bin
00453           __num_bins = std::min<_DifferenceType>(__s.TLB_size / 2, __num_bins);
00454 #endif
00455           __num_bins = __round_up_to_pow2(__num_bins);
00456 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
00457     }
00458 #endif
00459 
00460       int __num_bits = __rd_log2(__num_bins);
00461 
00462       if (__num_bins > 1)
00463     {
00464           _ValueType* __target =
00465         static_cast<_ValueType*>(::operator new(sizeof(_ValueType) * __n));
00466           _BinIndex* __oracles = new _BinIndex[__n];
00467           _DifferenceType* __dist0 = new _DifferenceType[__num_bins + 1],
00468                      * __dist1 = new _DifferenceType[__num_bins + 1];
00469 
00470           for (int __b = 0; __b < __num_bins + 1; ++__b)
00471             __dist0[__b] = 0;
00472 
00473           _RandomNumber __bitrng(__rng(0xFFFFFFFF));
00474 
00475           for (_DifferenceType __i = 0; __i < __n; ++__i)
00476             {
00477               _BinIndex __oracle = __random_number_pow2(__num_bits, __bitrng);
00478               __oracles[__i] = __oracle;
00479 
00480               // To allow prefix (partial) sum.
00481               ++(__dist0[__oracle + 1]);
00482             }
00483 
00484           // Sum up bins.
00485           __gnu_sequential::partial_sum(__dist0, __dist0 + __num_bins + 1,
00486                     __dist0);
00487 
00488           for (int __b = 0; __b < __num_bins + 1; ++__b)
00489             __dist1[__b] = __dist0[__b];
00490 
00491           // Distribute according to oracles.
00492           for (_DifferenceType __i = 0; __i < __n; ++__i)
00493             ::new(&(__target[(__dist0[__oracles[__i]])++])) 
00494         _ValueType(*(__begin + __i));
00495 
00496           for (int __b = 0; __b < __num_bins; ++__b)
00497         __sequential_random_shuffle(__target + __dist1[__b],
00498                     __target + __dist1[__b + 1], __rng);
00499 
00500           // Copy elements back.
00501           std::copy(__target, __target + __n, __begin);
00502 
00503           delete[] __dist0;
00504           delete[] __dist1;
00505           delete[] __oracles;
00506       
00507       for (_DifferenceType __i = 0; __i < __n; ++__i)
00508         __target[__i].~_ValueType();
00509           ::operator delete(__target);
00510     }
00511       else
00512     __gnu_sequential::random_shuffle(__begin, __end, __rng);
00513     }
00514 
00515   /** @brief Parallel random public call.
00516    *  @param __begin Begin iterator of sequence.
00517    *  @param __end End iterator of sequence.
00518    *  @param __rng Random number generator to use.
00519    */
00520   template<typename _RAIter, typename _RandomNumberGenerator>
00521     inline void
00522     __parallel_random_shuffle(_RAIter __begin, _RAIter __end,
00523                   _RandomNumberGenerator __rng = _RandomNumber())
00524     {
00525       typedef std::iterator_traits<_RAIter> _TraitsType;
00526       typedef typename _TraitsType::difference_type _DifferenceType;
00527       _DifferenceType __n = __end - __begin;
00528       __parallel_random_shuffle_drs(__begin, __end, __n,
00529                     __get_max_threads(), __rng);
00530     }
00531 }
00532 
00533 #endif /* _GLIBCXX_PARALLEL_RANDOM_SHUFFLE_H */