Merge branch 'master' of https://github.com/paboyle/Grid into scidac1_2

Conflicts:
	configure
	lib/qcd/action/Actions.h
	lib/qcd/action/fermion/WilsonKernels.h

lib/qcd/action/Actions.h

@@ -16,8 +16,6 @@
 #include <qcd/action/ActionBase.h>
 #include <qcd/action/ActionParams.h>
 
-
-
 ////////////////////////////////////////////
 // Utility functions
 ////////////////////////////////////////////
@@ -26,7 +24,6 @@
 #include <qcd/action/fermion/FermionOperator.h>
 #include <qcd/action/fermion/WilsonKernels.h>        //used by all wilson type fermions
 
-
 ////////////////////////////////////////////
 // Gauge Actions
 ////////////////////////////////////////////
@@ -55,12 +52,12 @@ typedef WilsonGaugeAction<LatticeGaugeFieldD>    WilsonGaugeActionD;
 
 #define FermOpTemplateInstantiate(A) \
   template class A<WilsonImplF>;		\
-  template class A<WilsonImplD>; 
-
-#define GparityFermOpTemplateInstantiate(A) \
-  template class A<GparityWilsonImplF>;	\
+  template class A<WilsonImplD>;		\
+  template class A<GparityWilsonImplF>;		\
   template class A<GparityWilsonImplD>;		
 
+#define GparityFermOpTemplateInstantiate(A) 
+
 ////////////////////////////////////////////
 // Fermion operators / actions
 ////////////////////////////////////////////

lib/qcd/action/fermion/FermionOperatorImpl.h

@@ -335,13 +335,33 @@ PARALLEL_FOR_LOOP
       }
 
       inline void InsertForce4D(GaugeField &mat, FermionField &Btilde, FermionField &A,int mu){
-	assert(0);
-	// Fixme
+	
+	// DhopDir provides U or Uconj depending on coor/flavour.
+	GaugeLinkField link(mat._grid);
+	// use lorentz for flavour as hack.
+	auto tmp = TraceIndex<SpinIndex>(outerProduct(Btilde,A));  
+PARALLEL_FOR_LOOP
+        for(auto ss=tmp.begin();ss<tmp.end();ss++){
+	  link[ss]() = tmp[ss](0,0) - conjugate(tmp[ss](1,1)) ;
+	}
+	PokeIndex<LorentzIndex>(mat,link,mu);
 	return;
       }
-      inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &A,int mu){
-	assert(0);
-	// Fixme
+      inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde,int mu){
+
+	int Ls=Btilde._grid->_fdimensions[0];
+
+	GaugeLinkField tmp(mat._grid);
+	tmp = zero;
+PARALLEL_FOR_LOOP
+	for(int ss=0;ss<tmp._grid->oSites();ss++){
+	  for(int s=0;s<Ls;s++){
+	    int sF = s+Ls*ss;
+	    auto ttmp = traceIndex<SpinIndex>(outerProduct(Btilde[sF],Atilde[sF]));
+	    tmp[ss]() = tmp[ss]()+ ttmp(0,0) + conjugate(ttmp(1,1));
+	  }
+	}
+	PokeIndex<LorentzIndex>(mat,tmp,mu);
 	return;
       }
     };

lib/qcd/action/fermion/WilsonFermion.cc

@@ -131,7 +131,7 @@ namespace QCD {
       // Flip gamma (1+g)<->(1-g) if dag
       ////////////////////////////////////////////////////////////////////////
       int gamma = mu;
-      if ( dag ) gamma+= Nd;
+      if ( !dag ) gamma+= Nd;
       
       ////////////////////////
       // Call the single hop
@@ -227,13 +227,15 @@ PARALLEL_FOR_LOOP
     DhopDir(in,out,dir,disp);
   }
   
+
   template<class Impl>
   void WilsonFermion<Impl>::DhopDir(const FermionField &in, FermionField &out,int dir,int disp){
     
     int skip = (disp==1) ? 0 : 1;
-    int dirdisp = dir+skip*4;
+    int dirdisp  = dir+skip*4;
+    int gamma    = dir+(1-skip)*4;
     
-    DhopDirDisp(in,out,dirdisp,dirdisp,DaggerNo);
+    DhopDirDisp(in,out,dirdisp,gamma,DaggerNo);
     
   };
   

lib/qcd/action/fermion/WilsonFermion5D.cc

