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 Comparator
>
56 : public std::binary_function
<std::pair
<T1
, T2
>, std::pair
<T1
, T2
>, bool>
62 lexicographic(Comparator
& _comp
) : comp(_comp
) { }
65 operator()(const std::pair
<T1
, T2
>& p1
,
66 const std::pair
<T1
, T2
>& p2
) const
68 if (comp(p1
.first
, p2
.first
))
71 if (comp(p2
.first
, p1
.first
))
75 return p1
.second
< p2
.second
;
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>
87 lexicographic_reverse(Comparator
& _comp
) : comp(_comp
) { }
90 operator()(const std::pair
<T1
, T2
>& p1
,
91 const std::pair
<T1
, T2
>& p2
) const
93 if (comp(p2
.first
, p1
.first
))
96 if (comp(p1
.first
, p2
.first
))
100 return p2
.second
< p1
.second
;
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>.
120 template<typename RanSeqs
, typename RankType
, typename RankIterator
,
123 multiseq_partition(RanSeqs begin_seqs
, RanSeqs end_seqs
,
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>
131 _GLIBCXX_CALL(end_seqs
- begin_seqs
)
133 typedef typename
std::iterator_traits
<RanSeqs
>::value_type::first_type
135 typedef typename
std::iterator_traits
<It
>::difference_type
137 typedef typename
std::iterator_traits
<It
>::value_type value_type
;
139 lexicographic
<value_type
, int, Comparator
> lcomp(comp
);
140 lexicographic_reverse
<value_type
, int, Comparator
> lrcomp(comp
);
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,
147 for (int i
= 0; i
< m
; i
++)
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);
156 for (int i
= 0; i
< m
; i
++)
157 begin_offsets
[i
] = begin_seqs
[i
].second
; // Very end.
162 _GLIBCXX_PARALLEL_ASSERT(m
!= 0);
163 _GLIBCXX_PARALLEL_ASSERT(N
!= 0);
164 _GLIBCXX_PARALLEL_ASSERT(rank
>= 0);
165 _GLIBCXX_PARALLEL_ASSERT(rank
< N
);
167 difference_type
* ns
= new difference_type
[m
];
168 difference_type
* a
= new difference_type
[m
];
169 difference_type
* b
= new difference_type
[m
];
172 ns
[0] = std::distance(begin_seqs
[0].first
, begin_seqs
[0].second
);
174 for (int i
= 0; i
< m
; i
++)
176 ns
[i
] = std::distance(begin_seqs
[i
].first
, begin_seqs
[i
].second
);
177 nmax
= std::max(nmax
, ns
[i
]);
180 r
= __log2(nmax
) + 1;
182 // Pad all lists to this length, at least as long as any ns[i],
183 // equality iff nmax = 2^k - 1.
186 // From now on, including padding.
189 for (int i
= 0; i
< m
; i
++)
197 // 0 <= a[i] <= ns[i], 0 <= b[i] <= l
199 #define S(i) (begin_seqs[i].first)
201 // Initial partition.
202 std::vector
<std::pair
<value_type
, int> > sample
;
204 for (int i
= 0; i
< m
; i
++)
205 if (n
< ns
[i
]) //sequence long enough
206 sample
.push_back(std::make_pair(S(i
)[n
], i
));
207 __gnu_sequential::sort(sample
.begin(), sample
.end(), lcomp
);
209 for (int i
= 0; i
< m
; i
++) //conceptual infinity
210 if (n
>= ns
[i
]) //sequence too short, conceptual infinity
211 sample
.push_back(std::make_pair(S(i
)[0] /*dummy element*/, i
));
213 difference_type localrank
= rank
* m
/ N
;
216 for (j
= 0; j
< localrank
&& ((n
+ 1) <= ns
[sample
[j
].second
]); ++j
)
217 a
[sample
[j
].second
] += n
+ 1;
219 b
[sample
[j
].second
] -= n
+ 1;
221 // Further refinement.
226 int lmax_seq
= -1; // to avoid warning
227 const value_type
* lmax
= NULL
; // impossible to avoid the warning?
228 for (int i
= 0; i
< m
; i
++)
234 lmax
= &(S(i
)[a
[i
] - 1]);
239 // Max, favor rear sequences.
240 if (!comp(S(i
)[a
[i
] - 1], *lmax
))
242 lmax
= &(S(i
)[a
[i
] - 1]);
250 for (i
= 0; i
< m
; i
++)
252 difference_type middle
= (b
[i
] + a
[i
]) / 2;
253 if (lmax
&& middle
< ns
[i
] &&
254 lcomp(std::make_pair(S(i
)[middle
], i
),
255 std::make_pair(*lmax
, lmax_seq
)))
256 a
[i
] = std::min(a
[i
] + n
+ 1, ns
[i
]);
261 difference_type leftsize
= 0, total
= 0;
262 for (int i
= 0; i
< m
; i
++)
264 leftsize
+= a
[i
] / (n
+ 1);
265 total
+= l
/ (n
+ 1);
268 difference_type skew
= static_cast<difference_type
>
269 (static_cast<uint64
>(total
) * rank
/ N
- leftsize
);
273 // Move to the left, find smallest.
274 std::priority_queue
<std::pair
<value_type
, int>,
275 std::vector
<std::pair
<value_type
, int> >,
276 lexicographic_reverse
<value_type
, int, Comparator
> >
279 for (int i
= 0; i
< m
; i
++)
281 pq
.push(std::make_pair(S(i
)[b
[i
]], i
));
283 for (; skew
!= 0 && !pq
.empty(); --skew
)
285 int source
= pq
.top().second
;
288 a
[source
] = std::min(a
[source
] + n
+ 1, ns
[source
]);
291 if (b
[source
] < ns
[source
])
292 pq
.push(std::make_pair(S(source
)[b
[source
]], source
));
297 // Move to the right, find greatest.
298 std::priority_queue
<std::pair
<value_type
, int>,
299 std::vector
<std::pair
<value_type
, int> >,
300 lexicographic
<value_type
, int, Comparator
> > pq(lcomp
);
302 for (int i
= 0; i
< m
; i
++)
304 pq
.push(std::make_pair(S(i
)[a
[i
] - 1], i
));
306 for (; skew
!= 0; ++skew
)
308 int source
= pq
.top().second
;
315 pq
.push(std::make_pair(S(source
)[a
[source
] - 1], source
));
321 // a[i] == b[i] in most cases, except when a[i] has been clamped
322 // because of having reached the boundary
324 // Now return the result, calculate the offset.
326 // Compare the keys on both edges of the border.
328 // Maximum of left edge, minimum of right edge.
329 value_type
* maxleft
= NULL
;
330 value_type
* minright
= NULL
;
331 for (int i
= 0; i
< m
; i
++)
336 maxleft
= &(S(i
)[a
[i
] - 1]);
339 // Max, favor rear sequences.
340 if (!comp(S(i
)[a
[i
] - 1], *maxleft
))
341 maxleft
= &(S(i
)[a
[i
] - 1]);
347 minright
= &(S(i
)[b
[i
]]);
350 // Min, favor fore sequences.
351 if (comp(S(i
)[b
[i
]], *minright
))
352 minright
= &(S(i
)[b
[i
]]);
358 for (int i
= 0; i
< m
; i
++)
359 begin_offsets
[i
] = S(i
) + a
[i
];
368 * @brief Selects the element at a certain global rank from several
371 * The sequences are passed via a sequence of random-access
372 * iterator pairs, none of the sequences may be empty.
373 * @param begin_seqs Begin of the sequence of iterator pairs.
374 * @param end_seqs End of the sequence of iterator pairs.
375 * @param rank The global rank to partition at.
376 * @param offset The rank of the selected element in the global
377 * subsequence of elements equal to the selected element. If the
378 * selected element is unique, this number is 0.
379 * @param comp The ordering functor, defaults to std::less.
381 template<typename T
, typename RanSeqs
, typename RankType
,
384 multiseq_selection(RanSeqs begin_seqs
, RanSeqs end_seqs
, RankType rank
,
385 RankType
& offset
, Comparator comp
= std::less
<T
>())
387 _GLIBCXX_CALL(end_seqs
- begin_seqs
)
389 typedef typename
std::iterator_traits
<RanSeqs
>::value_type::first_type
391 typedef typename
std::iterator_traits
<It
>::difference_type
394 lexicographic
<T
, int, Comparator
> lcomp(comp
);
395 lexicographic_reverse
<T
, int, Comparator
> lrcomp(comp
);
397 // Number of sequences, number of elements in total (possibly
398 // including padding).
399 difference_type m
= std::distance(begin_seqs
, end_seqs
);
400 difference_type N
= 0;
401 difference_type nmax
, n
, r
;
403 for (int i
= 0; i
< m
; i
++)
404 N
+= std::distance(begin_seqs
[i
].first
, begin_seqs
[i
].second
);
406 if (m
== 0 || N
== 0 || rank
< 0 || rank
>= N
)
408 // Result undefined when there is no data or rank is outside bounds.
409 throw std::exception();
413 difference_type
* ns
= new difference_type
[m
];
414 difference_type
* a
= new difference_type
[m
];
415 difference_type
* b
= new difference_type
[m
];
418 ns
[0] = std::distance(begin_seqs
[0].first
, begin_seqs
[0].second
);
420 for (int i
= 0; i
< m
; ++i
)
422 ns
[i
] = std::distance(begin_seqs
[i
].first
, begin_seqs
[i
].second
);
423 nmax
= std::max(nmax
, ns
[i
]);
426 r
= __log2(nmax
) + 1;
428 // Pad all lists to this length, at least as long as any ns[i],
429 // equality iff nmax = 2^k - 1
432 // From now on, including padding.
435 for (int i
= 0; i
< m
; ++i
)
443 // 0 <= a[i] <= ns[i], 0 <= b[i] <= l
445 #define S(i) (begin_seqs[i].first)
447 // Initial partition.
448 std::vector
<std::pair
<T
, int> > sample
;
450 for (int i
= 0; i
< m
; i
++)
452 sample
.push_back(std::make_pair(S(i
)[n
], i
));
453 __gnu_sequential::sort(sample
.begin(), sample
.end(),
454 lcomp
, sequential_tag());
456 // Conceptual infinity.
457 for (int i
= 0; i
< m
; i
++)
459 sample
.push_back(std::make_pair(S(i
)[0] /*dummy element*/, i
));
461 difference_type localrank
= rank
* m
/ N
;
464 for (j
= 0; j
< localrank
&& ((n
+ 1) <= ns
[sample
[j
].second
]); ++j
)
465 a
[sample
[j
].second
] += n
+ 1;
467 b
[sample
[j
].second
] -= n
+ 1;
469 // Further refinement.
474 const T
* lmax
= NULL
;
475 for (int i
= 0; i
< m
; ++i
)
480 lmax
= &(S(i
)[a
[i
] - 1]);
483 if (comp(*lmax
, S(i
)[a
[i
] - 1])) //max
484 lmax
= &(S(i
)[a
[i
] - 1]);
490 for (i
= 0; i
< m
; i
++)
492 difference_type middle
= (b
[i
] + a
[i
]) / 2;
493 if (lmax
&& middle
< ns
[i
] && comp(S(i
)[middle
], *lmax
))
494 a
[i
] = std::min(a
[i
] + n
+ 1, ns
[i
]);
499 difference_type leftsize
= 0, total
= 0;
500 for (int i
= 0; i
< m
; ++i
)
502 leftsize
+= a
[i
] / (n
+ 1);
503 total
+= l
/ (n
+ 1);
506 difference_type skew
= ((unsigned long long)total
* rank
/ N
511 // Move to the left, find smallest.
512 std::priority_queue
<std::pair
<T
, int>,
513 std::vector
<std::pair
<T
, int> >,
514 lexicographic_reverse
<T
, int, Comparator
> > pq(lrcomp
);
516 for (int i
= 0; i
< m
; ++i
)
518 pq
.push(std::make_pair(S(i
)[b
[i
]], i
));
520 for (; skew
!= 0 && !pq
.empty(); --skew
)
522 int source
= pq
.top().second
;
525 a
[source
] = std::min(a
[source
] + n
+ 1, ns
[source
]);
528 if (b
[source
] < ns
[source
])
529 pq
.push(std::make_pair(S(source
)[b
[source
]], source
));
534 // Move to the right, find greatest.
535 std::priority_queue
<std::pair
<T
, int>,
536 std::vector
<std::pair
<T
, int> >,
537 lexicographic
<T
, int, Comparator
> > pq(lcomp
);
539 for (int i
= 0; i
< m
; ++i
)
541 pq
.push(std::make_pair(S(i
)[a
[i
] - 1], i
));
543 for (; skew
!= 0; ++skew
)
545 int source
= pq
.top().second
;
552 pq
.push(std::make_pair(S(source
)[a
[source
] - 1], source
));
558 // a[i] == b[i] in most cases, except when a[i] has been clamped
559 // because of having reached the boundary
561 // Now return the result, calculate the offset.
563 // Compare the keys on both edges of the border.
565 // Maximum of left edge, minimum of right edge.
566 bool maxleftset
= false, minrightset
= false;
568 // Impossible to avoid the warning?
570 for (int i
= 0; i
< m
; ++i
)
576 maxleft
= S(i
)[a
[i
] - 1];
582 if (comp(maxleft
, S(i
)[a
[i
] - 1]))
583 maxleft
= S(i
)[a
[i
] - 1];
590 minright
= S(i
)[b
[i
]];
596 if (comp(S(i
)[b
[i
]], minright
))
597 minright
= S(i
)[b
[i
]];
602 // Minright is the splitter, in any case.
604 if (!maxleftset
|| comp(minright
, maxleft
))
606 // Good luck, everything is split unambiguously.
611 // We have to calculate an offset.
614 for (int i
= 0; i
< m
; ++i
)
616 difference_type lb
= std::lower_bound(S(i
), S(i
) + ns
[i
],
633 #endif /* _GLIBCXX_PARALLEL_MULTISEQ_SELECTION_H */