Staggered phases first cut, c1, c2, u0

benchmarks/Benchmark_staggered.cc

@@ -0,0 +1,133 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./benchmarks/Benchmark_wilson.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);
+
+  std::vector<int> latt_size   = GridDefaultLatt();
+  std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  std::vector<int> mpi_layout  = GridDefaultMpi();
+  GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
+
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+  std::cout<<GridLogMessage << "Grid floating point word size is REALF"<< sizeof(RealF)<<std::endl;
+  std::cout<<GridLogMessage << "Grid floating point word size is REALD"<< sizeof(RealD)<<std::endl;
+  std::cout<<GridLogMessage << "Grid floating point word size is REAL"<< sizeof(Real)<<std::endl;
+
+  std::vector<int> seeds({1,2,3,4});
+  GridParallelRNG          pRNG(&Grid);
+  pRNG.SeedFixedIntegers(seeds);
+  //  pRNG.SeedRandomDevice();
+
+  typedef typename ImprovedStaggeredFermionR::FermionField FermionField; 
+  typename ImprovedStaggeredFermionR::ImplParams params; 
+
+  FermionField src   (&Grid); random(pRNG,src);
+  FermionField result(&Grid); result=zero;
+  FermionField    ref(&Grid);    ref=zero;
+  FermionField    tmp(&Grid);    tmp=zero;
+  FermionField    err(&Grid);    tmp=zero;
+  LatticeGaugeField Umu(&Grid); random(pRNG,Umu);
+  std::vector<LatticeColourMatrix> U(4,&Grid);
+
+  double volume=1;
+  for(int mu=0;mu<Nd;mu++){
+    volume=volume*latt_size[mu];
+  }  
+
+  // Only one non-zero (y)
+#if 0
+  Umu=zero;
+  Complex cone(1.0,0.0);
+  for(int nn=0;nn<Nd;nn++){
+    random(pRNG,U[nn]);
+    if(1) {
+      if (nn!=2) { U[nn]=zero; std::cout<<GridLogMessage << "zeroing gauge field in dir "<<nn<<std::endl; }
+      //      else       { U[nn]= cone;std::cout<<GridLogMessage << "unit gauge field in dir "<<nn<<std::endl; }
+      else       { std::cout<<GridLogMessage << "random gauge field in dir "<<nn<<std::endl; }
+    }
+    PokeIndex<LorentzIndex>(Umu,U[nn],nn);
+  }
+#endif
+
+  for(int mu=0;mu<Nd;mu++){
+    U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
+  }
+  ref = zero;
+  /*  
+  { // Naive wilson implementation
+    ref = zero;
+    for(int mu=0;mu<Nd;mu++){
+      //    ref =  src + Gamma(Gamma::GammaX)* src ; // 1-gamma_x
+      tmp = U[mu]*Cshift(src,mu,1);
+      for(int i=0;i<ref._odata.size();i++){
+	ref._odata[i]+= tmp._odata[i] - Gamma(Gmu[mu])*tmp._odata[i]; ;
+      }
+
+      tmp =adj(U[mu])*src;
+      tmp =Cshift(tmp,mu,-1);
+      for(int i=0;i<ref._odata.size();i++){
+	ref._odata[i]+= tmp._odata[i] + Gamma(Gmu[mu])*tmp._odata[i]; ;
+      }
+    }
+  }
+  ref = -0.5*ref;
+  */
+
+  RealD mass=0.1;
+  RealD c1=9.0/8.0;
+  RealD c2=-1.0/24.0;
+  ImprovedStaggeredFermionR Ds(Umu,Grid,RBGrid,mass,c1,c2,params);
+  
+  std::cout<<GridLogMessage << "Calling Ds"<<std::endl;
+  int ncall=1000;
+  double t0=usecond();
+  for(int i=0;i<ncall;i++){
+    Ds.Dhop(src,result,0);
+  }
+  double t1=usecond();
+  double flops=(16*(3*(6+8+8)) + 15*3*2)*volume*ncall; // == 66*16 +  == 1146
+  
+  std::cout<<GridLogMessage << "Called Ds"<<std::endl;
+  std::cout<<GridLogMessage << "norm result "<< norm2(result)<<std::endl;
+  std::cout<<GridLogMessage << "norm ref    "<< norm2(ref)<<std::endl;
+  std::cout<<GridLogMessage << "mflop/s =   "<< flops/(t1-t0)<<std::endl;
+  err = ref-result; 
+  std::cout<<GridLogMessage << "norm diff   "<< norm2(err)<<std::endl;
+
+  Grid_finalize();
+}

lib/qcd/action/Actions.h

@@ -260,6 +260,9 @@ typedef MobiusFermion<GparityWilsonImplR> GparityMobiusFermionR;
 typedef MobiusFermion<GparityWilsonImplF> GparityMobiusFermionF;
 typedef MobiusFermion<GparityWilsonImplD> GparityMobiusFermionD;
 
