Checking in fixed adaptive WilsonFlow

.gitignore

@@ -88,7 +88,6 @@ Thumbs.db
 # build directory #
 ###################
 build*/*
-Documentation/_build
 
 # IDE related files #
 #####################

Grid/qcd/smearing/WilsonFlow.h

@@ -7,6 +7,7 @@ Source file: ./lib/qcd/modules/plaquette.h
 Copyright (C) 2017
 
 Author: Guido Cossu  <guido.cossu@ed.ac.uk>
+Author: Chulwoo Jung <chulwoo@bnl.gov>
 
 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
@@ -35,7 +36,7 @@ template <class Gimpl>
 class WilsonFlow: public Smear<Gimpl>{
   unsigned int Nstep;
   unsigned int measure_interval;
-  mutable RealD epsilon, taus;
+  mutable RealD epsilon, taus,tolerance;
 
 
   mutable WilsonGaugeAction<Gimpl> SG;
@@ -47,13 +48,15 @@ class WilsonFlow: public Smear<Gimpl>{
 public:
   INHERIT_GIMPL_TYPES(Gimpl)
 
-  explicit WilsonFlow(unsigned int Nstep, RealD epsilon, unsigned int interval = 1):
+  explicit WilsonFlow(unsigned int Nstep, RealD epsilon, unsigned int interval = 1, RealD tol = 1e-3):
   Nstep(Nstep),
     epsilon(epsilon),
+    tolerance(tol),
     measure_interval(interval),
     SG(WilsonGaugeAction<Gimpl>(3.0)) {
     // WilsonGaugeAction with beta 3.0
     assert(epsilon > 0.0);
+    assert(tolerance > 0.0);
     LogMessage();
   }
 
@@ -64,6 +67,8 @@ public:
 	      << "[WilsonFlow] epsilon : " << epsilon << std::endl;
     std::cout << GridLogMessage
 	      << "[WilsonFlow] full trajectory : " << Nstep * epsilon << std::endl;
+    std::cout << GridLogMessage
+	      << "[WilsonFlow] tolerance : " << tolerance << std::endl;
   }
 
   virtual void smear(GaugeField&, const GaugeField&) const;
@@ -106,11 +111,14 @@ void WilsonFlow<Gimpl>::evolve_step_adaptive(typename Gimpl::GaugeField &U, Real
   if (maxTau - taus < epsilon){
     epsilon = maxTau-taus;
   }
-  //std::cout << GridLogMessage << "Integration epsilon : " << epsilon << std::endl;
+  std::cout << GridLogMessage << "Integration epsilon : " << epsilon << std::endl;
   GaugeField Z(U.Grid());
   GaugeField Zprime(U.Grid());
-  GaugeField tmp(U.Grid()), Uprime(U.Grid());
+  GaugeField tmp(U.Grid()), Uprime(U.Grid()),Usave(U.Grid());
+
   Uprime = U;
+  Usave = U;
+
   SG.deriv(U, Z);
   Zprime = -Z;
   Z *= 0.25;                                  // Z0 = 1/4 * F(U)
@@ -128,18 +136,33 @@ void WilsonFlow<Gimpl>::evolve_step_adaptive(typename Gimpl::GaugeField &U, Real
   Z *= 3.0/4.0;                               // Z = 17/36*Z0 -8/9*Z1 +3/4*Z2
   Gimpl::update_field(Z, U, -2.0*epsilon);    // V(t+e) = exp(ep*Z)*W2
 
-  // Ramos 
+  // Ramos arXiv:1301.4388
   Gimpl::update_field(Zprime, Uprime, -2.0*epsilon); // V'(t+e) = exp(ep*Z')*W0
   // Compute distance as norm^2 of the difference
   GaugeField diffU = U - Uprime;
-  RealD diff = norm2(diffU);
+// Wrong
-  // adjust integration step
+//  RealD diff = norm2(diffU);
+//  std::cout << GridLogMessage << "norm2: " << diff << std::endl;
   
+//  RealD tol=1e-3;
+  
+  RealD diff = real(rankInnerMax(diffU,diffU));
+  diff = sqrt(diff)/18.; // distance defined in Ramos 
+
+  GridBase *grid = diffU.Grid();
+  std::cout << GridLogMessage << "max: " << diff << std::endl;
+  grid->GlobalMax(diff);
+  std::cout << GridLogMessage << "max: " << diff << std::endl;
+    
+  if(diff < tolerance) {
   taus += epsilon;
-  //std::cout << GridLogMessage << "Adjusting integration step with distance: " << diff << std::endl;
+//  std::cout << GridLogMessage << "Adjusting integration step with distance: " << diff << std::endl;
+  } else {
+    U = Usave;
+  }
     
-  epsilon = epsilon*0.95*std::pow(1e-4/diff,1./3.);
+  epsilon = epsilon*0.95*std::pow(tolerance/diff,1./3.);
-  //std::cout << GridLogMessage << "New epsilon : " << epsilon << std::endl;
+  std::cout << GridLogMessage << "Distance : "<<diff<<"New epsilon : " << epsilon << std::endl;
 
 }
 
@@ -184,8 +207,11 @@ void WilsonFlow<Gimpl>::smear(GaugeField& out, const GaugeField& in) const {
 template <class Gimpl>
 void WilsonFlow<Gimpl>::smear_adaptive(GaugeField& out, const GaugeField& in, RealD maxTau){
   out = in;
-  taus = epsilon;
+//  taus = epsilon;
+  taus = 0.;
   unsigned int step = 0;
+  double measTau = epsilon*measure_interval;
+  std::cout << GridLogMessage << "measTau :"<< measTau << std::endl;
   do{
     step++;
     //std::cout << GridLogMessage << "Evolution time :"<< taus << std::endl;
@@ -193,10 +219,12 @@ void WilsonFlow<Gimpl>::smear_adaptive(GaugeField& out, const GaugeField& in, Re
     std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : "
 		  << step << "  " << taus << "  "
 	      << energyDensityPlaquette(out) << std::endl;
-    if( step % measure_interval == 0){
+//    if( step % measure_interval == 0){
+    if( taus >  measTau ) {
       std::cout << GridLogMessage << "[WilsonFlow] Top. charge           : "
 		<< step << "  " 
 		<< WilsonLoops<PeriodicGimplR>::TopologicalCharge(out) << std::endl;
+      measTau += epsilon*measure_interval;
     }
   } while (taus < maxTau);
 

Grid/serialisation/BaseIO.cc

@@ -1,35 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/serialisation/BaseIO.h
-
-Copyright (C) 2015
-
-Author: Michael Marshall <michael.marshall@ed.ac.uk>
-
-This program 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.
-
-This program 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 this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-/*  END LEGAL */
-
-#include <Grid/GridCore.h>
-
-NAMESPACE_BEGIN(Grid)
-
-std::uint64_t EigenIO::EigenResizeCounter(0);
-
-NAMESPACE_END(Grid)

