Merge profiling from devel branch 10/5210/2 master
authorDavid <dhardy@ks.uiuc.edu>
Wed, 3 Jul 2019 22:52:49 +0000 (17:52 -0500)
committerDavid <dhardy@ks.uiuc.edu>
Mon, 8 Jul 2019 22:42:45 +0000 (17:42 -0500)
Includes NVIDIA NVTX, Charm++ Projections tracing user events, and
detailed in-house timer collections with optional histogramming.
The profiling here is toward understanding the CPU performance of the
numerical integration parts of the code in Sequencer and HomePatch.

To build with NVIDIA NVTX:
- define NAMD_NVTX_ENABLED (compiler flag -DNAMD_NVTX_ENABLED)
- link to libnvToolsExt (Make.config use "EXTRALINKLIBS = -lnvToolsExt")

The Charm++ tracing user events requires verion 6.9.0 of Charm++.
To build with Charm++ Projections tracing user events:
- build Charm++ with "--with-production --enable-tracing"
- define CMK_TRACE_ENABLED (compiler flag -DCMK_TRACE_ENABLED)
- define NAMD_CMK_TRACE_ENALBED (compiler flag -DNAMD_CMK_TRACE_ENABLED)
- build NAMD using "projections" target, e.g. "make projections"

To build with timer collections:
- define TIMER_COLLECTION (compiler flag -DTIMER_COLLECTION)
- optionally define TIMER_HISTOGRAM (compiler flag -DTIMER_HISTOGRAM)

Change-Id: Ia77ba9d4a6ca4ae81deb65103a8ee67ba74d485a

src/Debug.h
src/HomePatch.C
src/HomePatch.h
src/NamdEventsProfiling.def [new file with mode: 0644]
src/NamdEventsProfiling.h [new file with mode: 0644]
src/PatchTypes.h
src/Sequencer.C
src/SimParameters.C
src/SimParameters.h

index 1216cb5..ce970ae 100644 (file)
@@ -120,5 +120,41 @@ public:
 };
 #endif
 
