Cleanup in progress

Hadrons/Distil.hpp

@@ -38,39 +38,6 @@
 #include <Hadrons/A2AVectors.hpp>
 #include <Hadrons/DilutedNoise.hpp>
 
-/******************************************************************************
- 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
- ******************************************************************************/
-
-#if defined(__linux__)
-#  include <endian.h>
-#elif defined(__FreeBSD__) || defined(__NetBSD__)
-#  include <sys/endian.h>
-#elif defined(__OpenBSD__)
-#  include <sys/types.h>
-#  define be16toh(x) betoh16(x)
-#  define be32toh(x) betoh32(x)
-#  define be64toh(x) betoh64(x)
-#elif defined(__APPLE__)
-#include <libkern/OSByteOrder.h>
-#define htobe16(x) OSSwapHostToBigInt16(x)
-#define htole16(x) OSSwapHostToLittleInt16(x)
-#define be16toh(x) OSSwapBigToHostInt16(x)
-#define le16toh(x) OSSwapLittleToHostInt16(x)
-
-#define htobe32(x) OSSwapHostToBigInt32(x)
-#define htole32(x) OSSwapHostToLittleInt32(x)
-#define be32toh(x) OSSwapBigToHostInt32(x)
-#define le32toh(x) OSSwapLittleToHostInt32(x)
-
-#define htobe64(x) OSSwapHostToBigInt64(x)
-#define htole64(x) OSSwapHostToLittleInt64(x)
-#define be64toh(x) OSSwapBigToHostInt64(x)
-#define le64toh(x) OSSwapLittleToHostInt64(x)
-#endif
-
 /******************************************************************************
  This potentially belongs in CartesianCommunicator
  ******************************************************************************/
