(JEC) 9/4/15 timing function from Martin R.

Makefile

@@ -26,7 +26,10 @@ clean :
 #C++ common Objects
 CXXOBJ = $(OBJ)laguerreBuilder.o \
 	$(OBJ)laguerreTransform.o \
-	$(OBJ)lagSphericTransform.o
+	$(OBJ)lagSphericTransform.o \
+	$(OBJ)error_handling.o \
+	$(OBJ)walltimer.o
+
 
 #C++ common Headers
 CXXHDR = lagsht_exceptions.h \
@@ -34,7 +37,9 @@ CXXHDR = lagsht_exceptions.h \
 	lagsht_utils.h \
 	laguerreBuilder.h \
 	laguerreTransform.h \
-	lagSphericTransform.h
+	lagSphericTransform.h \
+	error_handling.h \
+	walltimer.h
 
 
 CXX = g++

README.md

@@ -12,7 +12,28 @@ git clone git@gitlab.in2p3.fr:campagne/LagSHT.git
 * extension of the computation of the weights (Laguerre function) up to N ~ 5500 by transfert of the exponential scale a the very beginning (see laguerreBuilder.cc)
 *  Start documentation (Doxygen)
 
-> doxygen doxydoc  
+`> doxygen doxydoc`  
+
+* Timing files (error_handling, walltimer) give nice Timing trees
+
+```
+Total wall clock time for 'MultiSphericalLaguerreTransform': 92.0287s
+|
++- processing                     : 92.05% (84.7126s)
+|  |
+|  +- processing part Analysis       : 51.18% (43.3588s)
+|  +- processing part Synthesis      : 48.82% (41.3539s)
+|  |  |
+|  |  +- SHT Synthesis                  : 68.20% (28.2026s)
+|  |  +- Laguerre MultiSynthesis        : 24.84% (10.2739s)
+|  |  +- <unaccounted>                  :  6.96% ( 2.8773s)
+|  |  
+|  +- <unaccounted>                  :  0.00% ( 0.0000s)
+|  
++- data input                     :  0.91% ( 0.8384s)
++- <unaccounted>                  :  7.04% ( 6.4777s)
+
+```
 
 # 8/4/15
 * creation d'un directory **Sophya** pour deplacer le code existant (v0.2) et pouvoir creer un code independant pour etre distribuable. 

error_handling.cc

+/*
+ *  This file is part of libcxxsupport.
+ *
+ *  libcxxsupport is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  libcxxsupport is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with libcxxsupport; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+/*
+ *  libcxxsupport is being developed at the Max-Planck-Institut fuer Astrophysik
+ *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
+ *  (DLR).
+ */
+
+/*
+ *  Utilities for error reporting
+ *
+ *  Copyright (C) 2003-2011 Max-Planck-Society
+ *  Author: Martin Reinecke
+ */
+
+#include "error_handling.h"
+
+using namespace std;
+
+PlanckError::PlanckError(const string &message) : msg (message) {}
+PlanckError::PlanckError(const char *message) : msg (message) {}
+
+//virtual
+PlanckError::~PlanckError() {}
+
+void planck_failure__(const char *file, int line, const char *func,
+  const string &msg)
+  {
+  cerr << "Error encountered at " << file << ", line " << line << endl;
+  if (func) cerr << "(function " << func << ")" << endl;
+  if (msg!="") cerr << endl << msg << endl;
+  cerr << endl;
+  }
+
+void planck_failure__(const char *file, int line, const char *func,
+  const char *msg)
+  { planck_failure__ (file,line,func,string(msg)); }
+
+void killjob__()
+  { throw; }

error_handling.h

