From 1c4bc7ed38c6e008972e205920d11e31880313d0 Mon Sep 17 00:00:00 2001 From: paboyle Date: Fri, 31 Mar 2017 14:41:48 +0900 Subject: [PATCH] Debugged staggered conventions --- .../fermion/ImprovedStaggeredFermion.cc | 61 ++- .../action/fermion/ImprovedStaggeredFermion.h | 9 + lib/qcd/action/fermion/StaggeredKernelsAsm.cc | 13 +- .../action/fermion/StaggeredKernelsHand.cc | 21 +- scripts/zmobius.sh | 35 ++ tests/qdpxx/Test_qdpxx_stag.cc | 364 ++++++++++++++++++ 6 files changed, 491 insertions(+), 12 deletions(-) create mode 100644 scripts/zmobius.sh create mode 100644 tests/qdpxx/Test_qdpxx_stag.cc diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc b/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc index ec5811e0..2ba4f4af 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc @@ -40,10 +40,10 @@ ImprovedStaggeredFermionStatic::displacements({1, 1, 1, 1, -1, -1, -1, -1, 3, 3, // Constructor and gauge import ///////////////////////////////// + template -ImprovedStaggeredFermion::ImprovedStaggeredFermion(GaugeField &_Uthin, GaugeField &_Ufat, GridCartesian &Fgrid, - GridRedBlackCartesian &Hgrid, RealD _mass, - RealD _c1, RealD _c2,RealD _u0, +ImprovedStaggeredFermion::ImprovedStaggeredFermion(GridCartesian &Fgrid, GridRedBlackCartesian &Hgrid, + RealD _mass, const ImplParams &p) : Kernels(p), _grid(&Fgrid), @@ -52,9 +52,6 @@ ImprovedStaggeredFermion::ImprovedStaggeredFermion(GaugeField &_Uthin, Gau 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), @@ -65,9 +62,29 @@ ImprovedStaggeredFermion::ImprovedStaggeredFermion(GaugeField &_Uthin, Gau UUUmuOdd(&Hgrid) , _tmp(&Hgrid) { - // Allocate the required comms buffer +} + +template +ImprovedStaggeredFermion::ImprovedStaggeredFermion(GaugeField &_Uthin, GaugeField &_Ufat, GridCartesian &Fgrid, + GridRedBlackCartesian &Hgrid, RealD _mass, + RealD _c1, RealD _c2,RealD _u0, + const ImplParams &p) + : ImprovedStaggeredFermion(Fgrid,Hgrid,_mass,p) +{ + c1=_c1; + c2=_c2; + u0=_u0; ImportGauge(_Uthin,_Ufat); } +template +ImprovedStaggeredFermion::ImprovedStaggeredFermion(GaugeField &_Uthin,GaugeField &_Utriple, GaugeField &_Ufat, GridCartesian &Fgrid, + GridRedBlackCartesian &Hgrid, RealD _mass, + const ImplParams &p) + : ImprovedStaggeredFermion(Fgrid,Hgrid,_mass,p) +{ + ImportGaugeSimple(_Utriple,_Ufat); +} + //////////////////////////////////////////////////////////// // Momentum space propagator should be @@ -86,6 +103,34 @@ void ImprovedStaggeredFermion::ImportGauge(const GaugeField &_Uthin) ImportGauge(_Uthin,_Uthin); }; template +void ImprovedStaggeredFermion::ImportGaugeSimple(const GaugeField &_Utriple,const GaugeField &_Ufat) +{ + ///////////////////////////////////////////////////////////////// + // Trivial import; phases and fattening and such like preapplied + ///////////////////////////////////////////////////////////////// + GaugeLinkField U(GaugeGrid()); + + for (int mu = 0; mu < Nd; mu++) { + + U = PeekIndex(_Utriple, mu); + PokeIndex(UUUmu, U, mu ); + + U = adj( Cshift(U, mu, -3)); + PokeIndex(UUUmu, -U, mu+4 ); + + U = PeekIndex(_Ufat, mu); + PokeIndex(Umu, U, mu); + + U = adj( Cshift(U, mu, -1)); + PokeIndex(Umu, -U, mu+4); + + } + pickCheckerboard(Even, UmuEven, Umu); + pickCheckerboard(Odd, UmuOdd , Umu); + pickCheckerboard(Even, UUUmuEven,UUUmu); + pickCheckerboard(Odd, UUUmuOdd, UUUmu); +} +template void ImprovedStaggeredFermion::ImportGauge(const GaugeField &_Uthin,const GaugeField &_Ufat) { GaugeLinkField U(GaugeGrid()); @@ -115,6 +160,8 @@ void ImprovedStaggeredFermion::ImportGauge(const GaugeField &_Uthin,const PokeIndex(UUUmu, U*(-0.5*c2/u0/u0/u0), mu+4); } + std::cout << " Umu " << Umu._odata[0]<, public ImprovedS RealD _c1=9.0/8.0, RealD _c2=-1.0/24.0,RealD _u0=1.0, const ImplParams &p = ImplParams()); + ImprovedStaggeredFermion(GaugeField &_Uthin, GaugeField &_Utriple, GaugeField &_Ufat, GridCartesian &Fgrid, + GridRedBlackCartesian &Hgrid, RealD _mass, + const ImplParams &p = ImplParams()); + + ImprovedStaggeredFermion(GridCartesian &Fgrid, GridRedBlackCartesian &Hgrid, RealD _mass, + const ImplParams &p = ImplParams()); + + // DoubleStore impl dependent + void ImportGaugeSimple(const GaugeField &_Utriple, const GaugeField &_Ufat); void ImportGauge(const GaugeField &_Uthin, const GaugeField &_Ufat); void ImportGauge(const GaugeField &_Uthin); diff --git a/lib/qcd/action/fermion/StaggeredKernelsAsm.cc b/lib/qcd/action/fermion/StaggeredKernelsAsm.cc index bfe13f07..fd881716 100644 --- a/lib/qcd/action/fermion/StaggeredKernelsAsm.cc +++ b/lib/qcd/action/fermion/StaggeredKernelsAsm.cc @@ -587,7 +587,6 @@ void StaggeredKernels::DhopSiteAsm(StencilImpl &st, LebesgueOrder &lo, int sU, const FermionField &in, FermionField &out) { assert(0); - }; @@ -905,9 +904,17 @@ template <> void StaggeredKernels::DhopSiteAsm(StencilImpl &st, #endif } +#define KERNEL_INSTANTIATE(CLASS,FUNC,IMPL) \ + template void CLASS::FUNC(StencilImpl &st, LebesgueOrder &lo, \ + DoubledGaugeField &U, \ + DoubledGaugeField &UUU, \ + SiteSpinor *buf, int LLs, \ + int sU, const FermionField &in, FermionField &out); -FermOpStaggeredTemplateInstantiate(StaggeredKernels); -FermOpStaggeredVec5dTemplateInstantiate(StaggeredKernels); +KERNEL_INSTANTIATE(StaggeredKernels,DhopSiteAsm,StaggeredImplD); +KERNEL_INSTANTIATE(StaggeredKernels,DhopSiteAsm,StaggeredImplF); +KERNEL_INSTANTIATE(StaggeredKernels,DhopSiteAsm,StaggeredVec5dImplD); +KERNEL_INSTANTIATE(StaggeredKernels,DhopSiteAsm,StaggeredVec5dImplF); }} diff --git a/lib/qcd/action/fermion/StaggeredKernelsHand.cc b/lib/qcd/action/fermion/StaggeredKernelsHand.cc index f5806657..7de8480c 100644 --- a/lib/qcd/action/fermion/StaggeredKernelsHand.cc +++ b/lib/qcd/action/fermion/StaggeredKernelsHand.cc @@ -299,7 +299,24 @@ void StaggeredKernels::DhopSiteDepthHand(StencilImpl &st, LebesgueOrder &l } -FermOpStaggeredTemplateInstantiate(StaggeredKernels); -FermOpStaggeredVec5dTemplateInstantiate(StaggeredKernels); +#define DHOP_SITE_HAND_INSTANTIATE(IMPL) \ + template void StaggeredKernels::DhopSiteHand(StencilImpl &st, LebesgueOrder &lo, \ + DoubledGaugeField &U,DoubledGaugeField &UUU, \ + SiteSpinor *buf, int LLs, \ + int sU, const FermionField &in, FermionField &out, int dag); + +#define DHOP_SITE_DEPTH_HAND_INSTANTIATE(IMPL) \ + template void StaggeredKernels::DhopSiteDepthHand(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, \ + SiteSpinor *buf, int sF, \ + int sU, const FermionField &in, SiteSpinor &out,int threeLink) ; +DHOP_SITE_HAND_INSTANTIATE(StaggeredImplD); +DHOP_SITE_HAND_INSTANTIATE(StaggeredImplF); +DHOP_SITE_HAND_INSTANTIATE(StaggeredVec5dImplD); +DHOP_SITE_HAND_INSTANTIATE(StaggeredVec5dImplF); + +DHOP_SITE_DEPTH_HAND_INSTANTIATE(StaggeredImplD); +DHOP_SITE_DEPTH_HAND_INSTANTIATE(StaggeredImplF); +DHOP_SITE_DEPTH_HAND_INSTANTIATE(StaggeredVec5dImplD); +DHOP_SITE_DEPTH_HAND_INSTANTIATE(StaggeredVec5dImplF); }} diff --git a/scripts/zmobius.sh b/scripts/zmobius.sh new file mode 100644 index 00000000..04b223d2 --- /dev/null +++ b/scripts/zmobius.sh @@ -0,0 +1,35 @@ +#!/bin/bash +fn=$1 + +grep "double zmobius_" $fn | +awk 'BEGIN{ m["zmobius_b_coeff"]=0; m["zmobius_c_coeff"]=1; }{ val[m[substr($2,0,15)]][substr($2,17)+0]=$4; }END{ + + ls=length(val[0])/2; + + print "ls = " ls + + bmc=-111; + + for (s=0;s(%.15g,%.15g) );\n",omegar[s],omegai[s]); + } + +}' diff --git a/tests/qdpxx/Test_qdpxx_stag.cc b/tests/qdpxx/Test_qdpxx_stag.cc new file mode 100644 index 00000000..a7563924 --- /dev/null +++ b/tests/qdpxx/Test_qdpxx_stag.cc @@ -0,0 +1,364 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/qdpxx/Test_qdpxx_munprec.cc + + Copyright (C) 2015 + +Author: Azusa Yamaguchi +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 + +double mq=0.1; + +typedef Grid::QCD::StaggeredImplR::FermionField FermionField; +typedef Grid::QCD::LatticeGaugeField GaugeField; + +void make_gauge (GaugeField & lat, FermionField &src); +void calc_grid (GaugeField & lat, GaugeField & uthin,GaugeField & ufat, FermionField &src, FermionField &res,int dag); +void calc_chroma (GaugeField & lat,GaugeField & uthin,GaugeField & ufat, FermionField &src, FermionField &res,int dag); + +#include +#include +#include + +namespace Chroma { + + +class ChromaWrapper { +public: + + typedef multi1d U; + typedef LatticeStaggeredFermion T4; + + static void ImportGauge(GaugeField & gr, + QDP::multi1d & ch) + { + Grid::QCD::LorentzColourMatrix LCM; + Grid::Complex cc; + QDP::ColorMatrix cm; + QDP::Complex c; + + std::vector x(4); + QDP::multi1d cx(4); + std::vector gd= gr._grid->GlobalDimensions(); + + for (x[0]=0;x[0] & ch) + { + Grid::QCD::LorentzColourMatrix LCM; + Grid::Complex cc; + QDP::ColorMatrix cm; + QDP::Complex c; + + std::vector x(4); + QDP::multi1d cx(4); + std::vector gd= gr._grid->GlobalDimensions(); + + for (x[0]=0;x[0] x(5); + QDP::multi1d cx(4); + std::vector gd= gr._grid->GlobalDimensions(); + + for (x[0]=0;x[0] x(5); + QDP::multi1d cx(4); + std::vector gd= gr._grid->GlobalDimensions(); + + for (x[0]=0;x[0] > GetLinOp (U &u,U &u_fat,U &u_triple) + { + QDP::Real _mq(mq); + QDP::multi1d bcs(QDP::Nd); + + bcs[0] = bcs[1] = bcs[2] = bcs[3] = 1; + + Chroma::AsqtadFermActParams p; + p.Mass = _mq; + p.u0 = Real(1.0); + + + Chroma::Handle > fbc(new Chroma::SimpleFermBC< T4, U, U >(bcs)); + Chroma::Handle > cfs( new Chroma::CreateSimpleFermState(fbc)); + Chroma::AsqtadFermAct S_f(cfs,p); + Chroma::Handle< Chroma::FermState > ffs( S_f.createState(u) ); + u_fat =ffs.cast()->getFatLinks(); + u_triple=ffs.cast()->getTripleLinks(); + return S_f.linOp(ffs); + } + +}; +} + +int main (int argc,char **argv ) +{ + + /******************************************************** + * Setup QDP + *********************************************************/ + Chroma::initialize(&argc,&argv); + Chroma::WilsonTypeFermActs4DEnv::registerAll(); + + /******************************************************** + * Setup Grid + *********************************************************/ + Grid::Grid_init(&argc,&argv); + Grid::GridCartesian * UGrid = Grid::QCD::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(), + Grid::GridDefaultSimd(Grid::QCD::Nd,Grid::vComplex::Nsimd()), + Grid::GridDefaultMpi()); + + std::vector gd = UGrid->GlobalDimensions(); + QDP::multi1d nrow(QDP::Nd); + for(int mu=0;mu<4;mu++) nrow[mu] = gd[mu]; + + QDP::Layout::setLattSize(nrow); + QDP::Layout::create(); + + GaugeField uthin (UGrid); + GaugeField ufat (UGrid); + GaugeField utriple(UGrid); + FermionField src(UGrid); + FermionField res_chroma(UGrid); + FermionField res_grid (UGrid); + + + { + + std::cout << "*****************************"< U; + + U u(4); + U ut(4); + U uf(4); + + // Chroma::HotSt(u); + Chroma::ChromaWrapper::ImportGauge(lat,u) ; + + QDP::LatticeStaggeredFermion check; + QDP::LatticeStaggeredFermion result; + QDP::LatticeStaggeredFermion tmp; + QDP::LatticeStaggeredFermion psi; + + Chroma::ChromaWrapper::ImportFermion(src,psi); + + auto linop =Chroma::ChromaWrapper::GetLinOp(u,uf,ut); + + Chroma::ChromaWrapper::ExportGauge(uthin,ut) ; + Chroma::ChromaWrapper::ExportGauge(ufat ,uf) ; + + enum Chroma::PlusMinus isign; + if ( dag ) { + isign=Chroma::MINUS; + } else { + isign=Chroma::PLUS; + } + + std::cout << "Calling Chroma Linop "<< std::endl; + linop->evenEvenLinOp(tmp,psi,isign); check[rb[0]] = tmp; + linop->oddOddLinOp (tmp,psi,isign); check[rb[1]] = tmp; + linop->evenOddLinOp(tmp,psi,isign) ; check[rb[0]]+= tmp; + linop->oddEvenLinOp(tmp,psi,isign) ; check[rb[1]]+= tmp; + + Chroma::ChromaWrapper::ExportFermion(res,check) ; +} + + +void make_gauge(GaugeField & Umu,FermionField &src) +{ + using namespace Grid; + using namespace Grid::QCD; + + std::vector seeds4({1,2,3,4}); + + Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu._grid; + Grid::GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + Grid::QCD::SU3::HotConfiguration(RNG4,Umu); + Grid::gaussian(RNG4,src); +} + +void calc_grid(GaugeField & Uthin, GaugeField & Utriple, GaugeField & Ufat, FermionField &src, FermionField &res,int dag) +{ + using namespace Grid; + using namespace Grid::QCD; + + Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Uthin._grid; + Grid::GridRedBlackCartesian * UrbGrid = Grid::QCD::SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + + Grid::QCD::ImprovedStaggeredFermionR Dstag(Uthin,Utriple,Ufat,*UGrid,*UrbGrid,mq*2.0); + + std::cout << Grid::GridLogMessage <<" Calling Grid staggered multiply "<