l4/pkg/libstdc++-v3/contrib/libstdc++-v3-4.3.3/include/parallel/random_shuffle.h

   1 // -*- C++ -*-
   2
   3 // Copyright (C) 2007, 2008 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 2, 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 // You should have received a copy of the GNU General Public License
  17 // along with this library; see the file COPYING.  If not, write to
  18 // the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
  19 // MA 02111-1307, USA.
  20
  21 // As a special exception, you may use this file as part of a free
  22 // software library without restriction.  Specifically, if other files
  23 // instantiate templates or use macros or inline functions from this
  24 // file, or you compile this file and link it with other files to
  25 // produce an executable, this file does not by itself cause the
  26 // resulting executable to be covered by the GNU General Public
  27 // License.  This exception does not however invalidate any other
  28 // reasons why the executable file might be covered by the GNU General
  29 // Public License.
  30
  31 /** @file parallel/random_shuffle.h
  32  *  @brief Parallel implementation of std::random_shuffle().
  33  *  This file is a GNU parallel extension to the Standard C++ Library.
  34  */
  35
  36 // Written by Johannes Singler.
  37
  38 #ifndef _GLIBCXX_PARALLEL_RANDOM_SHUFFLE_H
  39 #define _GLIBCXX_PARALLEL_RANDOM_SHUFFLE_H 1
  40
  41 #include <limits>
  42 #include <bits/stl_numeric.h>
  43 #include <parallel/parallel.h>
  44 #include <parallel/random_number.h>
  45
  46 namespace __gnu_parallel
  47 {
  48 /** @brief Type to hold the index of a bin.
  49   *
  50   *  Since many variables of this type are allocated, it should be
  51   *  chosen as small as possible.
  52   */
  53 typedef unsigned short bin_index;
  54
  55 /** @brief Data known to every thread participating in
  56     __gnu_parallel::parallel_random_shuffle(). */
  57 template<typename RandomAccessIterator>
  58   struct DRandomShufflingGlobalData
  59   {
  60     typedef std::iterator_traits<RandomAccessIterator> traits_type;
  61     typedef typename traits_type::value_type value_type;
  62     typedef typename traits_type::difference_type difference_type;
  63
  64     /** @brief Begin iterator of the source. */
  65     RandomAccessIterator& source;
  66
  67     /** @brief Temporary arrays for each thread. */
  68     value_type** temporaries;
  69
  70     /** @brief Two-dimensional array to hold the thread-bin distribution.
  71      *
  72      *  Dimensions (num_threads + 1) x (num_bins + 1). */
  73     difference_type** dist;
  74
  75     /** @brief Start indexes of the threads' chunks. */
  76     difference_type* starts;
  77
  78     /** @brief Number of the thread that will further process the
  79         corresponding bin. */
  80     thread_index_t* bin_proc;
  81
  82     /** @brief Number of bins to distribute to. */
  83     int num_bins;
  84
  85     /** @brief Number of bits needed to address the bins. */
  86     int num_bits;
  87
  88     /** @brief Constructor. */
  89     DRandomShufflingGlobalData(RandomAccessIterator& _source)
  90     : source(_source) { }
  91   };
  92
  93 /** @brief Local data for a thread participating in
  94     __gnu_parallel::parallel_random_shuffle().
  95   */
  96 template<typename RandomAccessIterator, typename RandomNumberGenerator>
  97   struct DRSSorterPU
  98   {
  99     /** @brief Number of threads participating in total. */
 100     int num_threads;
 101
 102     /** @brief Begin index for bins taken care of by this thread. */
 103     bin_index bins_begin;
 104
 105     /** @brief End index for bins taken care of by this thread. */
 106     bin_index bins_end;
 107
 108     /** @brief Random seed for this thread. */
 109     uint32 seed;
 110
 111     /** @brief Pointer to global data. */
 112     DRandomShufflingGlobalData<RandomAccessIterator>* sd;
 113   };
 114
 115 /** @brief Generate a random number in @c [0,2^logp).
 116   *  @param logp Logarithm (basis 2) of the upper range bound.
 117   *  @param rng Random number generator to use.
 118   */
 119 template<typename RandomNumberGenerator>
 120   inline int
 121   random_number_pow2(int logp, RandomNumberGenerator& rng)
 122   { return rng.genrand_bits(logp); }
 123
 124 /** @brief Random shuffle code executed by each thread.
 125   *  @param pus Array of thread-local data records. */
 126 template<typename RandomAccessIterator, typename RandomNumberGenerator>
 127   void
 128   parallel_random_shuffle_drs_pu(DRSSorterPU<RandomAccessIterator,
 129                                  RandomNumberGenerator>* pus)
 130   {
 131     typedef std::iterator_traits<RandomAccessIterator> traits_type;
 132     typedef typename traits_type::value_type value_type;
 133     typedef typename traits_type::difference_type difference_type;
 134
 135     thread_index_t iam = omp_get_thread_num();
 136     DRSSorterPU<RandomAccessIterator, RandomNumberGenerator>* d = &pus[iam];
 137     DRandomShufflingGlobalData<RandomAccessIterator>* sd = d->sd;
 138
 139     // Indexing: dist[bin][processor]
 140     difference_type length = sd->starts[iam + 1] - sd->starts[iam];
 141     bin_index* oracles = new bin_index[length];
 142     difference_type* dist = new difference_type[sd->num_bins + 1];
 143     bin_index* bin_proc = new bin_index[sd->num_bins];
 144     value_type** temporaries = new value_type*[d->num_threads];
 145
 146     // Compute oracles and count appearances.
 147     for (bin_index b = 0; b < sd->num_bins + 1; ++b)
 148       dist[b] = 0;
 149     int num_bits = sd->num_bits;
 150
 151     random_number rng(d->seed);
 152
 153     // First main loop.
 154     for (difference_type i = 0; i < length; ++i)
 155       {
 156         bin_index oracle = random_number_pow2(num_bits, rng);
 157         oracles[i] = oracle;
 158
 159         // To allow prefix (partial) sum.
 160         ++(dist[oracle + 1]);
 161       }
 162
 163     for (bin_index b = 0; b < sd->num_bins + 1; ++b)
 164       sd->dist[b][iam + 1] = dist[b];
 165
 166 #   pragma omp barrier
 167
 168 #   pragma omp single
 169     {
 170       // Sum up bins, sd->dist[s + 1][d->num_threads] now contains the
 171       // total number of items in bin s
 172       for (bin_index s = 0; s < sd->num_bins; ++s)
 173         __gnu_sequential::partial_sum(sd->dist[s + 1],
 174                                       sd->dist[s + 1] + d->num_threads + 1,
 175                                       sd->dist[s + 1]);
 176     }
 177
 178 #   pragma omp barrier
 179
 180     sequence_index_t offset = 0, global_offset = 0;
 181     for (bin_index s = 0; s < d->bins_begin; ++s)
 182       global_offset += sd->dist[s + 1][d->num_threads];
 183
 184 #   pragma omp barrier
 185
 186     for (bin_index s = d->bins_begin; s < d->bins_end; ++s)
 187       {
 188         for (int t = 0; t < d->num_threads + 1; ++t)
 189           sd->dist[s + 1][t] += offset;
 190         offset = sd->dist[s + 1][d->num_threads];
 191       }
 192
 193     sd->temporaries[iam] = static_cast<value_type*>(
 194       ::operator new(sizeof(value_type) * offset));
 195
 196 #   pragma omp barrier
 197
 198     // Draw local copies to avoid false sharing.
 199     for (bin_index b = 0; b < sd->num_bins + 1; ++b)
 200       dist[b] = sd->dist[b][iam];
 201     for (bin_index b = 0; b < sd->num_bins; ++b)
 202       bin_proc[b] = sd->bin_proc[b];
 203     for (thread_index_t t = 0; t < d->num_threads; ++t)
 204       temporaries[t] = sd->temporaries[t];
 205
 206     RandomAccessIterator source = sd->source;
 207     difference_type start = sd->starts[iam];
 208
 209     // Distribute according to oracles, second main loop.
 210     for (difference_type i = 0; i < length; ++i)
 211       {
 212         bin_index target_bin = oracles[i];
 213         thread_index_t target_p = bin_proc[target_bin];
 214
 215         // Last column [d->num_threads] stays unchanged.
 216         ::new(&(temporaries[target_p][dist[target_bin + 1]++]))
 217             value_type(*(source + i + start));
 218       }
 219
 220     delete[] oracles;
 221     delete[] dist;
 222     delete[] bin_proc;
 223     delete[] temporaries;
 224
 225 #   pragma omp barrier
 226
 227     // Shuffle bins internally.
 228     for (bin_index b = d->bins_begin; b < d->bins_end; ++b)
 229       {
 230         value_type* begin =
 231                     sd->temporaries[iam] +
 232                     ((b == d->bins_begin) ? 0 : sd->dist[b][d->num_threads]),
 233                   * end =
 234                     sd->temporaries[iam] + sd->dist[b + 1][d->num_threads];
 235         sequential_random_shuffle(begin, end, rng);
 236         std::copy(begin, end, sd->source + global_offset +
 237             ((b == d->bins_begin) ? 0 : sd->dist[b][d->num_threads]));
 238       }
 239
 240     ::operator delete(sd->temporaries[iam]);
 241   }
 242
 243 /** @brief Round up to the next greater power of 2.
 244   *  @param x Integer to round up */
 245 template<typename T>
 246   T
 247   round_up_to_pow2(T x)
 248   {
 249     if (x <= 1)
 250       return 1;
 251     else
 252       return (T)1 << (log2(x - 1) + 1);
 253   }
 254
 255 /** @brief Main parallel random shuffle step.
 256   *  @param begin Begin iterator of sequence.
 257   *  @param end End iterator of sequence.
 258   *  @param n Length of sequence.
 259   *  @param num_threads Number of threads to use.
 260   *  @param rng Random number generator to use.
 261   */
 262 template<typename RandomAccessIterator, typename RandomNumberGenerator>
 263   void
 264   parallel_random_shuffle_drs(RandomAccessIterator begin,
 265                               RandomAccessIterator end,
 266                               typename std::iterator_traits
 267                               <RandomAccessIterator>::difference_type n,
 268                               thread_index_t num_threads,
 269                               RandomNumberGenerator& rng)
 270   {
 271     typedef std::iterator_traits<RandomAccessIterator> traits_type;
 272     typedef typename traits_type::value_type value_type;
 273     typedef typename traits_type::difference_type difference_type;
 274
 275     _GLIBCXX_CALL(n)
 276
 277     const _Settings& __s = _Settings::get();
 278
 279     if (num_threads > n)
 280       num_threads = static_cast<thread_index_t>(n);
 281
 282     bin_index 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 = std::max<difference_type>(
 289         1, n / (__s.L1_cache_size_lb / sizeof(value_type)));
 290     num_bins_cache = round_up_to_pow2(num_bins_cache);
 291
 292     // No more buckets than TLB entries, power of 2
 293     // Power of 2 and at least one element per bin, at most the TLB size.
 294     num_bins = std::min<difference_type>(n, num_bins_cache);
 295
 296 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
 297     // 2 TLB entries needed per bin.
 298     num_bins = std::min<difference_type>(__s.TLB_size / 2, num_bins);
 299 #endif
 300     num_bins = round_up_to_pow2(num_bins);
 301
 302     if (num_bins < num_bins_cache)
 303       {
 304 #endif
 305         // Now try the L2 cache
 306         // Must fit into L2
 307         num_bins_cache = static_cast<bin_index>(std::max<difference_type>(
 308             1, n / (__s.L2_cache_size / sizeof(value_type))));
 309         num_bins_cache = round_up_to_pow2(num_bins_cache);
 310
 311         // No more buckets than TLB entries, power of 2.
 312         num_bins = static_cast<bin_index>(
 313             std::min(n, static_cast<difference_type>(num_bins_cache)));
 314         // Power of 2 and at least one element per bin, at most the TLB size.
 315 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
 316         // 2 TLB entries needed per bin.
 317         num_bins = std::min(
 318             static_cast<difference_type>(__s.TLB_size / 2), num_bins);
 319 #endif
 320           num_bins = round_up_to_pow2(num_bins);
 321 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
 322       }
 323 #endif
 324
 325     num_threads = std::min<bin_index>(num_threads, num_bins);
 326
 327     if (num_threads <= 1)
 328       return sequential_random_shuffle(begin, end, rng);
 329
 330     DRandomShufflingGlobalData<RandomAccessIterator> sd(begin);
 331     DRSSorterPU<RandomAccessIterator, random_number >* pus;
 332     difference_type* starts;
 333
 334 #   pragma omp parallel num_threads(num_threads)
 335       {
 336         thread_index_t num_threads = omp_get_num_threads();
 337 #       pragma omp single
 338           {
 339             pus = new DRSSorterPU<RandomAccessIterator, random_number>
 340                 [num_threads];
 341
 342             sd.temporaries = new value_type*[num_threads];
 343             sd.dist = new difference_type*[num_bins + 1];
 344             sd.bin_proc = new thread_index_t[num_bins];
 345             for (bin_index b = 0; b < num_bins + 1; ++b)
 346               sd.dist[b] = new difference_type[num_threads + 1];
 347             for (bin_index b = 0; b < (num_bins + 1); ++b)
 348               {
 349                 sd.dist[0][0] = 0;
 350                 sd.dist[b][0] = 0;
 351               }
 352             starts = sd.starts = new difference_type[num_threads + 1];
 353             int bin_cursor = 0;
 354             sd.num_bins = num_bins;
 355             sd.num_bits = log2(num_bins);
 356
 357             difference_type chunk_length = n / num_threads,
 358                             split = n % num_threads, start = 0;
 359             difference_type bin_chunk_length = num_bins / num_threads,
 360                             bin_split = num_bins % num_threads;
 361             for (thread_index_t i = 0; i < num_threads; ++i)
 362               {
 363                 starts[i] = start;
 364                 start += (i < split) ? (chunk_length + 1) : chunk_length;
 365                 int j = pus[i].bins_begin = bin_cursor;
 366
 367                 // Range of bins for this processor.
 368                 bin_cursor += (i < bin_split) ?
 369                     (bin_chunk_length + 1) : bin_chunk_length;
 370                 pus[i].bins_end = bin_cursor;
 371                 for (; j < bin_cursor; ++j)
 372                   sd.bin_proc[j] = i;
 373                 pus[i].num_threads = num_threads;
 374                 pus[i].seed = rng(std::numeric_limits<uint32>::max());
 375                 pus[i].sd = &sd;
 376               }
 377             starts[num_threads] = start;
 378           } //single
 379         // Now shuffle in parallel.
 380         parallel_random_shuffle_drs_pu(pus);
 381       }  // parallel
 382
 383     delete[] starts;
 384     delete[] sd.bin_proc;
 385     for (int s = 0; s < (num_bins + 1); ++s)
 386       delete[] sd.dist[s];
 387     delete[] sd.dist;
 388     delete[] sd.temporaries;
 389
 390     delete[] pus;
 391   }
 392
 393 /** @brief Sequential cache-efficient random shuffle.
 394  *  @param begin Begin iterator of sequence.
 395  *  @param end End iterator of sequence.
 396  *  @param rng Random number generator to use.
 397  */
 398 template<typename RandomAccessIterator, typename RandomNumberGenerator>
 399   void
 400   sequential_random_shuffle(RandomAccessIterator begin,
 401                             RandomAccessIterator end,
 402                             RandomNumberGenerator& rng)
 403   {
 404     typedef std::iterator_traits<RandomAccessIterator> traits_type;
 405     typedef typename traits_type::value_type value_type;
 406     typedef typename traits_type::difference_type difference_type;
 407
 408     difference_type n = end - begin;
 409     const _Settings& __s = _Settings::get();
 410
 411     bin_index num_bins, num_bins_cache;
 412
 413 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
 414     // Try the L1 cache first, must fit into L1.
 415     num_bins_cache =
 416         std::max<difference_type>
 417             (1, n / (__s.L1_cache_size_lb / sizeof(value_type)));
 418     num_bins_cache = round_up_to_pow2(num_bins_cache);
 419
 420     // No more buckets than TLB entries, power of 2
 421     // Power of 2 and at least one element per bin, at most the TLB size
 422     num_bins = std::min(n, (difference_type)num_bins_cache);
 423 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
 424     // 2 TLB entries needed per bin
 425     num_bins = std::min((difference_type)__s.TLB_size / 2, num_bins);
 426 #endif
 427     num_bins = round_up_to_pow2(num_bins);
 428
 429     if (num_bins < num_bins_cache)
 430       {
 431 #endif
 432         // Now try the L2 cache, must fit into L2.
 433         num_bins_cache =
 434             static_cast<bin_index>(std::max<difference_type>(
 435                 1, n / (__s.L2_cache_size / sizeof(value_type))));
 436         num_bins_cache = round_up_to_pow2(num_bins_cache);
 437
 438         // No more buckets than TLB entries, power of 2
 439         // Power of 2 and at least one element per bin, at most the TLB size.
 440         num_bins = static_cast<bin_index>
 441             (std::min(n, static_cast<difference_type>(num_bins_cache)));
 442
 443 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
 444         // 2 TLB entries needed per bin
 445         num_bins =
 446             std::min<difference_type>(__s.TLB_size / 2, num_bins);
 447 #endif
 448         num_bins = round_up_to_pow2(num_bins);
 449 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
 450       }
 451 #endif
 452
 453     int num_bits = log2(num_bins);
 454
 455     if (num_bins > 1)
 456       {
 457         value_type* target = static_cast<value_type*>(
 458           ::operator new(sizeof(value_type) * n));
 459         bin_index* oracles = new bin_index[n];
 460         difference_type* dist0 = new difference_type[num_bins + 1],
 461                        * dist1 = new difference_type[num_bins + 1];
 462
 463         for (int b = 0; b < num_bins + 1; ++b)
 464           dist0[b] = 0;
 465
 466         random_number bitrng(rng(0xFFFFFFFF));
 467
 468         for (difference_type i = 0; i < n; ++i)
 469           {
 470             bin_index oracle = random_number_pow2(num_bits, bitrng);
 471             oracles[i] = oracle;
 472
 473             // To allow prefix (partial) sum.
 474             ++(dist0[oracle + 1]);
 475           }
 476
 477         // Sum up bins.
 478         __gnu_sequential::partial_sum(dist0, dist0 + num_bins + 1, dist0);
 479
 480         for (int b = 0; b < num_bins + 1; ++b)
 481           dist1[b] = dist0[b];
 482
 483         // Distribute according to oracles.
 484         for (difference_type i = 0; i < n; ++i)
 485           ::new(&(target[(dist0[oracles[i]])++])) value_type(*(begin + i));
 486
 487         for (int b = 0; b < num_bins; ++b)
 488           {
 489             sequential_random_shuffle(target + dist1[b],
 490                                       target + dist1[b + 1],
 491                                       rng);
 492           }
 493
 494         // Copy elements back.
 495         std::copy(target, target + n, begin);
 496
 497         delete[] dist0;
 498         delete[] dist1;
 499         delete[] oracles;
 500         ::operator delete(target);
 501       }
 502     else
 503       __gnu_sequential::random_shuffle(begin, end, rng);
 504   }
 505
 506 /** @brief Parallel random public call.
 507  *  @param begin Begin iterator of sequence.
 508  *  @param end End iterator of sequence.
 509  *  @param rng Random number generator to use.
 510  */
 511 template<typename RandomAccessIterator, typename RandomNumberGenerator>
 512   inline void
 513   parallel_random_shuffle(RandomAccessIterator begin,
 514                           RandomAccessIterator end,
 515                           RandomNumberGenerator rng = random_number())
 516   {
 517     typedef std::iterator_traits<RandomAccessIterator> traits_type;
 518     typedef typename traits_type::difference_type difference_type;
 519     difference_type n = end - begin;
 520     parallel_random_shuffle_drs(begin, end, n, get_max_threads(), rng) ;
 521   }
 522
 523 }
 524
 525 #endif