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