Added index names to perambulator

Hadrons/Modules/MDistil/Distil.hpp

@@ -34,6 +34,10 @@
 #include <Hadrons/Module.hpp>
 #include <Hadrons/ModuleFactory.hpp>
 
+#ifndef COMMA
+#define COMMA ,
+#endif
+
 /******************************************************************************
  This potentially belongs in CartesianCommunicator
  Turns out I don't actually need this when running inside hadrons
@@ -199,10 +203,13 @@ inline GridCartesian * MakeLowerDimGrid( GridCartesian * gridHD )
 template<typename Scalar_, int NumIndices_>
 class Perambulator : public Eigen::Tensor<Scalar_, NumIndices_, Eigen::RowMajor | Eigen::DontAlign>
 {
+protected:
+public:
+  std::array<std::string,NumIndices_> IndexNames;
 public:
   template<typename... IndexTypes>
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Perambulator(Eigen::Index firstDimension, IndexTypes... otherDimensions)
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Perambulator(std::array<std::string,NumIndices_> &IndexNames_, Eigen::Index firstDimension, IndexTypes... otherDimensions)
-  : Eigen::Tensor<Scalar_, NumIndices_, Eigen::RowMajor | Eigen::DontAlign>(firstDimension, otherDimensions...)
+  : IndexNames{IndexNames_}, Eigen::Tensor<Scalar_, NumIndices_, Eigen::RowMajor | Eigen::DontAlign>(firstDimension, otherDimensions...)
   {
     // The number of dimensions used to construct a tensor must be equal to the rank of the tensor.
     EIGEN_STATIC_ASSERT(sizeof...(otherDimensions) + 1 == NumIndices_, YOU_MADE_A_PROGRAMMING_MISTAKE)
@@ -212,7 +219,7 @@ public:
 
   // Share data for timeslices we calculated with other nodes
   inline void SliceShare( GridCartesian * gridLowDim, GridCartesian * gridHighDim ) {
-    //Grid::SliceShare( gridLowDim, gridHighDim, data(), (int) (size() * sizeof(Scalar_)));
+    Grid::SliceShare( gridLowDim, gridHighDim, this->data(), (int) (this->size() * sizeof(Scalar_)));
   }
 };
 

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   = *env().template getObject<Perambulator<SpinVector,6> >(par().perambulator);
+    auto        &perambulator   = envGet(Perambulator<SpinVector COMMA 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

@@ -127,10 +127,17 @@ void TPerambLight<GImpl>::setup(void)
     // auto &noise = envGet(std::vector<std::vector<std::vector<SpinVector>>>, par().noise);
     const int nvec{par().nvec};
     const DistilParameters & Distil{par().Distil};
+    //const char * IndexNames[6] = {"Nt", "nvec", "LI", "nnoise", "Nt_inv", "SI"};
+    std::array<std::string,6> sIndexNames{"Nt", "nvec", "LI", "nnoise", "Nt_inv", "SI"};
+    //std::complex<double> z{0.6,-3.1};
 
-    //envCreate(Perambulator<SpinVector, 6>, getName() + "_perambulator_light", 1,
+    //envCreate(std::string, getName() + "_debug_delete_me", 1, "Bingonuts");
-    env().template createObject<Perambulator<SpinVector, 6> >(getName() + "_perambulator_light", Environment::Storage::object, 1,
+    //envCreate(std::complex<double>, getName() + "_debug_delete_me_2", 1, 0.6);
-              Distil.Nt,nvec,Distil.LI,Distil.nnoise,Distil.Nt_inv,Distil.SI);
+    //envCreate(std::complex<double>, getName() + "_debug_delete_me_3", 1, z);
+    //envCreate(std::complex<double>, getName() + "_debug_delete_me_4", 1, {0.6 COMMA -3.1});
+    //envCreate(std::array<std::string COMMA 3>, getName() + "_debug_delete_me_5", 1, {"One" COMMA "Two" COMMA "Three"});
+    envCreate(Perambulator<SpinVector COMMA 6>, getName() + "_perambulator_light", 1,
+              sIndexNames,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,9 +190,7 @@ 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>,
+    auto        &perambulator = envGet(Perambulator<SpinVector COMMA 6>, getName() + "_perambulator_light");
-                                *env().template getObject<Perambulator<SpinVector,6> >(
-                                    getName() + "_perambulator_light");
     auto        &epack   = envGet(Grid::Hadrons::EigenPack<LatticeColourVector>, par().eigenPack);
 
     envGetTmp(GaugeField, Umu);
@@ -252,9 +257,10 @@ void TPerambLight<GImpl>::execute(void)
             noise[inoise + nnoise*(t + Nt*(ivec+nvec*is))] = 1.;
             //noises[inoise][t][ivec]()(is)() = 1.;
           else{
-            random(sRNG,rn); 
+            random(sRNG,rn);
-            noise[inoise + nnoise*(t + Nt*(ivec+nvec*is))] = (rn-0.5 > 0) - (rn-0.5 < 0); //TODO: This could be 0 if rn==0.5!!
+            // We could use a greater number of complex roots of unity
-            //noises[inoise][t][ivec]()(is)() = (rn-0.5 > 0) - (rn-0.5 < 0); //TODO: This could be 0 if rn==0.5!!
+            // ... but this seems to work well
+            noise[inoise + nnoise*(t + Nt*(ivec+nvec*is))] = (rn > 0.5) ? -1 : 1;
           }
         }
       }

tests/hadrons/Test_hadrons_distil.cc

@@ -237,9 +237,8 @@ bool bNumber( int &ri, const char * & pstr, bool bGobbleWhiteSpace = true )
 
 typedef Eigen::Tensor<Complex,3,Eigen::RowMajor | Eigen::DontAlign> MyTensor;
 
-bool DebugEigenTest()
+void DebugShowTensor(MyTensor &x)
 {
-  MyTensor x(2,3,4);
   const MyTensor::Index s{x.size()};
   std::cout << "x.size() = " << s << std::endl;
   std::cout << "x.NumDimensions = " << x.NumDimensions << " (TensorBase)" << std::endl;
@@ -261,11 +260,11 @@ bool DebugEigenTest()
   // Initialise
   assert( d.size() == 3 );
   for( int i = 0 ; i < d[0] ; i++ )
-    for( int j = 0 ; j < d[1] ; j++ )
+  for( int j = 0 ; j < d[1] ; j++ )
-      for( int k = 0 ; k < d[2] ; k++ ) {
+  for( int k = 0 ; k < d[2] ; k++ ) {
-        x(i,j,k) = std::complex<double>(SizeCalculated, SizeCalculated);
+    x(i,j,k) = std::complex<double>(SizeCalculated, SizeCalculated);
-        SizeCalculated--;
+    SizeCalculated--;
-      }
+  }
   // Show raw data
   std::cout << "Data follow : " << std::endl;
   Complex * p = x.data();
@@ -274,6 +273,21 @@ bool DebugEigenTest()
     std::cout << "x.data()[" << i << "]=" << * p++;
   }
   std::cout << std::endl;
+}
+
+bool DebugEigenTest()
+{
+  MyTensor x(2,3,4);
+  DebugShowTensor(x);
+  // Test initialisation of an array of strings
+  std::array<std::string,3> as={"Alpha", "Beta", "Gamma"};
+  for( auto a : as )
+    std::cout << a << std::endl;
+  Grid::Hadrons::MDistil::Perambulator<Complex,3> p{as,2,7,2};
+  DebugShowTensor(p);
+  std::cout << "p.IndexNames follow" << std::endl;
+  for( auto a : p.IndexNames )
+    std::cout << a << std::endl;
   return true;
 }
 #endif
@@ -282,7 +296,7 @@ int main(int argc, char *argv[])
 {
 #ifdef DEBUG
   // Debug only - test of Eigen::Tensor
-  //if( DebugEigenTest() ) return 0;
+  if( DebugEigenTest() ) return 0;
 #endif
 
   // Decode command-line parameters. 1st one is which test to run