0d60813b23b7798549634c50d4827ea340111ef9
[charm.git] / examples / charm++ / cell / md / pairCompute.ci
1 module pairCompute {
2
3   // Include md_config.h on the accelerator (for physics constants)
4   accelblock { #include "md_config.h" };
5
6   array[2D] PairCompute {
7
8     entry PairCompute();
9
10     entry void init(int numParticlesPerPatch);
11
12     entry void patchData(int numParticles,
13                          MD_FLOAT particleX[numParticles],
14                          MD_FLOAT particleY[numParticles],
15                          MD_FLOAT particleZ[numParticles],
16                          MD_FLOAT particleQ[numParticles],
17                          int fromPatch
18                         );
19
20     entry [accel] void doCalc()[  readonly : int numParticles <impl_obj->numParticles>,
21                                   readonly : int thisIndex_x <impl_obj->thisIndex.x>,
22                                   readonly : int thisIndex_y <impl_obj->thisIndex.y>,
23                                   readonly : MD_FLOAT p0_x[numParticles] <impl_obj->particleX[0]>,
24                                   readonly : MD_FLOAT p1_x[numParticles] <impl_obj->particleX[1]>,
25                                   readonly : MD_FLOAT p0_y[numParticles] <impl_obj->particleY[0]>,
26                                   readonly : MD_FLOAT p1_y[numParticles] <impl_obj->particleY[1]>,
27                                   readonly : MD_FLOAT p0_z[numParticles] <impl_obj->particleZ[0]>,
28                                   readonly : MD_FLOAT p1_z[numParticles] <impl_obj->particleZ[1]>,
29                                   readonly : MD_FLOAT p0_q[numParticles] <impl_obj->particleQ[0]>,
30                                   readonly : MD_FLOAT p1_q[numParticles] <impl_obj->particleQ[1]>,
31                                  writeonly : MD_FLOAT f0_x[numParticles] <impl_obj->forceX[0]>,
32                                  writeonly : MD_FLOAT f1_x[numParticles] <impl_obj->forceX[1]>,
33                                  writeonly : MD_FLOAT f0_y[numParticles] <impl_obj->forceY[0]>,
34                                  writeonly : MD_FLOAT f1_y[numParticles] <impl_obj->forceY[1]>,
35                                  writeonly : MD_FLOAT f0_z[numParticles] <impl_obj->forceZ[0]>,
36                                  writeonly : MD_FLOAT f1_z[numParticles] <impl_obj->forceZ[1]>,
37                                  writeonly : unsigned int localFlopCount <impl_obj->localFlopCount>
38                                ] {
39
40       // DMK - DEBUG
41       #if (ENABLE_USER_EVENTS != 0) && defined(CMK_CELL) && (CMK_CELL == 0)
42         if (numParticles != 256) { CkPrintf("[DEBUG] :: PairCompute::doCalc() - Called with numParticles:%d...\n", numParticles); }
43         double __start_time__ = CmiWallTimer();
44       #endif
45
46       // Calculate the electrostatic force (coulumb's law) between the particles
47       //   F = (k_e * (q_0 * q_1)) / (r^2)
48
49       // DMK - DEBUG
50       #if COUNT_FLOPS != 0
51         localFlopCount = 0;
52       #endif
53
54       register MD_VEC* p1_x_vec = (MD_VEC*)p1_x;
55       register MD_VEC* p1_y_vec = (MD_VEC*)p1_y;
56       register MD_VEC* p1_z_vec = (MD_VEC*)p1_z;
57       register MD_VEC* p1_q_vec = (MD_VEC*)p1_q;
58       register MD_VEC* f1_x_vec = (MD_VEC*)f1_x;
59       register MD_VEC* f1_y_vec = (MD_VEC*)f1_y;
60       register MD_VEC* f1_z_vec = (MD_VEC*)f1_z;
61       register int i;
62       register int j;
63       register const int numParticlesByVecSize = numParticles / myvec_numElems;
64
65       // Zero out the force array for p1 (p0's force array will be zero'ed for each
66       //   particle in the outer loop)
67       MD_VEC zero_vec = vspread_MDF(zero);
68       for (j = 0; j < numParticlesByVecSize; j++) {
69         f1_x_vec[j] = zero_vec;
70         f1_y_vec[j] = zero_vec;
71         f1_z_vec[j] = zero_vec;
72       }
73
74       // Spread coulumb's constant across a vector
75       MD_VEC coulomb_vec = vspread_MDF(COULOMBS_CONSTANT);
76
77       // Outer-loop (p0)
78       for (i = 0; i < numParticles; i++) {
79
80         // Spread this particle's values out over vectors
81         MD_VEC p0_x_i_vec = vspread_MDF(p0_x[i]);
82         MD_VEC p0_y_i_vec = vspread_MDF(p0_y[i]);
83         MD_VEC p0_z_i_vec = vspread_MDF(p0_z[i]);
84         MD_VEC p0_q_i_vec = vspread_MDF(p0_q[i]);
85         MD_VEC f0_x_i_vec = vspread_MDF(zero);
86         MD_VEC f0_y_i_vec = vspread_MDF(zero);
87         MD_VEC f0_z_i_vec = vspread_MDF(zero);
88
89         #if 1
90
91         // Multi-stage software pipe
92         //   1) Load particle data
93         //   2) Calculate radius
94         //   3) Caluclate force vectors
95         //   4) Apply force vector
96
97         /// Loop Preamble ///
98
99         const int loopUpper = numParticlesByVecSize - 3;
100         // Setup Input for Stage 2 - Initially for j+2
101         const int jPlus2 = 2;
102         MD_VEC p1_x_j_vec = p1_x_vec[jPlus2];
103         MD_VEC p1_y_j_vec = p1_y_vec[jPlus2];
104         MD_VEC p1_z_j_vec = p1_z_vec[jPlus2];
105         MD_VEC p1_q_j_vec = p1_q_vec[jPlus2];
106
107         // Setup Input for Stage 3 - Initial for j+1
108         const int jPlus1 = 1;
109         MD_VEC p_x_diff_vec = p0_x_i_vec - p1_x_vec[jPlus1];
110         MD_VEC p_y_diff_vec = p0_y_i_vec - p1_y_vec[jPlus1];
111         MD_VEC p_z_diff_vec = p0_z_i_vec - p1_z_vec[jPlus1];
112         MD_VEC p1_q_j_vec_s2 = p1_q_vec[jPlus1];
113         MD_VEC r_2_vec = (p_x_diff_vec * p_x_diff_vec) + (p_y_diff_vec * p_y_diff_vec) + (p_z_diff_vec * p_z_diff_vec);
114         MD_VEC r_vec = vsqrt_MDF(r_2_vec);
115
116         // DMK - DEBUG
117         #if COUNT_FLOPS != 0
118           localFlopCount += (9 * 4);
119         #endif
120
121         // Setup Input for Stage 4 - Initially for j+0
122         const int jPlus0 = 0;
123         MD_VEC f_x_vec;
124         MD_VEC f_y_vec;
125         MD_VEC f_z_vec;
126         {
127           MD_VEC p_x_diff_vec = p0_x_i_vec - p1_x_vec[jPlus0];
128           MD_VEC p_y_diff_vec = p0_y_i_vec - p1_y_vec[jPlus0];
129           MD_VEC p_z_diff_vec = p0_z_i_vec - p1_z_vec[jPlus0];
130           MD_VEC r_2_vec = (p_x_diff_vec * p_x_diff_vec) + (p_y_diff_vec * p_y_diff_vec) + (p_z_diff_vec * p_z_diff_vec);
131           MD_VEC r_vec = vsqrt_MDF(r_2_vec);
132           MD_VEC r_recip_vec = vrecip_MDF(r_vec);
133           MD_VEC r_2_recip_vec = vrecip_MDF(r_2_vec);
134           MD_VEC p_x_diff_norm_vec = p_x_diff_vec * r_recip_vec;
135           MD_VEC p_y_diff_norm_vec = p_y_diff_vec * r_recip_vec;
136           MD_VEC p_z_diff_norm_vec = p_z_diff_vec * r_recip_vec;
137           MD_VEC f_mag_vec = coulomb_vec * ((p0_q_i_vec * p1_q_vec[jPlus0]) * r_2_recip_vec);
138           f_x_vec = p_x_diff_norm_vec * f_mag_vec;
139           f_y_vec = p_y_diff_norm_vec * f_mag_vec;
140           f_z_vec = p_z_diff_norm_vec * f_mag_vec;
141
142           // DMK - DEBUG
143           #if COUNT_FLOPS != 0
144             localFlopCount += (23 * 4);
145           #endif
146         }
147
148         /// Main Loop ///
149
150         for (j = 0; j < loopUpper; j++) {
151
152           #if 1
153
154           // Stage 4 - Apply force vector
155           f0_x_i_vec = vadd_MDF(f0_x_i_vec, f_x_vec);
156           f0_y_i_vec = vadd_MDF(f0_y_i_vec, f_y_vec);
157           f0_z_i_vec = vadd_MDF(f0_z_i_vec, f_z_vec);
158           f1_x_vec[j] = vsub_MDF(f1_x_vec[j], f_x_vec);
159           f1_y_vec[j] = vsub_MDF(f1_y_vec[j], f_y_vec);
160           f1_z_vec[j] = vsub_MDF(f1_z_vec[j], f_z_vec);
161
162           // Stage 3 - Calculate force vector
163           MD_VEC r_recip_vec = vrecip_MDF(r_vec);
164           MD_VEC r_2_recip_vec = vrecip_MDF(r_2_vec);
165           MD_VEC p_z_diff_norm_vec = vmul_MDF(p_z_diff_vec, r_recip_vec);
166           MD_VEC p_y_diff_norm_vec = vmul_MDF(p_y_diff_vec, r_recip_vec);
167           MD_VEC p_x_diff_norm_vec = vmul_MDF(p_x_diff_vec, r_recip_vec);
168           MD_VEC f_mag_vec = vmul_MDF(vmul_MDF(coulomb_vec, p0_q_i_vec), vmul_MDF(p1_q_j_vec_s2, r_2_recip_vec));
169           f_x_vec = vmul_MDF(p_x_diff_norm_vec, f_mag_vec);
170           f_y_vec = vmul_MDF(p_y_diff_norm_vec, f_mag_vec);
171           f_z_vec = vmul_MDF(p_z_diff_norm_vec, f_mag_vec);
172
173           // Stage 2 - Calculate radius
174           p1_q_j_vec_s2 = p1_q_j_vec;
175           p_x_diff_vec = vsub_MDF(p0_x_i_vec, p1_x_j_vec);
176           p_y_diff_vec = vsub_MDF(p0_y_i_vec, p1_y_j_vec);
177           p_z_diff_vec = vsub_MDF(p0_z_i_vec, p1_z_j_vec);
178           r_2_vec = vmadd_MDF(p_x_diff_vec, p_x_diff_vec, vmadd_MDF(p_y_diff_vec, p_y_diff_vec, vmul_MDF(p_z_diff_vec, p_z_diff_vec)));
179           r_vec = vsqrt_MDF(r_2_vec);
180
181           // Stage 1 - Load particle data
182           const int jIndexToLoad = j + 3;
183           p1_x_j_vec = p1_x_vec[jIndexToLoad];
184           p1_y_j_vec = p1_y_vec[jIndexToLoad];
185           p1_z_j_vec = p1_z_vec[jIndexToLoad];
186           p1_q_j_vec = p1_q_vec[jIndexToLoad];
187
188           #else
189
190           // Stage 4 - Apply force vector
191           f0_x_i_vec += f_x_vec;
192           f0_y_i_vec += f_y_vec;
193           f0_z_i_vec += f_z_vec;
194           f1_x_vec[j] -= f_x_vec;
195           f1_y_vec[j] -= f_y_vec;
196           f1_z_vec[j] -= f_z_vec;
197
198           // Stage 3 - Calculate force vector
199           MD_VEC r_recip_vec = vrecip_MDF(r_vec);
200           MD_VEC r_2_recip_vec = vrecip_MDF(r_2_vec);
201           MD_VEC p_x_diff_norm_vec = p_x_diff_vec * r_recip_vec;
202           MD_VEC p_y_diff_norm_vec = p_y_diff_vec * r_recip_vec;
203           MD_VEC p_z_diff_norm_vec = p_z_diff_vec * r_recip_vec;
204           MD_VEC f_mag_vec = coulomb_vec * p0_q_i_vec * p1_q_j_vec_s2 * r_2_recip_vec;  // Calc force magnitude
205           f_x_vec = p_x_diff_norm_vec * f_mag_vec;  // Multiply normalized vector by force magnitude
206           f_y_vec = p_y_diff_norm_vec * f_mag_vec;
207           f_z_vec = p_z_diff_norm_vec * f_mag_vec;
208
209           // Stage 2 - Calculate radius
210           p1_q_j_vec_s2 = p1_q_j_vec;
211           p_x_diff_vec = p0_x_i_vec - p1_x_j_vec;
212           p_y_diff_vec = p0_y_i_vec - p1_y_j_vec;
213           p_z_diff_vec = p0_z_i_vec - p1_z_j_vec;
214           r_2_vec = (p_x_diff_vec * p_x_diff_vec) + (p_y_diff_vec * p_y_diff_vec) + (p_z_diff_vec * p_z_diff_vec);
215           r_vec = vsqrt_MDF(r_2_vec);
216
217           // Stage 1 - Load particle data
218           const int jIndexToLoad = j + 3;
219           p1_x_j_vec = p1_x_vec[jIndexToLoad];
220           p1_y_j_vec = p1_y_vec[jIndexToLoad];
221           p1_z_j_vec = p1_z_vec[jIndexToLoad];
222           p1_q_j_vec = p1_q_vec[jIndexToLoad];
223
224           #endif
225
226           // DMK - DEBUG
227           #if COUNT_FLOPS != 0
228             localFlopCount += (31 * 4);
229           #endif
230         }
231
232         /// Loop Cleanup ///
233
234         // Complete Calc for Stage 3 Data - Apply Force Vector
235         {
236           f0_x_i_vec += f_x_vec;
237           f0_y_i_vec += f_y_vec;
238           f0_z_i_vec += f_z_vec;
239           f1_x_vec[loopUpper] -= f_x_vec;
240           f1_y_vec[loopUpper] -= f_y_vec;
241           f1_z_vec[loopUpper] -= f_z_vec;
242
243           // DMK - DEBUG
244           #if COUNT_FLOPS != 0
245             localFlopCount += (6 * 4);
246           #endif
247         }
248
249         // Complete Calc for Stage 2 Data - Calculate Force Vector, Apply Force Vector
250         {
251           MD_VEC r_recip_vec = vrecip_MDF(r_vec);
252           MD_VEC r_2_recip_vec = vrecip_MDF(r_2_vec);
253           MD_VEC p_x_diff_norm_vec = p_x_diff_vec * r_recip_vec;
254           MD_VEC p_y_diff_norm_vec = p_y_diff_vec * r_recip_vec;
255           MD_VEC p_z_diff_norm_vec = p_z_diff_vec * r_recip_vec;
256           MD_VEC f_mag_vec = coulomb_vec * ((p0_q_i_vec * p1_q_j_vec_s2) * r_2_recip_vec);
257           MD_VEC f_x_vec = p_x_diff_norm_vec * f_mag_vec;
258           MD_VEC f_y_vec = p_y_diff_norm_vec * f_mag_vec;
259           MD_VEC f_z_vec = p_z_diff_norm_vec * f_mag_vec;
260           f0_x_i_vec += f_x_vec;
261           f0_y_i_vec += f_y_vec;
262           f0_z_i_vec += f_z_vec;
263           f1_x_vec[loopUpper + 1] -= f_x_vec;
264           f1_y_vec[loopUpper + 1] -= f_y_vec;
265           f1_z_vec[loopUpper + 1] -= f_z_vec;
266
267           // DMK - DEBUG
268           #if COUNT_FLOPS != 0
269             localFlopCount += (19 * 4);
270           #endif
271         }
272
273         // Complete Calc for Stage 1 Data - Calculate Radius, Calculate Force Vector, Apply Force Vector
274         {
275           MD_VEC p_x_diff_vec = p0_x_i_vec - p1_x_j_vec;
276           MD_VEC p_y_diff_vec = p0_y_i_vec - p1_y_j_vec;
277           MD_VEC p_z_diff_vec = p0_z_i_vec - p1_z_j_vec;
278           MD_VEC r_2_vec = (p_x_diff_vec * p_x_diff_vec) + (p_y_diff_vec * p_y_diff_vec) + (p_z_diff_vec * p_z_diff_vec);
279           MD_VEC r_vec = vsqrt_MDF(r_2_vec);
280           MD_VEC r_recip_vec = vrecip_MDF(r_vec);
281           MD_VEC r_2_recip_vec = vrecip_MDF(r_2_vec);
282           MD_VEC p_x_diff_norm_vec = p_x_diff_vec * r_recip_vec;
283           MD_VEC p_y_diff_norm_vec = p_y_diff_vec * r_recip_vec;
284           MD_VEC p_z_diff_norm_vec = p_z_diff_vec * r_recip_vec;
285           MD_VEC f_mag_vec = coulomb_vec * ((p0_q_i_vec * p1_q_j_vec) * r_2_recip_vec);
286           MD_VEC f_x_vec = p_x_diff_norm_vec * f_mag_vec;
287           MD_VEC f_y_vec = p_y_diff_norm_vec * f_mag_vec;
288           MD_VEC f_z_vec = p_z_diff_norm_vec * f_mag_vec;
289           f0_x_i_vec += f_x_vec;
290           f0_y_i_vec += f_y_vec;
291           f0_z_i_vec += f_z_vec;
292           f1_x_vec[loopUpper + 2] -= f_x_vec;
293           f1_y_vec[loopUpper + 2] -= f_y_vec;
294           f1_z_vec[loopUpper + 2] -= f_z_vec;
295
296           // DMK - DEBUG
297           #if COUNT_FLOPS != 0
298             localFlopCount += (27 * 4);
299           #endif
300         }
301
302         #else
303
304         // Inner-loop (p1)
305         for (j = 0; j < numParticlesByVecSize; j++) {
306
307           // Load the particles' data
308           MD_VEC p1_x_j_vec = p1_x_vec[j];
309           MD_VEC p1_y_j_vec = p1_y_vec[j];
310           MD_VEC p1_z_j_vec = p1_z_vec[j];
311           MD_VEC p1_q_j_vec = p1_q_vec[j];
312
313           // Calculate the vector between the particles
314           MD_VEC p_x_diff_vec = p0_x_i_vec - p1_x_j_vec;
315           MD_VEC p_y_diff_vec = p0_y_i_vec - p1_y_j_vec;
316           MD_VEC p_z_diff_vec = p0_z_i_vec - p1_z_j_vec;
317
318           // Calculate r and r^2 between the particles
319           MD_VEC p_x_diff_2_vec = p_x_diff_vec * p_x_diff_vec;
320           MD_VEC p_y_diff_2_vec = p_y_diff_vec * p_y_diff_vec;
321           MD_VEC p_z_diff_2_vec = p_z_diff_vec * p_z_diff_vec;
322           MD_VEC r_2_vec = p_x_diff_2_vec + p_y_diff_2_vec + p_z_diff_2_vec;
323           MD_VEC r_vec = vsqrt_MDF(r_2_vec);
324
325           // Normalize the vector between the particles
326           MD_VEC p_x_diff_norm_vec = p_x_diff_vec / r_vec;
327           MD_VEC p_y_diff_norm_vec = p_y_diff_vec / r_vec;
328           MD_VEC p_z_diff_norm_vec = p_z_diff_vec / r_vec;
329
330           // Calculate the magnitude of the electrostatic force between the particles
331           MD_VEC f_mag_vec = coulomb_vec * ((p0_q_i_vec * p1_q_j_vec) / r_2_vec);
332
333           // Multiply the magnitude by the normalized postition difference vector to
334           //   create the force vector
335           MD_VEC f_x_vec = p_x_diff_norm_vec * f_mag_vec;
336           MD_VEC f_y_vec = p_y_diff_norm_vec * f_mag_vec;
337           MD_VEC f_z_vec = p_z_diff_norm_vec * f_mag_vec;
338
339           // Add the force to the outer-loop particle and subtract it from the inner-loop particles
340           f0_x_i_vec += f_x_vec;
341           f0_y_i_vec += f_y_vec;
342           f0_z_i_vec += f_z_vec;
343           f1_x_vec[j] -= f_x_vec;
344           f1_y_vec[j] -= f_y_vec;
345           f1_z_vec[j] -= f_z_vec;
346
347
348           // DMK - DEBUG
349           #if COUNT_FLOPS != 0
350             localFlopCount += (24 * 4);
351           #endif
352         }
353
354         #endif
355         // Write force values for the outer-loop particle back to memory
356 #ifdef USE_DOUBLE
357         f0_x[i] = vextract_MDF(f0_x_i_vec, 0) + vextract_MDF(f0_x_i_vec, 1);
358         f0_y[i] = vextract_MDF(f0_y_i_vec, 0) + vextract_MDF(f0_y_i_vec, 1);
359         f0_z[i] = vextract_MDF(f0_z_i_vec, 0) + vextract_MDF(f0_z_i_vec, 1);
360 #else
361
362         f0_x[i] = vextract_MDF(f0_x_i_vec, 0) + vextract_MDF(f0_x_i_vec, 1) + vextract_MDF(f0_x_i_vec, 2) + vextract_MDF(f0_x_i_vec, 3);
363         f0_y[i] = vextract_MDF(f0_y_i_vec, 0) + vextract_MDF(f0_y_i_vec, 1) + vextract_MDF(f0_y_i_vec, 2) + vextract_MDF(f0_y_i_vec, 3);
364         f0_z[i] = vextract_MDF(f0_z_i_vec, 0) + vextract_MDF(f0_z_i_vec, 1) + vextract_MDF(f0_z_i_vec, 2) + vextract_MDF(f0_z_i_vec, 3);
365 #endif
366 #if SANITY_CHECK
367         CkAssert(finite(f0_x[i]));
368         CkAssert(finite(f0_y[i]));
369         CkAssert(finite(f0_z[i]));
370
371 #endif  
372       }
373
374       // DMK - DEBUG
375       #if (ENABLE_USER_EVENTS != 0) && defined(CMK_CELL) && (CMK_CELL == 0)
376         double __end_time__ = CmiWallTimer();
377         double deltaTime = __end_time__ - __start_time__;
378         traceUserBracketEvent(PROJ_USER_EVENT_PAIRCOMPUTE_DOCALC_WORK, __start_time__, __end_time__);
379         if (deltaTime > 0.002) {
380           CkPrintf("[DEBUG] :: PairCompute[%d,%d]::doCalc() - Took %lf ms to execute... numParticles:%d\n"
381                    "             p0[0] = { %e, %e, %e : %e }, p1[0] = { %e, %e, %e : %e }...\n"
382                    "             f0[0] = { %e, %e, %e }, f1[0] = { %e, %e, %e }...\n",
383                    thisIndex_x, thisIndex_y, deltaTime * 100zero, numParticles,
384                    p0_x[0], p0_y[0], p0_z[0], p0_q[0], p1_x[0], p1_y[0], p1_z[0], p1_q[0],
385                    f0_x[0], f0_y[0], f0_z[0], f1_x[0], f1_y[0], f1_z[0]
386                   );
387         }
388       #endif
389
390     } doCalc_callback;
391
392   };
393
394 };