Branch data 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
|
