From eda4dd622e6470af3ee260a2d8568f8a2f40e985 Mon Sep 17 00:00:00 2001
From: Guido Cossu <guido.cossu@ed.ac.uk>
Date: Tue, 11 Oct 2016 15:45:20 +0100
Subject: [PATCH] Some more edit

---
 lib/parallelIO/BinaryIO.h      |  2 +-
 lib/qcd/QCD.h                  | 87 +++++++++++++++++++++++++++++-----
 lib/qcd/hmc/GenericHMCrunner.h | 11 +----
 lib/tensors/Tensor_outer.h     |  3 ++
 4 files changed, 82 insertions(+), 21 deletions(-)

diff --git a/lib/parallelIO/BinaryIO.h b/lib/parallelIO/BinaryIO.h
index f9a8ae3c..3265ba54 100644
--- a/lib/parallelIO/BinaryIO.h
+++ b/lib/parallelIO/BinaryIO.h
@@ -505,7 +505,7 @@ class BinaryIO {
     // available (how short sighted is that?)
     //////////////////////////////////////////////////////////
     Umu = zero;
-    static uint32_t csum=0;
+    static uint32_t csum; csum=0;//static for SHMEM
     fobj fileObj;
     static sobj siteObj; // Static to place in symmetric region for SHMEM
 
diff --git a/lib/qcd/QCD.h b/lib/qcd/QCD.h
index 74107a3e..61313f33 100644
--- a/lib/qcd/QCD.h
+++ b/lib/qcd/QCD.h
@@ -33,6 +33,71 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #define GRID_QCD_H
 namespace Grid{
 
+// First steps in the complete generalization of the Physics part
+namespace LatticeTheories {
+
+template <int Dimensions>
+struct LatticeTheory {
+  static const int Nd = Dimensions;
+  static const int Nds = Dimensions * 2;  // double stored field
+  template <typename vtype>
+  using iSinglet = iScalar<iScalar<iScalar<vtype> > >;
+};
+
+template <int Dimensions, int Colours>
+struct LatticeGaugeTheory : public LatticeTheory<Dimensions> {
+  static const int Nds = Dimensions * 2;
+  static const int Nd = Dimensions;
+  static const int Nc = Colours;
+
+  template <typename vtype>
+  using iColourMatrix = iScalar<iScalar<iMatrix<vtype, Nc> > >;
+  template <typename vtype>
+  using iLorentzColourMatrix = iVector<iScalar<iMatrix<vtype, Nc> >, Nd>;
+  template <typename vtype>
+  using iDoubleStoredColourMatrix = iVector<iScalar<iMatrix<vtype, Nc> >, Nds>;
+  template <typename vtype>
+  using iColourVector = iScalar<iScalar<iVector<vtype, Nc> > >;
+};
+
+template <int Dimensions, int Colours, int Spin>
+struct FermionicLatticeGaugeTheory
+    : public LatticeGaugeTheory<Dimensions, Colours> {
+  static const int Nd = Dimensions;
+  static const int Nds = Dimensions * 2;
+  static const int Nc = Colours;
+  static const int Ns = Spin;
+
+  template <typename vtype>
+  using iSpinMatrix = iScalar<iMatrix<iScalar<vtype>, Ns> >;
+  template <typename vtype>
+  using iSpinColourMatrix = iScalar<iMatrix<iMatrix<vtype, Nc>, Ns> >;
+  template <typename vtype>
+  using iSpinVector = iScalar<iVector<iScalar<vtype>, Ns> >;
+  template <typename vtype>
+  using iSpinColourVector = iScalar<iVector<iVector<vtype, Nc>, Ns> >;
+  // These 2 only is Spin is a multiple of 2
+  static const int Nhs = Spin / 2;
+  template <typename vtype>
+  using iHalfSpinVector = iScalar<iVector<iScalar<vtype>, Nhs> >;
+  template <typename vtype>
+  using iHalfSpinColourVector = iScalar<iVector<iVector<vtype, Nc>, Nhs> >;
+};
+
+struct QCD : public FermionicLatticeGaugeTheory<4, 3, 4> {
+    typedef FermionicLatticeGaugeTheory FLGT;
+    typedef FLGT::iSpinMatrix<Complex  >          SpinMatrix;
+    typedef FLGT::iSpinMatrix<ComplexF >          SpinMatrixF;
+    typedef FLGT::iSpinMatrix<ComplexD >          SpinMatrixD;
+
+};
+struct QED : public FermionicLatticeGaugeTheory<4, 1, 4> {};
+
+template <int Dimensions>
+struct Scalar : public LatticeTheory<Dimensions> {};
+
+}  // LatticeTheories
+
 namespace QCD {
 
 
@@ -355,36 +420,36 @@ namespace QCD {
     //////////////////////////////////////////////
     template<class vobj> 
       void pokeColour(Lattice<vobj> &lhs,
-		      const Lattice<decltype(peekIndex<ColourIndex>(lhs._odata[0],0))> & rhs,
-		      int i)
+              const Lattice<decltype(peekIndex<ColourIndex>(lhs._odata[0],0))> & rhs,
+              int i)
     {
       PokeIndex<ColourIndex>(lhs,rhs,i);
     }
     template<class vobj> 
       void pokeColour(Lattice<vobj> &lhs,
-		      const Lattice<decltype(peekIndex<ColourIndex>(lhs._odata[0],0,0))> & rhs,
-		      int i,int j)
+              const Lattice<decltype(peekIndex<ColourIndex>(lhs._odata[0],0,0))> & rhs,
+              int i,int j)
     {
       PokeIndex<ColourIndex>(lhs,rhs,i,j);
     }
     template<class vobj> 
       void pokeSpin(Lattice<vobj> &lhs,
-		      const Lattice<decltype(peekIndex<SpinIndex>(lhs._odata[0],0))> & rhs,
-		      int i)
+              const Lattice<decltype(peekIndex<SpinIndex>(lhs._odata[0],0))> & rhs,
+              int i)
     {
       PokeIndex<SpinIndex>(lhs,rhs,i);
     }
     template<class vobj> 
       void pokeSpin(Lattice<vobj> &lhs,
-		      const Lattice<decltype(peekIndex<SpinIndex>(lhs._odata[0],0,0))> & rhs,
-		      int i,int j)
+              const Lattice<decltype(peekIndex<SpinIndex>(lhs._odata[0],0,0))> & rhs,
+              int i,int j)
     {
       PokeIndex<SpinIndex>(lhs,rhs,i,j);
     }
     template<class vobj> 
       void pokeLorentz(Lattice<vobj> &lhs,
-		      const Lattice<decltype(peekIndex<LorentzIndex>(lhs._odata[0],0))> & rhs,
-		      int i)
+              const Lattice<decltype(peekIndex<LorentzIndex>(lhs._odata[0],0))> & rhs,
+              int i)
     {
       PokeIndex<LorentzIndex>(lhs,rhs,i);
     }
