WIP

Hadrons/Distil.hpp

@@ -47,9 +47,7 @@
 /******************************************************************************
  A consistent set of cross-platform methods for big endian <-> host byte ordering
  I imagine this exists already?
- 
+ This can be removed once the (deprecated) NamedTensor::ReadBinary & WriteBinary methods are deleted
- RELEASE NOTE:
- 
  ******************************************************************************/
 
 #if defined(__linux__)
@@ -81,7 +79,6 @@
 
 /******************************************************************************
  This potentially belongs in CartesianCommunicator
- Turns out I don't actually need this when running inside hadrons
  ******************************************************************************/
 
 BEGIN_MODULE_NAMESPACE(Grid)
@@ -138,15 +135,13 @@ inline void SliceShare( GridBase * gridLowDim, GridBase * gridHighDim, void * Bu
 }
 
 /*************************************************************************************
+
+ Not sure where the right home for this is? But presumably in Grid
  
- -Grad^2 (Peardon, 2009, pg 2, equation 3)
+ -Grad^2 (Peardon, 2009, pg 2, equation 3, https://arxiv.org/abs/0905.2160)
  Field      Type of field the operator will be applied to
  GaugeField Gauge field the operator will smear using
  
- TODO CANDIDATE for integration into laplacian operator
- should just require adding number of dimensions to act on to constructor,
- where the default=all dimensions, but we could specify 3 spatial dimensions
- 
  *************************************************************************************/
 
 template<typename Field, typename GaugeField=LatticeGaugeField>
@@ -202,15 +197,6 @@ public:
   }
 };
 
-
-class BFieldIO: Serializable{
-public:
-  using BaryonTensorSet = Eigen::Tensor<Complex, 7>;
-  GRID_SERIALIZABLE_CLASS_MEMBERS(BFieldIO,
-                                  BaryonTensorSet, BField
-		                  );
-};
-
 END_MODULE_NAMESPACE // Grid
 
 /******************************************************************************
@@ -221,29 +207,33 @@ BEGIN_HADRONS_NAMESPACE
 
 BEGIN_MODULE_NAMESPACE(MDistil)
 
-typedef Grid::Hadrons::EigenPack<LatticeColourVector>       DistilEP;
+using DistilEP = Grid::Hadrons::EigenPack<LatticeColourVector>;
-typedef std::vector<std::vector<std::vector<SpinVector> > > DistilNoises;
+
+// Noise vector index order: nnoise, nt, nvec, ns
+using NoiseTensor = Eigen::Tensor<Complex, 4, Eigen::RowMajor>;
+
+struct DistilParameters: Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(DistilParameters,
+                                  int, TI,
+                                  int, LI,
+                                  int, nnoise,
+                                  int, tsrc,
+                                  int, SI,
+                                  int, Ns,
+                                  int, Nt_inv)
+  DistilParameters() = default;
+  template <class ReaderClass> DistilParameters(Reader<ReaderClass>& Reader){read(Reader,"Distil",*this);}
+};
+
+class BFieldIO: Serializable{
+public:
+  using BaryonTensorSet = Eigen::Tensor<Complex, 7>;
+  GRID_SERIALIZABLE_CLASS_MEMBERS(BFieldIO, BaryonTensorSet, BField );
+};
 
 /******************************************************************************
- Make a lower dimensional grid
+ Default for distillation file operations. For now only used by NamedTensor
- ******************************************************************************/
+******************************************************************************/
-
-inline GridCartesian * MakeLowerDimGrid( GridCartesian * gridHD )
-{
-  //LOG(Message) << "MakeLowerDimGrid() begin" << std::endl;
-  int nd{static_cast<int>(gridHD->_ndimension)};
-  std::vector<int> latt_size   = gridHD->_gdimensions;
-  latt_size[nd-1] = 1;
-
-  std::vector<int> simd_layout = GridDefaultSimd(nd-1, vComplex::Nsimd());
-  simd_layout.push_back( 1 );
-
-  std::vector<int> mpi_layout  = gridHD->_processors;
-  mpi_layout[nd-1] = 1;
-  GridCartesian * gridLD = new GridCartesian(latt_size,simd_layout,mpi_layout,*gridHD);
-  //LOG(Message) << "MakeLowerDimGrid() end" << std::endl;
-  return gridLD;
-}
 
 #ifdef HAVE_HDF5
 using Default_Reader = Grid::Hdf5Reader;
@@ -348,6 +338,14 @@ template<typename T, typename V = void> struct is_named_tensor : public std::fal
 template<typename Scalar_, int NumIndices_, uint16_t Endian_Scalar_Size_> struct is_named_tensor<NamedTensor<Scalar_, NumIndices_, Endian_Scalar_Size_>> : public std::true_type {};
 template<typename T> struct is_named_tensor<T, typename std::enable_if<std::is_base_of<NamedTensor<typename T::Scalar, T::NumIndices, T::Endian_Scalar_Size_>, T>::value>::type> : public std::true_type {};
 
+/******************************************************************************
+ Perambulator object
+ Endian_Scalar_Size can be removed once (deprecated) NamedTensor::ReadBinary & WriteBinary methods deleted
+ ******************************************************************************/
+
+template<typename Scalar_, int NumIndices_, uint16_t Endian_Scalar_Size = sizeof(Scalar_)>
+using Perambulator = NamedTensor<Scalar_, NumIndices_, Endian_Scalar_Size>;
+
 /******************************************************************************
  Save NamedTensor binary format (NB: On-disk format is Big Endian)
  Assumes the Scalar_ objects are contiguous (no padding)
@@ -586,31 +584,25 @@ void NamedTensor<Scalar_, NumIndices_, Endian_Scalar_Size>::read(const char * fi
 }
 
 /******************************************************************************
- Perambulator object
+ Make a lower dimensional grid in preparation for local slice operations
  ******************************************************************************/
 
