Merge branch 'feature/distil' of github.com:mmphys/Grid into feature/distil

Grid/qcd/utils/A2Autils.h

@@ -41,12 +41,21 @@ public:
 			     const std::vector<ComplexField > &mom,
 			     int orthogdim);
 
+
+  static void NucleonFieldMom(Eigen::Tensor<ComplexD,6> &mat, 
+				     const FermionField *one,
+				     const FermionField *two,
+				     const FermionField *three,
+				     const std::vector<ComplexField > &mom,
+				     int parity,
+				     int orthogdim); 
+ 
   static void PionFieldXX(Eigen::Tensor<ComplexD,3> &mat, 
 			  const FermionField *wi,
 			  const FermionField *vj,
 			  int orthogdim,
 			  int g5);
-  
+ 
   static void PionFieldWV(Eigen::Tensor<ComplexD,3> &mat, 
 			  const FermionField *wi,
 			  const FermionField *vj,
@@ -101,6 +110,187 @@ public:
 #endif
 };
 
+template<class FImpl>
+void A2Autils<FImpl>::NucleonFieldMom(Eigen::Tensor<ComplexD,6> &mat, 
+				     const FermionField *one,
+				     const FermionField *two,
+				     const FermionField *three,
+				     const std::vector<ComplexField > &mom,
+				     int parity,
+				     int orthogdim) 
+{
+  assert(parity == 1 || parity == -1);
+
+  typedef typename FImpl::SiteSpinor vobj;
+
+  typedef typename vobj::scalar_object sobj;
+  typedef typename vobj::scalar_type scalar_type;
+  typedef typename vobj::vector_type vector_type;
+
+  typedef iSpinVector<vector_type> SpinVector_v;
+  typedef iSpinVector<scalar_type> SpinVector_s;
+  
+  int oneBlock = mat.dimension(2); 
+  int twoBlock = mat.dimension(3);
+  int threeBlock = mat.dimension(4);
+
+  GridBase *grid = wi[0]._grid;
+  
+  const int    nd = grid->_ndimension;
+  const int Nsimd = grid->Nsimd();
+
+  int Nt     = grid->GlobalDimensions()[orthogdim];
+  int Nmom   = mom.size();
+
+  int fd=grid->_fdimensions[orthogdim];
+  int ld=grid->_ldimensions[orthogdim];
+  int rd=grid->_rdimensions[orthogdim];
+
+  // will locally sum vectors first
+  // sum across these down to scalars
+  // splitting the SIMD
+  int MFrvol = rd*oneBlock*twoBlock*threeBlock*Nmom;
+  int MFlvol = ld*oneBlock*twoBlock*threeBlock*Nmom;
+
+  Vector<SpinVector_v > lvSum(MFrvol);
+  parallel_for (int r = 0; r < MFrvol; r++){
+    lvSum[r] = zero;
+  }
+
+  Vector<SpinVector_s > lsSum(MFlvol);             
+  parallel_for (int r = 0; r < MFlvol; r++){
+    lsSum[r]=scalar_type(0.0);
+  }
+
+  int e1=    grid->_slice_nblock[orthogdim];
+  int e2=    grid->_slice_block [orthogdim];
+  int stride=grid->_slice_stride[orthogdim];
+
+  parallel_for(int r=0;r<rd;r++){
+
+    int so=r*grid->_ostride[orthogdim]; // base offset for start of plane 
+
+    for(int n=0;n<e1;n++){
+      for(int b=0;b<e2;b++){
+
+	int ss= so+n*stride+b;
+
+	for(int i=0;i<oneBlock;i++){
+
+	  auto v1 = one[i]._odata[ss];
+	  auto pv1 = 0.5*(double)parity*(v1 + Gamma(Gamma::Algebra::GammaT)*v1); 
+
+       	  for(int j=0;j<twoBlock;j++){
+
+	    auto v2 = conjugate(two[j]._odata[ss]);
+
+	    for(int k=0;k<threeBlock;k++){
+
+	      auto v3 = three[k]._odata[ss];
+              // C = i gamma_2 gamma_4 => C gamma_5 = - i gamma_1 gamma_3
+	      auto gv3 = Gamma(Gamma::Algebra::SigmaXZ) * v3;
+	      SpinVector_v vv;
+            
+	      vv()()() =  pv1()()(0) * v2()()(1) * gv3()()(2)   //Cross product
+                -         pv1()()(0) * v2()()(2) * gv3()()(1)    
+                +         pv1()()(1) * v2()()(2) * gv3()()(0)    
+                -         pv1()()(1) * v2()()(0) * gv3()()(2)    
+                +         pv1()()(2) * v2()()(0) * gv3()()(1)    
+                -         pv1()()(2) * v2()()(1) * gv3()()(0);    
+
+	    
+	      // After getting the sitewise product do the mom phase loop
+	      int base = Nmom*i+Nmom*Lblock*j+Nmom*Lblock*Rblock*r;
+	      for ( int m=0;m<Nmom;m++){
+	        int idx = m+base;
+	        auto phase = mom[m]._odata[ss];
+	        mac(&lvSum[idx],&vv,&phase()()());
+              }
+	    }
+	  }
+	}
+      }
+    }
+  }
+
+
+  // Sum across simd lanes in the plane, breaking out orthog dir.
+  parallel_for(int rt=0;rt<rd;rt++){
+
+    std::vector<int> icoor(nd);
+    iScalar<vector_type> temp; 
+    std::vector<iScalar<SpinVector_s> > extracted(Nsimd);               
+
+    for(int i=0;i<oneBlock;i++){
+    for(int j=0;j<twoBlock;j++){
+    for(int k=0;k<threeBlock;k++){
+    for(int m=0;m<Nmom;m++){
+
+      int ij_rdx = m+Nmom*i + Nmom*oneBlock * j + Nmom*oneBlock * twoBlock * k + Nmom*oneBlock * twoBlock *threeBlock * rt;
+
+      temp._internal = lvSum[ij_rdx];
+      extract(temp,extracted);
+
+      for(int idx=0;idx<Nsimd;idx++){
+
+	grid->iCoorFromIindex(icoor,idx);
+
+	int ldx    = rt+icoor[orthogdim]*rd;
+
+	int ij_ldx = m+Nmom*i + Nmom*oneBlock * j + Nmom*oneBlock * twoBlock * k + Nmom*oneBlock * twoBlock *threeBlock * ldx;
+
+	lsSum[ij_ldx]=lsSum[ij_ldx]+extracted[idx]._internal;
+
+      }
+    }}}}
+  }
+
+  assert(mat.dimension(0) == Nmom);
+  assert(mat.dimension(1) == Nt);
+
+  int pd = grid->_processors[orthogdim];
+  int pc = grid->_processor_coor[orthogdim];
+  parallel_for_nest2(int lt=0;lt<ld;lt++)
+  {
+    for(int pt=0;pt<pd;pt++){
+      int t = lt + pt*ld;
+      if (pt == pc){
+        for(int i=0;i<oneBlock;i++){
+          for(int j=0;j<twoBlock;j++){
+            for(int k=0;k<threeBlock;k++){
+	      for(int m=0;m<Nmom;m++){
+	        int ij_dx = m+Nmom*i + Nmom*oneBlock * j + Nmom*oneBlock * twoBlock * k + Nmom*oneBlock * twoBlock *threeBlock * lt;
+	        for(int is=0;is<4;is++){
+	          mat(m,t,i,j,k,is) = lsSum[ij_dx]()(is)();
+                }
+              }
+	    }
+	  }
+	}
+      } else { 
+	const scalar_type zz(0.0);
+        for(int i=0;i<oneBlock;i++){
+          for(int j=0;j<twoBlock;j++){
+            for(int k=0;k<threeBlock;k++){
+	      for(int m=0;m<Nmom;m++){
+	        for(int is=0;is<4;is++){
+	          mat(m,t,i,j,k,is) =zz;
+                }
+              }
+	    }
+	  }
+	}
+      }
+    }
+  }
+
+  grid->GlobalSumVector(&mat(0,0,0,0,0,0),Nmom*Nt*oneBlock*twoBlock*threeBlock);
+}
+
+
+
+
+
 /*
 template <class FImpl>
 template <typename TensorType>
@@ -122,6 +312,8 @@ void A2Autils<FImpl>::BaryonField(TensorType &mat,
   typedef iSpinMatrix<vector_type> SpinMatrix_v;
   typedef iSpinMatrix<scalar_type> SpinMatrix_s;
   
+  typedef iSpinColourMatrix<vector_type> SpinColourMatrix_v;
+
   int oneBlock = mat.dimension(3); 
   int twoBlock = mat.dimension(4);
   int threeBlock = mat.dimension(5);
@@ -143,17 +335,23 @@ void A2Autils<FImpl>::BaryonField(TensorType &mat,
   // will locally sum vectors first
   // sum across these down to scalars
   // splitting the SIMD
-  int MFrvol = rd*oneBlock*twoBlock*threeBlock*Nmom;
+  int MFrvol = rd*twoBlock*threeBlock*Nmom;
-  int MFlvol = ld*oneBlock*twoBlock*threeBlock*Nmom;
+  int MFlvol = ld*twoBlock*threeBlock*Nmom;
 
-  Vector<SpinMatrix_v > lvSum(MFrvol);
+  Vector<Vector<SpinMatrix_v >> lvSum(3);
-  parallel_for (int r = 0; r < MFrvol; r++){
+  for (int ic=0;ic<3;ic++){
-    lvSum[r] = zero;
+    lvSum[ic].resize(MFrvol);
+    parallel_for (int r = 0; r < MFrvol; r++){
+      lvSum[ic][r] = zero;
+    }
   }
 
-  Vector<SpinMatrix_s > lsSum(MFlvol);             
+  Vector<Vector<SpinMatrix_s >> lsSum(3);
-  parallel_for (int r = 0; r < MFlvol; r++){
+  for (int ic=0;ic<3;ic++){
-    lsSum[r]=scalar_type(0.0);
+    lsSum[ic].resize(MFlvol);
+    parallel_for (int r = 0; r < MFlvol; r++){
+      lsSum[ic][r] = scalar_type(0.0);
+    }
   }
 
   int e1=    grid->_slice_nblock[orthogdim];
@@ -180,21 +378,26 @@ void A2Autils<FImpl>::BaryonField(TensorType &mat,
 
 	    auto three_k = three[j]._odata[ss];
 
-	    SpinMatrix_v vv;
+	    Vector<SpinMatrix_v > vv(3);
 
 	    for(int s1=0;s1<Ns;s1++){
 	    for(int s2=0;s2<Ns;s2++){
-	      vv()(s1,s2)() = two_j()(s2)(0) * three_k()(s1)(0)   //make this a colorMatrix for the diquark???
+	      vv[0]()(s1,s2)() = two_j()(s2)(1) * three_k()(s1)(2)           //ideal would be SpinMatrix but ColourVector...
-		+             two_j()(s2)(1) * three_k()(s1)(1)
+                -                two_j()(s2)(2) * three_k()(s1)(1);          //this is the cross product (two x three)^i
-		+             two_j()(s2)(2) * three_k()(s1)(2);
+	      vv[1]()(s1,s2)() = two_j()(s2)(2) * three_k()(s1)(0)           
+                -                two_j()(s2)(0) * three_k()(s1)(2);         
+	      vv[2]()(s1,s2)() = two_j()(s2)(0) * three_k()(s1)(1)           
+                -                two_j()(s2)(1) * three_k()(s1)(0);         
 	    }}
 	    
 	    // After getting the sitewise product do the mom phase loop
 	    int base = Nmom*i+Nmom*Lblock*j+Nmom*Lblock*Rblock*r;
 	    for ( int m=0;m<Nmom;m++){
-	      int idx = m+base;
+	      for ( int ic=0;ic<3;ic++){
-	      auto phase = mom[m]._odata[ss];
+	        int idx = m+base;
-	      mac(&lvSum[idx],&vv,&phase);
+	        auto phase = mom[m]._odata[ss];
+	        mac(&lvSum[ic][idx],&vv,&phase);
+	      }
 	    }
 	 
 	  }
@@ -204,19 +407,20 @@ void A2Autils<FImpl>::BaryonField(TensorType &mat,
   }
 
 
+  for ( int ic=0;ic<3;ic++){
   // Sum across simd lanes in the plane, breaking out orthog dir.
   parallel_for(int rt=0;rt<rd;rt++){
 
     std::vector<int> icoor(Nd);
     std::vector<SpinMatrix_s> extracted(Nsimd);               
 
-    for(int i=0;i<Lblock;i++){
+    for(int i=0;i<twoBlock;i++){
-    for(int j=0;j<Rblock;j++){
+    for(int j=0;j<threeBlock;j++){
     for(int m=0;m<Nmom;m++){
 
       int ij_rdx = m+Nmom*i+Nmom*Lblock*j+Nmom*Lblock*Rblock*rt;
 
-      extract(lvSum[ij_rdx],extracted);
+      extract(lvSum[ic][ij_rdx],extracted);
 
       for(int idx=0;idx<Nsimd;idx++){
 
@@ -226,16 +430,20 @@ void A2Autils<FImpl>::BaryonField(TensorType &mat,
 
 	int ij_ldx = m+Nmom*i+Nmom*Lblock*j+Nmom*Lblock*Rblock*ldx;
 
-	lsSum[ij_ldx]=lsSum[ij_ldx]+extracted[idx];
+	lsSum[ic][ij_ldx]=lsSum[ic][ij_ldx]+extracted[idx];
 
       }
     }}}
   }
+ 
+
   if (t_kernel) *t_kernel += usecond();
   assert(mat.dimension(0) == Nmom);
   assert(mat.dimension(1) == Ngamma);
   assert(mat.dimension(2) == Nt);
 
+  TensorType diquark; // Need this instead of mat!!!
+
   // ld loop and local only??
   int pd = grid->_processors[orthogdim];
   int pc = grid->_processor_coor[orthogdim];
@@ -244,21 +452,21 @@ void A2Autils<FImpl>::BaryonField(TensorType &mat,
     for(int pt=0;pt<pd;pt++){
       int t = lt + pt*ld;
       if (pt == pc){
-	for(int i=0;i<Lblock;i++){
+	for(int i=0;i<twoBlock;i++){
-	  for(int j=0;j<Rblock;j++){
+	  for(int j=0;j<threeBlock;j++){
 	    for(int m=0;m<Nmom;m++){
 	      int ij_dx = m+Nmom*i + Nmom*Lblock * j + Nmom*Lblock * Rblock * lt;
 	      for(int mu=0;mu<Ngamma;mu++){
 		// this is a bit slow
-		mat(m,mu,t,i,j) = trace(lsSum[ij_dx]*Gamma(gammas[mu]));
+		mat(m,mu,t,i,j) = trace(lsSum[ic][ij_dx]*Gamma(gammaB[mu]));
 	      }
 	    }
 	  }
 	}
       } else { 
 	const scalar_type zz(0.0);
-	for(int i=0;i<Lblock;i++){
+	for(int i=0;i<twoBlock;i++){
-	  for(int j=0;j<Rblock;j++){
+	  for(int j=0;j<threeBlock;j++){
 	    for(int mu=0;mu<Ngamma;mu++){
 	      for(int m=0;m<Nmom;m++){
 		mat(m,mu,t,i,j) =zz;
@@ -269,7 +477,7 @@ void A2Autils<FImpl>::BaryonField(TensorType &mat,
       }
     }
   }
-
+  }
   ////////////////////////////////////////////////////////////////////
   // This global sum is taking as much as 50% of time on 16 nodes
   // Vector size is 7 x 16 x 32 x 16 x 16 x sizeof(complex) = 2MB - 60MB depending on volume

Hadrons/Modules/MDistil/LapEvec.hpp

@@ -239,7 +239,7 @@ void TLapEvec<GImpl>::execute(void)
   auto &Umu = envGet(GaugeField, par().gauge);
   envGetTmp(GaugeField, Umu_smear);
   FieldMetaData header;
-  if((1)) {
+  if((0)) {
     const std::vector<int> seeds({1, 2, 3, 4, 5});
     GridParallelRNG pRNG4d(gridHD);
     pRNG4d.SeedFixedIntegers(seeds);

Hadrons/Modules/MDistil/PerambLight.hpp

@@ -212,7 +212,7 @@ void TPerambLight<FImpl>::execute(void)
 
     envGetTmp(GaugeField, Umu);
   FieldMetaData header;
-  if((1)){
+  if((0)){
     const std::vector<int> seeds({1, 2, 3, 4, 5});
     GridParallelRNG pRNG4d(grid4d);
     pRNG4d.SeedFixedIntegers(seeds);

tests/hadrons/Test_24.cc

@@ -0,0 +1,772 @@
+/*************************************************************************************
+ 
+ Grid physics library, www.github.com/paboyle/Grid
+ 
+ Source file: Tests/Hadrons/Test_hadrons_distil.cc
+ 
+ Copyright (C) 2015-2019
+ 
+ Author: Felix Erben <ferben@ed.ac.uk>
+ Author: Michael Marshall <Michael.Marshall@ed.ac.uk>
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+ 
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+ 
+ You should have received a copy of the GNU General Public License along
+ with this program; if not, write to the Free Software Foundation, Inc.,
+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ 
+ See the full license in the file "LICENSE" in the top level distribution directory
+ *************************************************************************************/
+/*  END LEGAL */
+
+#include <typeinfo>
+#include <Hadrons/Application.hpp>
+#include <Hadrons/Modules.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+
+/////////////////////////////////////////////////////////////
+// Test creation of laplacian eigenvectors
+/////////////////////////////////////////////////////////////
+int Nconf = 3160;
+
+void test_Global(Application &application)
+{
+  // global parameters
+  Application::GlobalPar globalPar;
+  globalPar.trajCounter.start = Nconf;
+  globalPar.trajCounter.end   = Nconf + 20;
+  globalPar.trajCounter.step  = 20;
+  globalPar.runId             = "test";
+  application.setPar(globalPar);
+} 
+
+/////////////////////////////////////////////////////////////
+// Test creation of laplacian eigenvectors
+/////////////////////////////////////////////////////////////
+
+void test_LapEvec(Application &application)
+{
+  const char szGaugeName[] = "gauge";
+  // gauge field
+  application.createModule<MGauge::Random>(szGaugeName);
+  // Now make an instance of the LapEvec object
+  MDistil::LapEvecPar p;
+  //p.ConfigFileDir="/home/dp008/dp008/dc-rich6/Scripts/ConfigsDeflQED/";
+  //p.ConfigFileName="ckpoint_lat.3000";
+  p.ConfigFileDir="/home/dp008/dp008/dc-rich6/Scripts/ConfigsDeflQED/";
+  p.ConfigFileName="ckpoint_lat." + std::to_string(Nconf);
+  p.gauge = szGaugeName;
+  //p.EigenPackName = "ePack";
+  //p.Distil.TI = 8;
+  //p.Distil.LI = 3;
+  //p.Distil.Nnoise = 2;
+  //p.Distil.tSrc = 0;
+  p.Stout.steps = 3;
+  p.Stout.parm = 0.2;
+  p.Cheby.PolyOrder = 11;
+  p.Cheby.alpha = 0.3;
+  p.Cheby.beta = 12.5;
+  p.Lanczos.Nvec = 50;
+  p.Lanczos.Nk = 60;
+  p.Lanczos.Np = 20;
+  p.Lanczos.MaxIt = 1000;
+  p.Lanczos.resid = 1e-8;
+  application.createModule<MDistil::LapEvec>("LapEvec",p);
+}
+
+/////////////////////////////////////////////////////////////
+// Perambulators
+/////////////////////////////////////////////////////////////
+
+void test_Perambulators(Application &application)
+{
+  // PerambLight parameters
+  MDistil::PerambLight::Par PerambPar;
+  PerambPar.eigenPack="LapEvec";
+  PerambPar.PerambFileName="peramb_" + std::to_string(Nconf) + ".bin";
+  PerambPar.ConfigFileDir="/home/dp008/dp008/dc-rich6/Scripts/ConfigsDeflQED/";
+  PerambPar.ConfigFileName="ckpoint_lat." + std::to_string(Nconf);
+  PerambPar.UniqueIdentifier="full_dilution";
+  PerambPar.Distil.tsrc = 0;
+  PerambPar.Distil.nnoise = 1;
+  PerambPar.Distil.LI=50;
+  PerambPar.Distil.SI=4;
+  PerambPar.Distil.TI=64;
+  PerambPar.nvec=50;
+  PerambPar.Distil.Ns=4;
+  PerambPar.Distil.Nt=64;
+  PerambPar.Distil.Nt_inv=1;
+  PerambPar.Solver.mass=0.005;
+  PerambPar.Solver.M5=1.8;
+  PerambPar.Ls=16;
+  PerambPar.Solver.CGPrecision=1e-7;
+  PerambPar.Solver.MaxIterations=10000;
+  application.createModule<MDistil::PerambLight>("Peramb",PerambPar);
+}
+/////////////////////////////////////////////////////////////
+// DistilVectors
+/////////////////////////////////////////////////////////////
+
+void test_DistilVectors(Application &application)
+{
+  // DistilVectors parameters
+  MDistil::DistilVectors::Par DistilVecPar;
+  DistilVecPar.noise="Peramb_noise";
+  DistilVecPar.perambulator="Peramb_perambulator_light";
+  DistilVecPar.eigenPack="LapEvec";
+  DistilVecPar.tsrc = 0;
+  DistilVecPar.nnoise = 1;
+  DistilVecPar.LI=50;
+  DistilVecPar.SI=4;
+  DistilVecPar.TI=64;
+  DistilVecPar.nvec=50;
+  DistilVecPar.Ns=4;
+  DistilVecPar.Nt=64;
+  DistilVecPar.Nt_inv=1;
+  application.createModule<MDistil::DistilVectors>("DistilVecs",DistilVecPar);
+}
+void test_PerambulatorsS(Application &application)
+{
+  // PerambLight parameters
+  MDistil::PerambLight::Par PerambPar;
+  PerambPar.eigenPack="LapEvec";
+  PerambPar.PerambFileName="perambS.bin";
+  PerambPar.ConfigFileDir="/home/dp008/dp008/paboyle/A2A/run/";
+  PerambPar.ConfigFileName="ckpoint_lat.IEEE64BIG.1100";
+  PerambPar.UniqueIdentifier="full_dilution";
+  PerambPar.Distil.tsrc = 0;
+  PerambPar.Distil.nnoise = 1;
+  PerambPar.Distil.LI=50;
+  PerambPar.Distil.SI=4;
+  PerambPar.Distil.TI=64;
+  PerambPar.nvec=50;
+  PerambPar.Distil.Ns=4;
+  PerambPar.Distil.Nt=64;
+  PerambPar.Distil.Nt_inv=1;
+  PerambPar.Solver.mass=0.04; //strange mass???
+  PerambPar.Solver.M5=1.8;
+  PerambPar.Ls=16;
+  PerambPar.Solver.CGPrecision=1e-8;
+  PerambPar.Solver.MaxIterations=10000;
+  application.createModule<MDistil::PerambLight>("PerambS",PerambPar);
+}
+/////////////////////////////////////////////////////////////
+// DistilVectors
+/////////////////////////////////////////////////////////////
+
+void test_DistilVectorsS(Application &application)
+{
+  // DistilVectors parameters
+  MDistil::DistilVectors::Par DistilVecPar;
+  DistilVecPar.noise="PerambS_noise";
+  DistilVecPar.perambulator="PerambS_perambulator_light";
+  DistilVecPar.eigenPack="LapEvec";
+  DistilVecPar.tsrc = 0;
+  DistilVecPar.nnoise = 1;
+  DistilVecPar.LI=50;
+  DistilVecPar.SI=4;
+  DistilVecPar.TI=64;
+  DistilVecPar.nvec=50;
+  DistilVecPar.Ns=4;
+  DistilVecPar.Nt=64;
+  DistilVecPar.Nt_inv=1;
+  application.createModule<MDistil::DistilVectors>("DistilVecsS",DistilVecPar);
+}
+/////////////////////////////////////////////////////////////
+// MesonSink
+/////////////////////////////////////////////////////////////
+
+void test_MesonSink(Application &application)
+{
+  // DistilVectors parameters
+  MContraction::A2AMesonField::Par A2AMesonFieldPar;
+  A2AMesonFieldPar.left="Peramb_unsmeared_sink";
+  A2AMesonFieldPar.right="Peramb_unsmeared_sink";
+  A2AMesonFieldPar.output="DistilFields";
+  A2AMesonFieldPar.gammas="all";
+  A2AMesonFieldPar.mom={"0 0 0"};
+  A2AMesonFieldPar.cacheBlock=2;
+  A2AMesonFieldPar.block=4;
+  application.createModule<MContraction::A2AMesonField>("DistilMesonSink",A2AMesonFieldPar);
+}
+/////////////////////////////////////////////////////////////
+// g5*unsmeared
+/////////////////////////////////////////////////////////////
+void test_g5_sinks(Application &application)
+{
+  MDistil::g5_multiply::Par g5_multiplyPar;
+  g5_multiplyPar.input="Peramb_unsmeared_sink";
+  g5_multiplyPar.nnoise = 1;
+  g5_multiplyPar.LI=50;
+  g5_multiplyPar.Ns=4;
+  g5_multiplyPar.Nt_inv=1;
+  application.createModule<MDistil::g5_multiply>("g5phi",g5_multiplyPar);
+}
+/////////////////////////////////////////////////////////////
+// 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);
+}
+/////////////////////////////////////////////////////////////
+// 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);
+}
+
+/////////////////////////////////////////////////////////////
+// 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)
+{
+   MContraction::A2AAslashField::Par A2AAslashFieldPar;
+   A2AAslashFieldPar.left="Peramb_unsmeared_sink";
+   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);
+}
+
+bool bNumber( int &ri, const char * & pstr, bool bGobbleWhiteSpace = true )
+{
+  if( bGobbleWhiteSpace )
+    while( std::isspace(static_cast<unsigned char>(*pstr)) )
+      pstr++;
+  const char * p = pstr;
+  bool bMinus = false;
+  char c = * p++;
+  if( c == '+' )
+    c = * p++;
+  else if( c == '-' ) {
+    bMinus = true;
+    c = * p++;
+  }
+  int n = c - '0';
+  if( n < 0 || n > 9 )
+    return false;
+  while( * p >= '0' && * p <= '9' ) {
+    n = n * 10 + ( * p ) - '0';
+    p++;
+  }
+  if( bMinus )
+    n *= -1;
+  ri = n;
+  pstr = p;
+  return true;
+}
+
+#ifdef DEBUG
+
+typedef Grid::Hadrons::MDistil::NamedTensor<Complex,3,sizeof(Real)> MyTensor;
+
+template<typename T>
+void DebugShowTensor(T &x, const char * n)
+{
+  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 << /*", " << x.IndexNames[i] << */"]=" << 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;
+}
+
+// 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.WriteBinary(pszTestFileName);
+  DebugShowTensor(x, "x");
+  // Test initialisation of an array of strings
+  for( auto a : as )
+    std::cout << a << std::endl;
+  Grid::Hadrons::MDistil::Perambulator<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.ReadBinary(pszTestFileName);
+  DebugShowTensor(y, "y");
+
+  // 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
+  << ", rank_non_trivial " << EigenIO::Traits<T>::rank_non_trivial
+  << ", count " << EigenIO::Traits<T>::count
+  << ", scalar_size " << EigenIO::Traits<T>::scalar_size
+  << ", size " << EigenIO::Traits<T>::size
+  << 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 <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;
+  DumpMemoryOrder( 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;
+  DumpMemoryOrder( 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:";
+  for_all( a, [&](TestScalar &c, float f, const std::array<size_t,T3::NumIndices> &Dims ){
+    c = TestScalar{f,-f};
+    std::cout << " a(" << Dims[0] << "," << Dims[1] << "," << Dims[2] << ")=" << c;
+  } );
+  std::cout << std::endl;
+  //std::cout << "Validating   a:";
+  float z = 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++ ) {
+        TestScalar w{z, -z};
+        //std::cout << " a(" << i << "," << j << "," << k << ")=" << w;
+        assert( a(i,j,k) == w );
+        z++;
+      }
+  //std::cout << std::endl;
+  //std::cout << "a initialised to:\n" << a << std::endl;
+  DumpMemoryOrder( a, "a" );
+  std::cout << "for_all(a):";
+  for_all( a, [&](TestScalar c, typename T3::Index n, const std::array<size_t,T3::NumIndices> &Dims ){
+    std::cout << " (" << Dims[0] << "," << Dims[1] << "," << Dims[2] << ")<" << n << ">=" << c;
+  } );
+  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;
+  DumpMemoryOrder( b, "b" );
+  T2 c(3,4);
+  c = a.chip(0,1);
+  std::cout << "c = a.chip(0,0):\n" << c << std::endl;
+  DumpMemoryOrder( 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;
+}
+#endif
+
+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) << std::endl;
+  if( DebugEigenTest() ) return 0;
+  if(DebugGridTensorTest()) return 0;
+#endif
+
+  // Decode command-line parameters. 1st one is which test to run
+  int iTestNum = -1;
+
+  for(int i = 1 ; i < argc ; i++ ) {
+    std::cout << "argv[" << i << "]=\"" << argv[i] << "\"" << std::endl;
+    const char * p = argv[i];
+    if( * p == '/' || * p == '-' ) {
+      p++;
+      char c = * p++;
+      switch(toupper(c)) {
+        case 'T':
+          if( bNumber( iTestNum, p ) ) {
+            std::cout << "Test " << iTestNum << " requested";
+            if( * p )
+              std::cout << " (ignoring trailer \"" << p << "\")";
+            std::cout << std::endl;
+          }
+          else
+            std::cout << "Invalid test \"" << &argv[i][2] << "\"" << std::endl;
+          break;
+        default:
+          std::cout << "Ignoring switch \"" << &argv[i][1] << "\"" << std::endl;
+          break;
+      }
+    }
+  }
+
+  // initialization //////////////////////////////////////////////////////////
+  Grid_init(&argc, &argv);
+  HadronsLogError.Active(GridLogError.isActive());
+  HadronsLogWarning.Active(GridLogWarning.isActive());
+  HadronsLogMessage.Active(GridLogMessage.isActive());
+  HadronsLogIterative.Active(GridLogIterative.isActive());
+  HadronsLogDebug.Active(GridLogDebug.isActive());
+  LOG(Message) << "Grid initialized" << std::endl;
+  
+  // run setup ///////////////////////////////////////////////////////////////
+  Application              application;
+
+  // 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) {
+    default:
+      test_Global( application );
+      test_LapEvec( application );
+      test_Perambulators( application );
+      test_MesonSink( application );
+      test_g5_sinks( application );
+      test_em( application );
+      test_Aslash( application );
+      test_DistilVectors( application );
+      test_MesonField( application );
+      test_MesonFieldRho( application );
+      break;
+  }
+  LOG(Message) << "====== XML creation for test " << iTestNum << " complete ======" << std::endl;
+
+  // execution
+  application.saveParameterFile("test_hadrons_distil.xml");
+  application.run();
+  
+  // epilogue
+  LOG(Message) << "Grid is finalizing now" << std::endl;
+  Grid_finalize();
+  
+  return EXIT_SUCCESS;
+}