Grid/serialisation/BaseIO.h

@@ -9,7 +9,6 @@
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 Author: Guido Cossu <guido.cossu@ed.ac.uk>
-Author: Michael Marshall <michael.marshall@ed.ac.uk>
 
     This program is free software; you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
@@ -31,7 +30,6 @@ Author: Michael Marshall <michael.marshall@ed.ac.uk>
 #ifndef GRID_SERIALISATION_ABSTRACT_READER_H
 #define GRID_SERIALISATION_ABSTRACT_READER_H
 
-#include <atomic>
 #include <type_traits>
 #include <Grid/tensors/Tensors.h>
 #include <Grid/serialisation/VectorUtils.h>
@@ -112,10 +110,6 @@ namespace Grid {
     template <typename ET>
     inline typename std::enable_if<is_tensor_of_container<ET>::value, typename Traits<ET>::scalar_type *>::type
     getFirstScalar(ET &eigenTensor) { return eigenTensor.data()->begin(); }
-
-    // Counter for resized EigenTensors (poor man's substitute for allocator)
-    // Defined in BinaryIO.cc
-    extern std::uint64_t EigenResizeCounter;
   }
 
   // Abstract writer/reader classes ////////////////////////////////////////////
@@ -503,14 +497,8 @@ namespace Grid {
   typename std::enable_if<EigenIO::is_tensor_variable<ETensor>::value, void>::type
   Reader<T>::Reshape(ETensor &t, const std::array<typename ETensor::Index, ETensor::NumDimensions> &dims )
   {
-#ifdef GRID_OMP
-    // The memory counter is the reason this must be done from the primary thread
-    assert(omp_in_parallel()==0 && "Deserialisation which resizes Eigen tensor must happen from primary thread");
-#endif
-    EigenIO::EigenResizeCounter -= static_cast<uint64_t>(t.size()) * sizeof(typename ETensor::Scalar);
     //t.reshape( dims );
     t.resize( dims );
-    EigenIO::EigenResizeCounter += static_cast<uint64_t>(t.size()) * sizeof(typename ETensor::Scalar);
   }
 
   template <typename T>

Grid/serialisation/Hdf5IO.cc

@@ -1,34 +1,3 @@
-/*************************************************************************************
- 
- Grid physics library, www.github.com/paboyle/Grid
- 
- Source file: ./Grid/serialisation/VectorUtils.h
- 
- Copyright (C) 2015
- 
- Author: Antonin Portelli <antonin.portelli@me.com>
- Author: Peter Boyle <paboyle@ed.ac.uk>
- Author: Guido Cossu <guido.cossu@ed.ac.uk>
- Author: Michael Marshall <michael.marshall@ed.ac.uk>
-
- This program 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.
- 
- This program 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 this program; if not, write to the Free Software Foundation, Inc.,
- 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- 
- See the full license in the file "LICENSE" in the top level distribution directory
- *************************************************************************************/
-/*  END LEGAL */
-
 #include <Grid/Grid.h>
 
 using namespace Grid;

Grid/serialisation/Hdf5IO.h

@@ -1,34 +1,3 @@
-/*************************************************************************************
- 
- Grid physics library, www.github.com/paboyle/Grid
- 
- Source file: ./Grid/serialisation/VectorUtils.h
- 
- Copyright (C) 2015
- 
- Author: Peter Boyle <paboyle@ed.ac.uk>
- Author: Antonin Portelli <antonin.portelli@me.com>
- Author: Guido Cossu <guido.cossu@ed.ac.uk>
- Author: Michael Marshall <michael.marshall@ed.ac.uk>
-
- This program 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.
- 
- This program 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 this program; if not, write to the Free Software Foundation, Inc.,
- 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- 
- See the full license in the file "LICENSE" in the top level distribution directory
- *************************************************************************************/
-/*  END LEGAL */
-
 #ifndef GRID_SERIALISATION_HDF5_H
 #define GRID_SERIALISATION_HDF5_H
 
@@ -65,13 +34,11 @@ namespace Grid
     template <typename U>
     void writeDefault(const std::string &s, const U &x);
     template <typename U>
-    void writeRagged(const std::string &s, const std::vector<U> &x);
+    typename std::enable_if<element<std::vector<U>>::is_number, void>::type
-    template <typename U>
-    typename std::enable_if<is_flattenable<std::vector<U>>::value>::type
     writeDefault(const std::string &s, const std::vector<U> &x);
     template <typename U>
-    typename std::enable_if<!is_flattenable<std::vector<U>>::value>::type
+    typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
-    writeDefault(const std::string &s, const std::vector<U> &x) { writeRagged(s, x); }
+    writeDefault(const std::string &s, const std::vector<U> &x);
     template <typename U>
     void writeMultiDim(const std::string &s, const std::vector<size_t> & Dimensions, const U * pDataRowMajor, size_t NumElements);
     H5NS::Group & getGroup(void);
@@ -97,13 +64,11 @@ namespace Grid
     template <typename U>
     void readDefault(const std::string &s, U &output);
     template <typename U>
-    void readRagged(const std::string &s, std::vector<U> &x);
+    typename std::enable_if<element<std::vector<U>>::is_number, void>::type
-    template <typename U>
-    typename std::enable_if<is_flattenable<std::vector<U>>::value>::type
     readDefault(const std::string &s, std::vector<U> &x);
     template <typename U>
-    typename std::enable_if<!is_flattenable<std::vector<U>>::value>::type
+    typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
-    readDefault(const std::string &s, std::vector<U> &x) { readRagged(s, x); }
+    readDefault(const std::string &s, std::vector<U> &x);
     template <typename U>
     void readMultiDim(const std::string &s, std::vector<U> &buf, std::vector<size_t> &dim);
     H5NS::Group & getGroup(void);
@@ -211,13 +176,11 @@ namespace Grid
   }
 
   template <typename U>