+
+#undef D_MSG
+#undef D_STR
+#undef D_INT
+#undef D_REAL
+#undef D_REXP
+#undef D_VEC
+
+#ifdef NL_DEBUG
+
+#define D_MSG(t) \
+  fprintf(stderr, "DEBUG: %s:%d: %s\n", __FILE__, __LINE__, t)
+#define D_STR(t) \
+  fprintf(stderr, "DEBUG: %s:%d: %s=\"%s\"\n", __FILE__, __LINE__, #t, t)
+#define D_INT(n) \
+  fprintf(stderr, "DEBUG: %s:%d: %s=%d\n", __FILE__, __LINE__, #n, n)
+#define D_REAL(r) \
+  fprintf(stderr, "DEBUG: %s:%d: %s=%g\n", __FILE__, __LINE__, #r, double(r))
+#define D_REXP(r) \
+  fprintf(stderr, "DEBUG: %s:%d: %s=%e\n", __FILE__, __LINE__, #r, double(r))
+#define D_VEC(v) \
+  fprintf(stderr, "DEBUG: %s:%d: %s=%g %g %g\n", __FILE__, __LINE__, \
+      #v, double(v.x), double(v.y), double(v.z))
+
+#else
+
+#define D_MSG(t)
+#define D_STR(t)
+#define D_INT(t)
+#define D_REAL(t)
+#define D_REXP(t)
+#define D_VEC(t)
+
+#endif // NL_DEBUG
+
+
 #endif /* DEBUG_H */
 
index f794d6e..86dcf9a 100644 (file)
 #include "ComputeQM.h"
 #include "ComputeQMMgr.decl.h"
 
+#include "NamdEventsProfiling.h"
+
 //#define PRINT_COMP
 #define TINY 1.0e-20;
 #define MAXHGS 10
 #define MIN_DEBUG_LEVEL 2
 //#define DEBUGM
+//#define NL_DEBUG
 #include "Debug.h"
 
 #include <vector>
@@ -88,6 +91,20 @@ void mollify(CompAtom *qtilde,const HGArrayVector &q0,const BigReal *lambda, HGA
 
 #endif
 
+#ifdef TIMER_COLLECTION
+const char *TimerSet::tlabel[TimerSet::NUMTIMERS] = {
+  "kick",
+  "maxmove",
+  "drift",
+  "piston",
+  "submithalf",
+  "velbbk1",
+  "velbbk2",
+  "rattle1",
+  "submitfull",
+  "submitcollect",
+};
+#endif
 
 HomePatch::HomePatch(PatchID pd, FullAtomList &al) : Patch(pd)
 // DMK - Atom Separation (water vs. non-water)
@@ -921,7 +938,13 @@ void HomePatch::positionsReady(int doMigration)
       doMarginCheck();
     }
   }
-  
+
+#if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED)
+  char prbuf[32];
+  sprintf(prbuf, "%s: %d", NamdProfileEventStr[NamdProfileEvent::POSITIONS_READY], this->getPatchID());
+  NAMD_EVENT_START_EX(1, NamdProfileEvent::POSITIONS_READY, prbuf);
+#endif
+
   if (doMigration && simParams->qmLSSOn)
       qmSwapAtoms();
 
@@ -1006,13 +1029,26 @@ void HomePatch::positionsReady(int doMigration)
   doMigration = doMigration || doAtomUpdate;
   doAtomUpdate = false;
 
-  // Workaround for oversize groups
+  // Workaround for oversize groups:
+  // reset nonbondedGroupSize (ngs) before force calculation,
+  // making sure that subset of hydrogen group starting with 
+  // parent atom are all within 0.5 * hgroupCutoff.
+  // XXX hydrogentGroupSize remains constant but is checked for nonzero
+  // XXX should be skipped for CUDA, ngs not used by CUDA kernels
+  // XXX should this also be skipped for KNL kernels?
+  // ngs used by ComputeNonbondedBase.h - CPU nonbonded kernels
+  // ngs used by ComputeGBIS.C - CPU GB nonbonded kernels
+#if ! defined(NAMD_CUDA)
   doGroupSizeCheck();
+#endif
 
   // Copy information needed by computes and proxys to Patch::p.
-  p.resize(numAtoms);
+  // Resize only if atoms were migrated
+  if (doMigration) {
+    p.resize(numAtoms);
+    pExt.resize(numAtoms);
+  }
   CompAtom *p_i = p.begin();
-  pExt.resize(numAtoms);
   CompAtomExt *pExt_i = pExt.begin();
   FullAtom *a_i = atom.begin();
   int i; int n = numAtoms;
@@ -1254,6 +1290,9 @@ void HomePatch::positionsReady(int doMigration)
 
   // gzheng
   Sync::Object()->PatchReady();
+
+  NAMD_EVENT_STOP(1, NamdProfileEvent::POSITIONS_READY);
+
 }
 
 void HomePatch::replaceForces(ExtForce *f)
@@ -2340,8 +2379,9 @@ int HomePatch::rattle1(const BigReal timestep, Tensor *virial,
     for (int i = 0; i < hgs; ++i ) {
       ref[i] = atom[ig+i].position;
       pos[i] = atom[ig+i].position;
-      if (!(fixedAtomsOn && atom[ig+i].atomFixed))
+      if (!(fixedAtomsOn && atom[ig+i].atomFixed)) {
         pos[i] += atom[ig+i].velocity * dt;
+      }
       refx[i] = ref[i].x;
       refy[i] = ref[i].y;
       refz[i] = ref[i].z;
@@ -2350,7 +2390,6 @@ int HomePatch::rattle1(const BigReal timestep, Tensor *virial,
       posz[i] = pos[i].z;
     }
 
-
     bool done;
     bool consFailure;
     if (icnt == 1) {
@@ -2400,9 +2439,7 @@ int HomePatch::rattle1(const BigReal timestep, Tensor *virial,
         << (atom[ig].id + 1) << "!\n" << endi;
       }
     }
-
   }
-
   // Finally, we have to go through atoms that are not involved in rattle just so that we have
   // their positions and velocities up-to-date in posNew and velNew
   for (int j=0;j < noconstList.size();++j) {
@@ -3497,13 +3534,27 @@ void HomePatch::submitLoadStats(int timestep)
 }
 
 
+//
+// XXX operates on CompAtom, not FullAtom
+//
 void HomePatch::doPairlistCheck()
 {
+#if 0
+#if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED)
+  char dpcbuf[32];
+  sprintf(dpcbuf, "%s: %d", NamdProfileEventStr[NamdProfileEvent::DO_PAIRLIST_CHECK], this->getPatchID());
+  NAMD_EVENT_START_EX(1, NamdProfileEvent::DO_PAIRLIST_CHECK, dpcbuf);
+#endif
+#endif
+
   SimParameters *simParams = Node::Object()->simParameters;
 
   if ( numAtoms == 0 || ! flags.usePairlists ) {
     flags.pairlistTolerance = 0.;
     flags.maxAtomMovement = 99999.;
+#if 0
+    NAMD_EVENT_STOP(1, NamdProfileEvent::DO_PAIRLIST_CHECK);
+#endif
     return;
   }
 
@@ -3518,6 +3569,9 @@ void HomePatch::doPairlistCheck()
     doPairlistCheck_positions.resize(numAtoms);
     CompAtom *psave_i = doPairlistCheck_positions.begin();
     for ( i=0; i<n; ++i ) { psave_i[i] = p_i[i]; }
+#if 0
+    NAMD_EVENT_STOP(1, NamdProfileEvent::DO_PAIRLIST_CHECK);
+#endif
     return;
   }
 
@@ -3569,7 +3623,9 @@ void HomePatch::doPairlistCheck()
   if ( max_tol > doPairlistCheck_newTolerance ) {
     doPairlistCheck_newTolerance = max_tol / (1. - simParams->pairlistTrigger);
   }
-
+#if 0
+  NAMD_EVENT_STOP(1, NamdProfileEvent::DO_PAIRLIST_CHECK);
+#endif
 }
 
 void HomePatch::doGroupSizeCheck()
@@ -3587,7 +3643,7 @@ void HomePatch::doGroupSizeCheck()
     const int hgs = p_i->hydrogenGroupSize;
     if ( ! hgs ) break;  // avoid infinite loop on bug
     int ngs = hgs;
-    if ( ngs > 5 ) ngs = 5;  // limit to at most 5 atoms per group
+    if ( ngs > 5 ) ngs = 5;  // XXX why? limit to at most 5 atoms per group
     BigReal x = p_i->position.x;
     BigReal y = p_i->position.y;
     BigReal z = p_i->position.z;
@@ -3726,9 +3782,19 @@ HomePatch::doAtomMigration()
   #endif
 
   while ( atom_i != atom_e ) {
+    // Even though this code iterates through all atoms successively
+    // it moves entire hydrogen/migration groups as follows:
+    // Only the parent atom of the hydrogen/migration group has
+    // nonzero migrationGroupSize.  Values determined for xdev,ydev,zdev
+    // will persist through the remaining group members so that each
+    // following atom will again be added to the same mList.
     if ( atom_i->migrationGroupSize ) {
       Position pos = atom_i->position;
       if ( atom_i->migrationGroupSize != atom_i->hydrogenGroupSize ) {
+        // If there are multiple hydrogen groups in a migration group
+        // (e.g. for supporting lone pairs)
+        // the following code takes the average position (midpoint)
+        // of their parents.
         int mgs = atom_i->migrationGroupSize;
         int c = 1;
         for ( int j=atom_i->hydrogenGroupSize; j<mgs;
@@ -3740,6 +3806,11 @@ HomePatch::doAtomMigration()
         // iout << "mgroup " << atom_i->id << " at " << pos << "\n" << endi;
       }
 
+      // Scaling the position below transforms space within patch from
+      // what could have been a rotated parallelepiped into
+      // orthogonal coordinates, where we can use minmax comparison
+      // to detect which of our nearest neighbors this
+      // parent atom might have entered.
       ScaledPosition s = lattice.scale(pos);
 
       // check if atom is within bounds
@@ -3785,6 +3856,8 @@ HomePatch::doAtomMigration()
 
 
     } else {
+      // By keeping track of delnum total being deleted from FullAtomList
+      // the else clause allows us to fill holes as we visit each atom.
 
       if ( delnum ) { *(atom_i-delnum) = *atom_i; }
 
@@ -3816,7 +3889,9 @@ HomePatch::doAtomMigration()
      depositMigration(msgbuf[i]);
   }
   numMlBuf = 0;
-     
+
+  NAMD_EVENT_STOP(1, NamdProfileEvent::ATOM_MIGRATIONS);
+
   if (!allMigrationIn) {
     DebugM(3,"All Migrations NOT in, we are suspending patch "<<patchID<<"\n");
     migrationSuspended = true;
@@ -3926,6 +4001,7 @@ HomePatch::depositMigration(MigrateAtomsMsg *msg)
 //   atom list (regardless of whether or not the atom list is has been
 //   sorted yet or not).
 void HomePatch::separateAtoms() {
+  SimParameters *simParams = Node::Object()->simParameters;
 
   // Basic Idea:  Iterate through all the atoms in the current list
   //   of atoms.  Pack the waters in the current atoms list and move
@@ -4003,6 +4079,7 @@ void HomePatch::separateAtoms() {
 // NOTE: This function applies the transformations to the incoming
 //   atoms as it is separating them.
 void HomePatch::mergeAtomList(FullAtomList &al) {
+  SimParameters *simParams = Node::Object()->simParameters;
 
   // Sanity check
   if (al.size() <= 0) return;  // Nothing to do
index 373a8b5..18b5129 100644 (file)
 #include <string>
 #include <map>
 
-//
-// DJH: NAMDLite array buffers are memory aligned for vector instructions.
-//
-//#include "nl_Array.h"
-//
-
 class RegisterProxyMsg;
 class UnregisterProxyMsg;
 class ProxyResultVarsizeMsg;
@@ -49,41 +43,232 @@ class ProxyNodeAwareSpanningTreeMsg;
 
 class ComputeQMMgr;
 
+
+#ifdef TIMER_COLLECTION
+
+#include <time.h>
+
+struct TimerMicrosecond {
+  struct timespec ts;
+  inline void start() {
+    clock_gettime(CLOCK_REALTIME, &ts);
+  }
+  inline double stop() {
+    struct timespec tsend;
+    clock_gettime(CLOCK_REALTIME, &tsend);
+    return( (tsend.tv_sec - ts.tv_sec) * 1e6      // sec to microsec
+        + (tsend.tv_nsec - ts.tv_nsec) * 1e-3 );  // nanosec to microsec
+  }
+};
+
+#define TIMER_SLOTS  101
+#define TIMER_SLOT_WIDTH  1
+
+/// Timer entry
+struct TimerEntry {
+  TimerMicrosecond tmicro;
+  double tcur;   ///< timer current entry
+  double tavg;   ///< timer average
+  double tvar;   ///< timer variance
+  double tstd;   ///< timer standard deviation
+  double tmin;   ///< timer minimum
+  double tmax;   ///< timer maximum
+  double tsum;   ///< timer summation
+#if defined(DEBUG_TIMER_COLLECTION)
+  double tcharm; ///< compare with timer from Charm++ for debugging
+#endif
+#if defined(TIMER_HISTOGRAM)
+  double slotwidth;
+  double inv_slotwidth;
+  int hist[TIMER_SLOTS];
+#endif
+  int count;
+  TimerEntry() { reset(); }
+  /// Reset accumulators.
+  inline void reset() {
+    memset(this, 0, sizeof(TimerEntry));
+  }
+  inline void init(double t = TIMER_SLOT_WIDTH) {
+    tmin = 1e20;
+#if defined(TIMER_HISTOGRAM)
+    slotwidth = t;
+    inv_slotwidth = (slotwidth > 0 ? 1./slotwidth : 0);
+#endif
+  }
+  /// Start timer, set tcur to seconds since program started.
+  inline void start() {
+#if defined(DEBUG_TIMER_COLLECTION)
+    tcharm = CkWallTimer();
+#endif
+    tmicro.start();
+  }
+  /// Stop timer, set tcur to time difference between now and start.
+  inline void stop() {
+    tcur = tmicro.stop();
+#if defined(DEBUG_TIMER_COLLECTION)
+    tcharm = CkWallTimer() - tcharm;  // find ellapsed time
+    tcharm *= 1e6;  // convert to microseconds
+#endif
+  }
+  /// Update by including tcur into running statistics.
+  /// tavg is running average
+  /// tvar is unscaled variance
+  /// tmin is minimum value so far
+  /// tmax is maximum value so far
+  inline void update() {
+    count++;
+    tsum += tcur;
+    double delta = tcur - tavg;
+    tavg = tavg + delta / count;
+    double delta2 = tcur - tavg;
+    tvar += delta * delta2;
+    if (tcur > tmax) tmax = tcur;
+    if (tcur < tmin) tmin = tcur;
+#if defined(TIMER_HISTOGRAM)
+    int index = int(floor(tcur * inv_slotwidth));
+    if (index >= TIMER_SLOTS) index = TIMER_SLOTS - 1;
+    hist[index]++;
+#endif
+  }
+  /// Finalize the statistics.
+  /// tvar becomes scaled by the count
+  /// tstd calculated from tvar
+  inline void finalize() {
+    if (count > 0) tvar /= count;
+    tstd = sqrt(tvar);
+    if (tmin > tmax) tmin = tmax;
+  }
+};
+
+struct TimerSet {
+  enum {
+    KICK,
+    MAXMOVE,
+    DRIFT,
+    PISTON,
+    SUBMITHALF,
+    VELBBK1,
+    VELBBK2,
+    RATTLE1,
+    SUBMITFULL,
+    SUBMITCOLLECT,
+    NUMTIMERS
+  };
+  TimerEntry t[NUMTIMERS];
+  static const char *tlabel[NUMTIMERS];
+};
+
+#define TIMER_INIT(T,TYPE) \
+  do { \
+    (T).t[TimerSet::TYPE].init(); \
+  } while(0)
+
+#define TIMER_INIT_WIDTH(T,TYPE,WIDTH) \
+  do { \
+    (T).t[TimerSet::TYPE].init(WIDTH); \
+  } while(0)
+
+#define TIMER_START(T,TYPE) \
+  do { \
+    (T).t[TimerSet::TYPE].start(); \
+  } while(0)
+
+#if defined(DEBUG_TIMER_COLLECTION)
+
+// For debugging, compare clock_gettime with CkWallTimer.
+// The concern is regarding these routines that are on average less than
+// 10 us (microseconds) but are sometimes an order of magnitude slower.
 //
-// DJH: Array buffers defined here for storing data from FullAtomList and
-// also forces into SOA (structure of arrays) layout. Also declare methods
-// for copying from AOS to SOA and back again.
-//
-// We will also remove derived constants from FullAtom (e.g. recipMass).
-// The idea is reduce the messaging footprint as much as possible.
-// Recalculate constants after atom migration.
-//
-//struct PatchDataSOA {
-//
-//namdlite::Array<float> gaussrand; // fill with Gaussian random numbers
-//
-//namdlite::Array<float> mass;
-//namdlite::Array<float> recipMass; // derived from mass
-//namdlite::Array<float> langevinParam;
-//namdlite::Array<float> langScalVelBBK2;  // derived from langevinParam
-//namdlite::Array<float> langScalRandBBK2; // from langevinParam and recipMass
+// Note:  After testing, everything with use of clock_gettime seems
+// to be working correctly.
 //
-//namdlite::Array<double> vel_x;  // Jim recommends double precision velocity
-//namdlite::Array<double> vel_y;
-//namdlite::Array<double> vel_z;
-//namdlite::Array<double> pos_x;
-//namdlite::Array<double> pos_y;
-//namdlite::Array<double> pos_z;
-//namdlite::Array<double> f_normal_x;
-//namdlite::Array<double> f_normal_y;
-//namdlite::Array<double> f_normal_z;
-//namdlite::Array<double> f_nbond_x;
-//namdlite::Array<double> f_nbond_y;
-//namdlite::Array<double> f_nbond_z;
-//namdlite::Array<double> f_slow_x;
-//namdlite::Array<double> f_slow_y;
-//namdlite::Array<double> f_slow_z;
-//};
+#define TIMER_STOP(T,TYPE)  \
+  do { \
+    (T).t[TimerSet::TYPE].stop(); \
+    (T).t[TimerSet::TYPE].update(); \
+    double tcur = (T).t[TimerSet::TYPE].tcur; \
+    int count = (T).t[TimerSet::TYPE].count; \
+    if (tcur >= 100 && patch->patchID == SPECIAL_PATCH_ID) { \
+      printf("*** %s timing: %g   count: %d  line: %d  charm: %g\n", \
+          (T).tlabel[TimerSet::TYPE], tcur, count, __LINE__, \
+          (T).t[TimerSet::TYPE].tcharm); \
+    } \
+  } while(0)
+
+#else   // no DEBUG
+
+#define TIMER_STOP(T,TYPE)  \
+  do { \
+    (T).t[TimerSet::TYPE].stop(); \
+    (T).t[TimerSet::TYPE].update(); \
+  } while(0)
+
+#endif  // DEBUG_TIMER_COLLECTION
+
+#define TIMER_DONE(T) \
+  do { \
+    for (int i=0;  i < TimerSet::NUMTIMERS;  i++) { \
+      (T).t[i].finalize(); \
+    } \
+  } while(0)
+
+#if defined(TIMER_HISTOGRAM)
+
+#define TIMER_REPORT(T) \
+  do { \
+    printf("%13s %11s %11s %8s %8s %11s %8s\n", \
+        "name", "avg", "std", "min", "max", "sum", "calls"); \
+    printf("---------------------------------------------------------------" \
+        "-------------\n"); \
+    for (int i=0;  i < TimerSet::NUMTIMERS;  i++) { \
+      printf("%13s %11g %11g %8g %8g %11g %8d\n", \
+          (T).tlabel[i], (T).t[i].tavg, (T).t[i].tstd, \
+          (T).t[i].tmin, (T).t[i].tmax, (T).t[i].tsum, (T).t[i].count); \
+    } \
+    printf("---------------------------------------------------------------" \
+        "-------------\n"); \
+    for (int i=0;  i < TimerSet::NUMTIMERS;  i++) { \
+      printf("%13s %8s %8s %8s\n", \
+          (T).tlabel[i], "slot", "time", "count"); \
+      for (int j=0;  j < TIMER_SLOTS;  j++) { \
+        printf("%13s %8d %8g %8d\n", \
+            " ", j, (j+1)*(T).t[i].slotwidth, (T).t[i].hist[j]); \
+      } \
+      printf("---------------------------------------------------------------" \
+          "-------------\n"); \
+    } \
+  } while(0)
+
+#else  // no HISTOGRAM
+
+#define TIMER_REPORT(T) \
+  do { \
+    printf("%13s %11s %11s %8s %8s %11s %8s\n", \
+        "name", "avg", "std", "min", "max", "sum", "calls"); \
+    printf("---------------------------------------------------------------" \
+        "-------------\n"); \
+    for (int i=0;  i < TimerSet::NUMTIMERS;  i++) { \
+      printf("%13s %11g %11g %8g %8g %11g %8d\n", \
+          (T).tlabel[i], (T).t[i].tavg, (T).t[i].tstd, \
+          (T).t[i].tmin, (T).t[i].tmax, (T).t[i].tsum, (T).t[i].count); \
+    } \
+    printf("---------------------------------------------------------------" \
+        "-------------\n"); \
+  } while(0)
+
+#endif  // TIMER_HISTOGRAM
+
+#else   // no TIMER
+
+#define TIMER_INIT(T,TYPE)   do { } while(0)
+#define TIMER_INIT_WIDTH(T,TYPE,WIDTH)  do{ } while(0)
+#define TIMER_START(T,TYPE)  do { } while(0)
+#define TIMER_STOP(T,TYPE)   do { } while(0)
+#define TIMER_DONE(T)        do { } while(0)
+#define TIMER_REPORT(T)      do { } while(0)
+
+#endif  // TIMER_COLLECTION
+
 
 class HomePatch : public Patch {
   friend class PatchMgr;
@@ -292,6 +477,10 @@ public:
 #endif
 
   LDObjHandle ldObjHandle;
+
+#ifdef TIMER_COLLECTION
+  TimerSet timerSet;
+#endif
 protected:
   virtual void boxClosed(int);
 
@@ -312,25 +501,6 @@ private:
 
   CudaAtomList cudaAtomList;
 
-  //
-  // DJH: SOA data structure declared here.
-  //
-  //PatchDataSOA patchDataSOA;
-  //
-  // Copy fields from FullAtom into SOA form.
-  //void copy_atoms_to_SOA();
-  //
-  // Copy forces into SOA form.
-  //void copy_forces_to_SOA();
-  //
-  // Calculate derived constants after atom migration.
-  //void calculate_derived_SOA();
-  //
-  // Copy the updated quantities, e.g., positions and velocities, from SOA
-  // back to AOS form.
-  //void copy_updates_to_AOS();
-  //
-
   // DMK - Atom Separation (water vs. non-water)
   #if NAMD_SeparateWaters != 0
     FullAtomList tempAtom;  // A temporary array used to sort waters
diff --git a/src/NamdEventsProfiling.def b/src/NamdEventsProfiling.def
new file mode 100644 (file)
index 0000000..8ab6b2d
--- /dev/null
@@ -0,0 +1,63 @@
+NAMD_PROFILE_EVENT(DUMMY_EVENT, "Dummy Event")
+
+#ifdef SEQUENCER_SOA
+NAMD_PROFILE_EVENT(INTEGRATE_SOA_1, "integrate_SOA_1")
+NAMD_PROFILE_EVENT(INTEGRATE_SOA_2, "integrate_SOA_2")
+NAMD_PROFILE_EVENT(INTEGRATE_SOA_3, "integrate_SOA_3")
+NAMD_PROFILE_EVENT(ADD_FORCE_TO_MOMENTUM_SOA, "addForceToMomentum_SOA")
+NAMD_PROFILE_EVENT(ADD_VELOCITY_TO_POSITION_SOA, "addVelocityToPosition_SOA")
+NAMD_PROFILE_EVENT(POSITIONS_READY_SOA, "positionsReady_SOA")
+NAMD_PROFILE_EVENT(SUBMIT_HALFSTEP_SOA, "submitHalfstep_SOA")
+NAMD_PROFILE_EVENT(SUBMIT_REDUCTIONS_SOA, "submitReductions_SOA")
+NAMD_PROFILE_EVENT(SUBMIT_COLLECTIONS_SOA, "submitCollections_SOA")
+NAMD_PROFILE_EVENT(MAXIMUM_MOVE_SOA, "maximumMove_SOA")
+NAMD_PROFILE_EVENT(LANGEVIN_VELOCITIES_BBK1_SOA, "langevinVelocitiesBBK1_SOA")
+NAMD_PROFILE_EVENT(LANGEVIN_VELOCITIES_BBK2_SOA, "langevinVelocitiesBBK2_SOA")
+NAMD_PROFILE_EVENT(RATTLE1_SOA, "rattle1_SOA")
+NAMD_PROFILE_EVENT(COPY_FORCES_TO_SOA, "copy_forces_to_SOA")
+#endif
+
+NAMD_PROFILE_EVENT(DUMMY_EVENT14, "Dummy Event14")
+NAMD_PROFILE_EVENT(DUMMY_EVENT15, "Dummy Event15")
+NAMD_PROFILE_EVENT(DUMMY_EVENT16, "Dummy Event16")
+NAMD_PROFILE_EVENT(DUMMY_EVENT17, "Dummy Event17")
+NAMD_PROFILE_EVENT(DUMMY_EVENT18, "Dummy Event18")
+NAMD_PROFILE_EVENT(DUMMY_EVENT19, "Dummy Event19")
+NAMD_PROFILE_EVENT(DUMMY_EVENT20, "Dummy Event20")
+NAMD_PROFILE_EVENT(DUMMY_EVENT21, "Dummy Event21")
+NAMD_PROFILE_EVENT(DUMMY_EVENT22, "Dummy Event22")
+
+NAMD_PROFILE_EVENT(INTEGRATE_1, "integrate 1")
+NAMD_PROFILE_EVENT(INTEGRATE_2, "integrate 2")
+NAMD_PROFILE_EVENT(INTEGRATE_3, "integrate 3")
+NAMD_PROFILE_EVENT(ADD_FORCE_TO_MOMENTUM, "addForceToMomentum")
+NAMD_PROFILE_EVENT(ADD_VELOCITY_TO_POSITION, "addVelocityToPosition")
+NAMD_PROFILE_EVENT(POSITIONS_READY, "positionsReady")
+NAMD_PROFILE_EVENT(SUBMIT_HALFSTEP, "submitHalfstep")
+NAMD_PROFILE_EVENT(SUBMIT_REDUCTIONS, "submitReductions")
+NAMD_PROFILE_EVENT(SUBMIT_COLLECTIONS, "submitCollections")
+NAMD_PROFILE_EVENT(MAXIMUM_MOVE, "maximumMove")
+NAMD_PROFILE_EVENT(LANGEVIN_VELOCITIES_BBK1, "langevinVelocitiesBBK1")
+NAMD_PROFILE_EVENT(LANGEVIN_VELOCITIES_BBK2, "langevinVelocitiesBBK2")
+NAMD_PROFILE_EVENT(RATTLE1, "rattle1")
+
+NAMD_PROFILE_EVENT(NEWTONIAN_VELOCITIES, "newtonianVelocities")
+NAMD_PROFILE_EVENT(ATOM_MIGRATIONS, "atomMigrations")
+NAMD_PROFILE_EVENT(DO_NONBONDED_CUDA_WORK, "calcNonBondedForces")
+NAMD_PROFILE_EVENT(DO_BONDED_CUDA_WORK, "calcBondedForces")
+NAMD_PROFILE_EVENT(COMPUTE_PME_CUDA, "computePMECuda")
+NAMD_PROFILE_EVENT(NAMD_STARTUP, "NAMD startup")
+NAMD_PROFILE_EVENT(MOLECULE_CONSTRUCTOR, "moleculeConstructor")
+NAMD_PROFILE_EVENT(SPREAD_CHARGE, "spreadCharge")
+NAMD_PROFILE_EVENT(GATHER_FORCE, "gatherForce")
+NAMD_PROFILE_EVENT(REDUCE_VIRIAL_ENERGY, "reduceVirialEnergy")
+NAMD_PROFILE_EVENT(COMPUTE_BONDED_CUDA_FINISH_PATCHES, "computeBonded_FinishPatchesOnPe")
+NAMD_PROFILE_EVENT(COMPUTE_BONDED_CUDA_FINISH_PATCHES_LOCK, "computeBonded_FinishPatchesOnPe_Lock")
+NAMD_PROFILE_EVENT(COMPUTE_BONDED_CUDA_OPEN_BOXES, "computeBonded_OpenBoxesOnPe")
+NAMD_PROFILE_EVENT(COMPUTE_NONBONDED_CUDA_COPY_ATOMS_SUBSET, "computeNonBonded_CopyAtomsSubset")
+NAMD_PROFILE_EVENT(COMPUTE_NONBONDED_CUDA_FINISH_PATCHES, "computeNonBonded_FinishPatchesOnPe")
+NAMD_PROFILE_EVENT(COMPUTE_NONBONDED_OPEN_BOXES, "computeNonBonded_OpenBox")
+NAMD_PROFILE_EVENT(SORT_SOLVENT_ATOMS, "sortSolventAtoms")
+NAMD_PROFILE_EVENT(DO_PAIRLIST_CHECK, "doPairlistCheck")
+NAMD_PROFILE_EVENT(COMPUTE_BONDED_LOAD_TUPLES, "computeBonded_loadTuplesOnPe")
+NAMD_PROFILE_EVENT(COMPUTE_BONDED_TUPLE_COPY, "computeBonded_tupleCopy")
diff --git a/src/NamdEventsProfiling.h b/src/NamdEventsProfiling.h
new file mode 100644 (file)
index 0000000..6018382
--- /dev/null
@@ -0,0 +1,254 @@
+/**
+***  Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by
+***  The Board of Trustees of the University of Illinois.
+***  All rights reserved.
+**/
+
+#ifndef NAMDEVENTSPROFILING_H
+#define NAMDEVENTSPROFILING_H
+
+#include "common.h"
+#include <map>
+
+//
+// Charm++ Projections requires user defined trace events to be registered
+// with event ID tags and string names.
+//
+// Use the NAMD_PROFILE_EVENT macro within NamdEventsProfiling.def
+// to define these event ID tags and string names.  The IDs generate
+// a set of enumerated types and corresponding array of string names.
+//
+struct NamdProfileEvent {
+  typedef enum {
+    #define NAMD_PROFILE_EVENT(a,b) a,
+    #include "NamdEventsProfiling.def"
+    #undef NAMD_PROFILE_EVENT
+    EventsCount
+  } Event;
+};
+
+char const* const NamdProfileEventStr[] = {
+  #define NAMD_PROFILE_EVENT(a,b) b,
+  #include "NamdEventsProfiling.def"
+  #undef NAMD_PROFILE_EVENT
+  0
+};
+
+#undef NAMD_PROFILE_START
+#undef NAMD_PROFILE_STOP
+#undef NAMD_REGISTER_EVENT
+#undef NAMD_EVENT_START
+#undef NAMD_EVENT_START_EX
+#undef NAMD_EVENT_STOP
+#undef NAMD_EVENT_RANGE
+#undef NAMD_EVENT_RANGE_2
+
+//
+// Enable NVTX instrumentation for nvvp or Nsight profiling
+// NAMD_CUDA build by defining NAMD_NVTX_ENABLED in Make.config
+//
+#if defined(NAMD_CUDA) && defined(NAMD_NVTX_ENABLED)
+
+#include <nvToolsExt.h>
+#include <cuda_profiler_api.h>
+
+// start profiling
+#define NAMD_PROFILE_START() \
+  do { \
+    cudaProfilerStart(); \
+  } while (0)  // must terminate with semi-colon
+
+// stop profiling
+#define NAMD_PROFILE_STOP() \
+  do { \
+    cudaProfilerStop(); \
+  } while (0)  // must terminate with semi-colon
+
+// C++ note: declaring const variables implies static (internal) linkage,
+// and you have to explicitly specify "extern" to get external linkage.
+const uint32_t NAMD_nvtx_colors[] = {
+  0x0000ff00,
+  0x000000ff,
+  0x00ffff00,
+  0x00ff00ff,
+  0x0000ffff,
+  0x00ff0000,
+  0x00006600,
+  0x00663300,
+  0x00000000,
+  0x007300e6,
+  0x00ff8c00,
+};
+const int NAMD_nvtx_colors_len = sizeof(NAMD_nvtx_colors)/sizeof(uint32_t);
+
+#define NAMD_REGISTER_EVENT(name,cid) \
+  do { } while(0)  // must terminate with semi-colon
+
+// start recording an event
+#define NAMD_EVENT_START(eon,id) \
+  do { \
+    if (eon) { \
+      int color_id = id; \
+      color_id = color_id % NAMD_nvtx_colors_len; \
+      nvtxEventAttributes_t eventAttrib = {0}; \
+      eventAttrib.version = NVTX_VERSION; \
+      eventAttrib.size = NVTX_EVENT_ATTRIB_STRUCT_SIZE; \
+      eventAttrib.colorType = NVTX_COLOR_ARGB; \
+      eventAttrib.color = NAMD_nvtx_colors[color_id]; \
+      eventAttrib.messageType = NVTX_MESSAGE_TYPE_ASCII; \
+      eventAttrib.message.ascii = NamdProfileEventStr[id]; \
+      nvtxRangePushEx(&eventAttrib); \
+    } \
+  } while(0)  // must terminate with semi-colon
+
+// start recording an event
+#define NAMD_EVENT_START_EX(eon,id,str) \
+  do { \
+    if (eon) { \
+      int color_id = id; \
+      color_id = color_id % NAMD_nvtx_colors_len; \
+      nvtxEventAttributes_t eventAttrib = {0}; \
+      eventAttrib.version = NVTX_VERSION; \
+      eventAttrib.size = NVTX_EVENT_ATTRIB_STRUCT_SIZE; \
+      eventAttrib.colorType = NVTX_COLOR_ARGB; \
+      eventAttrib.color = NAMD_nvtx_colors[color_id]; \
+      eventAttrib.messageType = NVTX_MESSAGE_TYPE_ASCII; \
+      eventAttrib.message.ascii = str; \
+      nvtxRangePushEx(&eventAttrib); \
+    } \
+  } while(0)  // must terminate with semi-colon
+
+// stop recording an event
+#define NAMD_EVENT_STOP(eon,id) \
+  do { \
+    if (eon) { \
+      nvtxRangePop(); \
+    } \
+  } while (0)  // must terminate with semi-colon
+
+// embed event recording in class to automatically pop when destroyed
+class NAMD_NVTX_Tracer {
+  protected:
+    int evon;  // is event on?
+  public:
+    NAMD_NVTX_Tracer(int eon, int id = 0) : evon(eon) {
+      NAMD_EVENT_START(eon, id);
+    }
+    ~NAMD_NVTX_Tracer() { NAMD_EVENT_STOP(evon, 0); }
+};
+
+// call NAMD_EVENT_RANGE at beginning of function to push event recording
+// destructor is automatically called on return to pop event recording
+#define NAMD_EVENT_RANGE(eon,id) \
+  NAMD_NVTX_Tracer namd_nvtx_tracer(eon,id)
+  // must terminate with semi-colon
+
+#if defined(NAMD_PROFILE_EVENT_LAYER_2)
+#define NAMD_EVENT_RANGE_2(eon,id) \
+  NAMD_EVENT_RANGE(eon,id)
+#else
+#define NAMD_EVENT_RANGE_2(eon,id) \
+  do { } while(0)  // must terminate with semi-colon
+#endif
+
+//
+// Enable Projections trace events when built against CMK_TRACE_ENABLED
+// version of Charm++ by defining NAMD_CMK_TRACE_ENABLED in Make.config
+//
+#elif defined(CMK_TRACE_ENABLED) && defined(NAMD_CMK_TRACE_ENABLED)
+
+// XXX should be more general than tracing Sequencer functions
+// XXX why offset event numbers by 150?
+#define SEQUENCER_EVENT_ID_START 150
+
+#define NAMD_PROFILE_START() \
+  do { } while(0)  // must terminate with semi-colon
+
+#define NAMD_PROFILE_STOP() \
+  do { } while(0)  // must terminate with semi-colon
+
+#define NAMD_REGISTER_EVENT(name,id) \
+  do { \
+    int eventID = SEQUENCER_EVENT_ID_START+id; \
+    traceRegisterUserEvent(name, eventID); \
+  } while(0) // must terminate with semi-colon
+
+#define NAMD_EVENT_START(eon,id) \
+  do {\
+    if (eon) { \
+      int eventID = SEQUENCER_EVENT_ID_START+id; \
+      traceBeginUserBracketEvent(eventID); \
+    } \
+  } while(0) // must terminate with semi-colon
+
+#define NAMD_EVENT_START_EX(eon,id,str) \
+  NAMD_EVENT_START(eon,id)
+
+#define NAMD_EVENT_STOP(eon,id) \
+  do { \
+    if (eon) { \
+      int eventID = SEQUENCER_EVENT_ID_START+id; \
+      traceEndUserBracketEvent(eventID); \
+    } \
+  } while(0)  // must terminate with semi-colon
+
+// XXX should be more general than Sequencer
+class NAMD_Sequencer_Events_Tracer {
+  int tEventID;
+  int tEventOn;
+  public:
+    NAMD_Sequencer_Events_Tracer(int eon, int id = 0)
+      : tEventOn(eon), tEventID(id) {
+      NAMD_EVENT_START(tEventOn, tEventID);
+    }
+    ~NAMD_Sequencer_Events_Tracer() {
+      NAMD_EVENT_STOP(tEventOn, tEventID);
+    }
+};
+
+// call NAMD_EVENT_RANGE at beginning of function to push event recording
+// destructor is automatically called on return to pop event recording
+#define NAMD_EVENT_RANGE(eon,id) \
+    NAMD_Sequencer_Events_Tracer namd_events_tracer(eon,id)
+    // must terminate with semi-colon
+
+#if defined(NAMD_PROFILE_EVENT_LAYER_2)
+#define NAMD_EVENT_RANGE_2(eon,id) \
+  NAMD_EVENT_RANGE(eon,id)
+#else
+#define NAMD_EVENT_RANGE_2(eon,id) \
+  do { } while(0)  // must terminate with semi-colon
+#endif
+
+#else
+
+//
+// Otherwise all profiling macros become no-ops.
+//
+#define NAMD_PROFILE_START() \
+  do { } while(0)  // must terminate with semi-colon
+
+#define NAMD_PROFILE_STOP() \
+  do { } while(0)  // must terminate with semi-colon
+
+#define NAMD_REGISTER_EVENT(name,cid) \
+  do { } while(0)  // must terminate with semi-colon
+
+#define NAMD_EVENT_START(eon,id) \
+  do { } while(0)  // must terminate with semi-colon
+
+#define NAMD_EVENT_START_EX(eon,id,str) \
+  do { } while(0)  // must terminate with semi-colon
+
+#define NAMD_EVENT_STOP(eon,id) \
+  do { } while(0)  // must terminate with semi-colon
+
+#define NAMD_EVENT_RANGE(eon,id) \
+  do { } while(0)  // must terminate with semi-colon
+
+#define NAMD_EVENT_RANGE_2(eon,id) \
+  do { } while(0)  // must terminate with semi-colon
+
+#endif // NAMD_CUDA && NAMD_NVTX_ENABLED
+
+#endif /* NAMDEVENTSPROFILING_H */
index aeb43b0..3fb031e 100644 (file)
@@ -31,6 +31,10 @@ public:
   int maxForceUsed;            // may ignore slower force classes
   int maxForceMerged;          // add this and faster to normal
 
+#if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED)
+  int event_on;  // true or false to control NVTX profiling
+#endif
+
   int usePairlists;
   int savePairlists;
   BigReal pairlistTolerance;
@@ -39,6 +43,10 @@ public:
 
   Lattice lattice;             // rather than shipping around separately
 
+#if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED)
+  Flags() : event_on(0) { }
+#endif
+
   inline Flags& operator= (const Flags &flags) 
     {
       CmiMemcpy (this, &flags, sizeof(Flags));
index cbe94ac..df0f9fe 100644 (file)
  * $Revision: 1.1230 $
  *****************************************************************************/
 
+// The UPPER_BOUND macro is used to eliminate all of the per atom
+// computation done for the numerical integration in Sequencer::integrate()
+// other than the actual force computation and atom migration.
+// The idea is to "turn off" the integration for doing performance
+// profiling in order to get an upper bound on the speedup available
+// by moving the integration parts to the GPU.
+//
+// Define it in the Make.config file, i.e. CXXOPTS += -DUPPER_BOUND
+// or simply uncomment the line below.
+//
+//#define UPPER_BOUND
+
 //for gbis debugging; print net force on each atom
 #define PRINT_FORCES 0
 
 #include "PatchMap.inl"
 #include "ComputeMgr.h"
 #include "ComputeGlobal.h"
+#include "NamdEventsProfiling.h"
 
 #define MIN_DEBUG_LEVEL 4
 //#define DEBUGM
+//
+// Define NL_DEBUG below to activate D_*() macros in integrate_SOA()
+// for debugging.
+//
+//#define NL_DEBUG
 #include "Debug.h"
 
 #if USE_HPM
@@ -43,6 +61,8 @@
 #define STOP_HPM_STEP   1500
 #endif
 
+#define SPECIAL_PATCH_ID  91
+
 Sequencer::Sequencer(HomePatch *p) :
        simParams(Node::Object()->simParameters),
        patch(p),
@@ -207,6 +227,20 @@ void Sequencer::integrate(int scriptTask) {
     //
     //patch->copy_all_to_SOA();
 
+#ifdef TIMER_COLLECTION
+    TimerSet& t = patch->timerSet;
+#endif
+    TIMER_INIT_WIDTH(t, KICK, simParams->timerBinWidth);
+    TIMER_INIT_WIDTH(t, MAXMOVE, simParams->timerBinWidth);
+    TIMER_INIT_WIDTH(t, DRIFT, simParams->timerBinWidth);
+    TIMER_INIT_WIDTH(t, PISTON, simParams->timerBinWidth);
+    TIMER_INIT_WIDTH(t, SUBMITHALF, simParams->timerBinWidth);
+    TIMER_INIT_WIDTH(t, VELBBK1, simParams->timerBinWidth);
+    TIMER_INIT_WIDTH(t, VELBBK2, simParams->timerBinWidth);
+    TIMER_INIT_WIDTH(t, RATTLE1, simParams->timerBinWidth);
+    TIMER_INIT_WIDTH(t, SUBMITFULL, simParams->timerBinWidth);
+    TIMER_INIT_WIDTH(t, SUBMITCOLLECT, simParams->timerBinWidth);
+
     int &step = patch->flags.step;
     step = simParams->firstTimestep;
 
@@ -247,6 +281,7 @@ void Sequencer::integrate(int scriptTask) {
                                        maxForceMerged = Results::slow;
     if ( doFullElectrostatics ) maxForceUsed = Results::slow;
 
+//#ifndef UPPER_BOUND
     const Bool accelMDOn = simParams->accelMDOn;
     const Bool accelMDdihe = simParams->accelMDdihe;
     const Bool accelMDdual = simParams->accelMDdual;
@@ -279,21 +314,29 @@ void Sequencer::integrate(int scriptTask) {
     // Do we need to return forces to TCL script or Colvar module?
     int doTcl = simParams->tclForcesOn;
        int doColvars = simParams->colvarsOn;
+//#endif
     ComputeGlobal *computeGlobal = Node::Object()->computeMgr->computeGlobalObject;
 
     // Bother to calculate energies?
     int &doEnergy = patch->flags.doEnergy;
     int energyFrequency = simParams->outputEnergies;
-
+#ifndef UPPER_BOUND
     const int reassignFreq = simParams->reassignFreq;
+#endif
 
     int &doVirial = patch->flags.doVirial;
     doVirial = 1;
 
   if ( scriptTask == SCRIPT_RUN ) {
 
+#ifndef UPPER_BOUND
 //    printf("Doing initial rattle\n");
+#ifndef UPPER_BOUND
+D_MSG("rattle1()");
+    TIMER_START(t, RATTLE1);
     rattle1(0.,0);  // enforce rigid bond constraints on initial positions
+    TIMER_STOP(t, RATTLE1);
+#endif
 
     if (simParams->lonepairs) {
       patch->atomMapper->registerIDsFullAtom(
@@ -303,47 +346,85 @@ void Sequencer::integrate(int scriptTask) {
     if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
        reassignVelocities(timestep,step);
     }
+#endif
 
     doEnergy = ! ( step % energyFrequency );
+#ifndef UPPER_BOUND
     if ( accelMDOn && !accelMDdihe ) doEnergy=1;
     //Update energy every timestep for adaptive tempering
     if ( adaptTempOn ) doEnergy=1;
+#endif
+D_MSG("runComputeObjects()");
     runComputeObjects(1,step<numberOfSteps); // must migrate here!
+#ifndef UPPER_BOUND
     rescaleaccelMD(step, doNonbonded, doFullElectrostatics); // for accelMD 
     adaptTempUpdate(step); // update adaptive tempering temperature
+#endif
 
+#ifndef UPPER_BOUND
     if ( staleForces || doTcl || doColvars ) {
       if ( doNonbonded ) saveForce(Results::nbond);
       if ( doFullElectrostatics ) saveForce(Results::slow);
     }
     if ( ! commOnly ) {
+D_MSG("newtonianVelocities()");
+      TIMER_START(t, KICK);
       newtonianVelocities(-0.5,timestep,nbondstep,slowstep,0,1,1);
+      TIMER_STOP(t, KICK);
     }
     minimizationQuenchVelocity();
+#ifndef UPPER_BOUND
+D_MSG("rattle1()");
+    TIMER_START(t, RATTLE1);
     rattle1(-timestep,0);
+    TIMER_STOP(t, RATTLE1);
+#endif
+D_MSG("submitHalfstep()");
+    TIMER_START(t, SUBMITHALF);
     submitHalfstep(step);
+    TIMER_STOP(t, SUBMITHALF);
     if ( ! commOnly ) {
+D_MSG("newtonianVelocities()");
+      TIMER_START(t, KICK);
       newtonianVelocities(1.0,timestep,nbondstep,slowstep,0,1,1);
+      TIMER_STOP(t, KICK);
     }
+D_MSG("rattle1()");
+    TIMER_START(t, RATTLE1);
     rattle1(timestep,1);
+    TIMER_STOP(t, RATTLE1);
     if (doTcl || doColvars)  // include constraint forces
       computeGlobal->saveTotalForces(patch);
+D_MSG("submitHalfstep()");
+    TIMER_START(t, SUBMITHALF);
     submitHalfstep(step);
+    TIMER_STOP(t, SUBMITHALF);
     if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
     if ( ! commOnly ) {
+D_MSG("newtonianVelocities()");
+      TIMER_START(t, KICK);
       newtonianVelocities(-0.5,timestep,nbondstep,slowstep,0,1,1);
+      TIMER_STOP(t, KICK);
     }
+#endif
+D_MSG("submitReductions()");
+    TIMER_START(t, SUBMITFULL);
     submitReductions(step);
+    TIMER_STOP(t, SUBMITFULL);
+#ifndef UPPER_BOUND
     if(0){ // if(traceIsOn()){
         traceUserEvent(eventEndOfTimeStep);
         sprintf(traceNote, "%s%d",tracePrefix,step); 
         traceUserSuppliedNote(traceNote);
     }
+#endif
     rebalanceLoad(step);
 
   } // scriptTask == SCRIPT_RUN
 
+#ifndef UPPER_BOUND
   bool doMultigratorRattle = false;
+#endif
 
   //
   // DJH: There are a lot of mod operations below and elsewhere to
@@ -351,10 +432,36 @@ void Sequencer::integrate(int scriptTask) {
   // Mod and integer division are expensive!
   // Might be better to replace with counters and test equality.
   //
+#if 0
+    for(int i = 0; i < NamdProfileEvent::EventsCount; i++)
+       CkPrintf("-------------- [%d] %s -------------\n", i, NamdProfileEventStr[i]);
+#endif
 
+#if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED)
+  int& eon = patch->flags.event_on;
+  int epid = (simParams->beginEventPatchID <= patch->getPatchID()
+      && patch->getPatchID() <= simParams->endEventPatchID);
+  int beginStep = simParams->beginEventStep;
+  int endStep = simParams->endEventStep;
+  bool controlProfiling = patch->getPatchID() == 0;
+#endif
     for ( ++step; step <= numberOfSteps; ++step )
     {
-      PUSH_RANGE("integrate 1", 0);
+#if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED)
+      eon = epid && (beginStep < step && step <= endStep);
+
+      if (controlProfiling && step == beginStep) {
+        NAMD_PROFILE_START();
+      }
+      if (controlProfiling && step == endStep) {
+        NAMD_PROFILE_STOP();
+      }
+      char buf[32];
+      sprintf(buf, "%s: %d", NamdProfileEventStr[NamdProfileEvent::INTEGRATE_1], patch->getPatchID());
+      NAMD_EVENT_START_EX(eon, NamdProfileEvent::INTEGRATE_1, buf);
+#endif
+#ifndef UPPER_BOUND
 
       rescaleVelocities(step);
       tcoupleVelocities(timestep,step);
@@ -362,7 +469,9 @@ void Sequencer::integrate(int scriptTask) {
       berendsenPressure(step);
 
       if ( ! commOnly ) {
+        TIMER_START(t, KICK);
         newtonianVelocities(0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics); 
+        TIMER_STOP(t, KICK);
       }
 
       // We do RATTLE here if multigrator thermostat was applied in the previous step
@@ -379,12 +488,18 @@ void Sequencer::integrate(int scriptTask) {
          rattle1(timestep,0);
          } */
 
+      TIMER_START(t, MAXMOVE);
       maximumMove(timestep);
+      TIMER_STOP(t, MAXMOVE);
 
-      POP_RANGE;  // integrate 1
+      NAMD_EVENT_STOP(eon, NamdProfileEvent::INTEGRATE_1);  // integrate 1
 
       if ( simParams->langevinPistonOn || (simParams->langevinOn && simParams->langevin_useBAOAB) ) {
-        if ( ! commOnly ) addVelocityToPosition(0.5*timestep);
+        if ( ! commOnly ) {
+          TIMER_START(t, DRIFT);
+          addVelocityToPosition(0.5*timestep);
+          TIMER_STOP(t, DRIFT);
+        }
         // We add an Ornstein-Uhlenbeck integration step for the case of BAOAB (Langevin)
         langevinVelocities(timestep);
 
@@ -393,22 +508,32 @@ void Sequencer::integrate(int scriptTask) {
         // so split profiling around this conditional block.
         langevinPiston(step);
 
-        if ( ! commOnly ) addVelocityToPosition(0.5*timestep);
+        if ( ! commOnly ) {
+          TIMER_START(t, DRIFT);
+          addVelocityToPosition(0.5*timestep);
+          TIMER_STOP(t, DRIFT);
+        }
       } else {
         // If Langevin is not used, take full time step directly instread of two half steps      
-        if ( ! commOnly ) addVelocityToPosition(timestep); 
+        if ( ! commOnly ) {
+          TIMER_START(t, DRIFT);
+          addVelocityToPosition(timestep); 
+          TIMER_STOP(t, DRIFT);
+        }
       }
 
-      PUSH_RANGE("integrate 2", 1);
+      NAMD_EVENT_START(eon, NamdProfileEvent::INTEGRATE_2);
 
       // impose hard wall potential for Drude bond length
       hardWallDrude(timestep, 1);
 
       minimizationQuenchVelocity();
+#endif // UPPER_BOUND
 
       doNonbonded = !(step%nonbondedFrequency);
       doFullElectrostatics = (dofull && !(step%fullElectFrequency));
 
+#ifndef UPPER_BOUND 
       if ( zeroMomentum && doFullElectrostatics ) {
         // There is a blocking receive inside of correctMomentum().
         correctMomentum(step,slowstep);
@@ -416,7 +541,9 @@ void Sequencer::integrate(int scriptTask) {
 
       // There are NO sends in submitHalfstep() just local summation 
       // into the Reduction struct.
+      TIMER_START(t, SUBMITHALF);
       submitHalfstep(step);
+      TIMER_STOP(t, SUBMITHALF);
 
       doMolly = simParams->mollyOn && doFullElectrostatics;
       // BEGIN LA
@@ -444,14 +571,15 @@ void Sequencer::integrate(int scriptTask) {
         doKineticEnergy = 1;
         doMomenta = 1;
       }
-
-      POP_RANGE;  // integrate 2
+#endif
+      NAMD_EVENT_STOP(eon, NamdProfileEvent::INTEGRATE_2);  // integrate 2
 
       // The current thread of execution will suspend in runComputeObjects().
       runComputeObjects(!(step%stepsPerCycle),step<numberOfSteps);
 
-      PUSH_RANGE("integrate 3", 2);
+      NAMD_EVENT_START(eon, NamdProfileEvent::INTEGRATE_3);
 
+#ifndef UPPER_BOUND
       rescaleaccelMD(step, doNonbonded, doFullElectrostatics); // for accelMD
      
       if ( staleForces || doTcl || doColvars ) {
@@ -467,30 +595,48 @@ void Sequencer::integrate(int scriptTask) {
       }
 
       if ( ! commOnly ) {
+        TIMER_START(t, VELBBK1);
         langevinVelocitiesBBK1(timestep);
+        TIMER_STOP(t, VELBBK1);
+        TIMER_START(t, KICK);
         newtonianVelocities(1.0,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
+        TIMER_STOP(t, KICK);
+        TIMER_START(t, VELBBK2);
         langevinVelocitiesBBK2(timestep);
+        TIMER_STOP(t, VELBBK2);
       }
 
       // add drag to each atom's positions
       if ( ! commOnly && movDragOn ) addMovDragToPosition(timestep);
       if ( ! commOnly && rotDragOn ) addRotDragToPosition(timestep);
 
+      TIMER_START(t, RATTLE1);
       rattle1(timestep,1);
+      TIMER_STOP(t, RATTLE1);
       if (doTcl || doColvars)  // include constraint forces
         computeGlobal->saveTotalForces(patch);
 
+      TIMER_START(t, SUBMITHALF);
       submitHalfstep(step);
+      TIMER_STOP(t, SUBMITHALF);
       if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
 
       if ( ! commOnly ) {
+        TIMER_START(t, KICK);
         newtonianVelocities(-0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
+        TIMER_STOP(t, KICK);
       }
 
        // rattle2(timestep,step);
+#endif
 
+        TIMER_START(t, SUBMITFULL);
        submitReductions(step);
+        TIMER_STOP(t, SUBMITFULL);
+        TIMER_START(t, SUBMITCOLLECT);
        submitCollections(step);
+        TIMER_STOP(t, SUBMITCOLLECT);
+#ifndef UPPER_BOUND
        //Update adaptive tempering temperature
         adaptTempUpdate(step);
 
@@ -501,7 +647,8 @@ void Sequencer::integrate(int scriptTask) {
       multigratorTemperature(step, 2);
       doMultigratorRattle = (simParams->multigratorOn && !(step % simParams->multigratorTemperatureFreq));
 
-      POP_RANGE;  // integrate 3
+      NAMD_EVENT_STOP(eon, NamdProfileEvent::INTEGRATE_3); // integrate 3
+#endif
 
 #if CYCLE_BARRIER
         cycleBarrier(!((step+1) % stepsPerCycle), step);
@@ -511,6 +658,7 @@ void Sequencer::integrate(int scriptTask) {
         cycleBarrier(1, step);
 #endif
 
+#ifndef UPPER_BOUND
         if(Node::Object()->specialTracing || simParams->statsOn){
                 int bstep = simParams->traceStartStep;
                 int estep = bstep + simParams->numTraceSteps;
@@ -534,6 +682,7 @@ void Sequencer::integrate(int scriptTask) {
             sprintf(traceNote, "%s%d",tracePrefix,step); 
             traceUserSuppliedNote(traceNote);
         }
+#endif // UPPER_BOUND
        rebalanceLoad(step);
 
 #if PME_BARRIER
@@ -551,6 +700,14 @@ void Sequencer::integrate(int scriptTask) {
 
     }
 
+  TIMER_DONE(t);
+#ifdef TIMER_COLLECTION
+  if (patch->patchID == SPECIAL_PATCH_ID) {
+    printf("Timer collection reporting in microseconds for "
+        "Patch %d\n", patch->patchID);
+    TIMER_REPORT(t);
+  }
+#endif // TIMER_COLLECTION
     //
     // DJH: Copy updates of SOA back into AOS.
     //
@@ -1140,7 +1297,8 @@ void Sequencer::newtonianVelocities(BigReal stepscale, const BigReal timestep,
                                     const int doNonbonded,
                                     const int doFullElectrostatics)
 {
-  RANGE("newtonianVelocities", 3);
+  NAMD_EVENT_RANGE_2(patch->flags.event_on,
+      NamdProfileEvent::NEWTONIAN_VELOCITIES);
 
   // Deterministic velocity update, account for multigrator
   if (staleForces || (doNonbonded && doFullElectrostatics)) {
@@ -1195,7 +1353,8 @@ void Sequencer::langevinVelocities(BigReal dt_fs)
 
 void Sequencer::langevinVelocitiesBBK1(BigReal dt_fs)
 {
-  RANGE("langevinVelocitiesBBK1", 7);
+  NAMD_EVENT_RANGE_2(patch->flags.event_on,
+      NamdProfileEvent::LANGEVIN_VELOCITIES_BBK1);
   if ( simParams->langevinOn && !simParams->langevin_useBAOAB )
   {
     FullAtom *a = patch->atom.begin();
@@ -1267,14 +1426,17 @@ void Sequencer::langevinVelocitiesBBK1(BigReal dt_fs)
 
 void Sequencer::langevinVelocitiesBBK2(BigReal dt_fs)
 {
-  RANGE("langevinVelocitiesBBK2", 8);
+  NAMD_EVENT_RANGE_2(patch->flags.event_on,
+      NamdProfileEvent::LANGEVIN_VELOCITIES_BBK2);
   if ( simParams->langevinOn && !simParams->langevin_useBAOAB ) 
   {
     //
     // DJH: This call is expensive. Avoid calling when gammas don't differ.
     // Set flag in SimParameters and make this call conditional.
     // 
+    TIMER_START(patch->timerSet, RATTLE1);
     rattle1(dt_fs,1);  // conserve momentum if gammas differ
+    TIMER_STOP(patch->timerSet, RATTLE1);
 
     FullAtom *a = patch->atom.begin();
     int numAtoms = patch->numAtoms;
@@ -1445,6 +1607,7 @@ void Sequencer::langevinPiston(int step)
    int numAtoms = patch->numAtoms;
    // Blocking receive for the updated lattice scaling factor.
    Tensor factor = broadcast->positionRescaleFactor.get(step);
+   TIMER_START(patch->timerSet, PISTON);
    // JCP FIX THIS!!!
    Vector velFactor(1/factor.xx,1/factor.yy,1/factor.zz);
    patch->lattice.rescale(factor);
@@ -1508,6 +1671,7 @@ void Sequencer::langevinPiston(int step)
       a[i].velocity.z *= velFactor.z;
     }
    }
+   TIMER_STOP(patch->timerSet, PISTON);
   }
 }
 
@@ -1740,6 +1904,8 @@ void Sequencer::saveForce(const int ftag)
 void Sequencer::addForceToMomentum(
     BigReal timestep, const int ftag, const int useSaved
     ) {
+  NAMD_EVENT_RANGE_2(patch->flags.event_on,
+      NamdProfileEvent::ADD_FORCE_TO_MOMENTUM);
 #if CMK_BLUEGENEL
   CmiNetworkProgressAfter (0);
 #endif
@@ -1755,6 +1921,8 @@ void Sequencer::addForceToMomentum3(
     const BigReal timestep2, const int ftag2, const int useSaved2,
     const BigReal timestep3, const int ftag3, const int useSaved3
     ) {
+  NAMD_EVENT_RANGE_2(patch->flags.event_on,
+      NamdProfileEvent::ADD_FORCE_TO_MOMENTUM);
 #if CMK_BLUEGENEL
   CmiNetworkProgressAfter (0);
 #endif
@@ -1774,7 +1942,8 @@ void Sequencer::addForceToMomentum3(
 
 void Sequencer::addVelocityToPosition(BigReal timestep)
 {
-  RANGE("addVelocityToPosition", 5);
+  NAMD_EVENT_RANGE_2(patch->flags.event_on,
+      NamdProfileEvent::ADD_VELOCITY_TO_POSITION);
 #if CMK_BLUEGENEL
   CmiNetworkProgressAfter (0);
 #endif
@@ -1800,7 +1969,7 @@ void Sequencer::hardWallDrude(BigReal dt, int pressure)
 
 void Sequencer::rattle1(BigReal dt, int pressure)
 {
-  RANGE("rattle1", 9);
+  NAMD_EVENT_RANGE_2(patch->flags.event_on, NamdProfileEvent::RATTLE1);
   if ( simParams->rigidBonds != RIGID_NONE ) {
     Tensor virial;
     Tensor *vp = ( pressure ? &virial : 0 );
@@ -1810,7 +1979,34 @@ void Sequencer::rattle1(BigReal dt, int pressure)
       Node::Object()->enableEarlyExit();
       terminate();
     }
+#if 0
+    printf("virial = %g %g %g  %g %g %g  %g %g %g\n",
+        virial.xx, virial.xy, virial.xz,
+        virial.yx, virial.yy, virial.yz,
+        virial.zx, virial.zy, virial.zz);
+#endif
     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
+#if 0
+    {
+      const FullAtom *a = patch->atom.const_begin();
+      for (int n=0;  n < patch->numAtoms;  n++) {
+        printf("pos[%d] =  %g %g %g\n", n,
+            a[n].position.x, a[n].position.y, a[n].position.z);
+      }
+      for (int n=0;  n < patch->numAtoms;  n++) {
+        printf("vel[%d] =  %g %g %g\n", n,
+            a[n].velocity.x, a[n].velocity.y, a[n].velocity.z);
+      }
+      if (pressure) {
+        for (int n=0;  n < patch->numAtoms;  n++) {
+          printf("force[%d] =  %g %g %g\n", n,
+              patch->f[Results::normal][n].x,
+              patch->f[Results::normal][n].y,
+              patch->f[Results::normal][n].z);
+        }
+      }
+    }
+#endif
   }
 }
 
@@ -1830,7 +2026,7 @@ void Sequencer::rattle1(BigReal dt, int pressure)
 
 void Sequencer::maximumMove(BigReal timestep)
 {
-  RANGE("maximumMove", 4);
+  NAMD_EVENT_RANGE_2(patch->flags.event_on, NamdProfileEvent::MAXIMUM_MOVE);
 
   FullAtom *a = patch->atom.begin();
   int numAtoms = patch->numAtoms;
@@ -1886,7 +2082,7 @@ void Sequencer::minimizationQuenchVelocity(void)
 
 void Sequencer::submitHalfstep(int step)
 {
-  RANGE("submitHalfstep", 6);
+  NAMD_EVENT_RANGE_2(patch->flags.event_on, NamdProfileEvent::SUBMIT_HALFSTEP);
 
   // velocity-dependent quantities *** ONLY ***
   // positions are not at half-step when called
@@ -2065,8 +2261,11 @@ void Sequencer::calcFixVirial(Tensor& fixVirialNormal, Tensor& fixVirialNbond, T
 
 void Sequencer::submitReductions(int step)
 {
-  RANGE("submitReductions", 10);
+#ifndef UPPER_BOUND
+  NAMD_EVENT_RANGE_2(patch->flags.event_on,
+      NamdProfileEvent::SUBMIT_REDUCTIONS);
   FullAtom *a = patch->atom.begin();
+#endif
   int numAtoms = patch->numAtoms;
 
 #if CMK_BLUEGENEL
@@ -2076,6 +2275,7 @@ void Sequencer::submitReductions(int step)
   reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
   reduction->item(REDUCTION_MARGIN_VIOLATIONS) += patch->marginViolations;
 
+#ifndef UPPER_BOUND
   // For non-Multigrator doKineticEnergy = 1 always
   if (doKineticEnergy || doMomenta || patch->flags.doVirial)
   {
@@ -2294,9 +2494,12 @@ void Sequencer::submitReductions(int step)
       ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_SLOW,fixForceSlow);
     }
   }
+#endif // UPPER_BOUND
 
   reduction->submit();
+#ifndef UPPER_BOUND
   if (pressureProfileReduction) pressureProfileReduction->submit();
+#endif
 }
 
 void Sequencer::submitMinimizeReductions(int step, BigReal fmax2)
@@ -2446,7 +2649,8 @@ void Sequencer::submitCollections(int step, int zeroVel)
   //
   //patch->copy_updates_to_AOS();
 
-  RANGE("submitCollections", 11);
+  NAMD_EVENT_RANGE_2(patch->flags.event_on,
+      NamdProfileEvent::SUBMIT_COLLECTIONS);
   int prec = Output::coordinateNeeded(step);
   if ( prec ) {
     collection->submitPositions(step,patch->atom,patch->lattice,prec);
index 23f2b50..50ae162 100644 (file)
@@ -651,6 +651,23 @@ void SimParameters::config_parser_basic(ParseOptions &opts) {
 #endif
    opts.optional("main", "waterModel", "Water model to use", PARSE_STRING);
    opts.optionalB("main", "LJcorrection", "Apply analytical tail corrections for energy and virial", &LJcorrection, FALSE);
+#ifdef TIMER_COLLECTION
+   opts.optional("main", "TimerBinWidth",
+       "Bin width of timer histogram collection in microseconds",
+       &timerBinWidth, 1.0);
+#endif
+#if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED)
+   // default NVTX or Projections profiling is up to the first 1000 patches
+   opts.optional("main", "beginEventPatchID","Beginning patch ID for profiling",
+       &beginEventPatchID, 0);
+   opts.optional("main", "endEventPatchID", "Ending patch ID for profiling",
+       &endEventPatchID, 5000);
+   // default NVTX or Projections profiling is up to the first 1000 time steps
+   opts.optional("main", "beginEventStep", "Beginning time step for profiling",
+       &beginEventStep, 0);
+   opts.optional("main", "endEventStep", "Ending time step for profiling",
+       &endEventStep, 1000);
+#endif
 }
 
 void SimParameters::config_parser_fileio(ParseOptions &opts) {
@@ -4386,7 +4403,8 @@ void SimParameters::check_config(ParseOptions &opts, ConfigList *config, char *&
       }
     }
 #endif
-}
+} // check_config()
+
 
 void SimParameters::print_config(ParseOptions &opts, ConfigList *config, char *&cwd) {
 
index 94b8df0..d352171 100644 (file)
@@ -105,6 +105,17 @@ public:
 //  MAKE SURE THAT THIS CLASS CAN BE BIT COPIED OR YOU WILL HAVE TO
 //  ADD SPECIAL CODE TO send_SimParameters() and receive_SimParameters()
 
+#if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED)
+  int beginEventPatchID;
+  int endEventPatchID;
+  int beginEventStep;
+  int endEventStep;
+#endif
+
+#ifdef TIMER_COLLECTION
+  double timerBinWidth;  // default 1
+#endif
+
   Bool lonepairs;  // enable lone pairs
   int watmodel; // integer code for the water model in use
                 // choices are defined in common.h