+/*
+ *  This file is part of libcxxsupport.
+ *
+ *  libcxxsupport is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  libcxxsupport is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with libcxxsupport; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+/*
+ *  libcxxsupport is being developed at the Max-Planck-Institut fuer Astrophysik
+ *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
+ *  (DLR).
+ */
+
+/*! \file error_handling.h
+ *  Utilities for error reporting
+ *
+ *  Copyright (C) 2003-2011 Max-Planck-Society
+ *  Authors: Reinhard Hell, Martin Reinecke
+ */
+
+#ifndef PLANCK_ERROR_HANDLING_H
+#define PLANCK_ERROR_HANDLING_H
+
+#include <string>
+#include <iostream>
+
+#if defined (__GNUC__)
+#define PLANCK_FUNC_NAME__ __PRETTY_FUNCTION__
+#else
+#define PLANCK_FUNC_NAME__ 0
+#endif
+
+void planck_failure__(const char *file, int line, const char *func,
+  const std::string &msg);
+void planck_failure__(const char *file, int line, const char *func,
+  const char *msg);
+void killjob__();
+
+class PlanckError
+  {
+  private:
+    std::string msg;
+
+  public:
+    explicit PlanckError(const std::string &message);
+    explicit PlanckError(const char *message);
+
+    virtual const char* what() const
+      { return msg.c_str(); }
+
+    virtual ~PlanckError();
+  };
+
+/*! \defgroup errorgroup Error handling */
+/*! \{ */
+
+/*! Writes diagnostic output and exits with an error status. */
+#define planck_fail(msg) \
+do { planck_failure__(__FILE__,__LINE__,PLANCK_FUNC_NAME__,msg); \
+throw PlanckError(msg); } while(0)
+
+/*! Throws a PlanckError without diagnostic message. */
+#define planck_fail_quietly(msg) \
+do { throw PlanckError(msg); } while(0)
+
+/*! Writes diagnostic output and exits with an error status if \a testval
+    is \a false. */
+#define planck_assert(testval,msg) \
+do { if (testval); else planck_fail(msg); } while(0)
+
+/*! Macro for improving error diagnostics. Should be placed immediately
+    after the opening brace of \c main(). Must be used in conjunction with
+    \c PLANCK_DIAGNOSIS_END. */
+#define PLANCK_DIAGNOSIS_BEGIN try {
+/*! Macro for improving error diagnostics. Should be placed immediately
+    before the closing brace of \c main(). Must be used in conjunction with
+    \c PLANCK_DIAGNOSIS_BEGIN. */
+#define PLANCK_DIAGNOSIS_END \
+} \
+catch (PlanckError &) \
+  { killjob__(); /* no need for further diagnostics; they were shown already */ } \
+catch (std::exception &e) \
+  { std::cerr << "std::exception: " << e.what() << std::endl; killjob__(); } \
+catch (...) \
+  { std::cerr << "Unknown exception" << std::endl; killjob__(); }
+
+/*! \} */
+
+#endif

lagSphericTransform.cc

@@ -3,6 +3,8 @@
 
 #include "lagSphericTransform.h"
 #include "lagsht_exceptions.h"
+#include "walltimer.h"    //timing
+
 
 #define DEBUG 0
 
@@ -56,8 +58,10 @@ void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn,
   cout << "SphLagTransform::Synthesis type <r_8> I/O .... Go..." << endl;
   cout << "Start   Inverse Laguerre...." << endl;
 #endif
-
+  
+  tstack_push("Laguerre MultiSynthesis");
   lagTrans_.MultiSynthesis(flmn, flmk, Nalm_);
+  tstack_pop("Laguerre MultiSynthesis");
 
 
 #if DEBUG >= 1
@@ -65,6 +69,7 @@ void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn,
   cout << "Start Inverse H.Spherical Transform " << endl;
 #endif
 
+  tstack_push("SHT Synthesis");
   int off7n, off7k;
   //perform the synthesis per shell
   for(int n =0; n<N_; n++){
@@ -78,6 +83,9 @@ void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn,
     sharpjob_.alm2map(&(flmk.data()[off7n].real()),&fijk.data()[off7k],false);
 
   }//loop on shell
+
+    tstack_pop("SHT Synthesis");
+
 #if DEBUG >= 1
   cout << "SphLagTransform::Synthesis.... Done" << endl;
 #endif

lagsht_testsuite.cc