-template<typename Scalar_, int NumIndices_, uint16_t Endian_Scalar_Size = sizeof(Scalar_)>
+inline GridCartesian * MakeLowerDimGrid( GridCartesian * gridHD )
-using Perambulator = NamedTensor<Scalar_, NumIndices_, Endian_Scalar_Size>;
+{
-
+  int nd{static_cast<int>(gridHD->_ndimension)};
-struct DistilParameters: Serializable {
+  std::vector<int> latt_size   = gridHD->_gdimensions;
-  GRID_SERIALIZABLE_CLASS_MEMBERS(DistilParameters,
+  latt_size[nd-1] = 1;
-                                  int, TI,
+  std::vector<int> simd_layout = GridDefaultSimd(nd-1, vComplex::Nsimd());
-                                  int, LI,
+  simd_layout.push_back( 1 );
-                                  int, nnoise,
+  std::vector<int> mpi_layout  = gridHD->_processors;
-                                  int, tsrc,
+  mpi_layout[nd-1] = 1;
-                                  int, SI,
+  GridCartesian * gridLD = new GridCartesian(latt_size,simd_layout,mpi_layout,*gridHD);
-                                  int, Ns,
+  return gridLD;
-                                  int, Nt,
+}
-                                  int, Nt_inv)
-  DistilParameters() = default;
-  template <class ReaderClass> DistilParameters(Reader<ReaderClass>& Reader){read(Reader,"Distil",*this);}
-};
 
 /*************************************************************************************
- 
  Rotate eigenvectors into our phase convention
  First component of first eigenvector is real and positive
- 
  *************************************************************************************/
 
 inline void RotateEigen(std::vector<LatticeColourVector> & evec)
@@ -642,7 +634,5 @@ inline void RotateEigen(std::vector<LatticeColourVector> & evec)
 }
 
 END_MODULE_NAMESPACE
-
 END_HADRONS_NAMESPACE
-
 #endif // Hadrons_MDistil_Distil_hpp_

Hadrons/Modules/MDistil/DistilVectors.hpp

@@ -60,13 +60,11 @@ public:
                                   std::string, sink,
                                   bool, multiFile,
                                   int, tsrc,
-                                  //int, nnoise,
                                   int, LI,
                                   int, SI,
                                   int, TI,
                                   int, nvec,
                                   int, Ns,
-                                  int, Nt,
                                   int, Nt_inv);
 };
 
