multiseq_selection.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2007, 2008, 2009 Free Software Foundation, Inc.
00004 //
00005 // This file is part of the GNU ISO C++ Library.  This library is free
00006 // software; you can redistribute it and/or modify it under the terms
00007 // of the GNU General Public License as published by the Free Software
00008 // Foundation; either version 3, or (at your option) any later
00009 // version.
00010 
00011 // This library is distributed in the hope that it will be useful, but
00012 // WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014 // General Public License for more details.
00015 
00016 // Under Section 7 of GPL version 3, you are granted additional
00017 // permissions described in the GCC Runtime Library Exception, version
00018 // 3.1, as published by the Free Software Foundation.
00019 
00020 // You should have received a copy of the GNU General Public License and
00021 // a copy of the GCC Runtime Library Exception along with this program;
00022 // see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
00023 // <http://www.gnu.org/licenses/>.
00024 
00025 /** @file parallel/multiseq_selection.h
00026  *  @brief Functions to find elements of a certain global rank in
00027  *  multiple sorted sequences.  Also serves for splitting such
00028  *  sequence sets.
00029  *
00030  *  The algorithm description can be found in 
00031  *
00032  *  P. J. Varman, S. D. Scheufler, B. R. Iyer, and G. R. Ricard.
00033  *  Merging Multiple Lists on Hierarchical-Memory Multiprocessors.
00034  *  Journal of Parallel and Distributed Computing, 12(2):171–177, 1991.
00035  *
00036  *  This file is a GNU parallel extension to the Standard C++ Library.
00037  */
00038 
00039 // Written by Johannes Singler.
00040 
00041 #ifndef _GLIBCXX_PARALLEL_MULTISEQ_SELECTION_H
00042 #define _GLIBCXX_PARALLEL_MULTISEQ_SELECTION_H 1
00043 
00044 #include <vector>
00045 #include <queue>
00046 
00047 #include <bits/stl_algo.h>
00048 
00049 #include <parallel/sort.h>
00050 
00051 namespace __gnu_parallel
00052 {
00053   /** @brief Compare a pair of types lexicographically, ascending. */
00054   template<typename T1, typename T2, typename Comparator>
00055     class lexicographic
00056     : public std::binary_function<std::pair<T1, T2>, std::pair<T1, T2>, bool>
00057     {
00058     private:
00059       Comparator& comp;
00060 
00061     public:
00062       lexicographic(Comparator& _comp) : comp(_comp) { }
00063 
00064       bool
00065       operator()(const std::pair<T1, T2>& p1,
00066          const std::pair<T1, T2>& p2) const
00067       {
00068     if (comp(p1.first, p2.first))
00069       return true;
00070 
00071     if (comp(p2.first, p1.first))
00072       return false;
00073 
00074     // Firsts are equal.
00075     return p1.second < p2.second;
00076       }
00077     };
00078 
00079   /** @brief Compare a pair of types lexicographically, descending. */
00080   template<typename T1, typename T2, typename Comparator>
00081     class lexicographic_reverse : public std::binary_function<T1, T2, bool>
00082     {
00083     private:
00084       Comparator& comp;
00085 
00086     public:
00087       lexicographic_reverse(Comparator& _comp) : comp(_comp) { }
00088 
00089       bool
00090       operator()(const std::pair<T1, T2>& p1,
00091          const std::pair<T1, T2>& p2) const
00092       {
00093     if (comp(p2.first, p1.first))
00094       return true;
00095 
00096     if (comp(p1.first, p2.first))
00097       return false;
00098 
00099     // Firsts are equal.
00100     return p2.second < p1.second;
00101       }
00102     };
00103 
00104   /** 
00105    *  @brief Splits several sorted sequences at a certain global rank,
00106    *  resulting in a splitting point for each sequence.
00107    *  The sequences are passed via a sequence of random-access
00108    *  iterator pairs, none of the sequences may be empty.  If there
00109    *  are several equal elements across the split, the ones on the
00110    *  left side will be chosen from sequences with smaller number.
00111    *  @param begin_seqs Begin of the sequence of iterator pairs.
00112    *  @param end_seqs End of the sequence of iterator pairs.
00113    *  @param rank The global rank to partition at.
00114    *  @param begin_offsets A random-access sequence begin where the
00115    *  result will be stored in. Each element of the sequence is an
00116    *  iterator that points to the first element on the greater part of
00117    *  the respective sequence.
00118    *  @param comp The ordering functor, defaults to std::less<T>. 
00119    */
00120   template<typename RanSeqs, typename RankType, typename RankIterator,
00121             typename Comparator>
00122     void
00123     multiseq_partition(RanSeqs begin_seqs, RanSeqs end_seqs,
00124                        RankType rank,
00125                        RankIterator begin_offsets,
00126                        Comparator comp = std::less<
00127                        typename std::iterator_traits<typename
00128                        std::iterator_traits<RanSeqs>::value_type::
00129                        first_type>::value_type>()) // std::less<T>
00130     {
00131       _GLIBCXX_CALL(end_seqs - begin_seqs)
00132 
00133       typedef typename std::iterator_traits<RanSeqs>::value_type::first_type
00134         It;
00135       typedef typename std::iterator_traits<It>::difference_type
00136            difference_type;
00137       typedef typename std::iterator_traits<It>::value_type value_type;
00138 
00139       lexicographic<value_type, int, Comparator> lcomp(comp);
00140       lexicographic_reverse<value_type, int, Comparator> lrcomp(comp);
00141 
00142       // Number of sequences, number of elements in total (possibly
00143       // including padding).
00144       difference_type m = std::distance(begin_seqs, end_seqs), N = 0,
00145                       nmax, n, r;
00146 
00147       for (int i = 0; i < m; i++)
00148         {
00149           N += std::distance(begin_seqs[i].first, begin_seqs[i].second);
00150           _GLIBCXX_PARALLEL_ASSERT(
00151             std::distance(begin_seqs[i].first, begin_seqs[i].second) > 0);
00152         }
00153 
00154       if (rank == N)
00155         {
00156           for (int i = 0; i < m; i++)
00157             begin_offsets[i] = begin_seqs[i].second; // Very end.
00158           // Return m - 1;
00159           return;
00160         }
00161 
00162       _GLIBCXX_PARALLEL_ASSERT(m != 0);
00163       _GLIBCXX_PARALLEL_ASSERT(N != 0);
00164       _GLIBCXX_PARALLEL_ASSERT(rank >= 0);
00165       _GLIBCXX_PARALLEL_ASSERT(rank < N);
00166 
00167       difference_type* ns = new difference_type[m];
00168       difference_type* a = new difference_type[m];
00169       difference_type* b = new difference_type[m];
00170       difference_type l;
00171 
00172       ns[0] = std::distance(begin_seqs[0].first, begin_seqs[0].second);
00173       nmax = ns[0];
00174       for (int i = 0; i < m; i++)
00175     {
00176       ns[i] = std::distance(begin_seqs[i].first, begin_seqs[i].second);
00177       nmax = std::max(nmax, ns[i]);
00178     }
00179 
00180       r = __log2(nmax) + 1;
00181 
00182       // Pad all lists to this length, at least as long as any ns[i],
00183       // equality iff nmax = 2^k - 1.
00184       l = (1ULL << r) - 1;
00185 
00186       // From now on, including padding.
00187       N = l * m;
00188 
00189       for (int i = 0; i < m; i++)
00190     {
00191       a[i] = 0;
00192       b[i] = l;
00193     }
00194       n = l / 2;
00195 
00196       // Invariants:
00197       // 0 <= a[i] <= ns[i], 0 <= b[i] <= l
00198 
00199 #define S(i) (begin_seqs[i].first)
00200 
00201       // Initial partition.
00202       std::vector<std::pair<value_type, int> > sample;
00203 
00204       for (int i = 0; i < m; i++)
00205     if (n < ns[i])  //sequence long enough
00206       sample.push_back(std::make_pair(S(i)[n], i));
00207       __gnu_sequential::sort(sample.begin(), sample.end(), lcomp);
00208 
00209       for (int i = 0; i < m; i++)   //conceptual infinity
00210     if (n >= ns[i]) //sequence too short, conceptual infinity
00211       sample.push_back(std::make_pair(S(i)[0] /*dummy element*/, i));
00212 
00213       difference_type localrank = rank * m / N ;
00214 
00215       int j;
00216       for (j = 0; j < localrank && ((n + 1) <= ns[sample[j].second]); ++j)
00217     a[sample[j].second] += n + 1;
00218       for (; j < m; j++)
00219     b[sample[j].second] -= n + 1;
00220       
00221       // Further refinement.
00222       while (n > 0)
00223     {
00224       n /= 2;
00225 
00226       int lmax_seq = -1;    // to avoid warning
00227       const value_type* lmax = NULL; // impossible to avoid the warning?
00228       for (int i = 0; i < m; i++)
00229         {
00230           if (a[i] > 0)
00231         {
00232           if (!lmax)
00233             {
00234               lmax = &(S(i)[a[i] - 1]);
00235               lmax_seq = i;
00236             }
00237           else
00238             {
00239               // Max, favor rear sequences.
00240               if (!comp(S(i)[a[i] - 1], *lmax))
00241             {
00242               lmax = &(S(i)[a[i] - 1]);
00243               lmax_seq = i;
00244             }
00245             }
00246         }
00247         }
00248 
00249       int i;
00250       for (i = 0; i < m; i++)
00251         {
00252           difference_type middle = (b[i] + a[i]) / 2;
00253           if (lmax && middle < ns[i] &&
00254           lcomp(std::make_pair(S(i)[middle], i),
00255             std::make_pair(*lmax, lmax_seq)))
00256         a[i] = std::min(a[i] + n + 1, ns[i]);
00257           else
00258         b[i] -= n + 1;
00259         }
00260 
00261       difference_type leftsize = 0, total = 0;
00262       for (int i = 0; i < m; i++)
00263         {
00264           leftsize += a[i] / (n + 1);
00265           total += l / (n + 1);
00266         }
00267       
00268       difference_type skew = static_cast<difference_type>
00269         (static_cast<uint64>(total) * rank / N - leftsize);
00270 
00271       if (skew > 0)
00272         {
00273           // Move to the left, find smallest.
00274           std::priority_queue<std::pair<value_type, int>,
00275         std::vector<std::pair<value_type, int> >,
00276         lexicographic_reverse<value_type, int, Comparator> >
00277         pq(lrcomp);
00278           
00279           for (int i = 0; i < m; i++)
00280         if (b[i] < ns[i])
00281           pq.push(std::make_pair(S(i)[b[i]], i));
00282 
00283           for (; skew != 0 && !pq.empty(); --skew)
00284         {
00285           int source = pq.top().second;
00286           pq.pop();
00287 
00288           a[source] = std::min(a[source] + n + 1, ns[source]);
00289           b[source] += n + 1;
00290 
00291           if (b[source] < ns[source])
00292             pq.push(std::make_pair(S(source)[b[source]], source));
00293         }
00294         }
00295       else if (skew < 0)
00296         {
00297           // Move to the right, find greatest.
00298           std::priority_queue<std::pair<value_type, int>,
00299         std::vector<std::pair<value_type, int> >,
00300         lexicographic<value_type, int, Comparator> > pq(lcomp);
00301 
00302           for (int i = 0; i < m; i++)
00303         if (a[i] > 0)
00304           pq.push(std::make_pair(S(i)[a[i] - 1], i));
00305 
00306           for (; skew != 0; ++skew)
00307         {
00308           int source = pq.top().second;
00309           pq.pop();
00310 
00311           a[source] -= n + 1;
00312           b[source] -= n + 1;
00313 
00314           if (a[source] > 0)
00315             pq.push(std::make_pair(S(source)[a[source] - 1], source));
00316         }
00317         }
00318     }
00319 
00320       // Postconditions:
00321       // a[i] == b[i] in most cases, except when a[i] has been clamped
00322       // because of having reached the boundary
00323 
00324       // Now return the result, calculate the offset.
00325 
00326       // Compare the keys on both edges of the border.
00327 
00328       // Maximum of left edge, minimum of right edge.
00329       value_type* maxleft = NULL;
00330       value_type* minright = NULL;
00331       for (int i = 0; i < m; i++)
00332     {
00333       if (a[i] > 0)
00334         {
00335           if (!maxleft)
00336         maxleft = &(S(i)[a[i] - 1]);
00337           else
00338         {
00339           // Max, favor rear sequences.
00340           if (!comp(S(i)[a[i] - 1], *maxleft))
00341             maxleft = &(S(i)[a[i] - 1]);
00342         }
00343         }
00344       if (b[i] < ns[i])
00345         {
00346           if (!minright)
00347         minright = &(S(i)[b[i]]);
00348           else
00349         {
00350           // Min, favor fore sequences.
00351           if (comp(S(i)[b[i]], *minright))
00352             minright = &(S(i)[b[i]]);
00353         }
00354         }
00355     }
00356 
00357       int seq = 0;
00358       for (int i = 0; i < m; i++)
00359     begin_offsets[i] = S(i) + a[i];
00360 
00361       delete[] ns;
00362       delete[] a;
00363       delete[] b;
00364     }
00365 
00366 
00367   /** 
00368    *  @brief Selects the element at a certain global rank from several
00369    *  sorted sequences.
00370    *
00371    *  The sequences are passed via a sequence of random-access
00372    *  iterator pairs, none of the sequences may be empty.
00373    *  @param begin_seqs Begin of the sequence of iterator pairs.
00374    *  @param end_seqs End of the sequence of iterator pairs.
00375    *  @param rank The global rank to partition at.
00376    *  @param offset The rank of the selected element in the global
00377    *  subsequence of elements equal to the selected element. If the
00378    *  selected element is unique, this number is 0.
00379    *  @param comp The ordering functor, defaults to std::less. 
00380    */
00381   template<typename T, typename RanSeqs, typename RankType,
00382        typename Comparator>
00383     T
00384     multiseq_selection(RanSeqs begin_seqs, RanSeqs end_seqs, RankType rank,
00385                RankType& offset, Comparator comp = std::less<T>())
00386     {
00387       _GLIBCXX_CALL(end_seqs - begin_seqs)
00388 
00389       typedef typename std::iterator_traits<RanSeqs>::value_type::first_type
00390     It;
00391       typedef typename std::iterator_traits<It>::difference_type
00392     difference_type;
00393 
00394       lexicographic<T, int, Comparator> lcomp(comp);
00395       lexicographic_reverse<T, int, Comparator> lrcomp(comp);
00396 
00397       // Number of sequences, number of elements in total (possibly
00398       // including padding).
00399       difference_type m = std::distance(begin_seqs, end_seqs);
00400       difference_type N = 0;
00401       difference_type nmax, n, r;
00402 
00403       for (int i = 0; i < m; i++)
00404     N += std::distance(begin_seqs[i].first, begin_seqs[i].second);
00405 
00406       if (m == 0 || N == 0 || rank < 0 || rank >= N)
00407     {
00408       // Result undefined when there is no data or rank is outside bounds.
00409       throw std::exception();
00410     }
00411 
00412 
00413       difference_type* ns = new difference_type[m];
00414       difference_type* a = new difference_type[m];
00415       difference_type* b = new difference_type[m];
00416       difference_type l;
00417 
00418       ns[0] = std::distance(begin_seqs[0].first, begin_seqs[0].second);
00419       nmax = ns[0];
00420       for (int i = 0; i < m; ++i)
00421     {
00422       ns[i] = std::distance(begin_seqs[i].first, begin_seqs[i].second);
00423       nmax = std::max(nmax, ns[i]);
00424     }
00425 
00426       r = __log2(nmax) + 1;
00427 
00428       // Pad all lists to this length, at least as long as any ns[i],
00429       // equality iff nmax = 2^k - 1
00430       l = pow2(r) - 1;
00431 
00432       // From now on, including padding.
00433       N = l * m;
00434 
00435       for (int i = 0; i < m; ++i)
00436     {
00437       a[i] = 0;
00438       b[i] = l;
00439     }
00440       n = l / 2;
00441 
00442       // Invariants:
00443       // 0 <= a[i] <= ns[i], 0 <= b[i] <= l
00444 
00445 #define S(i) (begin_seqs[i].first)
00446 
00447       // Initial partition.
00448       std::vector<std::pair<T, int> > sample;
00449 
00450       for (int i = 0; i < m; i++)
00451     if (n < ns[i])
00452       sample.push_back(std::make_pair(S(i)[n], i));
00453       __gnu_sequential::sort(sample.begin(), sample.end(),
00454                  lcomp, sequential_tag());
00455 
00456       // Conceptual infinity.
00457       for (int i = 0; i < m; i++)
00458     if (n >= ns[i])
00459       sample.push_back(std::make_pair(S(i)[0] /*dummy element*/, i));
00460 
00461       difference_type localrank = rank * m / N ;
00462 
00463       int j;
00464       for (j = 0; j < localrank && ((n + 1) <= ns[sample[j].second]); ++j)
00465     a[sample[j].second] += n + 1;
00466       for (; j < m; ++j)
00467     b[sample[j].second] -= n + 1;
00468 
00469       // Further refinement.
00470       while (n > 0)
00471     {
00472       n /= 2;
00473 
00474       const T* lmax = NULL;
00475       for (int i = 0; i < m; ++i)
00476         {
00477           if (a[i] > 0)
00478         {
00479           if (!lmax)
00480             lmax = &(S(i)[a[i] - 1]);
00481           else
00482             {
00483               if (comp(*lmax, S(i)[a[i] - 1]))  //max
00484             lmax = &(S(i)[a[i] - 1]);
00485             }
00486         }
00487         }
00488 
00489       int i;
00490       for (i = 0; i < m; i++)
00491         {
00492           difference_type middle = (b[i] + a[i]) / 2;
00493           if (lmax && middle < ns[i] && comp(S(i)[middle], *lmax))
00494         a[i] = std::min(a[i] + n + 1, ns[i]);
00495           else
00496         b[i] -= n + 1;
00497         }
00498 
00499       difference_type leftsize = 0, total = 0;
00500       for (int i = 0; i < m; ++i)
00501         {
00502           leftsize += a[i] / (n + 1);
00503           total += l / (n + 1);
00504         }
00505 
00506       difference_type skew = ((unsigned long long)total * rank / N
00507                   - leftsize);
00508 
00509       if (skew > 0)
00510         {
00511           // Move to the left, find smallest.
00512           std::priority_queue<std::pair<T, int>,
00513         std::vector<std::pair<T, int> >,
00514         lexicographic_reverse<T, int, Comparator> > pq(lrcomp);
00515 
00516           for (int i = 0; i < m; ++i)
00517         if (b[i] < ns[i])
00518           pq.push(std::make_pair(S(i)[b[i]], i));
00519 
00520           for (; skew != 0 && !pq.empty(); --skew)
00521         {
00522           int source = pq.top().second;
00523           pq.pop();
00524           
00525           a[source] = std::min(a[source] + n + 1, ns[source]);
00526           b[source] += n + 1;
00527           
00528           if (b[source] < ns[source])
00529             pq.push(std::make_pair(S(source)[b[source]], source));
00530         }
00531         }
00532       else if (skew < 0)
00533         {
00534           // Move to the right, find greatest.
00535           std::priority_queue<std::pair<T, int>,
00536         std::vector<std::pair<T, int> >,
00537         lexicographic<T, int, Comparator> > pq(lcomp);
00538 
00539           for (int i = 0; i < m; ++i)
00540         if (a[i] > 0)
00541           pq.push(std::make_pair(S(i)[a[i] - 1], i));
00542 
00543           for (; skew != 0; ++skew)
00544         {
00545           int source = pq.top().second;
00546           pq.pop();
00547 
00548           a[source] -= n + 1;
00549           b[source] -= n + 1;
00550 
00551           if (a[source] > 0)
00552             pq.push(std::make_pair(S(source)[a[source] - 1], source));
00553         }
00554         }
00555     }
00556 
00557       // Postconditions:
00558       // a[i] == b[i] in most cases, except when a[i] has been clamped
00559       // because of having reached the boundary
00560 
00561       // Now return the result, calculate the offset.
00562 
00563       // Compare the keys on both edges of the border.
00564 
00565       // Maximum of left edge, minimum of right edge.
00566       bool maxleftset = false, minrightset = false;
00567 
00568       // Impossible to avoid the warning?
00569       T maxleft, minright;
00570       for (int i = 0; i < m; ++i)
00571     {
00572       if (a[i] > 0)
00573         {
00574           if (!maxleftset)
00575         {
00576           maxleft = S(i)[a[i] - 1];
00577           maxleftset = true;
00578         }
00579           else
00580         {
00581           // Max.
00582           if (comp(maxleft, S(i)[a[i] - 1]))
00583             maxleft = S(i)[a[i] - 1];
00584         }
00585         }
00586       if (b[i] < ns[i])
00587         {
00588           if (!minrightset)
00589         {
00590           minright = S(i)[b[i]];
00591           minrightset = true;
00592         }
00593           else
00594         {
00595           // Min.
00596           if (comp(S(i)[b[i]], minright))
00597             minright = S(i)[b[i]];
00598         }
00599         }
00600       }
00601 
00602       // Minright is the splitter, in any case.
00603 
00604       if (!maxleftset || comp(minright, maxleft))
00605     {
00606       // Good luck, everything is split unambiguously.
00607       offset = 0;
00608     }
00609       else
00610     {
00611       // We have to calculate an offset.
00612       offset = 0;
00613 
00614       for (int i = 0; i < m; ++i)
00615         {
00616           difference_type lb = std::lower_bound(S(i), S(i) + ns[i],
00617                             minright,
00618                             comp) - S(i);
00619           offset += a[i] - lb;
00620         }
00621     }
00622 
00623       delete[] ns;
00624       delete[] a;
00625       delete[] b;
00626 
00627       return minright;
00628     }
00629 }
00630 
00631 #undef S
00632 
00633 #endif /* _GLIBCXX_PARALLEL_MULTISEQ_SELECTION_H */

Generated on Tue Apr 21 13:13:28 2009 for libstdc++ by  doxygen 1.5.8