3 {
4 using double4 = Point4D;
5 using double2 = Point2D;
6 using int2 = std::pair<int, int>;
7
8
9
10#pragma omp parallel
11 {
12 int i, j, ch;
13
14 double4 lu[2];
15
16 double2 mom0;
17 double2 mom1;
18 double2 mom2;
19 double2 mom3;
20 double2 mom4;
21 double2 mom5;
22 double2 mom6;
23 double2 mom7;
24 double2 mom8;
25 double2 mom9;
26 double2 mom10;
27 double2 mom11;
28 double2 mom12;
29 double2 mom13;
30 double2 mom14;
31 double2 cen;
32 double2 dr;
33
34 int m[2];
35 int cm;
36 const int nnodes = 2 * (int)object.size() - 1;
37 const int nbodies = (int)object.size();
38#pragma omp for schedule(dynamic,5)
39 for (int k = nbodies; k < nnodes; ++k)
40 {
41
42
43
44
45
46
47
48
49
50
51
52 j = 0;
53 cm = 0;
54
55
56 while (cm == 0)
57 {
58 j = 2;
59 int srt = indexSort[(nnodes - 1) - k];
60 int2 chdPair = child[srt];
61
62 for (i = 0; i < 2; i++) {
63 int chd = i * chdPair.second + (1 - i) * chdPair.first;
64
65 ch = (chd >= nbodies) ? (chd - nbodies) : ((nnodes - 1) - indexSortT[chd]);
66
67 if ((chd >= nbodies) || (mass[nnodes - 1 - ch] >= 0))
68 j--;
69 }
70
71 if (j == 0)
72 {
73
75
76
77
78 for (i = 0; i < 2; i++)
79 {
80 const int chd = i * chdPair.second + (1 - i) * chdPair.first;
81 if (chd >= nbodies)
82 {
83 ch = chd - nbodies;
84 const int sortedBody = mortonCodesIdx[ch];
85
86 double4 xyAB = gabForLeaves[sortedBody];
87 lu[i] = double4{
88 ::fmin(xyAB[0], xyAB[2]), ::fmin(xyAB[1], xyAB[3]),
89 ::fmax(xyAB[0], xyAB[2]), ::fmax(xyAB[1], xyAB[3])
90 };
91 }
92 else
93 {
94 ch = indexSortT[chd];
95
96 lu[i] = lowerupper[chd];
97 }
98 }
99
100
101 lowerupper[srt] = double4{
102 ::fmin(lu[0][0], lu[1][0]), ::fmin(lu[0][1], lu[1][1]),
103 ::fmax(lu[0][2], lu[1][2]), ::fmax(lu[0][3], lu[1][3])
104 };
105
106 cen = center[srt] = double2{ 0.5 * (lowerupper[srt][0] + lowerupper[srt][2]),
107 0.5 * (lowerupper[srt][1] + lowerupper[srt][3]) };
108
109 const double2 zero = { 0.0, 0.0 };
110 double2 momh0 = zero;
111 double2 momh1 = zero;
112 double2 momh2 = zero;
113 double2 momh3 = zero;
114 double2 momh4 = zero;
115 double2 momh5 = zero;
116 double2 momh6 = zero;
117 double2 momh7 = zero;
118 double2 momh8 = zero;
119 double2 momh9 = zero;
120 double2 momh10 = zero;
121 double2 momh11 = zero;
122 double2 momh12 = zero;
123 double2 momh13 = zero;
124 double2 momh14 = zero;
125
126 for (i = 0; i < 2; i++)
127 {
128 const int chd = i * chdPair.second + (1 - i) * chdPair.first;
129 if (chd >= nbodies)
130 {
131 ch = chd - nbodies;
132 const int sortedBody = mortonCodesIdx[ch];
133
135 {
136 mom0 = double2{ gamma[sortedBody], 0.0 };
137
138mom1 = mom2 = mom3 = mom4 = mom5 = mom6 = mom7 = mom8 = mom9 = mom10 = mom11 = mom12 = mom13 = mom14 = double2{ 0.0, 0.0 };
139 double2 pos = object[sortedBody];
140 dr = pos - cen;
141 m[i] = 1;
142 }
143
145 {
146
147 }
148
149 }
150 else
151 {
152 const int srtT = indexSortT[chd];
153 ch = (nnodes - 1) - srtT;
154
156
157 mom0 = double2{ moms[nch + 0][0], (double)0 };
158 mom1 = moms[nch + 1];
159 mom2 = moms[nch + 2];
160 mom3 = moms[nch + 3];
161 mom4 = moms[nch + 4];
162 mom5 = moms[nch + 5];
163 mom6 = moms[nch + 6];
164 mom7 = moms[nch + 7];
165 mom8 = moms[nch + 8];
166 mom9 = moms[nch + 9];
167 mom10 = moms[nch + 10];
168 mom11 = moms[nch + 11];
169 mom12 = moms[nch + 12];
170 mom13 = moms[nch + 13];
171 mom14 = moms[nch + 14];
172 dr = center[chd] - cen;
173 m[i] = mass[srtT];
174 }
175
176
177 momh0 += mom0;
178 momh1 += mom1;
179 momh2 += mom2;
180 momh3 += mom3;
181 momh4 += mom4;
182 momh5 += mom5;
183 momh6 += mom6;
184 momh7 += mom7;
185 momh8 += mom8;
186 momh9 += mom9;
187 momh10 += mom10;
188 momh11 += mom11;
189 momh12 += mom12;
190 momh13 += mom13;
191 momh14 += mom14;
192
193 double2 z = dr;
194
195 momh1 +=
multz(mom0, z);
196 momh2 += 2 *
multz(mom1, z);
197 momh3 += 3 *
multz(mom2, z);
198 momh4 += 4 *
multz(mom3, z);
199 momh5 += 5 *
multz(mom4, z);
200 momh6 += 6 *
multz(mom5, z);
201 momh7 += 7 *
multz(mom6, z);
202 momh8 += 8 *
multz(mom7, z);
203 momh9 += 9 *
multz(mom8, z);
204 momh10 += 10 *
multz(mom9, z);
205 momh11 += 11 *
multz(mom10, z);
206 momh12 += 12 *
multz(mom11, z);
207 momh13 += 13 *
multz(mom12, z);
208 momh14 += 14 *
multz(mom13, z);
209
211
212 momh2 +=
multz(mom0, z);
213 momh3 += 3 *
multz(mom1, z);
214 momh4 += 6 *
multz(mom2, z);
215 momh5 += 10 *
multz(mom3, z);
216 momh6 += 15 *
multz(mom4, z);
217 momh7 += 21 *
multz(mom5, z);
218 momh8 += 28 *
multz(mom6, z);
219 momh9 += 36 *
multz(mom7, z);
220 momh10 += 45 *
multz(mom8, z);
221 momh11 += 55 *
multz(mom9, z);
222 momh12 += 66 *
multz(mom10, z);
223 momh13 += 78 *
multz(mom11, z);
224 momh14 += 91 *
multz(mom12, z);
225
227
228 momh3 +=
multz(mom0, z);
229 momh4 += 4 *
multz(mom1, z);
230 momh5 += 10 *
multz(mom2, z);
231 momh6 += 20 *
multz(mom3, z);
232 momh7 += 35 *
multz(mom4, z);
233 momh8 += 56 *
multz(mom5, z);
234 momh9 += 84 *
multz(mom6, z);
235 momh10 += 120 *
multz(mom7, z);
236 momh11 += 165 *
multz(mom8, z);
237 momh12 += 220 *
multz(mom9, z);
238 momh13 += 286 *
multz(mom10, z);
239 momh14 += 364 *
multz(mom11, z);
240
242
243 momh4 +=
multz(mom0, z);
244 momh5 += 5 *
multz(mom1, z);
245 momh6 += 15 *
multz(mom2, z);
246 momh7 += 35 *
multz(mom3, z);
247 momh8 += 70 *
multz(mom4, z);
248 momh9 += 126 *
multz(mom5, z);
249 momh10 += 210 *
multz(mom6, z);
250 momh11 += 330 *
multz(mom7, z);
251 momh12 += 495 *
multz(mom8, z);
252 momh13 += 715 *
multz(mom9, z);
253 momh14 += 1001 *
multz(mom10, z);
254
256
257 momh5 +=
multz(mom0, z);
258 momh6 += 6 *
multz(mom1, z);
259 momh7 += 21 *
multz(mom2, z);
260 momh8 += 56 *
multz(mom3, z);
261 momh9 += 126 *
multz(mom4, z);
262 momh10 += 252 *
multz(mom5, z);
263 momh11 += 462 *
multz(mom6, z);
264 momh12 += 792 *
multz(mom7, z);
265 momh13 += 1287 *
multz(mom8, z);
266 momh14 += 2002 *
multz(mom9, z);
267
269
270 momh6 +=
multz(mom0, z);
271 momh7 += 7 *
multz(mom1, z);
272 momh8 += 28 *
multz(mom2, z);
273 momh9 += 84 *
multz(mom3, z);
274 momh10 += 210 *
multz(mom4, z);
275 momh11 += 462 *
multz(mom5, z);
276 momh12 += 924 *
multz(mom6, z);
277 momh13 += 1716 *
multz(mom7, z);
278 momh14 += 3003 *
multz(mom8, z);
279
281
282 momh7 +=
multz(mom0, z);
283 momh8 += 8 *
multz(mom1, z);
284 momh9 += 36 *
multz(mom2, z);
285 momh10 += 120 *
multz(mom3, z);
286 momh11 += 330 *
multz(mom4, z);
287 momh12 += 792 *
multz(mom5, z);
288 momh13 += 1716 *
multz(mom6, z);
289 momh14 += 3432 *
multz(mom7, z);
290
292
293 momh8 +=
multz(mom0, z);
294 momh9 += 9 *
multz(mom1, z);
295 momh10 += 45 *
multz(mom2, z);
296 momh11 += 165 *
multz(mom3, z);
297 momh12 += 495 *
multz(mom4, z);
298 momh13 += 1287 *
multz(mom5, z);
299 momh14 += 3003 *
multz(mom6, z);
300
302
303 momh9 +=
multz(mom0, z);
304 momh10 += 10 *
multz(mom1, z);
305 momh11 += 55 *
multz(mom2, z);
306 momh12 += 220 *
multz(mom3, z);
307 momh13 += 715 *
multz(mom4, z);
308 momh14 += 2002 *
multz(mom5, z);
309
311
312 momh10 +=
multz(mom0, z);
313 momh11 += 11 *
multz(mom1, z);
314 momh12 += 66 *
multz(mom2, z);
315 momh13 += 286 *
multz(mom3, z);
316 momh14 += 1001 *
multz(mom4, z);
317
319
320 momh11 +=
multz(mom0, z);
321 momh12 += 12 *
multz(mom1, z);
322 momh13 += 78 *
multz(mom2, z);
323 momh14 += 364 *
multz(mom3, z);
324
326
327 momh12 +=
multz(mom0, z);
328 momh13 += 13 *
multz(mom1, z);
329 momh14 += 91 *
multz(mom2, z);
330
332
333 momh13 +=
multz(mom0, z);
334 momh14 += 14 *
multz(mom1, z);
335
337
338 momh14 +=
multz(mom0, z);
339 }
340
341
342 momh0[1] = 0;
343 moms[kch + 0] = momh0;
344 moms[kch + 1] = momh1;
345 moms[kch + 2] = momh2;
346 moms[kch + 3] = momh3;
347 moms[kch + 4] = momh4;
348 moms[kch + 5] = momh5;
349 moms[kch + 6] = momh6;
350 moms[kch + 7] = momh7;
351 moms[kch + 8] = momh8;
352 moms[kch + 9] = momh9;
353 moms[kch + 10] = momh10;
354 moms[kch + 11] = momh11;
355 moms[kch + 12] = momh12;
356 moms[kch + 13] = momh13;
357 moms[kch + 14] = momh14;
358 cm = m[0] + m[1];
359 }
360
361#pragma omp flush
362
363 if (cm != 0)
364 {
365 mass[nnodes - 1 - k] = cm;
366
367
368 }
369 }
370 }
371 }
372 }
Point2D multz(const Point2D &a, const Point2D &b)
Умножение комплексных чисел
1.9.8