@@ -139,7 +106,7 @@ inline void SliceShare( GridBase * gridLowDim, GridBase * gridHighDim, void * Bu
  *************************************************************************************/
 
 template<typename Field, typename GaugeField=LatticeGaugeField>
-class LinOpPeardonNabla : public LinearOperatorBase<Field>, public LinearFunction<Field> {
+class Laplacian3D : public LinearOperatorBase<Field>, public LinearFunction<Field> {
   typedef typename GaugeField::vector_type vCoeff_t;
 protected: // I don't really mind if _gf is messed with ... so make this public?
   //GaugeField & _gf;
@@ -147,7 +114,7 @@ protected: // I don't really mind if _gf is messed with ... so make this public?
   std::vector<Lattice<iColourMatrix<vCoeff_t> > > U;
 public:
   // Construct this operator given a gauge field and the number of dimensions it should act on
-  LinOpPeardonNabla( GaugeField& gf, int dimSpatial = Tdir ) : /*_gf(gf),*/ nd{dimSpatial} {
+  Laplacian3D( GaugeField& gf, int dimSpatial = Tdir ) : /*_gf(gf),*/ nd{dimSpatial} {
     assert(dimSpatial>=1);
     for( int mu = 0 ; mu < nd ; mu++ )
       U.push_back(PeekIndex<LorentzIndex>(gf,mu));
@@ -178,12 +145,12 @@ public:
 };
 
 template<typename Field>
-class LinOpPeardonNablaHerm : public LinearFunction<Field> {
+class Laplacian3DHerm : public LinearFunction<Field> {
 public:
   OperatorFunction<Field>   & _poly;
   LinearOperatorBase<Field> &_Linop;
   
-  LinOpPeardonNablaHerm(OperatorFunction<Field> & poly,LinearOperatorBase<Field>& linop)
+  Laplacian3DHerm(OperatorFunction<Field> & poly,LinearOperatorBase<Field>& linop)
   : _poly{poly}, _Linop{linop} {}
   
   void operator()(const Field& in, Field& out) {
@@ -244,12 +211,6 @@ const bool full_tdil{ TI == Nt }; \
 const bool exact_distillation{ full_tdil && LI == nvec }; \
 const int Nt_inv{ full_tdil ? 1 : TI }
 
-class BFieldIO: Serializable{
-public:
-  using BaryonTensorSet = Eigen::Tensor<ComplexD, 6>;
-  GRID_SERIALIZABLE_CLASS_MEMBERS(BFieldIO, BaryonTensorSet, BField );
-};
-
 /******************************************************************************
  Default for distillation file operations. For now only used by NamedTensor
 ******************************************************************************/
@@ -268,9 +229,6 @@ static const char * FileExtension = ".dat";
  NamedTensor object
  This is an Eigen::Tensor of type Scalar_ and rank NumIndices_ (row-major order)
  They can be persisted to disk
- Scalar_ objects are assumed to be composite objects of size Endian_Scalar_Size.
- (Disable big-endian by setting Endian_Scalar_Size=1).
- NB: Endian_Scalar_Size will disappear when ReadBinary & WriteBinary retired
  IndexNames contains one name for each index, and IndexNames are validated on load.
  WHAT TO SAVE / VALIDATE ON LOAD (Override to warn instead of assert on load)
  Ensemble string
@@ -280,20 +238,18 @@ static const char * FileExtension = ".dat";
 
  ******************************************************************************/
 
-template<typename Scalar_, int NumIndices_, uint16_t Endian_Scalar_Size_ = sizeof(Scalar_)>
+template<typename Scalar_, int NumIndices_>
 class NamedTensor : Serializable
 {
 public:
   using Scalar = Scalar_;
   static constexpr int NumIndices = NumIndices_;
-  static constexpr uint16_t Endian_Scalar_Size = Endian_Scalar_Size_;
   using ET = Eigen::Tensor<Scalar_, NumIndices_, Eigen::RowMajor>;
   using Index = typename ET::Index;
   GRID_SERIALIZABLE_CLASS_MEMBERS(NamedTensor
                                   , ET, tensor
                                   , std::vector<std::string>, IndexNames
                                   );
-public:
   // Named tensors are intended to be a superset of Eigen tensor
   inline operator ET&() { return tensor; }
   template<typename... IndexTypes>
@@ -351,9 +307,6 @@ public:
   // Read/Write in default format, i.e. HDF5 if present, else binary
   inline void read (const char * filename, const char * pszTag = nullptr);
   inline void write(const char * filename, const char * pszTag = nullptr) const;
-  // Original I/O implementation. This will be removed when we're sure it's no longer needed
-  EIGEN_DEPRECATED inline void ReadBinary (const std::string filename); // To be removed
-  EIGEN_DEPRECATED inline void WriteBinary(const std::string filename); // To be removed
 
   // Case insensitive compare of two strings
   // Pesumably this exists already? Where should this go?
@@ -375,218 +328,31 @@ public:
 
 // Is this a named tensor
 template<typename T, typename V = void> struct is_named_tensor : public std::false_type {};
-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 {};
+template<typename Scalar_, int NumIndices_> struct is_named_tensor<NamedTensor<Scalar_, NumIndices_>> : 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>::value>::type> : public std::true_type {};
 
 /******************************************************************************
  PerambTensor 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 PerambTensor = NamedTensor<SpinVector, 6, sizeof(Real)>;
+using PerambTensor = NamedTensor<SpinVector, 6>;
 static const std::array<std::string, 6> PerambIndexNames{"nT", "nVec", "LI", "nNoise", "nT_inv", "SI"};
 
-/******************************************************************************
- Save NamedTensor binary format (NB: On-disk format is Big Endian)
- Assumes the Scalar_ objects are contiguous (no padding)
- ******************************************************************************/
-
-template<typename Scalar_, int NumIndices_, uint16_t Endian_Scalar_Size>
-void NamedTensor<Scalar_, NumIndices_, Endian_Scalar_Size>::WriteBinary(const std::string filename) {
-  LOG(Message) << "Writing NamedTensor to \"" << filename << "\"" << std::endl;
-  std::ofstream w(filename, std::ios::binary);
-  // Enforce assumption that the scalar is composed of fundamental elements of size Endian_Scalar_Size
-  assert((Endian_Scalar_Size == 1 || Endian_Scalar_Size == 2 || Endian_Scalar_Size == 4 || Endian_Scalar_Size == 8 )
-         && "NamedTensor error: Endian_Scalar_Size should be 1, 2, 4 or 8");
-  assert((sizeof(Scalar_) % Endian_Scalar_Size) == 0 && "NamedTensor error: Scalar_ is not composed of Endian_Scalar_Size" );
-  // Size of the data (in bytes)
-  const uint32_t Scalar_Size{sizeof(Scalar_)};
-  const auto NumElements = tensor.size();
-  const std::streamsize TotalDataSize{static_cast<std::streamsize>(NumElements * Scalar_Size)};
-  uint64_t u64 = htobe64(static_cast<uint64_t>(TotalDataSize));
-  w.write(reinterpret_cast<const char *>(&u64), sizeof(u64));
-  // Size of a Scalar_
-  uint32_t u32{htobe32(Scalar_Size)};
-  w.write(reinterpret_cast<const char *>(&u32), sizeof(u32));
-  // Endian_Scalar_Size
-  uint16_t u16{htobe16(Endian_Scalar_Size)};
-  w.write(reinterpret_cast<const char *>(&u16), sizeof(u16));
-  // number of dimensions which aren't 1
-  u16 = static_cast<uint16_t>(this->NumIndices);
-  for( auto dim : tensor.dimensions() )
-    if( dim == 1 )
-      u16--;
-  u16 = htobe16( u16 );
-  w.write(reinterpret_cast<const char *>(&u16), sizeof(u16));
-  // dimensions together with names
-  int d = 0;
-  for( auto dim : tensor.dimensions() ) {
-    if( dim != 1 ) {
-      // size of this dimension
-      u16 = htobe16( static_cast<uint16_t>( dim ) );
-      w.write(reinterpret_cast<const char *>(&u16), sizeof(u16));
-      // length of this dimension name
-      u16 = htobe16( static_cast<uint16_t>( IndexNames[d].size() ) );
-      w.write(reinterpret_cast<const char *>(&u16), sizeof(u16));
-      // dimension name
-      w.write(IndexNames[d].c_str(), IndexNames[d].size());
-    }
-    d++;
-  }
-  // Actual data
-  char * const pStart{reinterpret_cast<char *>(tensor.data())};
-  // Swap to network byte order in place (alternative is to copy memory - still slow)
-  void * const pEnd{pStart + TotalDataSize};
-  if(Endian_Scalar_Size == 8)
-    for(uint64_t * p = reinterpret_cast<uint64_t *>(pStart) ; p < pEnd ; p++ )
-      * p = htobe64( * p );
-  else if(Endian_Scalar_Size == 4)
-    for(uint32_t * p = reinterpret_cast<uint32_t *>(pStart) ; p < pEnd ; p++ )
-      * p = htobe32( * p );
-  else if(Endian_Scalar_Size == 2)
-    for(uint16_t * p = reinterpret_cast<uint16_t *>(pStart) ; p < pEnd ; p++ )
-      * p = htobe16( * p );
-  w.write(pStart, TotalDataSize);
-  // Swap back from network byte order
-  if(Endian_Scalar_Size == 8)
-    for(uint64_t * p = reinterpret_cast<uint64_t *>(pStart) ; p < pEnd ; p++ )
-      * p = be64toh( * p );
-  else if(Endian_Scalar_Size == 4)
-    for(uint32_t * p = reinterpret_cast<uint32_t *>(pStart) ; p < pEnd ; p++ )
-      * p = be32toh( * p );
-  else if(Endian_Scalar_Size == 2)
-    for(uint16_t * p = reinterpret_cast<uint16_t *>(pStart) ; p < pEnd ; p++ )
-      * p = be16toh( * p );
-  // checksum
-#ifdef USE_IPP
-  u32 = htobe32(GridChecksum::crc32c(tensor.data(), TotalDataSize));
-#else
-  u32 = htobe32(GridChecksum::crc32(tensor.data(), TotalDataSize));
-#endif
-  w.write(reinterpret_cast<const char *>(&u32), sizeof(u32));
-}
-
-/******************************************************************************
- Load NamedTensor binary format (NB: On-disk format is Big Endian)
- Assumes the Scalar_ objects are contiguous (no padding)
- ******************************************************************************/
-
-template<typename Scalar_, int NumIndices_, uint16_t Endian_Scalar_Size>
-void NamedTensor<Scalar_, NumIndices_, Endian_Scalar_Size>::ReadBinary(const std::string filename) {
-  LOG(Message) << "Reading NamedTensor from \"" << filename << "\"" << std::endl;
-  std::ifstream r(filename, std::ios::binary);
-  // Enforce assumption that the scalar is composed of fundamental elements of size Endian_Scalar_Size
-  assert((Endian_Scalar_Size == 1 || Endian_Scalar_Size == 2 || Endian_Scalar_Size == 4 || Endian_Scalar_Size == 8 )
-         && "NamedTensor error: Endian_Scalar_Size should be 1, 2, 4 or 8");
-  assert((sizeof(Scalar_) % Endian_Scalar_Size) == 0 && "NamedTensor error: Scalar_ is not composed of Endian_Scalar_Size" );
-  // Size of the data in bytes
-  const uint32_t Scalar_Size{sizeof(Scalar_)};
-  Index NumElements{tensor.size()};
-  std::streamsize TotalDataSize{static_cast<std::streamsize>(NumElements * Scalar_Size)};
-  uint64_t u64;
-  r.read(reinterpret_cast<char *>(&u64), sizeof(u64));
-  assert( TotalDataSize == 0 || TotalDataSize == be64toh( u64 ) && "NamedTensor error: Size of the data in bytes" );
-  // Size of a Scalar_
-  uint32_t u32;
-  r.read(reinterpret_cast<char *>(&u32), sizeof(u32));
-  assert( Scalar_Size == be32toh( u32 ) && "NamedTensor error: sizeof(Scalar_)");
-  // Endian_Scalar_Size
-  uint16_t u16;
-  r.read(reinterpret_cast<char *>(&u16), sizeof(u16));
-  assert( Endian_Scalar_Size == be16toh( u16 ) && "NamedTensor error: Scalar_Unit_size");
-  // number of dimensions which aren't 1
-  uint16_t NumFileDimensions;
-  r.read(reinterpret_cast<char *>(&NumFileDimensions), sizeof(NumFileDimensions));
-  NumFileDimensions = be16toh( NumFileDimensions );
-  /*for( auto dim : tensor.dimensions() )
-    if( dim == 1 )
-      u16++;*/
-  assert( ( TotalDataSize == 0 && this->NumIndices >= NumFileDimensions || this->NumIndices == NumFileDimensions )
-         && "NamedTensor error: number of dimensions which aren't 1" );
-  if( TotalDataSize == 0 ) {
-    // Read each dimension, using names to skip past dimensions == 1
-    std::array<Index,NumIndices_> NewDimensions;
-    for( Index &i : NewDimensions ) i = 1;
-    int d = 0;
-    for( int FileDimension = 0; FileDimension < NumFileDimensions; FileDimension++ ) {
-      // read dimension
-      uint16_t thisDim;
-      r.read(reinterpret_cast<char *>(&thisDim), sizeof(thisDim));
-      // read dimension name
-      r.read(reinterpret_cast<char *>(&u16), sizeof(u16));
-      size_t l = be16toh( u16 );
-      std::string s( l, '?' );
-      r.read(&s[0], l);
-      // skip forward to matching name
-      while( IndexNames[d].size() > 0 && !CompareCaseInsensitive( s, IndexNames[d] ) )
-        assert(++d < NumIndices && "NamedTensor error: dimension name" );
-      if( IndexNames[d].size() == 0 )
-        IndexNames[d] = s;
-      NewDimensions[d++] = be16toh( thisDim );
-    }
-    tensor.resize(NewDimensions);
-    NumElements = 1;
-    for( Index i : NewDimensions ) NumElements *= i;
-    TotalDataSize = NumElements * Scalar_Size;
-  } else {
-    // dimensions together with names
-    const auto & TensorDims = tensor.dimensions();
-    for( int d = 0; d < NumIndices_; d++ ) {
-      // size of dimension
-      r.read(reinterpret_cast<char *>(&u16), sizeof(u16));
-      u16 = be16toh( u16 );
-      assert( TensorDims[d] == u16 && "size of dimension" );
-      // length of dimension name
-      r.read(reinterpret_cast<char *>(&u16), sizeof(u16));
-      size_t l = be16toh( u16 );
-      assert( l == IndexNames[d].size() && "NamedTensor error: length of dimension name" );
-      // dimension name
-      std::string s( l, '?' );
-      r.read(&s[0], l);
-      assert( s == IndexNames[d] && "NamedTensor error: dimension name" );
-    }
-  }
-  // Actual data
-  char * const pStart{reinterpret_cast<char *>(tensor.data())};
-  void * const pEnd{pStart + TotalDataSize};
-  r.read(pStart,TotalDataSize);
-  // Swap back from network byte order
-  if(Endian_Scalar_Size == 8)
-    for(uint64_t * p = reinterpret_cast<uint64_t *>(pStart) ; p < pEnd ; p++ )
-      * p = be64toh( * p );
-  else if(Endian_Scalar_Size == 4)
-    for(uint32_t * p = reinterpret_cast<uint32_t *>(pStart) ; p < pEnd ; p++ )
-      * p = be32toh( * p );
-  else if(Endian_Scalar_Size == 2)
-    for(uint16_t * p = reinterpret_cast<uint16_t *>(pStart) ; p < pEnd ; p++ )
-      * p = be16toh( * p );
-  // checksum
-  r.read(reinterpret_cast<char *>(&u32), sizeof(u32));
-  u32 = be32toh( u32 );
-#ifdef USE_IPP
-  u32 -= GridChecksum::crc32c(tensor.data(), TotalDataSize);
-#else
-  u32 -= GridChecksum::crc32(tensor.data(), TotalDataSize);
-#endif
-  assert( u32 == 0 && "NamedTensor error: PerambTensor checksum invalid");
-}
-
 /******************************************************************************
  Write NamedTensor
  ******************************************************************************/
 
-template<typename Scalar_, int NumIndices_, uint16_t Endian_Scalar_Size>
+template<typename Scalar_, int NumIndices_>
 template<typename Writer>
-void NamedTensor<Scalar_, NumIndices_, Endian_Scalar_Size>::write(Writer &w, const char * pszTag)const{
+void NamedTensor<Scalar_, NumIndices_>::write(Writer &w, const char * pszTag)const{
   if( pszTag == nullptr )
     pszTag = "NamedTensor";
   LOG(Message) << "Writing NamedTensor to tag  " << pszTag << std::endl;
   write(w, pszTag, *this);
 }
 
-template<typename Scalar_, int NumIndices_, uint16_t Endian_Scalar_Size>
-void NamedTensor<Scalar_, NumIndices_, Endian_Scalar_Size>::write(const char * filename, const char * pszTag)const{
+template<typename Scalar_, int NumIndices_>
+void NamedTensor<Scalar_, NumIndices_>::write(const char * filename, const char * pszTag)const{
   std::string sFileName{filename};
   sFileName.append( MDistil::FileExtension );
   LOG(Message) << "Writing NamedTensor to file " << sFileName << std::endl;
@@ -598,8 +364,8 @@ void NamedTensor<Scalar_, NumIndices_, Endian_Scalar_Size>::write(const char * f
  Validate named tensor index names
  ******************************************************************************/
 
-template<typename Scalar_, int NumIndices_, uint16_t Endian_Scalar_Size>
-bool NamedTensor<Scalar_, NumIndices_, Endian_Scalar_Size>::ValidateIndexNames( int iNumNames, const std::string * MatchNames ) const {
+template<typename Scalar_, int NumIndices_>
+bool NamedTensor<Scalar_, NumIndices_>::ValidateIndexNames( int iNumNames, const std::string * MatchNames ) const {
   bool bSame{ iNumNames == NumIndices_ && IndexNames.size() == NumIndices_ };
   for( int i = 0; bSame && i < NumIndices_; i++ )
     bSame = CompareCaseInsensitive( MatchNames[i], IndexNames[i] );
@@ -610,9 +376,9 @@ bool NamedTensor<Scalar_, NumIndices_, Endian_Scalar_Size>::ValidateIndexNames(
  Read NamedTensor
  ******************************************************************************/
 
-template<typename Scalar_, int NumIndices_, uint16_t Endian_Scalar_Size>
+template<typename Scalar_, int NumIndices_>
 template<typename Reader>
-void NamedTensor<Scalar_, NumIndices_, Endian_Scalar_Size>::read(Reader &r, const char * pszTag) {
+void NamedTensor<Scalar_, NumIndices_>::read(Reader &r, const char * pszTag) {
   if( pszTag == nullptr )
     pszTag = "NamedTensor";
   // Grab index names and dimensions
@@ -626,8 +392,8 @@ void NamedTensor<Scalar_, NumIndices_, Endian_Scalar_Size>::read(Reader &r, cons
   assert( ValidateIndexNames( OldIndexNames.size(), &OldIndexNames[0] ) && "NamedTensor::load dimension name" );
 }
 
-template<typename Scalar_, int NumIndices_, uint16_t Endian_Scalar_Size>
-void NamedTensor<Scalar_, NumIndices_, Endian_Scalar_Size>::read(const char * filename, const char * pszTag) {
+template<typename Scalar_, int NumIndices_>
+void NamedTensor<Scalar_, NumIndices_>::read(const char * filename, const char * pszTag) {
   std::string sFileName{filename};
   sFileName.append( MDistil::FileExtension );
   LOG(Message) << "Reading NamedTensor from file " << sFileName << std::endl;
@@ -664,19 +430,14 @@ inline void RotateEigen(std::vector<LatticeColourVector> & evec)
   Coordinate siteFirst(grid->Nd(),0);
   peekSite(cv0, evec[0], siteFirst);
   Grid::Complex cplx0 = cv0()()(0);
-#ifdef GRID_NVCC
   if( cplx0.imag() == 0 )
-#else
-  if( std::imag(cplx0) == 0 )
-#endif
     std::cout << GridLogMessage << "RotateEigen() : Site 0 : " << cplx0 << " => already meets phase convention" << std::endl;
   else {
+    const Real cplx0_mag = Grid::abs(cplx0);
 #ifdef GRID_NVCC
-    const Real cplx0_mag = thrust::abs(cplx0);
     const Grid::Complex phase = thrust::conj(cplx0 / cplx0_mag);
     const Real argphase = thrust::arg(phase);
 #else
-    const Real cplx0_mag = std::abs(cplx0);
     const Grid::Complex phase = std::conj(cplx0 / cplx0_mag);
     const Real argphase = std::arg(phase);
 #endif

Hadrons/Modules/MDistil/LapEvec.hpp

@@ -219,7 +219,7 @@ void TLapEvec<GImpl>::execute(void)
   }
 
   ////////////////////////////////////////////////////////////////////////
-  // Invert Peardon Nabla operator separately on each time-slice
+  // Invert nabla operator separately on each time-slice
   ////////////////////////////////////////////////////////////////////////
   
   auto & eig4d = envGet(LapEvecs, getName() );
@@ -237,7 +237,7 @@ void TLapEvec<GImpl>::execute(void)
     
     // Construct smearing operator
     ExtractSliceLocal(UmuNoTime,Umu_smear,0,t,Tdir); // switch to 3d/4d objects
-    LinOpPeardonNabla<LatticeColourVector> PeardonNabla(UmuNoTime);
+    Laplacian3D<LatticeColourVector> Nabla(UmuNoTime);
     LOG(Debug) << "Chebyshev preconditioning to order " << ChebPar.PolyOrder
       << " with parameters (alpha,beta) = (" << ChebPar.alpha << "," << ChebPar.beta << ")" << std::endl;
     Chebyshev<LatticeColourVector> Cheb(ChebPar.alpha,ChebPar.beta,ChebPar.PolyOrder);
@@ -248,9 +248,9 @@ void TLapEvec<GImpl>::execute(void)
     nn = Grid::sqrt(nn);
     src = src * (1.0/nn);
 
-    LinOpPeardonNablaHerm<LatticeColourVector> PeardonNablaCheby(Cheb,PeardonNabla);
+    Laplacian3DHerm<LatticeColourVector> NablaCheby(Cheb,Nabla);
     ImplicitlyRestartedLanczos<LatticeColourVector>
-      IRL(PeardonNablaCheby,PeardonNabla,LPar.Nvec,LPar.Nk,LPar.Nk+LPar.Np,LPar.resid,LPar.MaxIt);
+      IRL(NablaCheby,Nabla,LPar.Nvec,LPar.Nk,LPar.Nk+LPar.Np,LPar.resid,LPar.MaxIt);
     int Nconv = 0;
     IRL.calc(eig[t].eval,eig[t].evec,src,Nconv);
     if( Nconv < LPar.Nvec ) {

Hadrons/Modules/MDistil/Noises.hpp

@@ -134,6 +134,7 @@ void TNoises<FImpl>::execute(void)
   }
   UniqueIdentifier.append( std::to_string( vm().getTrajectory() ) ); 
 
+  // We use our own seeds so we can specify different noises per quark
   GridSerialRNG sRNG;
   sRNG.SeedUniqueString(UniqueIdentifier);
   Real rn;

Hadrons/Modules/MDistil/Perambulator.hpp

@@ -30,7 +30,6 @@
 #ifndef Hadrons_MDistil_Perambulator_hpp_
 #define Hadrons_MDistil_Perambulator_hpp_
 
-// These are members of Distillation
 #include <Hadrons/Distil.hpp>
 
 BEGIN_HADRONS_NAMESPACE
@@ -154,7 +153,8 @@ void TPerambulator<FImpl>::setup(void)
 template <typename FImpl>
 void TPerambulator<FImpl>::Cleanup(void)
 {
-  if( grid3d != nullptr ) {
+    if( grid3d != nullptr )
+    {
         delete grid3d;
         grid3d = nullptr;
     }
@@ -166,23 +166,15 @@ template <typename FImpl>
 void TPerambulator<FImpl>::execute(void)
 {
     DISTIL_PARAMETERS_DEFINE( false );
-
     auto &solver=envGet(Solver, par().solver);
     auto &mat = solver.getFMat();
     envGetTmp(FermionField, v4dtmp);
     envGetTmp(FermionField, v5dtmp);
     envGetTmp(FermionField, v5dtmp_sol);
-
     auto &noise = envGet(NoiseTensor, sNoiseName);
     auto &perambulator = envGet(PerambTensor, getName());
     auto &epack = envGet(LapEvecs, sLapEvecName);
     auto &unsmeared_sink = envGet(std::vector<FermionField>, getName() + "_unsmeared_sink");
-
-
-    // Load perambulator if it exists on disk instead of creating it
-    // Not sure this is how we want it - rather specify an input flag 'read' 
-    // and assert that the file is there.
-
     envGetTmp(LatticeSpinColourVector, dist_source);
     envGetTmp(LatticeSpinColourVector, tmp2);
     envGetTmp(LatticeSpinColourVector, result);
@@ -191,26 +183,33 @@ void TPerambulator<FImpl>::execute(void)
     envGetTmp(LatticeColourVector, tmp3d_nospin);
     envGetTmp(LatticeColourVector, result_3d);
     envGetTmp(LatticeColourVector, evec3d);
-
     const int Ntlocal{grid4d->LocalDimensions()[3]};
     const int Ntfirst{grid4d->LocalStarts()[3]};
+    const std::string UnsmearedSinkFileName{ par().UnsmearedSinkFileName };
     
     {
-
         int t_inv;
-    for (int inoise = 0; inoise < nnoise; inoise++) {
-      for (int dk = 0; dk < LI; dk++) {
-        for (int dt = 0; dt < Nt_inv; dt++) {
-          for (int ds = 0; ds < SI; ds++) {
+        for (int inoise = 0; inoise < nnoise; inoise++)
+        {
+            for (int dk = 0; dk < LI; dk++)
+            {
+                for (int dt = 0; dt < Nt_inv; dt++)
+                {
+                    for (int ds = 0; ds < SI; ds++)
+                    {
                         LOG(Message) <<  "LapH source vector from noise " << inoise << " and dilution component (d_k,d_t,d_alpha) : (" << dk << ","<< dt << "," << ds << ")" << std::endl;
                         dist_source = 0;
                         tmp3d_nospin = 0;
                         evec3d = 0;
-            for (int it = dt; it < Nt; it += TI){
+                        for (int it = dt; it < Nt; it += TI)
+                        {
                             if (full_tdil) t_inv = tsrc; else t_inv = it;
-              if( t_inv >= Ntfirst && t_inv < Ntfirst + Ntlocal ) {
-                for (int ik = dk; ik < nvec; ik += LI){
-                  for (int is = ds; is < Ns; is += SI){ 
+                            if( t_inv >= Ntfirst && t_inv < Ntfirst + Ntlocal )
+                            {
+                                for (int ik = dk; ik < nvec; ik += LI)
+                                {
+                                    for (int is = ds; is < Ns; is += SI)
+                                    {
                                         ExtractSliceLocal(evec3d,epack.evec[ik],0,t_inv-Ntfirst,Tdir);
                                         tmp3d_nospin = evec3d * noise(inoise, t_inv, ik, is);
                                         tmp3d=0;
@@ -224,9 +223,12 @@ void TPerambulator<FImpl>::execute(void)
                         }
                         result=0;
                         v4dtmp = dist_source;
-	    if (Ls_ == 1){
+                        if (Ls_ == 1)
+                        {
                             solver(result, v4dtmp);
-	    } else {
+                        }
+                        else
+                        {
                             mat.ImportPhysicalFermionSource(v4dtmp, v5dtmp);
                             solver(v5dtmp_sol, v5dtmp);
                             mat.ExportPhysicalFermionSolution(v5dtmp_sol, v4dtmp);
@@ -234,11 +236,14 @@ void TPerambulator<FImpl>::execute(void)
                         }
                         if( !UnsmearedSinkFileName.empty() )
                             unsmeared_sink[inoise+nnoise*(dk+LI*(dt+Nt_inv*ds))] = result;
-            for (int is = 0; is < Ns; is++) {
+                        for (int is = 0; is < Ns; is++)
+                        {
                             result_nospin = peekSpin(result,is);
-              for (int t = Ntfirst; t < Ntfirst + Ntlocal; t++) {
+                            for (int t = Ntfirst; t < Ntfirst + Ntlocal; t++)
+                            {
                                 ExtractSliceLocal(result_3d,result_nospin,0,t-Ntfirst,Tdir);
-                for (int ivec = 0; ivec < nvec; ivec++) {
+                                for (int ivec = 0; ivec < nvec; ivec++)
+                                {
                                     ExtractSliceLocal(evec3d,epack.evec[ivec],0,t-Ntfirst,Tdir);
                                     pokeSpin(perambulator(t, ivec, dk, inoise,dt,ds),static_cast<Complex>(innerProduct(evec3d, result_3d)),is);
                                 }
@@ -252,7 +257,8 @@ void TPerambulator<FImpl>::execute(void)
     LOG(Message) <<  "perambulator done" << std::endl;
     perambulator.SliceShare( grid3d, grid4d );
     
-  if(grid4d->IsBoss()) {
+    if(grid4d->IsBoss())
+    {
         std::string sPerambName{par().PerambFileName};
         if( sPerambName.length() == 0 )
             sPerambName = getName();
@@ -261,8 +267,8 @@ void TPerambulator<FImpl>::execute(void)
         perambulator.write(sPerambName.c_str());
     }
     
-  const std::string UnsmearedSinkFileName{ par().UnsmearedSinkFileName };
-  if( !UnsmearedSinkFileName.empty() ) {
+    if( !UnsmearedSinkFileName.empty() )
+    {
         bool bMulti = ( Hadrons::MDistil::DistilParameters::ParameterDefault( par().UnsmearedSinkMultiFile, 1, false ) != 0 );
         LOG(Message) << "Writing unsmeared sink to " << UnsmearedSinkFileName << std::endl;
         A2AVectorsIo::write(UnsmearedSinkFileName, unsmeared_sink, bMulti, vm().getTrajectory());

tests/hadrons/Test_distil.cc

@@ -190,98 +190,22 @@ void test_Perambulators( Application &application, const char * pszSuffix = null
 // DistilVectors
 /////////////////////////////////////////////////////////////
 
-#define TEST_DISTIL_VECTORS_COMMON \
-std::string sModuleName{"DistilVecs"}; \
-if( pszSuffix ) \
-  sModuleName.append( pszSuffix ); \
-std::string sPerambName{"Peramb"}; \
-if( pszSuffix ) \
-  sPerambName.append( pszSuffix ); \
-MDistil::DistilVectors::Par DistilVecPar; \
-DistilVecPar.noise = sPerambName + "_noise"; \
-DistilVecPar.perambulator = sPerambName; \
-DistilVecPar.lapevec = "LapEvec"; \
-DistilVecPar.tsrc = 0; \
-if( pszNvec ) \
-  DistilVecPar.nvec = pszNvec
-
-#define TEST_DISTIL_VECTORS_COMMON_END \
-application.createModule<MDistil::DistilVectors>(sModuleName,DistilVecPar)
-
 void test_DistilVectors(Application &application, const char * pszSuffix = nullptr, const char * pszNvec = nullptr )
 {
-  TEST_DISTIL_VECTORS_COMMON;
-  TEST_DISTIL_VECTORS_COMMON_END;
-}
-
-void test_DistilVectorsSS(Application &application, const char * pszSink, const char * pszSource,
-                          const char * pszSuffix = nullptr, const char * pszNvec = nullptr )
-{
-  TEST_DISTIL_VECTORS_COMMON;
-  if( pszSink )
-    DistilVecPar.sink = pszSink;
-  if( pszSource )
-    DistilVecPar.source = pszSource;
-  TEST_DISTIL_VECTORS_COMMON_END;
-}
-
-/////////////////////////////////////////////////////////////
-// Multiple Perambulators
-/////////////////////////////////////////////////////////////
-
-void test_MultiPerambulators(Application &application)
-{
-  test_Perambulators( application, "5" );
-  MDistil::PerambFromSolve::Par SolvePar;
-  SolvePar.eigenPack="LapEvec";
-  SolvePar.PerambFileName="Peramb2";
-  SolvePar.solve = "Peramb5_unsmeared_sink";
-  SolvePar.Distil.nnoise = 1;
-  SolvePar.Distil.LI=5;
-  SolvePar.Distil.SI=4;
-  SolvePar.Distil.TI=8;
-  SolvePar.nvec=5;
-  SolvePar.nvec_reduced=2;
-  SolvePar.LI_reduced=2;
-  application.createModule<MDistil::PerambFromSolve>("Peramb2",SolvePar);
-  SolvePar.PerambFileName="Peramb3";
-  SolvePar.nvec_reduced=3;
-  SolvePar.LI_reduced=3;
-  application.createModule<MDistil::PerambFromSolve>("Peramb3",SolvePar);
-
-  test_DistilVectors( application, "2", "2" );
-  test_DistilVectors( application, "3", "3" );
-  test_DistilVectors( application, "5", "5" );
-
-  MContraction::A2AMesonField::Par A2AMesonFieldPar;
-  A2AMesonFieldPar.left="DistilVecs2_rho";
-  A2AMesonFieldPar.right="DistilVecs2_rho";
-  A2AMesonFieldPar.output="MesonSinksRho2";
-  A2AMesonFieldPar.gammas="Identity";
-  A2AMesonFieldPar.mom={"0 0 0"};
-  A2AMesonFieldPar.cacheBlock=2;
-  A2AMesonFieldPar.block=4;
-  application.createModule<MContraction::A2AMesonField>("DistilMesonFieldRho2",A2AMesonFieldPar);
-  A2AMesonFieldPar.left="DistilVecs2_phi";
-  A2AMesonFieldPar.right="DistilVecs2_phi";
-  A2AMesonFieldPar.output="MesonSinksPhi2";
-  application.createModule<MContraction::A2AMesonField>("DistilMesonFieldPhi2",A2AMesonFieldPar);
-  A2AMesonFieldPar.left="DistilVecs3_rho";
-  A2AMesonFieldPar.right="DistilVecs3_rho";
-  A2AMesonFieldPar.output="MesonSinksRho3";
-  application.createModule<MContraction::A2AMesonField>("DistilMesonFieldRho3",A2AMesonFieldPar);
-  A2AMesonFieldPar.left="DistilVecs3_phi";
-  A2AMesonFieldPar.right="DistilVecs3_phi";
-  A2AMesonFieldPar.output="MesonSinksPhi3";
-  application.createModule<MContraction::A2AMesonField>("DistilMesonFieldPhi3",A2AMesonFieldPar);
-  A2AMesonFieldPar.left="DistilVecs5_rho";
-  A2AMesonFieldPar.right="DistilVecs5_rho";
-  A2AMesonFieldPar.output="MesonSinksRho5";
-  application.createModule<MContraction::A2AMesonField>("DistilMesonFieldRho5",A2AMesonFieldPar);
-  A2AMesonFieldPar.left="DistilVecs5_phi";
-  A2AMesonFieldPar.right="DistilVecs5_phi";
-  A2AMesonFieldPar.output="MesonSinksPhi5";
-  application.createModule<MContraction::A2AMesonField>("DistilMesonFieldPhi5",A2AMesonFieldPar);
+  std::string sModuleName{"DistilVecs"};
+  if( pszSuffix )
+    sModuleName.append( pszSuffix );
+  std::string sPerambName{"Peramb"};
+  if( pszSuffix )
+    sPerambName.append( pszSuffix );
+  MDistil::DistilVectors::Par DistilVecPar;
+  DistilVecPar.noise = sPerambName + "_noise";
+  DistilVecPar.perambulator = sPerambName;
+  DistilVecPar.lapevec = "LapEvec";
+  DistilVecPar.tsrc = 0;
+  if( pszNvec )
+    DistilVecPar.nvec = pszNvec;
+  application.createModule<MDistil::DistilVectors>(sModuleName,DistilVecPar);
 }
 
 /////////////////////////////////////////////////////////////
@@ -330,163 +254,6 @@ void test_MesonField(Application &application, const char * pszFileSuffix,
   application.createModule<MContraction::A2AMesonField>(sObjectName, A2AMesonFieldPar);
 }
 
-/////////////////////////////////////////////////////////////
-// g5*unsmeared
-/////////////////////////////////////////////////////////////
-
-#ifdef DISTIL_PRE_RELEASE
-void test_g5_sinks(Application &application)
-{
-  // DistilVectors parameters
-  MDistil::g5_multiply::Par g5_multiplyPar;
-  g5_multiplyPar.input="Peramb_unsmeared_sink";
-  g5_multiplyPar.nnoise = 1;
-  g5_multiplyPar.LI=5;
-  g5_multiplyPar.Ns=4;
-  g5_multiplyPar.Nt_inv=1;
-  application.createModule<MDistil::g5_multiply>("g5phi",g5_multiplyPar);
-}
-
-/////////////////////////////////////////////////////////////
-// BaryonFields - phiphiphi - efficient
-/////////////////////////////////////////////////////////////
-
-void test_BaryonFieldPhi2(Application &application)
-{
-  // DistilVectors parameters
-  MDistil::BC2::Par BC2Par;
-  BC2Par.one="DistilVecs_phi";
-  BC2Par.two="DistilVecs_phi";
-  BC2Par.three="DistilVecs_phi";
-  BC2Par.output="BaryonFieldPhi2";
-  BC2Par.parity=1;
-  BC2Par.mom={"0 0 0"};
-  application.createModule<MDistil::BC2>("BaryonFieldPhi2",BC2Par);
-}
-/////////////////////////////////////////////////////////////
-// BaryonFields - rhorhorho - efficient
-/////////////////////////////////////////////////////////////
-
-void test_BaryonFieldRho2(Application &application)
-{
-  // DistilVectors parameters
-  MDistil::BC2::Par BC2Par;
-  BC2Par.one="DistilVecs_rho";
-  BC2Par.two="DistilVecs_rho";
-  BC2Par.three="DistilVecs_rho";
-  BC2Par.output="BaryonFieldRho2";
-  BC2Par.parity=1;
-  BC2Par.mom={"0 0 0"};
-  application.createModule<MDistil::BC2>("BaryonFieldRho2",BC2Par);
-}
-/////////////////////////////////////////////////////////////
-// BaryonFields - phiphiphi
-/////////////////////////////////////////////////////////////
-
-void test_BaryonFieldPhi(Application &application)
-{
-  // DistilVectors parameters
-  MDistil::BContraction::Par BContractionPar;
-  BContractionPar.one="DistilVecs_phi";
-  BContractionPar.two="DistilVecs_phi";
-  BContractionPar.three="DistilVecs_phi";
-  BContractionPar.output="BaryonFieldPhi";
-  BContractionPar.parity=1;
-  BContractionPar.mom={"0 0 0"};
-  application.createModule<MDistil::BContraction>("BaryonFieldPhi",BContractionPar);
-}
-/////////////////////////////////////////////////////////////
-// BaryonFields - rhorhorho
-/////////////////////////////////////////////////////////////
-
-void test_BaryonFieldRho(Application &application)
-{
-  // DistilVectors parameters
-  MDistil::BContraction::Par BContractionPar;
-  BContractionPar.one="DistilVecs_rho";
-  BContractionPar.two="DistilVecs_rho";
-  BContractionPar.three="DistilVecs_rho";
-  BContractionPar.output="BaryonFieldRho";
-  BContractionPar.parity=1;
-  BContractionPar.mom={"0 0 0"};
-  application.createModule<MDistil::BContraction>("BaryonFieldRho",BContractionPar);
-}
-/////////////////////////////////////////////////////////////
-// BaryonContraction
-/////////////////////////////////////////////////////////////
-
-void test_Baryon2pt(Application &application)
-{
-  // DistilVectors parameters
-  MDistil::Baryon2pt::Par Baryon2ptPar;
-  Baryon2ptPar.inputL="BaryonFieldPhi";
-  Baryon2ptPar.inputR="BaryonFieldRho";
-  Baryon2ptPar.quarksL="uud";
-  Baryon2ptPar.quarksR="uud";
-  Baryon2ptPar.output="C2_baryon";
-  application.createModule<MDistil::Baryon2pt>("C2_b",Baryon2ptPar);
-}
-#endif
-
-/////////////////////////////////////////////////////////////
-// emField
-/////////////////////////////////////////////////////////////
-void test_em(Application &application)
-{
-  MGauge::StochEm::Par StochEmPar;
-  StochEmPar.gauge=PhotonR::Gauge::feynman;
-  StochEmPar.zmScheme=PhotonR::ZmScheme::qedL;
-  application.createModule<MGauge::StochEm>("Em",StochEmPar);
-}
-
-/////////////////////////////////////////////////////////////
-// MesonA2ASlash
-/////////////////////////////////////////////////////////////
-
-void test_Aslash(Application &application)
-{
-  // DistilVectors parameters
-  MContraction::A2AAslashField::Par A2AAslashFieldPar;
-  A2AAslashFieldPar.left="g5phi";
-  //A2AAslashFieldPar.right="DistilVecs_phi";
-  A2AAslashFieldPar.right="Peramb_unsmeared_sink";
-  A2AAslashFieldPar.output="unsmeared_Aslash";
-  A2AAslashFieldPar.emField={"Em"};
-  A2AAslashFieldPar.cacheBlock=2;
-  A2AAslashFieldPar.block=4;
-  application.createModule<MContraction::A2AAslashField>("Aslash_field",A2AAslashFieldPar);
-}
-
-/////////////////////////////////////////////////////////////
-// MesonA2ASlashSequential
-/////////////////////////////////////////////////////////////
-
-void test_AslashSeq(Application &application)
-{
-  // DistilVectors parameters
-  MSolver::A2AAslashVectors::Par A2AAslashVectorsPar;
-  A2AAslashVectorsPar.vector="PerambS_unsmeared_sink";
-  A2AAslashVectorsPar.emField="Em";
-  A2AAslashVectorsPar.solver="CG_s";
-  A2AAslashVectorsPar.output="AslashSeq";
-  application.createModule<MSolver::A2AAslashVectors>("Aslash_seq",A2AAslashVectorsPar);
-}
-/////////////////////////////////////////////////////////////
-// Aslash_perambulators
-/////////////////////////////////////////////////////////////
-void test_PerambulatorsSolve(Application &application)
-{
-  // Perambulator parameters
-  MDistil::PerambFromSolve::Par PerambFromSolvePar;
-  PerambFromSolvePar.eigenPack="LapEvec";
-  PerambFromSolvePar.solve="Aslash_seq";
-  PerambFromSolvePar.PerambFileName="perambAslashS.bin";
-  PerambFromSolvePar.Distil.tsrc = 0;
-  PerambFromSolvePar.Distil.nnoise = 1;
-  PerambFromSolvePar.nvec=5;
-  application.createModule<MDistil::PerambFromSolve>("PerambAslashS",PerambFromSolvePar);
-}
-
 bool bNumber( int &ri, const char * & pstr, bool bGobbleWhiteSpace = true )
 {
   if( bGobbleWhiteSpace )
@@ -515,398 +282,8 @@ bool bNumber( int &ri, const char * & pstr, bool bGobbleWhiteSpace = true )
   return true;
 }
 
-#ifdef DEBUG
-
-typedef Grid::Hadrons::MDistil::NamedTensor<Complex,3,sizeof(Real)> MyTensor;
-
-template<typename T>
-typename std::enable_if<Grid::EigenIO::is_tensor<T>::value && !Grid::Hadrons::MDistil::is_named_tensor<T>::value>::type
-DebugShowTensor(T &x, const char * n, std::string * pIndexNames=nullptr)
-{
-  const MyTensor::Index s{x.size()};
-  std::cout << n << ".size() = " << s << std::endl;
-  std::cout << n << ".NumDimensions = " << x.NumDimensions << " (TensorBase)" << std::endl;
-  std::cout << n << ".NumIndices = " << x.NumIndices << std::endl;
-  const auto d{x.dimensions()};
-  //std::cout << n << ".dimensions().size() = " << d.size() << std::endl;
-  std::cout << "Dimensions are ";
-  for(auto i = 0; i < x.NumDimensions ; i++)
-    std::cout << "[" << d[i] << "]";
-  std::cout << std::endl;
-  MyTensor::Index SizeCalculated{1};
-  std::cout << "Dimensions again";
-  for(int i=0 ; i < x.NumDimensions ; i++ ) {
-    std::cout << " : [" << i;
-    if( pIndexNames )
-      std::cout << ", " << pIndexNames[i];
-    std::cout << "]=" << x.dimension(i);
-    SizeCalculated *= d[i];
-  }
-  std::cout << std::endl;
-  std::cout << "SizeCalculated = " << SizeCalculated << std::endl;\
-  assert( SizeCalculated == s );
-  // Initialise
-  assert( x.NumDimensions == 3 );
-  for( int i = 0 ; i < d[0] ; i++ )
-    for( int j = 0 ; j < d[1] ; j++ )
-      for( int k = 0 ; k < d[2] ; k++ ) {
-        x(i,j,k) = std::complex<double>(SizeCalculated, -SizeCalculated);
-        SizeCalculated--;
-      }
-  // Show raw data
-  std::cout << "Data follow : " << std::endl;
-  typename T::Scalar * p = x.data();
-  for( auto i = 0 ; i < s ; i++ ) {
-    if( i ) std::cout << ", ";
-    std::cout << n << ".data()[" << i << "]=" << * p++;
-  }
-  std::cout << std::endl;
-}
-
-template<typename T>
-typename std::enable_if<Grid::Hadrons::MDistil::is_named_tensor<T>::value>::type
-DebugShowTensor(T &x, const char * n)
-{
-  DebugShowTensor( x.tensor, n, &x.IndexNames[0] );
-}
-
-// Test whether typedef and underlying types are the same
-
-void DebugTestTypeEqualities(void)
-{
-  Real    r1;
-  RealD   r2;
-  double  r3;
-  const std::type_info &tr1{typeid(r1)};
-  const std::type_info &tr2{typeid(r2)};
-  const std::type_info &tr3{typeid(r3)};
-  if( tr1 == tr2 && tr2 == tr3 )
-    std::cout << "r1, r2 and r3 are the same type" << std::endl;
-  else
-    std::cout << "r1, r2 and r3 are different types" << std::endl;
-  std::cout << "r1 is a " << tr1.name() << std::endl;
-  std::cout << "r2 is a " << tr2.name() << std::endl;
-  std::cout << "r3 is a " << tr3.name() << std::endl;
-  
-  // These are the same
-  Complex             c1;
-  std::complex<Real>  c2;
-  const std::type_info &tc1{typeid(c1)};
-  const std::type_info &tc2{typeid(c2)};
-  const std::type_info &tc3{typeid(SpinVector::scalar_type)};
-  if( tc1 == tc2 && tc2 == tc3)
-    std::cout << "c1, c2 and SpinVector::scalar_type are the same type" << std::endl;
-  else
-    std::cout << "c1, c2 and SpinVector::scalar_type are different types" << std::endl;
-  std::cout << "c1                      is a " << tc1.name() << std::endl;
-  std::cout << "c2                      is a " << tc2.name() << std::endl;
-  std::cout << "SpinVector::scalar_type is a " << tc3.name() << std::endl;
-  
-  // These are the same
-  SpinVector                              s1;
-  iSpinVector<Complex >                   s2;
-  iScalar<iVector<iScalar<Complex>, Ns> > s3;
-  const std::type_info &ts1{typeid(s1)};
-  const std::type_info &ts2{typeid(s2)};
-  const std::type_info &ts3{typeid(s3)};
-  if( ts1 == ts2 && ts2 == ts3 )
-    std::cout << "s1, s2 and s3 are the same type" << std::endl;
-  else
-    std::cout << "s1, s2 and s3 are different types" << std::endl;
-  std::cout << "s1 is a " << ts1.name() << std::endl;
-  std::cout << "s2 is a " << ts2.name() << std::endl;
-  std::cout << "s3 is a " << ts3.name() << std::endl;
-  
-  // These are the same
-  SpinColourVector            sc1;
-  iSpinColourVector<Complex > sc2;
-  const std::type_info &tsc1{typeid(sc1)};
-  const std::type_info &tsc2{typeid(sc2)};
-  if( tsc1 == tsc2 )
-    std::cout << "sc1 and sc2 are the same type" << std::endl;
-  else
-    std::cout << "sc1 and sc2 are different types" << std::endl;
-  std::cout << "sc1 is a " << tsc1.name() << std::endl;
-  std::cout << "sc2 is a " << tsc2.name() << std::endl;
-}
-
-bool DebugEigenTest()
-{
-  {
-    Eigen::TensorFixedSize<std::complex<double>,Eigen::Sizes<3,4,5>> x;
-    DebugShowTensor(x, "fixed");
-  }
-  const char pszTestFileName[] = "test_tensor.bin";
-  std::array<std::string,3> as={"Alpha", "Beta", "Gamma"};
-  MyTensor x(as, 2,1,4);
-  DebugShowTensor(x, "x");
-  x.write(pszTestFileName);
-  // Test initialisation of an array of strings
-  for( auto a : as )
-    std::cout << a << std::endl;
-  Grid::Hadrons::MDistil::NamedTensor<Complex,3,sizeof(Real)> p{as,2,7,2};
-  DebugShowTensor(p, "p");
-  std::cout << "p.IndexNames follow" << std::endl;
-  for( auto a : p.IndexNames )
-    std::cout << a << std::endl;
-
-  // 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.read(pszTestFileName);
-  DebugShowTensor(y, "y");
-
-  // Now see whether we can read a tensor back from an hdf5 file
-  const char * pszFileName = "test";
-  y.write(pszFileName);
-  {
-    MyTensor z;
-    const char * pszName = "z1";
-    DebugShowTensor(z, pszName);
-    z.read(pszFileName);
-    DebugShowTensor(z, pszName);
-  }
-  {
-    MyTensor z(Names2,2,0,0);
-    const char * pszName = "z2";
-    DebugShowTensor(z, pszName);
-    z.read(pszFileName);
-    DebugShowTensor(z, pszName);
-  }
-  {
-    // Now see whether we can read a tensor back from an xml file
-    const char * pszXmlName = "test.xml";
-    {
-      XmlWriter w(pszXmlName);
-      y.write<XmlWriter>(w);
-    }
-    MyTensor z;
-    const char * pszName = "xml1";
-    DebugShowTensor(z, pszName);
-    XmlReader r(pszXmlName);
-    z.read<XmlReader>(r);
-    DebugShowTensor(z, pszName);
-  }
-  if((0)) // The following tests would fail
-  {
-    MyTensor z(Names2,2,0,78);
-    //std::array<std::string,3> NamesBad={"Alpha", "Gamma", "Kilo"};
-    //MyTensor z(NamesBad);
-    const char * pszName = "zFail";
-    DebugShowTensor(z, pszName);
-    z.read(pszFileName);
-    DebugShowTensor(z, pszName);
-  }
-
-  // Testing whether typedef produces the same type - yes it does
-
-  DebugTestTypeEqualities();
-  std::cout << std::endl;
-
-  // How to access members of SpinColourVector
-  SpinColourVector  sc;
-  for( int s = 0 ; s < Ns ; s++ ) {
-    auto cv{sc()(s)};
-    iVector<Complex,Nc> c2{sc()(s)};
-    std::cout << " cv is a " << typeid(cv).name() << std::endl;
-    std::cout << " c2 is a " << typeid(c2).name() << std::endl;
-    for( int c = 0 ; c < Nc ; c++ ) {
-      Complex & z{cv(c)};
-      std::cout << "  sc[spin=" << s << ", colour=" << c << "] = " << z << std::endl;
-    }
-  }
-  // We could have removed the Lorentz index independently, but much easier to do as we do above
-  iVector<iVector<Complex,Nc>,Ns> sc2{sc()};
-  std::cout << "sc() is a " << typeid(sc()).name() << std::endl;
-  std::cout << "sc2  is a " << typeid(sc2 ).name() << std::endl;
-
-  // Or you can access elements directly
-  std::complex<Real> z = sc()(0)(0);
-  std::cout << "z = " << z << std::endl;
-  sc()(3)(2) = std::complex<Real>{3.141,-3.141};
-  std::cout << "sc()(3)(2) = " << sc()(3)(2) << std::endl;
-
-  return true;
-}
-
-template <typename T>
-void DebugGridTensorTest_print( int i )
-{
-  // std::cout << i << " : " << EigenIO::is_tensor<T>::value
-  // << ", Rank " << EigenIO::Traits<T>::Rank
-  // << ", count " << EigenIO::Traits<T>::count
-  // << std::endl;
-}
-
-// begin() and end() are the minimum necessary to support range-for loops
-// should really turn this into an iterator ...
-template<typename T, int N>
-class TestObject {
-public:
-  using value_type = T;
-private:
-  value_type * m_p;
-public:
-  TestObject() {
-    m_p = reinterpret_cast<value_type *>(std::malloc(N * sizeof(value_type)));
-  }
-  ~TestObject() { std::free(m_p); }
-  inline value_type * begin(void) { return m_p; }
-  inline value_type * end(void) { return m_p + N; }
-};
-
-template<typename ET> typename std::enable_if<EigenIO::is_tensor<ET>::value>::type
-dump_tensor(const ET & et, const char * psz = nullptr) {
-  if( psz )
-    std::cout << psz << ": ";
-  else
-    std::cout << "Unnamed tensor: ";
-  Serializable::WriteMember( std::cout, et );
-}
-
-template <int Options>
-void EigenSliceExample()
-{
-  std::cout << "Eigen example, Options = " << Options << std::endl;
-  using T2 = Eigen::Tensor<int, 2, Options>;
-  T2 a(4, 3);
-  a.setValues({{0, 100, 200}, {300, 400, 500},
-    {600, 700, 800}, {900, 1000, 1100}});
-  std::cout << "a\n" << a << std::endl;
-  dump_tensor( a, "a" );
-  Eigen::array<typename T2::Index, 2> offsets = {0, 1};
-  Eigen::array<typename T2::Index, 2> extents = {4, 2};
-  T2 slice = a.slice(offsets, extents);
-  std::cout << "slice\n" << slice << std::endl;
-  dump_tensor( slice, "slice" );
-  std::cout << "\n========================================" << std::endl;
-}
-
-template <int Options>
-void EigenSliceExample2()
-{
-  using TestScalar = std::complex<float>;
-  using T3 = Eigen::Tensor<TestScalar, 3, Options>;
-  using T2 = Eigen::Tensor<TestScalar, 2, Options>;
-  T3 a(2,3,4);
-
-  std::cout << "Initialising a:";
-  TestScalar f{ 0 };
-  const TestScalar Inc{ 1, -1 };
-  for( auto &c : a ) {
-    c = f;
-    f += Inc;
-  }
-  std::cout << std::endl;
-  std::cout << "Validating   a (Eigen::" << ( ( Options & Eigen::RowMajor ) ? "Row" : "Col" ) << "Major):" << std::endl;
-  f = 0;
-  for( int i = 0 ; i < a.dimension(0) ; i++ )
-    for( int j = 0 ; j < a.dimension(1) ; j++ )
-      for( int k = 0 ; k < a.dimension(2) ; k++ ) {
-        std::cout << " a(" << i << "," << j << "," << k << ")=" << a(i,j,k) << std::endl;
-        assert( ( Options & Eigen::RowMajor ) == 0 || a(i,j,k) == f );
-        f += Inc;
-      }
-  //std::cout << std::endl;
-  //std::cout << "a initialised to:\n" << a << std::endl;
-  dump_tensor( a, "a" );
-  std::cout << std::endl;
-  Eigen::array<typename T3::Index, 3> offsets = {0,1,1};
-  Eigen::array<typename T3::Index, 3> extents = {1,2,2};
-  T3 b;
-  b = a.slice( offsets, extents );//.reshape(NewExtents);
-  std::cout << "b = a.slice( offsets, extents ):\n" << b << std::endl;
-  dump_tensor( b, "b" );
-  T2 c(3,4);
-  c = a.chip(0,1);
-  std::cout << "c = a.chip(0,0):\n" << c << std::endl;
-  dump_tensor( c, "c" );
-  //T2 d = b.reshape(extents);
-  //std::cout << "b.reshape(extents) is:\n" << d << std::endl;
-  std::cout << "\n========================================" << std::endl;
-}
-
-void DebugFelixTensorTest( void )
-{
-  unsigned int Nmom = 2;
-  unsigned int Nt = 2;
-  unsigned int N_1 = 2;
-  unsigned int N_2 = 2;
-  unsigned int N_3 = 2;
-  using BaryonTensorSet = Eigen::Tensor<Complex, 6, Eigen::RowMajor>;
-  BaryonTensorSet BField3(Nmom,4,Nt,N_1,N_2,N_3);
-  std::vector<Complex> Memory(Nmom * Nt * N_1 * N_2 * N_3 * 2);
-  using BaryonTensorMap = Eigen::TensorMap<BaryonTensorSet>;
-  BaryonTensorMap BField4 (&Memory[0], Nmom,4,Nt,N_1,N_2,N_3);
-
-  EigenSliceExample<Eigen::RowMajor>();
-  EigenSliceExample<0>();
-  EigenSliceExample2<Eigen::RowMajor>();
-  EigenSliceExample2<0>();
-}
-
-bool DebugGridTensorTest( void )
-{
-  DebugFelixTensorTest();
-  typedef Complex t1;
-  typedef iScalar<t1> t2;
-  typedef iVector<t1, Ns> t3;
-  typedef iMatrix<t1, Nc> t4;
-  typedef iVector<iMatrix<t1,1>,4> t5;
-  typedef iScalar<t5> t6;
-  typedef iMatrix<t6, 3> t7;
-  typedef iMatrix<iVector<iScalar<t7>,4>,2> t8;
-  int i = 1;
-  DebugGridTensorTest_print<t1>( i++ );
-  DebugGridTensorTest_print<t2>( i++ );
-  DebugGridTensorTest_print<t3>( i++ );
-  DebugGridTensorTest_print<t4>( i++ );
-  DebugGridTensorTest_print<t5>( i++ );
-  DebugGridTensorTest_print<t6>( i++ );
-  DebugGridTensorTest_print<t7>( i++ );
-  DebugGridTensorTest_print<t8>( i++ );
-
-  //using TOC7 = TestObject<std::complex<double>, 7>;
-  using TOC7 = t7;
-  TOC7 toc7;
-  constexpr std::complex<double> Inc{1,-1};
-  std::complex<double> Start{Inc};
-  for( auto &x : toc7 ) {
-    x = Start;
-    Start += Inc;
-  }
-  i = 0;
-  std::cout << "toc7:";
-  for( auto x : toc7 ) std::cout << " [" << i++ << "]=" << x;
-  std::cout << std::endl;
-
-  //t2 o2;
-  //auto a2 = TensorRemove(o2);
-  //t3 o3;
-  //t4 o4;
-  //auto a3 = TensorRemove(o3);
-  //auto a4 = TensorRemove(o4);
-  
-  return true;
-}
-
-bool ConvertPeramb(const char * pszSource, const char * pszDest) {
-  Grid::Hadrons::MDistil::PerambTensor p(Hadrons::MDistil::PerambIndexNames);
-  p.ReadBinary( pszSource );
-  p.write(pszDest);
-  return true;
-}
-#endif
-
 int main(int argc, char *argv[])
 {
-#ifdef DEBUG
-  // Debug only - test of Eigen::Tensor
-  //if( DebugEigenTest() ) return 0;
-  //if(DebugGridTensorTest()) return 0;
-  //if(ConvertPeramb("PerambL_100_tsrc0.3000","PerambL_100_tsrc0.3000")) return 0;
-#endif
-
   // Decode command-line parameters. 1st one is which test to run
   int iTestNum = -1;
 
@@ -948,31 +325,10 @@ int main(int argc, char *argv[])
 
   // For now perform free propagator test - replace this with distillation test(s)
   LOG(Message) << "====== Creating xml for test " << iTestNum << " ======" << std::endl;
-  //const unsigned int  nt    = GridDefaultLatt()[Tp];
   
   switch(iTestNum) {
-    case 0:
-      test_Global( application );
-      test_LapEvec( application );
-      break;
-    case 1:
-      test_Global( application );
-      test_LapEvec( application );
-      test_Perambulators( application );
-      break;
-    default: // 2
-      test_Global( application );
-      test_LapEvec( application );
-      test_Perambulators( application );
-      test_DistilVectors( application );
-      break;
-    case 3:
-      test_Global( application );
-      test_LapEvec( application );
-      test_LoadPerambulators( application );
-      test_DistilVectors( application );
-      break;
-    case 4:
+    default: // 0
+      LOG(Message) << "Computing Meson 2pt-function" << std::endl;
       test_Global( application );
       test_LapEvec( application );
       test_Perambulators( application );
@@ -980,7 +336,17 @@ int main(int argc, char *argv[])
       test_MesonField( application, "Phi", "_phi" );
       test_MesonField( application, "Rho", "_rho" );
       break;
-    case 5:
+    case 1:
+      LOG(Message) << "Computing Meson 2pt-function by loading perambulators" << std::endl;
+      test_Global( application );
+      test_LapEvec( application );
+      test_LoadPerambulators( application );
+      test_DistilVectors( application );
+      test_MesonField( application, "Phi", "_phi" );
+      test_MesonField( application, "Rho", "_rho" );
+      break;
+    case 2:
+      LOG(Message) << "Computing Meson 2pt-function for two quark flavours" << std::endl;
       test_Global( application );
       test_LapEvec( application );
       test_Perambulators( application );
@@ -990,69 +356,13 @@ int main(int argc, char *argv[])
       test_MesonField( application, "SPhi", "S_phi" );
       test_MesonField( application, "SRho", "S_rho" );
       break;
-#ifdef DISTIL_PRE_RELEASE
-    case 6: // 3
+    case 3:
+      LOG(Message) << "Computing Meson 2pt-function with current insertion" << std::endl;
       test_Global( application );
       test_LapEvec( application );
       test_Perambulators( application );
-      test_g5_sinks( application );
       test_MesonSink( application );
       break;
-    case 7: // 3
-      test_Global( application );
-      test_LapEvec( application );
-      test_Perambulators( application );
-      test_DistilVectors( application );
-      test_BaryonFieldPhi( application );
-      test_BaryonFieldRho( application );
-      break;
-#endif
-    case 8: // 3
-      test_Global( application );
-      test_LapEvec( application );
-      test_Perambulators( application );
-      test_DistilVectors( application );
-      test_MesonField( application, "Phi", "_phi" );
-      test_MesonField( application, "Rho", "_rho" );
-      break;
-#ifdef DISTIL_PRE_RELEASE
-    case 9: // 3
-      test_Global( application );
-      test_Solver( application );
-      test_Baryon2pt( application );
-      break;
-    case 10: // 3
-      test_Global( application );
-      test_LapEvec( application );
-      test_Perambulators( application );
-      test_g5_sinks( application );
-      test_em( application );
-      test_Aslash( application );
-      break;
-    case 11: // 3
-      test_Global( application );
-      test_LapEvec( application );
-      test_Perambulators( application );
-      test_DistilVectors( application );
-      test_BaryonFieldPhi2( application );
-      test_BaryonFieldRho2( application );
-      break;
-#endif
-    case 12: // 3
-      test_Global( application );
-      test_LapEvec( application );
-      test_Perambulators( application, "S" );
-      test_em( application );
-      test_AslashSeq( application );
-      test_PerambulatorsSolve( application );
-      test_DistilVectorsSS( application, "AslashSeq", nullptr, "S" );
-      test_MesonField( application, "AslashSeq" );
-      break;
-    case 13:
-      test_Global( application );
-      test_LapEvec( application );
-      test_MultiPerambulators( application );
-      break;
   }
   // execution
   static const char XmlFileName[] = "test_distil.xml";