Merge remote-tracking branch 'refs/remotes/paboyle/feature/hadrons' into feature/hadrons

.gitignore

@@ -47,7 +47,9 @@ Config.h.in
 config.log
 config.status
 .deps
-*.inc
+Make.inc
+eigen.inc
+Eigen.inc
 
 # http://www.gnu.org/software/autoconf #
 ########################################

Makefile.am

@@ -1,10 +1,12 @@
 # additional include paths necessary to compile the C++ library
-SUBDIRS = lib benchmarks tests programs
+SUBDIRS = lib benchmarks tests extras
 
-.PHONY: tests
+include $(top_srcdir)/doxygen.inc
 
 tests: all
 	$(MAKE) -C tests tests
 
+.PHONY: tests doxygen-run doxygen-doc $(DX_PS_GOAL) $(DX_PDF_GOAL)
+
 AM_CXXFLAGS += -I$(top_builddir)/include
 ACLOCAL_AMFLAGS = -I m4

README.md

@@ -116,13 +116,15 @@ If you want to build all the tests at once just use `make tests`.
 - `--with-fftw=<path>`: look for FFTW in the UNIX prefix `<path>`
 - `--enable-lapack[=<path>]`: enable LAPACK support in Lanczos eigensolver. A UNIX prefix containing the library can be specified (optional).
 - `--enable-mkl[=<path>]`: use Intel MKL for FFT (and LAPACK if enabled) routines. A UNIX prefix containing the library can be specified (optional).
-- `--enable-numa`: ???
+- `--enable-numa`: enable NUMA first touch optimisation
 - `--enable-simd=<code>`: setup Grid for the SIMD target `<code>` (default: `GEN`). A list of possible SIMD targets is detailed in a section below.
+- `--enable-gen-simd-width=<size>`: select the size (in bytes) of the generic SIMD vector type (default: 32 bytes).
 - `--enable-precision={single|double}`: set the default precision (default: `double`).
 - `--enable-precision=<comm>`: Use `<comm>` for message passing (default: `none`). A list of possible SIMD targets is detailed in a section below.
 - `--enable-rng={ranlux48|mt19937}`: choose the RNG (default: `ranlux48 `).
 - `--disable-timers`: disable system dependent high-resolution timers.
 - `--enable-chroma`: enable Chroma regression tests.
+- `--enable-doxygen-doc`: enable the Doxygen documentation generation (build with `make doxygen-doc`)
 
 ### Possible communication interfaces
 
@@ -136,7 +138,7 @@ The following options can be use with the `--enable-comms=` option to target dif
 | `mpi3l[-auto]` | MPI communications using MPI 3 shared memory and leader model |
 | `shmem `       | Cray SHMEM communications                                     |
 
-For the MPI interfaces the optional `-auto` suffix instructs the `configure` scripts to determine all the necessary compilation and linking flags. This is done by extracting the informations from the MPI wrapper specified in the environment variable `MPICXX` (if not specified `configure` will scan though a list of default names).
+For the MPI interfaces the optional `-auto` suffix instructs the `configure` scripts to determine all the necessary compilation and linking flags. This is done by extracting the informations from the MPI wrapper specified in the environment variable `MPICXX` (if not specified `configure` will scan though a list of default names). The `-auto` suffix is not supported by the Cray environment wrapper scripts. Use the standard versions instead.  
 
 ### Possible SIMD types
 
@@ -165,6 +167,7 @@ Alternatively, some CPU codenames can be directly used:
 - We currently support AVX512 only for the Intel compiler. Support for GCC and clang will appear in future versions of Grid when the AVX512 support within GCC and clang will be more advanced.
 - For BG/Q only [bgclang](http://trac.alcf.anl.gov/projects/llvm-bgq) is supported. We do not presently plan to support more compilers for this platform.
 - BG/Q performances are currently rather poor. This is being investigated for future versions.
+- The vector size for the `GEN` target can be specified with the `configure` script option `--enable-gen-simd-width`.
 
 ### Build setup for Intel Knights Landing platform
 

benchmarks/Benchmark_comms.cc

@@ -193,6 +193,7 @@ int main (int argc, char ** argv)
     }
   }  
 
+
   Nloop=100;
   std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
   std::cout<<GridLogMessage << "= Benchmarking concurrent STENCIL halo exchange in "<<nmu<<" dimensions"<<std::endl;
@@ -271,5 +272,90 @@ int main (int argc, char ** argv)
       std::cout<<GridLogMessage << lat<<"\t\t"<<Ls<<"\t\t"<<bytes<<"\t\t"<<xbytes/time<<"\t\t"<<bidibytes/time<<std::endl;
     }
   }    
+
+
+
+  Nloop=100;
+  std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
+  std::cout<<GridLogMessage << "= Benchmarking sequential STENCIL halo exchange in "<<nmu<<" dimensions"<<std::endl;
+  std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
+  std::cout<<GridLogMessage << "  L  "<<"\t\t"<<" Ls  "<<"\t\t"<<"bytes"<<"\t\t"<<"MB/s uni"<<"\t\t"<<"MB/s bidi"<<std::endl;
+
+  for(int lat=4;lat<=maxlat;lat+=2){
+    for(int Ls=1;Ls<=16;Ls*=2){
+
+      std::vector<int> latt_size  ({lat*mpi_layout[0],
+      				    lat*mpi_layout[1],
+      				    lat*mpi_layout[2],
+      				    lat*mpi_layout[3]});
+
+      GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
+
+      std::vector<HalfSpinColourVectorD *> xbuf(8);
+      std::vector<HalfSpinColourVectorD *> rbuf(8);
+      Grid.ShmBufferFreeAll();
+      for(int d=0;d<8;d++){
+	xbuf[d] = (HalfSpinColourVectorD *)Grid.ShmBufferMalloc(lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD));
+	rbuf[d] = (HalfSpinColourVectorD *)Grid.ShmBufferMalloc(lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD));
+      }
+
+      int ncomm;
+      int bytes=lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD);
+
+      double start=usecond();
+      for(int i=0;i<Nloop;i++){
+
+	std::vector<CartesianCommunicator::CommsRequest_t> requests;
+
+	ncomm=0;
+	for(int mu=0;mu<4;mu++){
+	
+	  if (mpi_layout[mu]>1 ) {
+	  
+	    ncomm++;
+	    int comm_proc=1;
+	    int xmit_to_rank;
+	    int recv_from_rank;
+	    
+	    Grid.ShiftedRanks(mu,comm_proc,xmit_to_rank,recv_from_rank);
+	    Grid.StencilSendToRecvFromBegin(requests,
+					    (void *)&xbuf[mu][0],
+					    xmit_to_rank,
+					    (void *)&rbuf[mu][0],
+					    recv_from_rank,
+					    bytes);
+	    //	    Grid.StencilSendToRecvFromComplete(requests);
+	    //	    requests.resize(0);
+
+	    comm_proc = mpi_layout[mu]-1;
+	  
+	    Grid.ShiftedRanks(mu,comm_proc,xmit_to_rank,recv_from_rank);
+	    Grid.StencilSendToRecvFromBegin(requests,
+					    (void *)&xbuf[mu+4][0],
+					    xmit_to_rank,
+					    (void *)&rbuf[mu+4][0],
+					    recv_from_rank,
+					    bytes);
+	    Grid.StencilSendToRecvFromComplete(requests);
+	    requests.resize(0);
+	  
+	  }
+	}
+	Grid.Barrier();
+
+      }
+      double stop=usecond();
+
+      double dbytes    = bytes;
+      double xbytes    = Nloop*dbytes*2.0*ncomm;
+      double rbytes    = xbytes;
+      double bidibytes = xbytes+rbytes;
+
+      double time = stop-start; // microseconds
+
+      std::cout<<GridLogMessage << lat<<"\t\t"<<Ls<<"\t\t"<<bytes<<"\t\t"<<xbytes/time<<"\t\t"<<bidibytes/time<<std::endl;
+    }
+  }    
+
   Grid_finalize();
 }

benchmarks/Benchmark_dwf.cc

@@ -57,7 +57,7 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
 
   std::vector<int> latt4 = GridDefaultLatt();
-  const int Ls=16;
+  const int Ls=8;
   GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
   GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
   GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
@@ -138,7 +138,7 @@ int main (int argc, char ** argv)
 
   int ncall =100;
   if (1) {
-
+    FGrid->Barrier();
     Dw.ZeroCounters();
     double t0=usecond();
     for(int i=0;i<ncall;i++){
@@ -147,6 +147,7 @@ int main (int argc, char ** argv)
       __SSC_STOP;
     }
     double t1=usecond();
+    FGrid->Barrier();
     
     double volume=Ls;  for(int mu=0;mu<Nd;mu++) volume=volume*latt4[mu];
     double flops=1344*volume*ncall;
@@ -158,7 +159,7 @@ int main (int argc, char ** argv)
     std::cout<<GridLogMessage << "mflop/s per rank =  "<< flops/(t1-t0)/NP<<std::endl;
     err = ref-result; 
     std::cout<<GridLogMessage << "norm diff   "<< norm2(err)<<std::endl;
-    assert (norm2(err)< 1.0e-5 );
+    assert (norm2(err)< 1.0e-4 );
     Dw.Report();
   }
 
@@ -193,6 +194,7 @@ int main (int argc, char ** argv)
       pokeSite(tmp,ssrc,site);
     }}}}}
     std::cout<<GridLogMessage<< "src norms "<< norm2(src)<<" " <<norm2(ssrc)<<std::endl;
+    FGrid->Barrier();
     double t0=usecond();
     sDw.ZeroCounters();
     for(int i=0;i<ncall;i++){
@@ -201,6 +203,7 @@ int main (int argc, char ** argv)
       __SSC_STOP;
     }
     double t1=usecond();
+    FGrid->Barrier();
     double volume=Ls;  for(int mu=0;mu<Nd;mu++) volume=volume*latt4[mu];
     double flops=1344*volume*ncall;
 
@@ -211,12 +214,12 @@ int main (int argc, char ** argv)
   
     if(0){
       for(int i=0;i< PerformanceCounter::NumTypes(); i++ ){
-  sDw.Dhop(ssrc,sresult,0);
+	sDw.Dhop(ssrc,sresult,0);
-  PerformanceCounter Counter(i);
+	PerformanceCounter Counter(i);
-  Counter.Start();
+	Counter.Start();
-  sDw.Dhop(ssrc,sresult,0);
+	sDw.Dhop(ssrc,sresult,0);
-  Counter.Stop();
+	Counter.Stop();
-  Counter.Report();
+	Counter.Report();
       }
     }
 
@@ -240,7 +243,7 @@ int main (int argc, char ** argv)
       }
     }}}}}
     std::cout<<GridLogMessage<<" difference between normal and simd is "<<sum<<std::endl;
-    assert (sum< 1.0e-5 );
+    assert (sum< 1.0e-4 );
 
 
     if (1) {
@@ -271,6 +274,7 @@ int main (int argc, char ** argv)
       if ( WilsonKernelsStatic::Opt == WilsonKernelsStatic::OptInlineAsm ) std::cout << GridLogMessage<< "* Using Asm Nc=3   WilsonKernels" <<std::endl;
       std::cout << GridLogMessage<< "*********************************************************" <<std::endl;
 
+      FGrid->Barrier();
       sDw.ZeroCounters();
       sDw.stat.init("DhopEO");
       double t0=usecond();
@@ -278,6 +282,7 @@ int main (int argc, char ** argv)
         sDw.DhopEO(ssrc_o, sr_e, DaggerNo);
       }
       double t1=usecond();
+      FGrid->Barrier();
       sDw.stat.print();
 
       double volume=Ls;  for(int mu=0;mu<Nd;mu++) volume=volume*latt4[mu];
@@ -301,7 +306,7 @@ int main (int argc, char ** argv)
 
       error+= norm2(ssrc_o);
       std::cout<<GridLogMessage << "sO norm diff   "<< norm2(ssrc_o)<< "  vec nrm"<<norm2(sr_o) <<std::endl;
-      if(error>1.0e-5) { 
+      if(error>1.0e-4) { 
 	setCheckerboard(ssrc,ssrc_o);
 	setCheckerboard(ssrc,ssrc_e);
 	std::cout<< ssrc << std::endl;
@@ -337,7 +342,7 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage << "norm ref    "<< norm2(ref)<<std::endl;
   err = ref-result; 
   std::cout<<GridLogMessage << "norm diff   "<< norm2(err)<<std::endl;
-  assert(norm2(err)<1.0e-5);
+  assert(norm2(err)<1.0e-4);
   LatticeFermion src_e (FrbGrid);
   LatticeFermion src_o (FrbGrid);
   LatticeFermion r_e   (FrbGrid);
@@ -363,11 +368,13 @@ int main (int argc, char ** argv)
   std::cout << GridLogMessage<< "*********************************************************" <<std::endl;
   {
     Dw.ZeroCounters();
+    FGrid->Barrier();
     double t0=usecond();
     for(int i=0;i<ncall;i++){
       Dw.DhopEO(src_o,r_e,DaggerNo);
     }
     double t1=usecond();
+    FGrid->Barrier();
     
     double volume=Ls;  for(int mu=0;mu<Nd;mu++) volume=volume*latt4[mu];
     double flops=(1344.0*volume*ncall)/2;
@@ -389,14 +396,14 @@ int main (int argc, char ** argv)
 
   err = r_eo-result; 
   std::cout<<GridLogMessage << "norm diff   "<< norm2(err)<<std::endl;
-  assert(norm2(err)<1.0e-5);
+  assert(norm2(err)<1.0e-4);
 
   pickCheckerboard(Even,src_e,err);
   pickCheckerboard(Odd,src_o,err);
   std::cout<<GridLogMessage << "norm diff even  "<< norm2(src_e)<<std::endl;
   std::cout<<GridLogMessage << "norm diff odd   "<< norm2(src_o)<<std::endl;
-  assert(norm2(src_e)<1.0e-5);
+  assert(norm2(src_e)<1.0e-4);
-  assert(norm2(src_o)<1.0e-5);
+  assert(norm2(src_o)<1.0e-4);
 
   Grid_finalize();
 }

benchmarks/Benchmark_dwf_sweep.cc

@@ -69,8 +69,8 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage << "Volume \t\t\tProcs \t Dw \t eoDw \t sDw \t eosDw (Mflop/s)  "<<std::endl;
   std::cout<<GridLogMessage << "=========================================================================="<<std::endl;
 
-  int Lmax=32;
+  int Lmax=16;
-  int dmin=0;
+  int dmin=2;
   if ( getenv("LMAX") ) Lmax=atoi(getenv("LMAX"));
   if ( getenv("DMIN") ) dmin=atoi(getenv("DMIN"));
   for (int L=8;L<=Lmax;L*=2){

benchmarks/Benchmark_mooee.cc

@@ -0,0 +1,183 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./benchmarks/Benchmark_dwf.cc
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: paboyle <paboyle@ph.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 <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+  std::vector<int> latt4 = GridDefaultLatt();
+  const int Ls=16;
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+  std::cout << GridLogMessage << "Making Vec5d innermost grids"<<std::endl;
+  GridCartesian         * sUGrid   = SpaceTimeGrid::makeFourDimDWFGrid(GridDefaultLatt(),GridDefaultMpi());
+  GridRedBlackCartesian * sUrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(sUGrid);
+  GridCartesian         * sFGrid   = SpaceTimeGrid::makeFiveDimDWFGrid(Ls,UGrid);
+  GridRedBlackCartesian * sFrbGrid = SpaceTimeGrid::makeFiveDimDWFRedBlackGrid(Ls,UGrid);
+
+  std::vector<int> seeds4({1,2,3,4});
+  std::vector<int> seeds5({5,6,7,8});
+
+  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
+  std::cout << GridLogMessage << "Seeded"<<std::endl;
+
+  LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
+
+  std::cout << GridLogMessage << "made random gauge fields"<<std::endl;
+
+  RealD mass=0.1;
+  RealD M5  =1.8;
+  RealD NP = UGrid->_Nprocessors;
+
+
+  if (1)
+  {
+    const int ncall=1000;
+
+    std::cout << GridLogMessage<< "*********************************************************" <<std::endl;
+    std::cout << GridLogMessage<< "* Benchmarking DomainWallFermionR::Dhop "<<std::endl;
+    std::cout << GridLogMessage<< "*********************************************************" <<std::endl;
+
+    GridParallelRNG RNG5(FGrid);
+    LatticeFermion src(FGrid); random(RNG5,src);
+    LatticeFermion result(FGrid);
+
+    DomainWallFermionR Dw(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
+    double t0,t1;
+
+    LatticeFermion r_eo(FGrid);
+    LatticeFermion src_e (FrbGrid);
+    LatticeFermion src_o (FrbGrid);
+    LatticeFermion r_e   (FrbGrid);
+    LatticeFermion r_o   (FrbGrid);
+    
+    pickCheckerboard(Even,src_e,src);
+    pickCheckerboard(Odd,src_o,src);
+    
+    setCheckerboard(r_eo,src_o);
+    setCheckerboard(r_eo,src_e);
+    
+    r_e = zero;
+    r_o = zero;
+
+
+#define BENCH_DW(A,in,out)			\
+    Dw.CayleyZeroCounters();			\
+    Dw. A (in,out);				\
+    FGrid->Barrier();				\
+    t0=usecond();				\
+    for(int i=0;i<ncall;i++){			\
+      Dw. A (in,out);				\
+    }						\
+    t1=usecond();				\
+    FGrid->Barrier();				\
+    Dw.CayleyReport();					\
+    std::cout<<GridLogMessage << "Called " #A " "<< (t1-t0)/ncall<<" us"<<std::endl;\
+    std::cout<<GridLogMessage << "******************"<<std::endl;
+
+#define BENCH_DW_MEO(A,in,out)			\
+    Dw.CayleyZeroCounters();			\
+    Dw. A (in,out,0);				\
+    FGrid->Barrier();				\
+    t0=usecond();				\
+    for(int i=0;i<ncall;i++){			\
+      Dw. A (in,out,0);				\
+    }						\
+    t1=usecond();				\
+    FGrid->Barrier();				\
+    Dw.CayleyReport();					\
+    std::cout<<GridLogMessage << "Called " #A " "<< (t1-t0)/ncall<<" us"<<std::endl;\
+    std::cout<<GridLogMessage << "******************"<<std::endl;
+
+    BENCH_DW_MEO(Dhop    ,src,result);
+    BENCH_DW_MEO(DhopEO  ,src_o,r_e);
+    BENCH_DW(Meooe   ,src_o,r_e);
+    BENCH_DW(Mooee   ,src_o,r_o);
+    BENCH_DW(MooeeInv,src_o,r_o);
+
+  }
+
+  if (1)
+  {
+    const int ncall=1000;
+
+    std::cout << GridLogMessage<< "*********************************************************" <<std::endl;
+    std::cout << GridLogMessage<< "* Benchmarking DomainWallFermionVec5dR::Dhop "<<std::endl;
+    std::cout << GridLogMessage<< "*********************************************************" <<std::endl;
+
+    GridParallelRNG RNG5(sFGrid);
+    LatticeFermion src(sFGrid); random(RNG5,src);
+    LatticeFermion sref(sFGrid);
+    LatticeFermion result(sFGrid);
+
+    std::cout<<GridLogMessage << "Constructing Vec5D Dw "<<std::endl;
+    DomainWallFermionVec5dR Dw(Umu,*sFGrid,*sFrbGrid,*sUGrid,*sUrbGrid,mass,M5);
+
+    std::cout<<GridLogMessage << "Calling Dhop "<<std::endl;
+    FGrid->Barrier();
+
+    double t0,t1;
+
+    LatticeFermion r_eo(sFGrid);
+    LatticeFermion src_e (sFrbGrid);
+    LatticeFermion src_o (sFrbGrid);
+    LatticeFermion r_e   (sFrbGrid);
+    LatticeFermion r_o   (sFrbGrid);
+    
+    pickCheckerboard(Even,src_e,src);
+    pickCheckerboard(Odd,src_o,src);
+    
+    setCheckerboard(r_eo,src_o);
+    setCheckerboard(r_eo,src_e);
+    
+    r_e = zero;
+    r_o = zero;
+
+    BENCH_DW_MEO(Dhop    ,src,result);
+    BENCH_DW_MEO(DhopEO  ,src_o,r_e);
+    BENCH_DW(Meooe   ,src_o,r_e);
+    BENCH_DW(Mooee   ,src_o,r_o);
+    BENCH_DW(MooeeInv,src_o,r_o);
+
+  }
+
+  Grid_finalize();
+}

benchmarks/Makefile.am

@@ -1 +1,11 @@
 include Make.inc
+
+simple: simple_su3_test.o simple_su3_expr.o simple_simd_test.o
+
+EXTRA_LIBRARIES = libsimple_su3_test.a libsimple_su3_expr.a libsimple_simd_test.a
+
+libsimple_su3_test_a_SOURCES = simple_su3_test.cc
+
+libsimple_su3_expr_a_SOURCES = simple_su3_expr.cc
+
+libsimple_simd_test_a_SOURCES = simple_simd_test.cc

benchmarks/simple_simd_test.cc

@@ -0,0 +1,11 @@
+#include <Grid/Grid.h>
+
+Grid::vRealD add(const Grid::vRealD &x, const Grid::vRealD &y)
+{
+	return x+y;
+}
+
+Grid::vRealD sub(const Grid::vRealD &x, const Grid::vRealD &y)
+{
+	return x-y;
+}

benchmarks/simple_su3_expr.cc

@@ -25,7 +25,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 using namespace std;
 using namespace Grid;

benchmarks/simple_su3_test.cc

@@ -25,7 +25,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 using namespace std;
 using namespace Grid;

configure.ac

@@ -149,8 +149,14 @@ CXXFLAGS=$CXXFLAGS_CPY
 LDFLAGS=$LDFLAGS_CPY
 
 ############### SIMD instruction selection
-AC_ARG_ENABLE([simd],[AC_HELP_STRING([--enable-simd=<code>],
+AC_ARG_ENABLE([simd],[AC_HELP_STRING([--enable-simd=code],
-	            [select SIMD target (cf. README.md)])], [ac_SIMD=${enable_simd}], [ac_SIMD=GEN])
+              [select SIMD target (cf. README.md)])], [ac_SIMD=${enable_simd}], [ac_SIMD=GEN])
+
+AC_ARG_ENABLE([gen-simd-width],
+            [AS_HELP_STRING([--enable-gen-simd-width=size],
+            [size (in bytes) of the generic SIMD vectors (default: 32)])],
+            [ac_gen_simd_width=$enable_gen_simd_width],
+            [ac_gen_simd_width=32])
 
 case ${ax_cv_cxx_compiler_vendor} in
   clang|gnu)
@@ -180,7 +186,10 @@ case ${ax_cv_cxx_compiler_vendor} in
         AC_DEFINE([AVX512],[1],[AVX512 intrinsics])
         SIMD_FLAGS='-march=knl';;
       GEN)
-        AC_DEFINE([GENERIC_VEC],[1],[generic vector code])
+        AC_DEFINE([GEN],[1],[generic vector code])
+        AC_DEFINE_UNQUOTED([GEN_SIMD_WIDTH],[$ac_gen_simd_width],
+                           [generic SIMD vector width (in bytes)])
+        SIMD_GEN_WIDTH_MSG=" (width= $ac_gen_simd_width)"
         SIMD_FLAGS='';;
       QPX|BGQ)
         AC_DEFINE([QPX],[1],[QPX intrinsics for BG/Q])
@@ -197,8 +206,8 @@ case ${ax_cv_cxx_compiler_vendor} in
         AC_DEFINE([AVX1],[1],[AVX intrinsics])
         SIMD_FLAGS='-mavx -xavx';;
       AVXFMA)
-        AC_DEFINE([AVXFMA],[1],[AVX intrinsics with FMA4])
+        AC_DEFINE([AVXFMA],[1],[AVX intrinsics with FMA3])
-        SIMD_FLAGS='-mavx -mfma';;
+        SIMD_FLAGS='-mavx -fma';;
       AVX2)
         AC_DEFINE([AVX2],[1],[AVX2 intrinsics])
         SIMD_FLAGS='-march=core-avx2 -xcore-avx2';;
@@ -212,7 +221,10 @@ case ${ax_cv_cxx_compiler_vendor} in
         AC_DEFINE([AVX512],[1],[AVX512 intrinsics for Knights Landing])
         SIMD_FLAGS='-xmic-avx512';;
       GEN)
-        AC_DEFINE([GENERIC_VEC],[1],[generic vector code])
+        AC_DEFINE([GEN],[1],[generic vector code])
+        AC_DEFINE_UNQUOTED([GEN_SIMD_WIDTH],[$ac_gen_simd_width],
+                           [generic SIMD vector width (in bytes)])
+        SIMD_GEN_WIDTH_MSG=" (width= $ac_gen_simd_width)"
         SIMD_FLAGS='';;
       *)
         AC_MSG_ERROR(["SIMD option ${ac_SIMD} not supported by the Intel compiler"]);;
@@ -278,7 +290,7 @@ esac
 case ${ac_COMMS} in
     *-auto)
         LX_FIND_MPI
-        if test "x$have_CXX_mpi" = 'xno'; then AC_MSG_ERROR(["MPI not found"]); fi
+        if test "x$have_CXX_mpi" = 'xno'; then AC_MSG_ERROR(["The configure could not find the MPI compilation flags. N.B. The -auto mode is not supported by Cray wrappers. Use the non -auto version in this case."]); fi
         AM_CXXFLAGS="$MPI_CXXFLAGS $AM_CXXFLAGS"
         AM_CFLAGS="$MPI_CFLAGS $AM_CFLAGS"
         AM_LDFLAGS="`echo $MPI_CXXLDFLAGS | sed -E 's/-l@<:@^ @:>@+//g'` $AM_LDFLAGS"
@@ -342,12 +354,17 @@ esac
 AM_CONDITIONAL(BUILD_CHROMA_REGRESSION,[ test "X${ac_CHROMA}X" == "XyesX" ])
 
 ############### Doxygen
-AC_PROG_DOXYGEN
+DX_DOXYGEN_FEATURE([OFF])
-
+DX_DOT_FEATURE([OFF])
-if test -n "$DOXYGEN"
+DX_HTML_FEATURE([ON])
-then
+DX_CHM_FEATURE([OFF])
-AC_CONFIG_FILES([docs/doxy.cfg])
+DX_CHI_FEATURE([OFF])
-fi
+DX_MAN_FEATURE([OFF])
+DX_RTF_FEATURE([OFF])
+DX_XML_FEATURE([OFF])
+DX_PDF_FEATURE([OFF])
+DX_PS_FEATURE([OFF])
+DX_INIT_DOXYGEN([$PACKAGE_NAME], [doxygen.cfg])
 
 ############### Ouput
 cwd=`pwd -P`; cd ${srcdir}; abs_srcdir=`pwd -P`; cd ${cwd}
@@ -364,12 +381,13 @@ AC_CONFIG_FILES(tests/IO/Makefile)
 AC_CONFIG_FILES(tests/core/Makefile)
 AC_CONFIG_FILES(tests/debug/Makefile)
 AC_CONFIG_FILES(tests/forces/Makefile)
+AC_CONFIG_FILES(tests/hadrons/Makefile)
 AC_CONFIG_FILES(tests/hmc/Makefile)
 AC_CONFIG_FILES(tests/solver/Makefile)
 AC_CONFIG_FILES(tests/qdpxx/Makefile)
 AC_CONFIG_FILES(benchmarks/Makefile)
-AC_CONFIG_FILES(programs/Makefile)
+AC_CONFIG_FILES(extras/Makefile)
-AC_CONFIG_FILES(programs/Hadrons/Makefile)
+AC_CONFIG_FILES(extras/Hadrons/Makefile)
 AC_OUTPUT
 
 echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -384,7 +402,7 @@ os (target)                 : $target_os
 compiler vendor             : ${ax_cv_cxx_compiler_vendor}
 compiler version            : ${ax_cv_gxx_version}
 ----- BUILD OPTIONS -----------------------------------
-SIMD                        : ${ac_SIMD}
+SIMD                        : ${ac_SIMD}${SIMD_GEN_WIDTH_MSG}
 Threading                   : ${ac_openmp} 
 Communications type         : ${comms_type}
 Default precision           : ${ac_PRECISION}
@@ -392,8 +410,7 @@ RNG choice                  : ${ac_RNG}
 GMP                         : `if test "x$have_gmp" = xtrue; then echo yes; else echo no; fi`
 LAPACK                      : ${ac_LAPACK}
 FFTW                        : `if test "x$have_fftw" = xtrue; then echo yes; else echo no; fi`
-build DOXYGEN documentation : `if test "x$enable_doc" = xyes; then echo yes; else echo no; fi`
+build DOXYGEN documentation : `if test "$DX_FLAG_doc" = '1'; then echo yes; else echo no; fi`
-graphs and diagrams         : `if test "x$enable_dot" = xyes; then echo yes; else echo no; fi`
 ----- BUILD FLAGS -------------------------------------
 CXXFLAGS:
 `echo ${AM_CXXFLAGS} ${CXXFLAGS} | tr ' ' '\n' | sed 's/^-/    -/g'`

doxygen.inc