+typedef ImprovedStaggeredFermion<StaggeredImplR> ImprovedStaggeredFermionR;
+typedef ImprovedStaggeredFermion<StaggeredImplF> ImprovedStaggeredFermionF;
+typedef ImprovedStaggeredFermion<StaggeredImplD> ImprovedStaggeredFermionD;
 
   }}
 ///////////////////////////////////////////////////////////////////////////////

lib/qcd/action/fermion/FermionOperatorImpl.h

@@ -520,14 +520,17 @@ PARALLEL_FOR_LOOP
 
     INHERIT_GIMPL_TYPES(Gimpl);
       
+    template <typename vtype> using iImplScalar            = iScalar<iScalar<iScalar<vtype> > >;
     template <typename vtype> using iImplSpinor            = iScalar<iScalar<iVector<vtype, Dimension> > >;
     template <typename vtype> using iImplHalfSpinor        = iVector<iScalar<iVector<vtype, Dimension> >, Ngp>;
     template <typename vtype> using iImplDoubledGaugeField = iVector<iScalar<iMatrix<vtype, Dimension> >, Nds>;
     
+    typedef iImplScalar<Simd>            SiteComplex;
     typedef iImplSpinor<Simd>            SiteSpinor;
     typedef iImplHalfSpinor<Simd>        SiteHalfSpinor;
     typedef iImplDoubledGaugeField<Simd> SiteDoubledGaugeField;
     
+    typedef Lattice<SiteComplex>           ComplexField;
     typedef Lattice<SiteSpinor>            FermionField;
     typedef Lattice<SiteDoubledGaugeField> DoubledGaugeField;
     
@@ -564,15 +567,46 @@ PARALLEL_FOR_LOOP
       conformable(Uds._grid, GaugeGrid);
       conformable(Umu._grid, GaugeGrid);
       GaugeLinkField U(GaugeGrid);
+      GaugeLinkField UU(GaugeGrid);
+      GaugeLinkField UUU(GaugeGrid);
+      GaugeLinkField Udag(GaugeGrid);
+      GaugeLinkField UUUdag(GaugeGrid);
       for (int mu = 0; mu < Nd; mu++) {
-	U = PeekIndex<LorentzIndex>(Umu, mu);
+
+	// Staggered Phase.
+	ComplexField coor(GaugeGrid);
+	ComplexField phases(GaugeGrid);
+	ComplexField x(GaugeGrid); LatticeCoordinate(x,0);
+	ComplexField y(GaugeGrid); LatticeCoordinate(y,1);
+	ComplexField z(GaugeGrid); LatticeCoordinate(z,2);
+	ComplexField t(GaugeGrid); LatticeCoordinate(t,3);
+
+	SiteComplex zz(0.0,0.0);
+	SiteComplex one(1.0,0.0);
+
+	phases = one;
+	if ( mu == 1 ) phases = where( mod(x    ,2)== zz, phases,-phases);
+	if ( mu == 2 ) phases = where( mod(x+y  ,2)== zz, phases,-phases);
+	if ( mu == 3 ) phases = where( mod(x+y+z,2)== zz, phases,-phases);
+
+	U  = PeekIndex<LorentzIndex>(Umu, mu);
+	UU = Gimpl::CovShiftForward(U,mu,U);
+	UUU= Gimpl::CovShiftForward(U,mu,UU);
+	
+	U   = U   *phases;
+	UUU = UUU *phases;
+
 	PokeIndex<LorentzIndex>(Uds, U, mu);
-	PokeIndex<LorentzIndex>(UUUds, U, mu);
+	PokeIndex<LorentzIndex>(UUUds, UUU, mu);
-	std::cout << GridLogMessage << " NOT created the treble links for staggered yet" <<std::endl;
+
-	std::cout << GridLogMessage << " Must do this and also apply the staggered phases which requires understanding the action conventions" <<std::endl;
+	std::cout << GridLogMessage << " Created the treble links for staggered Naik term" <<std::endl;
-	U = adj(Cshift(U, mu, -1));
+	std::cout << GridLogMessage << " Multiplied the staggered phases which requires understanding the action conventions" <<std::endl;
-	PokeIndex<LorentzIndex>(Uds, U, mu + 4);
+
-	PokeIndex<LorentzIndex>(UUUds, U, mu+4);
+	Udag   = adj( Cshift(U, mu, -1));
+	UUUdag = adj( Cshift(UUU, mu, -3));
+	PokeIndex<LorentzIndex>(Uds, Udag, mu + 4);
+	PokeIndex<LorentzIndex>(UUUds, UUUdag, mu+4);
+
       }
     }
 

lib/qcd/action/fermion/ImprovedStaggeredFermion.cc

