Line data    Source code

       1              : //
       2              : // Class Partitioner
       3              : //   Partition a domain into subdomains.
       4              : //
       5              : 
       6              : #include <algorithm>
       7              : #include <numeric>
       8              : #include <vector>
       9              : 
      10              : namespace ippl {
      11              :     namespace detail {
      12              : 
      13              :         template <unsigned Dim>
      14              :         template <typename view_type>
      15            0 :         void Partitioner<Dim>::split(const NDIndex<Dim>& domain, view_type& view,
      16              :                                      const std::array<bool, Dim>& isParallel, int nSplits) const {
      17              :             using NDIndex_t = NDIndex<Dim>;
      18              : 
      19              :             // Recursively split the domain until we have generated all the domains.
      20            0 :             std::vector<NDIndex_t> domains_c(nSplits);
      21            0 :             NDIndex_t leftDomain;
      22              : 
      23              :             // Start with the whole domain.
      24            0 :             domains_c[0] = domain;
      25              :             int v;
      26            0 :             unsigned int d = 0;
      27              : 
      28              :             int v1, v2, rm, vtot, vl, vr;
      29              :             double a, lmax, len;
      30              : 
      31            0 :             for (v = nSplits, rm = 0; v > 1; v /= 2) {
      32            0 :                 rm += (v % 2);
      33              :             }
      34              : 
      35            0 :             if (rm == 0) {
      36              :                 // nSplits is a power of 2
      37              : 
      38            0 :                 std::vector<NDIndex_t> copy_c(nSplits);
      39              : 
      40            0 :                 for (v = 1; v < nSplits; v *= 2) {
      41              :                     // Go to the next parallel dimension.
      42            0 :                     while (!isParallel[d])
      43            0 :                         if (++d == Dim)
      44            0 :                             d = 0;
      45              : 
      46              :                     // Split all the current nSplits.
      47              :                     int i, j;
      48            0 :                     for (i = 0, j = 0; i < v; ++i, j += 2) {
      49              :                         // Split to the left and to the right, saving both.
      50            0 :                         domains_c[i].split(copy_c[j], copy_c[j + 1], d);
      51              :                     }
      52              :                     // Copy back.
      53            0 :                     std::copy(copy_c.begin(), copy_c.begin() + v * 2, domains_c.begin());
      54              : 
      55              :                     // On to the next dimension.
      56            0 :                     if (++d == Dim)
      57            0 :                         d = 0;
      58              :                 }
      59              : 
      60            0 :             } else {
      61            0 :                 vtot = 1;  // count the number of nSplits to make sure that it worked
      62              :                            // nSplits is not a power of 2 so we need to do some fancy splitting
      63              :                            // sorry... this would be much cleaner with recursion
      64              :                 /*
      65              :                     The way this works is to recursively split on the longest dimension.
      66              :                     Suppose you request 11 nSplits.  It will split the longest dimension
      67              :                     in the ratio 5:6 and put the new domains in node 0 and node 5.  Then
      68              :                     it splits the longest dimension of the 0 domain and puts the results
      69              :                     in node 0 and node 2 and then splits the longest dimension of node 5
      70              :                     and puts the results in node 5 and node 8. etc.
      71              :                     The logic is kind of bizarre, but it works.
      72              :                 */
      73            0 :                 for (v = 1; v < 2 * nSplits; ++v) {
      74              :                     // kind of reverse the bits of v
      75            0 :                     for (v2 = v, v1 = 1; v2 > 1; v2 /= 2) {
      76            0 :                         v1 = 2 * v1 + (v2 % 2);
      77              :                     }
      78            0 :                     vl = 0;
      79            0 :                     vr = nSplits;
      80              : 
      81            0 :                     while (v1 > 1) {
      82            0 :                         if ((v1 % 2) == 1) {
      83            0 :                             vl = vl + (vr - vl) / 2;
      84              :                         } else {
      85            0 :                             vr = vl + (vr - vl) / 2;
      86              :                         }
      87            0 :                         v1 /= 2;
      88              :                     }
      89              : 
      90            0 :                     v2 = vl + (vr - vl) / 2;
      91              : 
      92            0 :                     if (v2 > vl) {
      93            0 :                         a = v2 - vl;
      94            0 :                         a /= vr - vl;
      95            0 :                         vr         = v2;
      96            0 :                         leftDomain = domains_c[vl];
      97            0 :                         lmax       = 0;
      98            0 :                         d          = std::numeric_limits<unsigned int>::max();
      99            0 :                         for (unsigned int dd = 0; dd < Dim; ++dd) {
     100            0 :                             if (isParallel[dd]) {
     101            0 :                                 if ((len = leftDomain[dd].length()) > lmax) {
     102            0 :                                     lmax = len;
     103            0 :                                     d    = dd;
     104              :                                 }
     105              :                             }
     106              :                         }
     107            0 :                         NDIndex_t temp;
     108            0 :                         domains_c[vl].split(temp, domains_c[vr], d, a);
     109            0 :                         domains_c[vl] = temp;
     110            0 :                         ++vtot;
     111              :                     }
     112              :                 }
     113              : 
     114            0 :                 v = vtot;
     115              :             }
     116              : 
     117              :             // Make sure v is the same number of nSplits at this stage.
     118            0 :             PAssert_EQ(v, nSplits);
     119              : 
     120            0 :             for (size_t i = 0; i < domains_c.size(); ++i) {
     121            0 :                 view(i) = domains_c[i];
     122              :             }
     123            0 :         }
     124              :     }  // namespace detail
     125              : }  // namespace ippl