@@ -0,0 +1,184 @@
+# Copyright (C) 2004 Oren Ben-Kiki
+# This file is distributed under the same terms as the Automake macro files.
+
+# Generate automatic documentation using Doxygen. Goals and variables values
+# are controlled by the various DX_COND_??? conditionals set by autoconf.
+#
+# The provided goals are:
+# doxygen-doc: Generate all doxygen documentation.
+# doxygen-run: Run doxygen, which will generate some of the documentation
+#              (HTML, CHM, CHI, MAN, RTF, XML) but will not do the post
+#              processing required for the rest of it (PS, PDF, and some MAN).
+# doxygen-man: Rename some doxygen generated man pages.
+# doxygen-ps: Generate doxygen PostScript documentation.
+# doxygen-pdf: Generate doxygen PDF documentation.
+#
+# Note that by default these are not integrated into the automake goals. If
+# doxygen is used to generate man pages, you can achieve this integration by
+# setting man3_MANS to the list of man pages generated and then adding the
+# dependency:
+#
+#   $(man3_MANS): doxygen-doc
+#
+# This will cause make to run doxygen and generate all the documentation.
+#
+# The following variable is intended for use in Makefile.am:
+#
+# DX_CLEANFILES = everything to clean.
+#
+# This is usually added to MOSTLYCLEANFILES.
+
+## --------------------------------- ##
+## Format-independent Doxygen rules. ##
+## --------------------------------- ##
+
+if DX_COND_doc
+
+## ------------------------------- ##
+## Rules specific for HTML output. ##
+## ------------------------------- ##
+
+if DX_COND_html
+
+DX_CLEAN_HTML = @DX_DOCDIR@/html
+
+endif DX_COND_html
+
+## ------------------------------ ##
+## Rules specific for CHM output. ##
+## ------------------------------ ##
+
+if DX_COND_chm
+
+DX_CLEAN_CHM = @DX_DOCDIR@/chm
+
+if DX_COND_chi
+
+DX_CLEAN_CHI = @DX_DOCDIR@/@PACKAGE@.chi
+
+endif DX_COND_chi
+
+endif DX_COND_chm
+
+## ------------------------------ ##
+## Rules specific for MAN output. ##
+## ------------------------------ ##
+
+if DX_COND_man
+
+DX_CLEAN_MAN = @DX_DOCDIR@/man
+
+endif DX_COND_man
+
+## ------------------------------ ##
+## Rules specific for RTF output. ##
+## ------------------------------ ##
+
+if DX_COND_rtf
+
+DX_CLEAN_RTF = @DX_DOCDIR@/rtf
+
+endif DX_COND_rtf
+
+## ------------------------------ ##
+## Rules specific for XML output. ##
+## ------------------------------ ##
+
+if DX_COND_xml
+
+DX_CLEAN_XML = @DX_DOCDIR@/xml
+
+endif DX_COND_xml
+
+## ----------------------------- ##
+## Rules specific for PS output. ##
+## ----------------------------- ##
+
+if DX_COND_ps
+
+DX_CLEAN_PS = @DX_DOCDIR@/@PACKAGE@.ps
+
+DX_PS_GOAL = doxygen-ps
+
+doxygen-ps: @DX_DOCDIR@/@PACKAGE@.ps
+
+@DX_DOCDIR@/@PACKAGE@.ps: @DX_DOCDIR@/@PACKAGE@.tag
+	cd @DX_DOCDIR@/latex; \
+	rm -f *.aux *.toc *.idx *.ind *.ilg *.log *.out; \
+	$(DX_LATEX) refman.tex; \
+	$(MAKEINDEX_PATH) refman.idx; \
+	$(DX_LATEX) refman.tex; \
+	countdown=5; \
+	while $(DX_EGREP) 'Rerun (LaTeX|to get cross-references right)' \
+	                  refman.log > /dev/null 2>&1 \
+	   && test $$countdown -gt 0; do \
+	    $(DX_LATEX) refman.tex; \
+	    countdown=`expr $$countdown - 1`; \
+	done; \
+	$(DX_DVIPS) -o ../@PACKAGE@.ps refman.dvi
+
+endif DX_COND_ps
+
+## ------------------------------ ##
+## Rules specific for PDF output. ##
+## ------------------------------ ##
+
+if DX_COND_pdf
+
+DX_CLEAN_PDF = @DX_DOCDIR@/@PACKAGE@.pdf
+
+DX_PDF_GOAL = doxygen-pdf
+
+doxygen-pdf: @DX_DOCDIR@/@PACKAGE@.pdf
+
+@DX_DOCDIR@/@PACKAGE@.pdf: @DX_DOCDIR@/@PACKAGE@.tag
+	cd @DX_DOCDIR@/latex; \
+	rm -f *.aux *.toc *.idx *.ind *.ilg *.log *.out; \
+	$(DX_PDFLATEX) refman.tex; \
+	$(DX_MAKEINDEX) refman.idx; \
+	$(DX_PDFLATEX) refman.tex; \
+	countdown=5; \
+	while $(DX_EGREP) 'Rerun (LaTeX|to get cross-references right)' \
+	                  refman.log > /dev/null 2>&1 \
+	   && test $$countdown -gt 0; do \
+	    $(DX_PDFLATEX) refman.tex; \
+	    countdown=`expr $$countdown - 1`; \
+	done; \
+	mv refman.pdf ../@PACKAGE@.pdf
+
+endif DX_COND_pdf
+
+## ------------------------------------------------- ##
+## Rules specific for LaTeX (shared for PS and PDF). ##
+## ------------------------------------------------- ##
+
+if DX_COND_latex
+
+DX_CLEAN_LATEX = @DX_DOCDIR@/latex
+
+endif DX_COND_latex
+
+.INTERMEDIATE: doxygen-run $(DX_PS_GOAL) $(DX_PDF_GOAL)
+
+doxygen-run: @DX_DOCDIR@/@PACKAGE@.tag
+
+doxygen-doc: doxygen-run $(DX_PS_GOAL) $(DX_PDF_GOAL)
+
+@DX_DOCDIR@/@PACKAGE@.tag: $(DX_CONFIG) $(pkginclude_HEADERS)
+	rm -rf @DX_DOCDIR@
+	$(DX_ENV) $(DX_DOXYGEN) $(srcdir)/$(DX_CONFIG)
+
+DX_CLEANFILES = \
+    @DX_DOCDIR@/@PACKAGE@.tag \
+    -r \
+    $(DX_CLEAN_HTML) \
+    $(DX_CLEAN_CHM) \
+    $(DX_CLEAN_CHI) \
+    $(DX_CLEAN_MAN) \
+    $(DX_CLEAN_RTF) \
+    $(DX_CLEAN_XML) \
+    $(DX_CLEAN_PS) \
+    $(DX_CLEAN_PDF) \
+    $(DX_CLEAN_LATEX)
+
+endif DX_COND_doc

extras/Hadrons/Application.cc

@@ -0,0 +1,312 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Application.cc
+
+Copyright (C) 2015
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 <Grid/Hadrons/Application.hpp>
+#include <Grid/Hadrons/GeneticScheduler.hpp>
+
+using namespace Grid;
+using namespace QCD;
+using namespace Hadrons;
+
+#define BIG_SEP "==============="
+#define SEP     "---------------"
+
+/******************************************************************************
+ *                       Application implementation                           *
+ ******************************************************************************/
+// constructors ////////////////////////////////////////////////////////////////
+Application::Application(void)
+: env_(Environment::getInstance())
+{
+    LOG(Message) << "Modules available:" << std::endl;
+    auto list = ModuleFactory::getInstance().getBuilderList();
+    for (auto &m: list)
+    {
+        LOG(Message) << "  " << m << std::endl;
+    }
+    auto dim = GridDefaultLatt(), mpi = GridDefaultMpi(), loc(dim);
+    locVol_ = 1;
+    for (unsigned int d = 0; d < dim.size(); ++d)
+    {
+        loc[d]  /= mpi[d];
+        locVol_ *= loc[d];
+    }
+    LOG(Message) << "Global lattice: " << dim << std::endl;
+    LOG(Message) << "MPI partition : " << mpi << std::endl;
+    LOG(Message) << "Local lattice : " << loc << std::endl;
+}
+
+Application::Application(const Application::GlobalPar &par)
+: Application()
+{
+    setPar(par);
+}
+
+Application::Application(const std::string parameterFileName)
+: Application()
+{
+    parameterFileName_ = parameterFileName;
+}
+
+// access //////////////////////////////////////////////////////////////////////
+void Application::setPar(const Application::GlobalPar &par)
+{
+    par_ = par;
+    env_.setSeed(strToVec<int>(par_.seed));
+}
+
+const Application::GlobalPar & Application::getPar(void)
+{
+    return par_;
+}
+
+// execute /////////////////////////////////////////////////////////////////////
+void Application::run(void)
+{
+    if (!parameterFileName_.empty() and (env_.getNModule() == 0))
+    {
+        parseParameterFile(parameterFileName_);
+    }
+    if (!scheduled_)
+    {
+        schedule();
+    }
+    printSchedule();
+    configLoop();
+}
+
+// parse parameter file ////////////////////////////////////////////////////////
+class ObjectId: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(ObjectId,
+                                    std::string, name,
+                                    std::string, type);
+};
+
+void Application::parseParameterFile(const std::string parameterFileName)
+{
+    XmlReader reader(parameterFileName);
+    GlobalPar par;
+    ObjectId  id;
+    
+    LOG(Message) << "Building application from '" << parameterFileName << "'..." << std::endl;
+    read(reader, "parameters", par);
+    setPar(par);
+    push(reader, "modules");
+    push(reader, "module");
+    do
+    {
+        read(reader, "id", id);
+        env_.createModule(id.name, id.type, reader);
+    } while (reader.nextElement("module"));
+    pop(reader);
+    pop(reader);
+}
+
+void Application::saveParameterFile(const std::string parameterFileName)
+{
+    XmlWriter          writer(parameterFileName);
+    ObjectId           id;
+    const unsigned int nMod = env_.getNModule();
+    
+    LOG(Message) << "Saving application to '" << parameterFileName << "'..." << std::endl;
+    write(writer, "parameters", getPar());
+    push(writer, "modules");
+    for (unsigned int i = 0; i < nMod; ++i)
+    {
+        push(writer, "module");
+        id.name = env_.getModuleName(i);
+        id.type = env_.getModule(i)->getRegisteredName();
+        write(writer, "id", id);
+        env_.getModule(i)->saveParameters(writer, "options");
+        pop(writer);
+    }
+    pop(writer);
+    pop(writer);
+}
+
+// schedule computation ////////////////////////////////////////////////////////
+#define MEM_MSG(size)\
+sizeString((size)*locVol_) << " (" << sizeString(size)  << "/site)"
+
+#define DEFINE_MEMPEAK \
+auto memPeak = [this](const std::vector<unsigned int> &program)\
+{\
+    unsigned int memPeak;\
+    bool         msg;\
+    \
+    msg = HadronsLogMessage.isActive();\
+    HadronsLogMessage.Active(false);\
+    env_.dryRun(true);\
+    memPeak = env_.executeProgram(program);\
+    env_.dryRun(false);\
+    env_.freeAll();\
+    HadronsLogMessage.Active(true);\
+    \
+    return memPeak;\
+}
+
+void Application::schedule(void)
+{
+    DEFINE_MEMPEAK;
+    
+    // build module dependency graph
+    LOG(Message) << "Building module graph..." << std::endl;
+    auto graph = env_.makeModuleGraph();
+    auto con = graph.getConnectedComponents();
+    
+    // constrained topological sort using a genetic algorithm
+    LOG(Message) << "Scheduling computation..." << std::endl;
+    LOG(Message) << "               #module= " << graph.size() << std::endl;
+    LOG(Message) << "       population size= " << par_.genetic.popSize << std::endl;
+    LOG(Message) << "       max. generation= " << par_.genetic.maxGen << std::endl;
+    LOG(Message) << "  max. cst. generation= " << par_.genetic.maxCstGen << std::endl;
+    LOG(Message) << "         mutation rate= " << par_.genetic.mutationRate << std::endl;
+    
+    unsigned int                               k = 0, gen, prevPeak, nCstPeak = 0;
+    std::random_device                         rd;
+    GeneticScheduler<unsigned int>::Parameters par;
+    
+    par.popSize      = par_.genetic.popSize;
+    par.mutationRate = par_.genetic.mutationRate;
+    par.seed         = rd();
+    memPeak_         = 0;
+    CartesianCommunicator::BroadcastWorld(0, &(par.seed), sizeof(par.seed));
+    for (unsigned int i = 0; i < con.size(); ++i)
+    {
+        GeneticScheduler<unsigned int> scheduler(con[i], memPeak, par);
+        
+        gen = 0;
+        do
+        {
+            LOG(Debug) << "Generation " << gen << ":" << std::endl;
+            scheduler.nextGeneration();
+            if (gen != 0)
+            {
+                if (prevPeak == scheduler.getMinValue())
+                {
+                    nCstPeak++;
+                }
+                else
+                {
+                    nCstPeak = 0;
+                }
+            }
+            
+            prevPeak = scheduler.getMinValue();
+            if (gen % 10 == 0)
+            {
+                LOG(Iterative) << "Generation " << gen << ": "
+                               << MEM_MSG(scheduler.getMinValue()) << std::endl;
+            }
+            
+            gen++;
+        } while ((gen < par_.genetic.maxGen)
+                 and (nCstPeak < par_.genetic.maxCstGen));
+        auto &t = scheduler.getMinSchedule();
+        if (scheduler.getMinValue() > memPeak_)
+        {
+            memPeak_ = scheduler.getMinValue();
+        }
+        for (unsigned int j = 0; j < t.size(); ++j)
+        {
+            program_.push_back(t[j]);
+        }
+    }
+    scheduled_ = true;
+}
+
+void Application::saveSchedule(const std::string filename)
+{
+    TextWriter               writer(filename);
+    std::vector<std::string> program;
+    
+    if (!scheduled_)
+    {
+        HADRON_ERROR("Computation not scheduled");
+    }
+    LOG(Message) << "Saving current schedule to '" << filename << "'..."
+                 << std::endl;
+    for (auto address: program_)
+    {
+        program.push_back(env_.getModuleName(address));
+    }
+    write(writer, "schedule", program);
+}
+
+void Application::loadSchedule(const std::string filename)
+{
+    DEFINE_MEMPEAK;
+    
+    TextReader               reader(filename);
+    std::vector<std::string> program;
+    
+    LOG(Message) << "Loading schedule from '" << filename << "'..."
+                 << std::endl;
+    read(reader, "schedule", program);
+    program_.clear();
+    for (auto &name: program)
+    {
+        program_.push_back(env_.getModuleAddress(name));
+    }
+    scheduled_ = true;
+    memPeak_   = memPeak(program_);
+}
+
+void Application::printSchedule(void)
+{
+    if (!scheduled_)
+    {
+        HADRON_ERROR("Computation not scheduled");
+    }
+    LOG(Message) << "Schedule (memory peak: " << MEM_MSG(memPeak_) << "):"
+                 << std::endl;
+    for (unsigned int i = 0; i < program_.size(); ++i)
+    {
+        LOG(Message) << std::setw(4) << i + 1 << ": "
+                     << env_.getModuleName(program_[i]) << std::endl;
+    }
+}
+
+// loop on configurations //////////////////////////////////////////////////////
+void Application::configLoop(void)
+{
+    auto range = par_.trajCounter;
+    
+    for (unsigned int t = range.start; t < range.end; t += range.step)
+    {
+        LOG(Message) << BIG_SEP << " Starting measurement for trajectory " << t
+                     << " " << BIG_SEP << std::endl;
+        env_.setTrajectory(t);
+        env_.executeProgram(program_);
+    }
+    LOG(Message) << BIG_SEP << " End of measurement " << BIG_SEP << std::endl;
+    env_.freeAll();
+}

extras/Hadrons/Application.hpp

@@ -0,0 +1,130 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Application.hpp
+
+Copyright (C) 2015
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 */
+
+#ifndef Hadrons_Application_hpp_
+#define Hadrons_Application_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Environment.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/Modules.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         Main program manager                               *
+ ******************************************************************************/
+class Application
+{
+public:
+    class TrajRange: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(TrajRange,
+                                        unsigned int, start,
+                                        unsigned int, end,
+                                        unsigned int, step);
+    };
+    class GeneticPar: Serializable
+    {
+    public:
+        GeneticPar(void):
+            popSize{20}, maxGen{1000}, maxCstGen{100}, mutationRate{.1} {};
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(GeneticPar,
+                                        unsigned int, popSize,
+                                        unsigned int, maxGen,
+                                        unsigned int, maxCstGen,
+                                        double      , mutationRate);
+    };
+    class GlobalPar: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(GlobalPar,
+                                        TrajRange,   trajCounter,
+                                        GeneticPar,  genetic,
+                                        std::string, seed);
+    };
+public:
+    // constructors
+    Application(void);
+    Application(const GlobalPar &par);
+    Application(const std::string parameterFileName);
+    // destructor
+    virtual ~Application(void) = default;
+    // access
+    void              setPar(const GlobalPar &par);
+    const GlobalPar & getPar(void);
+    // module creation
+    template <typename M>
+    void createModule(const std::string name);
+    template <typename M>
+    void createModule(const std::string name, const typename M::Par &par);
+    // execute
+    void run(void);
+    // XML parameter file I/O
+    void parseParameterFile(const std::string parameterFileName);
+    void saveParameterFile(const std::string parameterFileName);
+    // schedule computation
+    void schedule(void);
+    void saveSchedule(const std::string filename);
+    void loadSchedule(const std::string filename);
+    void printSchedule(void);
+    // loop on configurations
+    void configLoop(void);
+private:
+    long unsigned int         locVol_;
+    std::string               parameterFileName_{""};
+    GlobalPar                 par_;
+    Environment               &env_;
+    std::vector<unsigned int> program_;
+    Environment::Size         memPeak_;
+    bool                      scheduled_{false};
+};
+
+/******************************************************************************
+ *                     Application template implementation                    *
+ ******************************************************************************/
+// module creation /////////////////////////////////////////////////////////////
+template <typename M>
+void Application::createModule(const std::string name)
+{
+    env_.createModule<M>(name);
+}
+
+template <typename M>
+void Application::createModule(const std::string name,
+                               const typename M::Par &par)
+{
+    env_.createModule<M>(name, par);
+}
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Application_hpp_

extras/Hadrons/Environment.cc

@@ -1,9 +1,11 @@
-/*******************************************************************************
+/*************************************************************************************
+
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: programs/Hadrons/Environment.cc
+Source file: extras/Hadrons/Environment.cc
 
 Copyright (C) 2015
+Copyright (C) 2016
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 
@@ -21,13 +23,13 @@ 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 
+See the full license in the file "LICENSE" in the top level distribution directory
-directory.
+*************************************************************************************/
-*******************************************************************************/
+/*  END LEGAL */
 
-#include <Hadrons/Environment.hpp>
+#include <Grid/Hadrons/Environment.hpp>
-#include <Hadrons/Module.hpp>
+#include <Grid/Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
 
 using namespace Grid;
 using namespace QCD;
@@ -136,19 +138,19 @@ GridParallelRNG * Environment::get4dRng(void) const
 }
 
 // module management ///////////////////////////////////////////////////////////
