xref: /dflybsd-src/contrib/gcc-4.7/libstdc++-v3/include/parallel/multiseq_selection.h (revision 04febcfb30580676d3e95f58a16c5137ee478b32)
1*e4b17023SJohn Marino // -*- C++ -*-
2*e4b17023SJohn Marino 
3*e4b17023SJohn Marino // Copyright (C) 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
4*e4b17023SJohn Marino //
5*e4b17023SJohn Marino // This file is part of the GNU ISO C++ Library.  This library is free
6*e4b17023SJohn Marino // software; you can redistribute it and/or modify it under the terms
7*e4b17023SJohn Marino // of the GNU General Public License as published by the Free Software
8*e4b17023SJohn Marino // Foundation; either version 3, or (at your option) any later
9*e4b17023SJohn Marino // version.
10*e4b17023SJohn Marino 
11*e4b17023SJohn Marino // This library is distributed in the hope that it will be useful, but
12*e4b17023SJohn Marino // WITHOUT ANY WARRANTY; without even the implied warranty of
13*e4b17023SJohn Marino // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14*e4b17023SJohn Marino // General Public License for more details.
15*e4b17023SJohn Marino 
16*e4b17023SJohn Marino // Under Section 7 of GPL version 3, you are granted additional
17*e4b17023SJohn Marino // permissions described in the GCC Runtime Library Exception, version
18*e4b17023SJohn Marino // 3.1, as published by the Free Software Foundation.
19*e4b17023SJohn Marino 
20*e4b17023SJohn Marino // You should have received a copy of the GNU General Public License and
21*e4b17023SJohn Marino // a copy of the GCC Runtime Library Exception along with this program;
22*e4b17023SJohn Marino // see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
23*e4b17023SJohn Marino // <http://www.gnu.org/licenses/>.
24*e4b17023SJohn Marino 
25*e4b17023SJohn Marino /** @file parallel/multiseq_selection.h
26*e4b17023SJohn Marino  *  @brief Functions to find elements of a certain global __rank in
27*e4b17023SJohn Marino  *  multiple sorted sequences.  Also serves for splitting such
28*e4b17023SJohn Marino  *  sequence sets.
29*e4b17023SJohn Marino  *
30*e4b17023SJohn Marino  *  The algorithm description can be found in
31*e4b17023SJohn Marino  *
32*e4b17023SJohn Marino  *  P. J. Varman, S. D. Scheufler, B. R. Iyer, and G. R. Ricard.
33*e4b17023SJohn Marino  *  Merging Multiple Lists on Hierarchical-Memory Multiprocessors.
34*e4b17023SJohn Marino  *  Journal of Parallel and Distributed Computing, 12(2):171–177, 1991.
35*e4b17023SJohn Marino  *
36*e4b17023SJohn Marino  *  This file is a GNU parallel extension to the Standard C++ Library.
37*e4b17023SJohn Marino  */
38*e4b17023SJohn Marino 
39*e4b17023SJohn Marino // Written by Johannes Singler.
40*e4b17023SJohn Marino 
41*e4b17023SJohn Marino #ifndef _GLIBCXX_PARALLEL_MULTISEQ_SELECTION_H
42*e4b17023SJohn Marino #define _GLIBCXX_PARALLEL_MULTISEQ_SELECTION_H 1
43*e4b17023SJohn Marino 
44*e4b17023SJohn Marino #include <vector>
45*e4b17023SJohn Marino #include <queue>
46*e4b17023SJohn Marino 
47*e4b17023SJohn Marino #include <bits/stl_algo.h>
48*e4b17023SJohn Marino 
49*e4b17023SJohn Marino namespace __gnu_parallel
50*e4b17023SJohn Marino {
51*e4b17023SJohn Marino   /** @brief Compare __a pair of types lexicographically, ascending. */
52*e4b17023SJohn Marino   template<typename _T1, typename _T2, typename _Compare>
53*e4b17023SJohn Marino     class _Lexicographic
54*e4b17023SJohn Marino     : public std::binary_function<std::pair<_T1, _T2>,
55*e4b17023SJohn Marino 				  std::pair<_T1, _T2>, bool>
56*e4b17023SJohn Marino     {
57*e4b17023SJohn Marino     private:
58*e4b17023SJohn Marino       _Compare& _M_comp;
59*e4b17023SJohn Marino 
60*e4b17023SJohn Marino     public:
_Lexicographic(_Compare & __comp)61*e4b17023SJohn Marino       _Lexicographic(_Compare& __comp) : _M_comp(__comp) { }
62*e4b17023SJohn Marino 
63*e4b17023SJohn Marino       bool
operator()64*e4b17023SJohn Marino       operator()(const std::pair<_T1, _T2>& __p1,
65*e4b17023SJohn Marino                  const std::pair<_T1, _T2>& __p2) const
66*e4b17023SJohn Marino       {
67*e4b17023SJohn Marino         if (_M_comp(__p1.first, __p2.first))
68*e4b17023SJohn Marino           return true;
69*e4b17023SJohn Marino 
70*e4b17023SJohn Marino         if (_M_comp(__p2.first, __p1.first))
71*e4b17023SJohn Marino           return false;
72*e4b17023SJohn Marino 
73*e4b17023SJohn Marino         // Firsts are equal.
74*e4b17023SJohn Marino         return __p1.second < __p2.second;
75*e4b17023SJohn Marino       }
76*e4b17023SJohn Marino     };
77*e4b17023SJohn Marino 
78*e4b17023SJohn Marino   /** @brief Compare __a pair of types lexicographically, descending. */
79*e4b17023SJohn Marino   template<typename _T1, typename _T2, typename _Compare>
80*e4b17023SJohn Marino     class _LexicographicReverse : public std::binary_function<_T1, _T2, bool>
81*e4b17023SJohn Marino     {
82*e4b17023SJohn Marino     private:
83*e4b17023SJohn Marino       _Compare& _M_comp;
84*e4b17023SJohn Marino 
85*e4b17023SJohn Marino     public:
_LexicographicReverse(_Compare & __comp)86*e4b17023SJohn Marino       _LexicographicReverse(_Compare& __comp) : _M_comp(__comp) { }
87*e4b17023SJohn Marino 
88*e4b17023SJohn Marino       bool
operator()89*e4b17023SJohn Marino       operator()(const std::pair<_T1, _T2>& __p1,
90*e4b17023SJohn Marino                  const std::pair<_T1, _T2>& __p2) const
91*e4b17023SJohn Marino       {
92*e4b17023SJohn Marino         if (_M_comp(__p2.first, __p1.first))
93*e4b17023SJohn Marino           return true;
94*e4b17023SJohn Marino 
95*e4b17023SJohn Marino         if (_M_comp(__p1.first, __p2.first))
96*e4b17023SJohn Marino           return false;
97*e4b17023SJohn Marino 
98*e4b17023SJohn Marino         // Firsts are equal.
99*e4b17023SJohn Marino         return __p2.second < __p1.second;
100*e4b17023SJohn Marino       }
101*e4b17023SJohn Marino     };
102*e4b17023SJohn Marino 
103*e4b17023SJohn Marino   /**
104*e4b17023SJohn Marino    *  @brief Splits several sorted sequences at a certain global __rank,
105*e4b17023SJohn Marino    *  resulting in a splitting point for each sequence.
106*e4b17023SJohn Marino    *  The sequences are passed via a sequence of random-access
107*e4b17023SJohn Marino    *  iterator pairs, none of the sequences may be empty.  If there
108*e4b17023SJohn Marino    *  are several equal elements across the split, the ones on the
109*e4b17023SJohn Marino    *  __left side will be chosen from sequences with smaller number.
110*e4b17023SJohn Marino    *  @param __begin_seqs Begin of the sequence of iterator pairs.
111*e4b17023SJohn Marino    *  @param __end_seqs End of the sequence of iterator pairs.
112*e4b17023SJohn Marino    *  @param __rank The global rank to partition at.
113*e4b17023SJohn Marino    *  @param __begin_offsets A random-access __sequence __begin where the
114*e4b17023SJohn Marino    *  __result will be stored in. Each element of the sequence is an
115*e4b17023SJohn Marino    *  iterator that points to the first element on the greater part of
116*e4b17023SJohn Marino    *  the respective __sequence.
117*e4b17023SJohn Marino    *  @param __comp The ordering functor, defaults to std::less<_Tp>.
118*e4b17023SJohn Marino    */
119*e4b17023SJohn Marino   template<typename _RanSeqs, typename _RankType, typename _RankIterator,
120*e4b17023SJohn Marino             typename _Compare>
121*e4b17023SJohn Marino     void
122*e4b17023SJohn Marino     multiseq_partition(_RanSeqs __begin_seqs, _RanSeqs __end_seqs,
123*e4b17023SJohn Marino                        _RankType __rank,
124*e4b17023SJohn Marino                        _RankIterator __begin_offsets,
125*e4b17023SJohn Marino                        _Compare __comp = std::less<
126*e4b17023SJohn Marino                        typename std::iterator_traits<typename
127*e4b17023SJohn Marino                        std::iterator_traits<_RanSeqs>::value_type::
128*e4b17023SJohn Marino                        first_type>::value_type>()) // std::less<_Tp>
129*e4b17023SJohn Marino     {
130*e4b17023SJohn Marino       _GLIBCXX_CALL(__end_seqs - __begin_seqs)
131*e4b17023SJohn Marino 
132*e4b17023SJohn Marino       typedef typename std::iterator_traits<_RanSeqs>::value_type::first_type
133*e4b17023SJohn Marino         _It;
134*e4b17023SJohn Marino       typedef typename std::iterator_traits<_RanSeqs>::difference_type
135*e4b17023SJohn Marino         _SeqNumber;
136*e4b17023SJohn Marino       typedef typename std::iterator_traits<_It>::difference_type
137*e4b17023SJohn Marino                _DifferenceType;
138*e4b17023SJohn Marino       typedef typename std::iterator_traits<_It>::value_type _ValueType;
139*e4b17023SJohn Marino 
140*e4b17023SJohn Marino       _Lexicographic<_ValueType, _SeqNumber, _Compare> __lcomp(__comp);
141*e4b17023SJohn Marino       _LexicographicReverse<_ValueType, _SeqNumber, _Compare> __lrcomp(__comp);
142*e4b17023SJohn Marino 
143*e4b17023SJohn Marino       // Number of sequences, number of elements in total (possibly
144*e4b17023SJohn Marino       // including padding).
145*e4b17023SJohn Marino       _DifferenceType __m = std::distance(__begin_seqs, __end_seqs), __nn = 0,
146*e4b17023SJohn Marino                       __nmax, __n, __r;
147*e4b17023SJohn Marino 
148*e4b17023SJohn Marino       for (_SeqNumber __i = 0; __i < __m; __i++)
149*e4b17023SJohn Marino         {
150*e4b17023SJohn Marino           __nn += std::distance(__begin_seqs[__i].first,
151*e4b17023SJohn Marino                                __begin_seqs[__i].second);
152*e4b17023SJohn Marino           _GLIBCXX_PARALLEL_ASSERT(
153*e4b17023SJohn Marino             std::distance(__begin_seqs[__i].first,
154*e4b17023SJohn Marino                           __begin_seqs[__i].second) > 0);
155*e4b17023SJohn Marino         }
156*e4b17023SJohn Marino 
157*e4b17023SJohn Marino       if (__rank == __nn)
158*e4b17023SJohn Marino         {
159*e4b17023SJohn Marino           for (_SeqNumber __i = 0; __i < __m; __i++)
160*e4b17023SJohn Marino             __begin_offsets[__i] = __begin_seqs[__i].second; // Very end.
161*e4b17023SJohn Marino           // Return __m - 1;
162*e4b17023SJohn Marino           return;
163*e4b17023SJohn Marino         }
164*e4b17023SJohn Marino 
165*e4b17023SJohn Marino       _GLIBCXX_PARALLEL_ASSERT(__m != 0);
166*e4b17023SJohn Marino       _GLIBCXX_PARALLEL_ASSERT(__nn != 0);
167*e4b17023SJohn Marino       _GLIBCXX_PARALLEL_ASSERT(__rank >= 0);
168*e4b17023SJohn Marino       _GLIBCXX_PARALLEL_ASSERT(__rank < __nn);
169*e4b17023SJohn Marino 
170*e4b17023SJohn Marino       _DifferenceType* __ns = new _DifferenceType[__m];
171*e4b17023SJohn Marino       _DifferenceType* __a = new _DifferenceType[__m];
172*e4b17023SJohn Marino       _DifferenceType* __b = new _DifferenceType[__m];
173*e4b17023SJohn Marino       _DifferenceType __l;
174*e4b17023SJohn Marino 
175*e4b17023SJohn Marino       __ns[0] = std::distance(__begin_seqs[0].first, __begin_seqs[0].second);
176*e4b17023SJohn Marino       __nmax = __ns[0];
177*e4b17023SJohn Marino       for (_SeqNumber __i = 0; __i < __m; __i++)
178*e4b17023SJohn Marino         {
179*e4b17023SJohn Marino           __ns[__i] = std::distance(__begin_seqs[__i].first,
180*e4b17023SJohn Marino                                     __begin_seqs[__i].second);
181*e4b17023SJohn Marino           __nmax = std::max(__nmax, __ns[__i]);
182*e4b17023SJohn Marino         }
183*e4b17023SJohn Marino 
184*e4b17023SJohn Marino       __r = __rd_log2(__nmax) + 1;
185*e4b17023SJohn Marino 
186*e4b17023SJohn Marino       // Pad all lists to this length, at least as long as any ns[__i],
187*e4b17023SJohn Marino       // equality iff __nmax = 2^__k - 1.
188*e4b17023SJohn Marino       __l = (1ULL << __r) - 1;
189*e4b17023SJohn Marino 
190*e4b17023SJohn Marino       for (_SeqNumber __i = 0; __i < __m; __i++)
191*e4b17023SJohn Marino         {
192*e4b17023SJohn Marino           __a[__i] = 0;
193*e4b17023SJohn Marino           __b[__i] = __l;
194*e4b17023SJohn Marino         }
195*e4b17023SJohn Marino       __n = __l / 2;
196*e4b17023SJohn Marino 
197*e4b17023SJohn Marino       // Invariants:
198*e4b17023SJohn Marino       // 0 <= __a[__i] <= __ns[__i], 0 <= __b[__i] <= __l
199*e4b17023SJohn Marino 
200*e4b17023SJohn Marino #define __S(__i) (__begin_seqs[__i].first)
201*e4b17023SJohn Marino 
202*e4b17023SJohn Marino       // Initial partition.
203*e4b17023SJohn Marino       std::vector<std::pair<_ValueType, _SeqNumber> > __sample;
204*e4b17023SJohn Marino 
205*e4b17023SJohn Marino       for (_SeqNumber __i = 0; __i < __m; __i++)
206*e4b17023SJohn Marino         if (__n < __ns[__i])    //__sequence long enough
207*e4b17023SJohn Marino           __sample.push_back(std::make_pair(__S(__i)[__n], __i));
208*e4b17023SJohn Marino       __gnu_sequential::sort(__sample.begin(), __sample.end(), __lcomp);
209*e4b17023SJohn Marino 
210*e4b17023SJohn Marino       for (_SeqNumber __i = 0; __i < __m; __i++)       //conceptual infinity
211*e4b17023SJohn Marino         if (__n >= __ns[__i])   //__sequence too short, conceptual infinity
212*e4b17023SJohn Marino           __sample.push_back(
213*e4b17023SJohn Marino             std::make_pair(__S(__i)[0] /*__dummy element*/, __i));
214*e4b17023SJohn Marino 
215*e4b17023SJohn Marino       _DifferenceType __localrank = __rank / __l;
216*e4b17023SJohn Marino 
217*e4b17023SJohn Marino       _SeqNumber __j;
218*e4b17023SJohn Marino       for (__j = 0;
219*e4b17023SJohn Marino            __j < __localrank && ((__n + 1) <= __ns[__sample[__j].second]);
220*e4b17023SJohn Marino            ++__j)
221*e4b17023SJohn Marino         __a[__sample[__j].second] += __n + 1;
222*e4b17023SJohn Marino       for (; __j < __m; __j++)
223*e4b17023SJohn Marino         __b[__sample[__j].second] -= __n + 1;
224*e4b17023SJohn Marino 
225*e4b17023SJohn Marino       // Further refinement.
226*e4b17023SJohn Marino       while (__n > 0)
227*e4b17023SJohn Marino         {
228*e4b17023SJohn Marino           __n /= 2;
229*e4b17023SJohn Marino 
230*e4b17023SJohn Marino           _SeqNumber __lmax_seq = -1;  // to avoid warning
231*e4b17023SJohn Marino           const _ValueType* __lmax = 0; // impossible to avoid the warning?
232*e4b17023SJohn Marino           for (_SeqNumber __i = 0; __i < __m; __i++)
233*e4b17023SJohn Marino             {
234*e4b17023SJohn Marino               if (__a[__i] > 0)
235*e4b17023SJohn Marino                 {
236*e4b17023SJohn Marino                   if (!__lmax)
237*e4b17023SJohn Marino                     {
238*e4b17023SJohn Marino                       __lmax = &(__S(__i)[__a[__i] - 1]);
239*e4b17023SJohn Marino                       __lmax_seq = __i;
240*e4b17023SJohn Marino                     }
241*e4b17023SJohn Marino                   else
242*e4b17023SJohn Marino                     {
243*e4b17023SJohn Marino                       // Max, favor rear sequences.
244*e4b17023SJohn Marino                       if (!__comp(__S(__i)[__a[__i] - 1], *__lmax))
245*e4b17023SJohn Marino                         {
246*e4b17023SJohn Marino                           __lmax = &(__S(__i)[__a[__i] - 1]);
247*e4b17023SJohn Marino                           __lmax_seq = __i;
248*e4b17023SJohn Marino                         }
249*e4b17023SJohn Marino                     }
250*e4b17023SJohn Marino                 }
251*e4b17023SJohn Marino             }
252*e4b17023SJohn Marino 
253*e4b17023SJohn Marino           _SeqNumber __i;
254*e4b17023SJohn Marino           for (__i = 0; __i < __m; __i++)
255*e4b17023SJohn Marino             {
256*e4b17023SJohn Marino               _DifferenceType __middle = (__b[__i] + __a[__i]) / 2;
257*e4b17023SJohn Marino               if (__lmax && __middle < __ns[__i] &&
258*e4b17023SJohn Marino                   __lcomp(std::make_pair(__S(__i)[__middle], __i),
259*e4b17023SJohn Marino                         std::make_pair(*__lmax, __lmax_seq)))
260*e4b17023SJohn Marino                 __a[__i] = std::min(__a[__i] + __n + 1, __ns[__i]);
261*e4b17023SJohn Marino               else
262*e4b17023SJohn Marino                 __b[__i] -= __n + 1;
263*e4b17023SJohn Marino             }
264*e4b17023SJohn Marino 
265*e4b17023SJohn Marino           _DifferenceType __leftsize = 0;
266*e4b17023SJohn Marino           for (_SeqNumber __i = 0; __i < __m; __i++)
267*e4b17023SJohn Marino               __leftsize += __a[__i] / (__n + 1);
268*e4b17023SJohn Marino 
269*e4b17023SJohn Marino           _DifferenceType __skew = __rank / (__n + 1) - __leftsize;
270*e4b17023SJohn Marino 
271*e4b17023SJohn Marino           if (__skew > 0)
272*e4b17023SJohn Marino             {
273*e4b17023SJohn Marino               // Move to the left, find smallest.
274*e4b17023SJohn Marino               std::priority_queue<std::pair<_ValueType, _SeqNumber>,
275*e4b17023SJohn Marino                 std::vector<std::pair<_ValueType, _SeqNumber> >,
276*e4b17023SJohn Marino                 _LexicographicReverse<_ValueType, _SeqNumber, _Compare> >
277*e4b17023SJohn Marino                 __pq(__lrcomp);
278*e4b17023SJohn Marino 
279*e4b17023SJohn Marino               for (_SeqNumber __i = 0; __i < __m; __i++)
280*e4b17023SJohn Marino                 if (__b[__i] < __ns[__i])
281*e4b17023SJohn Marino                   __pq.push(std::make_pair(__S(__i)[__b[__i]], __i));
282*e4b17023SJohn Marino 
283*e4b17023SJohn Marino               for (; __skew != 0 && !__pq.empty(); --__skew)
284*e4b17023SJohn Marino                 {
285*e4b17023SJohn Marino                   _SeqNumber __source = __pq.top().second;
286*e4b17023SJohn Marino                   __pq.pop();
287*e4b17023SJohn Marino 
288*e4b17023SJohn Marino                   __a[__source]
289*e4b17023SJohn Marino                       = std::min(__a[__source] + __n + 1, __ns[__source]);
290*e4b17023SJohn Marino                   __b[__source] += __n + 1;
291*e4b17023SJohn Marino 
292*e4b17023SJohn Marino                   if (__b[__source] < __ns[__source])
293*e4b17023SJohn Marino                     __pq.push(
294*e4b17023SJohn Marino                       std::make_pair(__S(__source)[__b[__source]], __source));
295*e4b17023SJohn Marino                 }
296*e4b17023SJohn Marino             }
297*e4b17023SJohn Marino           else if (__skew < 0)
298*e4b17023SJohn Marino             {
299*e4b17023SJohn Marino               // Move to the right, find greatest.
300*e4b17023SJohn Marino               std::priority_queue<std::pair<_ValueType, _SeqNumber>,
301*e4b17023SJohn Marino                 std::vector<std::pair<_ValueType, _SeqNumber> >,
302*e4b17023SJohn Marino                 _Lexicographic<_ValueType, _SeqNumber, _Compare> >
303*e4b17023SJohn Marino                   __pq(__lcomp);
304*e4b17023SJohn Marino 
305*e4b17023SJohn Marino               for (_SeqNumber __i = 0; __i < __m; __i++)
306*e4b17023SJohn Marino                 if (__a[__i] > 0)
307*e4b17023SJohn Marino                   __pq.push(std::make_pair(__S(__i)[__a[__i] - 1], __i));
308*e4b17023SJohn Marino 
309*e4b17023SJohn Marino               for (; __skew != 0; ++__skew)
310*e4b17023SJohn Marino                 {
311*e4b17023SJohn Marino                   _SeqNumber __source = __pq.top().second;
312*e4b17023SJohn Marino                   __pq.pop();
313*e4b17023SJohn Marino 
314*e4b17023SJohn Marino                   __a[__source] -= __n + 1;
315*e4b17023SJohn Marino                   __b[__source] -= __n + 1;
316*e4b17023SJohn Marino 
317*e4b17023SJohn Marino                   if (__a[__source] > 0)
318*e4b17023SJohn Marino                     __pq.push(std::make_pair(
319*e4b17023SJohn Marino                         __S(__source)[__a[__source] - 1], __source));
320*e4b17023SJohn Marino                 }
321*e4b17023SJohn Marino             }
322*e4b17023SJohn Marino         }
323*e4b17023SJohn Marino 
324*e4b17023SJohn Marino       // Postconditions:
325*e4b17023SJohn Marino       // __a[__i] == __b[__i] in most cases, except when __a[__i] has been
326*e4b17023SJohn Marino       // clamped because of having reached the boundary
327*e4b17023SJohn Marino 
328*e4b17023SJohn Marino       // Now return the result, calculate the offset.
329*e4b17023SJohn Marino 
330*e4b17023SJohn Marino       // Compare the keys on both edges of the border.
331*e4b17023SJohn Marino 
332*e4b17023SJohn Marino       // Maximum of left edge, minimum of right edge.
333*e4b17023SJohn Marino       _ValueType* __maxleft = 0;
334*e4b17023SJohn Marino       _ValueType* __minright = 0;
335*e4b17023SJohn Marino       for (_SeqNumber __i = 0; __i < __m; __i++)
336*e4b17023SJohn Marino         {
337*e4b17023SJohn Marino           if (__a[__i] > 0)
338*e4b17023SJohn Marino             {
339*e4b17023SJohn Marino               if (!__maxleft)
340*e4b17023SJohn Marino                 __maxleft = &(__S(__i)[__a[__i] - 1]);
341*e4b17023SJohn Marino               else
342*e4b17023SJohn Marino                 {
343*e4b17023SJohn Marino                   // Max, favor rear sequences.
344*e4b17023SJohn Marino                   if (!__comp(__S(__i)[__a[__i] - 1], *__maxleft))
345*e4b17023SJohn Marino                     __maxleft = &(__S(__i)[__a[__i] - 1]);
346*e4b17023SJohn Marino                 }
347*e4b17023SJohn Marino             }
348*e4b17023SJohn Marino           if (__b[__i] < __ns[__i])
349*e4b17023SJohn Marino             {
350*e4b17023SJohn Marino               if (!__minright)
351*e4b17023SJohn Marino                 __minright = &(__S(__i)[__b[__i]]);
352*e4b17023SJohn Marino               else
353*e4b17023SJohn Marino                 {
354*e4b17023SJohn Marino                   // Min, favor fore sequences.
355*e4b17023SJohn Marino                   if (__comp(__S(__i)[__b[__i]], *__minright))
356*e4b17023SJohn Marino                     __minright = &(__S(__i)[__b[__i]]);
357*e4b17023SJohn Marino                 }
358*e4b17023SJohn Marino             }
359*e4b17023SJohn Marino         }
360*e4b17023SJohn Marino 
361*e4b17023SJohn Marino       _SeqNumber __seq = 0;
362*e4b17023SJohn Marino       for (_SeqNumber __i = 0; __i < __m; __i++)
363*e4b17023SJohn Marino         __begin_offsets[__i] = __S(__i) + __a[__i];
364*e4b17023SJohn Marino 
365*e4b17023SJohn Marino       delete[] __ns;
366*e4b17023SJohn Marino       delete[] __a;
367*e4b17023SJohn Marino       delete[] __b;
368*e4b17023SJohn Marino     }
369*e4b17023SJohn Marino 
370*e4b17023SJohn Marino 
371*e4b17023SJohn Marino   /**
372*e4b17023SJohn Marino    *  @brief Selects the element at a certain global __rank from several
373*e4b17023SJohn Marino    *  sorted sequences.
374*e4b17023SJohn Marino    *
375*e4b17023SJohn Marino    *  The sequences are passed via a sequence of random-access
376*e4b17023SJohn Marino    *  iterator pairs, none of the sequences may be empty.
377*e4b17023SJohn Marino    *  @param __begin_seqs Begin of the sequence of iterator pairs.
378*e4b17023SJohn Marino    *  @param __end_seqs End of the sequence of iterator pairs.
379*e4b17023SJohn Marino    *  @param __rank The global rank to partition at.
380*e4b17023SJohn Marino    *  @param __offset The rank of the selected element in the global
381*e4b17023SJohn Marino    *  subsequence of elements equal to the selected element. If the
382*e4b17023SJohn Marino    *  selected element is unique, this number is 0.
383*e4b17023SJohn Marino    *  @param __comp The ordering functor, defaults to std::less.
384*e4b17023SJohn Marino    */
385*e4b17023SJohn Marino   template<typename _Tp, typename _RanSeqs, typename _RankType,
386*e4b17023SJohn Marino            typename _Compare>
387*e4b17023SJohn Marino     _Tp
388*e4b17023SJohn Marino     multiseq_selection(_RanSeqs __begin_seqs, _RanSeqs __end_seqs,
389*e4b17023SJohn Marino                        _RankType __rank,
390*e4b17023SJohn Marino                        _RankType& __offset, _Compare __comp = std::less<_Tp>())
391*e4b17023SJohn Marino     {
392*e4b17023SJohn Marino       _GLIBCXX_CALL(__end_seqs - __begin_seqs)
393*e4b17023SJohn Marino 
394*e4b17023SJohn Marino       typedef typename std::iterator_traits<_RanSeqs>::value_type::first_type
395*e4b17023SJohn Marino         _It;
396*e4b17023SJohn Marino       typedef typename std::iterator_traits<_RanSeqs>::difference_type
397*e4b17023SJohn Marino         _SeqNumber;
398*e4b17023SJohn Marino       typedef typename std::iterator_traits<_It>::difference_type
399*e4b17023SJohn Marino         _DifferenceType;
400*e4b17023SJohn Marino 
401*e4b17023SJohn Marino       _Lexicographic<_Tp, _SeqNumber, _Compare> __lcomp(__comp);
402*e4b17023SJohn Marino       _LexicographicReverse<_Tp, _SeqNumber, _Compare> __lrcomp(__comp);
403*e4b17023SJohn Marino 
404*e4b17023SJohn Marino       // Number of sequences, number of elements in total (possibly
405*e4b17023SJohn Marino       // including padding).
406*e4b17023SJohn Marino       _DifferenceType __m = std::distance(__begin_seqs, __end_seqs);
407*e4b17023SJohn Marino       _DifferenceType __nn = 0;
408*e4b17023SJohn Marino       _DifferenceType __nmax, __n, __r;
409*e4b17023SJohn Marino 
410*e4b17023SJohn Marino       for (_SeqNumber __i = 0; __i < __m; __i++)
411*e4b17023SJohn Marino         __nn += std::distance(__begin_seqs[__i].first,
412*e4b17023SJohn Marino 			      __begin_seqs[__i].second);
413*e4b17023SJohn Marino 
414*e4b17023SJohn Marino       if (__m == 0 || __nn == 0 || __rank < 0 || __rank >= __nn)
415*e4b17023SJohn Marino         {
416*e4b17023SJohn Marino           // result undefined if there is no data or __rank is outside bounds
417*e4b17023SJohn Marino           throw std::exception();
418*e4b17023SJohn Marino         }
419*e4b17023SJohn Marino 
420*e4b17023SJohn Marino 
421*e4b17023SJohn Marino       _DifferenceType* __ns = new _DifferenceType[__m];
422*e4b17023SJohn Marino       _DifferenceType* __a = new _DifferenceType[__m];
423*e4b17023SJohn Marino       _DifferenceType* __b = new _DifferenceType[__m];
424*e4b17023SJohn Marino       _DifferenceType __l;
425*e4b17023SJohn Marino 
426*e4b17023SJohn Marino       __ns[0] = std::distance(__begin_seqs[0].first, __begin_seqs[0].second);
427*e4b17023SJohn Marino       __nmax = __ns[0];
428*e4b17023SJohn Marino       for (_SeqNumber __i = 0; __i < __m; ++__i)
429*e4b17023SJohn Marino         {
430*e4b17023SJohn Marino           __ns[__i] = std::distance(__begin_seqs[__i].first,
431*e4b17023SJohn Marino                                     __begin_seqs[__i].second);
432*e4b17023SJohn Marino           __nmax = std::max(__nmax, __ns[__i]);
433*e4b17023SJohn Marino         }
434*e4b17023SJohn Marino 
435*e4b17023SJohn Marino       __r = __rd_log2(__nmax) + 1;
436*e4b17023SJohn Marino 
437*e4b17023SJohn Marino       // Pad all lists to this length, at least as long as any ns[__i],
438*e4b17023SJohn Marino       // equality iff __nmax = 2^__k - 1
439*e4b17023SJohn Marino       __l = __round_up_to_pow2(__r) - 1;
440*e4b17023SJohn Marino 
441*e4b17023SJohn Marino       for (_SeqNumber __i = 0; __i < __m; ++__i)
442*e4b17023SJohn Marino         {
443*e4b17023SJohn Marino           __a[__i] = 0;
444*e4b17023SJohn Marino           __b[__i] = __l;
445*e4b17023SJohn Marino         }
446*e4b17023SJohn Marino       __n = __l / 2;
447*e4b17023SJohn Marino 
448*e4b17023SJohn Marino       // Invariants:
449*e4b17023SJohn Marino       // 0 <= __a[__i] <= __ns[__i], 0 <= __b[__i] <= __l
450*e4b17023SJohn Marino 
451*e4b17023SJohn Marino #define __S(__i) (__begin_seqs[__i].first)
452*e4b17023SJohn Marino 
453*e4b17023SJohn Marino       // Initial partition.
454*e4b17023SJohn Marino       std::vector<std::pair<_Tp, _SeqNumber> > __sample;
455*e4b17023SJohn Marino 
456*e4b17023SJohn Marino       for (_SeqNumber __i = 0; __i < __m; __i++)
457*e4b17023SJohn Marino         if (__n < __ns[__i])
458*e4b17023SJohn Marino           __sample.push_back(std::make_pair(__S(__i)[__n], __i));
459*e4b17023SJohn Marino       __gnu_sequential::sort(__sample.begin(), __sample.end(),
460*e4b17023SJohn Marino                              __lcomp, sequential_tag());
461*e4b17023SJohn Marino 
462*e4b17023SJohn Marino       // Conceptual infinity.
463*e4b17023SJohn Marino       for (_SeqNumber __i = 0; __i < __m; __i++)
464*e4b17023SJohn Marino         if (__n >= __ns[__i])
465*e4b17023SJohn Marino           __sample.push_back(
466*e4b17023SJohn Marino             std::make_pair(__S(__i)[0] /*__dummy element*/, __i));
467*e4b17023SJohn Marino 
468*e4b17023SJohn Marino       _DifferenceType __localrank = __rank / __l;
469*e4b17023SJohn Marino 
470*e4b17023SJohn Marino       _SeqNumber __j;
471*e4b17023SJohn Marino       for (__j = 0;
472*e4b17023SJohn Marino            __j < __localrank && ((__n + 1) <= __ns[__sample[__j].second]);
473*e4b17023SJohn Marino            ++__j)
474*e4b17023SJohn Marino         __a[__sample[__j].second] += __n + 1;
475*e4b17023SJohn Marino       for (; __j < __m; ++__j)
476*e4b17023SJohn Marino         __b[__sample[__j].second] -= __n + 1;
477*e4b17023SJohn Marino 
478*e4b17023SJohn Marino       // Further refinement.
479*e4b17023SJohn Marino       while (__n > 0)
480*e4b17023SJohn Marino         {
481*e4b17023SJohn Marino           __n /= 2;
482*e4b17023SJohn Marino 
483*e4b17023SJohn Marino           const _Tp* __lmax = 0;
484*e4b17023SJohn Marino           for (_SeqNumber __i = 0; __i < __m; ++__i)
485*e4b17023SJohn Marino             {
486*e4b17023SJohn Marino               if (__a[__i] > 0)
487*e4b17023SJohn Marino                 {
488*e4b17023SJohn Marino                   if (!__lmax)
489*e4b17023SJohn Marino                     __lmax = &(__S(__i)[__a[__i] - 1]);
490*e4b17023SJohn Marino                   else
491*e4b17023SJohn Marino                     {
492*e4b17023SJohn Marino                       if (__comp(*__lmax, __S(__i)[__a[__i] - 1]))      //max
493*e4b17023SJohn Marino                         __lmax = &(__S(__i)[__a[__i] - 1]);
494*e4b17023SJohn Marino                     }
495*e4b17023SJohn Marino                 }
496*e4b17023SJohn Marino             }
497*e4b17023SJohn Marino 
498*e4b17023SJohn Marino           _SeqNumber __i;
499*e4b17023SJohn Marino           for (__i = 0; __i < __m; __i++)
500*e4b17023SJohn Marino             {
501*e4b17023SJohn Marino               _DifferenceType __middle = (__b[__i] + __a[__i]) / 2;
502*e4b17023SJohn Marino               if (__lmax && __middle < __ns[__i]
503*e4b17023SJohn Marino                   && __comp(__S(__i)[__middle], *__lmax))
504*e4b17023SJohn Marino                 __a[__i] = std::min(__a[__i] + __n + 1, __ns[__i]);
505*e4b17023SJohn Marino               else
506*e4b17023SJohn Marino                 __b[__i] -= __n + 1;
507*e4b17023SJohn Marino             }
508*e4b17023SJohn Marino 
509*e4b17023SJohn Marino           _DifferenceType __leftsize = 0;
510*e4b17023SJohn Marino           for (_SeqNumber __i = 0; __i < __m; ++__i)
511*e4b17023SJohn Marino               __leftsize += __a[__i] / (__n + 1);
512*e4b17023SJohn Marino 
513*e4b17023SJohn Marino           _DifferenceType __skew = __rank / (__n + 1) - __leftsize;
514*e4b17023SJohn Marino 
515*e4b17023SJohn Marino           if (__skew > 0)
516*e4b17023SJohn Marino             {
517*e4b17023SJohn Marino               // Move to the left, find smallest.
518*e4b17023SJohn Marino               std::priority_queue<std::pair<_Tp, _SeqNumber>,
519*e4b17023SJohn Marino                 std::vector<std::pair<_Tp, _SeqNumber> >,
520*e4b17023SJohn Marino                 _LexicographicReverse<_Tp, _SeqNumber, _Compare> >
521*e4b17023SJohn Marino                   __pq(__lrcomp);
522*e4b17023SJohn Marino 
523*e4b17023SJohn Marino               for (_SeqNumber __i = 0; __i < __m; ++__i)
524*e4b17023SJohn Marino                 if (__b[__i] < __ns[__i])
525*e4b17023SJohn Marino                   __pq.push(std::make_pair(__S(__i)[__b[__i]], __i));
526*e4b17023SJohn Marino 
527*e4b17023SJohn Marino               for (; __skew != 0 && !__pq.empty(); --__skew)
528*e4b17023SJohn Marino                 {
529*e4b17023SJohn Marino                   _SeqNumber __source = __pq.top().second;
530*e4b17023SJohn Marino                   __pq.pop();
531*e4b17023SJohn Marino 
532*e4b17023SJohn Marino                   __a[__source]
533*e4b17023SJohn Marino                       = std::min(__a[__source] + __n + 1, __ns[__source]);
534*e4b17023SJohn Marino                   __b[__source] += __n + 1;
535*e4b17023SJohn Marino 
536*e4b17023SJohn Marino                   if (__b[__source] < __ns[__source])
537*e4b17023SJohn Marino                     __pq.push(
538*e4b17023SJohn Marino                       std::make_pair(__S(__source)[__b[__source]], __source));
539*e4b17023SJohn Marino                 }
540*e4b17023SJohn Marino             }
541*e4b17023SJohn Marino           else if (__skew < 0)
542*e4b17023SJohn Marino             {
543*e4b17023SJohn Marino               // Move to the right, find greatest.
544*e4b17023SJohn Marino               std::priority_queue<std::pair<_Tp, _SeqNumber>,
545*e4b17023SJohn Marino                 std::vector<std::pair<_Tp, _SeqNumber> >,
546*e4b17023SJohn Marino                 _Lexicographic<_Tp, _SeqNumber, _Compare> > __pq(__lcomp);
547*e4b17023SJohn Marino 
548*e4b17023SJohn Marino               for (_SeqNumber __i = 0; __i < __m; ++__i)
549*e4b17023SJohn Marino                 if (__a[__i] > 0)
550*e4b17023SJohn Marino                   __pq.push(std::make_pair(__S(__i)[__a[__i] - 1], __i));
551*e4b17023SJohn Marino 
552*e4b17023SJohn Marino               for (; __skew != 0; ++__skew)
553*e4b17023SJohn Marino                 {
554*e4b17023SJohn Marino                   _SeqNumber __source = __pq.top().second;
555*e4b17023SJohn Marino                   __pq.pop();
556*e4b17023SJohn Marino 
557*e4b17023SJohn Marino                   __a[__source] -= __n + 1;
558*e4b17023SJohn Marino                   __b[__source] -= __n + 1;
559*e4b17023SJohn Marino 
560*e4b17023SJohn Marino                   if (__a[__source] > 0)
561*e4b17023SJohn Marino                     __pq.push(std::make_pair(
562*e4b17023SJohn Marino                         __S(__source)[__a[__source] - 1], __source));
563*e4b17023SJohn Marino                 }
564*e4b17023SJohn Marino             }
565*e4b17023SJohn Marino         }
566*e4b17023SJohn Marino 
567*e4b17023SJohn Marino       // Postconditions:
568*e4b17023SJohn Marino       // __a[__i] == __b[__i] in most cases, except when __a[__i] has been
569*e4b17023SJohn Marino       // clamped because of having reached the boundary
570*e4b17023SJohn Marino 
571*e4b17023SJohn Marino       // Now return the result, calculate the offset.
572*e4b17023SJohn Marino 
573*e4b17023SJohn Marino       // Compare the keys on both edges of the border.
574*e4b17023SJohn Marino 
575*e4b17023SJohn Marino       // Maximum of left edge, minimum of right edge.
576*e4b17023SJohn Marino       bool __maxleftset = false, __minrightset = false;
577*e4b17023SJohn Marino 
578*e4b17023SJohn Marino       // Impossible to avoid the warning?
579*e4b17023SJohn Marino       _Tp __maxleft, __minright;
580*e4b17023SJohn Marino       for (_SeqNumber __i = 0; __i < __m; ++__i)
581*e4b17023SJohn Marino         {
582*e4b17023SJohn Marino           if (__a[__i] > 0)
583*e4b17023SJohn Marino             {
584*e4b17023SJohn Marino               if (!__maxleftset)
585*e4b17023SJohn Marino                 {
586*e4b17023SJohn Marino                   __maxleft = __S(__i)[__a[__i] - 1];
587*e4b17023SJohn Marino                   __maxleftset = true;
588*e4b17023SJohn Marino                 }
589*e4b17023SJohn Marino               else
590*e4b17023SJohn Marino                 {
591*e4b17023SJohn Marino                   // Max.
592*e4b17023SJohn Marino                   if (__comp(__maxleft, __S(__i)[__a[__i] - 1]))
593*e4b17023SJohn Marino                     __maxleft = __S(__i)[__a[__i] - 1];
594*e4b17023SJohn Marino                 }
595*e4b17023SJohn Marino             }
596*e4b17023SJohn Marino           if (__b[__i] < __ns[__i])
597*e4b17023SJohn Marino             {
598*e4b17023SJohn Marino               if (!__minrightset)
599*e4b17023SJohn Marino                 {
600*e4b17023SJohn Marino                   __minright = __S(__i)[__b[__i]];
601*e4b17023SJohn Marino                   __minrightset = true;
602*e4b17023SJohn Marino                 }
603*e4b17023SJohn Marino               else
604*e4b17023SJohn Marino                 {
605*e4b17023SJohn Marino                   // Min.
606*e4b17023SJohn Marino                   if (__comp(__S(__i)[__b[__i]], __minright))
607*e4b17023SJohn Marino                     __minright = __S(__i)[__b[__i]];
608*e4b17023SJohn Marino                 }
609*e4b17023SJohn Marino             }
610*e4b17023SJohn Marino       }
611*e4b17023SJohn Marino 
612*e4b17023SJohn Marino       // Minright is the __splitter, in any case.
613*e4b17023SJohn Marino 
614*e4b17023SJohn Marino       if (!__maxleftset || __comp(__minright, __maxleft))
615*e4b17023SJohn Marino         {
616*e4b17023SJohn Marino           // Good luck, everything is split unambiguously.
617*e4b17023SJohn Marino           __offset = 0;
618*e4b17023SJohn Marino         }
619*e4b17023SJohn Marino       else
620*e4b17023SJohn Marino         {
621*e4b17023SJohn Marino           // We have to calculate an offset.
622*e4b17023SJohn Marino           __offset = 0;
623*e4b17023SJohn Marino 
624*e4b17023SJohn Marino           for (_SeqNumber __i = 0; __i < __m; ++__i)
625*e4b17023SJohn Marino             {
626*e4b17023SJohn Marino               _DifferenceType lb
627*e4b17023SJohn Marino                 = std::lower_bound(__S(__i), __S(__i) + __ns[__i],
628*e4b17023SJohn Marino                                    __minright,
629*e4b17023SJohn Marino                                    __comp) - __S(__i);
630*e4b17023SJohn Marino               __offset += __a[__i] - lb;
631*e4b17023SJohn Marino             }
632*e4b17023SJohn Marino         }
633*e4b17023SJohn Marino 
634*e4b17023SJohn Marino       delete[] __ns;
635*e4b17023SJohn Marino       delete[] __a;
636*e4b17023SJohn Marino       delete[] __b;
637*e4b17023SJohn Marino 
638*e4b17023SJohn Marino       return __minright;
639*e4b17023SJohn Marino     }
640*e4b17023SJohn Marino }
641*e4b17023SJohn Marino 
642*e4b17023SJohn Marino #undef __S
643*e4b17023SJohn Marino 
644*e4b17023SJohn Marino #endif /* _GLIBCXX_PARALLEL_MULTISEQ_SELECTION_H */
645