From 63fefb8ea580d7b09c661b8f4611cdd469b524a1 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Wed, 19 Aug 2026 20:33:10 -0400 Subject: [PATCH] new tests for lex lattice --- tests/lexLattice/Makefile.am | 1 + tests/lexLattice/Test_GaugeAction_lex.cc | 128 +++++++++ tests/lexLattice/Test_RectPlaq_lex.cc | 216 +++++++++++++++ tests/lexLattice/Test_cshift_lex.cc | 175 +++++++++++++ tests/lexLattice/Test_dwf_lex.cc | 196 ++++++++++++++ .../lexLattice/Test_innerproduct_norm_lex.cc | 132 ++++++++++ tests/lexLattice/Test_io.cc | 194 ++++++++++++++ tests/lexLattice/Test_staggered_lex.cc | 246 ++++++++++++++++++ tests/lexLattice/Test_stencil_lex.cc | 122 +++++++++ tests/lexLattice/Test_wilson_lex.cc | 239 +++++++++++++++++ 10 files changed, 1649 insertions(+) create mode 100644 tests/lexLattice/Makefile.am create mode 100644 tests/lexLattice/Test_GaugeAction_lex.cc create mode 100644 tests/lexLattice/Test_RectPlaq_lex.cc create mode 100644 tests/lexLattice/Test_cshift_lex.cc create mode 100644 tests/lexLattice/Test_dwf_lex.cc create mode 100644 tests/lexLattice/Test_innerproduct_norm_lex.cc create mode 100644 tests/lexLattice/Test_io.cc create mode 100644 tests/lexLattice/Test_staggered_lex.cc create mode 100644 tests/lexLattice/Test_stencil_lex.cc create mode 100644 tests/lexLattice/Test_wilson_lex.cc diff --git a/tests/lexLattice/Makefile.am b/tests/lexLattice/Makefile.am new file mode 100644 index 000000000..60b82dd78 --- /dev/null +++ b/tests/lexLattice/Makefile.am @@ -0,0 +1 @@ +include Make.inc diff --git a/tests/lexLattice/Test_GaugeAction_lex.cc b/tests/lexLattice/Test_GaugeAction_lex.cc new file mode 100644 index 000000000..0b37b1e14 --- /dev/null +++ b/tests/lexLattice/Test_GaugeAction_lex.cc @@ -0,0 +1,128 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/lexLattice/Test_GaugeAction_lex.cc + + Copyright (C) 2026 + +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 */ + +// +// Wilson loops on a lexicographic gauge field. +// +// The configuration is generated on the vectorised grid, written as NERSC, +// and read back into the lexicographic layout: the reader recomputes the +// plaquette and link trace and checks them against the header written by +// the vectorised chart, so the crossing is validated inside Grid's own IO. +// Plaquette, link trace, rectangle and the staple identity are then compared +// between the two charts. +// +#include + +using namespace Grid; + +const RealD tol = 1.0e-10; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + typedef WilsonLoops vWL; + typedef WilsonLoops lexWL; + + Coordinate latt = GridDefaultLatt(); + Coordinate vsimd = GridDefaultSimd(Nd,vComplexD::Nsimd()); + Coordinate lsimd({1,1,1,1}); + Coordinate mpi = GridDefaultMpi(); + + GridCartesian vGrid(latt,vsimd,mpi); + GridCartesian lGrid(latt,lsimd,mpi); + + std::cout << GridLogMessage << "vectorised Nsimd = " << vGrid.Nsimd() << std::endl; + std::cout << GridLogMessage << "lexicographic Nsimd = " << lGrid.Nsimd() << std::endl; + GRID_ASSERT( lGrid.Nsimd() == 1 ); + + GridParallelRNG pRNG(&vGrid); + pRNG.SeedFixedIntegers(std::vector({1,2,3,4})); + + LatticeGaugeFieldD Umu(&vGrid); + SU::HotConfiguration(pRNG,Umu); + + std::string file("./ckpoint_lex.4000"); + NerscIO::writeConfiguration(Umu,file,0,0); + + lexLatticeGaugeFieldD Ulex(&lGrid); + FieldMetaData header; + NerscIO::readConfiguration(Ulex,header,file); + + ////////////////////////////////////////////////// + // Plaquette, link trace, rectangle in both charts + ////////////////////////////////////////////////// + RealD vplaq = vWL::avgPlaquette(Umu); + RealD lplaq = lexWL::avgPlaquette(Ulex); + + RealD vlink = vWL::linkTrace(Umu); + RealD llink = lexWL::linkTrace(Ulex); + + RealD vrect = vWL::avgRectangle(Umu); + RealD lrect = lexWL::avgRectangle(Ulex); + + std::cout << GridLogMessage << "plaquette simd " << vplaq << " lex " << lplaq + << " header " << header.plaquette << std::endl; + std::cout << GridLogMessage << "link trace simd " << vlink << " lex " << llink + << " header " << header.link_trace << std::endl; + std::cout << GridLogMessage << "rectangle simd " << vrect << " lex " << lrect << std::endl; + + GRID_ASSERT( fabs(vplaq-lplaq) < tol ); + GRID_ASSERT( fabs(vlink-llink) < tol ); + GRID_ASSERT( fabs(vrect-lrect) < tol ); + + ////////////////////////////////////////////////// + // Plaquette via staples, in the lexicographic chart + ////////////////////////////////////////////////// + { + RealD vol = lGrid.gSites(); + RealD stap_plaq = 0.0; + + lexLatticeColourMatrixD stap(&lGrid); + lexLatticeColourMatrixD Ul(&lGrid); + lexLatticeComplexD stap_tr(&lGrid); + + for(int mu=0;mu(Ulex,mu); + lexWL::Staple(stap,Ulex,mu); + stap_tr = trace(Ul*stap); + auto Ts = sum(stap_tr); + stap_plaq += real(TensorRemove(Ts)); + } + RealD StapScale = 1.0/vol/6.0/Nc/4.0; + RealD plaq_from_staples = stap_plaq*StapScale; + + std::cout << GridLogMessage << "plaquette via staples (lex) " << plaq_from_staples + << " direct " << lplaq << std::endl; + GRID_ASSERT( fabs(plaq_from_staples - lplaq) < 1.0e-8 ); + } + + std::cout << GridLogMessage << "Test_GaugeAction_lex: ALL PASS" << std::endl; + + Grid_finalize(); +} diff --git a/tests/lexLattice/Test_RectPlaq_lex.cc b/tests/lexLattice/Test_RectPlaq_lex.cc new file mode 100644 index 000000000..c52d269ee --- /dev/null +++ b/tests/lexLattice/Test_RectPlaq_lex.cc @@ -0,0 +1,216 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/lexLattice/Test_RectPlaq_lex.cc + + Copyright (C) 2026 + +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 */ + +// +// Plaquette and 2x1 rectangle built from covariant shifts, checked against +// WilsonLoops, plus link trace and a blockSum coarsening. +// +// The measurement is written once against a gauge implementation and run in +// both the vectorised and the lexicographic chart; every number must agree. +// The configuration crosses layouts through NERSC IO, whose reader verifies +// the header written by the other chart. +// +#include + +using namespace Grid; + +const RealD tol = 1.0e-10; + +template +struct Measurements +{ + RealD plaq_shift; + RealD plaq_loops; + RealD plaq_staple; + RealD rect_shift; + RealD rect_loops; + RealD link; + RealD coarse_plaq; +}; + +template +void Measure(typename Gimpl::Field &Umu,GridBase *coarse,Measurements &m) +{ + typedef typename Gimpl::LinkField LinkField; + typedef typename Gimpl::ComplexField ComplexField; + typedef WilsonLoops WL; + + GridBase *grid = Umu.Grid(); + RealD vol = grid->gSites(); + + std::vector U(Nd,grid); + for(int mu=0;mu(Umu,mu); + } + + /////////////////////////////////////////////////// + // Link trace + /////////////////////////////////////////////////// + ComplexField LinkTrace(grid); + LinkTrace = Zero(); + for(int mu=0;mu +void Report(std::string name,Measurements &m) +{ + std::cout << GridLogMessage << name << ": plaquette shifts " << m.plaq_shift + << " loops " << m.plaq_loops << " staples " << m.plaq_staple << std::endl; + std::cout << GridLogMessage << name << ": rectangle shifts " << m.rect_shift + << " loops " << m.rect_loops << std::endl; + std::cout << GridLogMessage << name << ": link trace " << m.link + << " coarsened plaquette " << m.coarse_plaq << std::endl; +} + +template +void SelfConsistent(Measurements &m) +{ + GRID_ASSERT( fabs(m.plaq_shift - m.plaq_loops) < 1.0e-8 ); + GRID_ASSERT( fabs(m.plaq_staple - m.plaq_loops) < 1.0e-8 ); + GRID_ASSERT( fabs(m.rect_shift - m.rect_loops) < 1.0e-8 ); + GRID_ASSERT( fabs(m.coarse_plaq - m.plaq_shift) < 1.0e-8 ); +} + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + Coordinate latt = GridDefaultLatt(); + Coordinate vsimd = GridDefaultSimd(Nd,vComplexD::Nsimd()); + Coordinate lsimd({1,1,1,1}); + Coordinate mpi = GridDefaultMpi(); + + Coordinate clatt(Nd); + for(int d=0;d vm; + Measurements lm; + + Measure (Umu ,&vCoarse,vm); + Measure(Ulex,&lCoarse,lm); + + Report("simd",vm); + Report("lex ",lm); + + SelfConsistent(vm); + SelfConsistent(lm); + + GRID_ASSERT( fabs(vm.plaq_shift - lm.plaq_shift ) < tol ); + GRID_ASSERT( fabs(vm.plaq_loops - lm.plaq_loops ) < tol ); + GRID_ASSERT( fabs(vm.plaq_staple - lm.plaq_staple) < tol ); + GRID_ASSERT( fabs(vm.rect_shift - lm.rect_shift ) < tol ); + GRID_ASSERT( fabs(vm.rect_loops - lm.rect_loops ) < tol ); + GRID_ASSERT( fabs(vm.link - lm.link ) < tol ); + GRID_ASSERT( fabs(vm.coarse_plaq - lm.coarse_plaq) < tol ); + + std::cout << GridLogMessage << "Test_RectPlaq_lex: ALL PASS" << std::endl; + + Grid_finalize(); +} diff --git a/tests/lexLattice/Test_cshift_lex.cc b/tests/lexLattice/Test_cshift_lex.cc new file mode 100644 index 000000000..cf98af2d9 --- /dev/null +++ b/tests/lexLattice/Test_cshift_lex.cc @@ -0,0 +1,175 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/lexLattice/Test_cshift_lex.cc + + Copyright (C) 2026 + +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 */ + +// +// Lexicographic (Nsimd()==1) lattice types: LatticeCoordinate, Cshift in all +// directions and shifts, and predicated where(). +// +// Cshift is checked by an identity rather than a reference loop: with +// d = Cshift(coor_mu,mu,shift) - coor_mu +// every site has d = shift (no wrap) or d = shift - L (wrapped), so +// (d - shift)*(d - shift + L) == 0 +// everywhere. +// +#include + +using namespace Grid; + +const RealD tol = 1.0e-5; + +int main(int argc, char** argv) +{ + Grid_init(&argc, &argv); + + Coordinate latt = GridDefaultLatt(); + Coordinate simd({1,1,1,1}); // lexicographic: no SIMD + Coordinate mpi = GridDefaultMpi(); + + GridCartesian grid(latt,simd,mpi); + + std::cout << GridLogMessage << "Lexicographic grid, Nsimd = " << grid.Nsimd() << std::endl; + GRID_ASSERT( grid.Nsimd() == 1 ); + + + + ////////////////////////////////////////////////// + // Cshift against the wrap identity + ////////////////////////////////////////////////// + for(int mu=0;mu tol ) { + std::cout << GridLogMessage << "FAIL mu " << mu << " shift " << shift + << " residual " << n << std::endl; + } + GRID_ASSERT( n < tol ); + } + std::cout << GridLogMessage << "Cshift mu = " << mu << " all shifts pass" << std::endl; + } + + ////////////////////////////////////////////////// + // Predicated where(): zero the second half in time + ////////////////////////////////////////////////// + { + int Tdir = Nd-1; + RealD T = latt[Tdir]; + RealD vol = grid.gSites(); + + lexLatticeInteger tcoor(&grid); + LatticeCoordinate(tcoor,Tdir); + + lexLatticeComplexD f(&grid); f = 2.0; + lexLatticeComplexD zz(&grid); zz = Zero(); + lexLatticeComplexD g(&grid); + + g = where( tcoor < Integer(T/2) , f , zz ); + + RealD n_g = norm2(g); + RealD expect = 4.0 * vol / 2.0; + + std::cout << GridLogMessage << "where(): norm2 = " << n_g + << " expect " << expect << std::endl; + DumpSliceNorm("where() slice norm",g,Tdir); + + GRID_ASSERT( fabs(n_g - expect) < tol*expect ); + } + + ////////////////////////////////////////////////// + // Cshift against sliceSum of a slice-zeroed random field. + // Exhaustive over direction, zeroed slice and shift: + // sliceSum(Cshift(f,mu,s))[t] == sliceSum(f)[(t+s)%L] + ////////////////////////////////////////////////// + { + GridParallelRNG pRNG(&grid); + pRNG.SeedFixedIntegers(std::vector({7,8,9,10})); + + typedef lexLatticeComplexD::vector_object::scalar_object sobj; + + lexLatticeComplexD rnd(&grid); random(pRNG,rnd); + lexLatticeComplexD zz(&grid); zz = Zero(); + + for(int mu=0;mu ref; + sliceSum(f,ref,mu); + + for(int shift=0;shift res; + sliceSum(g,res,mu); + + for(int t=0;t tol*(abs(expect)+1.0) ) { + std::cout << GridLogMessage << "FAIL mu " << mu << " t0 " << t0 + << " shift " << shift << " t " << t + << " got " << got << " expect " << expect << std::endl; + } + GRID_ASSERT( abs(got-expect) < tol*(abs(expect)+1.0) ); + } + } + } + std::cout << GridLogMessage << "sliceSum/Cshift consistency mu = " << mu + << " all zeroed slices and shifts pass" << std::endl; + } + } + + std::cout << GridLogMessage << "Test_cshift_lex: ALL PASS" << std::endl; + + Grid_finalize(); +} diff --git a/tests/lexLattice/Test_dwf_lex.cc b/tests/lexLattice/Test_dwf_lex.cc new file mode 100644 index 000000000..03dea9d14 --- /dev/null +++ b/tests/lexLattice/Test_dwf_lex.cc @@ -0,0 +1,196 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/lexLattice/Test_dwf_lex.cc + + Copyright (C) 2026 + +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 */ + +// +// Five dimensional fermions on a lexicographic (Nsimd()==1) lattice. +// +// Domain wall and Mobius, driven from one gauge field and one source in both +// charts and compared elementwise through the layout interchange. Mobius uses +// b=1.5 c=0.5, the production choice. +// +#include + +using namespace Grid; + +const RealD tol = 1.0e-10; + +//////////////////////////////////////////////////////////////////////// +// Layout interchange; both lattices share the same scalar_object. +//////////////////////////////////////////////////////////////////////// +template +void transfer(Lattice &out,const Lattice &in) +{ + typedef typename vobjIn::scalar_object sobj; + static_assert(std::is_same::value, + "transfer: lattices must share scalar_object"); + + GRID_ASSERT(out.Grid()->gSites() == in.Grid()->gSites()); + + std::vector buf; + unvectorizeToLexOrdArray(buf,in); + vectorizeFromLexOrdArray(buf,out); +} + +//////////////////////////////////////////////////////////////////////// +// Compare one operator application between the two charts +//////////////////////////////////////////////////////////////////////// +template +void Compare(std::string name,int dag, + vOp &vD, typename vOp::FermionField &vsrc, + lexOp &lD, typename lexOp::FermionField &lsrc) +{ + typename vOp::FermionField vres(vsrc.Grid()); + typename lexOp::FermionField lres(lsrc.Grid()); + typename lexOp::FermionField vres_lex(lsrc.Grid()); + typename lexOp::FermionField err(lsrc.Grid()); + + if ( dag == DaggerNo ) { vD.M (vsrc,vres); lD.M (lsrc,lres); } + else { vD.Mdag(vsrc,vres); lD.Mdag(lsrc,lres); } + + transfer(vres_lex,vres); + err = vres_lex - lres; + RealD n = norm2(err)/norm2(lres); + std::cout << GridLogMessage << name << " simd vs lex: relative " << n + << " |simd|^2 " << norm2(vres) << " |lex|^2 " << norm2(lres) << std::endl; + GRID_ASSERT( n < tol ); +} + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + const int Ls = 8; + + Coordinate latt = GridDefaultLatt(); + Coordinate vsimd = GridDefaultSimd(Nd,vComplexD::Nsimd()); + Coordinate lsimd({1,1,1,1}); + Coordinate mpi = GridDefaultMpi(); + + GridCartesian *vUGrid = new GridCartesian(latt,vsimd,mpi); + GridRedBlackCartesian *vUrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(vUGrid); + GridCartesian *vFGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,vUGrid); + GridRedBlackCartesian *vFrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,vUGrid); + + GridCartesian *lUGrid = new GridCartesian(latt,lsimd,mpi); + GridRedBlackCartesian *lUrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(lUGrid); + GridCartesian *lFGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,lUGrid); + GridRedBlackCartesian *lFrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,lUGrid); + + std::cout << GridLogMessage << "Ls = " << Ls + << " vectorised Nsimd(4d) = " << vUGrid->Nsimd() + << " Nsimd(5d) = " << vFGrid->Nsimd() << std::endl; + std::cout << GridLogMessage << " lexicographic Nsimd(4d) = " << lUGrid->Nsimd() + << " Nsimd(5d) = " << lFGrid->Nsimd() << std::endl; + GRID_ASSERT( lUGrid->Nsimd() == 1 ); + GRID_ASSERT( lFGrid->Nsimd() == 1 ); + + RealD mass = 0.05; + RealD M5 = 1.8; + + GridParallelRNG pRNG4(vUGrid); pRNG4.SeedFixedIntegers(std::vector({1,2,3,4})); + GridParallelRNG pRNG5(vFGrid); pRNG5.SeedFixedIntegers(std::vector({5,6,7,8})); + + ////////////////////////////////////////////////// + // One gauge field and one source, both charts + ////////////////////////////////////////////////// + LatticeGaugeFieldD Umu(vUGrid); SU::HotConfiguration(pRNG4,Umu); + lexLatticeGaugeFieldD Ulex(lUGrid); transfer(Ulex,Umu); + + typedef DomainWallFermion vDwf; + typedef DomainWallFermion lexDwf; + typedef MobiusFermion vMob; + typedef MobiusFermion lexMob; + + typename vDwf::FermionField src(vFGrid); random(pRNG5,src); + typename vDwf::FermionField phi(vFGrid); random(pRNG5,phi); + typename lexDwf::FermionField srclex(lFGrid); transfer(srclex,src); + typename lexDwf::FermionField philex(lFGrid); transfer(philex,phi); + + ////////////////////////////////////////////////// + // Domain wall + ////////////////////////////////////////////////// + { + vDwf Dv (Umu ,*vFGrid,*vFrbGrid,*vUGrid,*vUrbGrid,mass,M5); + lexDwf Dl (Ulex,*lFGrid,*lFrbGrid,*lUGrid,*lUrbGrid,mass,M5); + + Compare("DWF M ",DaggerNo ,Dv,src,Dl,srclex); + Compare("DWF Mdag",DaggerYes,Dv,src,Dl,srclex); + + // Adjoint identity within the lexicographic chart + typename lexDwf::FermionField Mphi(lFGrid),Mdagchi(lFGrid); + Dl.M (srclex,Mphi); + Dl.Mdag(philex,Mdagchi); + ComplexD lhs = innerProduct(philex,Mphi); + ComplexD rhs = innerProduct(Mdagchi,srclex); + RealD n = abs(lhs-rhs)/abs(lhs); + std::cout << GridLogMessage << "DWF lex = " << lhs + << " = " << rhs << " relative " << n << std::endl; + GRID_ASSERT( n < tol ); + + // Checkerboarded hopping term + typename vDwf::FermionField vo(vFrbGrid),ve(vFrbGrid); + typename lexDwf::FermionField lo(lFrbGrid),le(lFrbGrid); + vo.Checkerboard()=Odd; ve.Checkerboard()=Even; + lo.Checkerboard()=Odd; le.Checkerboard()=Even; + pickCheckerboard(Odd,vo,src); + pickCheckerboard(Odd,lo,srclex); + Dv.Meooe(vo,ve); + Dl.Meooe(lo,le); + RealD vn = norm2(ve), ln = norm2(le); + std::cout << GridLogMessage << "DWF Meooe |simd|^2 " << vn << " |lex|^2 " << ln + << " relative " << fabs(vn-ln)/vn << std::endl; + GRID_ASSERT( fabs(vn-ln)/vn < tol ); + } + + ////////////////////////////////////////////////// + // Mobius, production coefficients + ////////////////////////////////////////////////// + { + RealD b=1.5, c=0.5; + vMob Dv (Umu ,*vFGrid,*vFrbGrid,*vUGrid,*vUrbGrid,mass,M5,b,c); + lexMob Dl (Ulex,*lFGrid,*lFrbGrid,*lUGrid,*lUrbGrid,mass,M5,b,c); + + Compare("Mobius M ",DaggerNo ,Dv,src,Dl,srclex); + Compare("Mobius Mdag",DaggerYes,Dv,src,Dl,srclex); + + // MooeeInv is the Ls-direction solve; check it inverts Mooee in the lex chart + typename lexMob::FermionField lo(lFrbGrid),t1(lFrbGrid),t2(lFrbGrid),e(lFrbGrid); + lo.Checkerboard()=Odd; t1.Checkerboard()=Odd; t2.Checkerboard()=Odd; e.Checkerboard()=Odd; + pickCheckerboard(Odd,lo,srclex); + Dl.Mooee(lo,t1); + Dl.MooeeInv(t1,t2); + e = t2 - lo; + RealD n = norm2(e)/norm2(lo); + std::cout << GridLogMessage << "Mobius lex MooeeInv(Mooee) - 1: relative " << n << std::endl; + GRID_ASSERT( n < tol ); + } + + std::cout << GridLogMessage << "Test_dwf_lex: ALL PASS" << std::endl; + + Grid_finalize(); +} diff --git a/tests/lexLattice/Test_innerproduct_norm_lex.cc b/tests/lexLattice/Test_innerproduct_norm_lex.cc new file mode 100644 index 000000000..9530d1cbb --- /dev/null +++ b/tests/lexLattice/Test_innerproduct_norm_lex.cc @@ -0,0 +1,132 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/lexLattice/Test_innerproduct_norm_lex.cc + + Copyright (C) 2026 + +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 */ + +// +// Lexicographic (Nsimd()==1) lattice types: linear combination, expression +// templates, norm2, innerProduct, innerProductNorm and RNG fill. +// +#include + +using namespace Grid; + +const RealD tol = 1.0e-5; + +template +void BasicChecks(GridCartesian *grid,GridParallelRNG &pRNG,std::string precision) +{ + RealD vol = grid->gSites(); + + Field a(grid); a = 1.0; + Field b(grid); b = 2.0; + Field c(grid); c = 3.0; + Field d(grid); + + /////////////////////////////////////////////////// + // Linear combination through expression templates + /////////////////////////////////////////////////// + d = a + b - c; + RealD n_zero = norm2(d); + + RealD n_a = norm2(a); + RealD n_b = norm2(b); + + std::cout << GridLogMessage << precision << ": norm2(1+2-3) = " << n_zero << std::endl; + std::cout << GridLogMessage << precision << ": norm2(1) = " << n_a << " expect " << vol << std::endl; + std::cout << GridLogMessage << precision << ": norm2(2) = " << n_b << " expect " << 4.0*vol << std::endl; + + GRID_ASSERT( n_zero < tol ); + GRID_ASSERT( fabs(n_a - vol) < tol*vol ); + GRID_ASSERT( fabs(n_b - 4.0*vol) < tol*vol ); + + /////////////////////////////////////////////////// + // innerProduct of constant fields + /////////////////////////////////////////////////// + ComplexD ip = innerProduct(a,b); + std::cout << GridLogMessage << precision << ": innerProduct(1,2) = " << ip + << " expect (" << 2.0*vol << ",0)" << std::endl; + GRID_ASSERT( fabs(real(ip) - 2.0*vol) < tol*vol ); + GRID_ASSERT( fabs(imag(ip)) < tol*vol ); + + /////////////////////////////////////////////////// + // RNG fill; norm2(r+r) == 4 norm2(r) + /////////////////////////////////////////////////// + Field r(grid); random(pRNG,r); + Field s(grid); random(pRNG,s); + + RealD n_r = norm2(r); + Field rr(grid); rr = r + r; + RealD n_rr = norm2(rr); + + std::cout << GridLogMessage << precision << ": norm2(r) = " << n_r + << " norm2(r+r) = " << n_rr << " ratio " << n_rr/n_r << std::endl; + GRID_ASSERT( n_r > tol ); // RNG actually filled it + GRID_ASSERT( fabs(n_rr - 4.0*n_r) < tol*n_rr ); + + /////////////////////////////////////////////////// + // Fused innerProductNorm against separate calls + /////////////////////////////////////////////////// + ComplexD ip_ref = innerProduct(r,s); + RealD n2_ref = norm2(r); + + ComplexD ip_res; + RealD n2_res; + innerProductNorm(ip_res,n2_res,r,s); + + std::cout << GridLogMessage << precision << ": innerProductNorm ip diff " + << abs(ip_ref-ip_res) << " norm2 diff " << fabs(n2_ref-n2_res) << std::endl; + GRID_ASSERT( abs(ip_ref-ip_res) < tol*abs(ip_ref) + tol ); + GRID_ASSERT( fabs(n2_ref-n2_res) < tol*n2_ref ); + + std::cout << GridLogMessage << precision << ": all checks passed" << std::endl; +} + +int main(int argc, char** argv) +{ + Grid_init(&argc, &argv); + + Coordinate latt = GridDefaultLatt(); + Coordinate simd({1,1,1,1}); // lexicographic: no SIMD + Coordinate mpi = GridDefaultMpi(); + + GridCartesian grid(latt,simd,mpi); + + std::cout << GridLogMessage << "Lexicographic grid, Nsimd = " << grid.Nsimd() << std::endl; + GRID_ASSERT( grid.Nsimd() == 1 ); + + GridParallelRNG pRNG(&grid); + pRNG.SeedFixedIntegers(std::vector({1,2,3,4})); + + + + BasicChecks(&grid,pRNG,"Double"); + BasicChecks(&grid,pRNG,"Single"); + + std::cout << GridLogMessage << "Test_innerproduct_norm_lex: ALL PASS" << std::endl; + + Grid_finalize(); +} diff --git a/tests/lexLattice/Test_io.cc b/tests/lexLattice/Test_io.cc new file mode 100644 index 000000000..fcb9a4576 --- /dev/null +++ b/tests/lexLattice/Test_io.cc @@ -0,0 +1,194 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/lexLattice/Test_io.cc + + Copyright (C) 2026 + +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 */ + +// +// Interchange between vectorised and lexicographic (Nsimd()==1) lattices. +// +// - transfer() in both directions via an unvectorise/vectorise pair +// - RNG seeded identically on both layouts produces identical fields +// - binary I/O written from one layout and read into the other +// +#include + +using namespace Grid; + +const RealD tol = 1.0e-10; + +//////////////////////////////////////////////////////////////////////// +// Layout interchange. Both lattices share the same scalar_object, so a +// single lexicographic array mediates; works in either direction. +//////////////////////////////////////////////////////////////////////// +template +void transfer(Lattice &out,const Lattice &in) +{ + typedef typename vobjIn::scalar_object sobj; + static_assert(std::is_same::value, + "transfer: lattices must share scalar_object"); + + GRID_ASSERT(out.Grid()->gSites() == in.Grid()->gSites()); + + std::vector buf; + unvectorizeToLexOrdArray(buf,in); + vectorizeFromLexOrdArray(buf,out); +} + +//////////////////////////////////////////////////////////////////////// +// Binary write/read of a single lattice object +//////////////////////////////////////////////////////////////////////// +template +void writeField(const Field &f,std::string file) +{ + typedef typename Field::vector_object vobj; + typedef typename Field::scalar_object sobj; + + BinarySimpleMunger munge; + std::string format = getFormatString(); + uint64_t offset = 0; + uint32_t nersc_csum,scidac_csuma,scidac_csumb; + + BinaryIO::writeLatticeObject(const_cast(f),file,munge,offset,format, + nersc_csum,scidac_csuma,scidac_csumb); +} + +template +void readField(Field &f,std::string file) +{ + typedef typename Field::vector_object vobj; + typedef typename Field::scalar_object sobj; + + BinarySimpleMunger munge; + std::string format = getFormatString(); + uint64_t offset = 0; + uint32_t nersc_csum,scidac_csuma,scidac_csumb; + + BinaryIO::readLatticeObject(f,file,munge,offset,format, + nersc_csum,scidac_csuma,scidac_csumb); +} + +int main(int argc, char** argv) +{ + Grid_init(&argc, &argv); + + Coordinate latt = GridDefaultLatt(); + Coordinate mpi = GridDefaultMpi(); + Coordinate vsimd = GridDefaultSimd(Nd,vComplexD::Nsimd()); + Coordinate ssimd({1,1,1,1}); + + GridCartesian vGrid(latt,vsimd,mpi); + GridCartesian sGrid(latt,ssimd,mpi); + + std::cout << GridLogMessage << "vectorised grid Nsimd = " << vGrid.Nsimd() << std::endl; + std::cout << GridLogMessage << "lexicographic grid Nsimd = " << sGrid.Nsimd() << std::endl; + GRID_ASSERT( sGrid.Nsimd() == 1 ); + + typedef LatticeComplexD vField; + typedef lexLatticeComplexD sField; + + std::vector seeds({1,2,3,4}); + + GridParallelRNG vRNG(&vGrid); vRNG.SeedFixedIntegers(seeds); + GridParallelRNG sRNG(&sGrid); sRNG.SeedFixedIntegers(seeds); + + vField v(&vGrid); random(vRNG,v); + sField s(&sGrid); random(sRNG,s); + + RealD nv = norm2(v); + RealD ns = norm2(s); + std::cout << GridLogMessage << "norm2 vectorised = " << nv << std::endl; + std::cout << GridLogMessage << "norm2 lexicographic = " << ns << std::endl; + GRID_ASSERT( nv > tol ); + + //////////////////////////////////////////////////////// + // A: transfer round trip v -> s -> v + //////////////////////////////////////////////////////// + { + sField st(&sGrid); + vField vt(&vGrid); + + transfer(st,v); + transfer(vt,st); + + vField d(&vGrid); d = vt - v; + RealD n = norm2(d); + std::cout << GridLogMessage << "A: round trip v->s->v norm2(diff) = " << n << std::endl; + GRID_ASSERT( n < tol ); + + // and the transferred copy must carry the same norm + std::cout << GridLogMessage << "A: norm2 transferred = " << norm2(st) << std::endl; + GRID_ASSERT( fabs(norm2(st) - nv) < tol*nv ); + } + + //////////////////////////////////////////////////////// + // B: identically seeded RNGs agree across layouts + //////////////////////////////////////////////////////// + { + vField vs(&vGrid); + transfer(vs,s); + + vField d(&vGrid); d = vs - v; + RealD n = norm2(d); + std::cout << GridLogMessage << "B: same seed, both layouts norm2(diff) = " << n << std::endl; + GRID_ASSERT( n < tol ); + } + + //////////////////////////////////////////////////////// + // C: write vectorised, read lexicographic + //////////////////////////////////////////////////////// + { + sField sref(&sGrid); transfer(sref,v); // what the file should contain + sField sin(&sGrid); sin = Zero(); + + writeField(v,"nonsimd_io_v.bin"); + readField(sin,"nonsimd_io_v.bin"); + + sField d(&sGrid); d = sin - sref; + RealD n = norm2(d); + std::cout << GridLogMessage << "C: write simd / read lexicographic norm2(diff) = " << n << std::endl; + GRID_ASSERT( n < tol ); + } + + //////////////////////////////////////////////////////// + // D: write lexicographic, read vectorised + //////////////////////////////////////////////////////// + { + sField sout(&sGrid); transfer(sout,v); + vField vin(&vGrid); vin = Zero(); + + writeField(sout,"nonsimd_io_s.bin"); + readField(vin,"nonsimd_io_s.bin"); + + vField d(&vGrid); d = vin - v; + RealD n = norm2(d); + std::cout << GridLogMessage << "D: write lexicographic / read simd norm2(diff) = " << n << std::endl; + GRID_ASSERT( n < tol ); + } + + std::cout << GridLogMessage << "Test_io: ALL PASS" << std::endl; + + Grid_finalize(); +} diff --git a/tests/lexLattice/Test_staggered_lex.cc b/tests/lexLattice/Test_staggered_lex.cc new file mode 100644 index 000000000..f50c84258 --- /dev/null +++ b/tests/lexLattice/Test_staggered_lex.cc @@ -0,0 +1,246 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/lexLattice/Test_staggered_lex.cc + + Copyright (C) 2026 + +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 */ + +// +// Naive staggered fermion operator on a lexicographic (Nsimd()==1) lattice. +// +// Two independent checks: +// - lexicographic Dhop against a covariant-shift reference built in the +// same chart +// - both charts driven from the same gauge field and source, compared +// elementwise through the layout interchange +// +#include + +using namespace Grid; + +const RealD tol = 1.0e-10; + +typedef NaiveStaggeredFermion vStagOp; +typedef NaiveStaggeredFermion lexStagOp; + +//////////////////////////////////////////////////////////////////////// +// Layout interchange; both lattices share the same scalar_object. +//////////////////////////////////////////////////////////////////////// +template +void transfer(Lattice &out,const Lattice &in) +{ + typedef typename vobjIn::scalar_object sobj; + static_assert(std::is_same::value, + "transfer: lattices must share scalar_object"); + + GRID_ASSERT(out.Grid()->gSites() == in.Grid()->gSites()); + + std::vector buf; + unvectorizeToLexOrdArray(buf,in); + vectorizeFromLexOrdArray(buf,out); +} + +//////////////////////////////////////////////////////////////////////// +// Dhop from covariant shifts, in whichever chart Gimpl names +//////////////////////////////////////////////////////////////////////// +template +void ReferenceDhop(typename Impl::FermionField &ref, + const typename Impl::GaugeField &Umu, + const typename Impl::FermionField &src, + RealD c1,RealD u0) +{ + typedef typename Impl::GaugeLinkField LinkField; + typedef typename Impl::ComplexField ComplexField; + typedef typename Impl::FermionField FermionField; + + GridBase *grid = src.Grid(); + RealD c1tad = 0.5*c1/u0; + + std::vector U(Nd,grid); + for(int mu=0;mu(Umu,mu); + } + + typedef Lattice::Integerified> > IntField; + + IntField x(grid); LatticeCoordinate(x,0); + IntField y(grid); LatticeCoordinate(y,1); + IntField z(grid); LatticeCoordinate(z,2); + IntField lin_z(grid); lin_z = x+y; + IntField lin_t(grid); lin_t = x+y+z; + + FermionField tmp(grid); + ref = Zero(); + + for(int mu=0;mu = " << lhs + << " = " << rhs << " relative " << n << std::endl; + GRID_ASSERT( n < tol ); + } + + ////////////////////////////////////////////////// + // Checkerboarded hopping term, both charts + ////////////////////////////////////////////////// + { + typename vStagOp::FermionField vsrc_o(&vRBGrid); vsrc_o.Checkerboard()=Odd; + typename vStagOp::FermionField vres_e(&vRBGrid); vres_e.Checkerboard()=Even; + typename lexStagOp::FermionField lsrc_o(&lRBGrid); lsrc_o.Checkerboard()=Odd; + typename lexStagOp::FermionField lres_e(&lRBGrid); lres_e.Checkerboard()=Even; + + pickCheckerboard(Odd,vsrc_o,src); + pickCheckerboard(Odd,lsrc_o,srclex); + + Ds.Meooe (vsrc_o,vres_e); + Dslex.Meooe(lsrc_o,lres_e); + + RealD vn = norm2(vres_e); + RealD ln = norm2(lres_e); + RealD n = fabs(vn-ln)/vn; + std::cout << GridLogMessage << "Meooe |simd|^2 " << vn << " |lex|^2 " << ln + << " relative " << n << std::endl; + GRID_ASSERT( n < tol ); + } + + std::cout << GridLogMessage << "Test_staggered_lex: ALL PASS" << std::endl; + + Grid_finalize(); +} diff --git a/tests/lexLattice/Test_stencil_lex.cc b/tests/lexLattice/Test_stencil_lex.cc new file mode 100644 index 000000000..f8ec1dcfd --- /dev/null +++ b/tests/lexLattice/Test_stencil_lex.cc @@ -0,0 +1,122 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/lexLattice/Test_stencil_lex.cc + + Copyright (C) 2026 + +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 */ + +// +// CartesianStencil on a lexicographic (Nsimd()==1) lattice, checked against +// Cshift for every direction and displacement of either sign. +// +// On a unit simd grid no stencil entry may request a permute; the Grid_simd1 +// permute is an assert, so a request would abort rather than pass silently. +// +#include + +using namespace Grid; + +const RealD tol = 1.0e-10; + +int main(int argc, char ** argv) +{ + Grid_init(&argc, &argv); + + typedef lexLatticeComplexD Field; + typedef Field::vector_object vobj; + + Coordinate latt = GridDefaultLatt(); + Coordinate simd({1,1,1,1}); + Coordinate mpi = GridDefaultMpi(); + + GridCartesian Fine(latt,simd,mpi); + + std::cout << GridLogMessage << "Lexicographic grid, Nsimd = " << Fine.Nsimd() << std::endl; + GRID_ASSERT( Fine.Nsimd() == 1 ); + + GridParallelRNG fRNG(&Fine); + fRNG.SeedFixedIntegers(std::vector({1,2,3,4})); + + Field Foo(&Fine); random(fRNG,Foo); + Field Bar(&Fine); + Field Check(&Fine); + Field Diff(&Fine); + + typedef CartesianStencil Stencil; + SimpleStencilParams p; + + for(int dir=0;dir directions(npoint,dir); + std::vector displacements(npoint,disp); + + Stencil myStencil(&Fine,npoint,0,directions,displacements,p); + + SimpleCompressor compress; + myStencil.HaloExchange(Foo,compress); + + Bar = Cshift(Foo,dir,disp); + + { + autoView( check , Check, AcceleratorWrite); + autoView( foo , Foo, AcceleratorRead); + autoView( st_v , myStencil, AcceleratorRead); + auto CBp=myStencil.CommBuf(); + accelerator_for(i,Check.Grid()->oSites(), 1, { + + int permute_type; + StencilEntry *SE; + SE = st_v.GetEntry(permute_type,0,i); + + if ( SE->_is_local && SE->_permute ) + permute(check[i],foo[SE->_offset],permute_type); + else if (SE->_is_local) + check[i] = foo[SE->_offset]; + else + check[i] = CBp[SE->_offset]; + }); + } + + Diff = Check-Bar; + RealD nrm = norm2(Diff); + + if ( nrm > tol ) { + std::cout << GridLogMessage << "FAIL dir " << dir << " disp " << disp + << " norm2(stencil-cshift) = " << nrm << std::endl; + } + GRID_ASSERT( nrm < tol ); + } + std::cout << GridLogMessage << "Stencil == Cshift, dir " << dir + << ", all displacements " << -(L-1) << " .. " << L-1 << std::endl; + } + + std::cout << GridLogMessage << "Test_stencil_lex: ALL PASS" << std::endl; + + Grid_finalize(); +} diff --git a/tests/lexLattice/Test_wilson_lex.cc b/tests/lexLattice/Test_wilson_lex.cc new file mode 100644 index 000000000..4f2240c9c --- /dev/null +++ b/tests/lexLattice/Test_wilson_lex.cc @@ -0,0 +1,239 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/lexLattice/Test_wilson_lex.cc + + Copyright (C) 2026 + +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 */ + +// +// Wilson fermion operator on a lexicographic (Nsimd()==1) lattice. +// +// Exercises the spin-projected half spinor path -- WilsonCompressor and +// WilsonStencil -- which the staggered and stencil tests do not reach. +// +// - lexicographic Dhop against a covariant-shift reference in the same chart +// - both charts from one gauge field and source, compared elementwise +// - gamma5 hermiticity and the even-odd hopping term in the lex chart +// +#include + +using namespace Grid; + +const RealD tol = 1.0e-10; + +Gamma::Algebra Gmu [] = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT +}; + +typedef WilsonFermion vWilsonOp; +typedef WilsonFermion lexWilsonOp; + +//////////////////////////////////////////////////////////////////////// +// Layout interchange; both lattices share the same scalar_object. +//////////////////////////////////////////////////////////////////////// +template +void transfer(Lattice &out,const Lattice &in) +{ + typedef typename vobjIn::scalar_object sobj; + static_assert(std::is_same::value, + "transfer: lattices must share scalar_object"); + + GRID_ASSERT(out.Grid()->gSites() == in.Grid()->gSites()); + + std::vector buf; + unvectorizeToLexOrdArray(buf,in); + vectorizeFromLexOrdArray(buf,out); +} + +//////////////////////////////////////////////////////////////////////// +// Wilson hopping term from covariant shifts, in whichever chart Impl names +//////////////////////////////////////////////////////////////////////// +template +void ReferenceDhop(typename Impl::FermionField &ref, + const typename Impl::GaugeField &Umu, + const typename Impl::FermionField &src) +{ + typedef typename Impl::GaugeLinkField LinkField; + typedef typename Impl::FermionField FermionField; + + GridBase *grid = src.Grid(); + + std::vector U(Nd,grid); + for(int mu=0;mu(Umu,mu); + } + + FermionField tmp(grid); + ref = Zero(); + + for(int mu=0;mu