MkFitCore/src/radix_sort.cc

0001 #include "RecoTracker/MkFitCore/interface/radix_sort.h"
0002
0003 #include <array>
0004 #include <numeric>
0005
0006 namespace mkfit {
0007
0008   // --- Driver function
0009
0010   template <typename V, typename R>
0011   void radix_sort<V, R>::sort(const std::vector<V>& values, std::vector<R>& ranks) {
0012     if (values.empty()) {
0013       ranks.clear();
0014       return;
0015     }
0016
0017     rank_t histos[c_NBytes * 256] = {0};
0018
0019     histo_loop(values, histos);
0020     radix_loop(values, histos, ranks);
0021   }
0022
0023   // --- Histogramming
0024
0025   template <typename V, typename R>
0026   void radix_sort<V, R>::histo_loop(const std::vector<V>& values, rank_t* histos) {
0027     // Create histograms (counters). Counters for all passes are created in one run.
0028     const ubyte_t* p = (const ubyte_t*)values.data();
0029     const ubyte_t* pe = p + (values.size() * c_NBytes);
0030     std::array<rank_t*, c_NBytes> ha;
0031     for (rank_t j = 0; j < c_NBytes; ++j)
0032       ha[j] = &histos[j << 8];
0033     while (p != pe) {
0034       for (rank_t j = 0; j < c_NBytes; ++j)
0035         ha[j][*p++]++;
0036     }
0037   }
0038
0039   // --- Radix
0040
0041   template <typename V, typename R>
0042   void radix_sort<V, R>::radix_loop(const std::vector<V>& values, rank_t* histos, std::vector<R>& ranks) {
0043     const rank_t nb = values.size();
0044     rank_t* link[256];
0045     ranks.resize(nb);
0046     std::vector<rank_t> ranks2(nb);
0047     bool ranks_are_invalid = true;
0048     // Radix sort, j is the pass number (0=LSB, 3=MSB)
0049     for (rank_t j = 0; j < c_NBytes; j++) {
0050       // Shortcut to current counters
0051       rank_t* cur_count = &histos[j << 8];
0052
0053       // Get first byte - if that byte's counter equals nb, all values are the same.
0054       const ubyte_t first_entry_val = *(((const ubyte_t*)values.data()) + j);
0055       if (cur_count[first_entry_val] != nb) {
0056         // Create offsets
0057         link[0] = ranks2.data();
0058         for (rank_t i = 1; i < 256; i++)
0059           link[i] = link[i - 1] + cur_count[i - 1];
0060
0061         // Perform Radix Sort
0062         ubyte_t* input_bytes = (ubyte_t*)values.data();
0063         input_bytes += j;
0064         if (ranks_are_invalid) {
0065           for (rank_t i = 0; i < nb; i++)
0066             *link[input_bytes[i << 2]]++ = i;
0067           ranks_are_invalid = false;
0068         } else {
0069           rank_t* indices = ranks.data();
0070           rank_t* indices_end = indices + nb;
0071           while (indices != indices_end) {
0072             rank_t id = *indices++;
0073             *link[input_bytes[id << 2]]++ = id;
0074           }
0075         }
0076
0077         // Swap ranks - valid indices are in ranks after the swap.
0078         ranks.swap(ranks2);
0079       }
0080     }
0081     // If all values are equal, fill ranks with sequential integers.
0082     if (ranks_are_invalid)
0083       std::iota(ranks.begin(), ranks.end(), 0);
0084   }
0085
0086   // Instantiate supported sort types.
0087   template class radix_sort<unsigned int, unsigned int>;
0088   template class radix_sort<unsigned int, unsigned short>;
0089 }  // namespace mkfit