libstdc++
multiseq_selection.h
Go to the documentation of this file.
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/multiseq_selection.h
26  * @brief Functions to find elements of a certain global rank in
27  * multiple sorted sequences. Also serves for splitting such
28  * sequence sets.
29  *
30  * The algorithm description can be found in
31  *
32  * P. J. Varman, S. D. Scheufler, B. R. Iyer, and G. R. Ricard.
33  * Merging Multiple Lists on Hierarchical-Memory Multiprocessors.
34  * Journal of Parallel and Distributed Computing, 12(2):171–177, 1991.
35  *
36  * This file is a GNU parallel extension to the Standard C++ Library.
37  */
38 
39 // Written by Johannes Singler.
40 
41 #ifndef _GLIBCXX_PARALLEL_MULTISEQ_SELECTION_H
42 #define _GLIBCXX_PARALLEL_MULTISEQ_SELECTION_H 1
43 
44 #include <vector>
45 #include <queue>
46 
47 #include <bits/stl_algo.h>
48 
49 #include <parallel/sort.h>
50 
51 namespace __gnu_parallel
52 {
53  /** @brief Compare a pair of types lexicographically, ascending. */
54  template<typename T1, typename T2, typename Comparator>
56  : public std::binary_function<std::pair<T1, T2>, std::pair<T1, T2>, bool>
57  {
58  private:
59  Comparator& comp;
60 
61  public:
62  lexicographic(Comparator& _comp) : comp(_comp) { }
63 
64  bool
65  operator()(const std::pair<T1, T2>& p1,
66  const std::pair<T1, T2>& p2) const
67  {
68  if (comp(p1.first, p2.first))
69  return true;
70 
71  if (comp(p2.first, p1.first))
72  return false;
73 
74  // Firsts are equal.
75  return p1.second < p2.second;
76  }
77  };
78 
79  /** @brief Compare a pair of types lexicographically, descending. */
80  template<typename T1, typename T2, typename Comparator>
81  class lexicographic_reverse : public std::binary_function<T1, T2, bool>
82  {
83  private:
84  Comparator& comp;
85 
86  public:
87  lexicographic_reverse(Comparator& _comp) : comp(_comp) { }
88 
89  bool
90  operator()(const std::pair<T1, T2>& p1,
91  const std::pair<T1, T2>& p2) const
92  {
93  if (comp(p2.first, p1.first))
94  return true;
95 
96  if (comp(p1.first, p2.first))
97  return false;
98 
99  // Firsts are equal.
100  return p2.second < p1.second;
101  }
102  };
103 
104  /**
105  * @brief Splits several sorted sequences at a certain global rank,
106  * resulting in a splitting point for each sequence.
107  * The sequences are passed via a sequence of random-access
108  * iterator pairs, none of the sequences may be empty. If there
109  * are several equal elements across the split, the ones on the
110  * left side will be chosen from sequences with smaller number.
111  * @param begin_seqs Begin of the sequence of iterator pairs.
112  * @param end_seqs End of the sequence of iterator pairs.
113  * @param rank The global rank to partition at.
114  * @param begin_offsets A random-access sequence begin where the
115  * result will be stored in. Each element of the sequence is an
116  * iterator that points to the first element on the greater part of
117  * the respective sequence.
118  * @param comp The ordering functor, defaults to std::less<T>.
119  */
120  template<typename RanSeqs, typename RankType, typename RankIterator,
121  typename Comparator>
122  void
123  multiseq_partition(RanSeqs begin_seqs, RanSeqs end_seqs,
124  RankType rank,
125  RankIterator begin_offsets,
126  Comparator comp = std::less<
127  typename std::iterator_traits<typename
128  std::iterator_traits<RanSeqs>::value_type::
129  first_type>::value_type>()) // std::less<T>
130  {
131  _GLIBCXX_CALL(end_seqs - begin_seqs)
132 
134  It;
135  typedef typename std::iterator_traits<It>::difference_type
136  difference_type;
137  typedef typename std::iterator_traits<It>::value_type value_type;
138 
141 
142  // Number of sequences, number of elements in total (possibly
143  // including padding).
144  difference_type m = std::distance(begin_seqs, end_seqs), N = 0,
145  nmax, n, r;
146 
147  for (int i = 0; i < m; i++)
148  {
149  N += std::distance(begin_seqs[i].first, begin_seqs[i].second);
150  _GLIBCXX_PARALLEL_ASSERT(
151  std::distance(begin_seqs[i].first, begin_seqs[i].second) > 0);
152  }
153 
154  if (rank == N)
155  {
156  for (int i = 0; i < m; i++)
157  begin_offsets[i] = begin_seqs[i].second; // Very end.
158  // Return m - 1;
159  return;
160  }
161 
162  _GLIBCXX_PARALLEL_ASSERT(m != 0);
163  _GLIBCXX_PARALLEL_ASSERT(N != 0);
164  _GLIBCXX_PARALLEL_ASSERT(rank >= 0);
165  _GLIBCXX_PARALLEL_ASSERT(rank < N);
166 
167  difference_type* ns = new difference_type[m];
168  difference_type* a = new difference_type[m];
169  difference_type* b = new difference_type[m];
170  difference_type l;
171 
172  ns[0] = std::distance(begin_seqs[0].first, begin_seqs[0].second);
173  nmax = ns[0];
174  for (int i = 0; i < m; i++)
175  {
176  ns[i] = std::distance(begin_seqs[i].first, begin_seqs[i].second);
177  nmax = std::max(nmax, ns[i]);
178  }
179 
180  r = __log2(nmax) + 1;
181 
182  // Pad all lists to this length, at least as long as any ns[i],
183  // equality iff nmax = 2^k - 1.
184  l = (1ULL << r) - 1;
185 
186  for (int i = 0; i < m; i++)
187  {
188  a[i] = 0;
189  b[i] = l;
190  }
191  n = l / 2;
192 
193  // Invariants:
194  // 0 <= a[i] <= ns[i], 0 <= b[i] <= l
195 
196 #define S(i) (begin_seqs[i].first)
197 
198  // Initial partition.
200 
201  for (int i = 0; i < m; i++)
202  if (n < ns[i]) //sequence long enough
203  sample.push_back(std::make_pair(S(i)[n], i));
204  __gnu_sequential::sort(sample.begin(), sample.end(), lcomp);
205 
206  for (int i = 0; i < m; i++) //conceptual infinity
207  if (n >= ns[i]) //sequence too short, conceptual infinity
208  sample.push_back(std::make_pair(S(i)[0] /*dummy element*/, i));
209 
210  difference_type localrank = rank / l;
211 
212  int j;
213  for (j = 0; j < localrank && ((n + 1) <= ns[sample[j].second]); ++j)
214  a[sample[j].second] += n + 1;
215  for (; j < m; j++)
216  b[sample[j].second] -= n + 1;
217 
218  // Further refinement.
219  while (n > 0)
220  {
221  n /= 2;
222 
223  int lmax_seq = -1; // to avoid warning
224  const value_type* lmax = NULL; // impossible to avoid the warning?
225  for (int i = 0; i < m; i++)
226  {
227  if (a[i] > 0)
228  {
229  if (!lmax)
230  {
231  lmax = &(S(i)[a[i] - 1]);
232  lmax_seq = i;
233  }
234  else
235  {
236  // Max, favor rear sequences.
237  if (!comp(S(i)[a[i] - 1], *lmax))
238  {
239  lmax = &(S(i)[a[i] - 1]);
240  lmax_seq = i;
241  }
242  }
243  }
244  }
245 
246  int i;
247  for (i = 0; i < m; i++)
248  {
249  difference_type middle = (b[i] + a[i]) / 2;
250  if (lmax && middle < ns[i] &&
251  lcomp(std::make_pair(S(i)[middle], i),
252  std::make_pair(*lmax, lmax_seq)))
253  a[i] = std::min(a[i] + n + 1, ns[i]);
254  else
255  b[i] -= n + 1;
256  }
257 
258  difference_type leftsize = 0;
259  for (int i = 0; i < m; i++)
260  leftsize += a[i] / (n + 1);
261 
262  difference_type skew = rank / (n + 1) - leftsize;
263 
264  if (skew > 0)
265  {
266  // Move to the left, find smallest.
270  pq(lrcomp);
271 
272  for (int i = 0; i < m; i++)
273  if (b[i] < ns[i])
274  pq.push(std::make_pair(S(i)[b[i]], i));
275 
276  for (; skew != 0 && !pq.empty(); --skew)
277  {
278  int source = pq.top().second;
279  pq.pop();
280 
281  a[source] = std::min(a[source] + n + 1, ns[source]);
282  b[source] += n + 1;
283 
284  if (b[source] < ns[source])
285  pq.push(std::make_pair(S(source)[b[source]], source));
286  }
287  }
288  else if (skew < 0)
289  {
290  // Move to the right, find greatest.
294 
295  for (int i = 0; i < m; i++)
296  if (a[i] > 0)
297  pq.push(std::make_pair(S(i)[a[i] - 1], i));
298 
299  for (; skew != 0; ++skew)
300  {
301  int source = pq.top().second;
302  pq.pop();
303 
304  a[source] -= n + 1;
305  b[source] -= n + 1;
306 
307  if (a[source] > 0)
308  pq.push(std::make_pair(S(source)[a[source] - 1], source));
309  }
310  }
311  }
312 
313  // Postconditions:
314  // a[i] == b[i] in most cases, except when a[i] has been clamped
315  // because of having reached the boundary
316 
317  // Now return the result, calculate the offset.
318 
319  // Compare the keys on both edges of the border.
320 
321  // Maximum of left edge, minimum of right edge.
322  value_type* maxleft = NULL;
323  value_type* minright = NULL;
324  for (int i = 0; i < m; i++)
325  {
326  if (a[i] > 0)
327  {
328  if (!maxleft)
329  maxleft = &(S(i)[a[i] - 1]);
330  else
331  {
332  // Max, favor rear sequences.
333  if (!comp(S(i)[a[i] - 1], *maxleft))
334  maxleft = &(S(i)[a[i] - 1]);
335  }
336  }
337  if (b[i] < ns[i])
338  {
339  if (!minright)
340  minright = &(S(i)[b[i]]);
341  else
342  {
343  // Min, favor fore sequences.
344  if (comp(S(i)[b[i]], *minright))
345  minright = &(S(i)[b[i]]);
346  }
347  }
348  }
349 
350  int seq = 0;
351  for (int i = 0; i < m; i++)
352  begin_offsets[i] = S(i) + a[i];
353 
354  delete[] ns;
355  delete[] a;
356  delete[] b;
357  }
358 
359 
360  /**
361  * @brief Selects the element at a certain global rank from several
362  * sorted sequences.
363  *
364  * The sequences are passed via a sequence of random-access
365  * iterator pairs, none of the sequences may be empty.
366  * @param begin_seqs Begin of the sequence of iterator pairs.
367  * @param end_seqs End of the sequence of iterator pairs.
368  * @param rank The global rank to partition at.
369  * @param offset The rank of the selected element in the global
370  * subsequence of elements equal to the selected element. If the
371  * selected element is unique, this number is 0.
372  * @param comp The ordering functor, defaults to std::less.
373  */
374  template<typename T, typename RanSeqs, typename RankType,
375  typename Comparator>
376  T
377  multiseq_selection(RanSeqs begin_seqs, RanSeqs end_seqs, RankType rank,
378  RankType& offset, Comparator comp = std::less<T>())
379  {
380  _GLIBCXX_CALL(end_seqs - begin_seqs)
381 
383  It;
384  typedef typename std::iterator_traits<It>::difference_type
385  difference_type;
386 
389 
390  // Number of sequences, number of elements in total (possibly
391  // including padding).
392  difference_type m = std::distance(begin_seqs, end_seqs);
393  difference_type N = 0;
394  difference_type nmax, n, r;
395 
396  for (int i = 0; i < m; i++)
397  N += std::distance(begin_seqs[i].first, begin_seqs[i].second);
398 
399  if (m == 0 || N == 0 || rank < 0 || rank >= N)
400  {
401  // Result undefined when there is no data or rank is outside bounds.
402  throw std::exception();
403  }
404 
405 
406  difference_type* ns = new difference_type[m];
407  difference_type* a = new difference_type[m];
408  difference_type* b = new difference_type[m];
409  difference_type l;
410 
411  ns[0] = std::distance(begin_seqs[0].first, begin_seqs[0].second);
412  nmax = ns[0];
413  for (int i = 0; i < m; ++i)
414  {
415  ns[i] = std::distance(begin_seqs[i].first, begin_seqs[i].second);
416  nmax = std::max(nmax, ns[i]);
417  }
418 
419  r = __log2(nmax) + 1;
420 
421  // Pad all lists to this length, at least as long as any ns[i],
422  // equality iff nmax = 2^k - 1
423  l = pow2(r) - 1;
424 
425  for (int i = 0; i < m; ++i)
426  {
427  a[i] = 0;
428  b[i] = l;
429  }
430  n = l / 2;
431 
432  // Invariants:
433  // 0 <= a[i] <= ns[i], 0 <= b[i] <= l
434 
435 #define S(i) (begin_seqs[i].first)
436 
437  // Initial partition.
439 
440  for (int i = 0; i < m; i++)
441  if (n < ns[i])
442  sample.push_back(std::make_pair(S(i)[n], i));
443  __gnu_sequential::sort(sample.begin(), sample.end(),
444  lcomp, sequential_tag());
445 
446  // Conceptual infinity.
447  for (int i = 0; i < m; i++)
448  if (n >= ns[i])
449  sample.push_back(std::make_pair(S(i)[0] /*dummy element*/, i));
450 
451  difference_type localrank = rank / l;
452 
453  int j;
454  for (j = 0; j < localrank && ((n + 1) <= ns[sample[j].second]); ++j)
455  a[sample[j].second] += n + 1;
456  for (; j < m; ++j)
457  b[sample[j].second] -= n + 1;
458 
459  // Further refinement.
460  while (n > 0)
461  {
462  n /= 2;
463 
464  const T* lmax = NULL;
465  for (int i = 0; i < m; ++i)
466  {
467  if (a[i] > 0)
468  {
469  if (!lmax)
470  lmax = &(S(i)[a[i] - 1]);
471  else
472  {
473  if (comp(*lmax, S(i)[a[i] - 1])) //max
474  lmax = &(S(i)[a[i] - 1]);
475  }
476  }
477  }
478 
479  int i;
480  for (i = 0; i < m; i++)
481  {
482  difference_type middle = (b[i] + a[i]) / 2;
483  if (lmax && middle < ns[i] && comp(S(i)[middle], *lmax))
484  a[i] = std::min(a[i] + n + 1, ns[i]);
485  else
486  b[i] -= n + 1;
487  }
488 
489  difference_type leftsize = 0;
490  for (int i = 0; i < m; ++i)
491  leftsize += a[i] / (n + 1);
492 
493  difference_type skew = rank / (n + 1) - leftsize;
494 
495  if (skew > 0)
496  {
497  // Move to the left, find smallest.
501 
502  for (int i = 0; i < m; ++i)
503  if (b[i] < ns[i])
504  pq.push(std::make_pair(S(i)[b[i]], i));
505 
506  for (; skew != 0 && !pq.empty(); --skew)
507  {
508  int source = pq.top().second;
509  pq.pop();
510 
511  a[source] = std::min(a[source] + n + 1, ns[source]);
512  b[source] += n + 1;
513 
514  if (b[source] < ns[source])
515  pq.push(std::make_pair(S(source)[b[source]], source));
516  }
517  }
518  else if (skew < 0)
519  {
520  // Move to the right, find greatest.
524 
525  for (int i = 0; i < m; ++i)
526  if (a[i] > 0)
527  pq.push(std::make_pair(S(i)[a[i] - 1], i));
528 
529  for (; skew != 0; ++skew)
530  {
531  int source = pq.top().second;
532  pq.pop();
533 
534  a[source] -= n + 1;
535  b[source] -= n + 1;
536 
537  if (a[source] > 0)
538  pq.push(std::make_pair(S(source)[a[source] - 1], source));
539  }
540  }
541  }
542 
543  // Postconditions:
544  // a[i] == b[i] in most cases, except when a[i] has been clamped
545  // because of having reached the boundary
546 
547  // Now return the result, calculate the offset.
548 
549  // Compare the keys on both edges of the border.
550 
551  // Maximum of left edge, minimum of right edge.
552  bool maxleftset = false, minrightset = false;
553 
554  // Impossible to avoid the warning?
555  T maxleft, minright;
556  for (int i = 0; i < m; ++i)
557  {
558  if (a[i] > 0)
559  {
560  if (!maxleftset)
561  {
562  maxleft = S(i)[a[i] - 1];
563  maxleftset = true;
564  }
565  else
566  {
567  // Max.
568  if (comp(maxleft, S(i)[a[i] - 1]))
569  maxleft = S(i)[a[i] - 1];
570  }
571  }
572  if (b[i] < ns[i])
573  {
574  if (!minrightset)
575  {
576  minright = S(i)[b[i]];
577  minrightset = true;
578  }
579  else
580  {
581  // Min.
582  if (comp(S(i)[b[i]], minright))
583  minright = S(i)[b[i]];
584  }
585  }
586  }
587 
588  // Minright is the splitter, in any case.
589 
590  if (!maxleftset || comp(minright, maxleft))
591  {
592  // Good luck, everything is split unambiguously.
593  offset = 0;
594  }
595  else
596  {
597  // We have to calculate an offset.
598  offset = 0;
599 
600  for (int i = 0; i < m; ++i)
601  {
602  difference_type lb = std::lower_bound(S(i), S(i) + ns[i],
603  minright,
604  comp) - S(i);
605  offset += a[i] - lb;
606  }
607  }
608 
609  delete[] ns;
610  delete[] a;
611  delete[] b;
612 
613  return minright;
614  }
615 }
616 
617 #undef S
618 
619 #endif /* _GLIBCXX_PARALLEL_MULTISEQ_SELECTION_H */
#define _GLIBCXX_CALL(n)
Macro to produce log message when entering a function.
void multiseq_partition(RanSeqs begin_seqs, RanSeqs end_seqs, RankType rank, RankIterator begin_offsets, Comparator comp=std::less< typename std::iterator_traits< typename std::iterator_traits< RanSeqs >::value_type::first_type >::value_type >())
Splits several sorted sequences at a certain global rank, resulting in a splitting point for each seq...
void push_back(const value_type &__x)
Add data to the end of the vector.
Definition: stl_vector.h:733
const_reference top() const
Definition: stl_queue.h:494
T multiseq_selection(RanSeqs begin_seqs, RanSeqs end_seqs, RankType rank, RankType &offset, Comparator comp=std::less< T >())
Selects the element at a certain global rank from several sorted sequences.
void pop()
Removes first element.
Definition: stl_queue.h:544
Size __log2(Size n)
Calculates the rounded-down logarithm of n for base 2.
Definition: base.h:105
_T1 first
first is a copy of the first object
Definition: stl_pair.h:72
void push(const value_type &__x)
Add data to the queue.
Definition: stl_queue.h:509
Compare a pair of types lexicographically, descending.
A standard container automatically sorting its contents.
Definition: stl_queue.h:369
One of the comparison functors.
Definition: stl_function.h:226
Forces sequential execution at compile time.
Definition: tags.h:42
GNU parallel code for public use.
_T2 second
second is a copy of the second object
Definition: stl_pair.h:73
pair holds two objects of arbitrary type.
Definition: stl_pair.h:67
Compare a pair of types lexicographically, ascending.
A standard container which offers fixed time access to individual elements in any order...
Definition: stl_vector.h:170
const _Tp & max(const _Tp &, const _Tp &)
This does what you think it does.
Definition: stl_algobase.h:209
bool empty() const
Definition: stl_queue.h:481
const _Tp & min(const _Tp &, const _Tp &)
This does what you think it does.
Definition: stl_algobase.h:186
iterator_traits< _InputIterator >::difference_type distance(_InputIterator __first, _InputIterator __last)
A generalization of pointer arithmetic.
Parallel sorting algorithm switch. This file is a GNU parallel extension to the Standard C++ Library...
iterator end()
Definition: stl_vector.h:443
iterator begin()
Definition: stl_vector.h:425