Updated code to refect changes to simd.h. In particular, the removal of the vector...
authorDavid Kunzman <kunzman2@illinois.edu>
Thu, 28 May 2009 23:06:37 +0000 (23:06 +0000)
committerDavid Kunzman <kunzman2@illinois.edu>
Thu, 28 May 2009 23:06:37 +0000 (23:06 +0000)
examples/charm++/cell/md/main.C
examples/charm++/cell/md/main.h
examples/charm++/cell/md/md_config.h
examples/charm++/cell/md/pairCompute.ci
examples/charm++/cell/md/patch.C
examples/charm++/cell/md/patch.ci
examples/charm++/cell/md/patch.h
examples/charm++/cell/md/selfCompute.ci

index d9c6b3b463994d3dc6e5a8807d3872190959afb1..5a99b3631ad1e462c8b1ac11e59c189be3f6dd41 100644 (file)
@@ -83,6 +83,8 @@ Main::Main(CkArgMsg* msg) {
     traceRegisterUserEvent("Patch::integrate_callback()", PROJ_USER_EVENT_PATCH_INTEGRATE_CALLBACK);
     traceRegisterUserEvent("SelfCompute::doCalc_callback()", PROJ_USER_EVENT_SELFCOMPUTE_DOCALC_CALLBACK);
     traceRegisterUserEvent("PairCompute::doCalc_callback()", PROJ_USER_EVENT_PAIRCOMPUTE_DOCALC_CALLBACK);
+    traceRegisterUserEvent("SelfCompute::doCalc() - Work", PROJ_USER_EVENT_SELFCOMPUTE_DOCALC_WORK);
+    traceRegisterUserEvent("PairCompute::doCalc() - Work", PROJ_USER_EVENT_PAIRCOMPUTE_DOCALC_WORK);
     traceRegisterUserEvent("CmiMachineProgressImpl", PROJ_USER_EVENT_MACHINEPROGRESS);
   #endif
 
@@ -105,19 +107,65 @@ Main::Main(CkArgMsg* msg) {
   selfComputeArrayProxy = CProxy_SelfCompute::ckNew(numPatchesX, numPatchesY, numPatchesZ);
 
   // Create the pair compute array
+  #if ENABLE_STATIC_LOAD_BALANCING != 0
+    // NOTE : For now, this code has to be manually changed to match the nodelist file since there is no way to
+    //   pass this information into the program at runtime.  In the future, this is something the runtime system
+    //   take care of in the ideal case.
+    int numPEs = CkNumPes();
+    #define W_X86    ( 10)  //  10
+    #define W_BLADE  (125)  // 100
+    #define W_PS3    ( 96)  //  75
+    // NOTE: The peWeights should match the hetero nodelist file being used
+    //int peWeights[13] = { W_X86, W_BLADE, W_PS3, W_BLADE, W_PS3, W_BLADE, W_PS3, W_BLADE, W_PS3, W_BLADE, W_BLADE, W_BLADE, W_BLADE };
+    int peWeights[13] = { W_BLADE, W_PS3, W_BLADE, W_PS3, W_BLADE, W_PS3, W_BLADE, W_PS3, W_BLADE, W_BLADE, W_BLADE, W_BLADE, W_BLADE };
+    //int peWeights[14] = { W_X86, W_X86, W_BLADE, W_PS3, W_BLADE, W_PS3, W_BLADE, W_PS3, W_BLADE, W_PS3, W_BLADE, W_BLADE, W_BLADE, W_BLADE };
+    int peStats[13] = { 0 };
+    CkAssert(numPEs <= 13);
+    int rValLimit  = 0;
+    for (int i = 0; i < numPEs; i++) { rValLimit += peWeights[i]; }
+  #endif
   pairComputeArrayProxy = CProxy_PairCompute::ckNew();
   const int numPatches = numPatchesX * numPatchesY * numPatchesZ;
   for (int p0 = 0; p0 < numPatches; p0++) {
     for (int p1 = p0 + 1; p1 < numPatches; p1++) {
-      pairComputeArrayProxy(p0, p1).insert();
+      #if ENABLE_STATIC_LOAD_BALANCING != 0
+        int pe = 0;
+        int rVal = rand() % rValLimit;
+        for (int i = 0; i < numPEs; i++) { if (rVal < peWeights[i]) { pe = i; break; } rVal -= peWeights[i]; }
+        pairComputeArrayProxy(p0, p1).insert(pe);
+        peStats[pe]++;
+      #else
+        pairComputeArrayProxy(p0, p1).insert();
+      #endif
     }
   }
   pairComputeArrayProxy.doneInserting();
+  #if ENABLE_STATIC_LOAD_BALANCING != 0
+    int numPairComputes = 0;
+    for (int i = 0; i < numPEs; i++) { numPairComputes += peStats[i]; }
+    for (int i = 0; i < numPEs; i++) {
+      CkPrintf("[STATS] :: peStats[%d] = %6d (%5.2f%%)\n", i, peStats[i], ((float)peStats[i]) / ((float)numPairComputes) * 100.0f);
+    }
+  #endif
 
   // Start initialization (NOTE: Patch will initiate proxy patches directly if proxy patches are being used)
-  patchArrayProxy.init(numParticlesPerPatch);
   selfComputeArrayProxy.init(numParticlesPerPatch);
   pairComputeArrayProxy.init(numParticlesPerPatch);
+
+  #if USE_PROXY_PATCHES != 0
+    for (int x = 0; x < numPatchesX; x++) {
+      for (int y = 0; y < numPatchesY; y++) {
+        for (int z = 0; z < numPatchesZ; z++) {
+          int patchIndex = PATCH_XYZ_TO_I(x, y, z);
+          CProxy_ProxyPatch proxyPatchProxy = CProxy_ProxyPatch::ckNew(patchIndex);
+          proxyPatchProxy.init(numParticlesPerPatch);
+          patchArrayProxy(x, y, z).init(numParticlesPerPatch, proxyPatchProxy);
+       }
+      }
+    }
+  #else
+    patchArrayProxy.init(numParticlesPerPatch);
+  #endif
 }
 
 
@@ -140,7 +188,7 @@ void Main::initCheckIn() {
     numCheckedIn = 0;
 
     // Start timing
-    simStartTime = CkWallTimer();
+    simPrevTime = simStartTime = CkWallTimer();
 
     // One step for main (patches do many steps)
     const int numStepsToDo = (numStepsRemaining > STEPS_PER_PRINT) ? (STEPS_PER_PRINT) : (numStepsRemaining);
@@ -163,11 +211,15 @@ void Main::patchCheckIn() {
 
     // Check to see if there is another step, if so start it, otherwise exit
     numStepsRemaining--;
-    CkPrintf("Main::patchCheckIn() - Simulation (%d steps remaining)...\n", numStepsRemaining);
+    double curTime = CkWallTimer();
+    CkPrintf("Main::patchCheckIn() - Simulation (%d steps remaining)... (deltaTime: %lf sec)...\n",
+             numStepsRemaining, curTime - simPrevTime
+            );
+    simPrevTime = curTime;
     if (numStepsRemaining <= 0) {
 
       // Stop timing and display elapsed time
-      double simStopTime = CkWallTimer();
+      double simStopTime = curTime;
       CkPrintf("Elapsed Time: %lf sec\n", simStopTime - simStartTime);
 
       // DMK - DEBUG
index 913368d9cf26d636fa59a6518f81f95dd99014f8..c2cd002d09bebee5a5229b59a1ce7428b28f7422 100644 (file)
@@ -32,6 +32,7 @@ class Main : public CBase_Main {
     /// Member Variables ///
     int numStepsRemaining;
     double simStartTime;
+    double simPrevTime;
 
     /// Member Functions ///
     void parseCommandLine(int argc, char** argv);
index fe639551631f57a4429a9fcdf8f2332ed2c682d1..ca6521fff8337c5268240c99512eeef87779490d 100644 (file)
 ////////////////////////////////////////////////////////////////////////////////
 // Misc. Macros for Performance Testing
 
+// DMK - DEBUG
+#define ENABLE_STATIC_LOAD_BALANCING     (0)
+
+// DMK - DEBUG
+#define DUMP_INITIAL_PARTICLE_DATA       (0)
+
 // DMK - DEBUG
 #define COUNT_FLOPS                      (0)
 
 // DMK - DEBUG
-#define ENABLE_USER_EVENTS               (1)
+#define ENABLE_USER_EVENTS               (0)
 #define PROJ_USER_EVENT_PATCH_FORCECHECKIN_CALLBACK  (1120)
 #define PROJ_USER_EVENT_PATCH_INTEGRATE_CALLBACK     (1121)
 #define PROJ_USER_EVENT_SELFCOMPUTE_DOCALC_CALLBACK  (1130)
+#define PROJ_USER_EVENT_SELFCOMPUTE_DOCALC_WORK      (1131)
 #define PROJ_USER_EVENT_PAIRCOMPUTE_DOCALC_CALLBACK  (1140)
+#define PROJ_USER_EVENT_PAIRCOMPUTE_DOCALC_WORK      (1141)
 #define PROJ_USER_EVENT_MACHINEPROGRESS              (1150)
 
 // DMK - DEBUG
index 17babbc9c667eadafc24e2eb96f451db8166df01..3b414390955d607682c63e255f00d290da627b97 100644 (file)
@@ -37,6 +37,12 @@ module pairCompute {
                                  writeonly : unsigned int localFlopCount <impl_obj->localFlopCount>
                                ] {
 
+      // DMK - DEBUG
+      #if (ENABLE_USER_EVENTS != 0) && defined(CMK_CELL) && (CMK_CELL == 0)
+        if (numParticles != 256) { CkPrintf("[DEBUG] :: PairCompute::doCalc() - Called with numParticles:%d...\n", numParticles); }
+        double __start_time__ = CmiWallTimer();
+      #endif
+
       // Calculate the electrostatic force (coulumb's law) between the particles
       //   F = (k_e * (q_0 * q_1)) / (r^2)
 
@@ -45,39 +51,40 @@ module pairCompute {
         localFlopCount = 0;
       #endif
 
-      register vec4f* p1_x_vec = (vec4f*)p1_x;
-      register vec4f* p1_y_vec = (vec4f*)p1_y;
-      register vec4f* p1_z_vec = (vec4f*)p1_z;
-      register vec4f* p1_q_vec = (vec4f*)p1_q;
-      register vec4f* f1_x_vec = (vec4f*)f1_x;
-      register vec4f* f1_y_vec = (vec4f*)f1_y;
-      register vec4f* f1_z_vec = (vec4f*)f1_z;
+      register vecf* p1_x_vec = (vecf*)p1_x;
+      register vecf* p1_y_vec = (vecf*)p1_y;
+      register vecf* p1_z_vec = (vecf*)p1_z;
+      register vecf* p1_q_vec = (vecf*)p1_q;
+      register vecf* f1_x_vec = (vecf*)f1_x;
+      register vecf* f1_y_vec = (vecf*)f1_y;
+      register vecf* f1_z_vec = (vecf*)f1_z;
       register int i;
       register int j;
-      register const int numParticlesByVecSize = numParticles / (sizeof(vec4f) / sizeof(float));
+      register const int numParticlesByVecSize = numParticles / vecf_numElems;
 
       // Zero out the force array for p1 (p0's force array will be zero'ed for each
       //   particle in the outer loop)
+      vecf zero_vec = vspreadf(0.0f);
       for (j = 0; j < numParticlesByVecSize; j++) {
-        f1_x_vec[j] = vspread4f(0.0f);
-        f1_y_vec[j] = vspread4f(0.0f);
-        f1_z_vec[j] = vspread4f(0.0f);
+        f1_x_vec[j] = zero_vec;
+        f1_y_vec[j] = zero_vec;
+        f1_z_vec[j] = zero_vec;
       }
 
       // Spread coulumb's constant across a vector
-      vec4f coulomb_vec = vspread4f(COULOMBS_CONSTANT);
+      vecf coulomb_vec = vspreadf(COULOMBS_CONSTANT);
 
       // Outer-loop (p0)
       for (i = 0; i < numParticles; i++) {
 
         // Spread this particle's values out over vectors
-        vec4f p0_x_i_vec = vspread4f(p0_x[i]);
-        vec4f p0_y_i_vec = vspread4f(p0_y[i]);
-        vec4f p0_z_i_vec = vspread4f(p0_z[i]);
-        vec4f p0_q_i_vec = vspread4f(p0_q[i]);
-        vec4f f0_x_i_vec = vspread4f(0.0f);
-        vec4f f0_y_i_vec = vspread4f(0.0f);
-        vec4f f0_z_i_vec = vspread4f(0.0f);
+        vecf p0_x_i_vec = vspreadf(p0_x[i]);
+        vecf p0_y_i_vec = vspreadf(p0_y[i]);
+        vecf p0_z_i_vec = vspreadf(p0_z[i]);
+        vecf p0_q_i_vec = vspreadf(p0_q[i]);
+        vecf f0_x_i_vec = vspreadf(0.0f);
+        vecf f0_y_i_vec = vspreadf(0.0f);
+        vecf f0_z_i_vec = vspreadf(0.0f);
 
         #if 1
 
@@ -93,19 +100,19 @@ module pairCompute {
 
         // Setup Input for Stage 2 - Initially for j+2
         const int jPlus2 = 2;
-        vec4f p1_x_j_vec = p1_x_vec[jPlus2];
-        vec4f p1_y_j_vec = p1_y_vec[jPlus2];
-        vec4f p1_z_j_vec = p1_z_vec[jPlus2];
-        vec4f p1_q_j_vec = p1_q_vec[jPlus2];
+        vecf p1_x_j_vec = p1_x_vec[jPlus2];
+        vecf p1_y_j_vec = p1_y_vec[jPlus2];
+        vecf p1_z_j_vec = p1_z_vec[jPlus2];
+        vecf p1_q_j_vec = p1_q_vec[jPlus2];
 
         // Setup Input for Stage 3 - Initial for j+1
         const int jPlus1 = 1;
-        vec4f p_x_diff_vec = p0_x_i_vec - p1_x_vec[jPlus1];
-        vec4f p_y_diff_vec = p0_y_i_vec - p1_y_vec[jPlus1];
-        vec4f p_z_diff_vec = p0_z_i_vec - p1_z_vec[jPlus1];
-        vec4f p1_q_j_vec_s2 = p1_q_vec[jPlus1];
-        vec4f 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);
-        vec4f r_vec = vsqrt4f(r_2_vec);
+        vecf p_x_diff_vec = p0_x_i_vec - p1_x_vec[jPlus1];
+        vecf p_y_diff_vec = p0_y_i_vec - p1_y_vec[jPlus1];
+        vecf p_z_diff_vec = p0_z_i_vec - p1_z_vec[jPlus1];
+        vecf p1_q_j_vec_s2 = p1_q_vec[jPlus1];
+        vecf 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);
+        vecf r_vec = vsqrtf(r_2_vec);
 
         // DMK - DEBUG
         #if COUNT_FLOPS != 0
@@ -114,21 +121,21 @@ module pairCompute {
 
         // Setup Input for Stage 4 - Initially for j+0
         const int jPlus0 = 0;
-        vec4f f_x_vec;
-        vec4f f_y_vec;
-        vec4f f_z_vec;
+        vecf f_x_vec;
+        vecf f_y_vec;
+        vecf f_z_vec;
         {
-          vec4f p_x_diff_vec = p0_x_i_vec - p1_x_vec[jPlus0];
-          vec4f p_y_diff_vec = p0_y_i_vec - p1_y_vec[jPlus0];
-          vec4f p_z_diff_vec = p0_z_i_vec - p1_z_vec[jPlus0];
-          vec4f 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);
-          vec4f r_vec = vsqrt4f(r_2_vec);
-          vec4f r_recip_vec = vrecip4f(r_vec);
-          vec4f r_2_recip_vec = vrecip4f(r_2_vec);
-          vec4f p_x_diff_norm_vec = p_x_diff_vec * r_recip_vec;
-          vec4f p_y_diff_norm_vec = p_y_diff_vec * r_recip_vec;
-          vec4f p_z_diff_norm_vec = p_z_diff_vec * r_recip_vec;
-          vec4f f_mag_vec = coulomb_vec * ((p0_q_i_vec * p1_q_vec[jPlus0]) * r_2_recip_vec);
+          vecf p_x_diff_vec = p0_x_i_vec - p1_x_vec[jPlus0];
+          vecf p_y_diff_vec = p0_y_i_vec - p1_y_vec[jPlus0];
+          vecf p_z_diff_vec = p0_z_i_vec - p1_z_vec[jPlus0];
+          vecf 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);
+          vecf r_vec = vsqrtf(r_2_vec);
+          vecf r_recip_vec = vrecipf(r_vec);
+          vecf r_2_recip_vec = vrecipf(r_2_vec);
+          vecf p_x_diff_norm_vec = p_x_diff_vec * r_recip_vec;
+          vecf p_y_diff_norm_vec = p_y_diff_vec * r_recip_vec;
+          vecf p_z_diff_norm_vec = p_z_diff_vec * r_recip_vec;
+          vecf f_mag_vec = coulomb_vec * ((p0_q_i_vec * p1_q_vec[jPlus0]) * r_2_recip_vec);
           f_x_vec = p_x_diff_norm_vec * f_mag_vec;
           f_y_vec = p_y_diff_norm_vec * f_mag_vec;
           f_z_vec = p_z_diff_norm_vec * f_mag_vec;
@@ -143,6 +150,44 @@ module pairCompute {
 
         for (j = 0; j < loopUpper; j++) {
 
+          #if 1
+
+          // Stage 4 - Apply force vector
+          f0_x_i_vec = vaddf(f0_x_i_vec, f_x_vec);
+          f0_y_i_vec = vaddf(f0_y_i_vec, f_y_vec);
+          f0_z_i_vec = vaddf(f0_z_i_vec, f_z_vec);
+          f1_x_vec[j] = vsubf(f1_x_vec[j], f_x_vec);
+          f1_y_vec[j] = vsubf(f1_y_vec[j], f_y_vec);
+          f1_z_vec[j] = vsubf(f1_z_vec[j], f_z_vec);
+
+          // Stage 3 - Calculate force vector
+          vecf r_recip_vec = vrecipf(r_vec);
+          vecf r_2_recip_vec = vrecipf(r_2_vec);
+          vecf p_z_diff_norm_vec = vmulf(p_z_diff_vec, r_recip_vec);
+          vecf p_y_diff_norm_vec = vmulf(p_y_diff_vec, r_recip_vec);
+          vecf p_x_diff_norm_vec = vmulf(p_x_diff_vec, r_recip_vec);
+          vecf f_mag_vec = vmulf(vmulf(coulomb_vec, p0_q_i_vec), vmulf(p1_q_j_vec_s2, r_2_recip_vec));
+          f_x_vec = vmulf(p_x_diff_norm_vec, f_mag_vec);
+          f_y_vec = vmulf(p_y_diff_norm_vec, f_mag_vec);
+          f_z_vec = vmulf(p_z_diff_norm_vec, f_mag_vec);
+
+          // Stage 2 - Calculate radius
+          p1_q_j_vec_s2 = p1_q_j_vec;
+          p_x_diff_vec = vsubf(p0_x_i_vec, p1_x_j_vec);
+          p_y_diff_vec = vsubf(p0_y_i_vec, p1_y_j_vec);
+          p_z_diff_vec = vsubf(p0_z_i_vec, p1_z_j_vec);
+          r_2_vec = vmaddf(p_x_diff_vec, p_x_diff_vec, vmaddf(p_y_diff_vec, p_y_diff_vec, vmulf(p_z_diff_vec, p_z_diff_vec)));
+          r_vec = vsqrtf(r_2_vec);
+
+          // Stage 1 - Load particle data
+          const int jIndexToLoad = j + 3;
+          p1_x_j_vec = p1_x_vec[jIndexToLoad];
+          p1_y_j_vec = p1_y_vec[jIndexToLoad];
+          p1_z_j_vec = p1_z_vec[jIndexToLoad];
+          p1_q_j_vec = p1_q_vec[jIndexToLoad];
+
+          #else
+
           // Stage 4 - Apply force vector
           f0_x_i_vec += f_x_vec;
           f0_y_i_vec += f_y_vec;
@@ -152,12 +197,12 @@ module pairCompute {
           f1_z_vec[j] -= f_z_vec;
 
           // Stage 3 - Calculate force vector
-          vec4f r_recip_vec = vrecip4f(r_vec);
-          vec4f r_2_recip_vec = vrecip4f(r_2_vec);
-          vec4f p_x_diff_norm_vec = p_x_diff_vec * r_recip_vec;
-          vec4f p_y_diff_norm_vec = p_y_diff_vec * r_recip_vec;
-          vec4f p_z_diff_norm_vec = p_z_diff_vec * r_recip_vec;
-          vec4f f_mag_vec = coulomb_vec * p0_q_i_vec * p1_q_j_vec_s2 * r_2_recip_vec;  // Calc force magnitude
+          vecf r_recip_vec = vrecipf(r_vec);
+          vecf r_2_recip_vec = vrecipf(r_2_vec);
+          vecf p_x_diff_norm_vec = p_x_diff_vec * r_recip_vec;
+          vecf p_y_diff_norm_vec = p_y_diff_vec * r_recip_vec;
+          vecf p_z_diff_norm_vec = p_z_diff_vec * r_recip_vec;
+          vecf f_mag_vec = coulomb_vec * p0_q_i_vec * p1_q_j_vec_s2 * r_2_recip_vec;  // Calc force magnitude
           f_x_vec = p_x_diff_norm_vec * f_mag_vec;  // Multiply normalized vector by force magnitude
           f_y_vec = p_y_diff_norm_vec * f_mag_vec;
           f_z_vec = p_z_diff_norm_vec * f_mag_vec;
@@ -168,7 +213,7 @@ module pairCompute {
           p_y_diff_vec = p0_y_i_vec - p1_y_j_vec;
           p_z_diff_vec = p0_z_i_vec - p1_z_j_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);
-          r_vec = vsqrt4f(r_2_vec);
+          r_vec = vsqrtf(r_2_vec);
 
           // Stage 1 - Load particle data
           const int jIndexToLoad = j + 3;
@@ -177,6 +222,7 @@ module pairCompute {
           p1_z_j_vec = p1_z_vec[jIndexToLoad];
           p1_q_j_vec = p1_q_vec[jIndexToLoad];
 
+          #endif
 
           // DMK - DEBUG
           #if COUNT_FLOPS != 0
@@ -203,15 +249,15 @@ module pairCompute {
 
         // Complete Calc for Stage 2 Data - Calculate Force Vector, Apply Force Vector
         {
-          vec4f r_recip_vec = vrecip4f(r_vec);
-          vec4f r_2_recip_vec = vrecip4f(r_2_vec);
-          vec4f p_x_diff_norm_vec = p_x_diff_vec * r_recip_vec;
-          vec4f p_y_diff_norm_vec = p_y_diff_vec * r_recip_vec;
-          vec4f p_z_diff_norm_vec = p_z_diff_vec * r_recip_vec;
-          vec4f f_mag_vec = coulomb_vec * ((p0_q_i_vec * p1_q_j_vec_s2) * r_2_recip_vec);
-          vec4f f_x_vec = p_x_diff_norm_vec * f_mag_vec;
-          vec4f f_y_vec = p_y_diff_norm_vec * f_mag_vec;
-          vec4f f_z_vec = p_z_diff_norm_vec * f_mag_vec;
+          vecf r_recip_vec = vrecipf(r_vec);
+          vecf r_2_recip_vec = vrecipf(r_2_vec);
+          vecf p_x_diff_norm_vec = p_x_diff_vec * r_recip_vec;
+          vecf p_y_diff_norm_vec = p_y_diff_vec * r_recip_vec;
+          vecf p_z_diff_norm_vec = p_z_diff_vec * r_recip_vec;
+          vecf f_mag_vec = coulomb_vec * ((p0_q_i_vec * p1_q_j_vec_s2) * r_2_recip_vec);
+          vecf f_x_vec = p_x_diff_norm_vec * f_mag_vec;
+          vecf f_y_vec = p_y_diff_norm_vec * f_mag_vec;
+          vecf f_z_vec = p_z_diff_norm_vec * f_mag_vec;
           f0_x_i_vec += f_x_vec;
           f0_y_i_vec += f_y_vec;
           f0_z_i_vec += f_z_vec;
@@ -227,20 +273,20 @@ module pairCompute {
 
         // Complete Calc for Stage 1 Data - Calculate Radius, Calculate Force Vector, Apply Force Vector
         {
-          vec4f p_x_diff_vec = p0_x_i_vec - p1_x_j_vec;
-          vec4f p_y_diff_vec = p0_y_i_vec - p1_y_j_vec;
-          vec4f p_z_diff_vec = p0_z_i_vec - p1_z_j_vec;
-          vec4f 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);
-          vec4f r_vec = vsqrt4f(r_2_vec);
-          vec4f r_recip_vec = vrecip4f(r_vec);
-          vec4f r_2_recip_vec = vrecip4f(r_2_vec);
-          vec4f p_x_diff_norm_vec = p_x_diff_vec * r_recip_vec;
-          vec4f p_y_diff_norm_vec = p_y_diff_vec * r_recip_vec;
-          vec4f p_z_diff_norm_vec = p_z_diff_vec * r_recip_vec;
-          vec4f f_mag_vec = coulomb_vec * ((p0_q_i_vec * p1_q_j_vec) * r_2_recip_vec);
-          vec4f f_x_vec = p_x_diff_norm_vec * f_mag_vec;
-          vec4f f_y_vec = p_y_diff_norm_vec * f_mag_vec;
-          vec4f f_z_vec = p_z_diff_norm_vec * f_mag_vec;
+          vecf p_x_diff_vec = p0_x_i_vec - p1_x_j_vec;
+          vecf p_y_diff_vec = p0_y_i_vec - p1_y_j_vec;
+          vecf p_z_diff_vec = p0_z_i_vec - p1_z_j_vec;
+          vecf 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);
+          vecf r_vec = vsqrtf(r_2_vec);
+          vecf r_recip_vec = vrecipf(r_vec);
+          vecf r_2_recip_vec = vrecipf(r_2_vec);
+          vecf p_x_diff_norm_vec = p_x_diff_vec * r_recip_vec;
+          vecf p_y_diff_norm_vec = p_y_diff_vec * r_recip_vec;
+          vecf p_z_diff_norm_vec = p_z_diff_vec * r_recip_vec;
+          vecf f_mag_vec = coulomb_vec * ((p0_q_i_vec * p1_q_j_vec) * r_2_recip_vec);
+          vecf f_x_vec = p_x_diff_norm_vec * f_mag_vec;
+          vecf f_y_vec = p_y_diff_norm_vec * f_mag_vec;
+          vecf f_z_vec = p_z_diff_norm_vec * f_mag_vec;
           f0_x_i_vec += f_x_vec;
           f0_y_i_vec += f_y_vec;
           f0_z_i_vec += f_z_vec;
@@ -260,36 +306,36 @@ module pairCompute {
         for (j = 0; j < numParticlesByVecSize; j++) {
 
           // Load the particles' data
-          vec4f p1_x_j_vec = p1_x_vec[j];
-          vec4f p1_y_j_vec = p1_y_vec[j];
-          vec4f p1_z_j_vec = p1_z_vec[j];
-          vec4f p1_q_j_vec = p1_q_vec[j];
+          vecf p1_x_j_vec = p1_x_vec[j];
+          vecf p1_y_j_vec = p1_y_vec[j];
+          vecf p1_z_j_vec = p1_z_vec[j];
+          vecf p1_q_j_vec = p1_q_vec[j];
 
           // Calculate the vector between the particles
-          vec4f p_x_diff_vec = p0_x_i_vec - p1_x_j_vec;
-          vec4f p_y_diff_vec = p0_y_i_vec - p1_y_j_vec;
-          vec4f p_z_diff_vec = p0_z_i_vec - p1_z_j_vec;
+          vecf p_x_diff_vec = p0_x_i_vec - p1_x_j_vec;
+          vecf p_y_diff_vec = p0_y_i_vec - p1_y_j_vec;
+          vecf p_z_diff_vec = p0_z_i_vec - p1_z_j_vec;
 
           // Calculate r and r^2 between the particles
-          vec4f p_x_diff_2_vec = p_x_diff_vec * p_x_diff_vec;
-          vec4f p_y_diff_2_vec = p_y_diff_vec * p_y_diff_vec;
-          vec4f p_z_diff_2_vec = p_z_diff_vec * p_z_diff_vec;
-          vec4f r_2_vec = p_x_diff_2_vec + p_y_diff_2_vec + p_z_diff_2_vec;
-          vec4f r_vec = vsqrt4f(r_2_vec);
+          vecf p_x_diff_2_vec = p_x_diff_vec * p_x_diff_vec;
+          vecf p_y_diff_2_vec = p_y_diff_vec * p_y_diff_vec;
+          vecf p_z_diff_2_vec = p_z_diff_vec * p_z_diff_vec;
+          vecf r_2_vec = p_x_diff_2_vec + p_y_diff_2_vec + p_z_diff_2_vec;
+          vecf r_vec = vsqrtf(r_2_vec);
 
           // Normalize the vector between the particles
-          vec4f p_x_diff_norm_vec = p_x_diff_vec / r_vec;
-          vec4f p_y_diff_norm_vec = p_y_diff_vec / r_vec;
-          vec4f p_z_diff_norm_vec = p_z_diff_vec / r_vec;
+          vecf p_x_diff_norm_vec = p_x_diff_vec / r_vec;
+          vecf p_y_diff_norm_vec = p_y_diff_vec / r_vec;
+          vecf p_z_diff_norm_vec = p_z_diff_vec / r_vec;
 
           // Calculate the magnitude of the electrostatic force between the particles
-          vec4f f_mag_vec = coulomb_vec * ((p0_q_i_vec * p1_q_j_vec) / r_2_vec);
+          vecf f_mag_vec = coulomb_vec * ((p0_q_i_vec * p1_q_j_vec) / r_2_vec);
 
           // Multiply the magnitude by the normalized postition difference vector to
           //   create the force vector
-          vec4f f_x_vec = p_x_diff_norm_vec * f_mag_vec;
-          vec4f f_y_vec = p_y_diff_norm_vec * f_mag_vec;
-          vec4f f_z_vec = p_z_diff_norm_vec * f_mag_vec;
+          vecf f_x_vec = p_x_diff_norm_vec * f_mag_vec;
+          vecf f_y_vec = p_y_diff_norm_vec * f_mag_vec;
+          vecf f_z_vec = p_z_diff_norm_vec * f_mag_vec;
 
           // Add the force to the outer-loop particle and subtract it from the inner-loop particles
           f0_x_i_vec += f_x_vec;
@@ -309,11 +355,27 @@ module pairCompute {
         #endif
 
         // Write force values for the outer-loop particle back to memory
-        f0_x[i] += vextract4f(f0_x_i_vec, 0) + vextract4f(f0_x_i_vec, 1) + vextract4f(f0_x_i_vec, 2) + vextract4f(f0_x_i_vec, 3);
-        f0_y[i] += vextract4f(f0_y_i_vec, 0) + vextract4f(f0_y_i_vec, 1) + vextract4f(f0_y_i_vec, 2) + vextract4f(f0_y_i_vec, 3);
-        f0_z[i] += vextract4f(f0_z_i_vec, 0) + vextract4f(f0_z_i_vec, 1) + vextract4f(f0_z_i_vec, 2) + vextract4f(f0_z_i_vec, 3);
+        f0_x[i] = vextractf(f0_x_i_vec, 0) + vextractf(f0_x_i_vec, 1) + vextractf(f0_x_i_vec, 2) + vextractf(f0_x_i_vec, 3);
+        f0_y[i] = vextractf(f0_y_i_vec, 0) + vextractf(f0_y_i_vec, 1) + vextractf(f0_y_i_vec, 2) + vextractf(f0_y_i_vec, 3);
+        f0_z[i] = vextractf(f0_z_i_vec, 0) + vextractf(f0_z_i_vec, 1) + vextractf(f0_z_i_vec, 2) + vextractf(f0_z_i_vec, 3);
       }
 
+      // DMK - DEBUG
+      #if (ENABLE_USER_EVENTS != 0) && defined(CMK_CELL) && (CMK_CELL == 0)
+        double __end_time__ = CmiWallTimer();
+        double deltaTime = __end_time__ - __start_time__;
+        traceUserBracketEvent(PROJ_USER_EVENT_PAIRCOMPUTE_DOCALC_WORK, __start_time__, __end_time__);
+        if (deltaTime > 0.002) {
+          CkPrintf("[DEBUG] :: PairCompute[%d,%d]::doCalc() - Took %lf ms to execute... numParticles:%d\n"
+                   "             p0[0] = { %e, %e, %e : %e }, p1[0] = { %e, %e, %e : %e }...\n"
+                   "             f0[0] = { %e, %e, %e }, f1[0] = { %e, %e, %e }...\n",
+                   thisIndex_x, thisIndex_y, deltaTime * 1000.0f, numParticles,
+                   p0_x[0], p0_y[0], p0_z[0], p0_q[0], p1_x[0], p1_y[0], p1_z[0], p1_q[0],
+                   f0_x[0], f0_y[0], f0_z[0], f1_x[0], f1_y[0], f1_z[0]
+                  );
+        }
+      #endif
+
     } doCalc_callback;
 
   };
index 8235234ef5bbf097ce7d89ba39036a5098213a9b..195d7a8ecd25dfd17fae6be9c6d774ebb1ffa2a6 100644 (file)
@@ -16,17 +16,27 @@ void Patch::randomizeParticles() {
     velocityX[i] = 0.0f;
     velocityY[i] = 0.0f;
     velocityZ[i] = 0.0f;
+
+    // DMK - DEBUG
+    #if DUMP_INITIAL_PARTICLE_DATA != 0
+      CkPrintf("[INFO] :: Patch[%02d,%02d,%02d]::randomizeParticles() - particle[%04d] = { "
+               "px:%+6e, py:%+6e, pz:%+6e, q:%+6e, m:%+6e, vx:%+6e, vy:%+6e, vz:%+6e }\n",
+               thisIndex.x, thisIndex.y, thisIndex.z, i,
+               particleX[i], particleY[i], particleZ[i], particleQ[i], particleM[i],
+               velocityX[i], velocityY[i], velocityZ[i]
+              );
+    #endif
   }
 }
 
 
 float Patch::randf() {
-  return (((float)(rand() % 1483727)) / (1483727.0f));
+  const int mask = 0x7FFFFFFF;
+  return (((float)(rand() % mask)) / ((float)(mask)));
 }
 
 
 Patch::Patch() {
-
   particleX = particleY = particleZ = NULL;
   particleQ = particleM = NULL;
   forceSumX = forceSumY = forceSumZ = NULL;
@@ -57,15 +67,17 @@ Patch::~Patch() {
 
 
 void Patch::init(int numParticles) {
+  init_common(numParticles);
+}
 
-  // Also tell the proxy patches to initialize
+void Patch::init(int numParticles, CProxy_ProxyPatch proxyPatchProxy) {
   #if USE_PROXY_PATCHES != 0
-  {
-    const int patchIndex = PATCH_XYZ_TO_I(thisIndex.x, thisIndex.y, thisIndex.z);
-    proxyPatchProxy = CProxy_ProxyPatch::ckNew(patchIndex);
-    proxyPatchProxy.init(numParticles);
-  }
+    this->proxyPatchProxy = proxyPatchProxy;
   #endif
+  init_common(numParticles);
+}
+
+void Patch::init_common(int numParticles) {
 
   // Allocate memory for the particles
   this->numParticles = numParticles;
@@ -98,6 +110,7 @@ void Patch::init(int numParticles) {
     }
     pairComputeArraySection_lower = CProxySection_PairCompute::ckNew(pairComputeArrayProxy, pairIndexes_lower.getVec(), pairIndexes_lower.size());
 
+    // Upper section
     CkVec<CkArrayIndex2D> pairIndexes_upper;
     for (int i = patchIndex + 1; i < numPatches; i++) {
       pairIndexes_upper.push_back(CkArrayIndex2D(patchIndex, i));
@@ -129,11 +142,11 @@ void Patch::startIteration_common(int numIters) {
     memset(forceSumY, 0, sizeof(float) * numParticles);
     memset(forceSumZ, 0, sizeof(float) * numParticles);
   #else
-    register vec4f* fsx = (vec4f*)forceSumX;
-    register vec4f* fsy = (vec4f*)forceSumY;
-    register vec4f* fsz = (vec4f*)forceSumZ;
-    const vec4f zero_vec = vspread4f(0.0f);
-    register const int numParticles_vec = numParticles / (sizeof(vec4f) * sizeof(float));
+    register vecf* fsx = (vecf*)forceSumX;
+    register vecf* fsy = (vecf*)forceSumY;
+    register vecf* fsz = (vecf*)forceSumZ;
+    const vecf zero_vec = vspreadf(0.0f);
+    register const int numParticles_vec = numParticles / vecf_numElems;
     for (int i = 0; i < numParticles_vec; i++) {
       fsx[i] = zero_vec;
       fsy[i] = zero_vec;
@@ -194,18 +207,18 @@ void Patch::forceCheckIn(int numParticles, float* forceX, float* forceY, float*
 
   // Accumulate the force data
   #if 0
-    register vec4f* fsx = (vec4f*)forceSumX;
-    register vec4f* fsy = (vec4f*)forceSumY;
-    register vec4f* fsz = (vec4f*)forceSumZ;
-    register vec4f* fx = (vec4f*)forceX;
-    register vec4f* fy = (vec4f*)forceY;
-    register vec4f* fz = (vec4f*)forceZ;
-    register const int numParticles_vec = numParticles / (sizeof(vec4f) * sizeof(float));
+    register vecf* fsx = (vecf*)forceSumX;
+    register vecf* fsy = (vecf*)forceSumY;
+    register vecf* fsz = (vecf*)forceSumZ;
+    register vecf* fx = (vecf*)forceX;
+    register vecf* fy = (vecf*)forceY;
+    register vecf* fz = (vecf*)forceZ;
+    register const int numParticles_vec = numParticles / vecf_numElems;
     register int i;
     for (i = 0; i < numParticles_vec; i++) {
-      fsx[i] = vadd4f(fsx[i], fx[i]);
-      fsy[i] = vadd4f(fsy[i], fy[i]);
-      fsz[i] = vadd4f(fsz[i], fz[i]);
+      fsx[i] = vaddf(fsx[i], fx[i]);
+      fsy[i] = vaddf(fsy[i], fy[i]);
+      fsz[i] = vaddf(fsz[i], fz[i]);
     }
   #else
     for (int i = 0; i < numParticles; i++) {
@@ -352,13 +365,13 @@ void ProxyPatch::forceCheckIn(int numParticles, float* forceX, float* forceY, fl
 
   // Accumulate the force data
   #if USE_PROXY_PATCHES != 0  // Calls will be local and pointers will be aligned, so take advantage and vectorize the code
-    register vec4f* forceX_vec = (vec4f*)forceX;
-    register vec4f* forceY_vec = (vec4f*)forceY;
-    register vec4f* forceZ_vec = (vec4f*)forceZ;
-    register vec4f* forceSumX_vec = (vec4f*)forceSumX;
-    register vec4f* forceSumY_vec = (vec4f*)forceSumY;
-    register vec4f* forceSumZ_vec = (vec4f*)forceSumZ;
-    const int numParticles_vec = numParticles / (sizeof(vec4f) / sizeof(float));
+    register vecf* forceX_vec = (vecf*)forceX;
+    register vecf* forceY_vec = (vecf*)forceY;
+    register vecf* forceZ_vec = (vecf*)forceZ;
+    register vecf* forceSumX_vec = (vecf*)forceSumX;
+    register vecf* forceSumY_vec = (vecf*)forceSumY;
+    register vecf* forceSumZ_vec = (vecf*)forceSumZ;
+    const int numParticles_vec = numParticles / vecf_numElems;
     for (int i = 0; i < numParticles_vec; i++) {
       forceSumX_vec[i] += forceX_vec[i];
       forceSumY_vec[i] += forceY_vec[i];
index 032f77bc5da0742b6fbd49c3edaafe967513d369..5b23bf45c40db73703b737643a3be8ccf79f3a48 100644 (file)
@@ -8,27 +8,48 @@ module patch {
                     float* forceSumX, float* forceSumY, float* forceSumZ,
                     float* forceX, float* forceY, float* forceZ
                    ) {
-      register vec4f* fsx = (vec4f*)forceSumX;
-      register vec4f* fsy = (vec4f*)forceSumY;
-      register vec4f* fsz = (vec4f*)forceSumZ;
-      register vec4f* fx = (vec4f*)forceX;
-      register vec4f* fy = (vec4f*)forceY;
-      register vec4f* fz = (vec4f*)forceZ;
-      register const int numParticles_vec = numParticles / (sizeof(vec4f) * sizeof(float));
+      register vecf* fsx = (vecf*)forceSumX;
+      register vecf* fsy = (vecf*)forceSumY;
+      register vecf* fsz = (vecf*)forceSumZ;
+      register vecf* fx = (vecf*)forceX;
+      register vecf* fy = (vecf*)forceY;
+      register vecf* fz = (vecf*)forceZ;
+      register const int numParticles_vec = numParticles / vecf_numElems;
       register int i;
       for (i = 0; i < numParticles_vec; i++) {
-        fsx[i] = vadd4f(fsx[i], fx[i]);
-        fsy[i] = vadd4f(fsy[i], fy[i]);
-        fsz[i] = vadd4f(fsz[i], fz[i]);
+        fsx[i] = vaddf(fsx[i], fx[i]);
+        fsy[i] = vaddf(fsy[i], fy[i]);
+        fsz[i] = vaddf(fsz[i], fz[i]);
       }
     }
   };
 
+  group ProxyPatch {
+
+    entry ProxyPatch(int patchIndex);
+
+    entry void init(int numParticles);
+
+    entry void patchData(int numParticles,
+                         float particleX[numParticles],
+                         float particleY[numParticles],
+                         float particleZ[numParticles],
+                         float particleQ[numParticles]
+                        );
+
+    entry void forceCheckIn(int numParticles,
+                            float forceX[numParticles],
+                            float forceY[numParticles],
+                            float forceZ[numParticles]
+                           );
+  };
+
   array[3D] Patch {
 
     entry Patch();
 
     entry void init(int numParticles);
+    entry void init(int numParticles, CProxy_ProxyPatch proxyPatchProxy);
 
     entry void startIteration();
     entry void startIterations(int numIters);
@@ -67,28 +88,28 @@ module patch {
         localFlopCount = 0;
       #endif
 
-      register vec4f* p_x_vec = (vec4f*)p_x;
-      register vec4f* p_y_vec = (vec4f*)p_y;
-      register vec4f* p_z_vec = (vec4f*)p_z;
-      register vec4f* p_m_vec = (vec4f*)p_m;
-      register vec4f* f_x_vec = (vec4f*)f_x;
-      register vec4f* f_y_vec = (vec4f*)f_y;
-      register vec4f* f_z_vec = (vec4f*)f_z;
-      register vec4f* v_x_vec = (vec4f*)v_x;
-      register vec4f* v_y_vec = (vec4f*)v_y;
-      register vec4f* v_z_vec = (vec4f*)v_z;
+      register vecf* p_x_vec = (vecf*)p_x;
+      register vecf* p_y_vec = (vecf*)p_y;
+      register vecf* p_z_vec = (vecf*)p_z;
+      register vecf* p_m_vec = (vecf*)p_m;
+      register vecf* f_x_vec = (vecf*)f_x;
+      register vecf* f_y_vec = (vecf*)f_y;
+      register vecf* f_z_vec = (vecf*)f_z;
+      register vecf* v_x_vec = (vecf*)v_x;
+      register vecf* v_y_vec = (vecf*)v_y;
+      register vecf* v_z_vec = (vecf*)v_z;
       register int i;
-      register const int numParticles_vec = numParticles / (sizeof(vec4f) / sizeof(float));
+      register const int numParticles_vec = numParticles / vecf_numElems;
 
-      vec4f delta_time_vec = vspread4f(TIME_PER_STEP);
+      vecf delta_time_vec = vspreadf(TIME_PER_STEP);
 
       // Add the acceleration to the velocity and add the update velocity to the position
       for (i = 0; i < numParticles_vec; i++) {
 
-        // Update the velocity : v_new = v_old + (F/m)
-        vec4f new_v_x_vec = v_x_vec[i] + (delta_time_vec * (f_x_vec[i] / p_m_vec[i]));
-        vec4f new_v_y_vec = v_y_vec[i] + (delta_time_vec * (f_y_vec[i] / p_m_vec[i]));
-        vec4f new_v_z_vec = v_z_vec[i] + (delta_time_vec * (f_z_vec[i] / p_m_vec[i]));
+        // Update the velocity : v_new = v_old + (F/m)dt
+        vecf new_v_x_vec = v_x_vec[i] + (delta_time_vec * (f_x_vec[i] / p_m_vec[i]));
+        vecf new_v_y_vec = v_y_vec[i] + (delta_time_vec * (f_y_vec[i] / p_m_vec[i]));
+        vecf new_v_z_vec = v_z_vec[i] + (delta_time_vec * (f_z_vec[i] / p_m_vec[i]));
         v_x_vec[i] = new_v_x_vec;
         v_y_vec[i] = new_v_y_vec;
         v_z_vec[i] = new_v_z_vec;
@@ -108,25 +129,4 @@ module patch {
 
   };
 
-
-  group ProxyPatch {
-
-    entry ProxyPatch(int patchIndex);
-
-    entry void init(int numParticles);
-
-    entry void patchData(int numParticles,
-                         float particleX[numParticles],
-                         float particleY[numParticles],
-                         float particleZ[numParticles],
-                         float particleQ[numParticles]
-                        );
-
-    entry void forceCheckIn(int numParticles,
-                            float forceX[numParticles],
-                            float forceY[numParticles],
-                            float forceZ[numParticles]
-                           );
-  };
-
 };
index 5cdc8e77a75d38eab4ec04d957799f2318b473c6..eb8ed11028bc30397ffc3e31f001fba78b846958 100644 (file)
@@ -61,6 +61,8 @@ class Patch : public CBase_Patch {
 
     /// Entry Methods ///
     void init(int numParticles);
+    void init(int numParticles, CProxy_ProxyPatch proxyPatchProxy);
+    void init_common(int numParticles);
     void startIteration();
     void startIterations(int numIters);
     void forceCheckIn(int numParticles, float* forceX, float* forceY, float* forceZ);
index 3896f3bcb38d00a080d43f1c5a8e3249484cd752..e0c30d60f2195ee5564e51ece4131ea1f0626bff 100644 (file)
@@ -30,6 +30,11 @@ module selfCompute {
                                 writeonly : unsigned int localFlopCount <impl_obj->localFlopCount>
                                ] {
 
+      // DMK - DEBUG
+      #if (ENABLE_USER_EVENTS != 0) && defined(CMK_CELL) && (CMK_CELL == 0)
+        double __start_time__ = CmiWallTimer();
+      #endif
+
       // Calculate the electrostatic force (coulumb's law) between the particles
       //   F = (k_e * (q_0 * q_1)) / (r^2)
 
@@ -38,26 +43,27 @@ module selfCompute {
         localFlopCount = 0;
       #endif
 
-      register vec4f* p0_x_vec = (vec4f*)particleX;
-      register vec4f* p0_y_vec = (vec4f*)particleY;
-      register vec4f* p0_z_vec = (vec4f*)particleZ;
-      register vec4f* p0_q_vec = (vec4f*)particleQ;
-      register vec4f* f0_x_vec = (vec4f*)f0_x;
-      register vec4f* f0_y_vec = (vec4f*)f0_y;
-      register vec4f* f0_z_vec = (vec4f*)f0_z;
+      register vecf* p0_x_vec = (vecf*)particleX;
+      register vecf* p0_y_vec = (vecf*)particleY;
+      register vecf* p0_z_vec = (vecf*)particleZ;
+      register vecf* p0_q_vec = (vecf*)particleQ;
+      register vecf* f0_x_vec = (vecf*)f0_x;
+      register vecf* f0_y_vec = (vecf*)f0_y;
+      register vecf* f0_z_vec = (vecf*)f0_z;
       register int i;
       register int j;
-      register const int numParticlesByVecSize = numParticles / (sizeof(vec4f) / sizeof(float));
+      register const int numParticlesByVecSize = numParticles / vecf_numElems;
 
       // Zero out the force output
+      vecf zero_vec = vspreadf(0.0f);
       for (j = 0; j < numParticlesByVecSize; j++) {
-        f0_x_vec[j] = vspread4f(0.0f);
-        f0_y_vec[j] = vspread4f(0.0f);
-        f0_z_vec[j] = vspread4f(0.0f);
+        f0_x_vec[j] = zero_vec;
+        f0_y_vec[j] = zero_vec;
+        f0_z_vec[j] = zero_vec;
       }
 
       // Spread coulumb's constant across a vector
-      vec4f coulomb_vec = vspread4f(COULOMBS_CONSTANT);
+      vecf coulomb_vec = vspreadf(COULOMBS_CONSTANT);
 
       // Outer-loop
       for (i = 0; i < numParticles; i++) {
@@ -96,45 +102,45 @@ module selfCompute {
         }
 
         // Spread this particle's (p0[i]) values out over vectors
-        vec4f p0_x_i_vec = vspread4f(particleX[i]);
-        vec4f p0_y_i_vec = vspread4f(particleY[i]);
-        vec4f p0_z_i_vec = vspread4f(particleZ[i]);
-        vec4f p0_q_i_vec = vspread4f(particleQ[i]);
-        vec4f f0_x_i_vec = vspread4f(0.0f);
-        vec4f f0_y_i_vec = vspread4f(0.0f);
-        vec4f f0_z_i_vec = vspread4f(0.0f);
+        vecf p0_x_i_vec = vspreadf(particleX[i]);
+        vecf p0_y_i_vec = vspreadf(particleY[i]);
+        vecf p0_z_i_vec = vspreadf(particleZ[i]);
+        vecf p0_q_i_vec = vspreadf(particleQ[i]);
+        vecf f0_x_i_vec = vspreadf(0.0f);
+        vecf f0_y_i_vec = vspreadf(0.0f);
+        vecf f0_z_i_vec = vspreadf(0.0f);
 
         // Switch to vectorized loop for the remaining elements in the particle array
         for (j = j >> 2; j < numParticlesByVecSize; j++) {
 
           // Load the particles' data
-          vec4f p0_x_j_vec = p0_x_vec[j];
-          vec4f p0_y_j_vec = p0_y_vec[j];
-          vec4f p0_z_j_vec = p0_z_vec[j];
-          vec4f p0_q_j_vec = p0_q_vec[j];
+          vecf p0_x_j_vec = p0_x_vec[j];
+          vecf p0_y_j_vec = p0_y_vec[j];
+          vecf p0_z_j_vec = p0_z_vec[j];
+          vecf p0_q_j_vec = p0_q_vec[j];
 
           // Calculate the vector between the particles
-          vec4f p_x_diff_vec = p0_x_i_vec - p0_x_j_vec;
-          vec4f p_y_diff_vec = p0_y_i_vec - p0_y_j_vec;
-          vec4f p_z_diff_vec = p0_z_i_vec - p0_z_j_vec;
+          vecf p_x_diff_vec = p0_x_i_vec - p0_x_j_vec;
+          vecf p_y_diff_vec = p0_y_i_vec - p0_y_j_vec;
+          vecf p_z_diff_vec = p0_z_i_vec - p0_z_j_vec;
 
           // Calculate r and r^2 between the particles
-          vec4f 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);
-          vec4f r_vec = vsqrt4f(r_2_vec);
+          vecf 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);
+          vecf r_vec = vsqrtf(r_2_vec);
 
           // Normalize the vector between the particles
-          vec4f r_recip_vec = vrecip4f(r_vec);
-          vec4f r_2_recip_vec = vrecip4f(r_2_vec);
-          vec4f p_x_diff_norm_vec = p_x_diff_vec * r_recip_vec;
-          vec4f p_y_diff_norm_vec = p_y_diff_vec * r_recip_vec;
-          vec4f p_z_diff_norm_vec = p_z_diff_vec * r_recip_vec;
-          vec4f f_mag_vec = coulomb_vec * p0_q_i_vec * p0_q_j_vec * r_2_recip_vec;  // Calc force magnitude
+          vecf r_recip_vec = vrecipf(r_vec);
+          vecf r_2_recip_vec = vrecipf(r_2_vec);
+          vecf p_x_diff_norm_vec = p_x_diff_vec * r_recip_vec;
+          vecf p_y_diff_norm_vec = p_y_diff_vec * r_recip_vec;
+          vecf p_z_diff_norm_vec = p_z_diff_vec * r_recip_vec;
+          vecf f_mag_vec = coulomb_vec * p0_q_i_vec * p0_q_j_vec * r_2_recip_vec;  // Calc force magnitude
 
           // Multiply the magnitude by the normalized postition difference vector to
           //   create the force vector
-          vec4f f_x_vec = p_x_diff_norm_vec * f_mag_vec;
-          vec4f f_y_vec = p_y_diff_norm_vec * f_mag_vec;
-          vec4f f_z_vec = p_z_diff_norm_vec * f_mag_vec;
+          vecf f_x_vec = p_x_diff_norm_vec * f_mag_vec;
+          vecf f_y_vec = p_y_diff_norm_vec * f_mag_vec;
+          vecf f_z_vec = p_z_diff_norm_vec * f_mag_vec;
 
           // Add the force to the outer-loop particle and subtract it from the inner-loop particles
           f0_x_i_vec += f_x_vec;
@@ -151,11 +157,17 @@ module selfCompute {
         }
 
         // Add force values for p0[i] into f0
-        f0_x[i] += vextract4f(f0_x_i_vec, 0) + vextract4f(f0_x_i_vec, 1) + vextract4f(f0_x_i_vec, 2) + vextract4f(f0_x_i_vec, 3);
-        f0_y[i] += vextract4f(f0_y_i_vec, 0) + vextract4f(f0_y_i_vec, 1) + vextract4f(f0_y_i_vec, 2) + vextract4f(f0_y_i_vec, 3);
-        f0_z[i] += vextract4f(f0_z_i_vec, 0) + vextract4f(f0_z_i_vec, 1) + vextract4f(f0_z_i_vec, 2) + vextract4f(f0_z_i_vec, 3);
+        f0_x[i] += vextractf(f0_x_i_vec, 0) + vextractf(f0_x_i_vec, 1) + vextractf(f0_x_i_vec, 2) + vextractf(f0_x_i_vec, 3);
+        f0_y[i] += vextractf(f0_y_i_vec, 0) + vextractf(f0_y_i_vec, 1) + vextractf(f0_y_i_vec, 2) + vextractf(f0_y_i_vec, 3);
+        f0_z[i] += vextractf(f0_z_i_vec, 0) + vextractf(f0_z_i_vec, 1) + vextractf(f0_z_i_vec, 2) + vextractf(f0_z_i_vec, 3);
       }
 
+      // DMK - DEBUG
+      #if (ENABLE_USER_EVENTS != 0) && defined(CMK_CELL) && (CMK_CELL == 0)
+        double __end_time__ = CmiWallTimer();
+        traceUserBracketEvent(PROJ_USER_EVENT_SELFCOMPUTE_DOCALC_WORK, __start_time__, __end_time__);
+      #endif
+
     } doCalc_callback;
 
   };