1 // -*- C++ -*-
2
3 // Copyright (C) 2007-2023 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>
53 class _Lexicographic
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
127 std::iterator_traits<_RanSeqs>::value_type::
128 first_type>::value_type>()) // std::less<_Tp>
129 {
130 _GLIBCXX_CALL(__end_seqs - __begin_seqs)
131
132 typedef typename std::iterator_traits<_RanSeqs>::value_type::first_type
133 _It;
134 typedef typename std::iterator_traits<_RanSeqs>::difference_type
135 _SeqNumber;
136 typedef typename std::iterator_traits<_It>::difference_type
137 _DifferenceType;
138 typedef typename std::iterator_traits<_It>::value_type _ValueType;
139
140 _Lexicographic<_ValueType, _SeqNumber, _Compare> __lcomp(__comp);
141 _LexicographicReverse<_ValueType, _SeqNumber, _Compare> __lrcomp(__comp);
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.
203 std::vector<std::pair<_ValueType, _SeqNumber> > __sample;
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.
274 std::priority_queue<std::pair<_ValueType, _SeqNumber>,
275 std::vector<std::pair<_ValueType, _SeqNumber> >,
276 _LexicographicReverse<_ValueType, _SeqNumber, _Compare> >
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.
300 std::priority_queue<std::pair<_ValueType, _SeqNumber>,
301 std::vector<std::pair<_ValueType, _SeqNumber> >,
302 _Lexicographic<_ValueType, _SeqNumber, _Compare> >
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
394 typedef typename std::iterator_traits<_RanSeqs>::value_type::first_type
395 _It;
396 typedef typename std::iterator_traits<_RanSeqs>::difference_type
397 _SeqNumber;
398 typedef typename std::iterator_traits<_It>::difference_type
399 _DifferenceType;
400
401 _Lexicographic<_Tp, _SeqNumber, _Compare> __lcomp(__comp);
402 _LexicographicReverse<_Tp, _SeqNumber, _Compare> __lrcomp(__comp);
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.
454 std::vector<std::pair<_Tp, _SeqNumber> > __sample;
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.
518 std::priority_queue<std::pair<_Tp, _SeqNumber>,
519 std::vector<std::pair<_Tp, _SeqNumber> >,
520 _LexicographicReverse<_Tp, _SeqNumber, _Compare> >
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.
544 std::priority_queue<std::pair<_Tp, _SeqNumber>,
545 std::vector<std::pair<_Tp, _SeqNumber> >,
546 _Lexicographic<_Tp, _SeqNumber, _Compare> > __pq(__lcomp);
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 */