multiseq_selection.h

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

Generated on Sat Oct 25 05:09:05 2008 for libstdc++ by  doxygen 1.5.6