@@ -94,6 +94,7 @@ void WilsonFermion5D<Impl>::DhopDir(const FermionField &in, FermionField &out,in
   int skip = (disp==1) ? 0 : 1;
 
   int dirdisp = dir+skip*4;
+  int gamma   = dir+(1-skip)*4;
 
   assert(dirdisp<=7);
   assert(dirdisp>=0);
@@ -103,7 +104,7 @@ PARALLEL_FOR_LOOP
     for(int s=0;s<Ls;s++){
       int sU=ss;
       int sF = s+Ls*sU; 
-      Kernels::DiracOptDhopDir(Stencil,Umu,comm_buf,sF,sU,in,out,dirdisp,dirdisp);
+      Kernels::DiracOptDhopDir(Stencil,Umu,comm_buf,sF,sU,in,out,dirdisp,gamma);
     }
   }
 };
@@ -136,7 +137,7 @@ void WilsonFermion5D<Impl>::DerivInternal(StencilImpl & st,
     // Flip gamma if dag
     ////////////////////////////////////////////////////////////////////////
     int gamma = mu;
-    if ( dag ) gamma+= Nd;
+    if ( !dag ) gamma+= Nd;
 
     ////////////////////////
     // Call the single hop

lib/qcd/action/fermion/WilsonKernels.h

@@ -41,7 +41,7 @@ namespace Grid {
 #endif
 // doesn't seem to work with Gparity at the moment
 #undef HANDOPT
-//#define HANDOPT
+#if 1
      void DiracOptHandDhopSite(StencilImpl &st,DoubledGaugeField &U,
 			       std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
 			       int sF,int sU,const FermionField &in, FermionField &out);
@@ -49,6 +49,7 @@ namespace Grid {
      void DiracOptHandDhopSiteDag(StencilImpl &st,DoubledGaugeField &U,
 				  std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
 				  int sF,int sU,const FermionField &in, FermionField &out);
+#endif
 
      WilsonKernels(const ImplParams &p= ImplParams());
      

lib/qcd/action/fermion/WilsonKernelsHand.cc

@@ -280,7 +280,7 @@
 namespace Grid {
 namespace QCD {
 
-#ifdef HANDOPT
+
 template<class Impl>
 void WilsonKernels<Impl >::DiracOptHandDhopSiteDag(StencilImpl &st,DoubledGaugeField &U,
 					       std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
@@ -767,24 +767,43 @@ void WilsonKernels<Impl >::DiracOptHandDhopSite(StencilImpl &st,DoubledGaugeFiel
     vstream(ref()(3)(2),result_32*(-0.5));
   }
 }
-#else 
-template<class Impl>
-void WilsonKernels<Impl >::DiracOptHandDhopSite(StencilImpl &st,DoubledGaugeField &U,
+
+  ////////////////////////////////////////////////
+  // Specialise Gparity to simple implementation
+  ////////////////////////////////////////////////
+template<>
+void WilsonKernels<GparityWilsonImplF>::DiracOptHandDhopSite(StencilImpl &st,DoubledGaugeField &U,
 					       std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
-					       int ss,int sU,const FermionField &in, FermionField &out)
+					       int sF,int sU,const FermionField &in, FermionField &out)
 {
-  DiracOptDhopSite(st,U,buf,ss,sU,in,out); // will template override for Wilson Nc=3
+  DiracOptDhopSite(st,U,buf,sF,sU,in,out); // will template override for Wilson Nc=3
 }
 
-template<class Impl>
-void WilsonKernels<Impl >::DiracOptHandDhopSiteDag(StencilImpl &st,DoubledGaugeField &U,
+template<>
+void WilsonKernels<GparityWilsonImplF>::DiracOptHandDhopSiteDag(StencilImpl &st,DoubledGaugeField &U,
 						   std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
-						   int ss,int sU,const FermionField &in, FermionField &out)
+						   int sF,int sU,const FermionField &in, FermionField &out)
 {
-  DiracOptDhopSiteDag(st,U,buf,ss,sU,in,out); // will template override for Wilson Nc=3
+  DiracOptDhopSiteDag(st,U,buf,sF,sU,in,out); // will template override for Wilson Nc=3
 }
 
-#endif
+template<>
+void WilsonKernels<GparityWilsonImplD>::DiracOptHandDhopSite(StencilImpl &st,DoubledGaugeField &U,
+					       std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
+					       int sF,int sU,const FermionField &in, FermionField &out)
+{
+  DiracOptDhopSite(st,U,buf,sF,sU,in,out); // will template override for Wilson Nc=3
+}
+
+template<>
+void WilsonKernels<GparityWilsonImplD>::DiracOptHandDhopSiteDag(StencilImpl &st,DoubledGaugeField &U,
+						   std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
+						   int sF,int sU,const FermionField &in, FermionField &out)
+{
+  DiracOptDhopSiteDag(st,U,buf,sF,sU,in,out); // will template override for Wilson Nc=3
+}
+
+
 
 template void WilsonKernels<WilsonImplF>::DiracOptHandDhopSite(StencilImpl &st,DoubledGaugeField &U,
 							  std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,