Line data Source code
1 : //
2 : // Class Index
3 : // Define a slice in an array.
4 : //
5 : // This essentially defines a list of evenly spaced numbers.
6 : // Most commonly this list will be increasing (positive stride)
7 : // but it can also have negative stride and be decreasing.
8 : //
9 : // Index() --> A null interval with no elements.
10 : // Index(n) --> make an Index on [0..n-1]
11 : // Index(a,b) --> make an Index on [a..b]
12 : // Index(a,b,s) --> make an Index on [a..b] with stride s
13 : //
14 : // Example1:
15 : // --------
16 : // Index I(10); --> Index on [0..9]
17 : // Index Low(5); --> Index on [0..4]
18 : // Index High(5,9); --> Index on [5..9]
19 : // Index IOdd(1,9,2); --> Index on [1..9] stride 2
20 : // Index IEven(0,9,2); --> Index on [0..9] stride 2
21 : //
22 : // Given an Index I(a,n,s), and an integer j you can do the following:
23 : //
24 : // I+j : a+j+i*s for i in [0..n-1]
25 : // j-I : j-a-i*s
26 : // j*I : j*a + i*j*s
27 : // I/j : a/j + i*s/j
28 : //
29 : // j/I we don't do because it is not a uniform stride, and we don't
30 : // allow strides that are fractions.
31 : //
32 : #include "Utility/PAssert.h"
33 :
34 : namespace ippl {
35 :
36 749448 : KOKKOS_INLINE_FUNCTION Index::Index()
37 749448 : : first_m(0)
38 749448 : , stride_m(0)
39 749448 : , length_m(0) {}
40 :
41 4996 : KOKKOS_INLINE_FUNCTION Index::Index(size_t n)
42 4996 : : first_m(0)
43 4996 : , stride_m(1)
44 4996 : , length_m(n) {}
45 :
46 6912 : KOKKOS_INLINE_FUNCTION Index::Index(int f, int l)
47 6912 : : first_m(f)
48 6912 : , stride_m(1)
49 6912 : , length_m(l - f + 1) {
50 6912 : PAssert_GE(l - f + 1, 0);
51 6912 : }
52 :
53 0 : KOKKOS_INLINE_FUNCTION Index::Index(int f, int l, int s)
54 0 : : first_m(f)
55 0 : , stride_m(s) {
56 0 : PAssert_NE(s, 0);
57 0 : if (f == l) {
58 0 : length_m = 1;
59 0 : } else if ((l > f) ^ (s < 0)) {
60 0 : length_m = (l - f) / s + 1;
61 : } else {
62 0 : length_m = 0;
63 : }
64 0 : }
65 :
66 576 : KOKKOS_INLINE_FUNCTION Index::Index(int m, int a, const Index& b)
67 576 : : first_m(b.first_m * m + a)
68 576 : , stride_m(b.stride_m * m)
69 576 : , length_m(b.length_m) {}
70 :
71 : KOKKOS_INLINE_FUNCTION Index::Index(int f, int s, const Index* b)
72 : : first_m(f)
73 : , stride_m(s)
74 : , length_m(b->length_m) {}
75 :
76 31314524 : KOKKOS_INLINE_FUNCTION int Index::first() const noexcept {
77 31314524 : return first_m;
78 : }
79 :
80 408708 : KOKKOS_INLINE_FUNCTION int Index::stride() const noexcept {
81 408708 : return stride_m;
82 : }
83 :
84 : KOKKOS_INLINE_FUNCTION bool Index::empty() const noexcept {
85 : return length_m == 0;
86 : }
87 :
88 15025108 : KOKKOS_INLINE_FUNCTION size_t Index::length() const noexcept {
89 15025108 : return length_m;
90 : }
91 :
92 548052 : KOKKOS_INLINE_FUNCTION int Index::last() const noexcept {
93 548052 : return (length_m == 0) ? first_m : first_m + stride_m * (length_m - 1);
94 : }
95 :
96 152248 : KOKKOS_INLINE_FUNCTION int Index::min() const noexcept {
97 152248 : return (stride_m >= 0) ? first_m : first_m + stride_m * (length_m - 1);
98 : }
99 :
100 154072 : KOKKOS_INLINE_FUNCTION int Index::max() const noexcept {
101 154072 : return (stride_m >= 0) ? first_m + stride_m * (length_m - 1) : first_m;
102 : }
103 :
104 55248 : KOKKOS_INLINE_FUNCTION Index& Index::operator+=(int off) {
105 55248 : first_m += off;
106 55248 : return *this;
107 : }
108 :
109 55248 : KOKKOS_INLINE_FUNCTION Index& Index::operator-=(int off) {
110 55248 : first_m -= off;
111 55248 : return *this;
112 : }
113 :
114 480 : KOKKOS_INLINE_FUNCTION Index operator+(const Index& i, int off) {
115 480 : return Index(1, off, i);
116 : }
117 :
118 : KOKKOS_INLINE_FUNCTION Index operator+(int off, const Index& i) {
119 : return Index(1, off, i);
120 : }
121 :
122 96 : KOKKOS_INLINE_FUNCTION Index operator-(const Index& i, int off) {
123 96 : return Index(1, -off, i);
124 : }
125 :
126 : KOKKOS_INLINE_FUNCTION Index operator-(int off, const Index& i) {
127 : return Index(-1, off, i);
128 : }
129 :
130 : KOKKOS_INLINE_FUNCTION Index operator-(const Index& i) {
131 : return Index(-1, 0, i);
132 : }
133 :
134 : KOKKOS_INLINE_FUNCTION Index operator*(const Index& i, int m) {
135 : return Index(m, 0, i);
136 : }
137 :
138 : KOKKOS_INLINE_FUNCTION Index operator*(int m, const Index& i) {
139 : return Index(m, 0, i);
140 : }
141 :
142 : KOKKOS_INLINE_FUNCTION Index operator/(const Index& i, int d) {
143 : return Index(i.first_m / d, i.stride_m / d, &i);
144 : }
145 :
146 0 : KOKKOS_INLINE_FUNCTION Index Index::reverse() const {
147 0 : Index j;
148 0 : j.first_m = last();
149 0 : j.length_m = length_m;
150 0 : j.stride_m = -stride_m;
151 0 : return j;
152 : }
153 :
154 69072 : KOKKOS_INLINE_FUNCTION bool Index::touches(const Index& a) const {
155 69072 : return (min() <= a.max()) && (max() >= a.min());
156 : }
157 :
158 : KOKKOS_INLINE_FUNCTION bool Index::contains(const Index& a) const {
159 : return (min() <= a.min()) && (max() >= a.max());
160 : }
161 :
162 1324 : KOKKOS_INLINE_FUNCTION bool Index::split(Index& l, Index& r) const {
163 1324 : PAssert_EQ(stride_m, 1);
164 1324 : PAssert_GT(length_m, 1);
165 1324 : auto mid = first_m + length_m / 2 - 1;
166 1324 : l = Index(first_m, mid);
167 1324 : r = Index(mid + 1, first_m + length_m - 1);
168 1324 : return true;
169 : }
170 :
171 48 : KOKKOS_INLINE_FUNCTION bool Index::split(Index& l, Index& r, int i) const {
172 48 : PAssert_EQ(stride_m, 1);
173 48 : PAssert_GT(length_m, 1);
174 48 : if (i >= first_m + static_cast<int>(length_m)) {
175 0 : return false;
176 : }
177 48 : l = Index(first_m, i);
178 48 : r = Index(i + 1, first_m + length_m - 1);
179 48 : return true;
180 : }
181 :
182 0 : KOKKOS_INLINE_FUNCTION bool Index::split(Index& l, Index& r, double a) const {
183 0 : PAssert_EQ(stride_m, 1);
184 0 : PAssert_GT(length_m, 1);
185 0 : PAssert_LT(a, 1.0);
186 0 : PAssert_GT(a, 0.0);
187 0 : int mid = first_m + static_cast<int>(length_m * a + 0.5) - 1;
188 0 : l = Index(first_m, mid);
189 0 : r = Index(mid + 1, first_m + length_m - 1);
190 0 : return true;
191 : }
192 :
193 : //////////////////////////////////////////////////////////////////////
194 : // Calculate the lowest common multiple of s1 and s2.
195 : // put the result in s.
196 : // Also calculate m1 = s/s1 and m2 = s/s2.
197 : // This version is optimized for small s1 and s2 and
198 : // just uses an exhaustive search.
199 : //////////////////////////////////////////////////////////////////////
200 0 : KOKKOS_INLINE_FUNCTION void lcm(int s1, int s2, int& s, int& m1, int& m2) {
201 0 : PAssert_GT(s1, 0); // For simplicity, make some assumptions.
202 0 : PAssert_GT(s2, 0);
203 0 : int i1 = s1;
204 0 : int i2 = s2;
205 0 : int _m1 = 1;
206 0 : int _m2 = 1;
207 0 : if (i2 < i1) {
208 : while (true) {
209 0 : while (i2 < i1) {
210 0 : i2 += s2;
211 0 : ++_m2;
212 : }
213 0 : if (i1 == i2) {
214 0 : m1 = _m1;
215 0 : m2 = _m2;
216 0 : s = i1;
217 0 : return;
218 : }
219 0 : i1 += s1;
220 0 : ++_m1;
221 : }
222 : } else {
223 : while (true) {
224 0 : while (i1 < i2) {
225 0 : i1 += s1;
226 0 : ++_m1;
227 : }
228 0 : if (i1 == i2) {
229 0 : m1 = _m1;
230 0 : m2 = _m2;
231 0 : s = i1;
232 0 : return;
233 : }
234 0 : i2 += s2;
235 0 : ++_m2;
236 : }
237 : }
238 : }
239 :
240 : //
241 : // Intersect, with the code for the common case of
242 : // both strides equal to one.
243 : //
244 204312 : KOKKOS_INLINE_FUNCTION Index Index::intersect(const Index& rhs) const {
245 204312 : Index ret;
246 204312 : if ((stride() == 1) && (rhs.stride() == 1)) {
247 204312 : int lf = first();
248 204312 : int rf = rhs.first();
249 204312 : int ll = last();
250 204312 : int rl = rhs.last();
251 204312 : int f = lf > rf ? lf : rf;
252 204312 : int l = ll < rl ? ll : rl;
253 204312 : ret.first_m = f;
254 204312 : ret.length_m = ((l >= f) ? l - f + 1 : 0);
255 204312 : ret.stride_m = 1;
256 : } else
257 0 : ret = general_intersect(rhs);
258 204312 : return ret;
259 : }
260 :
261 140728 : KOKKOS_INLINE_FUNCTION Index Index::grow(int ncells) const {
262 140728 : Index index;
263 :
264 140728 : index.first_m = this->first_m - ncells;
265 140728 : index.length_m = this->length_m + 2 * ncells;
266 140728 : index.stride_m = this->stride_m;
267 :
268 140728 : return index;
269 : }
270 :
271 0 : KOKKOS_INLINE_FUNCTION static Index do_intersect(const Index& a, const Index& b) {
272 0 : PAssert_GT(a.stride(), 0); // This should be assured by the
273 0 : PAssert_GT(b.stride(), 0); // caller of this function.
274 :
275 : int newStride; // The stride for the new index is
276 : int a_mul, b_mul; // a_mul=newStride/a.stride() ...
277 0 : lcm(a.stride(), b.stride(), // The input strides...
278 : newStride, a_mul, b_mul); // the lcm of the strides of a and b.
279 :
280 : // Find the offset from a.first() in units of newStride
281 : // that puts the ranges close together.
282 0 : int a_i = (b.first() - a.first()) / a.stride();
283 0 : int a_off = a.first() + a_i * a.stride();
284 0 : if (a_off < b.first()) {
285 0 : a_i++;
286 0 : a_off += a.stride();
287 : }
288 :
289 0 : PAssert_GE(a_off, b.first()); // make sure I'm understanding this right...
290 :
291 : // Now do an exhaustive search for the first point in common.
292 : // Count over all possible offsets for a.
293 0 : for (int a_m = 0; (a_m < a_mul) && (a_i < (int)a.length());
294 0 : a_m++, a_i++, a_off += a.stride()) {
295 0 : int b_off = b.first();
296 : // Count over all possible offsets for b.
297 0 : for (int b_m = 0; (b_m < b_mul) && (b_m < (int)b.length()); b_m++, b_off += b.stride())
298 0 : if (a_off == b_off) { // If the offsets are the same, we found it!
299 0 : int am = a.max(); // Find the minimum maximum of a and b...
300 0 : int bm = b.max();
301 0 : int m = am < bm ? am : bm;
302 0 : return Index(a_off, m, newStride);
303 : }
304 : }
305 0 : return Index(0); // If we get to here there is no intersection.
306 : }
307 :
308 0 : KOKKOS_INLINE_FUNCTION Index Index::general_intersect(const Index& that) const {
309 : // If they just don't overlap, return null indexes.
310 0 : if ((min() > that.max()) || (that.min() > max()))
311 0 : return Index(0);
312 0 : if ((stride_m == 0) || (that.stride_m == 0))
313 0 : return Index(0);
314 :
315 : // If one or the other counts -ve, reverse it and intersect result.
316 0 : if (that.stride_m < 0)
317 0 : return intersect(that.reverse());
318 0 : if (stride_m < 0) {
319 0 : Index r;
320 0 : r = reverse().intersect(that).reverse();
321 0 : return r;
322 : }
323 :
324 : // Getting closer to the real thing: intersect them.
325 : // Pass the one that starts lower as the first argument.
326 0 : Index r;
327 0 : if (first_m < that.first_m)
328 0 : r = do_intersect(*this, that);
329 : else
330 0 : r = do_intersect(that, *this);
331 :
332 0 : return r;
333 : }
334 : } // namespace ippl
|