@@ -500,7 +565,7 @@ namespace QCD {
 #include <Grid/qcd/utils/LinalgUtils.h>
 #include <Grid/qcd/utils/CovariantCshift.h>
 
-// Include representations 	
+// Include representations  
 #include <Grid/qcd/utils/SUn.h>
 #include <Grid/qcd/utils/SUnAdjoint.h>
 #include <Grid/qcd/utils/SUnTwoIndex.h>
diff --git a/lib/qcd/hmc/GenericHMCrunner.h b/lib/qcd/hmc/GenericHMCrunner.h
index 8500cfb4..d627c47b 100644
--- a/lib/qcd/hmc/GenericHMCrunner.h
+++ b/lib/qcd/hmc/GenericHMCrunner.h
@@ -111,9 +111,11 @@ template <class IOCheckpointer = BinaryHmcCheckpointer<Implementation> >
     GridParallelRNG pRNG(UGrid);
     Field U(UGrid);
 
+    // This outside
     std::vector<int> SerSeed({1, 2, 3, 4, 5});
     std::vector<int> ParSeed({6, 7, 8, 9, 10});
 
+    // these decisions outside
     NoSmearing<Implementation> SmearingPolicy;
     typedef MinimumNorm2<Implementation, NoSmearing<Implementation>,
                          RepresentationsPolicy>
@@ -121,15 +123,6 @@ template <class IOCheckpointer = BinaryHmcCheckpointer<Implementation> >
     IntegratorParameters MDpar(20, 1.0);
     IntegratorType MDynamics(UGrid, MDpar, TheAction, SmearingPolicy);
 
-    // Checkpoint strategy   
-    /* 
-    int SaveInterval = 1;
-    std::string format = std::string("IEEE64BIG");
-    std::string conf_prefix = std::string("ckpoint_lat");
-    std::string rng_prefix = std::string("ckpoint_rng");
-    IOCheckpointer Checkpoint(conf_prefix, rng_prefix, SaveInterval, format);
-    */
-
     HMCparameters HMCpar;
     HMCpar.StartTrajectory = StartTraj;
     HMCpar.Trajectories = NumTraj;
diff --git a/lib/tensors/Tensor_outer.h b/lib/tensors/Tensor_outer.h
index 896fb8ed..6429a190 100644
--- a/lib/tensors/Tensor_outer.h
+++ b/lib/tensors/Tensor_outer.h
@@ -44,6 +44,8 @@ auto outerProduct (const iVector<l,N>& lhs,const iVector<r,N>& rhs) -> iMatrix<d
     }}
     return ret;
 }
+
+
 template<class l,class r> inline
 auto outerProduct (const iScalar<l>& lhs,const iScalar<r>& rhs) -> iScalar<decltype(outerProduct(lhs._internal,rhs._internal))>
 {
@@ -53,6 +55,7 @@ auto outerProduct (const iScalar<l>& lhs,const iScalar<r>& rhs) -> iScalar<declt
     return ret;
 }
 
+  
 inline ComplexF outerProduct(const ComplexF &l, const ComplexF& r)
 {
   std::cout << "outer product taking conj "<<r<<" "<<conj(r)<<std::endl;