CarpetLib, CycleClock: Move cycle counting clock into its own thorn

New thorn CycleClock, mostly taken out of CarpetLib. This thorn provides an (almost) cycle-accurate clock, taken from FFTW. This clock is both directly accessible, and is also wrapped in a Cactus clock.

Note: Applying this commit will require adding thorn CycleClock to thorn lists and parameter files.

14 files changed, 353 insertions, 222 deletions

diff --git a/Carpet/CarpetLib/configuration.ccl b/Carpet/CarpetLib/configuration.ccl
index a77f237fc..12ba9306c 100644
--- a/Carpet/CarpetLib/configuration.ccl
+++ b/Carpet/CarpetLib/configuration.ccl
@@ -6,7 +6,7 @@ PROVIDES CarpetLib
   LANG
 }
 
-REQUIRES MPI Vectors
+REQUIRES CycleClock MPI Vectors
 
 OPTIONAL LoopControl MPI
 {
diff --git a/Carpet/CarpetLib/interface.ccl b/Carpet/CarpetLib/interface.ccl
index 6495fdec3..0a2086c4f 100644
--- a/Carpet/CarpetLib/interface.ccl
+++ b/Carpet/CarpetLib/interface.ccl
@@ -2,6 +2,8 @@
 
 IMPLEMENTS: CarpetLib
 
+INHERITS: CycleClock
+
 includes header: mpi_string.hh in mpi_string.hh
 
 includes header: cacheinfo.hh in cacheinfo.hh
@@ -37,6 +39,7 @@ includes header: operators.hh in operators.hh
 uses include header: nompi.h
 
 uses include header: CarpetTimers.hh
+uses include header: cycleclock.h
 uses include header: loopcontrol.h
 uses include header: vectors.h
 
diff --git a/Carpet/CarpetLib/schedule.ccl b/Carpet/CarpetLib/schedule.ccl
index e719fc071..bf7dd05fc 100644
--- a/Carpet/CarpetLib/schedule.ccl
+++ b/Carpet/CarpetLib/schedule.ccl
@@ -1,11 +1,5 @@
 # Schedule definitions for thorn CarpetLib
 
-SCHEDULE CarpetLib_registercycleclock AT startup BEFORE Driver_Startup
-{
-  LANG: C
-  OPTIONS: global
-} "Register cycle based timer"
-
 SCHEDULE CarpetLib_test_prolongate_3d_rf2 AT paramcheck
 {
   LANG: C
diff --git a/Carpet/CarpetLib/src/timestat.cc b/Carpet/CarpetLib/src/timestat.cc
index e35fc3b43..08c02d270 100644
--- a/Carpet/CarpetLib/src/timestat.cc
+++ b/Carpet/CarpetLib/src/timestat.cc
@@ -34,53 +34,6 @@ namespace CarpetLib {
   
   
   
-  static
-  bool have_cputick = false;
-  
-  // CPU tick time in seconts
-  static
-  double cputick = 0.0;
-  
-  static
-  void
-  calculate_cputick ()
-  {
-    // Make a few warm-up measurements
-    getticks ();
-    getticks ();
-    getticks ();
-    
-#if 0
-    // Use usleep to calibrate the timer
-    for (int i=0; i<10; ++i) {
-      useconds_t const waittime = 1000 * 1000;
-      ticks const rstart = getticks ();
-      int const ierr = usleep (waittime);
-      ticks const rend = getticks ();
-      cputick = waittime / 1.0e6 / elapsed (rend, rstart);
-      if (not ierr) goto done;
-    }
-    CCTK_WARN (1, "Could not determine timer resolution");
-  done:;
-#endif
-    
-#if 1
-    // Use MPI_Wtime to calibrate the timer
-    ticks const rstart = getticks ();
-    double const wstart = MPI_Wtime ();
-    while (MPI_Wtime() < wstart + 1.0) {
-      // do nothing, just wait
-    }
-    ticks const rend = getticks ();
-    double const wend = MPI_Wtime ();
-    cputick = (wend - wstart) / elapsed (rend, rstart);
-#endif
-    
-    have_cputick = true;
-  }
-  
-  
-  
   // Call a timer
   static
   ticks
@@ -152,7 +105,6 @@ namespace CarpetLib {
     : timername (timername_)
   {
     assert (timername_);
-    if (not have_cputick) calculate_cputick ();
     resetstats ();
     if (not timerSet) timerSet = new TimerSet;
     timerSet->add (this);
@@ -280,11 +232,11 @@ namespace CarpetLib {
     
     os << timername << ":"
        << " cnt: " << count
-       << "   time: sum: " << cputick * wtime
-       << " avg: " << cputick * avg
-       << " stddev: " << cputick * stddev
-       << " min: " << cputick * wmin
-       << " max: " << cputick * wmax
+       << "   time: sum: " << seconds_per_tick() * wtime
+       << " avg: " << seconds_per_tick() * avg
+       << " stddev: " << seconds_per_tick() * stddev
+       << " min: " << seconds_per_tick() * wmin
+       << " max: " << seconds_per_tick() * wmax
        << "   bytes: sum: " << bytes
        << " avg: " << bavg
        << " stddev: " << bstddev
@@ -416,161 +368,4 @@ namespace CarpetLib {
     
   }
   
-  
-  
-  struct t_cycleclock {
-    double total;
-    double total_squared;
-    double min_total;
-    double max_total;
-    double count;
-    ticks last;
-    
-    t_cycleclock ()
-    {
-      reset();
-    }
-    
-    ~t_cycleclock ()
-    {
-    }
-    
-    void start ()
-    {
-      last = getticks();
-    }
-    
-    void stop ()
-    {
-      ticks const current = getticks();
-      double const difference = elapsed (current, last);
-      total += difference;
-      total_squared += pow (difference, 2);
-      min_total = min_total == 0.0 ? difference : min (min_total, difference);
-      max_total = max (min_total, difference);
-      count += 1;
-    }
-    
-    void reset ()
-    {
-      total         = 0.0;
-      total_squared = 0.0;
-      min_total     = 0.0;      // numeric_limits<double>::max();
-      max_total     = 0.0;
-      count         = 0.0;
-      // last          = 0.0;
-    }
-    
-  };
-  
-  
-  
-  void * cycleclock_create (int const timernum)
-  {
-    return new t_cycleclock;
-  }
-  
-  void cycleclock_destroy (int const timernum, void * const data)
-  {
-    if (data) {
-      delete static_cast<t_cycleclock*> (data);
-    }
-  }
-  
-  void cycleclock_start (int const timernum, void * const data)
-  {
-    static_cast<t_cycleclock*> (data) -> start();
-  }
-  
-  void cycleclock_stop (int const timernum, void * const data)
-  {
-    static_cast<t_cycleclock*> (data) -> stop();
-  }
-  
-  void cycleclock_reset (int const timernum, void * const data)
-  {
-    static_cast<t_cycleclock*> (data) -> reset();
-  }
-  
-  void cycleclock_get (int const timernum, void * const data_,
-                       cTimerVal * const vals)
-  {
-    t_cycleclock const & data = * static_cast<t_cycleclock const *> (data_);
-    
-    // Total time
-    vals[0].type       = val_double;
-    vals[0].heading    = "cycle";
-    vals[0].units      = "secs";
-    vals[0].val.d      = data.total;
-    vals[0].seconds    = cputick * vals[0].val.d;
-    vals[0].resolution = cputick;
-    
-    // Average
-    vals[1].type       = val_double;
-    vals[1].heading    = "cycle[avg]";
-    vals[1].units      = "secs";
-    vals[1].val.d      = data.count == 0.0 ? 0.0 : data.total / data.count;
-    vals[1].seconds    = cputick * vals[1].val.d;
-    vals[1].resolution = cputick;
-    
-    // Standard deviation
-    vals[2].type       = val_double;
-    vals[2].heading    = "cycle[stddev]";
-    vals[2].units      = "secs";
-    vals[2].val.d      = (data.count == 0.0
-                          ? 0.0
-                          : sqrt (abs (data.total_squared * data.count -
-                                       pow (data.total, 2)) / data.count));
-    vals[2].seconds    = cputick * vals[2].val.d;
-    vals[2].resolution = cputick;
-    
-    // Minimum
-    vals[3].type       = val_double;
-    vals[3].heading    = "cycle[min]";
-    vals[3].units      = "secs";
-    vals[3].val.d      = data.min_total;
-    vals[3].seconds    = cputick * vals[3].val.d;
-    vals[3].resolution = cputick;
-    
-    // Maximum
-    vals[4].type       = val_double;
-    vals[4].heading    = "cycle[max]";
-    vals[4].units      = "secs";
-    vals[4].val.d      = data.max_total;
-    vals[4].seconds    = cputick * vals[4].val.d;
-    vals[4].resolution = cputick;
-  }
-  
-  void cycleclock_set (int const timernum, void * const data_,
-                       cTimerVal * const vals)
-  {
-    t_cycleclock & data = * static_cast<t_cycleclock *> (data_);
-    
-    data.reset();               // punt
-    data.total = vals[0].val.d;
-  }
-  
-  extern "C" {
-    int CarpetLib_registercycleclock (void);
-  }
-  
-  int CarpetLib_registercycleclock (void)
-  {
-    if (not have_cputick) calculate_cputick ();
-    
-    cClockFuncs functions;
-    functions.n_vals  = 5;
-    functions.create  = cycleclock_create;
-    functions.destroy = cycleclock_destroy;
-    functions.start   = cycleclock_start;
-    functions.stop    = cycleclock_stop;
-    functions.reset   = cycleclock_reset;
-    functions.get     = cycleclock_get;
-    functions.set     = cycleclock_set;
-    
-    CCTK_ClockRegister("cycle", &functions);
-    
-    return 0;
-  }
-  
 } // namespace CarpetLib
diff --git a/Carpet/CarpetLib/src/timestat.hh b/Carpet/CarpetLib/src/timestat.hh
index e96d6aada..a773c4d2e 100644
--- a/Carpet/CarpetLib/src/timestat.hh
+++ b/Carpet/CarpetLib/src/timestat.hh
@@ -7,9 +7,7 @@
 
 #include <cctk.h>
 
-extern "C" {
-#include "cycle.h"
-}
+#include <cycleclock.h>
 
 
 
diff --git a/Carpet/CycleClock/configuration.ccl b/Carpet/CycleClock/configuration.ccl
new file mode 100644
index 000000000..a79586e82
--- /dev/null
+++ b/Carpet/CycleClock/configuration.ccl
@@ -0,0 +1,5 @@
+# Configuration definitions for thorn CycleClock
+
+PROVIDES CycleClock
+{
+}
diff --git a/Carpet/CycleClock/interface.ccl b/Carpet/CycleClock/interface.ccl
new file mode 100644
index 000000000..c52a14e8f
--- /dev/null
+++ b/Carpet/CycleClock/interface.ccl
@@ -0,0 +1,5 @@
+# Interface definition for thorn CycleClock
+
+IMPLEMENTS: CycleClock
+
+INCLUDES HEADER: cycleclock.h IN cycleclock.h
diff --git a/Carpet/CycleClock/param.ccl b/Carpet/CycleClock/param.ccl
new file mode 100644
index 000000000..887db2384
--- /dev/null
+++ b/Carpet/CycleClock/param.ccl
@@ -0,0 +1 @@
+# Parameter definitions for thorn CycleClock
diff --git a/Carpet/CycleClock/schedule.ccl b/Carpet/CycleClock/schedule.ccl
new file mode 100644
index 000000000..cde207d25
--- /dev/null
+++ b/Carpet/CycleClock/schedule.ccl
@@ -0,0 +1,6 @@
+# Schedule definitions for thorn CycleClock
+
+SCHEDULE CycleClock_Setup AT startup BEFORE Driver_Startup
+{
+  LANG: C
+} "Set up CycleClock"
diff --git a/Carpet/CycleClock/src/clock.cc b/Carpet/CycleClock/src/clock.cc
new file mode 100644
index 000000000..26e251ca1
--- /dev/null
+++ b/Carpet/CycleClock/src/clock.cc
@@ -0,0 +1,177 @@
+#include "cycleclock.h"
+
+#include <cctk.h>
+
+#include <cmath>
+
+
+
+namespace CycleClock {
+  
+  using namespace std;
+  
+  
+  
+  class cycleclock_t {
+    double sum;
+    double sum2;
+    double min;
+    double max;
+    double count;
+    
+    ticks last;
+    
+  public:
+    cycleclock_t()
+    {
+      reset();
+    }
+    
+    ~cycleclock_t()
+    {
+    }
+    
+    void start()
+    {
+      last = getticks();
+    }
+    
+    void stop()
+    {
+      ticks const current = getticks();
+      double const difference = elapsed(current, last);
+      sum += difference;
+      sum2 += pow(difference, 2.0);
+      min = min == 0.0 ? difference : fmin(min, difference);
+      max = fmax(min, difference);
+      count += 1.0;
+    }
+    
+    void reset()
+    {
+      sum   = 0.0;
+      sum2  = 0.0;
+      min   = 0.0;              // numeric_limits<double>::max();
+      max   = 0.0;
+      count = 0.0;
+    }
+    
+    void get(cTimerVal* restrict const vals) const
+    {
+      double const tick = seconds_per_tick();
+      
+      // Sum time
+      vals[0].type       = val_double;
+      vals[0].heading    = "cycle";
+      vals[0].units      = "secs";
+      vals[0].val.d      = sum;
+      vals[0].seconds    = tick * vals[0].val.d;
+      vals[0].resolution = tick;
+      
+      // Average
+      vals[1].type       = val_double;
+      vals[1].heading    = "cycle[avg]";
+      vals[1].units      = "secs";
+      vals[1].val.d      = count == 0.0 ? 0.0 : sum / count;
+      vals[1].seconds    = tick * vals[1].val.d;
+      vals[1].resolution = tick;
+      
+      // Standard deviation
+      vals[2].type       = val_double;
+      vals[2].heading    = "cycle[sdv]";
+      vals[2].units      = "secs";
+      vals[2].val.d      = (count == 0.0 ?
+                            0.0 :
+                            sqrt(fabs(sum2 * count - pow(sum, 2.0)) / count));
+      vals[2].seconds    = tick * vals[2].val.d;
+      vals[2].resolution = tick;
+      
+      // Minimum
+      vals[3].type       = val_double;
+      vals[3].heading    = "cycle[min]";
+      vals[3].units      = "secs";
+      vals[3].val.d      = min;
+      vals[3].seconds    = tick * vals[3].val.d;
+      vals[3].resolution = tick;
+      
+      // Maximum
+      vals[4].type       = val_double;
+      vals[4].heading    = "cycle[max]";
+      vals[4].units      = "secs";
+      vals[4].val.d      = max;
+      vals[4].seconds    = tick * vals[4].val.d;
+      vals[4].resolution = tick;
+    }
+    
+    void set(cTimerVal const* restrict const vals)
+    {
+      reset();                  // punt
+      sum = vals[0].val.d;
+    }
+  };
+  
+  
+  
+  void* cycleclock_create(int const timernum)
+  {
+    return new cycleclock_t;
+  }
+  
+  void cycleclock_destroy(int const timernum, void* const data)
+  {
+    if (!data) return;
+    delete static_cast<cycleclock_t*>(data);
+  }
+  
+  void cycleclock_start(int const timernum, void* const data)
+  {
+    static_cast<cycleclock_t*>(data)->start();
+  }
+  
+  void cycleclock_stop(int const timernum, void* const data)
+  {
+    static_cast<cycleclock_t*>(data)->stop();
+  }
+  
+  void cycleclock_reset(int const timernum, void* const data)
+  {
+    static_cast<cycleclock_t*>(data)->reset();
+  }
+  
+  void cycleclock_get(int const timernum, void* const data,
+                      cTimerVal* const vals)
+  {
+    static_cast<cycleclock_t const*>(data)->get(vals);
+  }
+  
+  void cycleclock_set(int const timernum, void* const data,
+                      cTimerVal* const vals)
+  {
+    static_cast<cycleclock_t*>(data)->set(vals);
+  }
+  
+  void cycleclock_register()
+  {
+    cClockFuncs functions;
+    functions.n_vals  = 5;
+    functions.create  = cycleclock_create;
+    functions.destroy = cycleclock_destroy;
+    functions.start   = cycleclock_start;
+    functions.stop    = cycleclock_stop;
+    functions.reset   = cycleclock_reset;
+    functions.get     = cycleclock_get;
+    functions.set     = cycleclock_set;
+    CCTK_ClockRegister("cycle", &functions);
+  }
+  
+  
+  
+  extern "C"
+  int CycleClock_Setup()
+  {
+    measure_tick();
+    cycleclock_register();
+    return 0;
+  }
+  
+} // namespace CycleClock
diff --git a/Carpet/CycleClock/src/cycle.cc b/Carpet/CycleClock/src/cycle.cc
new file mode 100644
index 000000000..3a6b48364
--- /dev/null
+++ b/Carpet/CycleClock/src/cycle.cc
@@ -0,0 +1,108 @@
+#include "cycleclock.h"
+
+#include <cctk.h>
+
+
+
+// Find a good wall clock timer
+#ifdef _OPENMP
+#  include <omp.h>
+#endif
+#ifdef HAVE_CAPABILITY_MPI
+#  include <mpi.h>
+#endif
+#ifdef HAVE_SYS_TIME_H
+#  include <sys/time.h>
+#endif
+
+
+
+static double cycleclock_tick = -1.0; // uninitialised
+
+
+
+#ifdef HAVE_SYS_TIME_H
+static double get_sys_time()
+{
+  timeval tp;
+  gettimeofday(&tp, NULL);
+  return tp.tv_sec + 1.0e-6 * tp.tv_usec;
+}
+#endif
+
+
+
+void measure_tick()
+{
+  // Make a few warm-up measurements
+  getticks();
+  getticks();
+  getticks();
+  
+#ifdef _OPENMP
+  if (cycleclock_tick < 0.0) {
+    // Use omp_get_sys_time to calibrate the timer
+    CCTK_INFO("Measuring CycleClock tick via OpenMP...");
+    ticks const rstart = getticks();
+    double const wstart = omp_get_wtime();
+    while (omp_get_wtime() < wstart + 1.0) {
+      // do nothing, just wait
+    }
+    ticks const rend = getticks();
+    double const wend = omp_get_wtime();
+    cycleclock_tick = (wend - wstart) / elapsed(rend, rstart);
+  }
+#endif
+  
+#ifdef HAVE_CAPABILITY_MPI
+  if (cycleclock_tick < 0.0) {
+    // Use MPI_Wtime to calibrate the timer
+    CCTK_INFO("Measuring CycleClock tick via MPI...");
+    ticks const rstart = getticks();
+    double const wstart = MPI_Wtime();
+    while (MPI_Wtime() < wstart + 1.0) {
+        // do nothing, just wait
+    }
+    ticks const rend = getticks();
+    double const wend = MPI_Wtime();
+    cycleclock_tick = (wend - wstart) / elapsed(rend, rstart);
+  }
+#endif
+  
+#ifdef HAVE_SYS_TIME_H
+  if (cycleclock_tick < 0.0) {
+    // Use gettimeofday to calibrate the timer
+    CCTK_INFO("Measuring CycleClock tick via gettimeofday...");
+    ticks const rstart = getticks();
+    double const wstart = get_sys_time();
+    while (get_sys_time() < wstart + 1.0) {
+      // do nothing, just wait
+    }
+    ticks const rend = getticks();
+    double const wend = get_sys_time();
+    cycleclock_tick = (wend - wstart) / elapsed(rend, rstart);
+  }
+#endif
+  
+  if (cycleclock_tick < 0.0) {
+    // Give up; just set the time scale to 1
+    CCTK_INFO("Could not measure CycleClock tick");
+    cycleclock_tick = 1.0;
+  }
+  
+  CCTK_VInfo(CCTK_THORNSTRING,
+             "Calibrated CycleClock: %g ns per clock tick (%g GHz)",
+             1.0e9 * cycleclock_tick, 1.0e-9 / cycleclock_tick);
+}
+
+
+
+double seconds_per_tick()
+{
+  // if (CCTK_BUILTIN_EXPECT(cycleclock_tick < 0.0, false)) {
+  //   CCTK_WARN(CCTK_WARN_ALERT,
+  //             "Thorn CycleClock has not been activated; the first measurements may be wrong");
+  //   measure_tick();
+  // }
+  return cycleclock_tick;
+}
diff --git a/Carpet/CarpetLib/src/cycle.h b/Carpet/CycleClock/src/cycle.h
index d80f2ede8..cbcf67f43 100644
--- a/Carpet/CarpetLib/src/cycle.h
+++ b/Carpet/CycleClock/src/cycle.h
@@ -222,9 +222,9 @@ typedef unsigned long long ticks;
 
 static __inline__ ticks getticks(void)
 {
-     unsigned a, d; 
+     ticks a, d;
      asm volatile("rdtsc" : "=a" (a), "=d" (d)); 
-     return ((ticks)a) | (((ticks)d) << 32); 
+     return a | (d << 32); 
 }
 
 INLINE_ELAPSED(__inline__)
diff --git a/Carpet/CycleClock/src/cycleclock.h b/Carpet/CycleClock/src/cycleclock.h
new file mode 100644
index 000000000..c985609f1
--- /dev/null
+++ b/Carpet/CycleClock/src/cycleclock.h
@@ -0,0 +1,32 @@
+#ifndef CYCLECLOCK_H
+#define CYCLECLOCK_H
+
+/* This defines:
+ *    typedef XXX ticks;
+ *    ticks getticks();
+ *    double elapsed(ticks t1, ticks t0);
+ * Use as:
+ *    #include <cycleclock.h>
+ *    ticks t0 = getticks();
+ *    ...
+ *    ticks t1 = getticks();
+ *    double elapsed_ticks = elapsed(t1, t0);
+ *    double elapsed_seconds = seconds_per_tick() * elapsed_ticks;
+*/
+
+//#include <cctk.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+  
+#include "cycle.h"
+  
+  double seconds_per_tick();
+  void measure_tick();
+  
+#ifdef __cplusplus
+}
+#endif
+
+#endif  /* CYCLECLOCK_H */
diff --git a/Carpet/CycleClock/src/make.code.defn b/Carpet/CycleClock/src/make.code.defn
new file mode 100644
index 000000000..eb3a889e6
--- /dev/null
+++ b/Carpet/CycleClock/src/make.code.defn
@@ -0,0 +1,7 @@
+# Main make.code.defn file for thorn CycleClock
+
+# Source files in this directory
+SRCS = clock.cc cycle.cc
+
+# Subdirectories containing source files
+SUBDIRS =