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