libstdc++
multiseq_selection.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 
3 // Copyright (C) 2007-2021 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 namespace __gnu_parallel
50 {
51  /** @brief Compare __a pair of types lexicographically, ascending. */
52  template<typename _T1, typename _T2, typename _Compare>
54  : public std::binary_function<std::pair<_T1, _T2>,
55  std::pair<_T1, _T2>, bool>
56  {
57  private:
58  _Compare& _M_comp;
59 
60  public:
61  _Lexicographic(_Compare& __comp) : _M_comp(__comp) { }
62 
63  bool
64  operator()(const std::pair<_T1, _T2>& __p1,
65  const std::pair<_T1, _T2>& __p2) const
66  {
67  if (_M_comp(__p1.first, __p2.first))
68  return true;
69 
70  if (_M_comp(__p2.first, __p1.first))
71  return false;
72 
73  // Firsts are equal.
74  return __p1.second < __p2.second;
75  }
76  };
77 
78  /** @brief Compare __a pair of types lexicographically, descending. */
79  template<typename _T1, typename _T2, typename _Compare>
80  class _LexicographicReverse : public std::binary_function<_T1, _T2, bool>
81  {
82  private:
83  _Compare& _M_comp;
84 
85  public:
86  _LexicographicReverse(_Compare& __comp) : _M_comp(__comp) { }
87 
88  bool
89  operator()(const std::pair<_T1, _T2>& __p1,
90  const std::pair<_T1, _T2>& __p2) const
91  {
92  if (_M_comp(__p2.first, __p1.first))
93  return true;
94 
95  if (_M_comp(__p1.first, __p2.first))
96  return false;
97 
98  // Firsts are equal.
99  return __p2.second < __p1.second;
100  }
101  };
102 
103  /**
104  * @brief Splits several sorted sequences at a certain global __rank,
105  * resulting in a splitting point for each sequence.
106  * The sequences are passed via a sequence of random-access
107  * iterator pairs, none of the sequences may be empty. If there
108  * are several equal elements across the split, the ones on the
109  * __left side will be chosen from sequences with smaller number.
110  * @param __begin_seqs Begin of the sequence of iterator pairs.
111  * @param __end_seqs End of the sequence of iterator pairs.
112  * @param __rank The global rank to partition at.
113  * @param __begin_offsets A random-access __sequence __begin where the
114  * __result will be stored in. Each element of the sequence is an
115  * iterator that points to the first element on the greater part of
116  * the respective __sequence.
117  * @param __comp The ordering functor, defaults to std::less<_Tp>.
118  */
119  template<typename _RanSeqs, typename _RankType, typename _RankIterator,
120  typename _Compare>
121  void
122  multiseq_partition(_RanSeqs __begin_seqs, _RanSeqs __end_seqs,
123  _RankType __rank,
124  _RankIterator __begin_offsets,
125  _Compare __comp = std::less<
126  typename std::iterator_traits<typename
128  first_type>::value_type>()) // std::less<_Tp>
129  {
130  _GLIBCXX_CALL(__end_seqs - __begin_seqs)
131 
133  _It;
135  _SeqNumber;
137  _DifferenceType;
138  typedef typename std::iterator_traits<_It>::value_type _ValueType;
139 
142 
143  // Number of sequences, number of elements in total (possibly
144  // including padding).
145  _DifferenceType __m = std::distance(__begin_seqs, __end_seqs), __nn = 0,
146  __nmax, __n, __r;
147 
148  for (_SeqNumber __i = 0; __i < __m; __i++)
149  {
150  __nn += std::distance(__begin_seqs[__i].first,
151  __begin_seqs[__i].second);
152  _GLIBCXX_PARALLEL_ASSERT(
153  std::distance(__begin_seqs[__i].first,
154  __begin_seqs[__i].second) > 0);
155  }
156 
157  if (__rank == __nn)
158  {
159  for (_SeqNumber __i = 0; __i < __m; __i++)
160  __begin_offsets[__i] = __begin_seqs[__i].second; // Very end.
161  // Return __m - 1;
162  return;
163  }
164 
165  _GLIBCXX_PARALLEL_ASSERT(__m != 0);
166  _GLIBCXX_PARALLEL_ASSERT(__nn != 0);
167  _GLIBCXX_PARALLEL_ASSERT(__rank >= 0);
168  _GLIBCXX_PARALLEL_ASSERT(__rank < __nn);
169 
170  _DifferenceType* __ns = new _DifferenceType[__m];
171  _DifferenceType* __a = new _DifferenceType[__m];
172  _DifferenceType* __b = new _DifferenceType[__m];
173  _DifferenceType __l;
174 
175  __ns[0] = std::distance(__begin_seqs[0].first, __begin_seqs[0].second);
176  __nmax = __ns[0];
177  for (_SeqNumber __i = 0; __i < __m; __i++)
178  {
179  __ns[__i] = std::distance(__begin_seqs[__i].first,
180  __begin_seqs[__i].second);
181  __nmax = std::max(__nmax, __ns[__i]);
182  }
183 
184  __r = __rd_log2(__nmax) + 1;
185 
186  // Pad all lists to this length, at least as long as any ns[__i],
187  // equality iff __nmax = 2^__k - 1.
188  __l = (1ULL << __r) - 1;
189 
190  for (_SeqNumber __i = 0; __i < __m; __i++)
191  {
192  __a[__i] = 0;
193  __b[__i] = __l;
194  }
195  __n = __l / 2;
196 
197  // Invariants:
198  // 0 <= __a[__i] <= __ns[__i], 0 <= __b[__i] <= __l
199 
200 #define __S(__i) (__begin_seqs[__i].first)
201 
202  // Initial partition.
204 
205  for (_SeqNumber __i = 0; __i < __m; __i++)
206  if (__n < __ns[__i]) //__sequence long enough
207  __sample.push_back(std::make_pair(__S(__i)[__n], __i));
208  __gnu_sequential::sort(__sample.begin(), __sample.end(), __lcomp);
209 
210  for (_SeqNumber __i = 0; __i < __m; __i++) //conceptual infinity
211  if (__n >= __ns[__i]) //__sequence too short, conceptual infinity
212  __sample.push_back(
213  std::make_pair(__S(__i)[0] /*__dummy element*/, __i));
214 
215  _DifferenceType __localrank = __rank / __l;
216 
217  _SeqNumber __j;
218  for (__j = 0;
219  __j < __localrank && ((__n + 1) <= __ns[__sample[__j].second]);
220  ++__j)
221  __a[__sample[__j].second] += __n + 1;
222  for (; __j < __m; __j++)
223  __b[__sample[__j].second] -= __n + 1;
224 
225  // Further refinement.
226  while (__n > 0)
227  {
228  __n /= 2;
229 
230  _SeqNumber __lmax_seq = -1; // to avoid warning
231  const _ValueType* __lmax = 0; // impossible to avoid the warning?
232  for (_SeqNumber __i = 0; __i < __m; __i++)
233  {
234  if (__a[__i] > 0)
235  {
236  if (!__lmax)
237  {
238  __lmax = &(__S(__i)[__a[__i] - 1]);
239  __lmax_seq = __i;
240  }
241  else
242  {
243  // Max, favor rear sequences.
244  if (!__comp(__S(__i)[__a[__i] - 1], *__lmax))
245  {
246  __lmax = &(__S(__i)[__a[__i] - 1]);
247  __lmax_seq = __i;
248  }
249  }
250  }
251  }
252 
253  _SeqNumber __i;
254  for (__i = 0; __i < __m; __i++)
255  {
256  _DifferenceType __middle = (__b[__i] + __a[__i]) / 2;
257  if (__lmax && __middle < __ns[__i] &&
258  __lcomp(std::make_pair(__S(__i)[__middle], __i),
259  std::make_pair(*__lmax, __lmax_seq)))
260  __a[__i] = std::min(__a[__i] + __n + 1, __ns[__i]);
261  else
262  __b[__i] -= __n + 1;
263  }
264 
265  _DifferenceType __leftsize = 0;
266  for (_SeqNumber __i = 0; __i < __m; __i++)
267  __leftsize += __a[__i] / (__n + 1);
268 
269  _DifferenceType __skew = __rank / (__n + 1) - __leftsize;
270 
271  if (__skew > 0)
272  {
273  // Move to the left, find smallest.
277  __pq(__lrcomp);
278 
279  for (_SeqNumber __i = 0; __i < __m; __i++)
280  if (__b[__i] < __ns[__i])
281  __pq.push(std::make_pair(__S(__i)[__b[__i]], __i));
282 
283  for (; __skew != 0 && !__pq.empty(); --__skew)
284  {
285  _SeqNumber __source = __pq.top().second;
286  __pq.pop();
287 
288  __a[__source]
289  = std::min(__a[__source] + __n + 1, __ns[__source]);
290  __b[__source] += __n + 1;
291 
292  if (__b[__source] < __ns[__source])
293  __pq.push(
294  std::make_pair(__S(__source)[__b[__source]], __source));
295  }
296  }
297  else if (__skew < 0)
298  {
299  // Move to the right, find greatest.
303  __pq(__lcomp);
304 
305  for (_SeqNumber __i = 0; __i < __m; __i++)
306  if (__a[__i] > 0)
307  __pq.push(std::make_pair(__S(__i)[__a[__i] - 1], __i));
308 
309  for (; __skew != 0; ++__skew)
310  {
311  _SeqNumber __source = __pq.top().second;
312  __pq.pop();
313 
314  __a[__source] -= __n + 1;
315  __b[__source] -= __n + 1;
316 
317  if (__a[__source] > 0)
318  __pq.push(std::make_pair(
319  __S(__source)[__a[__source] - 1], __source));
320  }
321  }
322  }
323 
324  // Postconditions:
325  // __a[__i] == __b[__i] in most cases, except when __a[__i] has been
326  // clamped because of having reached the boundary
327 
328  // Now return the result, calculate the offset.
329 
330  // Compare the keys on both edges of the border.
331 
332  // Maximum of left edge, minimum of right edge.
333  _ValueType* __maxleft = 0;
334  _ValueType* __minright = 0;
335  for (_SeqNumber __i = 0; __i < __m; __i++)
336  {
337  if (__a[__i] > 0)
338  {
339  if (!__maxleft)
340  __maxleft = &(__S(__i)[__a[__i] - 1]);
341  else
342  {
343  // Max, favor rear sequences.
344  if (!__comp(__S(__i)[__a[__i] - 1], *__maxleft))
345  __maxleft = &(__S(__i)[__a[__i] - 1]);
346  }
347  }
348  if (__b[__i] < __ns[__i])
349  {
350  if (!__minright)
351  __minright = &(__S(__i)[__b[__i]]);
352  else
353  {
354  // Min, favor fore sequences.
355  if (__comp(__S(__i)[__b[__i]], *__minright))
356  __minright = &(__S(__i)[__b[__i]]);
357  }
358  }
359  }
360 
361  _SeqNumber __seq = 0;
362  for (_SeqNumber __i = 0; __i < __m; __i++)
363  __begin_offsets[__i] = __S(__i) + __a[__i];
364 
365  delete[] __ns;
366  delete[] __a;
367  delete[] __b;
368  }
369 
370 
371  /**
372  * @brief Selects the element at a certain global __rank from several
373  * sorted sequences.
374  *
375  * The sequences are passed via a sequence of random-access
376  * iterator pairs, none of the sequences may be empty.
377  * @param __begin_seqs Begin of the sequence of iterator pairs.
378  * @param __end_seqs End of the sequence of iterator pairs.
379  * @param __rank The global rank to partition at.
380  * @param __offset The rank of the selected element in the global
381  * subsequence of elements equal to the selected element. If the
382  * selected element is unique, this number is 0.
383  * @param __comp The ordering functor, defaults to std::less.
384  */
385  template<typename _Tp, typename _RanSeqs, typename _RankType,
386  typename _Compare>
387  _Tp
388  multiseq_selection(_RanSeqs __begin_seqs, _RanSeqs __end_seqs,
389  _RankType __rank,
390  _RankType& __offset, _Compare __comp = std::less<_Tp>())
391  {
392  _GLIBCXX_CALL(__end_seqs - __begin_seqs)
393 
395  _It;
397  _SeqNumber;
399  _DifferenceType;
400 
403 
404  // Number of sequences, number of elements in total (possibly
405  // including padding).
406  _DifferenceType __m = std::distance(__begin_seqs, __end_seqs);
407  _DifferenceType __nn = 0;
408  _DifferenceType __nmax, __n, __r;
409 
410  for (_SeqNumber __i = 0; __i < __m; __i++)
411  __nn += std::distance(__begin_seqs[__i].first,
412  __begin_seqs[__i].second);
413 
414  if (__m == 0 || __nn == 0 || __rank < 0 || __rank >= __nn)
415  {
416  // result undefined if there is no data or __rank is outside bounds
417  throw std::exception();
418  }
419 
420 
421  _DifferenceType* __ns = new _DifferenceType[__m];
422  _DifferenceType* __a = new _DifferenceType[__m];
423  _DifferenceType* __b = new _DifferenceType[__m];
424  _DifferenceType __l;
425 
426  __ns[0] = std::distance(__begin_seqs[0].first, __begin_seqs[0].second);
427  __nmax = __ns[0];
428  for (_SeqNumber __i = 0; __i < __m; ++__i)
429  {
430  __ns[__i] = std::distance(__begin_seqs[__i].first,
431  __begin_seqs[__i].second);
432  __nmax = std::max(__nmax, __ns[__i]);
433  }
434 
435  __r = __rd_log2(__nmax) + 1;
436 
437  // Pad all lists to this length, at least as long as any ns[__i],
438  // equality iff __nmax = 2^__k - 1
439  __l = __round_up_to_pow2(__r) - 1;
440 
441  for (_SeqNumber __i = 0; __i < __m; ++__i)
442  {
443  __a[__i] = 0;
444  __b[__i] = __l;
445  }
446  __n = __l / 2;
447 
448  // Invariants:
449  // 0 <= __a[__i] <= __ns[__i], 0 <= __b[__i] <= __l
450 
451 #define __S(__i) (__begin_seqs[__i].first)
452 
453  // Initial partition.
455 
456  for (_SeqNumber __i = 0; __i < __m; __i++)
457  if (__n < __ns[__i])
458  __sample.push_back(std::make_pair(__S(__i)[__n], __i));
459  __gnu_sequential::sort(__sample.begin(), __sample.end(),
460  __lcomp, sequential_tag());
461 
462  // Conceptual infinity.
463  for (_SeqNumber __i = 0; __i < __m; __i++)
464  if (__n >= __ns[__i])
465  __sample.push_back(
466  std::make_pair(__S(__i)[0] /*__dummy element*/, __i));
467 
468  _DifferenceType __localrank = __rank / __l;
469 
470  _SeqNumber __j;
471  for (__j = 0;
472  __j < __localrank && ((__n + 1) <= __ns[__sample[__j].second]);
473  ++__j)
474  __a[__sample[__j].second] += __n + 1;
475  for (; __j < __m; ++__j)
476  __b[__sample[__j].second] -= __n + 1;
477 
478  // Further refinement.
479  while (__n > 0)
480  {
481  __n /= 2;
482 
483  const _Tp* __lmax = 0;
484  for (_SeqNumber __i = 0; __i < __m; ++__i)
485  {
486  if (__a[__i] > 0)
487  {
488  if (!__lmax)
489  __lmax = &(__S(__i)[__a[__i] - 1]);
490  else
491  {
492  if (__comp(*__lmax, __S(__i)[__a[__i] - 1])) //max
493  __lmax = &(__S(__i)[__a[__i] - 1]);
494  }
495  }
496  }
497 
498  _SeqNumber __i;
499  for (__i = 0; __i < __m; __i++)
500  {
501  _DifferenceType __middle = (__b[__i] + __a[__i]) / 2;
502  if (__lmax && __middle < __ns[__i]
503  && __comp(__S(__i)[__middle], *__lmax))
504  __a[__i] = std::min(__a[__i] + __n + 1, __ns[__i]);
505  else
506  __b[__i] -= __n + 1;
507  }
508 
509  _DifferenceType __leftsize = 0;
510  for (_SeqNumber __i = 0; __i < __m; ++__i)
511  __leftsize += __a[__i] / (__n + 1);
512 
513  _DifferenceType __skew = __rank / (__n + 1) - __leftsize;
514 
515  if (__skew > 0)
516  {
517  // Move to the left, find smallest.
521  __pq(__lrcomp);
522 
523  for (_SeqNumber __i = 0; __i < __m; ++__i)
524  if (__b[__i] < __ns[__i])
525  __pq.push(std::make_pair(__S(__i)[__b[__i]], __i));
526 
527  for (; __skew != 0 && !__pq.empty(); --__skew)
528  {
529  _SeqNumber __source = __pq.top().second;
530  __pq.pop();
531 
532  __a[__source]
533  = std::min(__a[__source] + __n + 1, __ns[__source]);
534  __b[__source] += __n + 1;
535 
536  if (__b[__source] < __ns[__source])
537  __pq.push(
538  std::make_pair(__S(__source)[__b[__source]], __source));
539  }
540  }
541  else if (__skew < 0)
542  {
543  // Move to the right, find greatest.
547 
548  for (_SeqNumber __i = 0; __i < __m; ++__i)
549  if (__a[__i] > 0)
550  __pq.push(std::make_pair(__S(__i)[__a[__i] - 1], __i));
551 
552  for (; __skew != 0; ++__skew)
553  {
554  _SeqNumber __source = __pq.top().second;
555  __pq.pop();
556 
557  __a[__source] -= __n + 1;
558  __b[__source] -= __n + 1;
559 
560  if (__a[__source] > 0)
561  __pq.push(std::make_pair(
562  __S(__source)[__a[__source] - 1], __source));
563  }
564  }
565  }
566 
567  // Postconditions:
568  // __a[__i] == __b[__i] in most cases, except when __a[__i] has been
569  // clamped because of having reached the boundary
570 
571  // Now return the result, calculate the offset.
572 
573  // Compare the keys on both edges of the border.
574 
575  // Maximum of left edge, minimum of right edge.
576  bool __maxleftset = false, __minrightset = false;
577 
578  // Impossible to avoid the warning?
579  _Tp __maxleft, __minright;
580  for (_SeqNumber __i = 0; __i < __m; ++__i)
581  {
582  if (__a[__i] > 0)
583  {
584  if (!__maxleftset)
585  {
586  __maxleft = __S(__i)[__a[__i] - 1];
587  __maxleftset = true;
588  }
589  else
590  {
591  // Max.
592  if (__comp(__maxleft, __S(__i)[__a[__i] - 1]))
593  __maxleft = __S(__i)[__a[__i] - 1];
594  }
595  }
596  if (__b[__i] < __ns[__i])
597  {
598  if (!__minrightset)
599  {
600  __minright = __S(__i)[__b[__i]];
601  __minrightset = true;
602  }
603  else
604  {
605  // Min.
606  if (__comp(__S(__i)[__b[__i]], __minright))
607  __minright = __S(__i)[__b[__i]];
608  }
609  }
610  }
611 
612  // Minright is the __splitter, in any case.
613 
614  if (!__maxleftset || __comp(__minright, __maxleft))
615  {
616  // Good luck, everything is split unambiguously.
617  __offset = 0;
618  }
619  else
620  {
621  // We have to calculate an offset.
622  __offset = 0;
623 
624  for (_SeqNumber __i = 0; __i < __m; ++__i)
625  {
626  _DifferenceType lb
627  = std::lower_bound(__S(__i), __S(__i) + __ns[__i],
628  __minright,
629  __comp) - __S(__i);
630  __offset += __a[__i] - lb;
631  }
632  }
633 
634  delete[] __ns;
635  delete[] __a;
636  delete[] __b;
637 
638  return __minright;
639  }
640 }
641 
642 #undef __S
643 
644 #endif /* _GLIBCXX_PARALLEL_MULTISEQ_SELECTION_H */
#define _GLIBCXX_CALL(__n)
Macro to produce log message when entering a function.
constexpr const _Tp & max(const _Tp &, const _Tp &)
This does what you think it does.
Definition: stl_algobase.h:254
constexpr const _Tp & min(const _Tp &, const _Tp &)
This does what you think it does.
Definition: stl_algobase.h:230
_RandomAccessIterator __sample(_InputIterator __first, _InputIterator __last, input_iterator_tag, _RandomAccessIterator __out, random_access_iterator_tag, _Size __n, _UniformRandomBitGenerator &&__g)
Reservoir sampling algorithm.
Definition: stl_algo.h:5748
GNU parallel code for public use.
_Tp multiseq_selection(_RanSeqs __begin_seqs, _RanSeqs __end_seqs, _RankType __rank, _RankType &__offset, _Compare __comp=std::less< _Tp >())
Selects the element at a certain global __rank from several sorted sequences.
_Tp __round_up_to_pow2(_Tp __x)
Round up to the next greater power of 2.
void multiseq_partition(_RanSeqs __begin_seqs, _RanSeqs __end_seqs, _RankType __rank, _RankIterator __begin_offsets, _Compare __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 s...
_Size __rd_log2(_Size __n)
Calculates the rounded-down logarithm of __n for base 2.
Definition: base.h:102
Traits class for iterators.
Base class for all library exceptions.
Definition: exception.h:61
One of the comparison functors.
Definition: stl_function.h:382
Struct holding two objects of arbitrary type.
Definition: stl_pair.h:202
_T1 first
The first member.
Definition: stl_pair.h:206
_T2 second
The second member.
Definition: stl_pair.h:207
A standard container automatically sorting its contents.
Definition: stl_queue.h:457
bool empty() const
Definition: stl_queue.h:612
void pop()
Removes first element.
Definition: stl_queue.h:675
const_reference top() const
Definition: stl_queue.h:625
void push(const value_type &__x)
Add data to the queue.
Definition: stl_queue.h:640
A standard container which offers fixed time access to individual elements in any order.
Definition: stl_vector.h:390
Compare __a pair of types lexicographically, ascending.
Compare __a pair of types lexicographically, descending.
Forces sequential execution at compile time.
Definition: tags.h:42