-  typename std::enable_if<is_flattenable<std::vector<U>>::value>::type
+  typename std::enable_if<element<std::vector<U>>::is_number, void>::type
   Hdf5Writer::writeDefault(const std::string &s, const std::vector<U> &x)
-  {
-    if (isRegularShape(x))
   {
     // alias to element type
-      using Scalar = typename is_flattenable<std::vector<U>>::type;
+    typedef typename element<std::vector<U>>::type Element;
     
     // flatten the vector and getting dimensions
     Flatten<std::vector<U>> flat(x);
@@ -225,16 +188,12 @@ namespace Grid
     const auto           &flatx = flat.getFlatVector();
     for (auto &d: flat.getDim())
       dim.push_back(d);
-      writeMultiDim<Scalar>(s, dim, &flatx[0], flatx.size());
+    writeMultiDim<Element>(s, dim, &flatx[0], flatx.size());
-    }
-    else
-    {
-      writeRagged(s, x);
-    }
   }
   
   template <typename U>
-  void Hdf5Writer::writeRagged(const std::string &s, const std::vector<U> &x)
+  typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
+  Hdf5Writer::writeDefault(const std::string &s, const std::vector<U> &x)
   {
     push(s);
     writeSingleAttribute(x.size(), HDF5_GRID_GUARD "vector_size",
@@ -270,7 +229,7 @@ namespace Grid
   void Hdf5Reader::readMultiDim(const std::string &s, std::vector<U> &buf, std::vector<size_t> &dim)
   {
     // alias to element type
-    using Scalar = typename is_flattenable<std::vector<U>>::type;
+    typedef typename element<std::vector<U>>::type Element;
     
     // read the dimensions
     H5NS::DataSpace       dataSpace;
@@ -301,33 +260,26 @@ namespace Grid
       H5NS::DataSet dataSet;
       
       dataSet = group_.openDataSet(s);
-      dataSet.read(buf.data(), Hdf5Type<Scalar>::type());
+      dataSet.read(buf.data(), Hdf5Type<Element>::type());
     }
     else
     {
       H5NS::Attribute attribute;
       
       attribute = group_.openAttribute(s);
-      attribute.read(Hdf5Type<Scalar>::type(), buf.data());
+      attribute.read(Hdf5Type<Element>::type(), buf.data());
     }
   }
 
   template <typename U>
-  typename std::enable_if<is_flattenable<std::vector<U>>::value>::type
+  typename std::enable_if<element<std::vector<U>>::is_number, void>::type
   Hdf5Reader::readDefault(const std::string &s, std::vector<U> &x)
-  {
-    if (H5Lexists        (group_.getId(), s.c_str(), H5P_DEFAULT) > 0
-     && H5Aexists_by_name(group_.getId(), s.c_str(), HDF5_GRID_GUARD "vector_size", H5P_DEFAULT ) > 0)
-    {
-      readRagged(s, x);
-    }
-    else
   {
     // alias to element type
-      using Scalar = typename is_flattenable<std::vector<U>>::type;
+    typedef typename element<std::vector<U>>::type Element;
 
     std::vector<size_t>   dim;
-      std::vector<Scalar>   buf;
+    std::vector<Element>  buf;
     readMultiDim( s, buf, dim );
 
     // reconstruct the multidimensional vector
@@ -335,10 +287,10 @@ namespace Grid
     
     x = r.getVector();
   }
-  }
   
   template <typename U>
-  void Hdf5Reader::readRagged(const std::string &s, std::vector<U> &x)
+  typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
+  Hdf5Reader::readDefault(const std::string &s, std::vector<U> &x)
   {
     uint64_t size;
     

Grid/serialisation/MacroMagic.h

@@ -118,13 +118,13 @@ static inline std::string SerialisableClassName(void) {return std::string(#cname
 static constexpr bool isEnum = false; \
 GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_MEMBER,__VA_ARGS__))\
 template <typename T>\
-static inline void write(::Grid::Writer<T> &WR,const std::string &s, const cname &obj){ \
+static inline void write(Writer<T> &WR,const std::string &s, const cname &obj){ \
   push(WR,s);\
   GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_WRITE_MEMBER,__VA_ARGS__))	\
   pop(WR);\
 }\
 template <typename T>\
-static inline void read(::Grid::Reader<T> &RD,const std::string &s, cname &obj){	\
+static inline void read(Reader<T> &RD,const std::string &s, cname &obj){	\
   if (!push(RD,s))\
   {\
     std::cout << ::Grid::GridLogWarning << "IO: Cannot open node '" << s << "'" << std::endl; \

Grid/serialisation/VectorUtils.h

@@ -9,7 +9,6 @@
  Author: Antonin Portelli <antonin.portelli@me.com>
  Author: Peter Boyle <paboyle@ph.ed.ac.uk>
  Author: paboyle <paboyle@ph.ed.ac.uk>
- Author: Michael Marshall <michael.marshall@ed.ac.uk>
  
  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
@@ -237,36 +236,21 @@ namespace Grid {
     }
   }
 
-  // is_flattenable<T>::value is true if T is a std::vector<> which can be flattened //////////////////////
+  // Vector element trait //////////////////////////////////////////////////////  
-  template <typename T, typename V = void>
+  template <typename T>
-  struct is_flattenable : std::false_type
+  struct element
   {
-    using type      = T;
+    typedef T type;
-    using grid_type = T;
+    static constexpr bool is_number = false;
-    static constexpr int vecRank = 0;
-    static constexpr bool isGridTensor = false;
-    static constexpr bool children_flattenable = std::is_arithmetic<T>::value or is_complex<T>::value;
   };
   
   template <typename T>
-  struct is_flattenable<T, typename std::enable_if<isGridTensor<T>::value>::type> : std::false_type
+  struct element<std::vector<T>>
   {
-    using type      = typename GridTypeMapper<T>::scalar_type;
+    typedef typename element<T>::type type;
-    using grid_type = T;
+    static constexpr bool is_number = std::is_arithmetic<T>::value
-    static constexpr int vecRank = 0;
+                                      or is_complex<T>::value
-    static constexpr bool isGridTensor = true;
+                                      or element<T>::is_number;
-    static constexpr bool children_flattenable = true;
-  };
-
-  template <typename T>
-  struct is_flattenable<std::vector<T>, typename std::enable_if<is_flattenable<T>::children_flattenable>::type>
-  : std::true_type
-  {
-    using type      = typename is_flattenable<T>::type;
-    using grid_type = typename is_flattenable<T>::grid_type;
-    static constexpr bool isGridTensor = is_flattenable<T>::isGridTensor;
-    static constexpr int vecRank = is_flattenable<T>::vecRank + 1;
-    static constexpr bool children_flattenable = true;
   };
   
   // Vector flattening utility class ////////////////////////////////////////////
@@ -275,29 +259,22 @@ namespace Grid {
   class Flatten
   {
   public:
-    using Scalar  = typename is_flattenable<V>::type;
+    typedef typename element<V>::type Element;
-    static constexpr bool isGridTensor = is_flattenable<V>::isGridTensor;
   public:
     explicit                     Flatten(const V &vector);
-    const V &                   getVector(void)     const { return vector_; }
+    const V &                    getVector(void);
-    const std::vector<Scalar> & getFlatVector(void) const { return flatVector_; }
+    const std::vector<Element> & getFlatVector(void);
-    const std::vector<size_t> & getDim(void)        const { return dim_; }
+    const std::vector<size_t>  & getDim(void);
   private:
-    template <typename W> typename std::enable_if<!is_flattenable<W>::value && !is_flattenable<W>::isGridTensor>::type
+    void accumulate(const Element &e);
-    accumulate(const W &e);
+    template <typename W>
-    template <typename W> typename std::enable_if<!is_flattenable<W>::value &&  is_flattenable<W>::isGridTensor>::type
+    void accumulate(const W &v);
-    accumulate(const W &e);
+    void accumulateDim(const Element &e);
-    template <typename W> typename std::enable_if< is_flattenable<W>::value>::type
+    template <typename W>
-    accumulate(const W &v);
+    void accumulateDim(const W &v);
-    template <typename W> typename std::enable_if<!is_flattenable<W>::value && !is_flattenable<W>::isGridTensor>::type
-    accumulateDim(const W &e) {} // Innermost is a scalar - do nothing
-    template <typename W> typename std::enable_if<!is_flattenable<W>::value &&  is_flattenable<W>::isGridTensor>::type
-    accumulateDim(const W &e);
-    template <typename W> typename std::enable_if< is_flattenable<W>::value>::type
-    accumulateDim(const W &v);
   private:
     const V              &vector_;
-    std::vector<Scalar> flatVector_;
+    std::vector<Element> flatVector_;
     std::vector<size_t>  dim_;
   };
   
@@ -306,32 +283,23 @@ namespace Grid {
   class Reconstruct
   {
   public:
-    using Scalar  = typename is_flattenable<V>::type;
+    typedef typename element<V>::type Element;
-    static constexpr bool isGridTensor = is_flattenable<V>::isGridTensor;
   public:
-    Reconstruct(const std::vector<Scalar> &flatVector,
+    Reconstruct(const std::vector<Element> &flatVector,
                 const std::vector<size_t> &dim);
-    const V &                   getVector(void)     const { return vector_; }
+    const V &                    getVector(void);
-    const std::vector<Scalar> & getFlatVector(void) const { return flatVector_; }
+    const std::vector<Element> & getFlatVector(void);
-    const std::vector<size_t> & getDim(void)        const { return dim_; }
+    const std::vector<size_t>  & getDim(void);
   private:
-    template <typename W> typename std::enable_if<!is_flattenable<W>::value && !is_flattenable<W>::isGridTensor>::type
+    void fill(std::vector<Element> &v);
-    fill(W &v);
+    template <typename W>
-    template <typename W> typename std::enable_if<!is_flattenable<W>::value &&  is_flattenable<W>::isGridTensor>::type
+    void fill(W &v);
-    fill(W &v);
+    void resize(std::vector<Element> &v, const unsigned int dim);
-    template <typename W> typename std::enable_if< is_flattenable<W>::value>::type
+    template <typename W>
-    fill(W &v);
+    void resize(W &v, const unsigned int dim);
-    template <typename W> typename std::enable_if< is_flattenable<W>::value &&  is_flattenable<W>::vecRank==1>::type
-    resize(W &v, const unsigned int dim);
-    template <typename W> typename std::enable_if< is_flattenable<W>::value && (is_flattenable<W>::vecRank>1)>::type
-    resize(W &v, const unsigned int dim);
-    template <typename W> typename std::enable_if<!is_flattenable<W>::isGridTensor>::type
-    checkInnermost(const W &e) {} // Innermost is a scalar - do nothing
-    template <typename W> typename std::enable_if< is_flattenable<W>::isGridTensor>::type
-    checkInnermost(const W &e);
   private:
     V                          vector_;
-    const std::vector<Scalar> &flatVector_;
+    const std::vector<Element> &flatVector_;
     std::vector<size_t>        dim_;
     size_t                     ind_{0};
     unsigned int               dimInd_{0};
@@ -339,24 +307,14 @@ namespace Grid {
 
   // Flatten class template implementation
   template <typename V>
-  template <typename W> typename std::enable_if<!is_flattenable<W>::value && !is_flattenable<W>::isGridTensor>::type
+  void Flatten<V>::accumulate(const Element &e)
-  Flatten<V>::accumulate(const W &e)
   {
     flatVector_.push_back(e);
   }
   
   template <typename V>
-  template <typename W> typename std::enable_if<!is_flattenable<W>::value && is_flattenable<W>::isGridTensor>::type
+  template <typename W>
-  Flatten<V>::accumulate(const W &e)
+  void Flatten<V>::accumulate(const W &v)
-  {
-    for (const Scalar &x: e) {
-      flatVector_.push_back(x);
-    }
-  }
-
-  template <typename V>
-  template <typename W> typename std::enable_if<is_flattenable<W>::value>::type
-  Flatten<V>::accumulate(const W &v)
   {
     for (auto &e: v)
     {
@@ -365,17 +323,11 @@ namespace Grid {
   }
   
   template <typename V>
-  template <typename W> typename std::enable_if<!is_flattenable<W>::value && is_flattenable<W>::isGridTensor>::type
+  void Flatten<V>::accumulateDim(const Element &e) {};
-  Flatten<V>::accumulateDim(const W &e)
-  {
-    using Traits = GridTypeMapper<typename is_flattenable<W>::grid_type>;
-    for (int rank=0; rank < Traits::Rank; ++rank)
-      dim_.push_back(Traits::Dimension(rank));
-  }
   
   template <typename V>
-  template <typename W> typename std::enable_if<is_flattenable<W>::value>::type
+  template <typename W>
-  Flatten<V>::accumulateDim(const W &v)
+  void Flatten<V>::accumulateDim(const W &v)
   {
     dim_.push_back(v.size());
     accumulateDim(v[0]);
@@ -385,26 +337,32 @@ namespace Grid {
   Flatten<V>::Flatten(const V &vector)
   : vector_(vector)
   {
-    accumulateDim(vector_);
-    std::size_t TotalSize{ dim_[0] };
-    for (int i = 1; i < dim_.size(); ++i) {
-      TotalSize *= dim_[i];
-    }
-    flatVector_.reserve(TotalSize);
     accumulate(vector_);
+    accumulateDim(vector_);
+  }
+  
+  template <typename V>
+  const V & Flatten<V>::getVector(void)
+  {
+    return vector_;
+  }
+  
+  template <typename V>
+  const std::vector<typename Flatten<V>::Element> &
+  Flatten<V>::getFlatVector(void)
+  {
+    return flatVector_;
+  }
+  
+  template <typename V>
+  const std::vector<size_t> & Flatten<V>::getDim(void)
+  {
+    return dim_;
   }
   
   // Reconstruct class template implementation
   template <typename V>
-  template <typename W> typename std::enable_if<!is_flattenable<W>::value && !is_flattenable<W>::isGridTensor>::type
+  void Reconstruct<V>::fill(std::vector<Element> &v)
-  Reconstruct<V>::fill(W &v)
-  {
-    v = flatVector_[ind_++];
-  }
-  
-  template <typename V>
-  template <typename W> typename std::enable_if<!is_flattenable<W>::value &&  is_flattenable<W>::isGridTensor>::type
-  Reconstruct<V>::fill(W &v)
   {
     for (auto &e: v)
     {
@@ -413,8 +371,8 @@ namespace Grid {
   }
   
   template <typename V>
-  template <typename W> typename std::enable_if<is_flattenable<W>::value>::type
+  template <typename W>
-  Reconstruct<V>::fill(W &v)
+  void Reconstruct<V>::fill(W &v)
   {
     for (auto &e: v)
     {
@@ -423,15 +381,14 @@ namespace Grid {
   }
   
   template <typename V>
-  template <typename W> typename std::enable_if<is_flattenable<W>::value && is_flattenable<W>::vecRank==1>::type
+  void Reconstruct<V>::resize(std::vector<Element> &v, const unsigned int dim)
-  Reconstruct<V>::resize(W &v, const unsigned int dim)
   {
     v.resize(dim_[dim]);
   }
   
   template <typename V>
-  template <typename W> typename std::enable_if<is_flattenable<W>::value && (is_flattenable<W>::vecRank>1)>::type
+  template <typename W>
-  Reconstruct<V>::resize(W &v, const unsigned int dim)
+  void Reconstruct<V>::resize(W &v, const unsigned int dim)
   {
     v.resize(dim_[dim]);
     for (auto &e: v)
@@ -441,31 +398,34 @@ namespace Grid {
   }
   
   template <typename V>
-  template <typename W> typename std::enable_if<is_flattenable<W>::isGridTensor>::type
+  Reconstruct<V>::Reconstruct(const std::vector<Element> &flatVector,
-  Reconstruct<V>::checkInnermost(const W &)
-  {
-    using Traits = GridTypeMapper<typename is_flattenable<W>::grid_type>;
-    const int gridRank{Traits::Rank};
-    const int dimRank{static_cast<int>(dim_.size())};
-    assert(dimRank >= gridRank && "Tensor rank too low for Grid tensor");
-    for (int i=0; i<gridRank; ++i) {
-      assert(dim_[dimRank - gridRank + i] == Traits::Dimension(i) && "Tensor dimension doesn't match Grid tensor");
-    }
-    dim_.resize(dimRank - gridRank);
-  }
-
-  template <typename V>
-  Reconstruct<V>::Reconstruct(const std::vector<Scalar> &flatVector,
                               const std::vector<size_t> &dim)
   : flatVector_(flatVector)
   , dim_(dim)
   {
-    checkInnermost(vector_);
-    assert(dim_.size() == is_flattenable<V>::vecRank && "Tensor rank doesn't match nested std::vector rank");
     resize(vector_, 0);
     fill(vector_);
   }
   
+  template <typename V>
+  const V & Reconstruct<V>::getVector(void)
+  {
+    return vector_;
+  }
+  
+  template <typename V>
+  const std::vector<typename Reconstruct<V>::Element> &
+  Reconstruct<V>::getFlatVector(void)
+  {
+    return flatVector_;
+  }
+  
+  template <typename V>
+  const std::vector<size_t> & Reconstruct<V>::getDim(void)
+  {
+    return dim_;
+  }
+
   // Vector IO utilities ///////////////////////////////////////////////////////
   // helper function to read space-separated values
   template <typename T>
@@ -499,64 +459,6 @@ namespace Grid {
     
     return os;
   }
-
-  // In general, scalar types are considered "flattenable" (regularly shaped)
-  template <typename T>
-  bool isRegularShapeHelper(const std::vector<T> &, std::vector<std::size_t> &, int, bool)
-  {
-    return true;
-  }
-
-  template <typename T>
-  bool isRegularShapeHelper(const std::vector<std::vector<T>> &v, std::vector<std::size_t> &Dims, int Depth, bool bFirst)
-  {
-    if( bFirst)
-    {
-      assert( Dims.size() == Depth     && "Bug: Delete this message after testing" );
-      Dims.push_back(v[0].size());
-      if (!Dims[Depth])
-        return false;
-    }
-    else
-    {
-      assert( Dims.size() >= Depth + 1 && "Bug: Delete this message after testing" );
-    }
-    for (std::size_t i = 0; i < v.size(); ++i)
-    {
-      if (v[i].size() != Dims[Depth] || !isRegularShapeHelper(v[i], Dims, Depth + 1, bFirst && i==0))
-      {
-        return false;
-      }
-    }
-    return true;
-  }
-
-  template <typename T>
-  bool isRegularShape(const T &t) { return true; }
-
-  template <typename T>
-  bool isRegularShape(const std::vector<T> &v) { return !v.empty(); }
-
-  // Return non-zero if all dimensions of this std::vector<std::vector<T>> are regularly shaped
-  template <typename T>
-  bool isRegularShape(const std::vector<std::vector<T>> &v)
-  {
-    if (v.empty() || v[0].empty())
-      return false;
-    // Make sure all of my rows are the same size
-    std::vector<std::size_t> Dims;
-    Dims.reserve(is_flattenable<T>::vecRank);
-    Dims.push_back(v.size());
-    Dims.push_back(v[0].size());
-    for (std::size_t i = 0; i < Dims[0]; ++i)
-    {
-      if (v[i].size() != Dims[1] || !isRegularShapeHelper(v[i], Dims, 2, i==0))
-      {
-        return false;
-      }
-    }
-    return true;
-  }
 }
 
 // helper function to read space-separated values

Grid/tensors/Tensor_class.h

@@ -417,7 +417,7 @@ public:
       stream << "{";
       for (int j = 0; j < N; j++) {
 	stream << o._internal[i][j];
-	if (j < N - 1) stream << ",";
+	if (i < N - 1) stream << ",";
       }
       stream << "}";
       if (i != N - 1) stream << "\n\t\t";

documentation/manual.rst

@@ -1787,7 +1787,7 @@ Hdf5Writer     Hdf5Reader       HDF5
 
 Write interfaces, similar to the XML facilities in QDP++ are presented. However,
 the serialisation routines are automatically generated by the macro, and a virtual
-reader and writer interface enables writing to any of a number of formats.
+reader adn writer interface enables writing to any of a number of formats.
 
 **Example**::
 
@@ -1814,91 +1814,6 @@ reader and writer interface enables writing to any of a number of formats.
    }
 
 
-Eigen tensor support -- added 2019H1
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-
-The Serialisation library was expanded in 2019 to support de/serialisation of
-Eigen tensors. De/serialisation of existing types was not changed. Data files
-without Eigen tensors remain compatible with earlier versions of Grid and other readers.
-Conversely, data files containing serialised Eigen tensors is a breaking change.
-
-Eigen tensor serialisation support was added to BaseIO, which was modified to provide a Traits class
-to recognise Eigen tensors with elements that are either: primitive scalars (arithmetic and complex types);
-or Grid tensors.
-
-**Traits determining de/serialisable scalars**::
-
-  // Is this an Eigen tensor
-  template<typename T> struct is_tensor : std::integral_constant<bool,
-    std::is_base_of<Eigen::TensorBase<T, Eigen::ReadOnlyAccessors>, T>::value> {};
-  // Is this an Eigen tensor of a supported scalar
-  template<typename T, typename V = void> struct is_tensor_of_scalar : public std::false_type {};
-  template<typename T> struct is_tensor_of_scalar<T, typename std::enable_if<is_tensor<T>::value && is_scalar<typename T::Scalar>::value>::type> : public std::true_type {};
-  // Is this an Eigen tensor of a supported container
-  template<typename T, typename V = void> struct is_tensor_of_container : public std::false_type {};
-  template<typename T> struct is_tensor_of_container<T, typename std::enable_if<is_tensor<T>::value && isGridTensor<typename T::Scalar>::value>::type> : public std::true_type {};
-
-
-Eigen tensors are regular, multidimensional objects, and each Reader/Writer
-was extended to support this new datatype. Where the Eigen tensor contains
-a Grid tensor, the dimensions of the data written are the dimensions of the
-Eigen tensor plus the dimensions of the underlying Grid scalar. Dimensions
-of size 1 are preserved.
-
-**New Reader/Writer methods for multi-dimensional data**::
-
-  template <typename U>
-  void readMultiDim(const std::string &s, std::vector<U> &buf, std::vector<size_t> &dim);
-  template <typename U>
-  void writeMultiDim(const std::string &s, const std::vector<size_t> & Dimensions, const U * pDataRowMajor, size_t NumElements);
-
-
-On readback, the Eigen tensor rank must match the data being read, but the tensor
-dimensions will be resized if necessary. Resizing is not possible for Eigen::TensorMap<T>
-because these tensors use a buffer provided at construction, and this buffer cannot be changed.
-Deserialisation failures cause Grid to assert.
-
-
-HDF5 Optimisations -- added June 2021
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-
-Grid serialisation is intended to be light, deterministic and provide a layer of abstraction over
-multiple file formats. HDF5 excels at handling multi-dimensional data, and the Grid HDF5Reader/HDF5Writer exploits this.
-When serialising nested ``std::vector<T>``, where ``T`` is an arithmetic or complex type,
-the Hdf5Writer writes the data as an Hdf5 DataSet object.
-
-However, nested ``std::vector<std::vector<...T>>`` might be "ragged", i.e. not necessarily regular. E.g. a 3d nested
-``std::vector`` might contain 2 rows, the first being a 2x2 block and the second row being a 1 x 2 block.
-A bug existed whereby this was not checked on write, so nested, ragged vectors
-were written as a regular dataset, with a buffer under/overrun and jumbled contents.
-
-Clearly this was not used in production, as the bug went undetected until now. Fixing this bug
-is an opportunity to further optimise the HDF5 file format.
-
-The goals of this change are to:
-
-* Make changes to the Hdf5 file format only -- i.e. do not impact other file formats
-
-* Implement file format changes in such a way that they are transparent to the Grid reader
-
-* Correct the bug for ragged vectors of numeric / complex types
-
-* Extend the support of nested std::vector<T> to arbitrarily nested Grid tensors
-
-
-The trait class ``element`` has been redefined to ``is_flattenable``, which is a trait class for
-potentially "flattenable" objects. These are (possibly nested) ``std::vector<T>`` where ``T`` is
-an arithmetic, complex or Grid tensor type. Flattenable objects are tested on write
-(with the function ``isRegularShape``) to see whether they actually are regular.
-
-Flattenable, regular objects are written to a multidimensional HDF5 DataSet.
-Otherwise, an Hdf5 sub group is created with the object "name", and each element of the outer dimension is
-recursively written to as object "name_n", where n is a 0-indexed number.
-
-On readback (by Grid)), the presence of a subgroup containing the attribute ``Grid_vector_size`` triggers a
-"ragged read", otherwise a read from a DataSet is attempted.
-
-
 Data parallel field IO
 -----------------------
 

tests/IO/Test_serialisation.cc

@@ -48,9 +48,7 @@ public:
                           std::vector<double>, array,
                           std::vector<std::vector<double> >, twodimarray,
 			  std::vector<std::vector<std::vector<std::complex<double>> > >, cmplx3darray,
-			  SpinColourMatrix, scm,
+			  SpinColourMatrix, scm
-                          std::vector<std::vector<std::vector<int> > >, ragged,
-                          std::vector<std::vector<SpinColourMatrix> >, vscm
                           );
   myclass() {}
   myclass(int i)
@@ -58,10 +56,6 @@ public:
   , twodimarray(3,std::vector<double>(5, 1.23456))
   , cmplx3darray(3,std::vector<std::vector<std::complex<double>>>(5, std::vector<std::complex<double>>(7, std::complex<double>(1.2, 3.4))))
   , ve(2, myenum::blue)
-  , ragged( {{{i+1},{i+2,i+3}}, // ragged
-            {{i+4,i+5,i+6,i+7},{i+8,i+9,i+10,i+11},{i+12,i+13,i+14,i+15}}, // block
-            {{i+16,i+17},{i+18,i+19,i+20}}} ) //ragged
-  , vscm(3, std::vector<SpinColourMatrix>(5))
   {
     e=myenum::red;
     x=i;
@@ -74,13 +68,6 @@ public:
     scm()(0, 2)(1, 1) = 6.336;
     scm()(2, 1)(2, 2) = 7.344;
     scm()(1, 1)(2, 0) = 8.3534;
-    int Counter = i;
-    for( auto & v : vscm ) {
-      for( auto & j : v ) {
-        j = std::complex<double>(Counter, -Counter);
-        Counter++;
-      }
-    }
   }
 };
 

tests/smearing/Test_WilsonFlow.cc

@@ -33,6 +33,7 @@ namespace Grid{
     GRID_SERIALIZABLE_CLASS_MEMBERS(WFParameters,
             int, steps,
             double, step_size,
+            double, tol,
             int, meas_interval,
             double, maxTau); // for the adaptive algorithm
        
@@ -82,13 +83,27 @@ int main(int argc, char **argv) {
   SU<Nc>::HotConfiguration(pRNG, Umu);
   
   typedef Grid::XmlReader       Serialiser;
-  Serialiser Reader("input.xml");
+//  Serialiser Reader("input.xml");
-  WFParameters WFPar(Reader);
+//  WFParameters WFPar(Reader);
-  ConfParameters CPar(Reader);
+//  ConfParameters CPar(Reader);
-  CheckpointerParameters CPPar(CPar.conf_prefix, CPar.rng_prefix);
+//  WFParameters WFPar;
+  int steps = 800;
+  double step_size=0.02;
+  double tol=1e-4;
+  int meas_interval=50;
+  double maxTau = 16;
+//  ConfParameters CPar;
+//  CPar. conf_prefix="configurations/ckpoint_lat";
+//  CPar. rng_prefix="rngs/ckpoint_rng";
+//  CPar. StartConfiguration=100,
+//  CPar. EndConfiguration=110,
+//  CPar. Skip=1;
+//  CheckpointerParameters CPPar(CPar.conf_prefix, CPar.rng_prefix);
+  CheckpointerParameters CPPar("configurations/ckpoint_lat","rngs/ckpoint_rng");
   BinaryHmcCheckpointer<PeriodicGimplR> CPBin(CPPar);
 
-  for (int conf = CPar.StartConfiguration; conf <= CPar.EndConfiguration; conf+= CPar.Skip){
+//  for (int conf = CPar.StartConfiguration; conf <= CPar.EndConfiguration; conf+= CPar.Skip){
+  for (int conf = 100; conf <= 110; conf+= 1){
 
   CPBin.CheckpointRestore(conf, Umu, sRNG, pRNG);
 
@@ -96,9 +111,10 @@ int main(int argc, char **argv) {
   std::cout << GridLogMessage << "Initial plaquette: "
     << WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu) << std::endl;
 
-  WilsonFlow<PeriodicGimplR> WF(WFPar.steps, WFPar.step_size, WFPar.meas_interval);
+  WilsonFlow<PeriodicGimplR> WF(steps, step_size, meas_interval);
 
-  WF.smear_adaptive(Uflow, Umu, WFPar.maxTau);
+//  WF.smear_adaptive(Uflow, Umu, maxTau);
+  WF.smear(Uflow, Umu);
 
   RealD WFlow_plaq = WilsonLoops<PeriodicGimplR>::avgPlaquette(Uflow);
   RealD WFlow_TC   = WilsonLoops<PeriodicGimplR>::TopologicalCharge(Uflow);