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