3 // Copyright (C) 2007 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 2, 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 // You should have received a copy of the GNU General Public License
17 // along with this library; see the file COPYING. If not, write to
18 // the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
19 // MA 02111-1307, USA.
21 // As a special exception, you may use this file as part of a free
22 // software library without restriction. Specifically, if other files
23 // instantiate templates or use macros or inline functions from this
24 // file, or you compile this file and link it with other files to
25 // produce an executable, this file does not by itself cause the
26 // resulting executable to be covered by the GNU General Public
27 // License. This exception does not however invalidate any other
28 // reasons why the executable file might be covered by the GNU General
31 /** @file parallel/multiseq_selection.h
32 * @brief Functions to find elements of a certain global rank in
33 * multiple sorted sequences. Also serves for splitting such
35 * This file is a GNU parallel extension to the Standard C++ Library.
38 // Written by Johannes Singler.
40 #ifndef _GLIBCXX_PARALLEL_MULTISEQ_SELECTION_H
41 #define _GLIBCXX_PARALLEL_MULTISEQ_SELECTION_H 1
46 #include <bits/stl_algo.h>
48 #include <parallel/sort.h>
50 namespace __gnu_parallel
52 /** @brief Compare a pair of types lexicographically, ascending. */
53 template<typename T1
, typename T2
, typename Comparator
>
54 class lexicographic
: public std::binary_function
<std::pair
<T1
, T2
>, std::pair
<T1
, T2
>, bool>
60 lexicographic(Comparator
& _comp
) : comp(_comp
) { }
64 operator()(const std::pair
<T1
, T2
>& p1
, const std::pair
<T1
, T2
>& p2
) const
66 if (comp(p1
.first
, p2
.first
))
69 if (comp(p2
.first
, p1
.first
))
73 return p1
.second
< p2
.second
;
77 /** @brief Compare a pair of types lexicographically, descending. */
78 template<typename T1
, typename T2
, typename Comparator
>
79 class lexicographic_reverse
: public std::binary_function
<T1
, T2
, bool>
85 lexicographic_reverse(Comparator
& _comp
) : comp(_comp
) { }
88 operator()(const std::pair
<T1
, T2
>& p1
, const std::pair
<T1
, T2
>& p2
) const
90 if (comp(p2
.first
, p1
.first
))
93 if (comp(p1
.first
, p2
.first
))
97 return p2
.second
< p1
.second
;
102 * @brief Splits several sorted sequences at a certain global rank,
103 * resulting in a splitting point for each sequence.
104 * The sequences are passed via a sequence of random-access
105 * iterator pairs, none of the sequences may be empty. If there
106 * are several equal elements across the split, the ones on the
107 * left side will be chosen from sequences with smaller number.
108 * @param begin_seqs Begin of the sequence of iterator pairs.
109 * @param end_seqs End of the sequence of iterator pairs.
110 * @param rank The global rank to partition at.
111 * @param begin_offsets A random-access sequence begin where the
112 * result will be stored in. Each element of the sequence is an
113 * iterator that points to the first element on the greater part of
114 * the respective sequence.
115 * @param comp The ordering functor, defaults to std::less<T>.
117 template<typename RanSeqs
, typename RankType
, typename RankIterator
, typename Comparator
>
119 multiseq_partition(RanSeqs begin_seqs
, RanSeqs end_seqs
, RankType rank
,
120 RankIterator begin_offsets
,
121 Comparator comp
= std::less
<
122 typename
std::iterator_traits
<typename
std::iterator_traits
<RanSeqs
>::value_type::first_type
>::value_type
>()) // std::less<T>
124 _GLIBCXX_CALL(end_seqs
- begin_seqs
)
126 typedef typename
std::iterator_traits
<RanSeqs
>::value_type::first_type It
;
127 typedef typename
std::iterator_traits
<It
>::difference_type difference_type
;
128 typedef typename
std::iterator_traits
<It
>::value_type T
;
130 lexicographic
<T
, int, Comparator
> lcomp(comp
);
131 lexicographic_reverse
<T
, int, Comparator
> lrcomp(comp
);
133 // Number of sequences, number of elements in total (possibly
134 // including padding).
135 difference_type m
= std::distance(begin_seqs
, end_seqs
), N
= 0, nmax
, n
, r
;
137 for (int i
= 0; i
< m
; i
++)
138 N
+= std::distance(begin_seqs
[i
].first
, begin_seqs
[i
].second
);
142 for (int i
= 0; i
< m
; i
++)
143 begin_offsets
[i
] = begin_seqs
[i
].second
; // Very end.
147 _GLIBCXX_PARALLEL_ASSERT(m
!= 0 && N
!= 0 && rank
>= 0 && rank
< N
);
149 difference_type
* ns
= new difference_type
[m
];
150 difference_type
* a
= new difference_type
[m
];
151 difference_type
* b
= new difference_type
[m
];
154 ns
[0] = std::distance(begin_seqs
[0].first
, begin_seqs
[0].second
);
156 for (int i
= 0; i
< m
; i
++)
158 ns
[i
] = std::distance(begin_seqs
[i
].first
, begin_seqs
[i
].second
);
159 nmax
= std::max(nmax
, ns
[i
]);
164 // Pad all lists to this length, at least as long as any ns[i],
165 // equality iff nmax = 2^k - 1.
168 // From now on, including padding.
171 for (int i
= 0; i
< m
; i
++)
179 // 0 <= a[i] <= ns[i], 0 <= b[i] <= l
181 #define S(i) (begin_seqs[i].first)
183 // Initial partition.
184 std::vector
<std::pair
<T
, int> > sample
;
186 for (int i
= 0; i
< m
; i
++)
187 if (n
< ns
[i
]) //sequence long enough
188 sample
.push_back(std::make_pair(S(i
)[n
], i
));
189 __gnu_sequential::sort(sample
.begin(), sample
.end(), lcomp
);
191 for (int i
= 0; i
< m
; i
++) //conceptual infinity
192 if (n
>= ns
[i
]) //sequence too short, conceptual infinity
193 sample
.push_back(std::make_pair(S(i
)[0] /*dummy element*/, i
));
195 difference_type localrank
= rank
* m
/ N
;
198 for (j
= 0; j
< localrank
&& ((n
+ 1) <= ns
[sample
[j
].second
]); j
++)
199 a
[sample
[j
].second
] += n
+ 1;
201 b
[sample
[j
].second
] -= n
+ 1;
203 // Further refinement.
208 int lmax_seq
= -1; // to avoid warning
209 const T
* lmax
= NULL
; // impossible to avoid the warning?
210 for (int i
= 0; i
< m
; i
++)
216 lmax
= &(S(i
)[a
[i
] - 1]);
221 // Max, favor rear sequences.
222 if (!comp(S(i
)[a
[i
] - 1], *lmax
))
224 lmax
= &(S(i
)[a
[i
] - 1]);
232 for (i
= 0; i
< m
; i
++)
234 difference_type middle
= (b
[i
] + a
[i
]) / 2;
235 if (lmax
&& middle
< ns
[i
] &&
236 lcomp(std::make_pair(S(i
)[middle
], i
), std::make_pair(*lmax
, lmax_seq
)))
237 a
[i
] = std::min(a
[i
] + n
+ 1, ns
[i
]);
242 difference_type leftsize
= 0, total
= 0;
243 for (int i
= 0; i
< m
; i
++)
245 leftsize
+= a
[i
] / (n
+ 1);
246 total
+= l
/ (n
+ 1);
249 difference_type skew
= static_cast<difference_type
>(static_cast<uint64
>(total
) * rank
/ N
- leftsize
);
253 // Move to the left, find smallest.
254 std::priority_queue
<std::pair
<T
, int>, std::vector
<std::pair
<T
, int> >, lexicographic_reverse
<T
, int, Comparator
> > pq(lrcomp
);
256 for (int i
= 0; i
< m
; i
++)
258 pq
.push(std::make_pair(S(i
)[b
[i
]], i
));
260 for (; skew
!= 0 && !pq
.empty(); skew
--)
262 int source
= pq
.top().second
;
265 a
[source
] = std::min(a
[source
] + n
+ 1, ns
[source
]);
268 if (b
[source
] < ns
[source
])
269 pq
.push(std::make_pair(S(source
)[b
[source
]], source
));
274 // Move to the right, find greatest.
275 std::priority_queue
<std::pair
<T
, int>, std::vector
<std::pair
<T
, int> >, lexicographic
<T
, int, Comparator
> > pq(lcomp
);
277 for (int i
= 0; i
< m
; i
++)
279 pq
.push(std::make_pair(S(i
)[a
[i
] - 1], i
));
281 for (; skew
!= 0; skew
++)
283 int source
= pq
.top().second
;
290 pq
.push(std::make_pair(S(source
)[a
[source
] - 1], source
));
296 // a[i] == b[i] in most cases, except when a[i] has been clamped
297 // because of having reached the boundary
299 // Now return the result, calculate the offset.
301 // Compare the keys on both edges of the border.
303 // Maximum of left edge, minimum of right edge.
304 bool maxleftset
= false, minrightset
= false;
305 T maxleft
, minright
; // Impossible to avoid the warning?
306 for (int i
= 0; i
< m
; i
++)
312 maxleft
= S(i
)[a
[i
] - 1];
317 // Max, favor rear sequences.
318 if (!comp(S(i
)[a
[i
] - 1], maxleft
))
319 maxleft
= S(i
)[a
[i
] - 1];
326 minright
= S(i
)[b
[i
]];
331 // Min, favor fore sequences.
332 if (comp(S(i
)[b
[i
]], minright
))
333 minright
= S(i
)[b
[i
]];
339 for (int i
= 0; i
< m
; i
++)
340 begin_offsets
[i
] = S(i
) + a
[i
];
350 * @brief Selects the element at a certain global rank from several
353 * The sequences are passed via a sequence of random-access
354 * iterator pairs, none of the sequences may be empty.
355 * @param begin_seqs Begin of the sequence of iterator pairs.
356 * @param end_seqs End of the sequence of iterator pairs.
357 * @param rank The global rank to partition at.
358 * @param offset The rank of the selected element in the global
359 * subsequence of elements equal to the selected element. If the
360 * selected element is unique, this number is 0.
361 * @param comp The ordering functor, defaults to std::less.
363 template<typename T
, typename RanSeqs
, typename RankType
, typename Comparator
>
365 multiseq_selection(RanSeqs begin_seqs
, RanSeqs end_seqs
, RankType rank
,
366 RankType
& offset
, Comparator comp
= std::less
<T
>())
368 _GLIBCXX_CALL(end_seqs
- begin_seqs
)
370 typedef typename
std::iterator_traits
<RanSeqs
>::value_type::first_type It
;
371 typedef typename
std::iterator_traits
<It
>::difference_type difference_type
;
373 lexicographic
<T
, int, Comparator
> lcomp(comp
);
374 lexicographic_reverse
<T
, int, Comparator
> lrcomp(comp
);
376 // Number of sequences, number of elements in total (possibly
377 // including padding).
378 difference_type m
= std::distance(begin_seqs
, end_seqs
);
379 difference_type N
= 0;
380 difference_type nmax
, n
, r
;
382 for (int i
= 0; i
< m
; i
++)
383 N
+= std::distance(begin_seqs
[i
].first
, begin_seqs
[i
].second
);
385 if (m
== 0 || N
== 0 || rank
< 0 || rank
>= N
)
387 // Result undefined when there is no data or rank is outside bounds.
388 throw std::exception();
392 difference_type
* ns
= new difference_type
[m
];
393 difference_type
* a
= new difference_type
[m
];
394 difference_type
* b
= new difference_type
[m
];
397 ns
[0] = std::distance(begin_seqs
[0].first
, begin_seqs
[0].second
);
399 for (int i
= 0; i
< m
; i
++)
401 ns
[i
] = std::distance(begin_seqs
[i
].first
, begin_seqs
[i
].second
);
402 nmax
= std::max(nmax
, ns
[i
]);
407 // Pad all lists to this length, at least as long as any ns[i],
408 // equality iff nmax = 2^k - 1
411 // From now on, including padding.
414 for (int i
= 0; i
< m
; i
++)
422 // 0 <= a[i] <= ns[i], 0 <= b[i] <= l
424 #define S(i) (begin_seqs[i].first)
426 // Initial partition.
427 std::vector
<std::pair
<T
, int> > sample
;
429 for (int i
= 0; i
< m
; i
++)
431 sample
.push_back(std::make_pair(S(i
)[n
], i
));
432 __gnu_sequential::sort(sample
.begin(), sample
.end(), lcomp
, sequential_tag());
434 // Conceptual infinity.
435 for (int i
= 0; i
< m
; i
++)
437 sample
.push_back(std::make_pair(S(i
)[0] /*dummy element*/, i
));
439 difference_type localrank
= rank
* m
/ N
;
442 for (j
= 0; j
< localrank
&& ((n
+ 1) <= ns
[sample
[j
].second
]); j
++)
443 a
[sample
[j
].second
] += n
+ 1;
445 b
[sample
[j
].second
] -= n
+ 1;
447 // Further refinement.
452 const T
* lmax
= NULL
;
453 for (int i
= 0; i
< m
; i
++)
459 lmax
= &(S(i
)[a
[i
] - 1]);
463 if (comp(*lmax
, S(i
)[a
[i
] - 1])) //max
464 lmax
= &(S(i
)[a
[i
] - 1]);
470 for (i
= 0; i
< m
; i
++)
472 difference_type middle
= (b
[i
] + a
[i
]) / 2;
473 if (lmax
&& middle
< ns
[i
] && comp(S(i
)[middle
], *lmax
))
474 a
[i
] = std::min(a
[i
] + n
+ 1, ns
[i
]);
479 difference_type leftsize
= 0, total
= 0;
480 for (int i
= 0; i
< m
; i
++)
482 leftsize
+= a
[i
] / (n
+ 1);
483 total
+= l
/ (n
+ 1);
486 difference_type skew
= (unsigned long long)total
* rank
/ N
- leftsize
;
490 // Move to the left, find smallest.
491 std::priority_queue
<std::pair
<T
, int>, std::vector
<std::pair
<T
, int> >, lexicographic_reverse
<T
, int, Comparator
> > pq(lrcomp
);
493 for (int i
= 0; i
< m
; i
++)
495 pq
.push(std::make_pair(S(i
)[b
[i
]], i
));
497 for (; skew
!= 0 && !pq
.empty(); skew
--)
499 int source
= pq
.top().second
;
502 a
[source
] = std::min(a
[source
] + n
+ 1, ns
[source
]);
505 if (b
[source
] < ns
[source
])
506 pq
.push(std::make_pair(S(source
)[b
[source
]], source
));
511 // Move to the right, find greatest.
512 std::priority_queue
<std::pair
<T
, int>, std::vector
<std::pair
<T
, int> >, lexicographic
<T
, int, Comparator
> > pq(lcomp
);
514 for (int i
= 0; i
< m
; i
++)
516 pq
.push(std::make_pair(S(i
)[a
[i
] - 1], i
));
518 for (; skew
!= 0; skew
++)
520 int source
= pq
.top().second
;
527 pq
.push(std::make_pair(S(source
)[a
[source
] - 1], source
));
533 // a[i] == b[i] in most cases, except when a[i] has been clamped
534 // because of having reached the boundary
536 // Now return the result, calculate the offset.
538 // Compare the keys on both edges of the border.
540 // Maximum of left edge, minimum of right edge.
541 bool maxleftset
= false, minrightset
= false;
543 // Impossible to avoid the warning?
545 for (int i
= 0; i
< m
; i
++)
551 maxleft
= S(i
)[a
[i
] - 1];
557 if (comp(maxleft
, S(i
)[a
[i
] - 1]))
558 maxleft
= S(i
)[a
[i
] - 1];
565 minright
= S(i
)[b
[i
]];
571 if (comp(S(i
)[b
[i
]], minright
))
572 minright
= S(i
)[b
[i
]];
577 // Minright is the splitter, in any case.
579 if (!maxleftset
|| comp(minright
, maxleft
))
581 // Good luck, everything is split unambigiously.
586 // We have to calculate an offset.
589 for (int i
= 0; i
< m
; i
++)
591 difference_type lb
= std::lower_bound(S(i
), S(i
) + ns
[i
], minright
,