libstdc++-v3/include/parallel/random_shuffle.h

   1 // -*- C++ -*-
   2
   3 // Copyright (C) 2007, 2008, 2009 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 _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>* _M_sd = d->_M_sd;
 132
 133     // Indexing: _M_dist[bin][processor]
 134     _DifferenceType __length = _M_sd->_M_starts[__iam + 1] -
 135                                _M_sd->_M_starts[__iam];
 136     _BinIndex* __oracles = new _BinIndex[__length];
 137     _DifferenceType* _M_dist = new _DifferenceType[_M_sd->_M_num_bins + 1];
 138     _BinIndex* _M_bin_proc = new _BinIndex[_M_sd->_M_num_bins];
 139     _ValueType** _M_temporaries = new _ValueType*[d->_M_num_threads];
 140
 141     // Compute oracles and count appearances.
 142     for (_BinIndex __b = 0; __b < _M_sd->_M_num_bins + 1; ++__b)
 143       _M_dist[__b] = 0;
 144     int _M_num_bits = _M_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(_M_num_bits, __rng);
 152         __oracles[__i] = __oracle;
 153
 154         // To allow prefix (partial) sum.
 155         ++(_M_dist[__oracle + 1]);
 156       }
 157
 158     for (_BinIndex __b = 0; __b < _M_sd->_M_num_bins + 1; ++__b)
 159       _M_sd->_M_dist[__b][__iam + 1] = _M_dist[__b];
 160
 161 #   pragma omp barrier
 162
 163 #   pragma omp single
 164     {
 165       // Sum up bins, _M_sd->_M_dist[__s + 1][d->_M_num_threads] now contains
 166       // the total number of items in bin __s
 167       for (_BinIndex __s = 0; __s < _M_sd->_M_num_bins; ++__s)
 168         __gnu_sequential::partial_sum(
 169           _M_sd->_M_dist[__s + 1],
 170           _M_sd->_M_dist[__s + 1] + d->_M_num_threads + 1,
 171           _M_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 += _M_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           _M_sd->_M_dist[__s + 1][__t] += __offset;
 186         __offset = _M_sd->_M_dist[__s + 1][d->_M_num_threads];
 187       }
 188
 189     _M_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 < _M_sd->_M_num_bins + 1; ++__b)
 196       _M_dist[__b] = _M_sd->_M_dist[__b][__iam];
 197     for (_BinIndex __b = 0; __b < _M_sd->_M_num_bins; ++__b)
 198       _M_bin_proc[__b] = _M_sd->_M_bin_proc[__b];
 199     for (_ThreadIndex __t = 0; __t < d->_M_num_threads; ++__t)
 200       _M_temporaries[__t] = _M_sd->_M_temporaries[__t];
 201
 202     _RAIter _M_source = _M_sd->_M_source;
 203     _DifferenceType __start = _M_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 = _M_bin_proc[target_bin];
 210
 211         // Last column [d->_M_num_threads] stays unchanged.
 212         ::new(&(_M_temporaries[target_p][_M_dist[target_bin + 1]++]))
 213             _ValueType(*(_M_source + __i + __start));
 214       }
 215
 216     delete[] __oracles;
 217     delete[] _M_dist;
 218     delete[] _M_bin_proc;
 219     delete[] _M_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                     _M_sd->_M_temporaries[__iam] +
 228                     ((__b == d->_M_bins_begin)
 229                       ? 0 : _M_sd->_M_dist[__b][d->_M_num_threads]),
 230                   * __end =
 231                     _M_sd->_M_temporaries[__iam] +
 232                       _M_sd->_M_dist[__b + 1][d->_M_num_threads];
 233         __sequential_random_shuffle(__begin, __end, __rng);
 234         std::copy(__begin, __end, _M_sd->_M_source + __global_offset +
 235                   ((__b == d->_M_bins_begin)
 236                   ? 0 : _M_sd->_M_dist[__b][d->_M_num_threads]));
 237       }
 238
 239     ::operator delete(_M_sd->_M_temporaries[__iam]);
 240   }
 241
 242 /** @brief Round up to the next greater power of 2.
 243   *  @param __x _Integer to round up */
 244 template<typename _Tp>
 245   _Tp
 246   __round_up_to_pow2(_Tp __x)
 247   {
 248     if (__x <= 1)
 249       return 1;
 250     else
 251       return (_Tp)1 << (__rd_log2(__x - 1) + 1);
 252   }
 253
 254 /** @brief Main parallel random shuffle step.
 255   *  @param __begin Begin iterator of sequence.
 256   *  @param __end End iterator of sequence.
 257   *  @param __n Length of sequence.
 258   *  @param __num_threads Number of threads to use.
 259   *  @param __rng Random number generator to use.
 260   */
 261 template<typename _RAIter, typename _RandomNumberGenerator>
 262   void
 263   __parallel_random_shuffle_drs(_RAIter __begin,
 264                               _RAIter __end,
 265                               typename std::iterator_traits
 266                               <_RAIter>::difference_type __n,
 267                               _ThreadIndex __num_threads,
 268                               _RandomNumberGenerator& __rng)
 269   {
 270     typedef std::iterator_traits<_RAIter> _TraitsType;
 271     typedef typename _TraitsType::value_type _ValueType;
 272     typedef typename _TraitsType::difference_type _DifferenceType;
 273
 274     _GLIBCXX_CALL(__n)
 275
 276     const _Settings& __s = _Settings::get();
 277
 278     if (__num_threads > __n)
 279       __num_threads = static_cast<_ThreadIndex>(__n);
 280
 281     _BinIndex _M_num_bins, __num_bins_cache;
 282
 283 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
 284     // Try the L1 cache first.
 285
 286     // Must fit into L1.
 287     __num_bins_cache = std::max<_DifferenceType>(
 288         1, __n / (__s.L1_cache_size_lb / sizeof(_ValueType)));
 289     __num_bins_cache = __round_up_to_pow2(__num_bins_cache);
 290
 291     // No more buckets than TLB entries, power of 2
 292     // Power of 2 and at least one element per bin, at most the TLB size.
 293     _M_num_bins = std::min<_DifferenceType>(__n, __num_bins_cache);
 294
 295 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
 296     // 2 TLB entries needed per bin.
 297     _M_num_bins = std::min<_DifferenceType>(__s.TLB_size / 2, _M_num_bins);
 298 #endif
 299     _M_num_bins = __round_up_to_pow2(_M_num_bins);
 300
 301     if (_M_num_bins < __num_bins_cache)
 302       {
 303 #endif
 304         // Now try the L2 cache
 305         // Must fit into L2
 306         __num_bins_cache = static_cast<_BinIndex>(std::max<_DifferenceType>(
 307             1, __n / (__s.L2_cache_size / sizeof(_ValueType))));
 308         __num_bins_cache = __round_up_to_pow2(__num_bins_cache);
 309
 310         // No more buckets than TLB entries, power of 2.
 311         _M_num_bins = static_cast<_BinIndex>(
 312             std::min(__n, static_cast<_DifferenceType>(__num_bins_cache)));
 313         // Power of 2 and at least one element per bin, at most the TLB size.
 314 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
 315         // 2 TLB entries needed per bin.
 316         _M_num_bins = std::min(
 317             static_cast<_DifferenceType>(__s.TLB_size / 2), _M_num_bins);
 318 #endif
 319           _M_num_bins = __round_up_to_pow2(_M_num_bins);
 320 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
 321       }
 322 #endif
 323
 324     __num_threads = std::min<_BinIndex>(__num_threads, _M_num_bins);
 325
 326     if (__num_threads <= 1)
 327       return __sequential_random_shuffle(__begin, __end, __rng);
 328
 329     _DRandomShufflingGlobalData<_RAIter> _M_sd(__begin);
 330     _DRSSorterPU<_RAIter, _RandomNumber >* __pus;
 331     _DifferenceType* _M_starts;
 332
 333 #   pragma omp parallel num_threads(__num_threads)
 334       {
 335         _ThreadIndex __num_threads = omp_get_num_threads();
 336 #       pragma omp single
 337           {
 338             __pus = new _DRSSorterPU<_RAIter, _RandomNumber>
 339                 [__num_threads];
 340
 341             _M_sd._M_temporaries = new _ValueType*[__num_threads];
 342             _M_sd._M_dist = new _DifferenceType*[_M_num_bins + 1];
 343             _M_sd._M_bin_proc = new _ThreadIndex[_M_num_bins];
 344             for (_BinIndex __b = 0; __b < _M_num_bins + 1; ++__b)
 345               _M_sd._M_dist[__b] = new _DifferenceType[__num_threads + 1];
 346             for (_BinIndex __b = 0; __b < (_M_num_bins + 1); ++__b)
 347               {
 348                 _M_sd._M_dist[0][0] = 0;
 349                 _M_sd._M_dist[__b][0] = 0;
 350               }
 351             _M_starts = _M_sd._M_starts
 352               = new _DifferenceType[__num_threads + 1];
 353             int bin_cursor = 0;
 354             _M_sd._M_num_bins = _M_num_bins;
 355             _M_sd._M_num_bits = __rd_log2(_M_num_bins);
 356
 357             _DifferenceType __chunk_length = __n / __num_threads,
 358                             __split = __n % __num_threads, __start = 0;
 359             _DifferenceType bin_chunk_length = _M_num_bins / __num_threads,
 360                             bin_split = _M_num_bins % __num_threads;
 361             for (_ThreadIndex __i = 0; __i < __num_threads; ++__i)
 362               {
 363                 _M_starts[__i] = __start;
 364                 __start += (__i < __split)
 365                            ? (__chunk_length + 1) : __chunk_length;
 366                 int __j = __pus[__i]._M_bins_begin = bin_cursor;
 367
 368                 // Range of bins for this processor.
 369                 bin_cursor += (__i < bin_split) ?
 370                     (bin_chunk_length + 1) : bin_chunk_length;
 371                 __pus[__i].__bins_end = bin_cursor;
 372                 for (; __j < bin_cursor; ++__j)
 373                   _M_sd._M_bin_proc[__j] = __i;
 374                 __pus[__i]._M_num_threads = __num_threads;
 375                 __pus[__i]._M_seed = __rng(std::numeric_limits<uint32>::max());
 376                 __pus[__i]._M_sd = &_M_sd;
 377               }
 378             _M_starts[__num_threads] = __start;
 379           } //single
 380         // Now shuffle in parallel.
 381         __parallel_random_shuffle_drs_pu(__pus);
 382       }  // parallel
 383
 384     delete[] _M_starts;
 385     delete[] _M_sd._M_bin_proc;
 386     for (int __s = 0; __s < (_M_num_bins + 1); ++__s)
 387       delete[] _M_sd._M_dist[__s];
 388     delete[] _M_sd._M_dist;
 389     delete[] _M_sd._M_temporaries;
 390
 391     delete[] __pus;
 392   }
 393
 394 /** @brief Sequential cache-efficient random shuffle.
 395  *  @param __begin Begin iterator of sequence.
 396  *  @param __end End iterator of sequence.
 397  *  @param __rng Random number generator to use.
 398  */
 399 template<typename _RAIter, typename _RandomNumberGenerator>
 400   void
 401   __sequential_random_shuffle(_RAIter __begin,
 402                             _RAIter __end,
 403                             _RandomNumberGenerator& __rng)
 404   {
 405     typedef std::iterator_traits<_RAIter> _TraitsType;
 406     typedef typename _TraitsType::value_type _ValueType;
 407     typedef typename _TraitsType::difference_type _DifferenceType;
 408
 409     _DifferenceType __n = __end - __begin;
 410     const _Settings& __s = _Settings::get();
 411
 412     _BinIndex _M_num_bins, __num_bins_cache;
 413
 414 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
 415     // Try the L1 cache first, must fit into L1.
 416     __num_bins_cache =
 417         std::max<_DifferenceType>
 418             (1, __n / (__s.L1_cache_size_lb / sizeof(_ValueType)));
 419     __num_bins_cache = __round_up_to_pow2(__num_bins_cache);
 420
 421     // No more buckets than TLB entries, power of 2
 422     // Power of 2 and at least one element per bin, at most the TLB size
 423     _M_num_bins = std::min(__n, (_DifferenceType)__num_bins_cache);
 424 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
 425     // 2 TLB entries needed per bin
 426     _M_num_bins = std::min((_DifferenceType)__s.TLB_size / 2, _M_num_bins);
 427 #endif
 428     _M_num_bins = __round_up_to_pow2(_M_num_bins);
 429
 430     if (_M_num_bins < __num_bins_cache)
 431       {
 432 #endif
 433         // Now try the L2 cache, must fit into L2.
 434         __num_bins_cache =
 435             static_cast<_BinIndex>(std::max<_DifferenceType>(
 436                 1, __n / (__s.L2_cache_size / sizeof(_ValueType))));
 437         __num_bins_cache = __round_up_to_pow2(__num_bins_cache);
 438
 439         // No more buckets than TLB entries, power of 2
 440         // Power of 2 and at least one element per bin, at most the TLB size.
 441         _M_num_bins = static_cast<_BinIndex>
 442             (std::min(__n, static_cast<_DifferenceType>(__num_bins_cache)));
 443
 444 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
 445         // 2 TLB entries needed per bin
 446         _M_num_bins =
 447             std::min<_DifferenceType>(__s.TLB_size / 2, _M_num_bins);
 448 #endif
 449         _M_num_bins = __round_up_to_pow2(_M_num_bins);
 450 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
 451       }
 452 #endif
 453
 454     int _M_num_bits = __rd_log2(_M_num_bins);
 455
 456     if (_M_num_bins > 1)
 457       {
 458         _ValueType* __target = static_cast<_ValueType*>(
 459           ::operator new(sizeof(_ValueType) * __n));
 460         _BinIndex* __oracles = new _BinIndex[__n];
 461         _DifferenceType* __dist0 = new _DifferenceType[_M_num_bins + 1],
 462                        * __dist1 = new _DifferenceType[_M_num_bins + 1];
 463
 464         for (int __b = 0; __b < _M_num_bins + 1; ++__b)
 465           __dist0[__b] = 0;
 466
 467         _RandomNumber bitrng(__rng(0xFFFFFFFF));
 468
 469         for (_DifferenceType __i = 0; __i < __n; ++__i)
 470           {
 471             _BinIndex __oracle = __random_number_pow2(_M_num_bits, bitrng);
 472             __oracles[__i] = __oracle;
 473
 474             // To allow prefix (partial) sum.
 475             ++(__dist0[__oracle + 1]);
 476           }
 477
 478         // Sum up bins.
 479         __gnu_sequential::
 480             partial_sum(__dist0, __dist0 + _M_num_bins + 1, __dist0);
 481
 482         for (int __b = 0; __b < _M_num_bins + 1; ++__b)
 483           __dist1[__b] = __dist0[__b];
 484
 485         // Distribute according to oracles.
 486         for (_DifferenceType __i = 0; __i < __n; ++__i)
 487           ::new(&(__target[(__dist0[__oracles[__i]])++]))
 488             _ValueType(*(__begin + __i));
 489
 490         for (int __b = 0; __b < _M_num_bins; ++__b)
 491           {
 492             __sequential_random_shuffle(__target + __dist1[__b],
 493                                       __target + __dist1[__b + 1],
 494                                       __rng);
 495           }
 496
 497         // Copy elements back.
 498         std::copy(__target, __target + __n, __begin);
 499
 500         delete[] __dist0;
 501         delete[] __dist1;
 502         delete[] __oracles;
 503         ::operator delete(__target);
 504       }
 505     else
 506       __gnu_sequential::random_shuffle(__begin, __end, __rng);
 507   }
 508
 509 /** @brief Parallel random public call.
 510  *  @param __begin Begin iterator of sequence.
 511  *  @param __end End iterator of sequence.
 512  *  @param __rng Random number generator to use.
 513  */
 514 template<typename _RAIter, typename _RandomNumberGenerator>
 515   inline void
 516   __parallel_random_shuffle(_RAIter __begin,
 517                           _RAIter __end,
 518                           _RandomNumberGenerator __rng = _RandomNumber())
 519   {
 520     typedef std::iterator_traits<_RAIter> _TraitsType;
 521     typedef typename _TraitsType::difference_type _DifferenceType;
 522     _DifferenceType __n = __end - __begin;
 523     __parallel_random_shuffle_drs(
 524       __begin, __end, __n, __get_max_threads(), __rng) ;
 525   }
 526
 527 }
 528
 529 #endif /* _GLIBCXX_PARALLEL_RANDOM_SHUFFLE_H */