Switched perambulator to sue Eigen::Tensor (file write temporarily excluded)

Hadrons/Modules/MDistil/Distil.hpp

@@ -179,7 +179,7 @@ inline GridCartesian * MakeLowerDimGrid( GridCartesian * gridHD )
 {
   //LOG(Message) << "MakeLowerDimGrid() begin" << std::endl;
   int nd{static_cast<int>(gridHD->_ndimension)};
-  std::vector<int> latt_size   = gridHD->_fdimensions;
+  std::vector<int> latt_size   = gridHD->_gdimensions;
   latt_size[nd-1] = 1;
 
   std::vector<int> simd_layout = GridDefaultSimd(nd-1, vComplex::Nsimd());
@@ -196,264 +196,23 @@ inline GridCartesian * MakeLowerDimGrid( GridCartesian * gridHD )
  Perambulator object
  ******************************************************************************/
 
-template<typename LatticeObj>
+template<typename Scalar_, int NumIndices_>
-class Perambulator : Serializable{
+class Perambulator : public Eigen::Tensor<Scalar_, NumIndices_, Eigen::RowMajor | Eigen::DontAlign>
-  // TODO: The next line makes friends across all combinations
+{
-  //     (not much of a problem given all public anyway ...)
-  //      FYI, the bug here was that I forgot that the friend is templated
-  template<typename T> friend std::ostream & operator<<(std::ostream &os, const Perambulator<T>& p);
-protected:
 public:
-  GRID_SERIALIZABLE_CLASS_MEMBERS( Perambulator,
+  template<typename... IndexTypes>
-                                  std::string,             ID,          // Allows owner to specialise
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Perambulator(Eigen::Index firstDimension, IndexTypes... otherDimensions)
-                                  std::string,             Provenance,  // For info only
+  : Eigen::Tensor<Scalar_, NumIndices_, Eigen::RowMajor | Eigen::DontAlign>(firstDimension, otherDimensions...)
-                                  std::vector<int>,        dimensions,
+  {
-                                  std::vector<LatticeObj>, perambulator,
+    // The number of dimensions used to construct a tensor must be equal to the rank of the tensor.
-                                  // Following items are redundant, but useful
+    EIGEN_STATIC_ASSERT(sizeof...(otherDimensions) + 1 == NumIndices_, YOU_MADE_A_PROGRAMMING_MISTAKE)
-                                  int,     nd,          // Number of dimensions
-                                  size_t,  NumElements); // Number of elements
-protected:
-  // Constructor common code
-  inline void ConstructCommon(const int * Dimensions) {
-    assert(nd > 0);
-    dimensions.resize(nd);
-    NumElements = 1;
-    for(int i = 0 ; i < nd ; i++) {
-      assert(Dimensions[i] > 0);
-      NumElements *= (size_t) Dimensions[i];
-      dimensions[i] = Dimensions[i];
-    }
-    //const LatticeObj perambulatorDefault;
-    perambulator.resize(NumElements);//,perambulatorDefault);
-  }
-public:
-  // Constructor with dimensions passed as std::vector<int>
-  inline Perambulator(const std::vector<int> & Dimensions)
-  : nd {(int) Dimensions.size()} {
-    ConstructCommon( &Dimensions[0] ); }
-  
-  // Constructor with dimensions passed as std::vector<int>
-  inline Perambulator(const std::vector<int> & Dimensions, const std::string sID)
-  : nd {(int) Dimensions.size()}, ID(sID) {
-    ConstructCommon( &Dimensions[0] ); }
-  
-  // Constructor with dimensions passed as std::vector<int>
-  inline Perambulator(const std::vector<int> & Dimensions, const std::string sID, const std::string sProvenance)
-  : nd {(int) Dimensions.size()}, ID(sID), Provenance(sProvenance) {
-    ConstructCommon( &Dimensions[0] ); }
-  
-  // Constructor with dimensions passed as individual parameters
-  // FYI: The caller is free to ignore the names and use the indices however they see fit
-  inline Perambulator(int NumNoise, int NumEvec=1, int NumTime=1, int NumSpin=1, int I_k=1, int I_t=1, int I_s=1) {
-    int Dimensions[]={NumNoise,NumEvec,NumTime,NumSpin,I_k,I_t,I_s};
-    nd = sizeof(Dimensions)/sizeof(Dimensions[0]);
-    while( nd > 1 && Dimensions[nd-1] == 1 )
-      nd--;
-    ConstructCommon( Dimensions );
   }
 
-  inline Perambulator(const std::string sID, int NumNoise, int NumEvec=1, int NumTime=1, int NumSpin=1, int I_k=1, int I_t=1, int I_s=1) : ID{sID} {
+  inline void WriteTemporary(const std::string &FileName){}
-    int Dimensions[]={NumNoise,NumEvec,NumTime,NumSpin,I_k,I_t,I_s};
-    nd = sizeof(Dimensions)/sizeof(Dimensions[0]);
-    while( nd > 1 && Dimensions[nd-1] == 1 )
-      nd--;
-    ConstructCommon( Dimensions );
-  }
-  
-  inline Perambulator(const std::string sID, const std::string sProvenance, int NumNoise, int NumEvec=1, int NumTime=1, int NumSpin=1, int I_k=1, int I_t=1, int I_s=1) : ID{sID}, Provenance{sProvenance} {
-    int Dimensions[]={NumNoise,NumEvec,NumTime,NumSpin,I_k,I_t,I_s};
-    nd = sizeof(Dimensions)/sizeof(Dimensions[0]);
-    while( nd > 1 && Dimensions[nd-1] == 1 )
-      nd--;
-    ConstructCommon( Dimensions );
-  }
-
-  inline LatticeObj & operator()(size_t count, const int * Coord) {
-    assert( count == nd );
-    assert( Coord );
-    size_t idx = 0;
-    // C memory order (???)
-    for( int d = 0 ; d < nd ; d++ ) {
-      assert( Coord[d] < dimensions[d] );
-      idx *= (size_t) dimensions[d];
-      idx += (size_t) Coord[d];
-    }
-    return perambulator[idx];
-  }
-  
-  inline LatticeObj & operator()(const std::vector<int> Coord) {
-    return operator()(Coord.size(), &Coord[0]);
-  }
-  
-  inline LatticeObj & operator()(int idxNoise, int idxEvec=0, int idxTime=0, int idxSpin=0, int I_k=0, int I_t=0, int I_s=0) {
-    int MyIndex[]={idxNoise,idxEvec,idxTime,idxSpin,I_k,I_t,I_s};
-    int i = sizeof(MyIndex)/sizeof(MyIndex[0]);
-    assert( i >= nd );
-    while( i > nd )
-      assert(MyIndex[--i] == 0);
-    return operator()(i, MyIndex);
-  }
 
   // Share data for timeslices we calculated with other nodes
   inline void SliceShare( GridCartesian * gridLowDim, GridCartesian * gridHighDim ) {
-    Grid::SliceShare( gridLowDim, gridHighDim, &perambulator[0],
+    //Grid::SliceShare( gridLowDim, gridHighDim, data(), (int) (size() * sizeof(Scalar_)));
-                                               (int) perambulator.size() * sizeof(perambulator[0]) );
-  }
-  
-  /*************************************************************************************
-   
-   Write/Read perambulator to/from disk
-   
-   Temporary version - keep the code running until such time as correct format written
-   
-   *************************************************************************************/
-  
-  inline void WriteTemporary(const std::string filename) const
-  {
-    std::cout << GridLogMessage << "Writing perambulator ID \"" << ID << "\" to " << filename << std::endl;
-    BinaryWriter myPhDThesis( filename + ".tmp" );
-    write( myPhDThesis, "Perambulator", *this );
-  }
-  
-  inline bool ReadTemporary(const std::string filename)
-  {
-    std::string _filename{filename};
-    _filename.append( ".tmp" );
-    bool bReturnValue = false;
-    std::fstream f;
-    f.open(_filename,std::ios_base::in);
-    if( !f.is_open() )
-      std::cout << GridLogMessage << "Cached perambulator file " << _filename << " does not exist" << std::endl;
-    else {
-      f.close();
-      bReturnValue = true;
-      auto MyID{ID};
-      std::cout << GridLogMessage << "Reading perambulator ID \"" << ID << "\" from " << _filename << std::endl;
-      BinaryReader reader( _filename );
-      read( reader, "Perambulator", *this );
-      std::cout << GridLogMessage << "Perambulator ID read from " << _filename << " was \"" << ID << "\"" << std::endl;
-      assert(MyID == ID);
-    }
-    return bReturnValue;
-  }
-  
-  /*************************************************************************************
-   
-   Write perambulator to disk
-   TODO 1. Ensure precision on disk can be specified independently of in memory
-   2. Validate format with Peter, Antonin et al
-   3. This object "works" for small lattii (lattices),
-   BUT, the final block is written as XML, and while write is fast, read is painfully slow
-   i.e. on a 24^3 x 64 lattice, I abandoned the perambulator read after 1 hour on Tesseract
-   
-   *************************************************************************************/
-  
-  inline void WritePoorly(LatticeGaugeField &field, const std::string filename)
-  {
-    std::cout << GridLogMessage << "Writing perambulator ID \"" << ID << "\" to " << filename << std::endl;
-    assert(nd>=2); // Really should be a little bigger
-    GridBase * gridHighDim = field._grid;
-    //ScidacWriterPerambulator binWriter(gridHighDim->IsBoss());
-    ScidacWriter binWriter(gridHighDim->IsBoss());
-    //makeFileDir(filename, gridHighDim); // Assume this makes directory ... but why pass it the grid?
-    binWriter.open(filename);
-    // Write the header
-    {
-      XmlWriter xmlWriter("", "perambulatorPar");
-      xmlWriter.pushXmlString("<ID>" + ID + "</ID>");
-      xmlWriter.pushXmlString("<Provenance>" + Provenance + "</Provenance>");
-      // TODO add all the perambulator parameters here
-      binWriter.writeLimeObject(1, 1, xmlWriter, "parameters", SCIDAC_FILE_XML);
-      std::cout << GridLogMessage << "Perambulator header written" << std::endl;
-    }
-    // Now write the local portion of the Perambulator
-    {
-      //binWriter.writeScidacPerambulatorRecord(field, *this);
-      //std::cout << GridLogMessage << "Perambulator body written" << std::endl;
-    }
-    {
-      ////////////////////////////////////////
-      // fill the Grid header
-      ////////////////////////////////////////
-      FieldMetaData header;
-      scidacRecord  _scidacRecord;
-      scidacFile    _scidacFile;
-      ScidacMetaData(field,header,_scidacRecord,_scidacFile);
-      
-      //////////////////////////////////////////////
-      // Fill the Lime file record by record
-      //////////////////////////////////////////////
-      constexpr auto precision = std::numeric_limits<Real>::digits10;
-      binWriter.writeLimeObject(1,0,header ,std::string("FieldMetaData"),std::string(GRID_FORMAT)); // Open message
-      binWriter.writeLimeObject(0, 0, _scidacRecord, _scidacRecord.SerialisableClassName(),
-                                std::string(SCIDAC_PRIVATE_RECORD_XML));
-      binWriter.writeLimeObject(0,1,*this,this->SerialisableClassName(),std::string("Perambulator"),precision);
-    }
-  }
-  
-  /*************************************************************************************
-   
-   Read perambulator from disk
-   TODO 1. Ensure precision on disk can be specified independently of in memory
-   2. Validate format with Peter
-   3. Abandoning for now because of synchronisation during write.
-   Object small enough to send to root and write from there
-   
-   *************************************************************************************/
-  
-  struct PerambHeader{
-    std::string ID, Provenance;
-  };
-  
-  inline bool ReadPoorly(LatticeGaugeField &field, const std::string filename)
-  {
-    assert(nd>=2); // Really should be a little bigger
-    bool bReturnValue = false;
-    std::fstream f;
-    f.open(filename,std::ios_base::in);
-    if( !f.is_open() )
-      std::cout << GridLogMessage << "Cached perambulator file " << filename << " does not exist" << std::endl;
-    else {
-      f.close();
-      ScidacReader binReader;
-      binReader.open(filename);
-      PerambHeader header;
-      // Read the header
-      {
-        std::string recordXml;
-        std::cout << GridLogMessage << "Reading perambulator header from " << filename << std::endl;
-        binReader.readLimeObject(recordXml, SCIDAC_FILE_XML);
-        XmlReader xmlReader(recordXml, true, "perambulatorPar");
-        xmlReader.push("perambulatorPar");
-        //xmlReader.readCurrentSubtree(header.ID);
-        xmlReader.readDefault("ID",header.ID);
-        std::cout << GridLogMessage << "Perambulator ID=" << header.ID << std::endl;
-        //xmlReader.nextElement();
-        //xmlReader.readCurrentSubtree(header.Provenance);
-        xmlReader.readDefault("Provenance",header.Provenance);
-        std::cout << GridLogMessage << "Perambulator Provenance=" << header.Provenance << std::endl;
-        assert( header.ID == ID );
-        bReturnValue = true;
-      }
-      // Now read the Perambulator
-      {
-        ////////////////////////////////////////
-        // fill the Grid header
-        ////////////////////////////////////////
-        FieldMetaData header;
-        scidacRecord  _scidacRecord;
-        
-        //////////////////////////////////////////////
-        // Fill the Lime file record by record
-        //////////////////////////////////////////////
-        binReader.readLimeObject(header ,std::string("FieldMetaData"),std::string(GRID_FORMAT)); // Open message
-        binReader.readLimeObject( _scidacRecord, _scidacRecord.SerialisableClassName(),
-                                 std::string(SCIDAC_PRIVATE_RECORD_XML));
-        binReader.readLimeObject(*this,this->SerialisableClassName(),std::string("Perambulator"));
-      }
-      // TODO Add validation that the field matches what we read in
-    }
-    return bReturnValue;
   }
 };
 

Hadrons/Modules/MDistil/DistilVectors.hpp

@@ -134,7 +134,7 @@ void TDistilVectors<FImpl>::execute(void)
    
     //auto        &noise     = envGet(std::vector<std::vector<std::vector<SpinVector>>>, par().noise);
     auto        &noise     = envGet(std::vector<Complex>, par().noise);
-    auto        &perambulator   = envGet(Perambulator<SpinVector>, par().perambulator);
+    auto        &perambulator   = *env().template getObject<Perambulator<SpinVector,6> >(par().perambulator);
     auto        &epack   = envGet(Grid::Hadrons::EigenPack<LatticeColourVector>, par().eigenPack);
     auto        &rho       = envGet(std::vector<FermionField>, getName() + "_rho");
     auto        &phi       = envGet(std::vector<FermionField>, getName() + "_phi");

Hadrons/Modules/MDistil/PerambLight.hpp

@@ -129,7 +129,8 @@ void TPerambLight<GImpl>::setup(void)
     const int nvec{par().nvec};
     const DistilParameters & Distil{par().Distil};
 
-    envCreate(Perambulator<SpinVector>, getName() + "_perambulator_light", 1,
+    //envCreate(Perambulator<SpinVector, 6>, getName() + "_perambulator_light", 1,
+    env().template createObject<Perambulator<SpinVector, 6> >(getName() + "_perambulator_light", Environment::Storage::object, 1,
               Distil.Nt,nvec,Distil.LI,Distil.nnoise,Distil.Nt_inv,Distil.SI);
     envCreate(std::vector<Complex>, getName() + "_noise", 1,
               nvec*Distil.Ns*Distil.Nt*Distil.nnoise);
@@ -183,7 +184,9 @@ void TPerambLight<GImpl>::execute(void)
 
     //auto        &noise     = envGet(std::vector<std::vector<std::vector<SpinVector>>>, par().noise);
     auto        &noise   = envGet(std::vector<Complex>, getName() + "_noise");
-    auto        &perambulator   = envGet(Perambulator<SpinVector>, getName() + "_perambulator_light");
+    auto        &perambulator = //envGet(Perambulator<SpinVector>,
+                                *env().template getObject<Perambulator<SpinVector,6> >(
+                                    getName() + "_perambulator_light");
     auto        &epack   = envGet(Grid::Hadrons::EigenPack<LatticeColourVector>, par().eigenPack);
 
     envGetTmp(GaugeField, Umu);