@@ -172,7 +170,7 @@ template <typename FImpl>
 void TDistilVectors<FImpl>::setup(void)
 {
   Cleanup();
-  auto &noise        = envGet(std::vector<Complex>, NoiseVectorName);
+  auto &noise        = envGet(NoiseTensor, NoiseVectorName);
   auto &perambulator = envGet(DistilPeramb, PerambulatorName);
 
   // We expect the perambulator to have been created with these indices
@@ -180,7 +178,8 @@ void TDistilVectors<FImpl>::setup(void)
   for(int i = 0; i < DistilPeramb::NumIndices; i++ )
     assert( sIndexNames[i] == perambulator.IndexNames[i] && "Perambulator indices bad" );
 
-  nnoise = static_cast<int>( noise.size() );
+  nnoise = static_cast<int>( noise.dimension(0) );
+  LOG(Message) << "NoiseTensor has " << nnoise << " noise vectors" << std::endl;
   int LI=par().LI;
   int Ns=par().Ns;
   int SI=par().SI;
@@ -225,7 +224,7 @@ void TDistilVectors<FImpl>::Cleanup(void)
 template <typename FImpl>
 void TDistilVectors<FImpl>::execute(void)
 {
-  auto &noise        = envGet(std::vector<Complex>, NoiseVectorName);
+  auto &noise        = envGet(NoiseTensor, NoiseVectorName);
   auto &perambulator = envGet(DistilPeramb, PerambulatorName);
   auto &epack        = envGet(Grid::Hadrons::EigenPack<LatticeColourVector>, LapEvecName);
   
@@ -244,7 +243,7 @@ void TDistilVectors<FImpl>::execute(void)
   int LI=par().LI;
   int Ns=par().Ns;
   int Nt_inv=par().Nt_inv; // TODO: No input, but define through Nt, TI
-  int Nt=par().Nt;
+  const int Nt{grid4d->GlobalDimensions()[Tdir]};
   int TI=par().TI;
   int nvec=par().nvec;
   int SI=par().SI;
@@ -255,10 +254,10 @@ void TDistilVectors<FImpl>::execute(void)
   int t_inv;
   if( bMakeSource ) {
     auto &rho = envGet(std::vector<FermionField>, SourceName);
-    for (int inoise = 0; inoise < nnoise; inoise++) {
+    for( int inoise = 0; inoise < nnoise; inoise++ ) {
-      for (int dk = 0; dk < LI; dk++) {
+      for( int dk = 0; dk < LI; dk++ ) {
-        for (int dt = 0; dt < Nt_inv; dt++) {
+        for( int dt = 0; dt < Nt_inv; dt++ ) {
-          for (int ds = 0; ds < SI; ds++) {
+          for( int ds = 0; ds < SI; ds++ ) {
             vecindex = inoise + nnoise * dk + nnoise * LI * ds + nnoise *LI * SI*dt;
             rho[vecindex] = zero;
             tmp3d_nospin = zero;
@@ -268,7 +267,8 @@ void TDistilVectors<FImpl>::execute(void)
                 for (int ik = dk; ik < nvec; ik += LI){
                   for (int is = ds; is < Ns; is += SI){
                     ExtractSliceLocal(evec3d,epack.evec[ik],0,t_inv,3);
-                    tmp3d_nospin = evec3d * noise[inoise + nnoise*(t_inv + Nt*(ik+nvec*is))];
+                    //tmp3d_nospin = evec3d * noise[inoise + nnoise*(t_inv + Nt*(ik+nvec*is))];
+                    tmp3d_nospin = evec3d * noise(inoise, t_inv, ik, is);
                     tmp3d=zero;
                     pokeSpin(tmp3d,tmp3d_nospin,is);
                     tmp2=zero;
@@ -285,10 +285,10 @@ void TDistilVectors<FImpl>::execute(void)
   }
   if( bMakeSink ) {
     auto &phi = envGet(std::vector<FermionField>, SinkName);
-    for (int inoise = 0; inoise < nnoise; inoise++) {
+    for( int inoise = 0; inoise < nnoise; inoise++ ) {
-      for (int dk = 0; dk < LI; dk++) {
+      for( int dk = 0; dk < LI; dk++ ) {
-        for (int dt = 0; dt < Nt_inv; dt++) {
+        for( int dt = 0; dt < Nt_inv; dt++ ) {
-          for (int ds = 0; ds < SI; ds++) {
+          for( int ds = 0; ds < SI; ds++ ) {
             vecindex = inoise + nnoise * dk + nnoise * LI * ds + nnoise *LI * SI*dt;
             phi[vecindex] = zero;
             for (int t = Ntfirst; t < Ntfirst + Ntlocal; t++) {
@@ -307,7 +307,5 @@ void TDistilVectors<FImpl>::execute(void)
 }
 
 END_MODULE_NAMESPACE
-
 END_HADRONS_NAMESPACE
-
 #endif // Hadrons_MDistil_DistilVectors_hpp_

Hadrons/Modules/MDistil/Noises.hpp

@@ -106,54 +106,53 @@ std::vector<std::string> TNoises<FImpl>::getOutput(void)
 template <typename FImpl>
 void TNoises<FImpl>::setup(void)
 {
-    const DistilParameters & Distil{par().Distil};
+  const int Nt{env().getGrid()->GlobalDimensions()[Tdir]};
-    const int nvec{par().nvec};
+  const int nvec{par().nvec};
- 
+  const DistilParameters & Distil{par().Distil};
-    envCreate(std::vector<Complex>, getName(), 1, nvec*Distil.Ns*Distil.Nt*Distil.nnoise);
+  //envCreate(std::vector<Complex>, getName(), 1, nvec*Distil.Ns*Distil.Nt*Distil.nnoise);
-
+  envCreate(NoiseTensor, getName(), 1, Distil.nnoise, Nt, nvec, Distil.Ns);
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 template <typename FImpl>
 void TNoises<FImpl>::execute(void)
 {
-       const std::string &UniqueIdentifier{par().UniqueIdentifier};
+  const std::string &UniqueIdentifier{par().UniqueIdentifier};
-       auto &noise   = envGet(std::vector<Complex>, getName());
+  auto &noise   = envGet(NoiseTensor, getName());
-       const int nvec{par().nvec};
+  const int nvec{par().nvec};
-       const DistilParameters & Distil{par().Distil};
+  const DistilParameters & Distil{par().Distil};
-       const int nnoise{Distil.nnoise};
+  const int nnoise{Distil.nnoise};
-       const int Nt{Distil.Nt};
+  const int Nt{env().getGrid()->GlobalDimensions()[Tdir]};
-       const int Ns{Distil.Ns};
+  const int Ns{Distil.Ns};
-       const int TI{Distil.TI};
+  const int TI{Distil.TI};
-       const int LI{Distil.LI};
+  const int LI{Distil.LI};
-       const bool full_tdil{TI==Nt};
+  const bool full_tdil{TI==Nt};
-       const bool exact_distillation{full_tdil && LI==nvec};
+  const bool exact_distillation{full_tdil && LI==nvec};
-
+  
-       GridSerialRNG sRNG; 
+  GridSerialRNG sRNG;
-       sRNG.SeedUniqueString(UniqueIdentifier + std::to_string(vm().getTrajectory())); //maybe add more??
+  sRNG.SeedUniqueString(UniqueIdentifier + std::to_string(vm().getTrajectory())); //maybe add more??
-       Real rn;
+  Real rn;
-		       
+  
-       for (int inoise=0;inoise<nnoise;inoise++) {
+  for( int inoise = 0; inoise < nnoise; inoise++ ) {
-         for (int t=0;t<Nt;t++) {
+    for( int t = 0; t < Nt; t++ ) {
-           for (int ivec=0;ivec<nvec;ivec++) {
+      for( int ivec = 0; ivec < nvec; ivec++ ) {
-             for (int is=0;is<Ns;is++) {
+        for( int is = 0; is < Ns; is++ ) {
-               if (exact_distillation)
+          if( exact_distillation )
-                 noise[inoise + nnoise*(t + Nt*(ivec+nvec*is))] = 1.;
+            //noise[inoise + nnoise*(t + Nt*(ivec+nvec*is))] = 1.;
-               else{
+            noise(inoise, t, ivec, is) = 1.;
-                 random(sRNG,rn);
+          else{
-                 // We could use a greater number of complex roots of unity
+            random(sRNG,rn);
-                 // ... but this seems to work well
+            // We could use a greater number of complex roots of unity
-                 noise[inoise + nnoise*(t + Nt*(ivec+nvec*is))] = (rn > 0.5) ? -1 : 1;
+            // ... but this seems to work well
-               }
+            //noise[inoise + nnoise*(t + Nt*(ivec+nvec*is))] = (rn > 0.5) ? -1 : 1;
-	     }
+            noise(inoise, t, ivec, is) = (rn > 0.5) ? -1 : 1;
-	   }
+          }
-	 }
+        }
-       }
+      }
-
+    }
+  }
 }
 
 END_MODULE_NAMESPACE
-
 END_HADRONS_NAMESPACE
-
 #endif // Hadrons_MDistil_Noises_hpp_

Hadrons/Modules/MDistil/PerambFromSolve.hpp

@@ -154,6 +154,7 @@ void TPerambFromSolve<FImpl>::setup(void)
 
     grid4d = env().getGrid();
         grid3d = MakeLowerDimGrid(grid4d);
+  const int Nt{grid4d->GlobalDimensions()[Tdir]};
 
 
 	const int nvec_reduced{par().nvec_reduced};
@@ -161,9 +162,9 @@ void TPerambFromSolve<FImpl>::setup(void)
     //envCreate(Perambulator<SpinVector COMMA 6 COMMA sizeof(Real)>, getName(), 1,
     //          sIndexNames,Distil.Nt,nvec,Distil.LI,Distil.nnoise,Distil.Nt_inv,Distil.SI);
     envCreate(Perambulator<SpinVector COMMA 6 COMMA sizeof(Real)>, getName(), 1,
-              sIndexNames,Distil.Nt,nvec_reduced,LI_reduced,Distil.nnoise,Distil.Nt_inv,Distil.SI);
+              sIndexNames,Nt,nvec_reduced,LI_reduced,Distil.nnoise,Distil.Nt_inv,Distil.SI);
     envCreate(std::vector<Complex>, getName() + "_noise", 1,
-              nvec*Distil.Ns*Distil.Nt*Distil.nnoise);
+              nvec*Distil.Ns*Nt*Distil.nnoise);
 
     envTmp(LatticeColourVector, "result_3d",1,LatticeColourVector(grid3d));
         envTmp(LatticeColourVector, "evec3d",1,LatticeColourVector(grid3d));
@@ -187,34 +188,32 @@ void TPerambFromSolve<FImpl>::Cleanup(void)
 template <typename FImpl>
 void TPerambFromSolve<FImpl>::execute(void)
 {
-	const int nvec_reduced{par().nvec_reduced};
+  GridCartesian * grid4d = env().getGrid();
-	const int LI_reduced{par().LI_reduced};
+  const int Nt{grid4d->GlobalDimensions()[Tdir]};
-    const int nvec{par().nvec};
+  const int Ntlocal{grid4d->LocalDimensions()[3]};
-    const DistilParameters & Distil{par().Distil};
+  const int Ntfirst{grid4d->LocalStarts()[3]};
-    const int LI{Distil.LI};
+
-       const int TI{Distil.TI};
+  const int nvec_reduced{par().nvec_reduced};
-       const int SI{Distil.SI};
+  const int LI_reduced{par().LI_reduced};
-    const int nnoise{Distil.nnoise};
+  const int nvec{par().nvec};
-    const int Nt{Distil.Nt};
+  const DistilParameters & Distil{par().Distil};
-    const int Nt_inv{Distil.Nt_inv};
+  const int LI{Distil.LI};
-    const int tsrc{Distil.tsrc};
+  const int TI{Distil.TI};
-    const int Ns{Distil.Ns};
+  const int SI{Distil.SI};
-    const bool full_tdil{TI==Nt};
+  const int nnoise{Distil.nnoise};
-    const bool exact_distillation{full_tdil && LI==nvec};
+  const int Nt_inv{Distil.Nt_inv};
-    auto        &perambulator = envGet(Perambulator<SpinVector COMMA 6 COMMA sizeof(Real)>,
+  const int tsrc{Distil.tsrc};
-                                       getName());
+  const int Ns{Distil.Ns};
- auto        &solve       = envGet(std::vector<FermionField>, par().solve);
+  const bool full_tdil{TI==Nt};
-    auto        &epack   = envGet(Grid::Hadrons::EigenPack<LatticeColourVector>, par().eigenPack);
+  const bool exact_distillation{full_tdil && LI==nvec};
+  auto &perambulator = envGet(Perambulator<SpinVector COMMA 6 COMMA sizeof(Real)>, getName());
+  auto &solve       = envGet(std::vector<FermionField>, par().solve);
+  auto &epack   = envGet(Grid::Hadrons::EigenPack<LatticeColourVector>, par().eigenPack);
 
   envGetTmp(LatticeColourVector, result_nospin);
   envGetTmp(LatticeColourVector, result_3d);
   envGetTmp(LatticeColourVector, evec3d);
 
-  GridCartesian * grid4d = env().getGrid();
-
-  int Ntlocal = grid4d->LocalDimensions()[3];
-  int Ntfirst = grid4d->LocalStarts()[3];
-
   for (int inoise = 0; inoise < nnoise; inoise++) {
     for (int dk = 0; dk < LI_reduced; dk++) {
       for (int dt = 0; dt < Nt_inv; dt++) {

Hadrons/Modules/MDistil/Perambulator.hpp

@@ -134,38 +134,36 @@ std::vector<std::string> TPerambulator<FImpl>::getOutput(void)
 template <typename FImpl>
 void TPerambulator<FImpl>::setup(void)
 {
-    Cleanup();
+  Cleanup();
-
+  grid4d = env().getGrid();
-    const int nvec{par().nvec};
+  grid3d = MakeLowerDimGrid(grid4d);
-    const DistilParameters & Distil{par().Distil};
+  const int Nt{grid4d->GlobalDimensions()[Tdir]};
-    const int LI{Distil.LI};
+  const int nvec{par().nvec};
-    const int nnoise{Distil.nnoise};
+  const DistilParameters & Distil{par().Distil};
-    const int Nt_inv{Distil.Nt_inv}; // TODO: PROBABLY BETTER: if (full_tdil) Nt_inv=1; else Nt_inv = TI;
+  const int LI{Distil.LI};
-    const int Ns{Distil.Ns};
+  const int nnoise{Distil.nnoise};
-    std::array<std::string,6> sIndexNames{"Nt", "nvec", "LI", "nnoise", "Nt_inv", "SI"};
+  const int Nt_inv{Distil.Nt_inv}; // TODO: PROBABLY BETTER: if (full_tdil) Nt_inv=1; else Nt_inv = TI;
-
+  const int Ns{Distil.Ns};
-    envCreate(Perambulator<SpinVector COMMA 6 COMMA sizeof(Real)>, getName(), 1,
+  std::array<std::string,6> sIndexNames{"Nt", "nvec", "LI", "nnoise", "Nt_inv", "SI"};
-              sIndexNames,Distil.Nt,nvec,Distil.LI,Distil.nnoise,Distil.Nt_inv,Distil.SI);
+  
-    envCreate(std::vector<FermionField>, getName() + "_unsmeared_sink", 1, 
+  envCreate(Perambulator<SpinVector COMMA 6 COMMA sizeof(Real)>, getName(), 1,
+            sIndexNames,Nt,nvec,Distil.LI,Distil.nnoise,Distil.Nt_inv,Distil.SI);
+  envCreate(std::vector<FermionField>, getName() + "_unsmeared_sink", 1,
             nnoise*LI*Ns*Nt_inv, envGetGrid(FermionField));
-
+  
-    grid4d = env().getGrid();
+  envTmpLat(LatticeSpinColourVector, "dist_source");
-    grid3d = MakeLowerDimGrid(grid4d);//new GridCartesian(latt_size,simd_layout_3,mpi_layout,*grid4d);
+  envTmpLat(LatticeSpinColourVector, "tmp2");
-
+  envTmpLat(LatticeSpinColourVector, "result");
-    envTmpLat(LatticeSpinColourVector, "dist_source");
+  envTmpLat(LatticeColourVector, "result_nospin");
-    envTmpLat(LatticeSpinColourVector, "tmp2");
+  envTmp(LatticeSpinColourVector, "tmp3d",1,LatticeSpinColourVector(grid3d));
-    envTmpLat(LatticeSpinColourVector, "result");
+  envTmp(LatticeColourVector, "tmp3d_nospin",1,LatticeColourVector(grid3d));
-    envTmpLat(LatticeColourVector, "result_nospin");
+  envTmp(LatticeColourVector, "result_3d",1,LatticeColourVector(grid3d));
-    envTmp(LatticeSpinColourVector, "tmp3d",1,LatticeSpinColourVector(grid3d));
+  envTmp(LatticeColourVector, "evec3d",1,LatticeColourVector(grid3d));
-    envTmp(LatticeColourVector, "tmp3d_nospin",1,LatticeColourVector(grid3d));
+  
-    envTmp(LatticeColourVector, "result_3d",1,LatticeColourVector(grid3d));
+  Ls_ = env().getObjectLs(par().solver);
-    envTmp(LatticeColourVector, "evec3d",1,LatticeColourVector(grid3d));
+  envTmpLat(FermionField, "v4dtmp");
-
+  envTmpLat(FermionField, "v5dtmp", Ls_);
-    Ls_ = env().getObjectLs(par().solver);
+  envTmpLat(FermionField, "v5dtmp_sol", Ls_);
-    envTmpLat(FermionField, "v4dtmp");
-    envTmpLat(FermionField, "v5dtmp", Ls_);
-    envTmpLat(FermionField, "v5dtmp_sol", Ls_);
-
 }
 
 // clean up any temporaries created by setup (that aren't stored in the environment)
@@ -183,30 +181,29 @@ void TPerambulator<FImpl>::Cleanup(void)
 template <typename FImpl>
 void TPerambulator<FImpl>::execute(void)
 {
+  const int Nt{grid4d->GlobalDimensions()[Tdir]};
     const int nvec{par().nvec};
     const DistilParameters & Distil{par().Distil};
     const int LI{Distil.LI};
     const int SI{Distil.SI};
     const int TI{Distil.TI};
     const int nnoise{Distil.nnoise};
-    const int Nt{Distil.Nt};
-    const int Nt_inv{Distil.Nt_inv}; // TODO: PROBABLY BETTER: if (full_tdil) Nt_inv=1; else Nt_inv = TI;
     const int tsrc{Distil.tsrc};
     const int Ns{Distil.Ns};
-  
+  const bool full_tdil{TI==Nt};
+  const bool exact_distillation{full_tdil && LI==nvec};
+  const int Nt_inv{full_tdil ? 1 : TI}; // TODO: PROBABLY BETTER: if (full_tdil) Nt_inv=1; else Nt_inv = TI;
+
     auto &solver=envGet(Solver, par().solver);
     auto &mat = solver.getFMat();
     envGetTmp(FermionField, v4dtmp);
     envGetTmp(FermionField, v5dtmp);
     envGetTmp(FermionField, v5dtmp_sol);
 
-
-    const bool full_tdil{TI==Nt};
-    const bool exact_distillation{full_tdil && LI==nvec};
-
     const std::string &UniqueIdentifier{par().UniqueIdentifier};
 
-    auto &noise = envGet(std::vector<Complex>, par().noise);
+    //auto &noise = envGet(std::vector<Complex>, par().noise);
+    auto &noise = envGet(NoiseTensor, par().noise);
     auto &perambulator = envGet(Perambulator<SpinVector COMMA 6 COMMA sizeof(Real)>, getName());
     auto &epack = envGet(Grid::Hadrons::EigenPack<LatticeColourVector>, par().eigenPack);
     auto &unsmeared_sink = envGet(std::vector<FermionField>, getName() + "_unsmeared_sink");
@@ -245,7 +242,8 @@ void TPerambulator<FImpl>::execute(void)
                 for (int ik = dk; ik < nvec; ik += LI){
                   for (int is = ds; is < Ns; is += SI){ 
                     ExtractSliceLocal(evec3d,epack.evec[ik],0,t_inv,3);
-                    tmp3d_nospin = evec3d * noise[inoise + nnoise*(t_inv + Nt*(ik+nvec*is))]; 
+                    //tmp3d_nospin = evec3d * noise[inoise + nnoise*(t_inv + Nt*(ik+nvec*is))];
+                    tmp3d_nospin = evec3d * noise(inoise, t_inv, ik, is);
                     tmp3d=zero;
                     pokeSpin(tmp3d,tmp3d_nospin,is);
                     tmp2=zero;

Hadrons/Modules/MIO/LoadPerambulator.hpp

@@ -101,17 +101,13 @@ std::vector<std::string> TLoadPerambulator<FImpl>::getOutput(void)
 template <typename FImpl>
 void TLoadPerambulator<FImpl>::setup(void)
 {
+  GridCartesian * grid4d = env().getGrid();
+  const int Nt{grid4d->GlobalDimensions()[Tdir]};
     const int nvec{par().nvec};
     const MDistil::DistilParameters & Distil{par().Distil};
-    const int LI{Distil.LI};
-    const int nnoise{Distil.nnoise};
-    const int Nt_inv{Distil.Nt_inv}; // TODO: PROBABLY BETTER: if (full_tdil) Nt_inv=1; else Nt_inv = TI;
-    const int Ns{Distil.Ns};
     std::array<std::string,6> sIndexNames{"Nt", "nvec", "LI", "nnoise", "Nt_inv", "SI"};
-
     envCreate(MDistil::Perambulator<SpinVector COMMA 6 COMMA sizeof(Real)>, getName(), 1,
-		              sIndexNames,Distil.Nt,nvec,Distil.LI,Distil.nnoise,Distil.Nt_inv,Distil.SI);
+		              sIndexNames,Nt,nvec,Distil.LI,Distil.nnoise,Distil.Nt_inv,Distil.SI);
-
 }
 
 // execution ///////////////////////////////////////////////////////////////////

tests/hadrons/Test_distil.cc

@@ -68,25 +68,22 @@ void test_Global(Application &application)
 
 void test_SolverS(Application &application)
 {
-    std::string boundary = "1 1 1 -1";
+  std::string boundary = "1 1 1 -1";
-
+  MAction::DWF::Par actionPar;
-
+  actionPar.gauge = "gauge";
-        MAction::DWF::Par actionPar;
+  actionPar.Ls    = 16;
-        actionPar.gauge = "gauge";
+  actionPar.M5    = 1.8;
-        actionPar.Ls    = 16;
+  actionPar.mass  = 0.005;
-        actionPar.M5    = 1.8;
+  actionPar.boundary = boundary;
-        actionPar.mass  = 0.005;
+  actionPar.twist = "0. 0. 0. 0.";
-        actionPar.boundary = boundary;
+  application.createModule<MAction::DWF>("DWF_s", actionPar);
-        actionPar.twist = "0. 0. 0. 0.";
+  MSolver::RBPrecCG::Par solverPar;
-        application.createModule<MAction::DWF>("DWF_s", actionPar);
+  solverPar.action       = "DWF_s";
-
+  solverPar.residual     = 1.0e-2;
-
+  solverPar.maxIteration = 10000;
-        MSolver::RBPrecCG::Par solverPar;
+  application.createModule<MSolver::RBPrecCG>("CG_s", solverPar);
-        solverPar.action       = "DWF_s";
-        solverPar.residual     = 1.0e-7;
-        solverPar.maxIteration = 10000;
-        application.createModule<MSolver::RBPrecCG>("CG_s", solverPar);
 }
+
 /////////////////////////////////////////////////////////////
 // Test creation of laplacian eigenvectors
 /////////////////////////////////////////////////////////////
@@ -138,7 +135,6 @@ void test_Perambulators(Application &application)
   PerambPar.Distil.TI=8;
   PerambPar.nvec=5;
   PerambPar.Distil.Ns=4;
-  PerambPar.Distil.Nt=8;
   PerambPar.Distil.Nt_inv=1;
   //PerambPar.Solver.mass=0.005;
   //PerambPar.Solver.M5=1.8;
@@ -166,7 +162,6 @@ void test_MultiPerambulators(Application &application)
   PerambPar.Distil.TI=8;
   PerambPar.nvec=5;
   PerambPar.Distil.Ns=4;
-  PerambPar.Distil.Nt=8;
   PerambPar.Distil.Nt_inv=1;
   application.createModule<MDistil::Peramb>("Peramb5",PerambPar);
   MDistil::PerambFromSolve::Par SolvePar;
@@ -181,7 +176,6 @@ void test_MultiPerambulators(Application &application)
   SolvePar.nvec_reduced=2;
   SolvePar.LI_reduced=2;
   SolvePar.Distil.Ns=4;
-  SolvePar.Distil.Nt=8;
   SolvePar.Distil.Nt_inv=1;
   application.createModule<MDistil::PerambFromSolve>("Peramb2",SolvePar);
   SolvePar.PerambFileName="Peramb3";
@@ -199,7 +193,6 @@ void test_MultiPerambulators(Application &application)
   DistilVecPar.TI=8;
   DistilVecPar.nvec=2;
   DistilVecPar.Ns=4;
-  DistilVecPar.Nt=8;
   DistilVecPar.Nt_inv=1;
   application.createModule<MDistil::DistilVectors>("DistilVecs2",DistilVecPar);
   DistilVecPar.perambulator="Peramb3";
@@ -215,7 +208,7 @@ void test_MultiPerambulators(Application &application)
   A2AMesonFieldPar.left="DistilVecs2_rho";
   A2AMesonFieldPar.right="DistilVecs2_rho";
   A2AMesonFieldPar.output="MesonSinksRho2";
-  A2AMesonFieldPar.gammas="all";
+  A2AMesonFieldPar.gammas="Identity";
   A2AMesonFieldPar.mom={"0 0 0"};
   A2AMesonFieldPar.cacheBlock=2;
   A2AMesonFieldPar.block=4;
@@ -247,11 +240,13 @@ void test_Noises(Application &application) {
   MDistil::NoisesPar NoisePar;
   NoisePar.UniqueIdentifier = "full_dilution";
   NoisePar.nvec = 5;
+  NoisePar.Distil.TI = 8;
+  NoisePar.Distil.SI = 4;
   NoisePar.Distil.LI = 5;
   NoisePar.Distil.nnoise = 1;
   NoisePar.Distil.Ns = 4;
-  NoisePar.Distil.Nt = 8;
+  NoisePar.Distil.Nt_inv = 1;
-  NoisePar.Distil.Nt = 1;
+  NoisePar.Distil.tsrc = 0;
   application.createModule<MDistil::Noises>("Peramb_noise",NoisePar);
 }
 
@@ -269,13 +264,11 @@ void test_DistilVectors(Application &application)
   DistilVecPar.perambulator="Peramb";
   DistilVecPar.lapevec="LapEvec";
   DistilVecPar.tsrc = 0;
-  //DistilVecPar.nnoise = 1;
   DistilVecPar.LI=5;
   DistilVecPar.SI=4;
   DistilVecPar.TI=8;
   DistilVecPar.nvec=5;
   DistilVecPar.Ns=4;
-  DistilVecPar.Nt=8;
   DistilVecPar.Nt_inv=1;
   application.createModule<MDistil::DistilVectors>("DistilVecs",DistilVecPar);
 }
@@ -294,7 +287,6 @@ void test_PerambulatorsS(Application &application)
   PerambPar.Distil.TI=8;
   PerambPar.nvec=5;
   PerambPar.Distil.Ns=4;
-  PerambPar.Distil.Nt=8;
   PerambPar.Distil.Nt_inv=1;
   //PerambPar.Solver.mass=0.005; //strange mass???
   //PerambPar.Solver.M5=1.8;
@@ -322,7 +314,6 @@ void test_DistilVectorsS(Application &application)
   DistilVecPar.TI=32;
   DistilVecPar.nvec=5;
   DistilVecPar.Ns=4;
-  DistilVecPar.Nt=8;
   DistilVecPar.Nt_inv=1;
   application.createModule<MDistil::DistilVectors>("DistilVecsS",DistilVecPar);
 }
@@ -338,12 +329,40 @@ void test_MesonSink(Application &application)
   A2AMesonFieldPar.left="g5phi";
   A2AMesonFieldPar.right="Peramb_unsmeared_sink";
   A2AMesonFieldPar.output="DistilFields";
-  A2AMesonFieldPar.gammas="all";
+  A2AMesonFieldPar.gammas="Identity";
   A2AMesonFieldPar.mom={"0 0 0"};
   A2AMesonFieldPar.cacheBlock=2;
   A2AMesonFieldPar.block=4;
   application.createModule<MContraction::A2AMesonField>("DistilMesonSink",A2AMesonFieldPar);
 }
+/////////////////////////////////////////////////////////////
+// MesonFields
+/////////////////////////////////////////////////////////////
+
+void test_MesonField(Application &application, const char * pszFileSuffix,
+                     const char * pszObjectLeft = nullptr, const char * pszObjectRight = nullptr )
+{
+  // DistilVectors parameters
+  if( pszObjectLeft == nullptr )
+    pszObjectLeft = pszFileSuffix;
+  if( pszObjectRight == nullptr )
+    pszObjectRight = pszObjectLeft;
+  MContraction::A2AMesonField::Par A2AMesonFieldPar;
+  A2AMesonFieldPar.left="DistilVecs";
+  A2AMesonFieldPar.right=A2AMesonFieldPar.left;
+  A2AMesonFieldPar.left.append( pszObjectLeft );
+  A2AMesonFieldPar.right.append( pszObjectRight );
+  A2AMesonFieldPar.output="MesonSinks";
+  A2AMesonFieldPar.output.append( pszFileSuffix );
+  A2AMesonFieldPar.gammas="Identity";
+  A2AMesonFieldPar.mom={"0 0 0"};
+  A2AMesonFieldPar.cacheBlock=2;
+  A2AMesonFieldPar.block=4;
+  std::string sObjectName{"DistilMesonField"};
+  sObjectName.append( pszFileSuffix );
+  application.createModule<MContraction::A2AMesonField>(sObjectName, A2AMesonFieldPar);
+}
+
 /////////////////////////////////////////////////////////////
 // g5*unsmeared
 /////////////////////////////////////////////////////////////
@@ -360,86 +379,11 @@ void test_g5_sinks(Application &application)
   g5_multiplyPar.Nt_inv=1;
   application.createModule<MDistil::g5_multiply>("g5phi",g5_multiplyPar);
 }
-#endif
 
-/////////////////////////////////////////////////////////////
-// MesonFields
-/////////////////////////////////////////////////////////////
-
-void test_MesonFieldSL(Application &application)
-{
-  // DistilVectors parameters
-  MContraction::A2AMesonField::Par A2AMesonFieldPar;
-  A2AMesonFieldPar.left="DistilVecsS_phi";
-  //A2AMesonFieldPar.right="DistilVecs_rho";
-  A2AMesonFieldPar.right="DistilVecs_phi";
-  A2AMesonFieldPar.output="DistilFieldsS";
-  A2AMesonFieldPar.gammas="all";
-  A2AMesonFieldPar.mom={"0 0 0"};
-  A2AMesonFieldPar.cacheBlock=2;
-  A2AMesonFieldPar.block=4;
-  application.createModule<MContraction::A2AMesonField>("DistilMesonFieldS",A2AMesonFieldPar);
-}
-/////////////////////////////////////////////////////////////
-// MesonFields - phiphi
-/////////////////////////////////////////////////////////////
-
-void test_MesonField(Application &application)
-{
-  // DistilVectors parameters
-  MContraction::A2AMesonField::Par A2AMesonFieldPar;
-  A2AMesonFieldPar.left="DistilVecs_phi";
-  //A2AMesonFieldPar.right="DistilVecs_rho";
-  A2AMesonFieldPar.right="DistilVecs_phi";
-  A2AMesonFieldPar.output="MesonSinksPhi";
-  A2AMesonFieldPar.gammas="all";
-  A2AMesonFieldPar.mom={"0 0 0"};
-  A2AMesonFieldPar.cacheBlock=2;
-  A2AMesonFieldPar.block=4;
-  application.createModule<MContraction::A2AMesonField>("DistilMesonField",A2AMesonFieldPar);
-}
-/////////////////////////////////////////////////////////////
-// MesonFields - rhorho
-/////////////////////////////////////////////////////////////
-
-void test_MesonFieldRho(Application &application)
-{
-  // DistilVectors parameters
-  MContraction::A2AMesonField::Par A2AMesonFieldPar;
-  A2AMesonFieldPar.left="DistilVecs_rho";
-  //A2AMesonFieldPar.right="DistilVecs_rho";
-  A2AMesonFieldPar.right="DistilVecs_rho";
-  A2AMesonFieldPar.output="MesonSinksRho";
-  A2AMesonFieldPar.gammas="all";
-  A2AMesonFieldPar.mom={"0 0 0"};
-  A2AMesonFieldPar.cacheBlock=2;
-  A2AMesonFieldPar.block=4;
-  application.createModule<MContraction::A2AMesonField>("DistilMesonFieldRho",A2AMesonFieldPar);
-}
-/////////////////////////////////////////////////////////////
-// MesonFields - rhorhoall
-/////////////////////////////////////////////////////////////
-
-void test_MesonFieldRhoAll(Application &application)
-{
-  // DistilVectors parameters
-  MContraction::A2AMesonField::Par A2AMesonFieldPar;
-  //A2AMesonFieldPar.left="DistilVecs_rho_all_tsrc";
-  //A2AMesonFieldPar.right="DistilVecs_rho_all_tsrc";
-  A2AMesonFieldPar.left="DistilVecs_rho";
-  A2AMesonFieldPar.right="DistilVecs_rho";
-  A2AMesonFieldPar.output="MesonSinksRhoAll";
-  A2AMesonFieldPar.gammas="all";
-  A2AMesonFieldPar.mom={"0 0 0"};
-  A2AMesonFieldPar.cacheBlock=2;
-  A2AMesonFieldPar.block=4;
-  application.createModule<MContraction::A2AMesonField>("DistilMesonFieldRhoAll",A2AMesonFieldPar);
-}
 /////////////////////////////////////////////////////////////
 // BaryonFields - phiphiphi - efficient
 /////////////////////////////////////////////////////////////
 
-#ifdef DISTIL_PRE_RELEASE
 void test_BaryonFieldPhi2(Application &application)
 {
   // DistilVectors parameters
@@ -577,7 +521,6 @@ void test_PerambulatorsSolve(Application &application)
   PerambFromSolvePar.Distil.TI=8;
   PerambFromSolvePar.nvec=5;
   PerambFromSolvePar.Distil.Ns=4;
-  PerambFromSolvePar.Distil.Nt=8;
   PerambFromSolvePar.Distil.Nt_inv=1;
   application.createModule<MDistil::PerambFromSolve>("PerambAslashS",PerambFromSolvePar);
 }
@@ -602,25 +545,6 @@ void test_DistilVectorsAslashSeq(Application &application)
   DistilSinkPar.Nt_inv=1;
   application.createModule<MDistil::DistilSink>("DistilVecsAslashSeq",DistilSinkPar);
 }
-/////////////////////////////////////////////////////////////
-// MesonFields - aslaaaash
-/////////////////////////////////////////////////////////////
-
-void test_MesonFieldAslashSeq(Application &application)
-{
-  // DistilVectors parameters
-  MContraction::A2AMesonField::Par A2AMesonFieldPar;
-  A2AMesonFieldPar.left="DistilVecsAslashSeq";
-  //A2AMesonFieldPar.right="DistilVecs_rho";
-  A2AMesonFieldPar.right="DistilVecsAslashSeq";
-  A2AMesonFieldPar.output="MesonSinksAslashSeq";
-  A2AMesonFieldPar.gammas="all";
-  A2AMesonFieldPar.mom={"0 0 0"};
-  A2AMesonFieldPar.cacheBlock=2;
-  A2AMesonFieldPar.block=4;
-  application.createModule<MContraction::A2AMesonField>("DistilMesonFieldAslashSeq",A2AMesonFieldPar);
-}
-
 bool bNumber( int &ri, const char * & pstr, bool bGobbleWhiteSpace = true )
 {
   if( bGobbleWhiteSpace )
@@ -774,7 +698,7 @@ bool DebugEigenTest()
   std::array<std::string,3> as={"Alpha", "Beta", "Gamma"};
   MyTensor x(as, 2,1,4);
   DebugShowTensor(x, "x");
-  x.WriteBinary(pszTestFileName);
+  x.write(pszTestFileName);
   // Test initialisation of an array of strings
   for( auto a : as )
     std::cout << a << std::endl;
@@ -787,7 +711,7 @@ bool DebugEigenTest()
   // Now see whether we can read a tensor back
   std::array<std::string,3> Names2={"Alpha", "Gamma", "Delta"};
   MyTensor y(Names2, 2,4,1);
-  y.ReadBinary(pszTestFileName);
+  y.read(pszTestFileName);
   DebugShowTensor(y, "y");
 
   // Now see whether we can read a tensor back from an hdf5 file
@@ -1037,15 +961,6 @@ int main(int argc, char *argv[])
 {
 #ifdef DEBUG
   // Debug only - test of Eigen::Tensor
-  std::cout << "sizeof(int) = " << sizeof(int)
-  << ", sizeof(long) = " << sizeof(long)
-  << ", sizeof(size_t) = " << sizeof(size_t)
-  << ", sizeof(std::size_t) = " << sizeof(std::size_t)
-  << ", sizeof(std::streamsize) = " << sizeof(std::streamsize)
-  << ", sizeof(Eigen::Index) = " << sizeof(Eigen::Index)
-  << ", sizeof(hsize_t) = " << sizeof(hsize_t)
-  << ", sizeof(unsigned long long) = " << sizeof(unsigned long long)
-  << std::endl;
   //if( DebugEigenTest() ) return 0;
   //if(DebugGridTensorTest()) return 0;
   //if(ConvertPeramb("PerambL_100_tsrc0.3000","PerambL_100_tsrc0.3000")) return 0;
@@ -1105,21 +1020,21 @@ int main(int argc, char *argv[])
       test_LapEvec( application );
       test_Perambulators( application );
       break;
-    case 3: // 3
+    default: // 3
       test_Global( application );
       test_SolverS( application );
       test_LapEvec( application );
       test_Perambulators( application );
       test_DistilVectors( application );
       break;
-    default: // 4
+    case 4:
       test_Global( application );
       test_SolverS( application );
       test_LapEvec( application );
       test_Perambulators( application );
       test_DistilVectors( application );
-      test_MesonField( application );
+      test_MesonField( application, "Phi", "_phi" );
-      test_MesonFieldRhoAll( application );
+      test_MesonField( application, "Rho", "_rho" );
       break;
     case 5: // 3
       test_Global( application );
@@ -1129,7 +1044,8 @@ int main(int argc, char *argv[])
       test_DistilVectors( application );
       test_PerambulatorsS( application );
       test_DistilVectorsS( application );
-      test_MesonFieldSL( application );
+      test_MesonField( application, "SPhi", "S_phi" );
+      test_MesonField( application, "SRho", "S_rho" );
       break;
 #ifdef DISTIL_PRE_RELEASE
     case 6: // 3
@@ -1156,7 +1072,8 @@ int main(int argc, char *argv[])
       test_LapEvec( application );
       test_Perambulators( application );
       test_DistilVectors( application );
-      test_MesonField( application );
+      test_MesonField( application, "Phi", "_phi" );
+      test_MesonField( application, "Rho", "_rho" );
       break;
 #ifdef DISTIL_PRE_RELEASE
     case 9: // 3
@@ -1192,7 +1109,7 @@ int main(int argc, char *argv[])
       test_AslashSeq( application );
       test_PerambulatorsSolve( application );
       test_DistilVectorsAslashSeq( application );
-      test_MesonFieldAslashSeq( application );
+      test_MesonField( application, "AslashSeq" );
       break;
     case 13:
       test_Global( application );
@@ -1201,10 +1118,9 @@ int main(int argc, char *argv[])
       test_MultiPerambulators( application );
       break;
   }
-  LOG(Message) << "====== XML creation for test " << iTestNum << " complete ======" << std::endl;
-
   // execution
   application.saveParameterFile("test_distil.xml");
+  LOG(Warning) << "These parameters are designed to run on a laptop usid --grid 4.4.4.8" << std::endl;
   application.run();
   
   // epilogue