-void Environment::createModule(const std::string name, const std::string type,
+void Environment::pushModule(Environment::ModPt &pt)
-                               XmlReader &reader)
 {
+    std::string name = pt->getName();
+    
     if (!hasModule(name))
     {
-        auto                      &factory = ModuleFactory::getInstance();
         std::vector<unsigned int> inputAddress;
+        unsigned int              address;
         ModuleInfo                m;
         
-        m.data = factory.create(type, name);
+        m.data = std::move(pt);
         m.type = typeIdPt(*m.data.get());
         m.name = name;
-        m.data->parseParameters(reader, "options");
         auto input  = m.data->getInput();
         for (auto &in: input)
         {
@@ -160,25 +162,26 @@ void Environment::createModule(const std::string name, const std::string type,
         }
         auto output = m.data->getOutput();
         module_.push_back(std::move(m));
-        moduleAddress_[name] = module_.size() - 1;
+        address              = static_cast<unsigned int>(module_.size() - 1);
+        moduleAddress_[name] = address;
         for (auto &out: output)
         {
             if (!hasObject(out))
             {
-                addObject(out , module_.size() - 1);
+                addObject(out, address);
             }
             else
             {
                 if (object_[objectAddress_[out]].module < 0)
                 {
-                    object_[objectAddress_[out]].module = module_.size() - 1;
+                    object_[objectAddress_[out]].module = address;
                 }
                 else
                 {
                     HADRON_ERROR("object '" + out
-                        + "' is already produced by module '"
+                                 + "' is already produced by module '"
-                        + module_[object_[getObjectAddress(out)].module].name
+                                 + module_[object_[getObjectAddress(out)].module].name
-                        + "' (while creating module '" + name + "')");
+                                 + "' (while pushing module '" + name + "')");
                 }
             }
         }
@@ -189,6 +192,21 @@ void Environment::createModule(const std::string name, const std::string type,
     }
 }
 
+unsigned int Environment::getNModule(void) const
+{
+    return module_.size();
+}
+
+void Environment::createModule(const std::string name, const std::string type,
+                               XmlReader &reader)
+{
+    auto &factory = ModuleFactory::getInstance();
+    auto pt       = factory.create(type, name);
+    
+    pt->parseParameters(reader, "options");
+    pushModule(pt);
+}
+
 ModuleBase * Environment::getModule(const unsigned int address) const
 {
     if (hasModule(address))
@@ -234,7 +252,7 @@ std::string Environment::getModuleType(const unsigned int address) const
 {
     if (hasModule(address))
     {
-        return module_[address].type->name();
+        return typeName(module_[address].type);
     }
     else
     {
@@ -278,9 +296,10 @@ Graph<unsigned int> Environment::makeModuleGraph(void) const
 #define MEM_MSG(size)\
 sizeString((size)*locVol_) << " (" << sizeString(size)  << "/site)"
 
-unsigned int Environment::executeProgram(const std::vector<unsigned int> &p)
+Environment::Size
+Environment::executeProgram(const std::vector<unsigned int> &p)
 {
-    unsigned int                        memPeak = 0, sizeBefore, sizeAfter;
+    Size                                memPeak = 0, sizeBefore, sizeAfter;
     std::vector<std::set<unsigned int>> freeProg;
     bool                                continueCollect, nothingFreed;
     
@@ -376,7 +395,7 @@ unsigned int Environment::executeProgram(const std::vector<unsigned int> &p)
     return memPeak;
 }
 
-unsigned int Environment::executeProgram(const std::vector<std::string> &p)
+Environment::Size Environment::executeProgram(const std::vector<std::string> &p)
 {
     std::vector<unsigned int> pAddress;
     
@@ -391,12 +410,19 @@ unsigned int Environment::executeProgram(const std::vector<std::string> &p)
 // general memory management ///////////////////////////////////////////////////
 void Environment::addObject(const std::string name, const int moduleAddress)
 {
-    ObjInfo info;
+    if (!hasObject(name))
+    {
+        ObjInfo info;
         
-    info.name   = name;
+        info.name   = name;
-    info.module = moduleAddress;
+        info.module = moduleAddress;
-    object_.push_back(std::move(info));
+        object_.push_back(std::move(info));
-    objectAddress_[name] = object_.size() - 1;
+        objectAddress_[name] = static_cast<unsigned int>(object_.size() - 1);
+    }
+    else
+    {
+        HADRON_ERROR("object '" + name + "' already exists");
+    }
 }
 
 void Environment::registerObject(const unsigned int address,
@@ -425,6 +451,10 @@ void Environment::registerObject(const unsigned int address,
 void Environment::registerObject(const std::string name,
                                  const unsigned int size, const unsigned int Ls)
 {
+    if (!hasObject(name))
+    {
+        addObject(name);
+    }
     registerObject(getObjectAddress(name), size, Ls);
 }
 
@@ -456,7 +486,7 @@ std::string Environment::getObjectType(const unsigned int address) const
 {
     if (hasRegisteredObject(address))
     {
-        return object_[address].type->name();
+        return typeName(object_[address].type);
     }
     else if (hasObject(address))
     {
@@ -474,7 +504,7 @@ std::string Environment::getObjectType(const std::string name) const
     return getObjectType(getObjectAddress(name));
 }
 
-unsigned int Environment::getObjectSize(const unsigned int address) const
+Environment::Size Environment::getObjectSize(const unsigned int address) const
 {
     if (hasRegisteredObject(address))
     {
@@ -491,7 +521,7 @@ unsigned int Environment::getObjectSize(const unsigned int address) const
     }
 }
 
-unsigned int Environment::getObjectSize(const std::string name) const
+Environment::Size Environment::getObjectSize(const std::string name) const
 {
     return getObjectSize(getObjectAddress(name));
 }
@@ -554,6 +584,30 @@ bool Environment::hasRegisteredObject(const std::string name) const
     }
 }
 
+bool Environment::hasCreatedObject(const unsigned int address) const
+{
+    if (hasObject(address))
+    {
+        return (object_[address].data != nullptr);
+    }
+    else
+    {
+        return false;
+    }
+}
+
+bool Environment::hasCreatedObject(const std::string name) const
+{
+    if (hasObject(name))
+    {
+        return hasCreatedObject(getObjectAddress(name));
+    }
+    else
+    {
+        return false;
+    }
+}
+
 bool Environment::isObject5d(const unsigned int address) const
 {
     return (getObjectLs(address) > 1);
@@ -564,9 +618,9 @@ bool Environment::isObject5d(const std::string name) const
     return (getObjectLs(name) > 1);
 }
 
-long unsigned int Environment::getTotalSize(void) const
+Environment::Size Environment::getTotalSize(void) const
 {
-    long unsigned int size = 0;
+    Environment::Size size = 0;
     
     for (auto &o: object_)
     {
@@ -628,7 +682,7 @@ bool Environment::freeObject(const unsigned int address)
 {
     if (!hasOwners(address))
     {
-        if (!isDryRun())
+        if (!isDryRun() and object_[address].isRegistered)
         {
             LOG(Message) << "Destroying object '" << object_[address].name
                          << "'" << std::endl;
@@ -671,14 +725,13 @@ void Environment::printContent(void)
     LOG(Message) << "Modules: " << std::endl;
     for (unsigned int i = 0; i < module_.size(); ++i)
     {
-        LOG(Message) << std::setw(4) << std::right << i << ": "
+        LOG(Message) << std::setw(4) << i << ": "
-                     << getModuleName(i) << " ("
+                     << getModuleName(i) << std::endl;
-                     << getModuleType(i) << ")" << std::endl;
     }
     LOG(Message) << "Objects: " << std::endl;
     for (unsigned int i = 0; i < object_.size(); ++i)
     {
-        LOG(Message) << std::setw(4) << std::right << i << ": "
+        LOG(Message) << std::setw(4) << i << ": "
                      << getObjectName(i) << std::endl;
     }
 }

extras/Hadrons/Environment.hpp

@@ -1,9 +1,11 @@
-/*******************************************************************************
+/*************************************************************************************
+
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: programs/Hadrons/Environment.hpp
+Source file: extras/Hadrons/Environment.hpp
 
 Copyright (C) 2015
+Copyright (C) 2016
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 
@@ -21,15 +23,19 @@ 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 
+See the full license in the file "LICENSE" in the top level distribution directory
-directory.
+*************************************************************************************/
-*******************************************************************************/
+/*  END LEGAL */
 
 #ifndef Hadrons_Environment_hpp_
 #define Hadrons_Environment_hpp_
 
-#include <Hadrons/Global.hpp>
+#include <Grid/Hadrons/Global.hpp>
-#include <Hadrons/Graph.hpp>
+#include <Grid/Hadrons/Graph.hpp>
+
+#ifndef SITE_SIZE_TYPE
+#define SITE_SIZE_TYPE unsigned int
+#endif
 
 BEGIN_HADRONS_NAMESPACE
 
@@ -64,26 +70,24 @@ class Environment
 {
     SINGLETON(Environment);
 public:
-    typedef std::unique_ptr<ModuleBase>                 ModPt;
+    typedef SITE_SIZE_TYPE                         Size;
-    typedef std::unique_ptr<GridCartesian>              GridPt;
+    typedef std::unique_ptr<ModuleBase>            ModPt;
-    typedef std::unique_ptr<GridRedBlackCartesian>      GridRbPt;
+    typedef std::unique_ptr<GridCartesian>         GridPt;
-    typedef FermionOperator<WilsonImplR>                FMat;
+    typedef std::unique_ptr<GridRedBlackCartesian> GridRbPt;
-    typedef std::unique_ptr<FMat>                       FMatPt;
+    typedef std::unique_ptr<GridParallelRNG>       RngPt;
-    typedef std::function<void(LatticeFermion &,
+    typedef std::unique_ptr<LatticeBase>           LatticePt;
-                               const LatticeFermion &)> Solver;
-    typedef std::unique_ptr<GridParallelRNG>            RngPt;
-    typedef std::unique_ptr<LatticeBase>                LatticePt;
 private:
     struct ModuleInfo
     {
-        const std::type_info        *type{nullptr};
+        const std::type_info      *type{nullptr};
-        std::string                 name;
+        std::string               name;
-        std::unique_ptr<ModuleBase> data{nullptr};
+        ModPt                     data{nullptr};
-        std::vector<unsigned int>   input;
+        std::vector<unsigned int> input;
     };
     struct ObjInfo
     {
-        unsigned int            size{0}, Ls{0};
+        Size                    size{0};
+        unsigned int            Ls{0};
         bool                    isRegistered{false};
         const std::type_info    *type{nullptr};
         std::string             name;
@@ -106,9 +110,16 @@ public:
     void                    setSeed(const std::vector<int> &seed);
     GridParallelRNG *       get4dRng(void) const;
     // module management
+    void                    pushModule(ModPt &pt);
+    template <typename M>
+    void                    createModule(const std::string name);
+    template <typename M>
+    void                    createModule(const std::string name,
+                                         const typename M::Par &par);
     void                    createModule(const std::string name,
                                          const std::string type,
                                          XmlReader &reader);
+    unsigned int            getNModule(void) const;
     ModuleBase *            getModule(const unsigned int address) const;
     ModuleBase *            getModule(const std::string name) const;
     template <typename M>
@@ -122,11 +133,11 @@ public:
     bool                    hasModule(const unsigned int address) const;
     bool                    hasModule(const std::string name) const;
     Graph<unsigned int>     makeModuleGraph(void) const;
-    unsigned int            executeProgram(const std::vector<unsigned int> &p);
+    Size                    executeProgram(const std::vector<unsigned int> &p);
-    unsigned int            executeProgram(const std::vector<std::string> &p);
+    Size                    executeProgram(const std::vector<std::string> &p);
     // general memory management
     void                    addObject(const std::string name,
-                                      const int moduleAddress);
+                                      const int moduleAddress = -1);
     void                    registerObject(const unsigned int address,
                                            const unsigned int size,
                                            const unsigned int Ls = 1);
@@ -157,17 +168,19 @@ public:
     std::string             getObjectName(const unsigned int address) const;
     std::string             getObjectType(const unsigned int address) const;
     std::string             getObjectType(const std::string name) const;
-    unsigned int            getObjectSize(const unsigned int address) const;
+    Size                    getObjectSize(const unsigned int address) const;
-    unsigned int            getObjectSize(const std::string name) const;
+    Size                    getObjectSize(const std::string name) const;
     unsigned int            getObjectLs(const unsigned int address) const;
     unsigned int            getObjectLs(const std::string name) const;
     bool                    hasObject(const unsigned int address) const;
     bool                    hasObject(const std::string name) const;
     bool                    hasRegisteredObject(const unsigned int address) const;
     bool                    hasRegisteredObject(const std::string name) const;
+    bool                    hasCreatedObject(const unsigned int address) const;
+    bool                    hasCreatedObject(const std::string name) const;
     bool                    isObject5d(const unsigned int address) const;
     bool                    isObject5d(const std::string name) const;
-    long unsigned int       getTotalSize(void) const;
+    Environment::Size       getTotalSize(void) const;
     void                    addOwnership(const unsigned int owner,
                                          const unsigned int property);
     void                    addOwnership(const std::string owner,
@@ -194,24 +207,21 @@ private:
     std::map<std::string, unsigned int>    moduleAddress_;
     // lattice store
     std::map<unsigned int, LatticePt>      lattice_;
-    // fermion matrix store
-    std::map<unsigned int, FMatPt>         fMat_;
-    // solver store & solver/action map
-    std::map<unsigned int, Solver>         solver_;
-    std::map<std::string, std::string>     solverAction_;
     // object store
     std::vector<ObjInfo>                   object_;
     std::map<std::string, unsigned int>    objectAddress_;
 };
 
 /******************************************************************************
- *                        template implementation                             *
+ *                       Holder template implementation                       *
  ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
 template <typename T>
 Holder<T>::Holder(T *pt)
 : objPt_(pt)
 {}
 
+// access //////////////////////////////////////////////////////////////////////
 template <typename T>
 T & Holder<T>::get(void) const
 {
@@ -230,6 +240,28 @@ void Holder<T>::reset(T *pt)
     objPt_.reset(pt);
 }
 
+/******************************************************************************
+ *                     Environment template implementation                    *
+ ******************************************************************************/
+// module management ///////////////////////////////////////////////////////////
+template <typename M>
+void Environment::createModule(const std::string name)
+{
+    ModPt pt(new M(name));
+    
+    pushModule(pt);
+}
+
+template <typename M>
+void Environment::createModule(const std::string name,
+                               const typename M::Par &par)
+{
+    ModPt pt(new M(name));
+    
+    static_cast<M *>(pt.get())->setPar(par);
+    pushModule(pt);
+}
+
 template <typename M>
 M * Environment::getModule(const unsigned int address) const
 {

extras/Hadrons/Factory.hpp

@@ -1,8 +1,10 @@
-/*******************************************************************************
+/*************************************************************************************
+
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: programs/Hadrons/Factory.hpp
+Source file: extras/Hadrons/Factory.hpp
 
+Copyright (C) 2015
 Copyright (C) 2016
 
 Author: Antonin Portelli <antonin.portelli@me.com>
@@ -21,14 +23,14 @@ 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 
+See the full license in the file "LICENSE" in the top level distribution directory
-directory.
+*************************************************************************************/
-*******************************************************************************/
+/*  END LEGAL */
 
 #ifndef Hadrons_Factory_hpp_
 #define Hadrons_Factory_hpp_
 
-#include <Hadrons/Global.hpp>
+#include <Grid/Hadrons/Global.hpp>
 
 BEGIN_HADRONS_NAMESPACE
 

extras/Hadrons/GeneticScheduler.hpp

@@ -0,0 +1,329 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/GeneticScheduler.hpp
+
+Copyright (C) 2015
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 */
+
+#ifndef Hadrons_GeneticScheduler_hpp_
+#define Hadrons_GeneticScheduler_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Graph.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                   Scheduler based on a genetic algorithm                   *
+ ******************************************************************************/
+template <typename T>
+class GeneticScheduler
+{
+public:
+    typedef std::vector<T>                   Gene;
+    typedef std::pair<Gene *, Gene *>        GenePair;
+    typedef std::function<int(const Gene &)> ObjFunc;
+    struct Parameters
+    {
+        double       mutationRate;
+        unsigned int popSize, seed;
+    };
+public:
+    // constructor
+    GeneticScheduler(Graph<T> &graph, const ObjFunc &func,
+                     const Parameters &par);
+    // destructor
+    virtual ~GeneticScheduler(void) = default;
+    // access
+    const Gene & getMinSchedule(void);
+    int          getMinValue(void);
+    // breed a new generation
+    void nextGeneration(void);
+    // heuristic benchmarks
+    void benchmarkCrossover(const unsigned int nIt);
+    // print population
+    friend std::ostream & operator<<(std::ostream &out,
+                                     const GeneticScheduler<T> &s)
+    {
+        out << "[";
+        for (auto &p: s.population_)
+        {
+            out << p.first << ", ";
+        }
+        out << "\b\b]";
+        
+        return out;
+    }
+private:
+    // evolution steps
+    void initPopulation(void);
+    void doCrossover(void);
+    void doMutation(void);
+    // genetic operators
+    GenePair selectPair(void);
+    void     crossover(Gene &c1, Gene &c2, const Gene &p1, const Gene &p2);
+    void     mutation(Gene &m, const Gene &c);
+    
+private:
+    Graph<T>                 &graph_;
+    const ObjFunc            &func_;
+    const Parameters         par_;
+    std::multimap<int, Gene> population_;
+    std::mt19937             gen_;
+};
+
+/******************************************************************************
+ *                       template implementation                              *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename T>
+GeneticScheduler<T>::GeneticScheduler(Graph<T> &graph, const ObjFunc &func,
+                                      const Parameters &par)
+: graph_(graph)
+, func_(func)
+, par_(par)
+{
+    gen_.seed(par_.seed);
+}
+
+// access //////////////////////////////////////////////////////////////////////
+template <typename T>
+const typename GeneticScheduler<T>::Gene &
+GeneticScheduler<T>::getMinSchedule(void)
+{
+    return population_.begin()->second;
+}
+
+template <typename T>
+int GeneticScheduler<T>::getMinValue(void)
+{
+    return population_.begin()->first;
+}
+
+// breed a new generation //////////////////////////////////////////////////////
+template <typename T>
+void GeneticScheduler<T>::nextGeneration(void)
+{
+    // random initialization of the population if necessary
+    if (population_.size() != par_.popSize)
+    {
+        initPopulation();
+    }
+    LOG(Debug) << "Starting population:\n" << *this << std::endl;
+    
+    // random mutations
+    //PARALLEL_FOR_LOOP
+    for (unsigned int i = 0; i < par_.popSize; ++i)
+    {
+        doMutation();
+    }
+    LOG(Debug) << "After mutations:\n" << *this << std::endl;
+    
+    // mating
+    //PARALLEL_FOR_LOOP
+    for (unsigned int i = 0; i < par_.popSize/2; ++i)
+    {
+        doCrossover();
+    }
+    LOG(Debug) << "After mating:\n" << *this << std::endl;
+    
+    // grim reaper
+    auto it = population_.begin();
+    
+    std::advance(it, par_.popSize);
+    population_.erase(it, population_.end());
+    LOG(Debug) << "After grim reaper:\n" << *this << std::endl;
+}
+
+// evolution steps /////////////////////////////////////////////////////////////
+template <typename T>
+void GeneticScheduler<T>::initPopulation(void)
+{
+    population_.clear();
+    for (unsigned int i = 0; i < par_.popSize; ++i)
+    {
+        auto p = graph_.topoSort(gen_);
+        
+        population_.emplace(func_(p), p);
+    }
+}
+
+template <typename T>
+void GeneticScheduler<T>::doCrossover(void)
+{
+    auto p = selectPair();
+    Gene &p1 = *(p.first), &p2 = *(p.second);
+    Gene c1, c2;
+    
+    crossover(c1, c2, p1, p2);
+    PARALLEL_CRITICAL
+    {
+        population_.emplace(func_(c1), c1);
+        population_.emplace(func_(c2), c2);
+    }
+}
+
+template <typename T>
+void GeneticScheduler<T>::doMutation(void)
+{
+    std::uniform_real_distribution<double>      mdis(0., 1.);
+    std::uniform_int_distribution<unsigned int> pdis(0, population_.size() - 1);
+    
+    if (mdis(gen_) < par_.mutationRate)
+    {
+        Gene m;
+        auto it = population_.begin();
+        
+        std::advance(it, pdis(gen_));
+        mutation(m, it->second);
+        PARALLEL_CRITICAL
+        {
+            population_.emplace(func_(m), m);
+        }
+    }
+}
+
+// genetic operators ///////////////////////////////////////////////////////////
+template <typename T>
+typename GeneticScheduler<T>::GenePair GeneticScheduler<T>::selectPair(void)
+{
+    std::vector<double> prob;
+    unsigned int        ind;
+    Gene                *p1, *p2;
+    
+    for (auto &c: population_)
+    {
+        prob.push_back(1./c.first);
+    }
+    do
+    {
+        double probCpy;
+        
+        std::discrete_distribution<unsigned int> dis1(prob.begin(), prob.end());
+        auto rIt = population_.begin();
+        ind = dis1(gen_);
+        std::advance(rIt, ind);
+        p1 = &(rIt->second);
+        probCpy   = prob[ind];
+        prob[ind] = 0.;
+        std::discrete_distribution<unsigned int> dis2(prob.begin(), prob.end());
+        rIt = population_.begin();
+        std::advance(rIt, dis2(gen_));
+        p2 = &(rIt->second);
+        prob[ind] = probCpy;
+    } while (p1 == p2);
+    
+    return std::make_pair(p1, p2);
+}
+
+template <typename T>
+void GeneticScheduler<T>::crossover(Gene &c1, Gene &c2, const Gene &p1,
+                                    const Gene &p2)
+{
+    Gene                                        buf;
+    std::uniform_int_distribution<unsigned int> dis(0, p1.size() - 1);
+    unsigned int                                cut = dis(gen_);
+    
+    c1.clear();
+    buf = p2;
+    for (unsigned int i = 0; i < cut; ++i)
+    {
+        c1.push_back(p1[i]);
+        buf.erase(std::find(buf.begin(), buf.end(), p1[i]));
+    }
+    for (unsigned int i = 0; i < buf.size(); ++i)
+    {
+        c1.push_back(buf[i]);
+    }
+    c2.clear();
+    buf = p2;
+    for (unsigned int i = cut; i < p1.size(); ++i)
+    {
+        buf.erase(std::find(buf.begin(), buf.end(), p1[i]));
+    }
+    for (unsigned int i = 0; i < buf.size(); ++i)
+    {
+        c2.push_back(buf[i]);
+    }
+    for (unsigned int i = cut; i < p1.size(); ++i)
+    {
+        c2.push_back(p1[i]);
+    }
+}
+
+template <typename T>
+void GeneticScheduler<T>::mutation(Gene &m, const Gene &c)
+{
+    Gene                                        buf;
+    std::uniform_int_distribution<unsigned int> dis(0, c.size() - 1);
+    unsigned int                                cut = dis(gen_);
+    Graph<T>                                    g1 = graph_, g2 = graph_;
+    
+    for (unsigned int i = 0; i < cut; ++i)
+    {
+        g1.removeVertex(c[i]);
+    }
+    for (unsigned int i = cut; i < c.size(); ++i)
+    {
+        g2.removeVertex(c[i]);
+    }
+    if (g1.size() > 0)
+    {
+        buf = g1.topoSort(gen_);
+    }
+    if (g2.size() > 0)
+    {
+        m = g2.topoSort(gen_);
+    }
+    for (unsigned int i = cut; i < c.size(); ++i)
+    {
+        m.push_back(buf[i - cut]);
+    }
+}
+
+template <typename T>
+void GeneticScheduler<T>::benchmarkCrossover(const unsigned int nIt)
+{
+    Gene   p1, p2, c1, c2;
+    double neg = 0., eq = 0., pos = 0., total;
+    int    improvement;
+    
+    LOG(Message) << "Benchmarking crossover..." << std::endl;
+    for (unsigned int i = 0; i < nIt; ++i)
+    {
+        p1 = graph_.topoSort(gen_);
+        p2 = graph_.topoSort(gen_);
+        crossover(c1, c2, p1, p2);
+        improvement = (func_(c1) + func_(c2) - func_(p1) - func_(p2))/2;
+        if (improvement < 0) neg++; else if (improvement == 0) eq++; else pos++;
+    }
+    total = neg + eq + pos;
+    LOG(Message) << "  -: " << neg/total << "  =: " << eq/total
+                 << "  +: " << pos/total << std::endl;
+}
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_GeneticScheduler_hpp_

extras/Hadrons/Global.cc

@@ -1,9 +1,11 @@
-/*******************************************************************************
+/*************************************************************************************
+
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: programs/Hadrons/Global.cc
+Source file: extras/Hadrons/Global.cc
 
 Copyright (C) 2015
+Copyright (C) 2016
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 
@@ -21,11 +23,11 @@ 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 
+See the full license in the file "LICENSE" in the top level distribution directory
-directory.
+*************************************************************************************/
-*******************************************************************************/
+/*  END LEGAL */
 
-#include <Hadrons/Global.hpp>
+#include <Grid/Hadrons/Global.hpp>
 
 using namespace Grid;
 using namespace QCD;
@@ -63,3 +65,18 @@ std::string Hadrons::sizeString(long unsigned int bytes)
     
     return std::string(buf);
 }
+
+// type utilities //////////////////////////////////////////////////////////////
+constexpr unsigned int maxNameSize = 1024u;
+
+std::string Hadrons::typeName(const std::type_info *info)
+{
+    char        *buf;
+    std::string name;
+    
+    buf  = abi::__cxa_demangle(info->name(), nullptr, nullptr, nullptr);
+    name = buf;
+    free(buf);
+    
+    return name;
+}

extras/Hadrons/Global.hpp

@@ -1,9 +1,11 @@
-/*******************************************************************************
+/*************************************************************************************
+
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: programs/Hadrons/Global.hpp
+Source file: extras/Hadrons/Global.hpp
 
 Copyright (C) 2015
+Copyright (C) 2016
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 
@@ -21,9 +23,9 @@ 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 
+See the full license in the file "LICENSE" in the top level distribution directory
-directory.
+*************************************************************************************/
-*******************************************************************************/
+/*  END LEGAL */
 
 #ifndef Hadrons_Global_hpp_
 #define Hadrons_Global_hpp_
@@ -31,6 +33,7 @@ directory.
 #include <set>
 #include <stack>
 #include <Grid/Grid.h>
+#include <cxxabi.h>
 
 #define BEGIN_HADRONS_NAMESPACE \
 namespace Grid {\
@@ -39,12 +42,33 @@ namespace Hadrons {\
 using Grid::operator<<;
 #define END_HADRONS_NAMESPACE }}
 
+#define BEGIN_MODULE_NAMESPACE(name)\
+namespace name {\
+using Grid::operator<<;
+#define END_MODULE_NAMESPACE }
+
 /* the 'using Grid::operator<<;' statement prevents a very nasty compilation
- * error with GCC (clang compiles fine without it).
+ * error with GCC 5 (clang & GCC 6 compile fine without it).
  */
 
+// FIXME: find a way to do that in a more general fashion
+#ifndef FIMPL
+#define FIMPL WilsonImplR
+#endif
+
 BEGIN_HADRONS_NAMESPACE
 
+// type aliases
+#define TYPE_ALIASES(FImpl, suffix)\
+typedef FermionOperator<FImpl>                       FMat##suffix;             \
+typedef typename FImpl::FermionField                 FermionField##suffix;     \
+typedef typename FImpl::PropagatorField              PropagatorField##suffix;  \
+typedef typename FImpl::SitePropagator               SitePropagator##suffix;   \
+typedef typename FImpl::DoubledGaugeField            DoubledGaugeField##suffix;\
+typedef std::function<void(FermionField##suffix &,                             \
+                      const FermionField##suffix &)> SolverFn##suffix;
+
+// logger
 class HadronsLogger: public Logger
 {
 public:
@@ -94,34 +118,33 @@ private:\
 // pretty size formating
 std::string sizeString(long unsigned int bytes);
 
-template <typename T>
+// type utilities
-std::string typeName(const T &x)
-{
-    std::string name(typeid(x).name());
-
-    return name;
-}
-
-template <typename T>
-std::string typeName(void)
-{
-    std::string name(typeid(T).name());
-
-    return name;
-}
-
 template <typename T>
 const std::type_info * typeIdPt(const T &x)
 {
     return &typeid(x);
 }
 
+std::string typeName(const std::type_info *info);
+
 template <typename T>
-const std::type_info * typeName(void)
+const std::type_info * typeIdPt(void)
 {
     return &typeid(T);
 }
 
+template <typename T>
+std::string typeName(const T &x)
+{
+    return typeName(typeIdPt(x));
+}
+
+template <typename T>
+std::string typeName(void)
+{
+    return typeName(typeIdPt<T>());
+}
+
 END_HADRONS_NAMESPACE
 
 #endif // Hadrons_Global_hpp_

extras/Hadrons/Graph.hpp

@@ -1,9 +1,11 @@
-/*******************************************************************************
+/*************************************************************************************
+
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: programs/Hadrons/Graph.hpp
+Source file: extras/Hadrons/Graph.hpp
 
 Copyright (C) 2015
+Copyright (C) 2016
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 
@@ -21,14 +23,14 @@ 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 
+See the full license in the file "LICENSE" in the top level distribution directory
-directory.
+*************************************************************************************/
-*******************************************************************************/
+/*  END LEGAL */
 
 #ifndef Hadrons_Graph_hpp_
 #define Hadrons_Graph_hpp_
 
-#include <Hadrons/Global.hpp>
+#include <Grid/Hadrons/Global.hpp>
 
 BEGIN_HADRONS_NAMESPACE
 

extras/Hadrons/HadronsXmlRun.cc

@@ -0,0 +1,80 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/HadronsXmlRun.cc
+
+Copyright (C) 2015
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 <Grid/Hadrons/Application.hpp>
+
+using namespace Grid;
+using namespace QCD;
+using namespace Hadrons;
+
+int main(int argc, char *argv[])
+{
+    // parse command line
+    std::string parameterFileName, scheduleFileName = "";
+    
+    if (argc < 2)
+    {
+        std::cerr << "usage: " << argv[0] << " <parameter file> [<precomputed schedule>] [Grid options]";
+        std::cerr << std::endl;
+        std::exit(EXIT_FAILURE);
+    }
+    parameterFileName = argv[1];
+    if (argc > 2)
+    {
+        if (argv[2][0] != '-')
+        {
+            scheduleFileName = argv[2];
+        }
+    }
+    
+    // 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;
+    
+    // execution
+    Application application(parameterFileName);
+    
+    application.parseParameterFile(parameterFileName);
+    if (!scheduleFileName.empty())
+    {
+        application.loadSchedule(scheduleFileName);
+    }
+    application.run();
+    
+    // epilogue
+    LOG(Message) << "Grid is finalizing now" << std::endl;
+    Grid_finalize();
+    
+    return EXIT_SUCCESS;
+}

extras/Hadrons/HadronsXmlSchedule.cc

@@ -1,9 +1,11 @@
-/*******************************************************************************
+/*************************************************************************************
+
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: programs/Hadrons/Hadrons.cc
+Source file: extras/Hadrons/HadronsXmlSchedule.cc
 
 Copyright (C) 2015
+Copyright (C) 2016
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 
@@ -21,11 +23,11 @@ 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 
+See the full license in the file "LICENSE" in the top level distribution directory
-directory.
+*************************************************************************************/
-*******************************************************************************/
+/*  END LEGAL */
 
-#include <Hadrons/Application.hpp>
+#include <Grid/Hadrons/Application.hpp>
 
 using namespace Grid;
 using namespace QCD;
@@ -34,15 +36,16 @@ using namespace Hadrons;
 int main(int argc, char *argv[])
 {
     // parse command line
-    std::string parameterFileName;
+    std::string parameterFileName, scheduleFileName;
     
-    if (argc < 2)
+    if (argc < 3)
     {
-        std::cerr << "usage: " << argv[0] << " <parameter file> [Grid options]";
+        std::cerr << "usage: " << argv[0] << " <parameter file> <schedule output> [Grid options]";
         std::cerr << std::endl;
         std::exit(EXIT_FAILURE);
     }
     parameterFileName = argv[1];
+    scheduleFileName  = argv[2];
     
     // initialization
     Grid_init(&argc, &argv);
@@ -54,9 +57,12 @@ int main(int argc, char *argv[])
     LOG(Message) << "Grid initialized" << std::endl;
     
     // execution
-    Application application(parameterFileName);
+    Application application;
     
-    application.run();
+    application.parseParameterFile(parameterFileName);
+    application.schedule();
+    application.printSchedule();
+    application.saveSchedule(scheduleFileName);
     
     // epilogue
     LOG(Message) << "Grid is finalizing now" << std::endl;

extras/Hadrons/Makefile.am

@@ -0,0 +1,29 @@
+lib_LIBRARIES = libHadrons.a
+bin_PROGRAMS  = HadronsXmlRun HadronsXmlSchedule
+
+include modules.inc
+
+libHadrons_a_SOURCES = \
+    $(modules_cc)      \
+    Application.cc     \
+    Environment.cc     \
+    Global.cc          \
+    Module.cc
+libHadrons_adir = $(pkgincludedir)/Hadrons
+nobase_libHadrons_a_HEADERS = \
+	$(modules_hpp)            \
+	Application.hpp           \
+	Environment.hpp           \
+	Factory.hpp               \
+	GeneticScheduler.hpp      \
+	Global.hpp                \
+	Graph.hpp                 \
+	Module.hpp                \
+	Modules.hpp               \
+	ModuleFactory.hpp
+
+HadronsXmlRun_SOURCES = HadronsXmlRun.cc
+HadronsXmlRun_LDADD   = libHadrons.a -lGrid
+
+HadronsXmlSchedule_SOURCES = HadronsXmlSchedule.cc
+HadronsXmlSchedule_LDADD   = libHadrons.a -lGrid

extras/Hadrons/Module.cc

@@ -1,9 +1,11 @@
-/*******************************************************************************
+/*************************************************************************************
+
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: programs/Hadrons/Module.cc
+Source file: extras/Hadrons/Module.cc
 
 Copyright (C) 2015
+Copyright (C) 2016
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 
@@ -21,11 +23,11 @@ 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 
+See the full license in the file "LICENSE" in the top level distribution directory
-directory.
+*************************************************************************************/
-*******************************************************************************/
+/*  END LEGAL */
 
-#include <Hadrons/Module.hpp>
+#include <Grid/Hadrons/Module.hpp>
 
 using namespace Grid;
 using namespace QCD;
@@ -51,6 +53,13 @@ Environment & ModuleBase::env(void) const
     return env_;
 }
 
+// get factory registration name if available
+std::string ModuleBase::getRegisteredName(void)
+{
+    HADRON_ERROR("module '" + getName() + "' has a type not registered"
+                 + " in the factory");
+}
+
 // execution ///////////////////////////////////////////////////////////////////
 void ModuleBase::operator()(void)
 {

extras/Hadrons/Module.hpp

@@ -1,9 +1,11 @@
-/*******************************************************************************
+/*************************************************************************************
+
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: programs/Hadrons/Module.hpp
+Source file: extras/Hadrons/Module.hpp
 
 Copyright (C) 2015
+Copyright (C) 2016
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 
@@ -21,20 +23,30 @@ 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 
+See the full license in the file "LICENSE" in the top level distribution directory
-directory.
+*************************************************************************************/
-*******************************************************************************/
+/*  END LEGAL */
 
 #ifndef Hadrons_Module_hpp_
 #define Hadrons_Module_hpp_
 
-#include <Hadrons/Global.hpp>
+#include <Grid/Hadrons/Global.hpp>
-#include <Hadrons/Environment.hpp>
+#include <Grid/Hadrons/Environment.hpp>
 
 BEGIN_HADRONS_NAMESPACE
 
-// module registration macro
+// module registration macros
-#define MODULE_REGISTER(mod)\
+#define MODULE_REGISTER(mod, base)\
+class mod: public base\
+{\
+public:\
+    typedef base Base;\
+    using Base::Base;\
+    virtual std::string getRegisteredName(void)\
+    {\
+        return std::string(#mod);\
+    }\
+};\
 class mod##ModuleRegistrar\
 {\
 public:\
@@ -49,16 +61,39 @@ public:\
 };\
 static mod##ModuleRegistrar mod##ModuleRegistrarInstance;
 
+#define MODULE_REGISTER_NS(mod, base, ns)\
+class mod: public base\
+{\
+public:\
+    typedef base Base;\
+    using Base::Base;\
+    virtual std::string getRegisteredName(void)\
+    {\
+        return std::string(#ns "::" #mod);\
+    }\
+};\
+class ns##mod##ModuleRegistrar\
+{\
+public:\
+    ns##mod##ModuleRegistrar(void)\
+    {\
+        ModuleFactory &modFac = ModuleFactory::getInstance();\
+        modFac.registerBuilder(#ns "::" #mod, [&](const std::string name)\
+                              {\
+                                  return std::unique_ptr<ns::mod>(new ns::mod(name));\
+                              });\
+    }\
+};\
+static ns##mod##ModuleRegistrar ns##mod##ModuleRegistrarInstance;
+
+#define ARG(...) __VA_ARGS__
+
 /******************************************************************************
  *                            Module class                                    *
  ******************************************************************************/
 // base class
 class ModuleBase
 {
-public:
-    // convenient type shortcuts
-    typedef Environment::FMat   FMat;
-    typedef Environment::Solver Solver;
 public:
     // constructor
     ModuleBase(const std::string name);
@@ -67,11 +102,14 @@ public:
     // access
     std::string getName(void) const;
     Environment &env(void) const;
+    // get factory registration name if available
+    virtual std::string getRegisteredName(void);
     // dependencies/products
     virtual std::vector<std::string> getInput(void) = 0;
     virtual std::vector<std::string> getOutput(void) = 0;
     // parse parameters
     virtual void parseParameters(XmlReader &reader, const std::string name) = 0;
+    virtual void saveParameters(XmlWriter &writer, const std::string name) = 0;
     // setup
     virtual void setup(void) {};
     // execution
@@ -95,6 +133,7 @@ public:
     virtual ~Module(void) = default;
     // parse parameters
     virtual void parseParameters(XmlReader &reader, const std::string name);
+    virtual void saveParameters(XmlWriter &writer, const std::string name);
     // parameter access
     const P & par(void) const;
     void      setPar(const P &par);
@@ -115,6 +154,11 @@ public:
     virtual ~Module(void) = default;
     // parse parameters (do nothing)
     virtual void parseParameters(XmlReader &reader, const std::string name) {};
+    virtual void saveParameters(XmlWriter &writer, const std::string name)
+    {
+        push(writer, "options");
+        pop(writer);
+    };
 };
 
 /******************************************************************************
@@ -131,6 +175,12 @@ void Module<P>::parseParameters(XmlReader &reader, const std::string name)
     read(reader, name, par_);
 }
 
+template <typename P>
+void Module<P>::saveParameters(XmlWriter &writer, const std::string name)
+{
+    write(writer, name, par_);
+}
+
 template <typename P>
 const P & Module<P>::par(void) const
 {

extras/Hadrons/ModuleFactory.hpp

@@ -1,9 +1,11 @@
-/*******************************************************************************
+/*************************************************************************************
+
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: programs/Hadrons/ModuleFactory.hpp
+Source file: extras/Hadrons/ModuleFactory.hpp
 
 Copyright (C) 2015
+Copyright (C) 2016
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 
@@ -21,16 +23,16 @@ 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 
+See the full license in the file "LICENSE" in the top level distribution directory
-directory.
+*************************************************************************************/
-*******************************************************************************/
+/*  END LEGAL */
 
 #ifndef Hadrons_ModuleFactory_hpp_
 #define Hadrons_ModuleFactory_hpp_
 
-#include <Hadrons/Global.hpp>
+#include <Grid/Hadrons/Global.hpp>
-#include <Hadrons/Factory.hpp>
+#include <Grid/Hadrons/Factory.hpp>
-#include <Hadrons/Module.hpp>
+#include <Grid/Hadrons/Module.hpp>
 
 BEGIN_HADRONS_NAMESPACE
 

extras/Hadrons/Modules.hpp

@@ -0,0 +1,40 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules.hpp
+
+Copyright (C) 2015
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 <Grid/Hadrons/Modules/MAction/DWF.hpp>
+#include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
+#include <Grid/Hadrons/Modules/MContraction/Baryon.hpp>
+#include <Grid/Hadrons/Modules/MContraction/Meson.hpp>
+#include <Grid/Hadrons/Modules/MGauge/Load.hpp>
+#include <Grid/Hadrons/Modules/MGauge/Random.hpp>
+#include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
+#include <Grid/Hadrons/Modules/MSolver/RBPrecCG.hpp>
+#include <Grid/Hadrons/Modules/MSource/Point.hpp>
+#include <Grid/Hadrons/Modules/MSource/SeqGamma.hpp>
+#include <Grid/Hadrons/Modules/MSource/Z2.hpp>
+#include <Grid/Hadrons/Modules/Quark.hpp>

extras/Hadrons/Modules/MAction/DWF.hpp

@@ -0,0 +1,134 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MAction/DWF.hpp
+
+Copyright (C) 2015
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 */
+
+#ifndef Hadrons_DWF_hpp_
+#define Hadrons_DWF_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                     Domain wall quark action                               *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MAction)
+
+class DWFPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(DWFPar,
+                                    std::string, gauge,
+                                    unsigned int, Ls,
+                                    double      , mass,
+                                    double      , M5);
+};
+
+template <typename FImpl>
+class TDWF: public Module<DWFPar>
+{
+public:
+    TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TDWF(const std::string name);
+    // destructor
+    virtual ~TDWF(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(DWF, TDWF<FIMPL>, MAction);
+
+/******************************************************************************
+ *                        DWF template implementation                         *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TDWF<FImpl>::TDWF(const std::string name)
+: Module<DWFPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TDWF<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().gauge};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TDWF<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TDWF<FImpl>::setup(void)
+{
+    unsigned int size;
+    
+    size = 2*env().template lattice4dSize<typename FImpl::DoubledGaugeField>();
+    env().registerObject(getName(), size, par().Ls);
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TDWF<FImpl>::execute(void)
+{
+    LOG(Message) << "Setting up domain wall fermion matrix with m= "
+                 << par().mass << ", M5= " << par().M5 << " and Ls= "
+                 << par().Ls << " using gauge field '" << par().gauge << "'"
+                 << std::endl;
+    env().createGrid(par().Ls);
+    auto &U      = *env().template getObject<LatticeGaugeField>(par().gauge);
+    auto &g4     = *env().getGrid();
+    auto &grb4   = *env().getRbGrid();
+    auto &g5     = *env().getGrid(par().Ls);
+    auto &grb5   = *env().getRbGrid(par().Ls);
+    FMat *fMatPt = new DomainWallFermion<FImpl>(U, g5, grb5, g4, grb4,
+                                                par().mass, par().M5);
+    env().setObject(getName(), fMatPt);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_DWF_hpp_

extras/Hadrons/Modules/MAction/Wilson.hpp

@@ -0,0 +1,126 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MAction/Wilson.hpp
+
+Copyright (C) 2015
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 */
+
+#ifndef Hadrons_Wilson_hpp_
+#define Hadrons_Wilson_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                            TWilson quark action                            *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MAction)
+
+class WilsonPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonPar,
+                                    std::string, gauge,
+                                    double     , mass);
+};
+
+template <typename FImpl>
+class TWilson: public Module<WilsonPar>
+{
+public:
+    TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TWilson(const std::string name);
+    // destructor
+    virtual ~TWilson(void) = default;
+    // dependencies/products
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(Wilson, TWilson<FIMPL>, MAction);
+
+/******************************************************************************
+ *                     TWilson template implementation                        *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TWilson<FImpl>::TWilson(const std::string name)
+: Module<WilsonPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TWilson<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().gauge};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TWilson<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TWilson<FImpl>::setup(void)
+{
+    unsigned int size;
+    
+    size = 2*env().template lattice4dSize<typename FImpl::DoubledGaugeField>();
+    env().registerObject(getName(), size);
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TWilson<FImpl>::execute()
+{
+    LOG(Message) << "Setting up TWilson fermion matrix with m= " << par().mass
+                 << " using gauge field '" << par().gauge << "'" << std::endl;
+    auto &U      = *env().template getObject<LatticeGaugeField>(par().gauge);
+    auto &grid   = *env().getGrid();
+    auto &gridRb = *env().getRbGrid();
+    FMat *fMatPt = new WilsonFermion<FImpl>(U, grid, gridRb, par().mass);
+    env().setObject(getName(), fMatPt);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Wilson_hpp_

extras/Hadrons/Modules/MContraction/Baryon.hpp

@@ -0,0 +1,131 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MContraction/Baryon.hpp
+
+Copyright (C) 2015
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 */
+
+#ifndef Hadrons_Baryon_hpp_
+#define Hadrons_Baryon_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                               Baryon                                       *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MContraction)
+
+class BaryonPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(BaryonPar,
+                                    std::string, q1,
+                                    std::string, q2,
+                                    std::string, q3,
+                                    std::string, output);
+};
+
+template <typename FImpl1, typename FImpl2, typename FImpl3>
+class TBaryon: public Module<BaryonPar>
+{
+public:
+    TYPE_ALIASES(FImpl1, 1);
+    TYPE_ALIASES(FImpl2, 2);
+    TYPE_ALIASES(FImpl3, 3);
+    class Result: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
+                                        std::vector<std::vector<std::vector<Complex>>>, corr);
+    };
+public:
+    // constructor
+    TBaryon(const std::string name);
+    // destructor
+    virtual ~TBaryon(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(Baryon, ARG(TBaryon<FIMPL, FIMPL, FIMPL>), MContraction);
+
+/******************************************************************************
+ *                         TBaryon implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2, typename FImpl3>
+TBaryon<FImpl1, FImpl2, FImpl3>::TBaryon(const std::string name)
+: Module<BaryonPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2, typename FImpl3>
+std::vector<std::string> TBaryon<FImpl1, FImpl2, FImpl3>::getInput(void)
+{
+    std::vector<std::string> input = {par().q1, par().q2, par().q3};
+    
+    return input;
+}
+
+template <typename FImpl1, typename FImpl2, typename FImpl3>
+std::vector<std::string> TBaryon<FImpl1, FImpl2, FImpl3>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2, typename FImpl3>
+void TBaryon<FImpl1, FImpl2, FImpl3>::execute(void)
+{
+    LOG(Message) << "Computing baryon contractions '" << getName() << "' using"
+                 << " quarks '" << par().q1 << "', '" << par().q2 << "', and '"
+                 << par().q3 << "'" << std::endl;
+    
+    XmlWriter             writer(par().output);
+    PropagatorField1      &q1 = *env().template getObject<PropagatorField1>(par().q1);
+    PropagatorField2      &q2 = *env().template getObject<PropagatorField2>(par().q2);
+    PropagatorField3      &q3 = *env().template getObject<PropagatorField3>(par().q2);
+    LatticeComplex        c(env().getGrid());
+    Result                result;
+    
+    // FIXME: do contractions
+    
+    write(writer, "meson", result);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Baryon_hpp_

extras/Hadrons/Modules/MContraction/Meson.hpp

@@ -1,9 +1,11 @@
-/*******************************************************************************
+/*************************************************************************************
+
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: programs/Hadrons/CMeson.cc
+Source file: extras/Hadrons/Modules/MContraction/Meson.hpp
 
 Copyright (C) 2015
+Copyright (C) 2016
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 
@@ -21,33 +23,79 @@ 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 
+See the full license in the file "LICENSE" in the top level distribution directory
-directory.
+*************************************************************************************/
-*******************************************************************************/
+/*  END LEGAL */
 
-#include <Hadrons/Modules/CMeson.hpp>
+#ifndef Hadrons_Meson_hpp_
+#define Hadrons_Meson_hpp_
 
-using namespace Grid;
+#include <Grid/Hadrons/Global.hpp>
-using namespace QCD;
+#include <Grid/Hadrons/Module.hpp>
-using namespace Hadrons;
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
 
 /******************************************************************************
- *                          CMeson implementation                             *
+ *                                TMeson                                       *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MContraction)
+
+class MesonPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(MesonPar,
+                                    std::string, q1,
+                                    std::string, q2,
+                                    std::string, output);
+};
+
+template <typename FImpl1, typename FImpl2>
+class TMeson: public Module<MesonPar>
+{
+public:
+    TYPE_ALIASES(FImpl1, 1);
+    TYPE_ALIASES(FImpl2, 2);
+    class Result: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
+                                        std::vector<std::vector<std::vector<Complex>>>, corr);
+    };
+public:
+    // constructor
+    TMeson(const std::string name);
+    // destructor
+    virtual ~TMeson(void) = default;
+    // dependencies/products
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(Meson, ARG(TMeson<FIMPL, FIMPL>), MContraction);
+
+/******************************************************************************
+ *                           TMeson implementation                            *
  ******************************************************************************/
 // constructor /////////////////////////////////////////////////////////////////
-CMeson::CMeson(const std::string name)
+template <typename FImpl1, typename FImpl2>
-: Module<CMesonPar>(name)
+TMeson<FImpl1, FImpl2>::TMeson(const std::string name)
+: Module<MesonPar>(name)
 {}
 
 // dependencies/products ///////////////////////////////////////////////////////
-std::vector<std::string> CMeson::getInput(void)
+template <typename FImpl1, typename FImpl2>
+std::vector<std::string> TMeson<FImpl1, FImpl2>::getInput(void)
 {
     std::vector<std::string> input = {par().q1, par().q2};
     
     return input;
 }
 
-std::vector<std::string> CMeson::getOutput(void)
+template <typename FImpl1, typename FImpl2>
+std::vector<std::string> TMeson<FImpl1, FImpl2>::getOutput(void)
 {
     std::vector<std::string> output = {getName()};
     
@@ -55,20 +103,20 @@ std::vector<std::string> CMeson::getOutput(void)
 }
 
 // execution ///////////////////////////////////////////////////////////////////
-void CMeson::execute(void)
+template <typename FImpl1, typename FImpl2>
+void TMeson<FImpl1, FImpl2>::execute(void)
 {
-    LOG(Message) << "Computing meson contraction '" << getName() << "' using"
+    LOG(Message) << "Computing meson contractions '" << getName() << "' using"
                  << " quarks '" << par().q1 << "' and '" << par().q2 << "'"
                  << std::endl;
     
     XmlWriter             writer(par().output);
-    LatticePropagator     &q1 = *env().getObject<LatticePropagator>(par().q1);
+    PropagatorField1      &q1 = *env().template getObject<PropagatorField1>(par().q1);
-    LatticePropagator     &q2 = *env().getObject<LatticePropagator>(par().q2);
+    PropagatorField2      &q2 = *env().template getObject<PropagatorField2>(par().q2);
     LatticeComplex        c(env().getGrid());
     SpinMatrix            g[Ns*Ns], g5;
     std::vector<TComplex> buf;
     Result                result;
-    unsigned int          nt = env().getGrid()->GlobalDimensions()[Tp];
     
     g5 = makeGammaProd(Ns*Ns - 1);
     result.corr.resize(Ns*Ns);
@@ -92,3 +140,9 @@ void CMeson::execute(void)
     }
     write(writer, "meson", result);
 }
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Meson_hpp_

extras/Hadrons/Modules/MGauge/Load.cc

@@ -1,8 +1,10 @@
-/*******************************************************************************
+/*************************************************************************************
+
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: programs/Hadrons/GLoad.cc
+Source file: extras/Hadrons/Modules/MGauge/Load.cc
 
+Copyright (C) 2015
 Copyright (C) 2016
 
 Author: Antonin Portelli <antonin.portelli@me.com>
@@ -21,32 +23,33 @@ 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 
+See the full license in the file "LICENSE" in the top level distribution directory
-directory.
+*************************************************************************************/
-*******************************************************************************/
+/*  END LEGAL */
 
-#include <Hadrons/Modules/GLoad.hpp>
+#include <Grid/Hadrons/Modules/MGauge/Load.hpp>
 
 using namespace Grid;
 using namespace Hadrons;
+using namespace MGauge;
 
 /******************************************************************************
-*                          GLoad implementation                               *
+*                           TLoad implementation                               *
 ******************************************************************************/
 // constructor /////////////////////////////////////////////////////////////////
-GLoad::GLoad(const std::string name)
+TLoad::TLoad(const std::string name)
-: Module<GLoadPar>(name)
+: Module<LoadPar>(name)
 {}
 
 // dependencies/products ///////////////////////////////////////////////////////
-std::vector<std::string> GLoad::getInput(void)
+std::vector<std::string> TLoad::getInput(void)
 {
     std::vector<std::string> in;
     
     return in;
 }
 
-std::vector<std::string> GLoad::getOutput(void)
+std::vector<std::string> TLoad::getOutput(void)
 {
     std::vector<std::string> out = {getName()};
     
@@ -54,13 +57,13 @@ std::vector<std::string> GLoad::getOutput(void)
 }
 
 // setup ///////////////////////////////////////////////////////////////////////
-void GLoad::setup(void)
+void TLoad::setup(void)
 {
     env().registerLattice<LatticeGaugeField>(getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
-void GLoad::execute(void)
+void TLoad::execute(void)
 {
     NerscField  header;
     std::string fileName = par().file + "."

extras/Hadrons/Modules/MGauge/Load.hpp

@@ -1,8 +1,10 @@
-/*******************************************************************************
+/*************************************************************************************
+
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: programs/Hadrons/GLoad.hpp
+Source file: extras/Hadrons/Modules/MGauge/Load.hpp
 
+Copyright (C) 2015
 Copyright (C) 2016
 
 Author: Antonin Portelli <antonin.portelli@me.com>
@@ -21,36 +23,38 @@ 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 
+See the full license in the file "LICENSE" in the top level distribution directory
-directory.
+*************************************************************************************/
-*******************************************************************************/
+/*  END LEGAL */
 
-#ifndef Hadrons_GLoad_hpp_
+#ifndef Hadrons_Load_hpp_
-#define Hadrons_GLoad_hpp_
+#define Hadrons_Load_hpp_
 
-#include <Hadrons/Global.hpp>
+#include <Grid/Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
+#include <Grid/Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
 
 BEGIN_HADRONS_NAMESPACE
 
 /******************************************************************************
  *                         Load a NERSC configuration                         *
  ******************************************************************************/
-class GLoadPar: Serializable
+BEGIN_MODULE_NAMESPACE(MGauge)
+
+class LoadPar: Serializable
 {
 public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(GLoadPar,
+    GRID_SERIALIZABLE_CLASS_MEMBERS(LoadPar,
                                     std::string, file);
 };
 
-class GLoad: public Module<GLoadPar>
+class TLoad: public Module<LoadPar>
 {
 public:
     // constructor
-    GLoad(const std::string name);
+    TLoad(const std::string name);
     // destructor
-    virtual ~GLoad(void) = default;
+    virtual ~TLoad(void) = default;
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -60,8 +64,10 @@ public:
     virtual void execute(void);
 };
 
-MODULE_REGISTER(GLoad);
+MODULE_REGISTER_NS(Load, TLoad, MGauge);
+
+END_MODULE_NAMESPACE
 
 END_HADRONS_NAMESPACE
 
-#endif // Hadrons_GLoad_hpp_
+#endif // Hadrons_Load_hpp_

extras/Hadrons/Modules/MGauge/Random.cc

@@ -1,8 +1,10 @@
-/*******************************************************************************
+/*************************************************************************************
+
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: programs/Hadrons/GRandom.cc
+Source file: extras/Hadrons/Modules/MGauge/Random.cc
 
+Copyright (C) 2015
 Copyright (C) 2016
 
 Author: Antonin Portelli <antonin.portelli@me.com>
@@ -21,30 +23,31 @@ 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 
+See the full license in the file "LICENSE" in the top level distribution directory
-directory.
+*************************************************************************************/
-*******************************************************************************/
+/*  END LEGAL */
 
-#include <Hadrons/Modules/GRandom.hpp>
+#include <Grid/Hadrons/Modules/MGauge/Random.hpp>
 
 using namespace Grid;
 using namespace Hadrons;
+using namespace MGauge;
 
 /******************************************************************************
-*                          GRandom implementation                             *
+*                           TRandom implementation                            *
 ******************************************************************************/
 // constructor /////////////////////////////////////////////////////////////////
-GRandom::GRandom(const std::string name)
+TRandom::TRandom(const std::string name)
 : Module<NoPar>(name)
 {}
 
 // dependencies/products ///////////////////////////////////////////////////////
-std::vector<std::string> GRandom::getInput(void)
+std::vector<std::string> TRandom::getInput(void)
 {
     return std::vector<std::string>();
 }
 
-std::vector<std::string> GRandom::getOutput(void)
+std::vector<std::string> TRandom::getOutput(void)
 {
     std::vector<std::string> out = {getName()};
     
@@ -52,13 +55,13 @@ std::vector<std::string> GRandom::getOutput(void)
 }
 
 // setup ///////////////////////////////////////////////////////////////////////
-void GRandom::setup(void)
+void TRandom::setup(void)
 {
     env().registerLattice<LatticeGaugeField>(getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
-void GRandom::execute(void)
+void TRandom::execute(void)
 {
     LOG(Message) << "Generating random gauge configuration" << std::endl;
     LatticeGaugeField &U = *env().createLattice<LatticeGaugeField>(getName());

extras/Hadrons/Modules/MGauge/Random.hpp

@@ -1,8 +1,10 @@
-/*******************************************************************************
+/*************************************************************************************
+
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: programs/Hadrons/GRandom.hpp
+Source file: extras/Hadrons/Modules/MGauge/Random.hpp
 
+Copyright (C) 2015
 Copyright (C) 2016
 
 Author: Antonin Portelli <antonin.portelli@me.com>
@@ -21,29 +23,31 @@ 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 
+See the full license in the file "LICENSE" in the top level distribution directory
-directory.
+*************************************************************************************/
-*******************************************************************************/
+/*  END LEGAL */
 
-#ifndef Hadrons_GRandom_hpp_
+#ifndef Hadrons_Random_hpp_
-#define Hadrons_GRandom_hpp_
+#define Hadrons_Random_hpp_
 
-#include <Hadrons/Global.hpp>
+#include <Grid/Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
+#include <Grid/Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
 
 BEGIN_HADRONS_NAMESPACE
 
 /******************************************************************************
  *                             Random gauge                                   *
  ******************************************************************************/
-class GRandom: public Module<NoPar>
+BEGIN_MODULE_NAMESPACE(MGauge)
+
+class TRandom: public Module<NoPar>
 {
 public:
     // constructor
-    GRandom(const std::string name);
+    TRandom(const std::string name);
     // destructor
-    virtual ~GRandom(void) = default;
+    virtual ~TRandom(void) = default;
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -53,8 +57,10 @@ public:
     virtual void execute(void);
 };
 
-MODULE_REGISTER(GRandom);
+MODULE_REGISTER_NS(Random, TRandom, MGauge);
+
+END_MODULE_NAMESPACE
 
 END_HADRONS_NAMESPACE
 
-#endif // Hadrons_GRandom_hpp_
+#endif // Hadrons_Random_hpp_

extras/Hadrons/Modules/MGauge/Unit.cc

@@ -1,8 +1,10 @@
-/*******************************************************************************
+/*************************************************************************************
+
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: programs/Hadrons/GUnit.cc
+Source file: extras/Hadrons/Modules/MGauge/Unit.cc
 
+Copyright (C) 2015
 Copyright (C) 2016
 
 Author: Antonin Portelli <antonin.portelli@me.com>
@@ -21,30 +23,31 @@ 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 
+See the full license in the file "LICENSE" in the top level distribution directory
-directory.
+*************************************************************************************/
-*******************************************************************************/
+/*  END LEGAL */
 
-#include <Hadrons/Modules/GUnit.hpp>
+#include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
 
 using namespace Grid;
 using namespace Hadrons;
+using namespace MGauge;
 
 /******************************************************************************
-*                            GUnit implementation                             *
+*                            TUnit implementation                             *
 ******************************************************************************/
 // constructor /////////////////////////////////////////////////////////////////
-GUnit::GUnit(const std::string name)
+TUnit::TUnit(const std::string name)
 : Module<NoPar>(name)
 {}
 
 // dependencies/products ///////////////////////////////////////////////////////
-std::vector<std::string> GUnit::getInput(void)
+std::vector<std::string> TUnit::getInput(void)
 {
     return std::vector<std::string>();
 }
 
-std::vector<std::string> GUnit::getOutput(void)
+std::vector<std::string> TUnit::getOutput(void)
 {
     std::vector<std::string> out = {getName()};
     
@@ -52,13 +55,13 @@ std::vector<std::string> GUnit::getOutput(void)
 }
 
 // setup ///////////////////////////////////////////////////////////////////////
-void GUnit::setup(void)
+void TUnit::setup(void)
 {
     env().registerLattice<LatticeGaugeField>(getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
-void GUnit::execute(void)
+void TUnit::execute(void)
 {
     LOG(Message) << "Creating unit gauge configuration" << std::endl;
     LatticeGaugeField &U = *env().createLattice<LatticeGaugeField>(getName());

extras/Hadrons/Modules/MGauge/Unit.hpp

@@ -1,8 +1,10 @@
-/*******************************************************************************
+/*************************************************************************************
+
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: programs/Hadrons/GUnit.hpp
+Source file: extras/Hadrons/Modules/MGauge/Unit.hpp
 
+Copyright (C) 2015
 Copyright (C) 2016
 
 Author: Antonin Portelli <antonin.portelli@me.com>
@@ -21,29 +23,31 @@ 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 
+See the full license in the file "LICENSE" in the top level distribution directory
-directory.
+*************************************************************************************/
-*******************************************************************************/
+/*  END LEGAL */
 
-#ifndef Hadrons_GUnit_hpp_
+#ifndef Hadrons_Unit_hpp_
-#define Hadrons_GUnit_hpp_
+#define Hadrons_Unit_hpp_
 
-#include <Hadrons/Global.hpp>
+#include <Grid/Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
+#include <Grid/Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
 
 BEGIN_HADRONS_NAMESPACE
 
 /******************************************************************************
  *                              Unit gauge                                    *
  ******************************************************************************/
-class GUnit: public Module<NoPar>
+BEGIN_MODULE_NAMESPACE(MGauge)
+
+class TUnit: public Module<NoPar>
 {
 public:
     // constructor
-    GUnit(const std::string name);
+    TUnit(const std::string name);
     // destructor
-    virtual ~GUnit(void) = default;
+    virtual ~TUnit(void) = default;
     // dependencies/products
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -53,8 +57,10 @@ public:
     virtual void execute(void);
 };
 
-MODULE_REGISTER(GUnit);
+MODULE_REGISTER_NS(Unit, TUnit, MGauge);
+
+END_MODULE_NAMESPACE
 
 END_HADRONS_NAMESPACE
 
-#endif // Hadrons_GUnit_hpp_
+#endif // Hadrons_Unit_hpp_

extras/Hadrons/Modules/MSolver/RBPrecCG.hpp

@@ -0,0 +1,132 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MSolver/RBPrecCG.hpp
+
+Copyright (C) 2015
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 */
+
+#ifndef Hadrons_RBPrecCG_hpp_
+#define Hadrons_RBPrecCG_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                     Schur red-black preconditioned CG                      *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MSolver)
+
+class RBPrecCGPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(RBPrecCGPar,
+                                    std::string, action,
+                                    double     , residual);
+};
+
+template <typename FImpl>
+class TRBPrecCG: public Module<RBPrecCGPar>
+{
+public:
+    TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TRBPrecCG(const std::string name);
+    // destructor
+    virtual ~TRBPrecCG(void) = default;
+    // dependencies/products
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(RBPrecCG, TRBPrecCG<FIMPL>, MSolver);
+
+/******************************************************************************
+ *                      TRBPrecCG template implementation                     *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TRBPrecCG<FImpl>::TRBPrecCG(const std::string name)
+: Module(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TRBPrecCG<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().action};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TRBPrecCG<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TRBPrecCG<FImpl>::setup(void)
+{
+    auto Ls = env().getObjectLs(par().action);
+    
+    env().registerObject(getName(), 0, Ls);
+    env().addOwnership(getName(), par().action);
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TRBPrecCG<FImpl>::execute(void)
+{
+    auto &mat   = *(env().template getObject<FMat>(par().action));
+    auto solver = [&mat, this](FermionField &sol, const FermionField &source)
+    {
+        ConjugateGradient<FermionField>           cg(par().residual, 10000);
+        SchurRedBlackDiagMooeeSolve<FermionField> schurSolver(cg);
+        
+        schurSolver(mat, source, sol);
+    };
+    
+    LOG(Message) << "setting up Schur red-black preconditioned CG for"
+                 << " action '" << par().action << "' with residual "
+                 << par().residual << std::endl;
+    env().setObject(getName(), new SolverFn(solver));
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_RBPrecCG_hpp_

extras/Hadrons/Modules/MSource/Point.hpp

@@ -0,0 +1,135 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MSource/Point.hpp
+
+Copyright (C) 2015
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 */
+
+#ifndef Hadrons_Point_hpp_
+#define Hadrons_Point_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/*
+ 
+ Point source
+ ------------
+ * src_x = delta_x,position
+ 
+ * options:
+ - position: space-separated integer sequence (e.g. "0 1 1 0")
+ 
+ */
+
+/******************************************************************************
+ *                                  TPoint                                     *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MSource)
+
+class PointPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(PointPar,
+                                    std::string, position);
+};
+
+template <typename FImpl>
+class TPoint: public Module<PointPar>
+{
+public:
+    TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TPoint(const std::string name);
+    // destructor
+    virtual ~TPoint(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(Point, TPoint<FIMPL>, MSource);
+
+/******************************************************************************
+ *                       TPoint template implementation                       *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TPoint<FImpl>::TPoint(const std::string name)
+: Module<PointPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TPoint<FImpl>::getInput(void)
+{
+    std::vector<std::string> in;
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TPoint<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TPoint<FImpl>::setup(void)
+{
+    env().template registerLattice<PropagatorField>(getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TPoint<FImpl>::execute(void)
+{
+    std::vector<int> position = strToVec<int>(par().position);
+    typename SitePropagator::scalar_object id;
+    
+    LOG(Message) << "Creating point source at position [" << par().position
+                 << "]" << std::endl;
+    PropagatorField &src = *env().template createLattice<PropagatorField>(getName());
+    id  = 1.;
+    src = zero;
+    pokeSite(id, src, position);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Point_hpp_

extras/Hadrons/Modules/MSource/SeqGamma.hpp

@@ -0,0 +1,164 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MSource/SeqGamma.hpp
+
+Copyright (C) 2015
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 */
+
+#ifndef Hadrons_SeqGamma_hpp_
+#define Hadrons_SeqGamma_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/*
+ 
+ Sequential source
+ -----------------------------
+ * src_x = q_x * theta(x_3 - tA) * theta(tB - x_3) * gamma * exp(i x.mom)
+ 
+ * options:
+ - q: input propagator (string)
+ - tA: begin timeslice (integer)
+ - tB: end timesilce (integer)
+ - gamma: gamma product to insert (integer)
+ - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.")
+ 
+ */
+
+/******************************************************************************
+ *                         SeqGamma                                 *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MSource)
+
+class SeqGammaPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(SeqGammaPar,
+                                    std::string,  q,
+                                    unsigned int, tA,
+                                    unsigned int, tB,
+                                    unsigned int, gamma,
+                                    std::string,  mom);
+};
+
+template <typename FImpl>
+class TSeqGamma: public Module<SeqGammaPar>
+{
+public:
+    TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TSeqGamma(const std::string name);
+    // destructor
+    virtual ~TSeqGamma(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(SeqGamma, TSeqGamma<FIMPL>, MSource);
+
+/******************************************************************************
+ *                         TSeqGamma implementation                           *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TSeqGamma<FImpl>::TSeqGamma(const std::string name)
+: Module<SeqGammaPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TSeqGamma<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().q};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TSeqGamma<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TSeqGamma<FImpl>::setup(void)
+{
+    env().template registerLattice<PropagatorField>(getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TSeqGamma<FImpl>::execute(void)
+{
+    if (par().tA == par().tB)
+    {
+        LOG(Message) << "Generating gamma_" << par().gamma
+                     << " sequential source at t= " << par().tA << std::endl;
+    }
+    else
+    {
+        LOG(Message) << "Generating gamma_" << par().gamma
+                     << " sequential source for "
+                     << par().tA << " <= t <= " << par().tB << std::endl;
+    }
+    PropagatorField &src = *env().template createLattice<PropagatorField>(getName());
+    PropagatorField &q   = *env().template getObject<PropagatorField>(par().q);
+    Lattice<iScalar<vInteger>> t(env().getGrid());
+    LatticeComplex             ph(env().getGrid()), coor(env().getGrid());
+    SpinMatrix                 g;
+    std::vector<Real>          p;
+    Complex                    i(0.0,1.0);
+    
+    g  = makeGammaProd(par().gamma);
+    p  = strToVec<Real>(par().mom);
+    ph = zero;
+    for(unsigned int mu = 0; mu < Nd; mu++)
+    {
+        LatticeCoordinate(coor, mu);
+        ph = ph + p[mu]*coor;
+    }
+    ph = exp(i*ph);
+    LatticeCoordinate(t, Tp);
+    src = where((t >= par().tA) and (t <= par().tB), g*ph*q, 0.*q);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_SeqGamma_hpp_

extras/Hadrons/Modules/MSource/Z2.hpp

@@ -0,0 +1,151 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MSource/Z2.hpp
+
+Copyright (C) 2015
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 */
+
+#ifndef Hadrons_Z2_hpp_
+#define Hadrons_Z2_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/*
+ 
+ Z_2 stochastic source
+ -----------------------------
+ * src_x = eta_x * theta(x_3 - tA) * theta(tB - x_3)
+ 
+ the eta_x are independent uniform random numbers in {+/- 1 +/- i}
+ 
+ * options:
+ - tA: begin timeslice (integer)
+ - tB: end timesilce (integer)
+ 
+ */
+ 
+/******************************************************************************
+ *                          Z2 stochastic source                              *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MSource)
+
+class Z2Par: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(Z2Par,
+                                    unsigned int, tA,
+                                    unsigned int, tB);
+};
+
+template <typename FImpl>
+class TZ2: public Module<Z2Par>
+{
+public:
+    TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TZ2(const std::string name);
+    // destructor
+    virtual ~TZ2(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(Z2, TZ2<FIMPL>, MSource);
+
+/******************************************************************************
+ *                       TZ2 template implementation                          *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TZ2<FImpl>::TZ2(const std::string name)
+: Module<Z2Par>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TZ2<FImpl>::getInput(void)
+{
+    std::vector<std::string> in;
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TZ2<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TZ2<FImpl>::setup(void)
+{
+    env().template registerLattice<PropagatorField>(getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TZ2<FImpl>::execute(void)
+{
+    Lattice<iScalar<vInteger>> t(env().getGrid());
+    LatticeComplex             eta(env().getGrid());
+    Complex                    shift(1., 1.);
+    
+    if (par().tA == par().tB)
+    {
+        LOG(Message) << "Generating Z_2 wall source at t= " << par().tA
+        << std::endl;
+    }
+    else
+    {
+        LOG(Message) << "Generating Z_2 band for " << par().tA << " <= t <= "
+        << par().tB << std::endl;
+    }
+    PropagatorField &src = *env().template createLattice<PropagatorField>(getName());
+    LatticeCoordinate(t, Tp);
+    bernoulli(*env().get4dRng(), eta);
+    eta = (2.*eta - shift)*(1./::sqrt(2.));
+    eta = where((t >= par().tA) and (t <= par().tB), eta, 0.*eta);
+    src = 1.;
+    src = src*eta;
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Z2_hpp_

extras/Hadrons/Modules/Quark.hpp

@@ -1,9 +1,11 @@
-/*******************************************************************************
+/*************************************************************************************
+
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: programs/Hadrons/MQuark.cc
+Source file: extras/Hadrons/Modules/Quark.hpp
 
 Copyright (C) 2015
+Copyright (C) 2016
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 
@@ -21,33 +23,74 @@ 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 
+See the full license in the file "LICENSE" in the top level distribution directory
-directory.
+*************************************************************************************/
-*******************************************************************************/
+/*  END LEGAL */
 
-#include <Hadrons/Modules/MQuark.hpp>
+#ifndef Hadrons_Quark_hpp_
+#define Hadrons_Quark_hpp_
 
-using namespace Grid;
+#include <Grid/Hadrons/Global.hpp>
-using namespace QCD;
+#include <Grid/Hadrons/Module.hpp>
-using namespace Hadrons;
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
 
 /******************************************************************************
- *                          MQuark implementation                             *
+ *                               TQuark                                       *
+ ******************************************************************************/
+class QuarkPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(QuarkPar,
+                                    std::string, source,
+                                    std::string, solver);
+};
+
+template <typename FImpl>
+class TQuark: public Module<QuarkPar>
+{
+public:
+    TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TQuark(const std::string name);
+    // destructor
+    virtual ~TQuark(void) = default;
+    // dependencies/products
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+private:
+    unsigned int Ls_;
+    SolverFn     *solver_{nullptr};
+};
+
+MODULE_REGISTER(Quark, TQuark<FIMPL>);
+
+/******************************************************************************
+ *                          TQuark implementation                             *
  ******************************************************************************/
 // constructor /////////////////////////////////////////////////////////////////
-MQuark::MQuark(const std::string name)
+template <typename FImpl>
+TQuark<FImpl>::TQuark(const std::string name)
 : Module(name)
 {}
 
 // dependencies/products ///////////////////////////////////////////////////////
-std::vector<std::string> MQuark::getInput(void)
+template <typename FImpl>
+std::vector<std::string> TQuark<FImpl>::getInput(void)
 {
     std::vector<std::string> in = {par().source, par().solver};
     
     return in;
 }
 
-std::vector<std::string> MQuark::getOutput(void)
+template <typename FImpl>
+std::vector<std::string> TQuark<FImpl>::getOutput(void)
 {
     std::vector<std::string> out = {getName(), getName() + "_5d"};
     
@@ -55,31 +98,33 @@ std::vector<std::string> MQuark::getOutput(void)
 }
 
 // setup ///////////////////////////////////////////////////////////////////////
-void MQuark::setup(void)
+template <typename FImpl>
+void TQuark<FImpl>::setup(void)
 {
     Ls_ = env().getObjectLs(par().solver);
-    env().registerLattice<LatticePropagator>(getName());
+    env().template registerLattice<PropagatorField>(getName());
     if (Ls_ > 1)
     {
-        env().registerLattice<LatticePropagator>(getName() + "_5d", Ls_);
+        env().template registerLattice<PropagatorField>(getName() + "_5d", Ls_);
     }
 }
 
 // execution ///////////////////////////////////////////////////////////////////
-void MQuark::execute(void)
+template <typename FImpl>
+void TQuark<FImpl>::execute(void)
 {
-    LatticeFermion source(env().getGrid(Ls_)), sol(env().getGrid(Ls_)),
-                   tmp(env().getGrid());
-    std::string    propName = (Ls_ == 1) ? getName() : (getName() + "_5d");
-    
     LOG(Message) << "Computing quark propagator '" << getName() << "'"
                  << std::endl;
-    LatticePropagator   &prop    = *env().createLattice<LatticePropagator>(propName);
+    
-    LatticePropagator   &fullSrc = *env().getObject<LatticePropagator>(par().source);
+    FermionField    source(env().getGrid(Ls_)), sol(env().getGrid(Ls_)),
-    Environment::Solver &solver  = *env().getObject<Environment::Solver>(par().solver);
+                    tmp(env().getGrid());
+    std::string     propName = (Ls_ == 1) ? getName() : (getName() + "_5d");
+    PropagatorField &prop    = *env().template createLattice<PropagatorField>(propName);
+    PropagatorField &fullSrc = *env().template getObject<PropagatorField>(par().source);
+    SolverFn        &solver  = *env().template getObject<SolverFn>(par().solver);
     if (Ls_ > 1)
     {
-        env().createLattice<LatticePropagator>(getName());
+        env().template createLattice<PropagatorField>(getName());
     }
     
     LOG(Message) << "Inverting using solver '" << par().solver
@@ -88,7 +133,7 @@ void MQuark::execute(void)
     for (unsigned int c = 0; c < Nc; ++c)
     {
         LOG(Message) << "Inversion for spin= " << s << ", color= " << c
-                     << std::endl;
+        << std::endl;
         // source conversion for 4D sources
         if (!env().isObject5d(par().source))
         {
@@ -124,7 +169,8 @@ void MQuark::execute(void)
         // create 4D propagators from 5D one if necessary
         if (Ls_ > 1)
         {
-            LatticePropagator &p4d = *env().getObject<LatticePropagator>(getName());
+            PropagatorField &p4d =
+                *env().template getObject<PropagatorField>(getName());
             
             axpby_ssp_pminus(sol, 0., sol, 1., sol, 0, 0);
             axpby_ssp_pplus(sol, 0., sol, 1., sol, 0, Ls_-1);
@@ -133,3 +179,7 @@ void MQuark::execute(void)
         }
     }
 }
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Quark_hpp_

extras/Hadrons/Modules/templates/Module.cc.template

@@ -1,25 +1,25 @@
-#include <Hadrons/___FILEBASENAME___.hpp>
+#include <Grid/Hadrons/Modules/___FILEBASENAME___.hpp>
 
 using namespace Grid;
 using namespace Hadrons;
 
 /******************************************************************************
-*                  ___FILEBASENAME___ implementation                             *
+*                  T___FILEBASENAME___ implementation                             *
 ******************************************************************************/
 // constructor /////////////////////////////////////////////////////////////////
-___FILEBASENAME___::___FILEBASENAME___(const std::string name)
+T___FILEBASENAME___::T___FILEBASENAME___(const std::string name)
 : Module<___FILEBASENAME___Par>(name)
 {}
 
 // dependencies/products ///////////////////////////////////////////////////////
-std::vector<std::string> ___FILEBASENAME___::getInput(void)
+std::vector<std::string> T___FILEBASENAME___::getInput(void)
 {
     std::vector<std::string> in;
     
     return in;
 }
 
-std::vector<std::string> ___FILEBASENAME___::getOutput(void)
+std::vector<std::string> T___FILEBASENAME___::getOutput(void)
 {
     std::vector<std::string> out = {getName()};
     
@@ -27,13 +27,13 @@ std::vector<std::string> ___FILEBASENAME___::getOutput(void)
 }
 
 // setup ///////////////////////////////////////////////////////////////////////
-void ___FILEBASENAME___::setup(void)
+void T___FILEBASENAME___::setup(void)
 {
 
 }
 
 // execution ///////////////////////////////////////////////////////////////////
-void ___FILEBASENAME___::execute(void)
+void T___FILEBASENAME___::execute(void)
 {
 
 }

extras/Hadrons/Modules/templates/Module.hpp.template

@@ -1,9 +1,9 @@
 #ifndef Hadrons____FILEBASENAME____hpp_
 #define Hadrons____FILEBASENAME____hpp_
 
-#include <Hadrons/Global.hpp>
+#include <Grid/Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
+#include <Grid/Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
 
 BEGIN_HADRONS_NAMESPACE
 
@@ -17,13 +17,13 @@ public:
                                     unsigned int, i);
 };
 
-class ___FILEBASENAME___: public Module<___FILEBASENAME___Par>
+class T___FILEBASENAME___: public Module<___FILEBASENAME___Par>
 {
 public:
     // constructor
-    ___FILEBASENAME___(const std::string name);
+    T___FILEBASENAME___(const std::string name);
     // destructor
-    virtual ~___FILEBASENAME___(void) = default;
+    virtual ~T___FILEBASENAME___(void) = default;
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -33,7 +33,7 @@ public:
     virtual void execute(void);
 };
 
-MODULE_REGISTER(___FILEBASENAME___);
+MODULE_REGISTER(___FILEBASENAME___, T___FILEBASENAME___);
 
 END_HADRONS_NAMESPACE
 

extras/Hadrons/Modules/templates/Module_in_NS.cc.template

@@ -0,0 +1,40 @@
+#include <Grid/Hadrons/Modules/___NAMESPACE___/___FILEBASENAME___.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace ___NAMESPACE___;
+
+/******************************************************************************
+*                  T___FILEBASENAME___ implementation                             *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+T___FILEBASENAME___::T___FILEBASENAME___(const std::string name)
+: Module<___FILEBASENAME___Par>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+std::vector<std::string> T___FILEBASENAME___::getInput(void)
+{
+    std::vector<std::string> in;
+    
+    return in;
+}
+
+std::vector<std::string> T___FILEBASENAME___::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+void T___FILEBASENAME___::setup(void)
+{
+
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+void T___FILEBASENAME___::execute(void)
+{
+
+}

extras/Hadrons/Modules/templates/Module_in_NS.hpp.template

@@ -0,0 +1,44 @@
+#ifndef Hadrons____FILEBASENAME____hpp_
+#define Hadrons____FILEBASENAME____hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         ___FILEBASENAME___                                 *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(___NAMESPACE___)
+
+class ___FILEBASENAME___Par: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(___FILEBASENAME___Par,
+                                    unsigned int, i);
+};
+
+class T___FILEBASENAME___: public Module<___FILEBASENAME___Par>
+{
+public:
+    // constructor
+    T___FILEBASENAME___(const std::string name);
+    // destructor
+    virtual ~T___FILEBASENAME___(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(___FILEBASENAME___, T___FILEBASENAME___, ___NAMESPACE___);
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons____FILEBASENAME____hpp_

extras/Hadrons/Modules/templates/Module_tmp.hpp.template

@@ -0,0 +1,81 @@
+#ifndef Hadrons____FILEBASENAME____hpp_
+#define Hadrons____FILEBASENAME____hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         ___FILEBASENAME___                                 *
+ ******************************************************************************/
+class ___FILEBASENAME___Par: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(___FILEBASENAME___Par,
+                                    unsigned int, i);
+};
+
+template <typename FImpl>
+class T___FILEBASENAME___: public Module<___FILEBASENAME___Par>
+{
+public:
+    // constructor
+    T___FILEBASENAME___(const std::string name);
+    // destructor
+    virtual ~T___FILEBASENAME___(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER(___FILEBASENAME___, T___FILEBASENAME___<FIMPL>);
+
+/******************************************************************************
+ *                 T___FILEBASENAME___ implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+T___FILEBASENAME___<FImpl>::T___FILEBASENAME___(const std::string name)
+: Module<___FILEBASENAME___Par>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> T___FILEBASENAME___<FImpl>::getInput(void)
+{
+    std::vector<std::string> in;
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> T___FILEBASENAME___<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void T___FILEBASENAME___<FImpl>::setup(void)
+{
+    
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void T___FILEBASENAME___<FImpl>::execute(void)
+{
+    
+}
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons____FILEBASENAME____hpp_

extras/Hadrons/Modules/templates/Module_tmp_in_NS.hpp.template

@@ -0,0 +1,85 @@
+#ifndef Hadrons____FILEBASENAME____hpp_
+#define Hadrons____FILEBASENAME____hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         ___FILEBASENAME___                                 *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(___NAMESPACE___)
+
+class ___FILEBASENAME___Par: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(___FILEBASENAME___Par,
+                                    unsigned int, i);
+};
+
+template <typename FImpl>
+class T___FILEBASENAME___: public Module<___FILEBASENAME___Par>
+{
+public:
+    // constructor
+    T___FILEBASENAME___(const std::string name);
+    // destructor
+    virtual ~T___FILEBASENAME___(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(___FILEBASENAME___, T___FILEBASENAME___<FIMPL>, ___NAMESPACE___);
+
+/******************************************************************************
+ *                 T___FILEBASENAME___ implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+T___FILEBASENAME___<FImpl>::T___FILEBASENAME___(const std::string name)
+: Module<___FILEBASENAME___Par>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> T___FILEBASENAME___<FImpl>::getInput(void)
+{
+    std::vector<std::string> in;
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> T___FILEBASENAME___<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void T___FILEBASENAME___<FImpl>::setup(void)
+{
+    
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void T___FILEBASENAME___<FImpl>::execute(void)
+{
+    
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons____FILEBASENAME____hpp_

extras/Hadrons/add_module.sh

@@ -0,0 +1,31 @@
+#!/usr/bin/env bash
+
+if (( $# != 1 && $# != 2)); then
+    echo "usage: `basename $0` <module name> [<namespace>]" 1>&2
+    exit 1
+fi
+NAME=$1
+NS=$2
+
+if (( $# == 1 )); then
+	if [ -e "Modules/${NAME}.cc" ] || [ -e "Modules/${NAME}.hpp" ]; then
+	    echo "error: files Modules/${NAME}.* already exists" 1>&2
+	    exit 1
+	fi
+	sed "s/___FILEBASENAME___/${NAME}/g" Modules/templates/Module.cc.template > Modules/${NAME}.cc
+	sed "s/___FILEBASENAME___/${NAME}/g" Modules/templates/Module.hpp.template > Modules/${NAME}.hpp
+elif (( $# == 2 )); then
+	mkdir -p Modules/${NS}
+	if [ -e "Modules/${NS}/${NAME}.cc" ] || [ -e "Modules/${NS}/${NAME}.hpp" ]; then
+	    echo "error: files Modules/${NS}/${NAME}.* already exists" 1>&2
+	    exit 1
+	fi
+	TMPCC=".${NS}.${NAME}.tmp.cc"
+	TMPHPP=".${NS}.${NAME}.tmp.hpp"
+	sed "s/___FILEBASENAME___/${NAME}/g" Modules/templates/Module_in_NS.cc.template  > ${TMPCC}
+	sed "s/___FILEBASENAME___/${NAME}/g" Modules/templates/Module_in_NS.hpp.template > ${TMPHPP}
+	sed "s/___NAMESPACE___/${NS}/g" ${TMPCC}  > Modules/${NS}/${NAME}.cc
+	sed "s/___NAMESPACE___/${NS}/g" ${TMPHPP} > Modules/${NS}/${NAME}.hpp
+	rm -f ${TMPCC} ${TMPHPP}
+fi
+./make_module_list.sh

extras/Hadrons/add_module_template.sh

@@ -0,0 +1,28 @@
+#!/usr/bin/env bash
+
+if (( $# != 1 && $# != 2)); then
+    echo "usage: `basename $0` <module name> [<namespace>]" 1>&2
+    exit 1
+fi
+NAME=$1
+NS=$2
+
+if (( $# == 1 )); then
+	if [ -e "Modules/${NAME}.cc" ] || [ -e "Modules/${NAME}.hpp" ]; then
+	    echo "error: files Modules/${NAME}.* already exists" 1>&2
+	    exit 1
+	fi
+	sed "s/___FILEBASENAME___/${NAME}/g" Modules/templates/Module_tmp.hpp.template > Modules/${NAME}.hpp
+elif (( $# == 2 )); then
+	mkdir -p Modules/${NS}
+	if [ -e "Modules/${NS}/${NAME}.cc" ] || [ -e "Modules/${NS}/${NAME}.hpp" ]; then
+	    echo "error: files Modules/${NS}/${NAME}.* already exists" 1>&2
+	    exit 1
+	fi
+	TMPCC=".${NS}.${NAME}.tmp.cc"
+	TMPHPP=".${NS}.${NAME}.tmp.hpp"
+	sed "s/___FILEBASENAME___/${NAME}/g" Modules/templates/Module_tmp_in_NS.hpp.template > ${TMPHPP}
+	sed "s/___NAMESPACE___/${NS}/g" ${TMPHPP} > Modules/${NS}/${NAME}.hpp
+	rm -f ${TMPCC} ${TMPHPP}
+fi
+./make_module_list.sh

extras/Hadrons/make_module_list.sh

@@ -0,0 +1,12 @@
+#!/usr/bin/env bash
+
+echo 'modules_cc =\' > modules.inc
+find Modules -name '*.cc' -type f -print | sed 's/^/  /;$q;s/$/ \\/' >> modules.inc
+echo '' >> modules.inc
+echo 'modules_hpp =\' >> modules.inc
+find Modules -name '*.hpp' -type f -print | sed 's/^/  /;$q;s/$/ \\/' >> modules.inc
+echo '' >> modules.inc
+rm -f Modules.hpp
+for f in `find Modules -name '*.hpp'`; do
+	echo "#include <Grid/Hadrons/${f}>" >> Modules.hpp
+done

extras/Hadrons/modules.inc

@@ -0,0 +1,19 @@
+modules_cc =\
+  Modules/MGauge/Load.cc \
+  Modules/MGauge/Random.cc \
+  Modules/MGauge/Unit.cc
+
+modules_hpp =\
+  Modules/MAction/DWF.hpp \
+  Modules/MAction/Wilson.hpp \
+  Modules/MContraction/Baryon.hpp \
+  Modules/MContraction/Meson.hpp \
+  Modules/MGauge/Load.hpp \
+  Modules/MGauge/Random.hpp \
+  Modules/MGauge/Unit.hpp \
+  Modules/MSolver/RBPrecCG.hpp \
+  Modules/MSource/Point.hpp \
+  Modules/MSource/SeqGamma.hpp \
+  Modules/MSource/Z2.hpp \
+  Modules/Quark.hpp
+

lib/FFT.h

@@ -244,7 +244,10 @@ namespace Grid {
             pokeLocalSite(s,pgbuf,cbuf);
           }
         }
-        result = Cshift(result,dim,L);
+        if (p != processors[dim] - 1)
+        {
+          result = Cshift(result,dim,L);
+        }
       }
       
       // Loop over orthog coords
@@ -287,10 +290,10 @@ namespace Grid {
           cgbuf = clbuf;
           cgbuf[dim] = clbuf[dim]+L*pc;
           peekLocalSite(s,pgbuf,cgbuf);
-          s = s * div;
           pokeLocalSite(s,result,clbuf);
         }
       }
+      result = result*div;
       
       // destroying plan
       FFTW<scalar>::fftw_destroy_plan(p);

lib/Hadrons

@@ -0,0 +1 @@
+../extras/Hadrons

lib/Log.cc

@@ -93,7 +93,7 @@ void GridLogConfigure(std::vector<std::string> &logstreams) {
 ////////////////////////////////////////////////////////////
 void Grid_quiesce_nodes(void) {
   int me = 0;
-#if defined(GRID_COMMS_MPI) || defined(GRID_COMMS_MPI3)
+#if defined(GRID_COMMS_MPI) || defined(GRID_COMMS_MPI3) || defined(GRID_COMMS_MPI3L)
   MPI_Comm_rank(MPI_COMM_WORLD, &me);
 #endif
 #ifdef GRID_COMMS_SHMEM

lib/Log.h

@@ -110,8 +110,8 @@ public:
   friend std::ostream& operator<< (std::ostream& stream, Logger& log){
 
     if ( log.active ) {
-      stream << log.background()<< log.topName << log.background()<< " : ";
+      stream << log.background()<< std::setw(10) << std::left << log.topName << log.background()<< " : ";
-      stream << log.colour() <<std::setw(14) << std::left << log.name << log.background() << " : ";
+      stream << log.colour() << std::setw(14) << std::left << log.name << log.background() << " : ";
       if ( log.timestamp ) {
 	StopWatch.Stop();
 	GridTime now = StopWatch.Elapsed();

lib/Threads.h

@@ -46,11 +46,13 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #endif
 #define PARALLEL_NESTED_LOOP2 _Pragma("omp parallel for collapse(2)")
 #define PARALLEL_REGION       _Pragma("omp parallel")
+#define PARALLEL_CRITICAL     _Pragma("omp critical")
 #else
 #define PARALLEL_FOR_LOOP
 #define PARALLEL_FOR_LOOP_INTERN
 #define PARALLEL_NESTED_LOOP2
 #define PARALLEL_REGION
+#define PARALLEL_CRITICAL
 #endif
 
 namespace Grid {

lib/algorithms/iterative/ConjugateGradient.h

@@ -154,7 +154,7 @@ class ConjugateGradient : public OperatorFunction<Field> {
                   << LinalgTimer.Elapsed();
         std::cout << std::endl;
 
-        if (ErrorOnNoConverge) assert(true_residual / Tolerance < 1000.0);
+        if (ErrorOnNoConverge) assert(true_residual / Tolerance < 10000.0);
 
         return;
       }

lib/algorithms/iterative/ImplicitlyRestartedLanczos.h

@@ -1080,10 +1080,10 @@ say con = 2
 **/
 
 template<class T>
-static void Lock(DenseMatrix<T> &H, 	///Hess mtx	
+static void Lock(DenseMatrix<T> &H, 	// Hess mtx	
-		 DenseMatrix<T> &Q, 	///Lock Transform
+		 DenseMatrix<T> &Q, 	// Lock Transform
-		 T val, 		///value to be locked
+		 T val, 		// value to be locked
-		 int con, 	///number already locked
+		 int con, 	// number already locked
 		 RealD small,
 		 int dfg,
 		 bool herm)

lib/qcd/action/Actions.h

@@ -195,6 +195,7 @@ typedef WilsonTMFermion<WilsonImplD> WilsonTMFermionD;
 typedef DomainWallFermion<WilsonImplR> DomainWallFermionR;
 typedef DomainWallFermion<WilsonImplF> DomainWallFermionF;
 typedef DomainWallFermion<WilsonImplD> DomainWallFermionD;
+
 typedef MobiusFermion<WilsonImplR> MobiusFermionR;
 typedef MobiusFermion<WilsonImplF> MobiusFermionF;
 typedef MobiusFermion<WilsonImplD> MobiusFermionD;
@@ -203,6 +204,20 @@ typedef ZMobiusFermion<ZWilsonImplR> ZMobiusFermionR;
 typedef ZMobiusFermion<ZWilsonImplF> ZMobiusFermionF;
 typedef ZMobiusFermion<ZWilsonImplD> ZMobiusFermionD;
 
+// Ls vectorised 
+typedef DomainWallFermion<DomainWallVec5dImplR> DomainWallFermionVec5dR;
+typedef DomainWallFermion<DomainWallVec5dImplF> DomainWallFermionVec5dF;
+typedef DomainWallFermion<DomainWallVec5dImplD> DomainWallFermionVec5dD;
+
+typedef MobiusFermion<DomainWallVec5dImplR> MobiusFermionVec5dR;
+typedef MobiusFermion<DomainWallVec5dImplF> MobiusFermionVec5dF;
+typedef MobiusFermion<DomainWallVec5dImplD> MobiusFermionVec5dD;
+
+typedef ZMobiusFermion<ZDomainWallVec5dImplR> ZMobiusFermionVec5dR;
+typedef ZMobiusFermion<ZDomainWallVec5dImplF> ZMobiusFermionVec5dF;
+typedef ZMobiusFermion<ZDomainWallVec5dImplD> ZMobiusFermionVec5dD;
+
+
 typedef ScaledShamirFermion<WilsonImplR> ScaledShamirFermionR;
 typedef ScaledShamirFermion<WilsonImplF> ScaledShamirFermionF;
 typedef ScaledShamirFermion<WilsonImplD> ScaledShamirFermionD;
@@ -254,6 +269,7 @@ typedef MobiusFermion<GparityWilsonImplF> GparityMobiusFermionF;
 typedef MobiusFermion<GparityWilsonImplD> GparityMobiusFermionD;
 
 
+
   }}
 ///////////////////////////////////////////////////////////////////////////////
 // G5 herm -- this has to live in QCD since dirac matrix is not in the broader sector of code

lib/qcd/action/fermion/CayleyFermion5D.cc

@@ -62,6 +62,50 @@ void CayleyFermion5D<Impl>::Dminus(const FermionField &psi, FermionField &chi)
     axpby_ssp(chi,Coeff_t(1.0),psi,-cs[s],tmp,s,s);// chi = (1-c[s] D_W) psi
   }
 }
+
+
+template<class Impl> void CayleyFermion5D<Impl>::CayleyReport(void)
+{
+  this->Report();
+  std::vector<int> latt = GridDefaultLatt();          
+  RealD volume = this->Ls;  for(int mu=0;mu<Nd;mu++) volume=volume*latt[mu];
+  RealD NP     = this->_FourDimGrid->_Nprocessors;
+  if ( M5Dcalls > 0 ) {
+    std::cout << GridLogMessage << "#### M5D calls report " << std::endl;
+    std::cout << GridLogMessage << "CayleyFermion5D Number of M5D Calls     : " << M5Dcalls   << std::endl;
+    std::cout << GridLogMessage << "CayleyFermion5D ComputeTime/Calls       : " << M5Dtime / M5Dcalls << " us" << std::endl;
+
+    // Flops = 6.0*(Nc*Ns) *Ls*vol
+    RealD mflops = 6.0*12*volume*M5Dcalls/M5Dtime/2; // 2 for red black counting
+    std::cout << GridLogMessage << "Average mflops/s per call                : " << mflops << std::endl;
+    std::cout << GridLogMessage << "Average mflops/s per call per rank       : " << mflops/NP << std::endl;
+  }
+
+  if ( MooeeInvCalls > 0 ) {
+
+    std::cout << GridLogMessage << "#### MooeeInv calls report " << std::endl;
+    std::cout << GridLogMessage << "CayleyFermion5D Number of MooeeInv Calls     : " << MooeeInvCalls   << std::endl;
+    std::cout << GridLogMessage << "CayleyFermion5D ComputeTime/Calls            : " << MooeeInvTime / MooeeInvCalls << " us" << std::endl;
+
+    // Flops = 9*12*Ls*vol/2
+    RealD mflops = 9.0*12*volume*MooeeInvCalls/MooeeInvTime/2; // 2 for red black counting
+    std::cout << GridLogMessage << "Average mflops/s per call                : " << mflops << std::endl;
+    std::cout << GridLogMessage << "Average mflops/s per call per rank       : " << mflops/NP << std::endl;
+  }
+
+}
+template<class Impl> void CayleyFermion5D<Impl>::CayleyZeroCounters(void)
+{
+  this->ZeroCounters();
+  M5Dflops=0;
+  M5Dcalls=0;
+  M5Dtime=0;
+  MooeeInvFlops=0;
+  MooeeInvCalls=0;
+  MooeeInvTime=0;
+}
+
+
 template<class Impl>  
 void CayleyFermion5D<Impl>::DminusDag(const FermionField &psi, FermionField &chi)
 {

lib/qcd/action/fermion/CayleyFermion5D.h

@@ -121,6 +121,18 @@ namespace Grid {
 		      RealD _mass,RealD _M5,const ImplParams &p= ImplParams());
 
       
+
+     void CayleyReport(void);
+     void CayleyZeroCounters(void);
+
+     double M5Dflops;
+     double M5Dcalls;
+     double M5Dtime;
+
+     double MooeeInvFlops;
+     double MooeeInvCalls;
+     double MooeeInvTime;
+
     protected:
       void SetCoefficientsZolotarev(RealD zolohi,Approx::zolotarev_data *zdata,RealD b,RealD c);
       void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c);

lib/qcd/action/fermion/CayleyFermion5Dcache.cc

@@ -51,6 +51,9 @@ void CayleyFermion5D<Impl>::M5D(const FermionField &psi,
   GridBase *grid=psi._grid;
   assert(phi.checkerboard == psi.checkerboard);
   chi.checkerboard=psi.checkerboard;
+  // Flops = 6.0*(Nc*Ns) *Ls*vol
+  M5Dcalls++;
+  M5Dtime-=usecond();
 PARALLEL_FOR_LOOP
   for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
     for(int s=0;s<Ls;s++){
@@ -76,6 +79,7 @@ PARALLEL_FOR_LOOP
       }
     }
   }
+  M5Dtime+=usecond();
 }
 
 template<class Impl>  
@@ -91,6 +95,9 @@ void CayleyFermion5D<Impl>::M5Ddag(const FermionField &psi,
   assert(phi.checkerboard == psi.checkerboard);
   chi.checkerboard=psi.checkerboard;
 
+  // Flops = 6.0*(Nc*Ns) *Ls*vol
+  M5Dcalls++;
+  M5Dtime-=usecond();
 PARALLEL_FOR_LOOP
   for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
     auto tmp = psi._odata[0];
@@ -116,6 +123,7 @@ PARALLEL_FOR_LOOP
       }
     }
   }
+  M5Dtime+=usecond();
 }
 
 template<class Impl>
@@ -126,10 +134,14 @@ void CayleyFermion5D<Impl>::MooeeInv    (const FermionField &psi, FermionField &
 
   chi.checkerboard=psi.checkerboard;
 
+  MooeeInvCalls++;
+  MooeeInvTime-=usecond();
+
 PARALLEL_FOR_LOOP
   for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
     auto tmp = psi._odata[0];
 
+    // flops = 12*2*Ls + 12*2*Ls + 3*12*Ls + 12*2*Ls  = 12*Ls * (9) = 108*Ls flops
     // Apply (L^{\prime})^{-1}
     chi[ss]=psi[ss]; // chi[0]=psi[0]
     for(int s=1;s<Ls;s++){
@@ -155,6 +167,9 @@ PARALLEL_FOR_LOOP
       chi[ss+s] = chi[ss+s] - uee[s]*tmp;
     }
   }
+
+  MooeeInvTime+=usecond();
+
 }
 
 template<class Impl>
@@ -166,6 +181,8 @@ void CayleyFermion5D<Impl>::MooeeInvDag (const FermionField &psi, FermionField &
   assert(psi.checkerboard == psi.checkerboard);
   chi.checkerboard=psi.checkerboard;
 
+  MooeeInvCalls++;
+  MooeeInvTime-=usecond();
 
 PARALLEL_FOR_LOOP
   for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
@@ -197,6 +214,9 @@ PARALLEL_FOR_LOOP
       chi[ss+s] = chi[ss+s] - lee[s]*tmp;
     }
   }
+
+  MooeeInvTime+=usecond();
+
 }
 
 #ifdef CAYLEY_DPERP_CACHE

lib/qcd/action/fermion/CayleyFermion5Dvec.cc

@@ -60,7 +60,7 @@ void CayleyFermion5D<Impl>::M5D(const FermionField &psi,
   GridBase *grid=psi._grid;
   int Ls   = this->Ls;
   int LLs  = grid->_rdimensions[0];
-  int nsimd= Simd::Nsimd();
+  const int nsimd= Simd::Nsimd();
 
   Vector<iSinglet<Simd> > u(LLs);
   Vector<iSinglet<Simd> > l(LLs);
@@ -86,35 +86,138 @@ void CayleyFermion5D<Impl>::M5D(const FermionField &psi,
     d_p[ss] = diag[s];
   }}
 
+
+  M5Dcalls++;
+  M5Dtime-=usecond();
+
+  assert(Nc==3);
+
 PARALLEL_FOR_LOOP
   for(int ss=0;ss<grid->oSites();ss+=LLs){ // adds LLs
+#if 0
+      alignas(64) SiteHalfSpinor hp;
+      alignas(64) SiteHalfSpinor hm;
+      alignas(64) SiteSpinor fp;
+      alignas(64) SiteSpinor fm;
 
-    alignas(64) SiteHalfSpinor hp;
+      for(int v=0;v<LLs;v++){
-    alignas(64) SiteHalfSpinor hm;
-    alignas(64) SiteSpinor fp;
-    alignas(64) SiteSpinor fm;
 
-    for(int v=0;v<LLs;v++){
+	int vp=(v+1)%LLs;
+	int vm=(v+LLs-1)%LLs;
 
-      int vp=(v+1)%LLs;
+	spProj5m(hp,psi[ss+vp]);
-      int vm=(v+LLs-1)%LLs;
+	spProj5p(hm,psi[ss+vm]);
 
-      spProj5m(hp,psi[ss+vp]);
+	if ( vp<=v ) rotate(hp,hp,1);
-      spProj5p(hm,psi[ss+vm]);
+	if ( vm>=v ) rotate(hm,hm,nsimd-1);
 	
-      if ( vp<=v ) rotate(hp,hp,1);
+	hp=0.5*hp;
-      if ( vm>=v ) rotate(hm,hm,nsimd-1);
+        hm=0.5*hm;
 
-      hp=hp*0.5;
+	spRecon5m(fp,hp);
-      hm=hm*0.5;
+	spRecon5p(fm,hm);
-      spRecon5m(fp,hp);
-      spRecon5p(fm,hm);
 
-      chi[ss+v] = d[v]*phi[ss+v]+u[v]*fp;
+	chi[ss+v] = d[v]*phi[ss+v];
-      chi[ss+v] = chi[ss+v]     +l[v]*fm;
+	chi[ss+v] = chi[ss+v]     +u[v]*fp;
+	chi[ss+v] = chi[ss+v]     +l[v]*fm;
 
-    }
+      }
+#else
+      for(int v=0;v<LLs;v++){
+
+	vprefetch(psi[ss+v+LLs]);
+	//	vprefetch(phi[ss+v+LLs]);
+
+	int vp= (v==LLs-1) ? 0     : v+1;
+	int vm= (v==0    ) ? LLs-1 : v-1;
+	
+	Simd hp_00 = psi[ss+vp]()(2)(0); 
+	Simd hp_01 = psi[ss+vp]()(2)(1); 
+	Simd hp_02 = psi[ss+vp]()(2)(2); 
+	Simd hp_10 = psi[ss+vp]()(3)(0); 
+	Simd hp_11 = psi[ss+vp]()(3)(1); 
+	Simd hp_12 = psi[ss+vp]()(3)(2); 
+	
+	Simd hm_00 = psi[ss+vm]()(0)(0); 
+	Simd hm_01 = psi[ss+vm]()(0)(1); 
+	Simd hm_02 = psi[ss+vm]()(0)(2); 
+	Simd hm_10 = psi[ss+vm]()(1)(0); 
+	Simd hm_11 = psi[ss+vm]()(1)(1); 
+	Simd hm_12 = psi[ss+vm]()(1)(2); 
+
+	//	if ( ss==0) std::cout << " hp_00 " <<hp_00<<std::endl;
+	//	if ( ss==0) std::cout << " hm_00 " <<hm_00<<std::endl;
+
+	if ( vp<=v ) {
+	  hp_00.v = Optimization::Rotate::tRotate<2>(hp_00.v);
+	  hp_01.v = Optimization::Rotate::tRotate<2>(hp_01.v);
+	  hp_02.v = Optimization::Rotate::tRotate<2>(hp_02.v);
+	  hp_10.v = Optimization::Rotate::tRotate<2>(hp_10.v);
+	  hp_11.v = Optimization::Rotate::tRotate<2>(hp_11.v);
+	  hp_12.v = Optimization::Rotate::tRotate<2>(hp_12.v);
+	}
+	if ( vm>=v ) {
+	  hm_00.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_00.v);
+	  hm_01.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_01.v);
+	  hm_02.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_02.v);
+	  hm_10.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_10.v);
+	  hm_11.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_11.v);
+	  hm_12.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_12.v);
+	}
+
+	/*
+	if ( ss==0) std::cout << " dphi_00 " <<d[v]()()() * phi[ss+v]()(0)(0) <<std::endl;
+	if ( ss==0) std::cout << " dphi_10 " <<d[v]()()() * phi[ss+v]()(1)(0) <<std::endl;
+	if ( ss==0) std::cout << " dphi_20 " <<d[v]()()() * phi[ss+v]()(2)(0) <<std::endl;
+	if ( ss==0) std::cout << " dphi_30 " <<d[v]()()() * phi[ss+v]()(3)(0) <<std::endl;
+	*/	
+	Simd p_00  = d[v]()()() * phi[ss+v]()(0)(0)  + l[v]()()()*hm_00; 
+	Simd p_01  = d[v]()()() * phi[ss+v]()(0)(1)  + l[v]()()()*hm_01; 
+	Simd p_02  = d[v]()()() * phi[ss+v]()(0)(2)  + l[v]()()()*hm_02; 
+	Simd p_10  = d[v]()()() * phi[ss+v]()(1)(0)  + l[v]()()()*hm_10; 
+	Simd p_11  = d[v]()()() * phi[ss+v]()(1)(1)  + l[v]()()()*hm_11; 
+	Simd p_12  = d[v]()()() * phi[ss+v]()(1)(2)  + l[v]()()()*hm_12; 
+	Simd p_20  = d[v]()()() * phi[ss+v]()(2)(0)  + u[v]()()()*hp_00; 
+	Simd p_21  = d[v]()()() * phi[ss+v]()(2)(1)  + u[v]()()()*hp_01; 
+	Simd p_22  = d[v]()()() * phi[ss+v]()(2)(2)  + u[v]()()()*hp_02;  
+	Simd p_30  = d[v]()()() * phi[ss+v]()(3)(0)  + u[v]()()()*hp_10; 
+	Simd p_31  = d[v]()()() * phi[ss+v]()(3)(1)  + u[v]()()()*hp_11; 
+	Simd p_32  = d[v]()()() * phi[ss+v]()(3)(2)  + u[v]()()()*hp_12; 
+
+	
+	//	if ( ss==0){
+	/*
+	std::cout << ss<<" "<< v<< " good "<< chi[ss+v]()(0)(0) << " bad "<<p_00<<" diff "<<chi[ss+v]()(0)(0)-p_00<<std::endl;
+	std::cout << ss<<" "<< v<< " good "<< chi[ss+v]()(0)(1) << " bad "<<p_01<<" diff "<<chi[ss+v]()(0)(1)-p_01<<std::endl;
+	std::cout << ss<<" "<< v<< " good "<< chi[ss+v]()(0)(2) << " bad "<<p_02<<" diff "<<chi[ss+v]()(0)(2)-p_02<<std::endl;
+	std::cout << ss<<" "<< v<< " good "<< chi[ss+v]()(1)(0) << " bad "<<p_10<<" diff "<<chi[ss+v]()(1)(0)-p_10<<std::endl;
+	std::cout << ss<<" "<< v<< " good "<< chi[ss+v]()(1)(1) << " bad "<<p_11<<" diff "<<chi[ss+v]()(1)(1)-p_11<<std::endl;
+	std::cout << ss<<" "<< v<< " good "<< chi[ss+v]()(1)(2) << " bad "<<p_12<<" diff "<<chi[ss+v]()(1)(2)-p_12<<std::endl;
+	std::cout << ss<<" "<< v<< " good "<< chi[ss+v]()(2)(0) << " bad "<<p_20<<" diff "<<chi[ss+v]()(2)(0)-p_20<<std::endl;
+	std::cout << ss<<" "<< v<< " good "<< chi[ss+v]()(2)(1) << " bad "<<p_21<<" diff "<<chi[ss+v]()(2)(1)-p_21<<std::endl;
+	std::cout << ss<<" "<< v<< " good "<< chi[ss+v]()(2)(2) << " bad "<<p_22<<" diff "<<chi[ss+v]()(2)(2)-p_22<<std::endl;
+	std::cout << ss<<" "<< v<< " good "<< chi[ss+v]()(3)(0) << " bad "<<p_30<<" diff "<<chi[ss+v]()(3)(0)-p_30<<std::endl;
+	std::cout << ss<<" "<< v<< " good "<< chi[ss+v]()(3)(1) << " bad "<<p_31<<" diff "<<chi[ss+v]()(3)(1)-p_31<<std::endl;
+	std::cout << ss<<" "<< v<< " good "<< chi[ss+v]()(3)(2) << " bad "<<p_32<<" diff "<<chi[ss+v]()(3)(2)-p_32<<std::endl;
+	}
+	*/
+	vstream(chi[ss+v]()(0)(0),p_00);
+	vstream(chi[ss+v]()(0)(1),p_01);
+	vstream(chi[ss+v]()(0)(2),p_02);
+	vstream(chi[ss+v]()(1)(0),p_10);
+	vstream(chi[ss+v]()(1)(1),p_11);
+	vstream(chi[ss+v]()(1)(2),p_12);
+	vstream(chi[ss+v]()(2)(0),p_20);
+	vstream(chi[ss+v]()(2)(1),p_21);
+	vstream(chi[ss+v]()(2)(2),p_22);
+	vstream(chi[ss+v]()(3)(0),p_30);
+	vstream(chi[ss+v]()(3)(1),p_31);
+	vstream(chi[ss+v]()(3)(2),p_32);
+
+      }
+#endif
   }
+  M5Dtime+=usecond();
 }
 
 template<class Impl>  
@@ -154,6 +257,8 @@ void CayleyFermion5D<Impl>::M5Ddag(const FermionField &psi,
     d_p[ss] = diag[s];
   }}
 
+  M5Dcalls++;
+  M5Dtime-=usecond();
 PARALLEL_FOR_LOOP
   for(int ss=0;ss<grid->oSites();ss+=LLs){ // adds LLs
 
@@ -183,8 +288,8 @@ PARALLEL_FOR_LOOP
 
     }
   }
+  M5Dtime+=usecond();
 }
-
 template<class Impl>
 void CayleyFermion5D<Impl>::MooeeInternal(const FermionField &psi, FermionField &chi,int dag, int inv)
 {
@@ -250,13 +355,11 @@ void CayleyFermion5D<Impl>::MooeeInternal(const FermionField &psi, FermionField
     }
   }
   
+  MooeeInvCalls++;
+  MooeeInvTime-=usecond();
   // Dynamic allocate on stack to get per thread without serialised heap acces
-PARALLEL_FOR_LOOP
+#pragma omp parallel  
-  for(auto site=0;site<vol;site++){
+  {
-    
-    //    SiteHalfSpinor *SitePplus =(SiteHalfSpinor *) alloca(LLs*sizeof(SiteHalfSpinor));
-    //    SiteHalfSpinor *SitePminus=(SiteHalfSpinor *) alloca(LLs*sizeof(SiteHalfSpinor));
-    //    SiteSpinor     *SiteChi   =(SiteSpinor *)     alloca(LLs*sizeof(SiteSpinor));
 
     Vector<SiteHalfSpinor> SitePplus(LLs);
     Vector<SiteHalfSpinor> SitePminus(LLs);
@@ -267,6 +370,9 @@ PARALLEL_FOR_LOOP
     SiteHalfSpinor BcastP;
     SiteHalfSpinor BcastM;
 
+#pragma omp for 
+  for(auto site=0;site<vol;site++){
+
     for(int s=0;s<LLs;s++){
       int lex = s+LLs*site;
       spProj5p(SitePplus[s] ,psi[lex]);
@@ -294,6 +400,8 @@ PARALLEL_FOR_LOOP
       chi[lex] = SiteChi[s]*0.5;
     }
   }
+  }
+  MooeeInvTime+=usecond();
 }
 
 INSTANTIATE_DPERP(DomainWallVec5dImplD);

lib/qcd/action/fermion/FermionOperatorImpl.h

@@ -48,8 +48,10 @@ namespace QCD {
   //    typedef typename XXX         GaugeField;
   //    typedef typename XXX      GaugeActField;
   //    typedef typename XXX       FermionField;
+  //    typedef typename XXX    PropagatorField;
   //    typedef typename XXX  DoubledGaugeField;
   //    typedef typename XXX         SiteSpinor;
+  //    typedef typename XXX     SitePropagator;
   //    typedef typename XXX     SiteHalfSpinor;	
   //    typedef typename XXX         Compressor;	
   //
@@ -95,13 +97,15 @@ namespace QCD {
   
 #define INHERIT_FIMPL_TYPES(Impl)\
   typedef typename Impl::FermionField           FermionField;		\
+  typedef typename Impl::PropagatorField     PropagatorField;		\
   typedef typename Impl::DoubledGaugeField DoubledGaugeField;		\
   typedef typename Impl::SiteSpinor               SiteSpinor;		\
+  typedef typename Impl::SitePropagator       SitePropagator;		\
   typedef typename Impl::SiteHalfSpinor       SiteHalfSpinor;		\
   typedef typename Impl::Compressor               Compressor;		\
   typedef typename Impl::StencilImpl             StencilImpl;		\
-  typedef typename Impl::ImplParams ImplParams;				\
+  typedef typename Impl::ImplParams               ImplParams;	        \
-  typedef typename Impl::Coeff_t       Coeff_t;
+  typedef typename Impl::Coeff_t                     Coeff_t;           \
   
 #define INHERIT_IMPL_TYPES(Base) \
   INHERIT_GIMPL_TYPES(Base)	 \
@@ -127,14 +131,17 @@ namespace QCD {
     INHERIT_GIMPL_TYPES(Gimpl);
       
     template <typename vtype> using iImplSpinor            = iScalar<iVector<iVector<vtype, Dimension>, Ns> >;
+    template <typename vtype> using iImplPropagator        = iScalar<iMatrix<iMatrix<vtype, Dimension>, Ns> >;
     template <typename vtype> using iImplHalfSpinor        = iScalar<iVector<iVector<vtype, Dimension>, Nhs> >;
     template <typename vtype> using iImplDoubledGaugeField = iVector<iScalar<iMatrix<vtype, Dimension> >, Nds>;
     
     typedef iImplSpinor<Simd>            SiteSpinor;
+    typedef iImplPropagator<Simd>        SitePropagator;
     typedef iImplHalfSpinor<Simd>        SiteHalfSpinor;
     typedef iImplDoubledGaugeField<Simd> SiteDoubledGaugeField;
     
     typedef Lattice<SiteSpinor>            FermionField;
+    typedef Lattice<SitePropagator>        PropagatorField;
     typedef Lattice<SiteDoubledGaugeField> DoubledGaugeField;
     
     typedef WilsonCompressor<SiteHalfSpinor, SiteSpinor> Compressor;
@@ -216,14 +223,17 @@ class DomainWallVec5dImpl :  public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres
   INHERIT_GIMPL_TYPES(Gimpl);
   
   template <typename vtype> using iImplSpinor            = iScalar<iVector<iVector<vtype, Nrepresentation>, Ns> >;
+  template <typename vtype> using iImplPropagator        = iScalar<iMatrix<iMatrix<vtype, Nrepresentation>, Ns> >;
   template <typename vtype> using iImplHalfSpinor        = iScalar<iVector<iVector<vtype, Nrepresentation>, Nhs> >;
   template <typename vtype> using iImplDoubledGaugeField = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nds>;
   template <typename vtype> using iImplGaugeField        = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nd>;
   template <typename vtype> using iImplGaugeLink         = iScalar<iScalar<iMatrix<vtype, Nrepresentation> > >;
   
   typedef iImplSpinor<Simd> SiteSpinor;
+  typedef iImplPropagator<Simd> SitePropagator;
   typedef iImplHalfSpinor<Simd> SiteHalfSpinor;
   typedef Lattice<SiteSpinor> FermionField;
+  typedef Lattice<SitePropagator> PropagatorField;
   
   // Make the doubled gauge field a *scalar*
   typedef iImplDoubledGaugeField<typename Simd::scalar_type>  SiteDoubledGaugeField;  // This is a scalar
@@ -315,14 +325,17 @@ class GparityWilsonImpl : public ConjugateGaugeImpl<GaugeImplTypes<S, Nrepresent
  INHERIT_GIMPL_TYPES(Gimpl);
       
  template <typename vtype> using iImplSpinor            = iVector<iVector<iVector<vtype, Nrepresentation>, Ns>, Ngp>;
+ template <typename vtype> using iImplPropagator        = iVector<iMatrix<iMatrix<vtype, Nrepresentation>, Ns>, Ngp >;
  template <typename vtype> using iImplHalfSpinor        = iVector<iVector<iVector<vtype, Nrepresentation>, Nhs>, Ngp>;
  template <typename vtype> using iImplDoubledGaugeField = iVector<iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nds>, Ngp>;
       
  typedef iImplSpinor<Simd> SiteSpinor;
+ typedef iImplPropagator<Simd> SitePropagator;
  typedef iImplHalfSpinor<Simd> SiteHalfSpinor;
  typedef iImplDoubledGaugeField<Simd> SiteDoubledGaugeField;
  
  typedef Lattice<SiteSpinor> FermionField;
+ typedef Lattice<SitePropagator> PropagatorField;
  typedef Lattice<SiteDoubledGaugeField> DoubledGaugeField;
  
  typedef WilsonCompressor<SiteHalfSpinor, SiteSpinor> Compressor;

lib/qcd/action/fermion/WilsonFermion5D.cc

@@ -194,6 +194,11 @@ void WilsonFermion5D<Impl>::Report(void)
     std::cout << GridLogMessage << "Average mflops/s per call                : " << mflops << std::endl;
     std::cout << GridLogMessage << "Average mflops/s per call per rank       : " << mflops/NP << std::endl;
 
+    RealD Fullmflops = 1344*volume*DhopCalls/(DhopComputeTime+DhopCommTime)/2; // 2 for red black counting
+    std::cout << GridLogMessage << "Average mflops/s per call (full)         : " << Fullmflops << std::endl;
+    std::cout << GridLogMessage << "Average mflops/s per call per rank (full): " << Fullmflops/NP << std::endl;
+
+
    }
 
   if ( DerivCalls > 0 ) {
@@ -209,12 +214,15 @@ void WilsonFermion5D<Impl>::Report(void)
     RealD mflops = 144*volume*DerivCalls/DerivDhopComputeTime;
     std::cout << GridLogMessage << "Average mflops/s per call                : " << mflops << std::endl;
     std::cout << GridLogMessage << "Average mflops/s per call per node       : " << mflops/NP << std::endl;
-  }
+
+    RealD Fullmflops = 144*volume*DerivCalls/(DerivDhopComputeTime+DerivCommTime)/2; // 2 for red black counting
+    std::cout << GridLogMessage << "Average mflops/s per call (full)         : " << Fullmflops << std::endl;
+    std::cout << GridLogMessage << "Average mflops/s per call per node (full): " << Fullmflops/NP << std::endl;  }
 
   if (DerivCalls > 0 || DhopCalls > 0){
-    std::cout << GridLogMessage << "WilsonFermion5D Stencil"<<std::endl;  Stencil.Report();
+    std::cout << GridLogMessage << "WilsonFermion5D Stencil"    <<std::endl;  Stencil.Report();
     std::cout << GridLogMessage << "WilsonFermion5D StencilEven"<<std::endl;  StencilEven.Report();
-    std::cout << GridLogMessage << "WilsonFermion5D StencilOdd"<<std::endl;  StencilOdd.Report();
+    std::cout << GridLogMessage << "WilsonFermion5D StencilOdd" <<std::endl;  StencilOdd.Report();
   }
 }
 

lib/qcd/hmc/HmcRunner.h

@@ -116,7 +116,7 @@ class NerscHmcRunnerTemplate {
     NoSmearing<Gimpl> SmearingPolicy;
     typedef MinimumNorm2<GaugeField, NoSmearing<Gimpl>, RepresentationsPolicy >
         IntegratorType;  // change here to change the algorithm
-    IntegratorParameters MDpar(20, 1.0);
+    IntegratorParameters MDpar(40, 1.0);
     IntegratorType MDynamics(UGrid, MDpar, TheAction, SmearingPolicy);
 
     // Checkpoint strategy

lib/simd/Grid_avx.h

@@ -167,7 +167,7 @@ namespace Optimization {
     }
     //Integer
     inline __m256i operator()(__m256i a, __m256i b){
-#if defined (AVX1) || defined (AVXFMA4)
+#if defined (AVX1) || defined (AVXFMA) || defined (AVXFMA4)
           __m128i a0,a1;
           __m128i b0,b1;
           a0 = _mm256_extractf128_si256(a,0);
@@ -195,7 +195,7 @@ namespace Optimization {
     }
     //Integer
     inline __m256i operator()(__m256i a, __m256i b){
-#if defined (AVX1) || defined (AVXFMA4)
+#if defined (AVX1) || defined (AVXFMA) || defined (AVXFMA4)
           __m128i a0,a1;
           __m128i b0,b1;
           a0 = _mm256_extractf128_si256(a,0);
@@ -233,7 +233,7 @@ namespace Optimization {
       a_imag = _mm256_mul_ps( a_imag,tmp  );  // (Ai, Ai) * (Bi, Br) = Ai Bi, Ai Br
       return _mm256_maddsub_ps( a_real, b, a_imag ); // Ar Br , Ar Bi   +- Ai Bi             = ArBr-AiBi , ArBi+AiBr
 #endif
-#if defined (AVX2)
+#if defined (AVX2)  || defined (AVXFMA)
       __m256 a_real = _mm256_moveldup_ps( a ); // Ar Ar
       __m256 a_imag = _mm256_movehdup_ps( a ); // Ai Ai
       a_imag = _mm256_mul_ps( a_imag, _mm256_shuffle_ps( b,b, _MM_SELECT_FOUR_FOUR(2,3,0,1) ));  // (Ai, Ai) * (Bi, Br) = Ai Bi, Ai Br
@@ -279,7 +279,7 @@ namespace Optimization {
       a_imag = _mm256_mul_pd( a_imag, _mm256_permute_pd( b, 0x5 ) );  // (Ai, Ai) * (Bi, Br) = Ai Bi, Ai Br
       return _mm256_maddsub_pd( a_real, b, a_imag ); // Ar Br , Ar Bi   +- Ai Bi             = ArBr-AiBi , ArBi+AiBr
 #endif
-#if defined (AVX2)
+#if defined (AVX2) || defined (AVXFMA)
       __m256d a_real = _mm256_movedup_pd( a ); // Ar Ar
       __m256d a_imag = _mm256_shuffle_pd(a,a,0xF);//aiai
       a_imag = _mm256_mul_pd( a_imag, _mm256_permute_pd( b, 0x5 ) );  // (Ai, Ai) * (Bi, Br) = Ai Bi, Ai Br
@@ -320,7 +320,7 @@ namespace Optimization {
 #if defined (AVXFMA4)
       a= _mm256_macc_ps(b,c,a);
 #endif
-#if defined (AVX2)
+#if defined (AVX2) || defined (AVXFMA)
       a= _mm256_fmadd_ps( b, c, a);
 #endif
     }
@@ -332,7 +332,7 @@ namespace Optimization {
 #if defined (AVXFMA4)
       a= _mm256_macc_pd(b,c,a);
 #endif
-#if defined (AVX2)
+#if defined (AVX2) || defined (AVXFMA)
       a= _mm256_fmadd_pd( b, c, a);
 #endif
     }
@@ -347,7 +347,7 @@ namespace Optimization {
     }
     // Integer
     inline __m256i operator()(__m256i a, __m256i b){
-#if defined (AVX1) 
+#if defined (AVX1) || defined (AVXFMA)
       __m128i a0,a1;
       __m128i b0,b1;
       a0 = _mm256_extractf128_si256(a,0);

lib/simd/Grid_avx512.h

@@ -27,15 +27,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-//----------------------------------------------------------------------
-/*! @file Grid_knc.h
-  @brief Optimization libraries for AVX512 instructions set for KNC
-
-  Using intrinsics
-*/
-// Time-stamp: <2015-06-09 14:27:28 neo>
-//----------------------------------------------------------------------
-
 #include <immintrin.h>
 
 
@@ -95,13 +86,13 @@ namespace Optimization {
   struct Vstream{
     //Float
     inline void operator()(float * a, __m512 b){
-      //_mm512_stream_ps(a,b);
+      _mm512_stream_ps(a,b);
-      _mm512_store_ps(a,b);
+      //      _mm512_store_ps(a,b);
     }
     //Double
     inline void operator()(double * a, __m512d b){
-      //_mm512_stream_pd(a,b);
+      _mm512_stream_pd(a,b);
-      _mm512_store_pd(a,b);
+      //      _mm512_store_pd(a,b);
     }
 
   };

lib/simd/Grid_generic.h

@@ -6,8 +6,7 @@
 
     Copyright (C) 2015
 
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Antonin Portelli <antonin.portelli@me.com>
-Author: neo <cossu@post.kek.jp>
 
     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
@@ -27,133 +26,352 @@ Author: neo <cossu@post.kek.jp>
     *************************************************************************************/
     /*  END LEGAL */
 
+static_assert(GEN_SIMD_WIDTH % 16u == 0, "SIMD vector size is not an integer multiple of 16 bytes");
+
+//#define VECTOR_LOOPS
+
+// playing with compiler pragmas
+#ifdef VECTOR_LOOPS
+#ifdef __clang__
+#define VECTOR_FOR(i, w, inc)\
+_Pragma("clang loop unroll(full) vectorize(enable) interleave(enable) vectorize_width(w)")\
+for (unsigned int i = 0; i < w; i += inc)
+#elif defined __INTEL_COMPILER
+#define VECTOR_FOR(i, w, inc)\
+_Pragma("simd vectorlength(w*8)")\
+for (unsigned int i = 0; i < w; i += inc)
+#else
+#define VECTOR_FOR(i, w, inc)\
+for (unsigned int i = 0; i < w; i += inc)
+#endif
+#else
+#define VECTOR_FOR(i, w, inc)\
+for (unsigned int i = 0; i < w; i += inc)
+#endif
+
 namespace Grid {
 namespace Optimization {
 
-  template<class vtype>
+  // type traits giving the number of elements for each vector type
-  union uconv {
+  template <typename T> struct W;
-    float f;
+  template <> struct W<double> {
-    vtype v;
+    constexpr static unsigned int c = GEN_SIMD_WIDTH/16u;
+    constexpr static unsigned int r = GEN_SIMD_WIDTH/8u;
+  };
+  template <> struct W<float> {
+    constexpr static unsigned int c = GEN_SIMD_WIDTH/8u;
+    constexpr static unsigned int r = GEN_SIMD_WIDTH/4u;
   };
   
-  union u128f {
+  // SIMD vector types
-    float v;
+  template <typename T>
-    float f[4];
+  struct vec {
-  };
+    alignas(GEN_SIMD_WIDTH) T v[W<T>::r];
-  union u128d {
-    double v;
-    double f[2];
   };
   
+  typedef vec<float>   vecf;
+  typedef vec<double>  vecd;
+  
   struct Vsplat{
-    //Complex float
+    // Complex
-    inline u128f operator()(float a, float b){
+    template <typename T>
-      u128f out; 
+    inline vec<T> operator()(T a, T b){
-      out.f[0] = a;
+      vec<T> out;
-      out.f[1] = b;
+      
-      out.f[2] = a;
+      VECTOR_FOR(i, W<T>::r, 2)
-      out.f[3] = b;
+      {
+        out.v[i]   = a;
+        out.v[i+1] = b;
+      }
+
       return out;
     }
-    // Real float
+    
-    inline u128f operator()(float a){
+    // Real
-      u128f out; 
+    template <typename T>
-      out.f[0] = a;
+    inline vec<T> operator()(T a){
-      out.f[1] = a;
+      vec<T> out;
-      out.f[2] = a;
+      
-      out.f[3] = a;
+      VECTOR_FOR(i, W<T>::r, 1)
+      {
+        out.v[i] = a;
+      }
+      
       return out;
     }
-    //Complex double
+    
-    inline u128d operator()(double a, double b){
+    // Integer
-      u128d out; 
-      out.f[0] = a;
-      out.f[1] = b;
-      return out;
-    }
-    //Real double
-    inline u128d operator()(double a){
-      u128d out; 
-      out.f[0] = a;
-      out.f[1] = a;
-      return out;
-    }
-    //Integer
     inline int operator()(Integer a){
       return a;
     }
   };
 
   struct Vstore{
-    //Float 
+    // Real
-    inline void operator()(u128f a, float* F){
+    template <typename T>
-      memcpy(F,a.f,4*sizeof(float));
+    inline void operator()(vec<T> a, T *D){
-    }
+      *((vec<T> *)D) = a;
-    //Double
-    inline void operator()(u128d a, double* D){
-      memcpy(D,a.f,2*sizeof(double));
     }
     //Integer
-    inline void operator()(int a, Integer* I){
+    inline void operator()(int a, Integer *I){
-      I[0] = a;
+      *I = a;
     }
 
   };
 
   struct Vstream{
-    //Float
+    // Real
-    inline void operator()(float * a, u128f b){
+    template <typename T>
-      memcpy(a,b.f,4*sizeof(float));
+    inline void operator()(T * a, vec<T> b){
+      *((vec<T> *)a) = b;
     }
-    //Double
-    inline void operator()(double * a, u128d b){
-      memcpy(a,b.f,2*sizeof(double));
-    }
-
-
   };
 
   struct Vset{
-    // Complex float 
+    // Complex
-    inline u128f operator()(Grid::ComplexF *a){
+    template <typename T>
-      u128f out; 
+    inline vec<T> operator()(std::complex<T> *a){
-      out.f[0] = a[0].real();
+      vec<T> out;
-      out.f[1] = a[0].imag();
+      
-      out.f[2] = a[1].real();
+      VECTOR_FOR(i, W<T>::c, 1)
-      out.f[3] = a[1].imag();
+      {
+        out.v[2*i]   = a[i].real();
+        out.v[2*i+1] = a[i].imag();
+      }
+      
       return out;
     }
-    // Complex double 
+    
-    inline u128d operator()(Grid::ComplexD *a){
+    // Real
-      u128d out; 
+    template <typename T>
-      out.f[0] = a[0].real();
+    inline vec<T> operator()(T *a){
-      out.f[1] = a[0].imag();
+      vec<T> out;
-      return out;
+      
-    }
+      out = *((vec<T> *)a);
-    // Real float 
+      
-    inline u128f operator()(float *a){
-      u128f out; 
-      out.f[0] = a[0];
-      out.f[1] = a[1];
-      out.f[2] = a[2];
-      out.f[3] = a[3];
-      return out;
-    }
-    // Real double
-    inline u128d operator()(double *a){
-      u128d out; 
-      out.f[0] = a[0];
-      out.f[1] = a[1];
       return out;
     }
+
     // Integer
     inline int operator()(Integer *a){
-      return a[0];
+      return *a;
+    }
+  };
+
+  /////////////////////////////////////////////////////
+  // Arithmetic operations
+  /////////////////////////////////////////////////////
+  struct Sum{
+    // Complex/Real
+    template <typename T>
+    inline vec<T> operator()(vec<T> a, vec<T> b){
+      vec<T> out;
+      
+      VECTOR_FOR(i, W<T>::r, 1)
+      {
+        out.v[i] = a.v[i] + b.v[i];
+      }
+      
+      return out;
     }
     
-
+    //I nteger
+    inline int operator()(int a, int b){
+      return a + b;
+    }
   };
 
+  struct Sub{
+    // Complex/Real
+    template <typename T>
+    inline vec<T> operator()(vec<T> a, vec<T> b){
+      vec<T> out;
+      
+      VECTOR_FOR(i, W<T>::r, 1)
+      {
+        out.v[i] = a.v[i] - b.v[i];
+      }
+      
+      return out;
+    }
+    
+    //Integer
+    inline int operator()(int a, int b){
+      return a-b;
+    }
+  };
+
+  struct Mult{
+    // Real
+    template <typename T>
+    inline vec<T> operator()(vec<T> a, vec<T> b){
+      vec<T> out;
+      
+      VECTOR_FOR(i, W<T>::r, 1)
+      {
+        out.v[i] = a.v[i]*b.v[i];
+      }
+      
+      return out;
+    }
+    
+    // Integer
+    inline int operator()(int a, int b){
+      return a*b;
+    }
+  };
+  
+  #define cmul(a, b, c, i)\
+  c[i]   = a[i]*b[i]   - a[i+1]*b[i+1];\
+  c[i+1] = a[i]*b[i+1] + a[i+1]*b[i];
+  
+  struct MultComplex{
+    // Complex
+    template <typename T>
+    inline vec<T> operator()(vec<T> a, vec<T> b){
+      vec<T> out;
+      
+      VECTOR_FOR(i, W<T>::c, 1)
+      {
+        cmul(a.v, b.v, out.v, 2*i);
+      }      
+      
+      return out;
+    }
+  };
+  
+  #undef cmul
+
+  struct Div{
+    // Real
+    template <typename T>
+    inline vec<T> operator()(vec<T> a, vec<T> b){
+      vec<T> out;
+      
+      VECTOR_FOR(i, W<T>::r, 1)
+      {
+        out.v[i] = a.v[i]/b.v[i];
+      }
+      
+      return out;
+    }
+  };
+  
+  #define conj(a, b, i)\
+  b[i]   = a[i];\
+  b[i+1] = -a[i+1];
+  
+  struct Conj{
+    // Complex
+    template <typename T>
+    inline vec<T> operator()(vec<T> a){
+      vec<T> out;
+      
+      VECTOR_FOR(i, W<T>::c, 1)
+      {
+        conj(a.v, out.v, 2*i);
+      }
+      
+      return out;
+    }
+  };
+  
+  #undef conj
+
+  #define timesmi(a, b, i)\
+  b[i]   = a[i+1];\
+  b[i+1] = -a[i];
+  
+  struct TimesMinusI{
+    // Complex
+    template <typename T>
+    inline vec<T> operator()(vec<T> a, vec<T> b){
+      vec<T> out;
+      
+      VECTOR_FOR(i, W<T>::c, 1)
+      {
+        timesmi(a.v, out.v, 2*i);
+      }
+      
+      return out;
+    }
+  };
+
+  #undef timesmi
+  
+  #define timesi(a, b, i)\
+  b[i]   = -a[i+1];\
+  b[i+1] = a[i];
+  
+  struct TimesI{
+    // Complex
+    template <typename T>
+    inline vec<T> operator()(vec<T> a, vec<T> b){
+      vec<T> out;
+      
+      VECTOR_FOR(i, W<T>::c, 1)
+      {
+        timesi(a.v, out.v, 2*i);
+      }
+      
+      return out;
+    }
+  };
+  
+  #undef timesi
+
+  //////////////////////////////////////////////
+  // Some Template specialization
+  #define perm(a, b, n, w)\
+  unsigned int _mask = w >> (n + 1);\
+  VECTOR_FOR(i, w, 1)\
+  {\
+    b[i] = a[i^_mask];\
+  }
+  
+  #define DECL_PERMUTE_N(n)\
+  template <typename T>\
+  static inline vec<T> Permute##n(vec<T> in) {\
+    vec<T> out;\
+    perm(in.v, out.v, n, W<T>::r);\
+    return out;\
+  }
+  
+  struct Permute{
+    DECL_PERMUTE_N(0);
+    DECL_PERMUTE_N(1);
+    DECL_PERMUTE_N(2);
+    DECL_PERMUTE_N(3);
+  };
+  
+  #undef perm
+  #undef DECL_PERMUTE_N
+  
+  #define rot(a, b, n, w)\
+  VECTOR_FOR(i, w, 1)\
+  {\
+    b[i] = a[(i + n)%w];\
+  }
+  
+  struct Rotate{
+    template <typename T>
+    static inline vec<T> rotate(vec<T> in, int n){
+      vec<T> out;
+      
+      rot(in.v, out.v, n, W<T>::r);
+      
+      return out;
+    }
+  };
+
+  #undef rot
+  
+  #define acc(v, a, off, step, n)\
+  for (unsigned int i = off; i < n; i += step)\
+  {\
+    a += v[i];\
+  }
+  
   template <typename Out_type, typename In_type>
   struct Reduce{
     //Need templated class to overload output type
@@ -165,315 +383,66 @@ namespace Optimization {
     }
   };
   
-  /////////////////////////////////////////////////////
-  // Arithmetic operations
-  /////////////////////////////////////////////////////
-  struct Sum{
-    //Complex/Real float
-    inline u128f operator()(u128f a, u128f b){
-      u128f out;
-      out.f[0] = a.f[0] + b.f[0];
-      out.f[1] = a.f[1] + b.f[1];
-      out.f[2] = a.f[2] + b.f[2];
-      out.f[3] = a.f[3] + b.f[3];
-      return out;
-    }
-    //Complex/Real double
-    inline u128d operator()(u128d a, u128d b){
-      u128d out;
-      out.f[0] = a.f[0] + b.f[0];
-      out.f[1] = a.f[1] + b.f[1];
-      return out;
-    }
-    //Integer
-    inline int operator()(int a, int b){
-      return a + b;
-    }
-  };
-
-  struct Sub{
-    //Complex/Real float
-    inline u128f operator()(u128f a, u128f b){
-      u128f out;
-      out.f[0] = a.f[0] - b.f[0];
-      out.f[1] = a.f[1] - b.f[1];
-      out.f[2] = a.f[2] - b.f[2];
-      out.f[3] = a.f[3] - b.f[3];
-      return out;
-    }
-    //Complex/Real double
-    inline u128d operator()(u128d a, u128d b){
-      u128d out;
-      out.f[0] = a.f[0] - b.f[0];
-      out.f[1] = a.f[1] - b.f[1];
-      return out;
-    }
-    //Integer
-    inline int operator()(int a, int b){
-      return a-b;
-    }
-  };
-
-  struct MultComplex{
-    // Complex float
-    inline u128f operator()(u128f a, u128f b){
-      u128f out;
-      out.f[0] = a.f[0]*b.f[0] - a.f[1]*b.f[1];
-      out.f[1] = a.f[0]*b.f[1] + a.f[1]*b.f[0];
-      out.f[2] = a.f[2]*b.f[2] - a.f[3]*b.f[3];
-      out.f[3] = a.f[2]*b.f[3] + a.f[3]*b.f[2];
-      return out;
-    }
-    // Complex double
-    inline u128d operator()(u128d a, u128d b){
-      u128d out;
-      out.f[0] = a.f[0]*b.f[0] - a.f[1]*b.f[1];
-      out.f[1] = a.f[0]*b.f[1] + a.f[1]*b.f[0];
-      return out;
-    }
-  };
-
-  struct Mult{
-    //CK: Appear unneeded
-    // inline float  mac(float a, float b,double c){
-    //   return 0;
-    // }
-    // inline double mac(double a, double b,double c){
-    //   return 0;
-    // }
-
-    // Real float
-    inline u128f operator()(u128f a, u128f b){
-      u128f out;
-      out.f[0] = a.f[0]*b.f[0];
-      out.f[1] = a.f[1]*b.f[1];
-      out.f[2] = a.f[2]*b.f[2];
-      out.f[3] = a.f[3]*b.f[3];
-      return out;
-    }
-    // Real double
-    inline u128d operator()(u128d a, u128d b){
-      u128d out;
-      out.f[0] = a.f[0]*b.f[0];
-      out.f[1] = a.f[1]*b.f[1];
-      return out;
-    }
-    // Integer
-    inline int operator()(int a, int b){
-      return a*b;
-    }
-  };
-
-  struct Conj{
-    // Complex single
-    inline u128f operator()(u128f in){
-      u128f out;
-      out.f[0] = in.f[0];
-      out.f[1] = -in.f[1];
-      out.f[2] = in.f[2];
-      out.f[3] = -in.f[3];
-      return out;
-    }
-    // Complex double
-    inline u128d operator()(u128d in){
-      u128d out;
-      out.f[0] = in.f[0];
-      out.f[1] = -in.f[1];
-      return out;
-    }
-    // do not define for integer input
-  };
-
-  struct TimesMinusI{
-    //Complex single
-    inline u128f operator()(u128f in, u128f ret){ //note ret is ignored
-      u128f out;
-      out.f[0] = in.f[1];
-      out.f[1] = -in.f[0];
-      out.f[2] = in.f[3];
-      out.f[3] = -in.f[2];
-      return out;
-    }
-    //Complex double
-    inline u128d operator()(u128d in, u128d ret){
-      u128d out;
-      out.f[0] = in.f[1];
-      out.f[1] = -in.f[0];
-      return out;
-    }
-  };
-
-  struct TimesI{
-    //Complex single
-    inline u128f operator()(u128f in, u128f ret){ //note ret is ignored
-      u128f out;
-      out.f[0] = -in.f[1];
-      out.f[1] = in.f[0];
-      out.f[2] = -in.f[3];
-      out.f[3] = in.f[2];
-      return out;
-    }
-    //Complex double
-    inline u128d operator()(u128d in, u128d ret){
-      u128d out;
-      out.f[0] = -in.f[1];
-      out.f[1] = in.f[0];
-      return out;
-    }
-  };
-
-  //////////////////////////////////////////////
-  // Some Template specialization
-  struct Permute{
-    //We just have to mirror the permutes of Grid_sse4.h
-    static inline u128f Permute0(u128f in){ //AB CD -> CD AB
-      u128f out;
-      out.f[0] = in.f[2];
-      out.f[1] = in.f[3];
-      out.f[2] = in.f[0];
-      out.f[3] = in.f[1];
-      return out;
-    };
-    static inline u128f Permute1(u128f in){ //AB CD -> BA DC
-      u128f out;
-      out.f[0] = in.f[1];
-      out.f[1] = in.f[0];
-      out.f[2] = in.f[3];
-      out.f[3] = in.f[2];
-      return out;
-    };
-    static inline u128f Permute2(u128f in){
-      return in;
-    };
-    static inline u128f Permute3(u128f in){
-      return in;
-    };
-
-    static inline u128d Permute0(u128d in){ //AB -> BA
-      u128d out;
-      out.f[0] = in.f[1];
-      out.f[1] = in.f[0];
-      return out;      
-    };
-    static inline u128d Permute1(u128d in){
-      return in;
-    };
-    static inline u128d Permute2(u128d in){
-      return in;
-    };
-    static inline u128d Permute3(u128d in){
-      return in;
-    };
-
-  };
-  
-  template < typename vtype > 
-    void permute(vtype &a, vtype b, int perm) {
-   };
-    
-  struct Rotate{
-
-    static inline u128f rotate(u128f in,int n){
-      u128f out;
-      switch(n){
-      case 0:
-        out.f[0] = in.f[0];
-        out.f[1] = in.f[1];
-        out.f[2] = in.f[2];
-        out.f[3] = in.f[3];
-        break;
-      case 1:
-        out.f[0] = in.f[1];
-        out.f[1] = in.f[2];
-        out.f[2] = in.f[3];
-        out.f[3] = in.f[0];
-        break;
-      case 2:
-        out.f[0] = in.f[2];
-        out.f[1] = in.f[3];
-        out.f[2] = in.f[0];
-        out.f[3] = in.f[1];
-        break;
-      case 3:
-        out.f[0] = in.f[3];
-        out.f[1] = in.f[0];
-        out.f[2] = in.f[1];
-        out.f[3] = in.f[2];
-        break;
-      default: assert(0);
-      }
-      return out;
-    }
-    static inline u128d rotate(u128d in,int n){
-      u128d out;
-      switch(n){
-      case 0:
-        out.f[0] = in.f[0];
-        out.f[1] = in.f[1];
-        break;
-      case 1:
-        out.f[0] = in.f[1];
-        out.f[1] = in.f[0];
-        break;
-      default: assert(0);
-      }
-      return out;
-    }
-  };
-
   //Complex float Reduce
-  template<>
+  template <>
-  inline Grid::ComplexF Reduce<Grid::ComplexF, u128f>::operator()(u128f in){ //2 complex
+  inline Grid::ComplexF Reduce<Grid::ComplexF, vecf>::operator()(vecf in){
-    return Grid::ComplexF(in.f[0] + in.f[2], in.f[1] + in.f[3]);
+    float a = 0.f, b = 0.f;
+    
+    acc(in.v, a, 0, 2, W<float>::r);
+    acc(in.v, b, 1, 2, W<float>::r);
+    
+    return Grid::ComplexF(a, b);
   }
+  
   //Real float Reduce
   template<>
-  inline Grid::RealF Reduce<Grid::RealF, u128f>::operator()(u128f in){ //4 floats
+  inline Grid::RealF Reduce<Grid::RealF, vecf>::operator()(vecf in){
-    return in.f[0] + in.f[1] + in.f[2] + in.f[3];
+    float a = 0.;
-  }
     
+    acc(in.v, a, 0, 1, W<float>::r);
+    
+    return a;
+  }
   
   //Complex double Reduce
   template<>
-  inline Grid::ComplexD Reduce<Grid::ComplexD, u128d>::operator()(u128d in){ //1 complex
+  inline Grid::ComplexD Reduce<Grid::ComplexD, vecd>::operator()(vecd in){
-    return Grid::ComplexD(in.f[0],in.f[1]);
+    double a = 0., b = 0.;
+    
+    acc(in.v, a, 0, 2, W<double>::r);
+    acc(in.v, b, 1, 2, W<double>::r);
+    
+    return Grid::ComplexD(a, b);
   }
   
   //Real double Reduce
   template<>
-  inline Grid::RealD Reduce<Grid::RealD, u128d>::operator()(u128d in){ //2 doubles
+  inline Grid::RealD Reduce<Grid::RealD, vecd>::operator()(vecd in){
-    return in.f[0] + in.f[1];
+    double a = 0.f;
+    
+    acc(in.v, a, 0, 1, W<double>::r);
+    
+    return a;
   }
 
   //Integer Reduce
   template<>
   inline Integer Reduce<Integer, int>::operator()(int in){
-    // FIXME unimplemented
+    return in;
-   printf("Reduce : Missing integer implementation -> FIX\n");
-    assert(0);
   }
 }
 
 //////////////////////////////////////////////////////////////////////////////////////
 // Here assign types 
 
-  typedef Optimization::u128f SIMD_Ftype;  // Single precision type
+  typedef Optimization::vecf SIMD_Ftype; // Single precision type
-  typedef Optimization::u128d SIMD_Dtype; // Double precision type
+  typedef Optimization::vecd SIMD_Dtype; // Double precision type
   typedef int SIMD_Itype; // Integer type
 
   // prefetch utilities
   inline void v_prefetch0(int size, const char *ptr){};
   inline void prefetch_HINT_T0(const char *ptr){};
 
-
-
-  // Gpermute function
-  template < typename VectorSIMD > 
-    inline void Gpermute(VectorSIMD &y,const VectorSIMD &b, int perm ) {
-    Optimization::permute(y.v,b.v,perm);
-  }
-
-
   // Function name aliases
   typedef Optimization::Vsplat   VsplatSIMD;
   typedef Optimization::Vstore   VstoreSIMD;
@@ -481,16 +450,13 @@ namespace Optimization {
   typedef Optimization::Vstream  VstreamSIMD;
   template <typename S, typename T> using ReduceSIMD = Optimization::Reduce<S,T>;
 
- 
-
-
   // Arithmetic operations
   typedef Optimization::Sum         SumSIMD;
   typedef Optimization::Sub         SubSIMD;
+  typedef Optimization::Div         DivSIMD;
   typedef Optimization::Mult        MultSIMD;
   typedef Optimization::MultComplex MultComplexSIMD;
   typedef Optimization::Conj        ConjSIMD;
   typedef Optimization::TimesMinusI TimesMinusISIMD;
   typedef Optimization::TimesI      TimesISIMD;
-
 }

lib/simd/Grid_imci.h

@@ -26,14 +26,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-//----------------------------------------------------------------------
-/*! @file Grid_knc.h
-  @brief Optimization libraries for AVX512 instructions set for KNC
-
-  Using intrinsics
-*/
-// Time-stamp: <2015-06-09 14:27:28 neo>
-//----------------------------------------------------------------------
 
 #include <immintrin.h>
 #include <zmmintrin.h>

lib/simd/Grid_qpx.h

@@ -245,6 +245,21 @@ namespace Optimization {
     }
   };
 
+  struct Div{
+    // Real double
+    inline vector4double operator()(vector4double a, vector4double b){
+      return vec_swdiv(a, b);
+    }
+
+    // Real float
+    FLOAT_WRAP_2(operator(), inline)
+
+    // Integer
+    inline int operator()(int a, int b){
+      return a/b;
+    }
+  };
+
   struct Conj{
     // Complex double
     inline vector4double operator()(vector4double v){
@@ -413,6 +428,7 @@ template <typename S, typename T> using ReduceSIMD = Optimization::Reduce<S,T>;
 typedef Optimization::Sum         SumSIMD;
 typedef Optimization::Sub         SubSIMD;
 typedef Optimization::Mult        MultSIMD;
+typedef Optimization::Div         DivSIMD;
 typedef Optimization::MultComplex MultComplexSIMD;
 typedef Optimization::Conj        ConjSIMD;
 typedef Optimization::TimesMinusI TimesMinusISIMD;

lib/simd/Grid_vector_types.h

@@ -38,13 +38,13 @@ directory
 #ifndef GRID_VECTOR_TYPES
 #define GRID_VECTOR_TYPES
 
-#ifdef GENERIC_VEC
+#ifdef GEN
 #include "Grid_generic.h"
 #endif
 #ifdef SSE4
 #include "Grid_sse4.h"
 #endif
-#if defined(AVX1) || defined(AVX2) || defined(AVXFMA4)
+#if defined(AVX1) || defined (AVXFMA) || defined(AVX2) || defined(AVXFMA4)
 #include "Grid_avx.h"
 #endif
 #if defined AVX512
@@ -130,7 +130,7 @@ class Grid_simd {
 
   Vector_type v;
 
-  static inline int Nsimd(void) {
+  static inline constexpr int Nsimd(void) {
     return sizeof(Vector_type) / sizeof(Scalar_type);
   }
 

lib/stencil/Lebesgue.cc

@@ -65,7 +65,7 @@ void LebesgueOrder::CartesianBlocking(void)
 {
   _LebesgueReorder.resize(0);
 
-  std::cout << GridLogDebug << " CartesianBlocking ";
+  //    std::cout << GridLogDebug << " CartesianBlocking ";
   //    for(int d=0;d<Block.size();d++) std::cout <<Block[d]<<" ";
   //    std::cout<<std::endl; 
 

m4/ac_prog_doxygen.m4

@@ -1,54 +0,0 @@
-dnl Check for doxygen to create API docs
-dnl
-AC_DEFUN([AC_PROG_DOXYGEN],
-[
-AC_ARG_ENABLE(doxygen,
-AS_HELP_STRING([--enable-doxygen],[enable documentation generation with doxygen (auto)]))
-
-AC_ARG_ENABLE(dot,
-AS_HELP_STRING([--enable-dot],[use 'dot' to generate graphs in doxygen (auto)]))              
-AC_ARG_ENABLE(html-docs,
-AS_HELP_STRING([--enable-html-docs],[enable HTML generation with doxygen (yes)]),
-[],[ enable_html_docs=yes])              
-AC_ARG_ENABLE(latex-docs,
-AS_HELP_STRING([--enable-latex-docs],
-[enable LaTeX documentation generation with doxygen (no)]),[],[enable_latex_docs=no])              
-
-if test "x$enable_doxygen" = xno; then
- enable_doc=no	
-else 
- AC_CHECK_PROG(DOXYGEN, doxygen, doxygen)
- if test x$DOXYGEN = x; then
-  if test "x$enable_doxygen" = xyes; then
-   AC_MSG_ERROR([could not find doxygen])
-  fi
-  enable_doc=no
- else
-  doxy_ver=`doxygen --version`
-  doxy_major=`expr "$doxy_ver" : '\(@<:@0-9@:>@\)\..*'`
-  doxy_minor=`expr "$doxy_ver" : '@<:@0-9@:>@\.\(@<:@0-9@:>@\).*'`
-  if test $doxy_major -eq "1" -a $doxy_minor -ge "3" ; then
-   enable_doc=yes
-   AC_CHECK_PROG(DOT, dot, dot)
-  else
-   AC_MSG_WARN([doxygen version $doxy_ver too old, doxygen will not be used.])
-   enable_doc=no
-  fi
- fi
-fi
-
-AM_CONDITIONAL(DOXYGEN_DOC, test x$enable_doc = xyes)
-
-if test x$DOT = x; then
- if test "x$enable_dot" = xyes; then
-  AC_MSG_ERROR([could not find dot])
- fi
- enable_dot=no
-else
- enable_dot=yes
-fi
-
-AC_SUBST(enable_dot)
-AC_SUBST(enable_html_docs)
-AC_SUBST(enable_latex_docs)
-])

programs/Hadrons/Application.cc

@@ -1,195 +0,0 @@
-/*******************************************************************************
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: programs/Hadrons/Application.cc
-
-Copyright (C) 2015
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-
-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.
-*******************************************************************************/
-
-#include <Hadrons/Application.hpp>
-#include <Hadrons/GeneticScheduler.hpp>
-
-using namespace Grid;
-using namespace QCD;
-using namespace Hadrons;
-
-#define BIG_SEP "==============="
-#define SEP     "---------------"
-
-/******************************************************************************
- *                       Application implementation                           *
- ******************************************************************************/
-// constructor /////////////////////////////////////////////////////////////////
-Application::Application(const std::string parameterFileName)
-: parameterFileName_(parameterFileName)
-, env_(Environment::getInstance())
-{
-    LOG(Message) << "Modules available:" << std::endl;
-    auto list = ModuleFactory::getInstance().getBuilderList();
-    for (auto &m: list)
-    {
-        LOG(Message) << "  " << m << std::endl;
-    }
-    auto dim = GridDefaultLatt(), mpi = GridDefaultMpi(), loc(dim);
-    locVol_ = 1;
-    for (unsigned int d = 0; d < dim.size(); ++d)
-    {
-        loc[d]  /= mpi[d];
-        locVol_ *= loc[d];
-    }
-    LOG(Message) << "Global lattice: " << dim << std::endl;
-    LOG(Message) << "MPI partition : " << mpi << std::endl;
-    LOG(Message) << "Local lattice : " << loc << std::endl;
-}
-
-// destructor //////////////////////////////////////////////////////////////////
-Application::~Application(void)
-{}
-
-// execute /////////////////////////////////////////////////////////////////////
-void Application::run(void)
-{
-    parseParameterFile();
-    schedule();
-    configLoop();
-}
-
-// parse parameter file ////////////////////////////////////////////////////////
-class ObjectId: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(ObjectId,
-                                    std::string, name,
-                                    std::string, type);
-};
-
-void Application::parseParameterFile(void)
-{
-    XmlReader reader(parameterFileName_);
-    ObjectId  id;
-    
-    LOG(Message) << "Reading '" << parameterFileName_ << "'..." << std::endl;
-    read(reader, "parameters", par_);
-    push(reader, "modules");
-    push(reader, "module");
-    do
-    {
-        read(reader, "id", id);
-        env_.createModule(id.name, id.type, reader);
-    } while (reader.nextElement("module"));
-    pop(reader);
-    pop(reader);
-    env_.setSeed(strToVec<int>(par_.seed));
-    env_.printContent();
-}
-
-// schedule computation ////////////////////////////////////////////////////////
-#define MEM_MSG(size)\
-sizeString((size)*locVol_) << " (" << sizeString(size)  << "/site)"
-
-void Application::schedule(void)
-{
-    // memory peak function
-    auto memPeak = [this](const std::vector<unsigned int> &program)
-    {
-        unsigned int memPeak;
-        bool         msg;
-        
-        msg = HadronsLogMessage.isActive();
-        HadronsLogMessage.Active(false);
-        env_.dryRun(true);
-        memPeak = env_.executeProgram(program);
-        env_.dryRun(false);
-        env_.freeAll();
-        HadronsLogMessage.Active(true);
-        
-        return memPeak;
-    };
-    
-    // constrained topological sort using a genetic algorithm
-    LOG(Message) << "Scheduling computation..." << std::endl;
-    constexpr unsigned int maxGen = 200, maxCstGen = 50;
-    unsigned int           k = 0, gen, prevPeak, nCstPeak = 0;
-    auto                   graph = env_.makeModuleGraph();
-    auto                   con = graph.getConnectedComponents();
-    std::random_device     rd;
-    GeneticScheduler<unsigned int>::Parameters par;
-
-    par.popSize      = 20;
-    par.mutationRate = .1;
-    par.seed         = rd();
-    CartesianCommunicator::BroadcastWorld(0, &(par.seed), sizeof(par.seed));
-    for (unsigned int i = 0; i < con.size(); ++i)
-    {
-        GeneticScheduler<unsigned int> scheduler(con[i], memPeak, par);
-        
-        gen = 0;
-        do
-        {
-            scheduler.nextGeneration();
-            if (gen != 0)
-            {
-                if (prevPeak == scheduler.getMinValue())
-                {
-                    nCstPeak++;
-                }
-                else
-                {
-                    nCstPeak = 0;
-                }
-            }
-            prevPeak = scheduler.getMinValue();
-            if (gen % 10 == 0)
-            {
-                LOG(Iterative) << "Generation " << gen << ": "
-                               << MEM_MSG(scheduler.getMinValue()) << std::endl;
-            }
-            gen++;
-        } while ((gen < maxGen) and (nCstPeak < maxCstGen));
-        auto &t = scheduler.getMinSchedule();
-        LOG(Message) << "Program " << i + 1 << " (memory peak: "
-                     << MEM_MSG(scheduler.getMinValue()) << "):" << std::endl;
-        for (unsigned int j = 0; j < t.size(); ++j)
-        {
-            program_.push_back(t[j]);
-            LOG(Message) << std::setw(4) << std::right << k + 1 << ": "
-                         << env_.getModuleName(program_[k]) << std::endl;
-            k++;
-        }
-    }
-}
-
-// loop on configurations //////////////////////////////////////////////////////
-void Application::configLoop(void)
-{
-    auto range = par_.trajCounter;
-    
-    for (unsigned int t = range.start; t < range.end; t += range.step)
-    {
-        LOG(Message) << BIG_SEP << " Starting measurement for trajectory " << t
-                     << " " << BIG_SEP << std::endl;
-        env_.setTrajectory(t);
-        env_.executeProgram(program_);
-        env_.freeAll();
-    }
-    LOG(Message) << BIG_SEP << " End of measurement " << BIG_SEP << std::endl;
-}

programs/Hadrons/Application.hpp

@@ -1,85 +0,0 @@
-/*******************************************************************************
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: programs/Hadrons/Application.hpp
-
-Copyright (C) 2015
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-
-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.
-*******************************************************************************/
-
-#ifndef Hadrons_Application_hpp_
-#define Hadrons_Application_hpp_
-
-#include <Hadrons/Global.hpp>
-#include <Hadrons/Environment.hpp>
-#include <Hadrons/ModuleFactory.hpp>
-
-BEGIN_HADRONS_NAMESPACE
-
-class TrajRange: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(TrajRange,
-                                    unsigned int, start,
-                                    unsigned int, end,
-                                    unsigned int, step);
-};
-
-class GlobalPar: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(GlobalPar,
-                                    TrajRange,   trajCounter,
-                                    std::string, seed);
-};
-
-/******************************************************************************
- *                         Main program manager                               *
- ******************************************************************************/
-class Application
-{
-public:
-
-public:
-    // constructor
-    Application(const std::string parameterFileName);
-    // destructor
-    virtual ~Application(void);
-    // execute
-    void run(void);
-private:
-    // parse parameter file
-    void parseParameterFile(void);
-    // schedule computation
-    void schedule(void);
-    // loop on configurations
-    void configLoop(void);
-private:
-    long unsigned int                               locVol_;
-    std::string                                     parameterFileName_;
-    GlobalPar                                       par_;
-    Environment                                     &env_;
-    std::vector<unsigned int>                       program_;
-};
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_Application_hpp_

programs/Hadrons/GeneticScheduler.hpp

@@ -1,236 +0,0 @@
-/*******************************************************************************
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: programs/Hadrons/GeneticScheduler.hpp
-
-Copyright (C) 2016
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-
-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.
-*******************************************************************************/
-
-#ifndef Hadrons_GeneticScheduler_hpp_
-#define Hadrons_GeneticScheduler_hpp_
-
-#include <Hadrons/Global.hpp>
-#include <Hadrons/Graph.hpp>
-
-BEGIN_HADRONS_NAMESPACE
-
-/******************************************************************************
- *                   Scheduler based on a genetic algorithm                   *
- ******************************************************************************/
-template <typename T>
-class GeneticScheduler
-{
-public:
-    typedef std::function<int(const std::vector<T> &)> ObjFunc;
-    struct Parameters
-    {
-        double       mutationRate;
-        unsigned int popSize, seed;
-    };
-public:
-    // constructor
-    GeneticScheduler(Graph<T> &graph, const ObjFunc &func,
-                     const Parameters &par);
-    // destructor
-    virtual ~GeneticScheduler(void) = default;
-    // access
-    const std::vector<T> & getMinSchedule(void);
-    int                    getMinValue(void);
-    // breed a new generation
-    void nextGeneration(void);
-    // print population
-    friend std::ostream & operator<<(std::ostream &out,
-                                     const GeneticScheduler<T> &s)
-    {
-        for (auto &p: s.population_)
-        {
-            out << p.second << ": " << p.first << std::endl;
-        }
-        
-        return out;
-    }
-private:
-    // randomly initialize population
-    void initPopulation(void);
-    // genetic operators
-    const std::vector<T> & selection(void);
-    void                   crossover(const std::vector<T> &c1,
-                                     const std::vector<T> &c2);
-    void                   mutation(std::vector<T> &c);
-private:
-    Graph<T>                           &graph_;
-    const ObjFunc                      &func_;
-    const Parameters                   par_;
-    std::multimap<int, std::vector<T>> population_;
-    std::mt19937                       gen_;
-};
-
-/******************************************************************************
- *                       template implementation                              *
- ******************************************************************************/
-// constructor /////////////////////////////////////////////////////////////////
-template <typename T>
-GeneticScheduler<T>::GeneticScheduler(Graph<T> &graph, const ObjFunc &func,
-                                      const Parameters &par)
-: graph_(graph)
-, func_(func)
-, par_(par)
-{
-    gen_.seed(par_.seed);
-}
-
-// access //////////////////////////////////////////////////////////////////////
-template <typename T>
-const std::vector<T> & GeneticScheduler<T>::getMinSchedule(void)
-{
-    return population_.begin()->second;
-}
-
-template <typename T>
-int GeneticScheduler<T>::getMinValue(void)
-{
-    return population_.begin()->first;
-}
-
-// breed a new generation //////////////////////////////////////////////////////
-template <typename T>
-void GeneticScheduler<T>::nextGeneration(void)
-{
-    std::uniform_real_distribution<double> dis(0., 1.);
-    
-    // random initialization of the population if necessary
-    if (population_.size() != par_.popSize)
-    {
-        initPopulation();
-    }
-
-    // mating
-    for (unsigned int i = 0; i < par_.popSize/2; ++i)
-    {
-        auto &p1 = selection(), &p2 = selection();
-        crossover(p1, p2);
-    }
-    
-    // random mutations
-    auto buf = population_;
-    population_.clear();
-    for (auto &c: buf)
-    {
-        if (dis(gen_) < par_.mutationRate)
-        {
-            mutation(c.second);
-        }
-        population_.emplace(func_(c.second), c.second);
-    }
-    
-    // grim reaper
-    auto it = population_.begin();
-    
-    std::advance(it, par_.popSize);
-    population_.erase(it, population_.end());
-}
-
-// randomly initialize population //////////////////////////////////////////////
-template <typename T>
-void GeneticScheduler<T>::initPopulation(void)
-{
-    population_.clear();
-    for (unsigned int i = 0; i < par_.popSize; ++i)
-    {
-        auto p = graph_.topoSort(gen_);
-        
-        population_.emplace(func_(p), p);
-    }
-}
-
-// genetic operators ///////////////////////////////////////////////////////////
-template <typename T>
-const std::vector<T> & GeneticScheduler<T>::selection(void)
-{
-    std::vector<double> prob;
-    
-    for (auto &c: population_)
-    {
-        prob.push_back(1./c.first);
-    }
-    std::discrete_distribution<unsigned int> dis(prob.begin(), prob.end());
-    auto rIt = population_.begin();
-    std::advance(rIt, dis(gen_));
-    
-    return rIt->second;
-}
-
-template <typename T>
-void GeneticScheduler<T>::crossover(const std::vector<T> &p1,
-                                    const std::vector<T> &p2)
-{
-    std::uniform_int_distribution<unsigned int> dis(0, p1.size() - 1);
-    unsigned int                                cut = dis(gen_);
-    std::vector<T>                              c1, c2, buf;
-    
-    auto cross = [&buf, cut](std::vector<T> &c, const std::vector<T> &p1,
-                             const std::vector<T> &p2)
-    {
-        buf = p2;
-        for (unsigned int i = 0; i < cut; ++i)
-        {
-            c.push_back(p1[i]);
-            buf.erase(std::find(buf.begin(), buf.end(), p1[i]));
-        }
-        for (unsigned int i = 0; i < buf.size(); ++i)
-        {
-            c.push_back(buf[i]);
-        }
-    };
-    
-    cross(c1, p1, p2);
-    cross(c2, p2, p1);
-    population_.emplace(func_(c1), c1);
-    population_.emplace(func_(c2), c2);
-}
-
-template <typename T>
-void GeneticScheduler<T>::mutation(std::vector<T> &c)
-{
-    std::uniform_int_distribution<unsigned int> dis(0, c.size() - 1);
-    unsigned int                                cut = dis(gen_);
-    Graph<T>                                    g = graph_;
-    std::vector<T>                              buf;
-    
-    for (unsigned int i = cut; i < c.size(); ++i)
-    {
-        g.removeVertex(c[i]);
-    }
-    if (g.size() > 0)
-    {
-        buf = g.topoSort(gen_);
-    }
-    for (unsigned int i = cut; i < c.size(); ++i)
-    {
-        buf.push_back(c[i]);
-    }
-    c = buf;
-}
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_GeneticScheduler_hpp_

programs/Hadrons/Make.inc

@@ -1,4 +0,0 @@
-bin_PROGRAMS = Hadrons
-
-Hadrons_SOURCES = Hadrons.cc
-Hadrons_LDADD   = -lGrid

programs/Hadrons/Makefile.am

@@ -1,15 +0,0 @@
-AM_CXXFLAGS += -I$(top_srcdir)/programs -I../$(top_srcdir)/programs
-
-bin_PROGRAMS = Hadrons
-
-include modules.inc
-
-Hadrons_SOURCES =   \
-    $(modules)      \
-    Application.cc  \
-    Environment.cc  \
-    Global.cc       \
-    Hadrons.cc      \
-    Module.cc
-
-Hadrons_LDADD   = -lGrid

programs/Hadrons/Modules/ADWF.cc

@@ -1,81 +0,0 @@
-/*******************************************************************************
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: programs/Hadrons/ADWF.cc
-
-Copyright (C) 2016
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-
-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.
-*******************************************************************************/
-
-#include <Hadrons/Modules/ADWF.hpp>
-
-using namespace Grid;
-using namespace Hadrons;
-
-/******************************************************************************
-*                          ADWF implementation                                *
-******************************************************************************/
-// constructor /////////////////////////////////////////////////////////////////
-ADWF::ADWF(const std::string name)
-: Module<ADWFPar>(name)
-{}
-
-// dependencies/products ///////////////////////////////////////////////////////
-std::vector<std::string> ADWF::getInput(void)
-{
-    std::vector<std::string> in = {par().gauge};
-    
-    return in;
-}
-
-std::vector<std::string> ADWF::getOutput(void)
-{
-    std::vector<std::string> out = {getName()};
-    
-    return out;
-}
-
-// setup ///////////////////////////////////////////////////////////////////////
-void ADWF::setup(void)
-{
-    unsigned int size;
-    
-    size = 3*env().lattice4dSize<WilsonFermionR::DoubledGaugeField>();
-    env().registerObject(getName(), size, par().Ls);
-}
-
-// execution ///////////////////////////////////////////////////////////////////
-void ADWF::execute(void)
-{
-    LOG(Message) << "Setting up domain wall fermion matrix with m= "
-                 << par().mass << ", M5= " << par().M5 << " and Ls= "
-                 << par().Ls << " using gauge field '" << par().gauge << "'"
-                 << std::endl;
-    env().createGrid(par().Ls);
-    auto &U      = *env().getObject<LatticeGaugeField>(par().gauge);
-    auto &g4     = *env().getGrid();
-    auto &grb4   = *env().getRbGrid();
-    auto &g5     = *env().getGrid(par().Ls);
-    auto &grb5   = *env().getRbGrid(par().Ls);
-    FMat *fMatPt = new DomainWallFermionR(U, g5, grb5, g4, grb4, par().mass,
-                                          par().M5);
-    env().setObject(getName(), fMatPt);
-}

programs/Hadrons/Modules/ADWF.hpp

@@ -1,70 +0,0 @@
-/*******************************************************************************
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: programs/Hadrons/ADWF.hpp
-
-Copyright (C) 2016
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-
-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.
-*******************************************************************************/
-
-#ifndef Hadrons_ADWF_hpp_
-#define Hadrons_ADWF_hpp_
-
-#include <Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
-
-BEGIN_HADRONS_NAMESPACE
-
-/******************************************************************************
- *                     Domain wall quark action                               *
- ******************************************************************************/
-class ADWFPar: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(ADWFPar, 
-                                    std::string, gauge,
-                                    unsigned int, Ls,
-                                    double      , mass,
-                                    double      , M5);
-};
-
-class ADWF: public Module<ADWFPar>
-{
-public:
-    // constructor
-    ADWF(const std::string name);
-    // destructor
-    virtual ~ADWF(void) = default;
-    // dependency relation
-    virtual std::vector<std::string> getInput(void);
-    virtual std::vector<std::string> getOutput(void);
-    // setup
-    virtual void setup(void);
-    // execution
-    virtual void execute(void);
-};
-
-MODULE_REGISTER(ADWF);
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_ADWF_hpp_

programs/Hadrons/Modules/AWilson.cc

@@ -1,75 +0,0 @@
-/*******************************************************************************
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: programs/Hadrons/AWilson.cc
-
-Copyright (C) 2016
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-
-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.
-*******************************************************************************/
-
-#include <Hadrons/Modules/AWilson.hpp>
-
-using namespace Grid;
-using namespace Hadrons;
-
-/******************************************************************************
-*                         AWilson implementation                              *
-******************************************************************************/
-// constructor /////////////////////////////////////////////////////////////////
-AWilson::AWilson(const std::string name)
-: Module<AWilsonPar>(name)
-{}
-
-// dependencies/products ///////////////////////////////////////////////////////
-std::vector<std::string> AWilson::getInput(void)
-{
-    std::vector<std::string> in = {par().gauge};
-    
-    return in;
-}
-
-std::vector<std::string> AWilson::getOutput(void)
-{
-    std::vector<std::string> out = {getName()};
-    
-    return out;
-}
-
-// setup ///////////////////////////////////////////////////////////////////////
-void AWilson::setup(void)
-{
-    unsigned int size;
-    
-    size = 3*env().lattice4dSize<WilsonFermionR::DoubledGaugeField>();
-    env().registerObject(getName(), size);
-}
-
-// execution ///////////////////////////////////////////////////////////////////
-void AWilson::execute()
-{
-    LOG(Message) << "Setting up Wilson fermion matrix with m= " << par().mass
-                 << " using gauge field '" << par().gauge << "'" << std::endl;
-    auto &U      = *env().getObject<LatticeGaugeField>(par().gauge);
-    auto &grid   = *env().getGrid();
-    auto &gridRb = *env().getRbGrid();
-    FMat *fMatPt = new WilsonFermionR(U, grid, gridRb, par().mass);
-    env().setObject(getName(), fMatPt);
-}

programs/Hadrons/Modules/AWilson.hpp

@@ -1,68 +0,0 @@
-/*******************************************************************************
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: programs/Hadrons/AWilson.hpp
-
-Copyright (C) 2016
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-
-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.
-*******************************************************************************/
-
-#ifndef Hadrons_AWilson_hpp_
-#define Hadrons_AWilson_hpp_
-
-#include <Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
-
-BEGIN_HADRONS_NAMESPACE
-
-/******************************************************************************
- *                            Wilson quark action                             *
- ******************************************************************************/
-class AWilsonPar: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(AWilsonPar,
-                                    std::string, gauge,
-                                    double     , mass);
-};
-
-class AWilson: public Module<AWilsonPar>
-{
-public:
-    // constructor
-    AWilson(const std::string name);
-    // destructor
-    virtual ~AWilson(void) = default;
-    // dependencies/products
-    virtual std::vector<std::string> getInput(void);
-    virtual std::vector<std::string> getOutput(void);
-    // setup
-    virtual void setup(void);
-    // execution
-    virtual void execute(void);
-};
-
-MODULE_REGISTER(AWilson);
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_AWilson_hpp_

programs/Hadrons/Modules/CMeson.hpp

@@ -1,74 +0,0 @@
-/*******************************************************************************
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: programs/Hadrons/CMeson.hpp
-
-Copyright (C) 2015
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-
-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.
-*******************************************************************************/
-
-#ifndef Hadrons_CMeson_hpp_
-#define Hadrons_CMeson_hpp_
-
-#include <Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
-
-BEGIN_HADRONS_NAMESPACE
-
-/******************************************************************************
- *                               CMeson                                       *
- ******************************************************************************/
-class CMesonPar: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(CMesonPar,
-                                    std::string, q1,
-                                    std::string, q2,
-                                    std::string, output);
-};
-
-class CMeson: public Module<CMesonPar>
-{
-public:
-    class Result: Serializable
-    {
-    public:
-        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
-                        std::vector<std::vector<std::vector<Complex>>>, corr);
-    };
-public:
-    // constructor
-    CMeson(const std::string name);
-    // destructor
-    virtual ~CMeson(void) = default;
-    // dependencies/products
-    virtual std::vector<std::string> getInput(void);
-    virtual std::vector<std::string> getOutput(void);
-    // execution
-    virtual void execute(void);
-};
-
-MODULE_REGISTER(CMeson);
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_CMeson_hpp_

programs/Hadrons/Modules/MQuark.hpp

@@ -1,71 +0,0 @@
-/*******************************************************************************
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: programs/Hadrons/MQuark.hpp
-
-Copyright (C) 2015
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-
-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.
-*******************************************************************************/
-
-#ifndef Hadrons_MQuark_hpp_
-#define Hadrons_MQuark_hpp_
-
-#include <Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
-
-BEGIN_HADRONS_NAMESPACE
-
-/******************************************************************************
- *                               MQuark                                       *
- ******************************************************************************/
-class MQuarkPar: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(MQuarkPar,
-                                    std::string, source,
-                                    std::string, solver);
-};
-
-class MQuark: public Module<MQuarkPar>
-{
-public:
-    // constructor
-    MQuark(const std::string name);
-    // destructor
-    virtual ~MQuark(void) = default;
-    // dependencies/products
-    virtual std::vector<std::string> getInput(void);
-    virtual std::vector<std::string> getOutput(void);
-    // setup
-    virtual void setup(void);
-    // execution
-    virtual void execute(void);
-private:
-    unsigned int        Ls_;
-    Environment::Solver *solver_{nullptr};
-};
-
-MODULE_REGISTER(MQuark);
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_MQuark_hpp_

programs/Hadrons/Modules/SolRBPrecCG.cc

@@ -1,83 +0,0 @@
-/*******************************************************************************
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: programs/Hadrons/SolRBPrecCG.cc
-
-Copyright (C) 2016
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-
-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.
-*******************************************************************************/
-
-#include <Hadrons/Modules/SolRBPrecCG.hpp>
-
-using namespace Grid;
-using namespace QCD;
-using namespace Hadrons;
-
-/******************************************************************************
-*                       SolRBPrecCG implementation                            *
-******************************************************************************/
-// constructor /////////////////////////////////////////////////////////////////
-SolRBPrecCG::SolRBPrecCG(const std::string name)
-: Module(name)
-{}
-
-// dependencies/products ///////////////////////////////////////////////////////
-std::vector<std::string> SolRBPrecCG::getInput(void)
-{
-    std::vector<std::string> in = {par().action};
-    
-    return in;
-}
-
-std::vector<std::string> SolRBPrecCG::getOutput(void)
-{
-    std::vector<std::string> out = {getName()};
-    
-    return out;
-}
-
-// setup ///////////////////////////////////////////////////////////////////////
-void SolRBPrecCG::setup(void)
-{
-    auto Ls = env().getObjectLs(par().action);
-    
-    env().registerObject(getName(), 0, Ls);
-    env().addOwnership(getName(), par().action);
-}
-
-// execution ///////////////////////////////////////////////////////////////////
-void SolRBPrecCG::execute(void)
-{
-    auto &mat   = *(env().getObject<Environment::FMat>(par().action));
-    auto solver = [&mat, this](LatticeFermion &sol,
-                               const LatticeFermion &source)
-    {
-        ConjugateGradient<LatticeFermion>           cg(par().residual, 10000);
-        SchurRedBlackDiagMooeeSolve<LatticeFermion> schurSolver(cg);
-        
-        schurSolver(mat, source, sol);
-    };
-    
-    LOG(Message) << "setting up Schur red-black preconditioned CG for"
-                 << " action '" << par().action << "' with residual "
-                 << par().residual << std::endl;
-    env().setObject(getName(), new Environment::Solver(solver));
-}

programs/Hadrons/Modules/SolRBPrecCG.hpp

@@ -1,68 +0,0 @@
-/*******************************************************************************
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: programs/Hadrons/SolRBPrecCG.hpp
-
-Copyright (C) 2016
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-
-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.
-*******************************************************************************/
-
-#ifndef Hadrons_SolRBPrecCG_hpp_
-#define Hadrons_SolRBPrecCG_hpp_
-
-#include <Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
-
-BEGIN_HADRONS_NAMESPACE
-
-/******************************************************************************
- *                     Schur red-black preconditioned CG                      *
- ******************************************************************************/
-class SolRBPrecCGPar: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(SolRBPrecCGPar,
-                                    std::string, action,
-                                    double     , residual);
-};
-
-class SolRBPrecCG: public Module<SolRBPrecCGPar>
-{
-public:
-    // constructor
-    SolRBPrecCG(const std::string name);
-    // destructor
-    virtual ~SolRBPrecCG(void) = default;
-    // dependencies/products
-    virtual std::vector<std::string> getInput(void);
-    virtual std::vector<std::string> getOutput(void);
-    // setup
-    virtual void setup(void);
-    // execution
-    virtual void execute(void);
-};
-
-MODULE_REGISTER(SolRBPrecCG);
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_SolRBPrecCG_hpp_

programs/Hadrons/Modules/SrcPoint.cc

@@ -1,74 +0,0 @@
-/*******************************************************************************
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: programs/Hadrons/SrcPoint.cc
-
-Copyright (C) 2016
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-
-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.
-*******************************************************************************/
-
-#include <Hadrons/Modules/SrcPoint.hpp>
-
-using namespace Grid;
-using namespace Hadrons;
-
-/******************************************************************************
-*                         SrcPoint implementation                             *
-******************************************************************************/
-// constructor /////////////////////////////////////////////////////////////////
-SrcPoint::SrcPoint(const std::string name)
-: Module<SrcPointPar>(name)
-{}
-
-// dependencies/products ///////////////////////////////////////////////////////
-std::vector<std::string> SrcPoint::getInput(void)
-{
-    std::vector<std::string> in;
-    
-    return in;
-}
-
-std::vector<std::string> SrcPoint::getOutput(void)
-{
-    std::vector<std::string> out = {getName()};
-    
-    return out;
-}
-
-// setup ///////////////////////////////////////////////////////////////////////
-void SrcPoint::setup(void)
-{
-    env().registerLattice<LatticePropagator>(getName());
-}
-
-// execution ///////////////////////////////////////////////////////////////////
-void SrcPoint::execute(void)
-{
-    std::vector<int> position = strToVec<int>(par().position);
-    SpinColourMatrix id;
-    
-    LOG(Message) << "Creating point source at position [" << par().position
-                 << "]" << std::endl;
-    LatticePropagator &src = *env().createLattice<LatticePropagator>(getName());
-    id  = 1.;
-    src = zero;
-    pokeSite(id, src, position);
-}

programs/Hadrons/Modules/SrcPoint.hpp

@@ -1,78 +0,0 @@
-/*******************************************************************************
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: programs/Hadrons/SrcPoint.hpp
-
-Copyright (C) 2016
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-
-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.
-*******************************************************************************/
-
-#ifndef Hadrons_SrcPoint_hpp_
-#define Hadrons_SrcPoint_hpp_
-
-#include <Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
-
-BEGIN_HADRONS_NAMESPACE
-
-/*
- 
- Point source
- ------------
- * src_x = delta_x,o
- 
- * options: o
- - position: space-separated integer sequence (e.g. "0 1 1 0")
- 
- */
-
-/******************************************************************************
- *                                SrcPoint                                    *
- ******************************************************************************/
-class SrcPointPar: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(SrcPointPar,
-                                    std::string, position);
-};
-
-class SrcPoint: public Module<SrcPointPar>
-{
-public:
-    // constructor
-    SrcPoint(const std::string name);
-    // destructor
-    virtual ~SrcPoint(void) = default;
-    // dependency relation
-    virtual std::vector<std::string> getInput(void);
-    virtual std::vector<std::string> getOutput(void);
-    // setup
-    virtual void setup(void);
-    // execution
-    virtual void execute(void);
-};
-
-MODULE_REGISTER(SrcPoint);
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_SrcPoint_hpp_

programs/Hadrons/Modules/SrcZ2.cc

@@ -1,87 +0,0 @@
-/*******************************************************************************
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: programs/Hadrons/SrcZ2.cc
-
-Copyright (C) 2016
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-
-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.
-*******************************************************************************/
-
-#include <Hadrons/Modules/SrcZ2.hpp>
-
-using namespace Grid;
-using namespace Hadrons;
-
-/******************************************************************************
-*                           SrcZ2 implementation                              *
-******************************************************************************/
-// constructor /////////////////////////////////////////////////////////////////
-SrcZ2::SrcZ2(const std::string name)
-: Module<SrcZ2Par>(name)
-{}
-
-// dependencies/products ///////////////////////////////////////////////////////
-std::vector<std::string> SrcZ2::getInput(void)
-{
-    std::vector<std::string> in;
-    
-    return in;
-}
-
-std::vector<std::string> SrcZ2::getOutput(void)
-{
-    std::vector<std::string> out = {getName()};
-    
-    return out;
-}
-
-// setup ///////////////////////////////////////////////////////////////////////
-void SrcZ2::setup(void)
-{
-    env().registerLattice<LatticePropagator>(getName());
-}
-
-// execution ///////////////////////////////////////////////////////////////////
-void SrcZ2::execute(void)
-{
-    Lattice<iScalar<vInteger>> t(env().getGrid());
-    LatticeComplex             eta(env().getGrid());
-    LatticeFermion             phi(env().getGrid());
-    Complex                    shift(1., 1.);
-    
-    if (par().tA == par().tB)
-    {
-        LOG(Message) << "Generating Z_2 wall source at t= " << par().tA
-                     << std::endl;
-    }
-    else
-    {
-        LOG(Message) << "Generating Z_2 band for " << par().tA << " <= t <= "
-                     << par().tB << std::endl;
-    }
-    LatticePropagator &src = *env().createLattice<LatticePropagator>(getName());
-    LatticeCoordinate(t, Tp);
-    bernoulli(*env().get4dRng(), eta);
-    eta = (2.*eta - shift)*(1./::sqrt(2.));
-    eta = where((t >= par().tA) and (t <= par().tB), eta, 0.*eta);
-    src = 1.;
-    src = src*eta;
-}

programs/Hadrons/Modules/SrcZ2.hpp

@@ -1,80 +0,0 @@
-/*******************************************************************************
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: programs/Hadrons/SrcZ2.hpp
-
-Copyright (C) 2016
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-
-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.
-*******************************************************************************/
-
-#ifndef Hadrons_SrcZ2_hpp_
-#define Hadrons_SrcZ2_hpp_
-
-#include <Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
-
-BEGIN_HADRONS_NAMESPACE
-
-/*
- 
- Z_2 stochastic source
- -----------------------------
- * src_x = eta_x * theta(x_3 - ta) * theta(tb - x_3)
- 
- * options:
- - tA: begin timeslice (integer)
- - tB: end timesilce (integer)
- 
- */
- 
-/******************************************************************************
- *                               SrcZ2                                        *
- ******************************************************************************/
-class SrcZ2Par: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(SrcZ2Par,
-                                    unsigned int, tA,
-                                    unsigned int, tB);
-};
-
-class SrcZ2: public Module<SrcZ2Par>
-{
-public:
-    // constructor
-    SrcZ2(const std::string name);
-    // destructor
-    virtual ~SrcZ2(void) = default;
-    // dependency relation
-    virtual std::vector<std::string> getInput(void);
-    virtual std::vector<std::string> getOutput(void);
-    // setup
-    virtual void setup(void);
-    // execution
-    virtual void execute(void);
-};
-
-MODULE_REGISTER(SrcZ2);
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_SrcZ2_hpp_

programs/Hadrons/add_module.sh

@@ -1,15 +0,0 @@
-#!/usr/bin/env bash
-
-if (( $# != 1 )); then
-    echo "usage: `basename $0` <module name>" 1>&2
-    exit 1
-fi
-NAME=$1
-
-if [ -e "Modules/${NAME}.cc" ] || [ -e "Modules/${NAME}.hpp" ]; then
-    echo "error: files Modules/${NAME}.* already exists" 1>&2
-    exit 1
-fi
-sed "s/___FILEBASENAME___/${NAME}/g" Module.cc.template > Modules/${NAME}.cc
-sed "s/___FILEBASENAME___/${NAME}/g" Module.hpp.template > Modules/${NAME}.hpp
-./make_module_list.sh

programs/Hadrons/make_module_list.sh

@@ -1,4 +0,0 @@
-#!/usr/bin/env bash
-
-echo 'modules =\' > modules.inc
-find Modules -name '*.cc' -type f -print | sed 's/^/  /;$q;s/$/ \\/' >> modules.inc