@@ -43,6 +43,7 @@ ImprovedStaggeredFermionStatic::displacements({1, 1, 1, 1, -1, -1, -1, -1, 3, 3,
 template <class Impl>
 ImprovedStaggeredFermion<Impl>::ImprovedStaggeredFermion(GaugeField &_Umu, GridCartesian &Fgrid,
 							 GridRedBlackCartesian &Hgrid, RealD _mass,
+							 RealD _c1, RealD _c2,RealD _u0,
 							 const ImplParams &p)
     : Kernels(p),
       _grid(&Fgrid),
@@ -51,6 +52,9 @@ ImprovedStaggeredFermion<Impl>::ImprovedStaggeredFermion(GaugeField &_Umu, GridC
       StencilEven(&Hgrid, npoint, Even, directions, displacements),  // source is Even
       StencilOdd(&Hgrid, npoint, Odd, directions, displacements),  // source is Odd
       mass(_mass),
+      c1(_c1),
+      c2(_c2),
+      u0(_u0),
       Lebesgue(_grid),
       LebesgueEvenOdd(_cbgrid),
       Umu(&Fgrid),
@@ -63,15 +67,52 @@ ImprovedStaggeredFermion<Impl>::ImprovedStaggeredFermion(GaugeField &_Umu, GridC
   ImportGauge(_Umu);
 }
 
+  ////////////////////////////////////////////////////////////
+  // Momentum space propagator should be 
+  // https://arxiv.org/pdf/hep-lat/9712010.pdf
+  //
+  // mom space action.
+  //   gamma_mu i ( c1 sin pmu + c2 sin 3 pmu ) + m
+  //
+  // must track through staggered flavour/spin reduction in literature to 
+  // turn to free propagator for the one component chi field, a la page 4/5
+  // of above link to implmement fourier based solver.
+  ////////////////////////////////////////////////////////////
 template <class Impl>
 void ImprovedStaggeredFermion<Impl>::ImportGauge(const GaugeField &_Umu) {
-  GaugeField HUmu(_Umu._grid);
+
-  HUmu = _Umu * (-0.5);
+  GaugeLinkField U(GaugeGrid);
-  Impl::DoubleStore(GaugeGrid(), Umu, UUUmu, HUmu);
+
+  ////////////////////////////////////////////////////////
+  // Double Store should take two fields for Naik and one hop separately.
+  ////////////////////////////////////////////////////////
+  Impl::DoubleStore(GaugeGrid(), Umu, UUUmu, _Umu);
+
+
+  ////////////////////////////////////////////////////////
+  // Apply scale factors to get the right fermion Kinetic term
+  // 
+  // 0.5 ( U p(x+mu) - Udag(x-mu) p(x-mu) ) 
+  ////////////////////////////////////////////////////////
+  for (int mu = 0; mu < Nd; mu++) {
+
+    U = PeekIndex<LorentzIndex>(Umu, mu);
+    PokeIndex<LorentzIndex>(Umu, U*( 0.5*c1/u0), mu );
+    
+    U = PeekIndex<LorentzIndex>(Umu, mu+4);
+    PokeIndex<LorentzIndex>(Umu, U*(-0.5*c1/u0), mu+4);
+
+    U = PeekIndex<LorentzIndex>(UUUmu, mu);
+    PokeIndex<LorentzIndex>(UUUmu, U*( 0.5*c2/u0/u0/u0), mu );
+    
+    U = PeekIndex<LorentzIndex>(UUUmu, mu+4);
+    PokeIndex<LorentzIndex>(UUUmu, U*(-0.5*c2/u0/u0/u0), mu+4);
+  }
+
   pickCheckerboard(Even, UmuEven, Umu);
-  pickCheckerboard(Odd, UmuOdd, Umu);
+  pickCheckerboard(Odd,  UmuOdd , Umu);
   pickCheckerboard(Even, UUUmuEven, UUUmu);
-  pickCheckerboard(Odd, UUUmuOdd, UUUmu);
+  pickCheckerboard(Odd,   UUUmuOdd, UUUmu);
 }
 
 /////////////////////////////

lib/qcd/action/fermion/ImprovedStaggeredFermion.h

@@ -108,6 +108,7 @@ class ImprovedStaggeredFermion : public StaggeredKernels<Impl>, public ImprovedS
   // Constructor
   ImprovedStaggeredFermion(GaugeField &_Umu, GridCartesian &Fgrid,
 			   GridRedBlackCartesian &Hgrid, RealD _mass,
+			   RealD _c1=9.0/8.0, RealD _c2=-1.0/24.0,RealD _u0,
 			   const ImplParams &p = ImplParams());
 
   // DoubleStore impl dependent
@@ -122,6 +123,9 @@ class ImprovedStaggeredFermion : public StaggeredKernels<Impl>, public ImprovedS
   // any other parameters of action ???
 
   RealD mass;
+  RealD u0;
+  RealD c1;
+  RealD c2;
 
   GridBase *_grid;
   GridBase *_cbgrid;