From 77fb25fb292878f4bdeae5442f9270b6a38d2e46 Mon Sep 17 00:00:00 2001 From: Azusa Yamaguchi Date: Tue, 29 Nov 2016 13:43:56 +0000 Subject: [PATCH] Push 5d tests --- .../fermion/ImprovedStaggeredFermion5D.cc | 344 ++++++++++++++++++ .../fermion/ImprovedStaggeredFermion5D.h | 164 +++++++++ tests/core/Test_staggered.cc | 291 +++++++++++++++ tests/core/Test_staggered5D.cc | 314 ++++++++++++++++ 4 files changed, 1113 insertions(+) create mode 100644 lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc create mode 100644 lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h create mode 100644 tests/core/Test_staggered.cc create mode 100644 tests/core/Test_staggered5D.cc diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc new file mode 100644 index 00000000..71a6bf06 --- /dev/null +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc @@ -0,0 +1,344 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc + + Copyright (C) 2015 + +Author: Azusa Yamaguchi +Author: Peter Boyle + + 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 +#include + +namespace Grid { +namespace QCD { + +// S-direction is INNERMOST and takes no part in the parity. +const std::vector +ImprovedStaggeredFermion5DStatic::directions({1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4}); +const std::vector +ImprovedStaggeredFermion5DStatic::displacements({1, 1, 1, 1, -1, -1, -1, -1, 3, 3, 3, 3, -3, -3, -3, -3}); + + // 5d lattice for DWF. +template +ImprovedStaggeredFermion5D::ImprovedStaggeredFermion5D(GaugeField &_Uthin,GaugeField &_Ufat, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass, + RealD _c1,RealD _c2, RealD _u0, + const ImplParams &p) : + Kernels(p), + _FiveDimGrid (&FiveDimGrid), + _FiveDimRedBlackGrid(&FiveDimRedBlackGrid), + _FourDimGrid (&FourDimGrid), + _FourDimRedBlackGrid(&FourDimRedBlackGrid), + Stencil (_FiveDimGrid,npoint,Even,directions,displacements), + StencilEven(_FiveDimRedBlackGrid,npoint,Even,directions,displacements), // source is Even + StencilOdd (_FiveDimRedBlackGrid,npoint,Odd ,directions,displacements), // source is Odd + mass(_mass), + c1(_c1), + c2(_c2), + u0(_u0), + Umu(_FourDimGrid), + UmuEven(_FourDimRedBlackGrid), + UmuOdd (_FourDimRedBlackGrid), + UUUmu(_FourDimGrid), + UUUmuEven(_FourDimRedBlackGrid), + UUUmuOdd(_FourDimRedBlackGrid), + Lebesgue(_FourDimGrid), + LebesgueEvenOdd(_FourDimRedBlackGrid) +{ + // some assertions + assert(FiveDimGrid._ndimension==5); + assert(FourDimGrid._ndimension==4); + assert(FiveDimRedBlackGrid._ndimension==5); + assert(FourDimRedBlackGrid._ndimension==4); + assert(FiveDimRedBlackGrid._checker_dim==1); + + // Dimension zero of the five-d is the Ls direction + Ls=FiveDimGrid._fdimensions[0]; + assert(FiveDimRedBlackGrid._fdimensions[0]==Ls); + assert(FiveDimRedBlackGrid._processors[0] ==1); + assert(FiveDimRedBlackGrid._simd_layout[0]==1); + assert(FiveDimGrid._processors[0] ==1); + assert(FiveDimGrid._simd_layout[0] ==1); + + // Other dimensions must match the decomposition of the four-D fields + for(int d=0;d<4;d++){ + assert(FourDimRedBlackGrid._fdimensions[d] ==FourDimGrid._fdimensions[d]); + assert(FiveDimRedBlackGrid._fdimensions[d+1]==FourDimGrid._fdimensions[d]); + + assert(FourDimRedBlackGrid._processors[d] ==FourDimGrid._processors[d]); + assert(FiveDimRedBlackGrid._processors[d+1] ==FourDimGrid._processors[d]); + + assert(FourDimRedBlackGrid._simd_layout[d] ==FourDimGrid._simd_layout[d]); + assert(FiveDimRedBlackGrid._simd_layout[d+1]==FourDimGrid._simd_layout[d]); + + assert(FiveDimGrid._fdimensions[d+1] ==FourDimGrid._fdimensions[d]); + assert(FiveDimGrid._processors[d+1] ==FourDimGrid._processors[d]); + assert(FiveDimGrid._simd_layout[d+1] ==FourDimGrid._simd_layout[d]); + } + + // Allocate the required comms buffer + ImportGauge(_Uthin,_Ufat); +} + +template +void ImprovedStaggeredFermion5D::ImportGauge(const GaugeField &_Uthin) +{ + ImportGauge(_Uthin,_Uthin); +}; +template +void ImprovedStaggeredFermion5D::ImportGauge(const GaugeField &_Uthin,const GaugeField &_Ufat) +{ + GaugeLinkField U(GaugeGrid()); + + //////////////////////////////////////////////////////// + // Double Store should take two fields for Naik and one hop separately. + //////////////////////////////////////////////////////// + Impl::DoubleStore(GaugeGrid(), UUUmu, Umu, _Uthin, _Ufat ); + + //////////////////////////////////////////////////////// + // Apply scale factors to get the right fermion Kinetic term + // Could pass coeffs into the double store to save work. + // 0.5 ( U p(x+mu) - Udag(x-mu) p(x-mu) ) + //////////////////////////////////////////////////////// + for (int mu = 0; mu < Nd; mu++) { + + U = PeekIndex(Umu, mu); + PokeIndex(Umu, U*( 0.5*c1/u0), mu ); + + U = PeekIndex(Umu, mu+4); + PokeIndex(Umu, U*(-0.5*c1/u0), mu+4); + + U = PeekIndex(UUUmu, mu); + PokeIndex(UUUmu, U*( 0.5*c2/u0/u0/u0), mu ); + + U = PeekIndex(UUUmu, mu+4); + PokeIndex(UUUmu, U*(-0.5*c2/u0/u0/u0), mu+4); + } + + pickCheckerboard(Even, UmuEven, Umu); + pickCheckerboard(Odd, UmuOdd , Umu); + pickCheckerboard(Even, UUUmuEven, UUUmu); + pickCheckerboard(Odd, UUUmuOdd, UUUmu); +} +template +void ImprovedStaggeredFermion5D::DhopDir(const FermionField &in, FermionField &out,int dir5,int disp) +{ + int dir = dir5-1; // Maps to the ordering above in "directions" that is passed to stencil + // we drop off the innermost fifth dimension + + Compressor compressor; + Stencil.HaloExchange(in,compressor); + + PARALLEL_FOR_LOOP + for(int ss=0;ssoSites();ss++){ + for(int s=0;s +void ImprovedStaggeredFermion5D::DerivInternal(StencilImpl & st, + DoubledGaugeField & U, + DoubledGaugeField & UUU, + GaugeField &mat, + const FermionField &A, + const FermionField &B, + int dag) +{ + // No force terms in multi-rhs solver staggered + assert(0); +} + +template +void ImprovedStaggeredFermion5D::DhopDeriv(GaugeField &mat, + const FermionField &A, + const FermionField &B, + int dag) +{ + assert(0); +} + +template +void ImprovedStaggeredFermion5D::DhopDerivEO(GaugeField &mat, + const FermionField &A, + const FermionField &B, + int dag) +{ + assert(0); +} + + +template +void ImprovedStaggeredFermion5D::DhopDerivOE(GaugeField &mat, + const FermionField &A, + const FermionField &B, + int dag) +{ + assert(0); +} + +template +void ImprovedStaggeredFermion5D::DhopInternal(StencilImpl & st, LebesgueOrder &lo, + DoubledGaugeField & U,DoubledGaugeField & UUU, + const FermionField &in, FermionField &out,int dag) +{ + Compressor compressor; + + int LLs = in._grid->_rdimensions[0]; + + st.HaloExchange(in,compressor); + + // Dhop takes the 4d grid from U, and makes a 5d index for fermion + if (dag == DaggerYes) { + PARALLEL_FOR_LOOP + for (int ss = 0; ss < U._grid->oSites(); ss++) { + for(int s=0;soSites(); ss++) { + for(int s=0;s +void ImprovedStaggeredFermion5D::DhopOE(const FermionField &in, FermionField &out,int dag) +{ + conformable(in._grid,FermionRedBlackGrid()); // verifies half grid + conformable(in._grid,out._grid); // drops the cb check + + assert(in.checkerboard==Even); + out.checkerboard = Odd; + + DhopInternal(StencilEven,LebesgueEvenOdd,UmuOdd,UUUmuOdd,in,out,dag); +} +template +void ImprovedStaggeredFermion5D::DhopEO(const FermionField &in, FermionField &out,int dag) +{ + conformable(in._grid,FermionRedBlackGrid()); // verifies half grid + conformable(in._grid,out._grid); // drops the cb check + + assert(in.checkerboard==Odd); + out.checkerboard = Even; + + DhopInternal(StencilOdd,LebesgueEvenOdd,UmuEven,UUUmuEven,in,out,dag); +} +template +void ImprovedStaggeredFermion5D::Dhop(const FermionField &in, FermionField &out,int dag) +{ + conformable(in._grid,FermionGrid()); // verifies full grid + conformable(in._grid,out._grid); + + out.checkerboard = in.checkerboard; + + DhopInternal(Stencil,Lebesgue,Umu,UUUmu,in,out,dag); +} + + +///////////////////////////////////////////////////////////////////////// +// Implement the general interface. Here we use SAME mass on all slices +///////////////////////////////////////////////////////////////////////// +template +void ImprovedStaggeredFermion5D::Mdir(const FermionField &in, FermionField &out, int dir, int disp) { + DhopDir(in, out, dir, disp); +} +template +RealD ImprovedStaggeredFermion5D::M(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + Dhop(in, out, DaggerNo); + return axpy_norm(out, mass, in, out); +} + +template +RealD ImprovedStaggeredFermion5D::Mdag(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + Dhop(in, out, DaggerYes); + return axpy_norm(out, mass, in, out); +} + +template +void ImprovedStaggeredFermion5D::Meooe(const FermionField &in, FermionField &out) { + if (in.checkerboard == Odd) { + DhopEO(in, out, DaggerNo); + } else { + DhopOE(in, out, DaggerNo); + } +} +template +void ImprovedStaggeredFermion5D::MeooeDag(const FermionField &in, FermionField &out) { + if (in.checkerboard == Odd) { + DhopEO(in, out, DaggerYes); + } else { + DhopOE(in, out, DaggerYes); + } +} + +template +void ImprovedStaggeredFermion5D::Mooee(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + typename FermionField::scalar_type scal(mass); + out = scal * in; +} + +template +void ImprovedStaggeredFermion5D::MooeeDag(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + Mooee(in, out); +} + +template +void ImprovedStaggeredFermion5D::MooeeInv(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + out = (1.0 / (mass)) * in; +} + +template +void ImprovedStaggeredFermion5D::MooeeInvDag(const FermionField &in, + FermionField &out) { + out.checkerboard = in.checkerboard; + MooeeInv(in, out); +} + + + +FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion5D); + +}} + + + diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h new file mode 100644 index 00000000..c3502229 --- /dev/null +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h @@ -0,0 +1,164 @@ + + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h + + Copyright (C) 2015 + +Author: Peter Boyle +Author: AzusaYamaguchi + + 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 GRID_QCD_IMPROVED_STAGGERED_FERMION_5D_H +#define GRID_QCD_IMPROVED_STAGGERED_FERMION_5D_H + +namespace Grid { +namespace QCD { + + //////////////////////////////////////////////////////////////////////////////// + // This is the 4d red black case appropriate to support + //////////////////////////////////////////////////////////////////////////////// + + class ImprovedStaggeredFermion5DStatic { + public: + // S-direction is INNERMOST and takes no part in the parity. + static const std::vector directions; + static const std::vector displacements; + const int npoint = 16; + }; + + template + class ImprovedStaggeredFermion5D : public StaggeredKernels, public ImprovedStaggeredFermion5DStatic + { + public: + INHERIT_IMPL_TYPES(Impl); + typedef StaggeredKernels Kernels; + + /////////////////////////////////////////////////////////////// + // Implement the abstract base + /////////////////////////////////////////////////////////////// + GridBase *GaugeGrid(void) { return _FourDimGrid ;} + GridBase *GaugeRedBlackGrid(void) { return _FourDimRedBlackGrid ;} + GridBase *FermionGrid(void) { return _FiveDimGrid;} + GridBase *FermionRedBlackGrid(void) { return _FiveDimRedBlackGrid;} + + // full checkerboard operations; leave unimplemented as abstract for now + RealD M (const FermionField &in, FermionField &out); + RealD Mdag (const FermionField &in, FermionField &out); + + // half checkerboard operations + void Meooe (const FermionField &in, FermionField &out); + void Mooee (const FermionField &in, FermionField &out); + void MooeeInv (const FermionField &in, FermionField &out); + + void MeooeDag (const FermionField &in, FermionField &out); + void MooeeDag (const FermionField &in, FermionField &out); + void MooeeInvDag (const FermionField &in, FermionField &out); + + void Mdir (const FermionField &in, FermionField &out,int dir,int disp); + void DhopDir(const FermionField &in, FermionField &out,int dir,int disp); + + // These can be overridden by fancy 5d chiral action + void DhopDeriv (GaugeField &mat,const FermionField &U,const FermionField &V,int dag); + void DhopDerivEO(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); + void DhopDerivOE(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); + + // Implement hopping term non-hermitian hopping term; half cb or both + void Dhop (const FermionField &in, FermionField &out,int dag); + void DhopOE(const FermionField &in, FermionField &out,int dag); + void DhopEO(const FermionField &in, FermionField &out,int dag); + + + /////////////////////////////////////////////////////////////// + // New methods added + /////////////////////////////////////////////////////////////// + void DerivInternal(StencilImpl & st, + DoubledGaugeField & U, + DoubledGaugeField & UUU, + GaugeField &mat, + const FermionField &A, + const FermionField &B, + int dag); + + void DhopInternal(StencilImpl & st, + LebesgueOrder &lo, + DoubledGaugeField &U, + DoubledGaugeField &UUU, + const FermionField &in, + FermionField &out, + int dag); + + // Constructors + ImprovedStaggeredFermion5D(GaugeField &_Uthin, + GaugeField &_Ufat, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + double _mass, + RealD _c1=9.0/8.0, RealD _c2=-1.0/24.0,RealD _u0=1.0, + const ImplParams &p= ImplParams()); + + // DoubleStore + void ImportGauge(const GaugeField &_U); + void ImportGauge(const GaugeField &_Uthin,const GaugeField &_Ufat); + + /////////////////////////////////////////////////////////////// + // Data members require to support the functionality + /////////////////////////////////////////////////////////////// + public: + + GridBase *_FourDimGrid; + GridBase *_FourDimRedBlackGrid; + GridBase *_FiveDimGrid; + GridBase *_FiveDimRedBlackGrid; + + RealD mass; + RealD c1; + RealD c2; + RealD u0; + int Ls; + + //Defines the stencils for even and odd + StencilImpl Stencil; + StencilImpl StencilEven; + StencilImpl StencilOdd; + + // Copy of the gauge field , with even and odd subsets + DoubledGaugeField Umu; + DoubledGaugeField UmuEven; + DoubledGaugeField UmuOdd; + + DoubledGaugeField UUUmu; + DoubledGaugeField UUUmuEven; + DoubledGaugeField UUUmuOdd; + + LebesgueOrder Lebesgue; + LebesgueOrder LebesgueEvenOdd; + + // Comms buffer + std::vector > comm_buf; + + }; + +}} + +#endif diff --git a/tests/core/Test_staggered.cc b/tests/core/Test_staggered.cc new file mode 100644 index 00000000..89055fc7 --- /dev/null +++ b/tests/core/Test_staggered.cc @@ -0,0 +1,291 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./benchmarks/Benchmark_wilson.cc + + Copyright (C) 2015 + +Author: Peter Boyle +Author: paboyle + + 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 + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector 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< seeds({1,2,3,4}); + GridParallelRNG pRNG(&Grid); + pRNG.SeedFixedIntegers(seeds); + // pRNG.SeedRandomDevice(); + + typedef typename ImprovedStaggeredFermionR::FermionField FermionField; + typedef typename ImprovedStaggeredFermionR::ComplexField ComplexField; + 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; + FermionField phi (&Grid); random(pRNG,phi); + FermionField chi (&Grid); random(pRNG,chi); + LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu); + std::vector U(4,&Grid); + + + double volume=1; + for(int mu=0;mu(Umu,mu); + /* Debug force unit + U[mu] = 1.0; + PokeIndex(Umu,U[mu],mu); + */ + } + + ref = zero; + + RealD mass=0.1; + RealD c1=9.0/8.0; + RealD c2=-1.0/24.0; + RealD u0=1.0; + + { // Simple improved staggered implementation + ref = zero; + RealD c1tad = 0.5*c1/u0; + RealD c2tad = 0.5*c2/u0/u0/u0; + + Lattice > coor(&Grid); + + Lattice > x(&Grid); LatticeCoordinate(x,0); + Lattice > y(&Grid); LatticeCoordinate(y,1); + Lattice > z(&Grid); LatticeCoordinate(z,2); + Lattice > t(&Grid); LatticeCoordinate(t,3); + + Lattice > lin_z(&Grid); lin_z=x+y; + Lattice > lin_t(&Grid); lin_t=x+y+z; + + for(int mu=0;mu * = < chi | Deo^dag| phi> "< HermOpEO(Ds); + HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2); + HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2); + + HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2); + HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2); + + pDce = innerProduct(phi_e,dchi_e); + pDco = innerProduct(phi_o,dchi_o); + cDpe = innerProduct(chi_e,dphi_e); + cDpo = innerProduct(chi_o,dphi_o); + + std::cout< +Author: paboyle + + 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 + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + + std::cout << GridLogMessage << "Making s innermost grids"< seeds({1,2,3,4}); + GridParallelRNG pRNG4(UGrid); + GridParallelRNG pRNG5(FGrid); + pRNG4.SeedFixedIntegers(seeds); + pRNG5.SeedFixedIntegers(seeds); + + typedef typename ImprovedStaggeredFermion5DR::FermionField FermionField; + typedef typename ImprovedStaggeredFermion5DR::ComplexField ComplexField; + typename ImprovedStaggeredFermion5DR::ImplParams params; + + FermionField src (FGrid); random(pRNG5,src); + FermionField result(FGrid); result=zero; + FermionField ref(FGrid); ref=zero; + FermionField tmp(FGrid); tmp=zero; + FermionField err(FGrid); tmp=zero; + FermionField phi (FGrid); random(pRNG5,phi); + FermionField chi (FGrid); random(pRNG5,chi); + + LatticeGaugeField Umu(UGrid); SU3::ColdConfiguration(pRNG4,Umu); + + double volume=Ls; + for(int mu=0;muoSites();ss++){ + for(int s=0;s U(4,FGrid); + for(int mu=0;mu(Umu5d,mu); + if ( mu!=0 ) U[mu]=zero; + PokeIndex(Umu5d,U[mu],mu); + } + + std::vector Ua(4,UGrid); + for(int mu=0;mu(Umu,mu); + if ( mu!=0 ) { + Ua[mu]=zero; + } + PokeIndex(Umu,Ua[mu],mu); + } + + RealD mass=0.1; + RealD c1=9.0/8.0; + RealD c2=-1.0/24.0; + RealD u0=1.0; + + { // Simple improved staggered implementation + ref = zero; + RealD c1tad = 0.5*c1/u0; + RealD c2tad = 0.5*c2/u0/u0/u0; + + Lattice > coor(FGrid); + + Lattice > x(FGrid); LatticeCoordinate(x,1); // s innermost + Lattice > y(FGrid); LatticeCoordinate(y,2); + Lattice > z(FGrid); LatticeCoordinate(z,3); + Lattice > t(FGrid); LatticeCoordinate(t,4); + + Lattice > lin_z(FGrid); lin_z=x+y; + Lattice > lin_t(FGrid); lin_t=x+y+z; + + for(int mu=0;mu * = < chi | Deo^dag| phi> "< HermOpEO(Ds); + HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2); + HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2); + + HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2); + HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2); + + pDce = innerProduct(phi_e,dchi_e); + pDco = innerProduct(phi_o,dchi_o); + cDpe = innerProduct(chi_e,dphi_e); + cDpo = innerProduct(chi_o,dphi_o); + + std::cout<