@@ -10,6 +10,7 @@ using namespace std;
 
 #include "lagsht_exceptions.h"
 #include "lagsht_utils.h" //getMemorySize
+#include "walltimer.h"    //timing
 
 
 #include "lagSphericTransform.h"
@@ -169,6 +170,8 @@ void MultiSphericalLaguerreTransform() {
 
   cout << " ___________  MultiSphericalLaguerreTransform   TEST _____________ " << endl;
 
+  tstack_push("data input");
+
   int state = 1234567 + 8912 ; //random seed
 
 
@@ -198,7 +201,10 @@ void MultiSphericalLaguerreTransform() {
   std::copy(flmnOrig.begin(), flmnOrig.end(), std::ostream_iterator<char>(std::cout, " "));
 #endif
       
-  //  Timer tm("Test 4 r_8 function");
+  
+  tstack_pop("data input");
+  tstack_push("processing");
+  tstack_push("processing part Synthesis");
 
 #if DEBUG >= 1
   cout << "Main:Synthesis r_8 function..." << endl;
@@ -210,16 +216,19 @@ void MultiSphericalLaguerreTransform() {
     cout << "fijk("<<i<<"): " << fijk(i) << endl;
   }
 #endif
-  //  tm.Split("SHTL r_8/end synthesis");
-
-
+  
+  tstack_pop("processing part Synthesis");
+  tstack_push("processing part Analysis");
+  
 #if DEBUG >= 1
   cout << "Main:Analysis r_8 function..." << endl;
 #endif
   vector< complex<r_8> > flmn(Ntot);
   sphlagtrans.Analysis(fijk,flmn);
     
-  //  tm.Split("SHTL r_8/end analysis");
+  tstack_pop("processing part Analysis");
+  tstack_pop("processing");  
+
 
 
   cout << "Error analysis..." << endl;
@@ -324,7 +333,11 @@ int main(int narg, char *arg[]) {
       MultiLaguerreTransformSdtAlone();
       break;
     case 3:
+      tstack_push("MultiSphericalLaguerreTransform");
       MultiSphericalLaguerreTransform();
+      tstack_pop("MultiSphericalLaguerreTransform");
+      tstack_report("MultiSphericalLaguerreTransform");
+  
       break;
     default:
       throw LagSHTError("EROOR Test type unkwown");

walltimer.cc

+/*
+ *  This file is part of libcxxsupport.
+ *
+ *  libcxxsupport is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  libcxxsupport is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with libcxxsupport; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+/*
+ *  libcxxsupport is being developed at the Max-Planck-Institut fuer Astrophysik
+ *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
+ *  (DLR).
+ */
+
+/*
+ *  This file contains functionality related to wall-clock timers
+ *
+ *  Copyright (C) 2010, 2011, 2012 Max-Planck-Society
+ *  Author: Martin Reinecke
+ */
+
+#include <iostream>
+#include <utility>
+#include <cstdio>
+#include <cmath>
+#include <algorithm>
+#include "walltimer.h"
+#include "walltime_c.h"
+#include "error_handling.h"
+
+using namespace std;
+
+void wallTimer::start()
+  { start(wallTime()); }
+void wallTimer::stop()
+  { stop(wallTime()); }
+double wallTimer::acc() const
+  { return acc(wallTime()); }
+
+int wallTimerSet::getIndex(const string &name)
+  {
+  maptype::const_iterator it = lut.find(name);
+  if (it!=lut.end())
+    return it->second;
+  timer.push_back(wallTimer());
+  lut[name]=timer.size()-1;
+  return timer.size()-1;
+  }
+
+void wallTimerSet::start(int index)
+  { timer[index].start(); }
+void wallTimerSet::stop(int index)
+  { timer[index].stop(); }
+void wallTimerSet::stopstart(int index1, int index2)
+  { double t=wallTime(); timer[index1].stop(t); timer[index2].start(t); }
+void wallTimerSet::reset(int index)
+  { timer[index].reset(); }
+double wallTimerSet::acc(int index)
+  { return timer[index].acc(); }
+void wallTimerSet::start(const string &name)
+  { start(getIndex(name)); }
+void wallTimerSet::stop(const string &name)
+  { stop(getIndex(name)); }
+void wallTimerSet::stopstart(const string &name1, const string &name2)
+  { stopstart(getIndex(name1),getIndex(name2)); }
+void wallTimerSet::reset(const string &name)
+  { reset(getIndex(name)); }
+double wallTimerSet::acc(const string &name)
+  { return acc(getIndex(name)); }
+
+void wallTimerSet::report() const
+  {
+  cout << "\nWall clock timer report:" << endl;
+  for (maptype::const_iterator it=lut.begin(); it!=lut.end(); ++it)
+    printf("  %-15s: %10.5fs\n", it->first.c_str(), timer[it->second].acc());
+  cout << "End wall clock timer report\n" << endl;
+  }
+
+wallTimerSet wallTimers;
+
+namespace {
+
+class tstack_node;
+
+typedef map<string,tstack_node>::iterator Ti;
+typedef map<string,tstack_node>::const_iterator Tci;
+typedef pair<Tci,double> Tipair;
+
+class tstack_node
+  {
+  public:
+    tstack_node *parent;
+    wallTimer wt;
+    string name;
+    map<string,tstack_node> child;
+
+    tstack_node(const string &name_, tstack_node *parent_)
+      : parent(parent_), name(name_) {}
+
+    int max_namelen() const
+      {
+      int res=name.length();
+      for (Tci it=child.begin(); it!=child.end(); ++it)
+        res=max(res,it->second.max_namelen());
+      return res;
+      }
+  };
+
+tstack_node tstack_root("root",0);
+tstack_node *curnode=0;
+double overhead=0.;
+
+struct timecomp
+  {
+  bool operator() (const Tipair &a, const Tipair &b) const
+    { return a.second>b.second; }
+  };
+
+void tstack_report(const tstack_node &node, const string &indent, int twidth,
+  int slen)
+  {
+  double total=node.wt.acc();
+  vector<Tipair> tmp;
+  for (Tci it=node.child.begin(); it!=node.child.end(); ++it)
+    tmp.push_back(make_pair(it,it->second.wt.acc()));
+
+  if (tmp.size()>0)
+    {
+    sort(tmp.begin(),tmp.end(),timecomp());
+    double tsum=0;
+    printf("%s|\n", indent.c_str());
+    for (unsigned i=0; i<tmp.size(); ++i)
+      {
+      printf("%s+- %-*s:%6.2f%% (%*.4fs)\n",indent.c_str(),slen,
+        (tmp[i].first->first).c_str(), 100*tmp[i].second/total,twidth,
+        tmp[i].second);
+      tstack_report(tmp[i].first->second,indent+"|  ",twidth,slen);
+      tsum+=tmp[i].second;
+      }
+    printf("%s+- %-*s:%6.2f%% (%*.4fs)\n%s\n",indent.c_str(),slen,
+      "<unaccounted>",100*(total-tsum)/total,twidth,total-tsum,indent.c_str());
+    }
+  }
+
+} // unnamed namespace
+
+void tstack_push(const string &name)
+  {
+  double t0=wallTime();
+  if (curnode==0) curnode=&tstack_root;
+  Ti it=curnode->child.find(name);
+  if (it==curnode->child.end())
+    it=curnode->child.insert (make_pair(name,tstack_node(name,curnode))).first;
+  curnode=&(it->second);
+  double t1=wallTime();
+  curnode->wt.start(0.5*(t0+t1));
+  overhead+=t1-t0;
+  }
+void tstack_pop(const string &name)
+  {
+  double t0=wallTime();
+  planck_assert(curnode && (curnode->name==name), "invalid tstack operation");
+  double t1=wallTime();
+  curnode->wt.stop(0.5*(t0+t1));
+  curnode=curnode->parent;
+  overhead+=t1-t0;
+  }
+void tstack_pop()
+  {
+  double t0=wallTime();
+  planck_assert(curnode, "invalid tstack operation");
+  double t1=wallTime();
+  curnode->wt.stop(0.5*(t0+t1));
+  curnode=curnode->parent;
+  overhead+=t1-t0;
+  }
+void tstack_replace(const string &name2)
+  {
+  double t0=wallTime();
+  planck_assert(curnode, "invalid tstack operation");
+  tstack_node *savenode=curnode;
+  curnode=curnode->parent;
+  Ti it=curnode->child.find(name2);
+  if (it==curnode->child.end())
+    it=curnode->child.insert(make_pair(name2,tstack_node(name2,curnode))).first;
+  curnode=&(it->second);
+  double t1=wallTime();
+  double t=0.5*(t0+t1);
+  savenode->wt.stop(t);
+  curnode->wt.start(t);
+  overhead+=t1-t0;
+  }
+void tstack_replace(const string &name1, const string &name2)
+  {
+  planck_assert(curnode && (curnode->name==name1), "invalid tstack operation");
+  tstack_replace(name2);
+  }
+
+void tstack_report(const string &stem)
+  {
+  const tstack_node *ptr = 0;
+  for (Tci it=tstack_root.child.begin(); it!=tstack_root.child.end(); ++it)
+    if (it->first==stem) ptr=&(it->second);
+  planck_assert(ptr,"invalid stem");
+  int slen=string("<unaccounted>").size();
+  slen = max(slen,ptr->max_namelen());
+
+  double total=ptr->wt.acc();
+  printf("\nTotal wall clock time for '%s': %1.4fs\n",stem.c_str(),total);
+
+  int logtime=max(1,int(log10(total)+1));
+  tstack_report(*ptr,"",logtime+5,slen);
+
+  printf("\nAccumulated timing overhead: approx. %1.4fs\n",overhead);
+  }

walltimer.h

+/*
+ *  This file is part of libcxxsupport.
+ *
+ *  libcxxsupport is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  libcxxsupport is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with libcxxsupport; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+/*
+ *  libcxxsupport is being developed at the Max-Planck-Institut fuer Astrophysik
+ *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
+ *  (DLR).
+ */
+
+/*! \file walltimer.h
+ *  Functionality related to wall-clock timers
+ *
+ *  Copyright (C) 2010, 2011 Max-Planck-Society
+ *  \author Martin Reinecke
+ */
+
+#ifndef PLANCK_WALLTIMER_H
+#define PLANCK_WALLTIMER_H
+
+#include <string>
+#include <map>
+#include <vector>
+
+class wallTimer
+  {
+  private:
+    double t_acc, t_started;
+    bool running;
+
+  public:
+    wallTimer() : t_acc(0.), t_started(0.), running(false) {}
+    void start(double wtime_now)
+      { if (!running) { t_started=wtime_now; running=true; } }
+    void start();
+    void stop(double wtime_now)
+      { if (running) { t_acc+=wtime_now-t_started; running=false; } }
+    void stop();
+    void reset() { t_acc=t_started=0.; running=false;}
+    double acc(double wtime_now) const
+      { return running ? t_acc+wtime_now-t_started : t_acc; }
+    double acc() const;
+  };
+
+class wallTimerSet
+  {
+  public:
+    typedef std::map<std::string,int> maptype;
+
+  private:
+    maptype lut;
+    std::vector<wallTimer> timer;
+
+  public:
+    int getIndex(const std::string &name);
+    void start(int index);
+    void stop(int index);
+    void stopstart(int index1, int index2);
+    void reset(int index);
+    double acc(int index);
+    void start(const std::string &name);
+    void stop(const std::string &name);
+    void stopstart(const std::string &name1, const std::string &name2);
+    void reset(const std::string &name);
+    double acc(const std::string &name);
+
+    void report() const;
+
+    const maptype &table() const { return lut; }
+  };
+
+extern wallTimerSet wallTimers;
+
+void tstack_push(const std::string &name);
+void tstack_pop(const std::string &name);
+void tstack_pop();
+void tstack_replace(const std::string &name1, const std::string &name2);
+void tstack_replace(const std::string &name);
+void tstack_report(const std::string &stem);
+