3 // Copyright (C) 2007, 2008, 2009 Free Software Foundation, Inc.
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
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.
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.
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/>.
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
30 * The algorithm description can be found in
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.
36 * This file is a GNU parallel extension to the Standard C++ Library.
39 // Written by Johannes Singler.
41 #ifndef _GLIBCXX_PARALLEL_MULTISEQ_SELECTION_H
42 #define _GLIBCXX_PARALLEL_MULTISEQ_SELECTION_H 1
47 #include <bits/stl_algo.h>
49 #include <parallel/sort.h>
51 namespace __gnu_parallel
53 /** @brief Compare __a pair of types lexicographically, ascending. */
54 template<typename _T1
, typename _T2
, typename _Compare
>
56 : public std::binary_function
<
57 std::pair
<_T1
, _T2
>, std::pair
<_T1
, _T2
>, bool>
63 _Lexicographic(_Compare
& __comp
) : _M_comp(__comp
) { }
66 operator()(const std::pair
<_T1
, _T2
>& __p1
,
67 const std::pair
<_T1
, _T2
>& __p2
) const
69 if (_M_comp(__p1
.first
, __p2
.first
))
72 if (_M_comp(__p2
.first
, __p1
.first
))
76 return __p1
.second
< __p2
.second
;
80 /** @brief Compare __a pair of types lexicographically, descending. */
81 template<typename _T1
, typename _T2
, typename _Compare
>
82 class _LexicographicReverse
: public std::binary_function
<_T1
, _T2
, bool>
88 _LexicographicReverse(_Compare
& __comp
) : _M_comp(__comp
) { }
91 operator()(const std::pair
<_T1
, _T2
>& __p1
,
92 const std::pair
<_T1
, _T2
>& __p2
) const
94 if (_M_comp(__p2
.first
, __p1
.first
))
97 if (_M_comp(__p1
.first
, __p2
.first
))
101 return __p2
.second
< __p1
.second
;
106 * @brief Splits several sorted sequences at __a certain global __rank,
107 * resulting in a splitting point for each sequence.
108 * The sequences are passed via __a __sequence of random-access
109 * iterator pairs, none of the sequences may be empty. If there
110 * are several equal elements across the split, the ones on the
111 * __left side will be chosen from sequences with smaller number.
112 * @param __begin_seqs Begin of the sequence of iterator pairs.
113 * @param __end_seqs End of the sequence of iterator pairs.
114 * @param __rank The global __rank to partition at.
115 * @param __begin_offsets A random-access __sequence __begin where the
116 * __result will be stored in. Each element of the sequence is an
117 * iterator that points to the first element on the greater part of
118 * the respective __sequence.
119 * @param __comp The ordering functor, defaults to std::less<_Tp>.
121 template<typename _RanSeqs
, typename _RankType
, typename _RankIterator
,
124 multiseq_partition(_RanSeqs __begin_seqs
, _RanSeqs __end_seqs
,
126 _RankIterator __begin_offsets
,
127 _Compare __comp
= std::less
<
128 typename
std::iterator_traits
<typename
129 std::iterator_traits
<_RanSeqs
>::value_type::
130 first_type
>::value_type
>()) // std::less<_Tp>
132 _GLIBCXX_CALL(__end_seqs
- __begin_seqs
)
134 typedef typename
std::iterator_traits
<_RanSeqs
>::value_type::first_type
136 typedef typename
std::iterator_traits
<_It
>::difference_type
138 typedef typename
std::iterator_traits
<_It
>::value_type _ValueType
;
140 _Lexicographic
<_ValueType
, int, _Compare
> __lcomp(__comp
);
141 _LexicographicReverse
<_ValueType
, int, _Compare
> __lrcomp(__comp
);
143 // Number of sequences, number of elements in total (possibly
144 // including padding).
145 _DifferenceType __m
= std::distance(__begin_seqs
, __end_seqs
), __N
= 0,
148 for (int __i
= 0; __i
< __m
; __i
++)
150 __N
+= 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);
159 for (int __i
= 0; __i
< __m
; __i
++)
160 __begin_offsets
[__i
] = __begin_seqs
[__i
].second
; // Very end.
165 _GLIBCXX_PARALLEL_ASSERT(__m
!= 0);
166 _GLIBCXX_PARALLEL_ASSERT(__N
!= 0);
167 _GLIBCXX_PARALLEL_ASSERT(__rank
>= 0);
168 _GLIBCXX_PARALLEL_ASSERT(__rank
< __N
);
170 _DifferenceType
* __ns
= new _DifferenceType
[__m
];
171 _DifferenceType
* __a
= new _DifferenceType
[__m
];
172 _DifferenceType
* __b
= new _DifferenceType
[__m
];
175 __ns
[0] = std::distance(__begin_seqs
[0].first
, __begin_seqs
[0].second
);
177 for (int __i
= 0; __i
< __m
; __i
++)
179 __ns
[__i
] = std::distance(__begin_seqs
[__i
].first
,
180 __begin_seqs
[__i
].second
);
181 __nmax
= std::max(__nmax
, __ns
[__i
]);
184 __r
= __log2(__nmax
) + 1;
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;
190 // From now on, including padding.
193 for (int __i
= 0; __i
< __m
; __i
++)
201 // 0 <= __a[__i] <= __ns[__i], 0 <= __b[__i] <= __l
203 #define __S(__i) (__begin_seqs[__i].first)
205 // Initial partition.
206 std::vector
<std::pair
<_ValueType
, int> > __sample
;
208 for (int __i
= 0; __i
< __m
; __i
++)
209 if (__n
< __ns
[__i
]) //__sequence long enough
210 __sample
.push_back(std::make_pair(__S(__i
)[__n
], __i
));
211 __gnu_sequential::sort(__sample
.begin(), __sample
.end(), __lcomp
);
213 for (int __i
= 0; __i
< __m
; __i
++) //conceptual infinity
214 if (__n
>= __ns
[__i
]) //__sequence too short, conceptual infinity
216 std::make_pair(__S(__i
)[0] /*__dummy element*/, __i
));
218 _DifferenceType localrank
= __rank
* __m
/ __N
;
222 __j
< localrank
&& ((__n
+ 1) <= __ns
[__sample
[__j
].second
]); ++__j
)
223 __a
[__sample
[__j
].second
] += __n
+ 1;
224 for (; __j
< __m
; __j
++)
225 __b
[__sample
[__j
].second
] -= __n
+ 1;
227 // Further refinement.
232 int __lmax_seq
= -1; // to avoid warning
233 const _ValueType
* __lmax
= NULL
; // impossible to avoid the warning?
234 for (int __i
= 0; __i
< __m
; __i
++)
240 __lmax
= &(__S(__i
)[__a
[__i
] - 1]);
245 // Max, favor rear sequences.
246 if (!__comp(__S(__i
)[__a
[__i
] - 1], *__lmax
))
248 __lmax
= &(__S(__i
)[__a
[__i
] - 1]);
256 for (__i
= 0; __i
< __m
; __i
++)
258 _DifferenceType __middle
= (__b
[__i
] + __a
[__i
]) / 2;
259 if (__lmax
&& __middle
< __ns
[__i
] &&
260 __lcomp(std::make_pair(__S(__i
)[__middle
], __i
),
261 std::make_pair(*__lmax
, __lmax_seq
)))
262 __a
[__i
] = std::min(__a
[__i
] + __n
+ 1, __ns
[__i
]);
267 _DifferenceType __leftsize
= 0, __total
= 0;
268 for (int __i
= 0; __i
< __m
; __i
++)
270 __leftsize
+= __a
[__i
] / (__n
+ 1);
271 __total
+= __l
/ (__n
+ 1);
274 _DifferenceType __skew
= static_cast<_DifferenceType
>
275 (static_cast<uint64
>(__total
) * __rank
/ __N
- __leftsize
);
279 // Move to the left, find smallest.
280 std::priority_queue
<std::pair
<_ValueType
, int>,
281 std::vector
<std::pair
<_ValueType
, int> >,
282 _LexicographicReverse
<_ValueType
, int, _Compare
> >
285 for (int __i
= 0; __i
< __m
; __i
++)
286 if (__b
[__i
] < __ns
[__i
])
287 __pq
.push(std::make_pair(__S(__i
)[__b
[__i
]], __i
));
289 for (; __skew
!= 0 && !__pq
.empty(); --__skew
)
291 int source
= __pq
.top().second
;
294 __a
[source
] = std::min(__a
[source
] + __n
+ 1, __ns
[source
]);
295 __b
[source
] += __n
+ 1;
297 if (__b
[source
] < __ns
[source
])
299 std::make_pair(__S(source
)[__b
[source
]], source
));
304 // Move to the right, find greatest.
305 std::priority_queue
<std::pair
<_ValueType
, int>,
306 std::vector
<std::pair
<_ValueType
, int> >,
307 _Lexicographic
<_ValueType
, int, _Compare
> > __pq(__lcomp
);
309 for (int __i
= 0; __i
< __m
; __i
++)
311 __pq
.push(std::make_pair(__S(__i
)[__a
[__i
] - 1], __i
));
313 for (; __skew
!= 0; ++__skew
)
315 int source
= __pq
.top().second
;
318 __a
[source
] -= __n
+ 1;
319 __b
[source
] -= __n
+ 1;
323 std::make_pair(__S(source
)[__a
[source
] - 1], source
));
329 // __a[__i] == __b[__i] in most cases, except when __a[__i] has been
330 // clamped because of having reached the boundary
332 // Now return the result, calculate the offset.
334 // Compare the keys on both edges of the border.
336 // Maximum of left edge, minimum of right edge.
337 _ValueType
* __maxleft
= NULL
;
338 _ValueType
* __minright
= NULL
;
339 for (int __i
= 0; __i
< __m
; __i
++)
344 __maxleft
= &(__S(__i
)[__a
[__i
] - 1]);
347 // Max, favor rear sequences.
348 if (!__comp(__S(__i
)[__a
[__i
] - 1], *__maxleft
))
349 __maxleft
= &(__S(__i
)[__a
[__i
] - 1]);
352 if (__b
[__i
] < __ns
[__i
])
355 __minright
= &(__S(__i
)[__b
[__i
]]);
358 // Min, favor fore sequences.
359 if (__comp(__S(__i
)[__b
[__i
]], *__minright
))
360 __minright
= &(__S(__i
)[__b
[__i
]]);
366 for (int __i
= 0; __i
< __m
; __i
++)
367 __begin_offsets
[__i
] = __S(__i
) + __a
[__i
];
376 * @brief Selects the element at __a certain global __rank from several
379 * The sequences are passed via __a __sequence of random-access
380 * iterator pairs, none of the sequences may be empty.
381 * @param __begin_seqs Begin of the sequence of iterator pairs.
382 * @param __end_seqs End of the sequence of iterator pairs.
383 * @param __rank The global __rank to partition at.
384 * @param __offset The rank of the selected element in the global
385 * subsequence of elements equal to the selected element. If the
386 * selected element is unique, this number is 0.
387 * @param __comp The ordering functor, defaults to std::less.
389 template<typename _Tp
, typename _RanSeqs
, typename _RankType
,
392 multiseq_selection(_RanSeqs __begin_seqs
, _RanSeqs __end_seqs
,
394 _RankType
& __offset
, _Compare __comp
= std::less
<_Tp
>())
396 _GLIBCXX_CALL(__end_seqs
- __begin_seqs
)
398 typedef typename
std::iterator_traits
<_RanSeqs
>::value_type::first_type
400 typedef typename
std::iterator_traits
<_It
>::difference_type
403 _Lexicographic
<_Tp
, int, _Compare
> __lcomp(__comp
);
404 _LexicographicReverse
<_Tp
, int, _Compare
> __lrcomp(__comp
);
406 // Number of sequences, number of elements in total (possibly
407 // including padding).
408 _DifferenceType __m
= std::distance(__begin_seqs
, __end_seqs
);
409 _DifferenceType __N
= 0;
410 _DifferenceType __nmax
, __n
, __r
;
412 for (int __i
= 0; __i
< __m
; __i
++)
413 __N
+= std::distance(__begin_seqs
[__i
].first
,
414 __begin_seqs
[__i
].second
);
416 if (__m
== 0 || __N
== 0 || __rank
< 0 || __rank
>= __N
)
418 // result undefined if there is no data or __rank is outside bounds
419 throw std::exception();
423 _DifferenceType
* __ns
= new _DifferenceType
[__m
];
424 _DifferenceType
* __a
= new _DifferenceType
[__m
];
425 _DifferenceType
* __b
= new _DifferenceType
[__m
];
428 __ns
[0] = std::distance(__begin_seqs
[0].first
, __begin_seqs
[0].second
);
430 for (int __i
= 0; __i
< __m
; ++__i
)
432 __ns
[__i
] = std::distance(__begin_seqs
[__i
].first
,
433 __begin_seqs
[__i
].second
);
434 __nmax
= std::max(__nmax
, __ns
[__i
]);
437 __r
= __log2(__nmax
) + 1;
439 // Pad all lists to this length, at least as long as any ns[__i],
440 // equality iff __nmax = 2^__k - 1
443 // From now on, including padding.
446 for (int __i
= 0; __i
< __m
; ++__i
)
454 // 0 <= __a[__i] <= __ns[__i], 0 <= __b[__i] <= __l
456 #define __S(__i) (__begin_seqs[__i].first)
458 // Initial partition.
459 std::vector
<std::pair
<_Tp
, int> > __sample
;
461 for (int __i
= 0; __i
< __m
; __i
++)
463 __sample
.push_back(std::make_pair(__S(__i
)[__n
], __i
));
464 __gnu_sequential::sort(__sample
.begin(), __sample
.end(),
465 __lcomp
, sequential_tag());
467 // Conceptual infinity.
468 for (int __i
= 0; __i
< __m
; __i
++)
469 if (__n
>= __ns
[__i
])
471 std::make_pair(__S(__i
)[0] /*__dummy element*/, __i
));
473 _DifferenceType localrank
= __rank
* __m
/ __N
;
477 __j
< localrank
&& ((__n
+ 1) <= __ns
[__sample
[__j
].second
]); ++__j
)
478 __a
[__sample
[__j
].second
] += __n
+ 1;
479 for (; __j
< __m
; ++__j
)
480 __b
[__sample
[__j
].second
] -= __n
+ 1;
482 // Further refinement.
487 const _Tp
* __lmax
= NULL
;
488 for (int __i
= 0; __i
< __m
; ++__i
)
493 __lmax
= &(__S(__i
)[__a
[__i
] - 1]);
496 if (__comp(*__lmax
, __S(__i
)[__a
[__i
] - 1])) //max
497 __lmax
= &(__S(__i
)[__a
[__i
] - 1]);
503 for (__i
= 0; __i
< __m
; __i
++)
505 _DifferenceType __middle
= (__b
[__i
] + __a
[__i
]) / 2;
506 if (__lmax
&& __middle
< __ns
[__i
]
507 && __comp(__S(__i
)[__middle
], *__lmax
))
508 __a
[__i
] = std::min(__a
[__i
] + __n
+ 1, __ns
[__i
]);
513 _DifferenceType __leftsize
= 0, __total
= 0;
514 for (int __i
= 0; __i
< __m
; ++__i
)
516 __leftsize
+= __a
[__i
] / (__n
+ 1);
517 __total
+= __l
/ (__n
+ 1);
520 _DifferenceType __skew
= ((unsigned long long)__total
* __rank
/ __N
525 // Move to the left, find smallest.
526 std::priority_queue
<std::pair
<_Tp
, int>,
527 std::vector
<std::pair
<_Tp
, int> >,
528 _LexicographicReverse
<_Tp
, int, _Compare
> > __pq(__lrcomp
);
530 for (int __i
= 0; __i
< __m
; ++__i
)
531 if (__b
[__i
] < __ns
[__i
])
532 __pq
.push(std::make_pair(__S(__i
)[__b
[__i
]], __i
));
534 for (; __skew
!= 0 && !__pq
.empty(); --__skew
)
536 int source
= __pq
.top().second
;
539 __a
[source
] = std::min(__a
[source
] + __n
+ 1, __ns
[source
]);
540 __b
[source
] += __n
+ 1;
542 if (__b
[source
] < __ns
[source
])
544 std::make_pair(__S(source
)[__b
[source
]], source
));
549 // Move to the right, find greatest.
550 std::priority_queue
<std::pair
<_Tp
, int>,
551 std::vector
<std::pair
<_Tp
, int> >,
552 _Lexicographic
<_Tp
, int, _Compare
> > __pq(__lcomp
);
554 for (int __i
= 0; __i
< __m
; ++__i
)
556 __pq
.push(std::make_pair(__S(__i
)[__a
[__i
] - 1], __i
));
558 for (; __skew
!= 0; ++__skew
)
560 int source
= __pq
.top().second
;
563 __a
[source
] -= __n
+ 1;
564 __b
[source
] -= __n
+ 1;
568 std::make_pair(__S(source
)[__a
[source
] - 1], source
));
574 // __a[__i] == __b[__i] in most cases, except when __a[__i] has been
575 // clamped because of having reached the boundary
577 // Now return the result, calculate the offset.
579 // Compare the keys on both edges of the border.
581 // Maximum of left edge, minimum of right edge.
582 bool __maxleftset
= false, __minrightset
= false;
584 // Impossible to avoid the warning?
585 _Tp __maxleft
, __minright
;
586 for (int __i
= 0; __i
< __m
; ++__i
)
592 __maxleft
= __S(__i
)[__a
[__i
] - 1];
598 if (__comp(__maxleft
, __S(__i
)[__a
[__i
] - 1]))
599 __maxleft
= __S(__i
)[__a
[__i
] - 1];
602 if (__b
[__i
] < __ns
[__i
])
606 __minright
= __S(__i
)[__b
[__i
]];
607 __minrightset
= true;
612 if (__comp(__S(__i
)[__b
[__i
]], __minright
))
613 __minright
= __S(__i
)[__b
[__i
]];
618 // Minright is the splitter, in any case.
620 if (!__maxleftset
|| __comp(__minright
, __maxleft
))
622 // Good luck, everything is split unambiguously.
627 // We have to calculate an offset.
630 for (int __i
= 0; __i
< __m
; ++__i
)
633 = std::lower_bound(__S(__i
), __S(__i
) + __ns
[__i
],
636 __offset
+= __a
[__i
] - lb
;
650 #endif /* _GLIBCXX_PARALLEL_MULTISEQ_SELECTION_H */