examples/charm++: Always use CkWallTimer, not CmiWallTimer (tested)
[charm.git] / examples / charm++ / cell / md / selfCompute.ci
1 module selfCompute {
2
3   // Include md_config.h on the accelerator (for physics constants)
4   accelblock { #include "md_config.h" };
5
6   array[3D] SelfCompute {
7
8     entry SelfCompute();
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                         );
18
19     entry [accel] void doCalc()[
20                                  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 particleX[numParticles] <impl_obj->particleX>,
24                                  readonly : MD_FLOAT particleY[numParticles] <impl_obj->particleY>,
25                                  readonly : MD_FLOAT particleZ[numParticles] <impl_obj->particleZ>,
26                                  readonly : MD_FLOAT particleQ[numParticles] <impl_obj->particleQ>,
27                                 writeonly : MD_FLOAT f0_x[numParticles] <impl_obj->forceX>,
28                                 writeonly : MD_FLOAT f0_y[numParticles] <impl_obj->forceY>,
29                                 writeonly : MD_FLOAT f0_z[numParticles] <impl_obj->forceZ>,
30                                 writeonly : unsigned int localFlopCount <impl_obj->localFlopCount>
31                                ] {
32
33       // DMK - DEBUG
34       #if (ENABLE_USER_EVENTS != 0) && defined(CMK_CELL) && (CMK_CELL == 0)
35         double __start_time__ = CkWallTimer();
36       #endif
37
38       // Calculate the electrostatic force (coulumb's law) between the particles
39       //   F = (k_e * (q_0 * q_1)) / (r^2)
40
41       // DMK - DEBUG
42       #if COUNT_FLOPS != 0
43         localFlopCount = 0;
44       #endif
45       const int vec_numElems= myvec_numElems;
46       const int VEC_MASK= vec_numElems-1;
47       const int shifter= sizeof(MD_FLOAT)==4 ? 2 :1;
48       register MD_VEC* p0_x_vec = (MD_VEC*)particleX;
49       register MD_VEC* p0_y_vec = (MD_VEC*)particleY;
50       register MD_VEC* p0_z_vec = (MD_VEC*)particleZ;
51       register MD_VEC* p0_q_vec = (MD_VEC*)particleQ;
52       register MD_VEC* f0_x_vec = (MD_VEC*)f0_x;
53       register MD_VEC* f0_y_vec = (MD_VEC*)f0_y;
54       register MD_VEC* f0_z_vec = (MD_VEC*)f0_z;
55       register int i;
56       register int j;
57       register const int numParticlesByVecSize = numParticles / myvec_numElems;
58
59       // Zero out the force output
60       MD_VEC zero_vec = vspread_MDF(zero);
61       for (j = 0; j < numParticlesByVecSize; j++) {
62         f0_x_vec[j] = zero_vec;
63         f0_y_vec[j] = zero_vec;
64         f0_z_vec[j] = zero_vec;
65       }
66
67       // Spread coulumb's constant across a vector
68       MD_VEC coulomb_vec = vspread_MDF(COULOMBS_CONSTANT);
69
70       // Outer-loop
71       for (i = 0; i < numParticles; i++) {
72
73         // Interact with all particles (one-by-one) until the vector boundary
74         for (j = i + 1; ((j & VEC_MASK) != 0x0) && (j < numParticles); j++) {
75
76           MD_FLOAT p_x_diff = particleX[i] - particleX[j];
77           MD_FLOAT p_y_diff = particleY[i] - particleY[j];
78           MD_FLOAT p_z_diff = particleZ[i] - particleZ[j];
79           MD_FLOAT r_2 = (p_x_diff * p_x_diff) + (p_y_diff * p_y_diff) + (p_z_diff * p_z_diff);
80           MD_FLOAT r = sqrt(r_2);
81
82           MD_FLOAT p_x_diff_norm = p_x_diff / r;
83           MD_FLOAT p_y_diff_norm = p_y_diff / r;
84           MD_FLOAT p_z_diff_norm = p_z_diff / r;
85
86           MD_FLOAT p_q_prod = particleQ[i] * particleQ[j];
87           MD_FLOAT f_mag = COULOMBS_CONSTANT * (p_q_prod / r_2);
88
89           MD_FLOAT f_x = f_mag * p_x_diff_norm;
90           MD_FLOAT f_y = f_mag * p_y_diff_norm;
91           MD_FLOAT f_z = f_mag * p_z_diff_norm;
92
93           f0_x[i] += f_x;
94           f0_y[i] += f_y;
95           f0_z[i] += f_z;
96           f0_x[j] -= f_x;
97           f0_y[j] -= f_y;
98           f0_z[j] -= f_z;
99
100           // DMK - DEBUG
101           #if COUNT_FLOPS != 0
102             localFlopCount += 26;
103           #endif
104         }
105
106         // Spread this particle's (p0[i]) values out over vectors
107         MD_VEC p0_x_i_vec = vspread_MDF(particleX[i]);
108         MD_VEC p0_y_i_vec = vspread_MDF(particleY[i]);
109         MD_VEC p0_z_i_vec = vspread_MDF(particleZ[i]);
110         MD_VEC p0_q_i_vec = vspread_MDF(particleQ[i]);
111         MD_VEC f0_x_i_vec = vspread_MDF(zero);
112         MD_VEC f0_y_i_vec = vspread_MDF(zero);
113         MD_VEC f0_z_i_vec = vspread_MDF(zero);
114
115         // Switch to vectorized loop for the remaining elements in the particle array
116         
117         for (j = j >>shifter; j < numParticlesByVecSize; j++) {
118
119           // Load the particles' data
120           MD_VEC p0_x_j_vec = p0_x_vec[j];
121           MD_VEC p0_y_j_vec = p0_y_vec[j];
122           MD_VEC p0_z_j_vec = p0_z_vec[j];
123           MD_VEC p0_q_j_vec = p0_q_vec[j];
124
125           // Calculate the vector between the particles
126           MD_VEC p_x_diff_vec = p0_x_i_vec - p0_x_j_vec;
127           MD_VEC p_y_diff_vec = p0_y_i_vec - p0_y_j_vec;
128           MD_VEC p_z_diff_vec = p0_z_i_vec - p0_z_j_vec;
129           // Calculate r and r^2 between the particles
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
133           // Normalize the vector between the particles
134           MD_VEC r_recip_vec = vrecip_MDF(r_vec);
135           MD_VEC r_2_recip_vec = vrecip_MDF(r_2_vec);
136 #if SANITY_CHECK
137           CkAssert(visfinite_MDF(r_2_vec));
138           CkAssert(visfinite_MDF(r_vec));
139 #endif
140 //              CkAssert(&p0_x_vec !=);
141 /*
142         { // how the hell does this not flake out?
143                 CkPrintf("i %d j %d px_i %+6e px_j %+6e py_i %+6e py_j %+6e pz_i %+6e pz_j%+6e\n",i,j,vextract_MDF(p0_x_i_vec,0),vextract_MDF(p0_x_j_vec,0),vextract_MDF(p0_y_i_vec,0),vextract_MDF(p0_y_j_vec,0),vextract_MDF(p0_z_i_vec,0),vextract_MDF(p0_z_j_vec,0));
144         }
145 */
146           MD_VEC p_x_diff_norm_vec = p_x_diff_vec * r_recip_vec;
147           MD_VEC p_y_diff_norm_vec = p_y_diff_vec * r_recip_vec;
148           MD_VEC p_z_diff_norm_vec = p_z_diff_vec * r_recip_vec;
149           MD_VEC f_mag_vec = coulomb_vec * p0_q_i_vec * p0_q_j_vec * r_2_recip_vec;  // Calc force magnitude
150
151           // Multiply the magnitude by the normalized postition difference vector to
152           //   create the force vector
153           MD_VEC f_x_vec = p_x_diff_norm_vec * f_mag_vec;
154           MD_VEC f_y_vec = p_y_diff_norm_vec * f_mag_vec;
155           MD_VEC f_z_vec = p_z_diff_norm_vec * f_mag_vec;
156
157           // Add the force to the outer-loop particle and subtract it from the inner-loop particles
158           f0_x_i_vec += f_x_vec;
159           f0_y_i_vec += f_y_vec;
160           f0_z_i_vec += f_z_vec;
161           f0_x_vec[j] -= f_x_vec;
162           f0_y_vec[j] -= f_y_vec;
163           f0_z_vec[j] -= f_z_vec;
164
165           // DMK - DEBUG
166           #if COUNT_FLOPS != 0
167             localFlopCount += (31 * 4);
168           #endif
169         }
170
171         // Add force values for p0[i] into f0
172         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);
173         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);
174         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);
175       }
176
177       // DMK - DEBUG
178       #if (ENABLE_USER_EVENTS != 0) && defined(CMK_CELL) && (CMK_CELL == 0)
179         double __end_time__ = CkWallTimer();
180         traceUserBracketEvent(PROJ_USER_EVENT_SELFCOMPUTE_DOCALC_WORK, __start_time__, __end_time__);
181       #endif
182
183     } doCalc_callback;
184
185   };
186
187 };