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