diff --git a/benchmarks/Benchmark_dwf.cc b/benchmarks/Benchmark_dwf.cc index 0d3e2653..bc9ab708 100644 --- a/benchmarks/Benchmark_dwf.cc +++ b/benchmarks/Benchmark_dwf.cc @@ -86,18 +86,6 @@ int main (int argc, char ** argv) LatticeFermion tmp(FGrid); LatticeFermion err(FGrid); - /* src=zero; - std::vector origin(5,0); - SpinColourVector f=zero; - for(int sp=0;sp<4;sp++){ - for(int co=0;co<3;co++){ - f()(sp)(co)=Complex(1.0,0.0); - }} - pokeSite(f,src,origin); - */ - - ColourMatrix cm = Complex(1.0,0.0); - LatticeGaugeField Umu(UGrid); random(RNG4,Umu); @@ -144,7 +132,8 @@ int main (int argc, char ** argv) DomainWallFermionR Dw(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); - std::cout< 1.0e-6 ) { - std::cout << "site "<SendToRecvFrom(xmit,to,rcv,from,bytes); @@ -304,23 +303,23 @@ double calls; void ZeroCounters(void) { - gathertime=0; - jointime=0; - commtime=0; - halogtime=0; - mergetime=0; - spintime=0; - gathermtime=0; - splicetime=0; - nosplicetime=0; - comms_bytes=0; - calls=0; + gathertime = 0.; + jointime = 0.; + commtime = 0.; + halogtime = 0.; + mergetime = 0.; + spintime = 0.; + gathermtime = 0.; + splicetime = 0.; + nosplicetime = 0.; + comms_bytes = 0.; + calls = 0.; }; void Report(void) { #define PRINTIT(A) \ std::cout << GridLogMessage << " Stencil " << #A << " "<< A/calls< 0 ) { + if ( calls > 0. ) { std::cout << GridLogMessage << " Stencil calls "< 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 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. +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. +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 */ +See the full license in the file "LICENSE" in the top level distribution +directory +*************************************************************************************/ +/* END LEGAL */ #ifndef GRID_CONJUGATE_GRADIENT_H #define GRID_CONJUGATE_GRADIENT_H namespace Grid { - ///////////////////////////////////////////////////////////// - // Base classes for iterative processes based on operators - // single input vec, single output vec. - ///////////////////////////////////////////////////////////// +///////////////////////////////////////////////////////////// +// Base classes for iterative processes based on operators +// single input vec, single output vec. +///////////////////////////////////////////////////////////// - template - class ConjugateGradient : public OperatorFunction { -public: - bool ErrorOnNoConverge; //throw an assert when the CG fails to converge. Defaults true. - RealD Tolerance; - Integer MaxIterations; - ConjugateGradient(RealD tol,Integer maxit, bool err_on_no_conv = true) : Tolerance(tol), MaxIterations(maxit), ErrorOnNoConverge(err_on_no_conv){ - }; +template +class ConjugateGradient : public OperatorFunction { + public: + bool ErrorOnNoConverge; // throw an assert when the CG fails to converge. + // Defaults true. + RealD Tolerance; + Integer MaxIterations; + ConjugateGradient(RealD tol, Integer maxit, bool err_on_no_conv = true) + : Tolerance(tol), + MaxIterations(maxit), + ErrorOnNoConverge(err_on_no_conv){}; + void operator()(LinearOperatorBase &Linop, const Field &src, + Field &psi) { + psi.checkerboard = src.checkerboard; + conformable(psi, src); - void operator() (LinearOperatorBase &Linop,const Field &src, Field &psi){ + RealD cp, c, a, d, b, ssq, qq, b_pred; - psi.checkerboard = src.checkerboard; - conformable(psi,src); + Field p(src); + Field mmp(src); + Field r(src); - RealD cp,c,a,d,b,ssq,qq,b_pred; - - Field p(src); - Field mmp(src); - Field r(src); - - //Initial residual computation & set up - RealD guess = norm2(psi); - assert(std::isnan(guess)==0); + // Initial residual computation & set up + RealD guess = norm2(psi); + assert(std::isnan(guess) == 0); - Linop.HermOpAndNorm(psi,mmp,d,b); - - r= src-mmp; - p= r; - - a =norm2(p); - cp =a; - ssq=norm2(src); + + Linop.HermOpAndNorm(psi, mmp, d, b); + - std::cout< WilsonFermion5DStatic::displacements({1,1,1,1,-1,-1,-1,-1 // 5d lattice for DWF. template WilsonFermion5D::WilsonFermion5D(GaugeField &_Umu, - GridCartesian &FiveDimGrid, - GridRedBlackCartesian &FiveDimRedBlackGrid, - GridCartesian &FourDimGrid, - GridRedBlackCartesian &FourDimRedBlackGrid, - RealD _M5,const ImplParams &p) : + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _M5,const ImplParams &p) : Kernels(p), _FiveDimGrid (&FiveDimGrid), _FiveDimRedBlackGrid(&FiveDimRedBlackGrid), @@ -135,10 +135,10 @@ WilsonFermion5D::WilsonFermion5D(GaugeField &_Umu, /* template WilsonFermion5D::WilsonFermion5D(int simd,GaugeField &_Umu, - GridCartesian &FiveDimGrid, - GridRedBlackCartesian &FiveDimRedBlackGrid, - GridCartesian &FourDimGrid, - RealD _M5,const ImplParams &p) : + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + RealD _M5,const ImplParams &p) : { int nsimd = Simd::Nsimd(); @@ -178,26 +178,64 @@ WilsonFermion5D::WilsonFermion5D(int simd,GaugeField &_Umu, template void WilsonFermion5D::Report(void) { - if ( Calls > 0 ) { - std::cout << GridLogMessage << "WilsonFermion5D Dhop Calls " < latt = GridDefaultLatt(); + RealD volume = Ls; for(int mu=0;mu_Nprocessors; - std::cout << GridLogMessage << "WilsonFermion5D Stencil"< 0 ) { + std::cout << GridLogMessage << "#### Dhop calls report " << std::endl; + std::cout << GridLogMessage << "WilsonFermion5D Number of Dhop Calls : " << DhopCalls << std::endl; + std::cout << GridLogMessage << "WilsonFermion5D Total Communication time : " << DhopCommTime + << " us" << std::endl; + std::cout << GridLogMessage << "WilsonFermion5D CommTime/Calls : " + << DhopCommTime / DhopCalls << " us" << std::endl; + std::cout << GridLogMessage << "WilsonFermion5D Total Compute time : " + << DhopComputeTime << " us" << std::endl; + std::cout << GridLogMessage << "WilsonFermion5D ComputeTime/Calls : " + << DhopComputeTime / DhopCalls << " us" << std::endl; - std::cout << GridLogMessage << "WilsonFermion5D StencilEven"< 0 ) { + std::cout << GridLogMessage << "#### Deriv calls report "<< std::endl; + std::cout << GridLogMessage << "WilsonFermion5D Number of Deriv Calls : " < 0 || DhopCalls > 0){ + std::cout << GridLogMessage << "WilsonFermion5D Stencil"< void WilsonFermion5D::ZeroCounters(void) { - Calls=0; - CommTime=0; - ComputeTime=0; + DhopCalls = 0; + DhopCommTime = 0; + DhopComputeTime = 0; + + DerivCalls = 0; + DerivCommTime = 0; + DerivComputeTime = 0; + DerivDhopComputeTime = 0; + Stencil.ZeroCounters(); StencilEven.ZeroCounters(); StencilOdd.ZeroCounters(); @@ -244,12 +282,13 @@ PARALLEL_FOR_LOOP template void WilsonFermion5D::DerivInternal(StencilImpl & st, - DoubledGaugeField & U, - GaugeField &mat, - const FermionField &A, - const FermionField &B, - int dag) + DoubledGaugeField & U, + GaugeField &mat, + const FermionField &A, + const FermionField &B, + int dag) { + DerivCalls++; assert((dag==DaggerNo) ||(dag==DaggerYes)); conformable(st._grid,A._grid); @@ -260,51 +299,53 @@ void WilsonFermion5D::DerivInternal(StencilImpl & st, FermionField Btilde(B._grid); FermionField Atilde(B._grid); + DerivCommTime-=usecond(); st.HaloExchange(B,compressor); + DerivCommTime+=usecond(); Atilde=A; - for(int mu=0;muoSites();sss++){ - for(int s=0;soSites(); sss++) { + for (int s = 0; s < Ls; s++) { + int sU = sss; + int sF = s + Ls * sU; - assert ( sF< B._grid->oSites()); - assert ( sU< U._grid->oSites()); + assert(sF < B._grid->oSites()); + assert(sU < U._grid->oSites()); - Kernels::DiracOptDhopDir(st,U,st.comm_buf,sF,sU,B,Btilde,mu,gamma); - - //////////////////////////// - // spin trace outer product - //////////////////////////// + Kernels::DiracOptDhopDir(st, U, st.comm_buf, sF, sU, B, Btilde, mu, + gamma); + //////////////////////////// + // spin trace outer product + //////////////////////////// } - } - - Impl::InsertForce5D(mat,Btilde,Atilde,mu); - + DerivDhopComputeTime += usecond(); + Impl::InsertForce5D(mat, Btilde, Atilde, mu); } + DerivComputeTime += usecond(); } template void WilsonFermion5D::DhopDeriv( GaugeField &mat, - const FermionField &A, - const FermionField &B, - int dag) + const FermionField &A, + const FermionField &B, + int dag) { conformable(A._grid,FermionGrid()); conformable(A._grid,B._grid); @@ -317,9 +358,9 @@ void WilsonFermion5D::DhopDeriv( GaugeField &mat, template void WilsonFermion5D::DhopDerivEO(GaugeField &mat, - const FermionField &A, - const FermionField &B, - int dag) + const FermionField &A, + const FermionField &B, + int dag) { conformable(A._grid,FermionRedBlackGrid()); conformable(GaugeRedBlackGrid(),mat._grid); @@ -335,9 +376,9 @@ void WilsonFermion5D::DhopDerivEO(GaugeField &mat, template void WilsonFermion5D::DhopDerivOE(GaugeField &mat, - const FermionField &A, - const FermionField &B, - int dag) + const FermionField &A, + const FermionField &B, + int dag) { conformable(A._grid,FermionRedBlackGrid()); conformable(GaugeRedBlackGrid(),mat._grid); @@ -352,37 +393,39 @@ void WilsonFermion5D::DhopDerivOE(GaugeField &mat, template void WilsonFermion5D::DhopInternal(StencilImpl & st, LebesgueOrder &lo, - DoubledGaugeField & U, - const FermionField &in, FermionField &out,int dag) + DoubledGaugeField & U, + const FermionField &in, FermionField &out,int dag) { - Calls++; + DhopCalls++; // assert((dag==DaggerNo) ||(dag==DaggerYes)); Compressor compressor(dag); int LLs = in._grid->_rdimensions[0]; - CommTime-=usecond(); + DhopCommTime-=usecond(); st.HaloExchange(in,compressor); - CommTime+=usecond(); + DhopCommTime+=usecond(); - ComputeTime-=usecond(); + DhopComputeTime-=usecond(); // Dhop takes the 4d grid from U, and makes a 5d index for fermion - if ( dag == DaggerYes ) { -PARALLEL_FOR_LOOP - for(int ss=0;ssoSites();ss++){ - int sU=ss; - int sF=LLs*sU; - Kernels::DiracOptDhopSiteDag(st,lo,U,st.comm_buf,sF,sU,LLs,1,in,out); + if (dag == DaggerYes) { + PARALLEL_FOR_LOOP + for (int ss = 0; ss < U._grid->oSites(); ss++) { + int sU = ss; + int sF = LLs * sU; + Kernels::DiracOptDhopSiteDag(st, lo, U, st.comm_buf, sF, sU, LLs, 1, in, + out); } } else { -PARALLEL_FOR_LOOP - for(int ss=0;ssoSites();ss++){ - int sU=ss; - int sF=LLs*sU; - Kernels::DiracOptDhopSite(st,lo,U,st.comm_buf,sF,sU,LLs,1,in,out); + PARALLEL_FOR_LOOP + for (int ss = 0; ss < U._grid->oSites(); ss++) { + int sU = ss; + int sF = LLs * sU; + Kernels::DiracOptDhopSite(st, lo, U, st.comm_buf, sF, sU, LLs, 1, in, + out); } } - ComputeTime+=usecond(); + DhopComputeTime+=usecond(); } diff --git a/lib/qcd/action/fermion/WilsonFermion5D.h b/lib/qcd/action/fermion/WilsonFermion5D.h index 35016cb2..b9c35b7c 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.h +++ b/lib/qcd/action/fermion/WilsonFermion5D.h @@ -63,9 +63,14 @@ namespace Grid { void Report(void); void ZeroCounters(void); - double Calls; - double CommTime; - double ComputeTime; + double DhopCalls; + double DhopCommTime; + double DhopComputeTime; + + double DerivCalls; + double DerivCommTime; + double DerivComputeTime; + double DerivDhopComputeTime; /////////////////////////////////////////////////////////////// // Implement the abstract base diff --git a/lib/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h b/lib/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h index edb6beaa..5e3b80d9 100644 --- a/lib/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h +++ b/lib/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h @@ -131,9 +131,11 @@ namespace Grid{ Vpc.MpcDag(PhiOdd,Y); // Y= Vdag phi X=zero; ActionSolver(Mpc,Y,X); // X= (MdagM)^-1 Vdag phi - Mpc.Mpc(X,Y); // Y= Mdag^-1 Vdag phi + //Mpc.Mpc(X,Y); // Y= Mdag^-1 Vdag phi + // Multiply by Ydag + RealD action = real(innerProduct(Y,X)); - RealD action = norm2(Y); + //RealD action = norm2(Y); // The EE factorised block; normally can replace with zero if det is constant (gauge field indept) // Only really clover term that creates this. Leave the EE portion as a future to do to make most diff --git a/tests/hmc/Test_hmc_EODWFRatio.cc b/tests/hmc/Test_hmc_EODWFRatio.cc index 78f81dd9..20ef7db6 100644 --- a/tests/hmc/Test_hmc_EODWFRatio.cc +++ b/tests/hmc/Test_hmc_EODWFRatio.cc @@ -76,6 +76,12 @@ public: TheAction.push_back(Level1); Run(argc,argv); + + std::cout << GridLogMessage << "Numerator report, Pauli-Villars term : " << std::endl; + NumOp.Report(); + std::cout << GridLogMessage << "Denominator report, Dw(m) term (includes CG) : " << std::endl; + DenOp.Report(); + }; }; diff --git a/tests/solver/Test_dwf_cg_prec.cc b/tests/solver/Test_dwf_cg_prec.cc index c03998ff..d248c560 100644 --- a/tests/solver/Test_dwf_cg_prec.cc +++ b/tests/solver/Test_dwf_cg_prec.cc @@ -1,87 +1,105 @@ - /************************************************************************************* +/************************************************************************************* - Grid physics library, www.github.com/paboyle/Grid +Grid physics library, www.github.com/paboyle/Grid - Source file: ./tests/Test_dwf_cg_prec.cc +Source file: ./tests/Test_dwf_cg_prec.cc - Copyright (C) 2015 +Copyright (C) 2015 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 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. +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. +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 */ +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; -template +template struct scal { d internal; }; - Gamma::GammaMatrix Gmu [] = { - Gamma::GammaX, - Gamma::GammaY, - Gamma::GammaZ, - Gamma::GammaT - }; +Gamma::GammaMatrix Gmu[] = {Gamma::GammaX, Gamma::GammaY, Gamma::GammaZ, + Gamma::GammaT}; -int main (int argc, char ** argv) -{ - Grid_init(&argc,&argv); +int main(int argc, char** argv) { + Grid_init(&argc, &argv); - const int Ls=8; + const int Ls = 16; - GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); - GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); - GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); - GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid( + GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()), + GridDefaultMpi()); + GridRedBlackCartesian* UrbGrid = + SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridCartesian* FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid); + GridRedBlackCartesian* FrbGrid = + SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid); - std::vector seeds4({1,2,3,4}); - std::vector seeds5({5,6,7,8}); - GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); - GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + std::vector seeds4({1, 2, 3, 4}); + std::vector seeds5({5, 6, 7, 8}); + GridParallelRNG RNG5(FGrid); + RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); + RNG4.SeedFixedIntegers(seeds4); - LatticeFermion src(FGrid); random(RNG5,src); - LatticeFermion result(FGrid); result=zero; - LatticeGaugeField Umu(UGrid); + LatticeFermion src(FGrid); + random(RNG5, src); + LatticeFermion result(FGrid); + result = zero; + LatticeGaugeField Umu(UGrid); - SU3::HotConfiguration(RNG4,Umu); + SU3::HotConfiguration(RNG4, Umu); - std::vector U(4,UGrid); - for(int mu=0;mu(Umu,mu); + std::cout << GridLogMessage << "Lattice dimensions: " << GridDefaultLatt() + << " Ls: " << Ls << std::endl; + + std::vector U(4, UGrid); + for (int mu = 0; mu < Nd; mu++) { + U[mu] = PeekIndex(Umu, mu); } - - RealD mass=0.1; - RealD M5=1.8; - DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); - LatticeFermion src_o(FrbGrid); + RealD mass = 0.01; + RealD M5 = 1.8; + DomainWallFermionR Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5); + + LatticeFermion src_o(FrbGrid); LatticeFermion result_o(FrbGrid); - pickCheckerboard(Odd,src_o,src); - result_o=zero; + pickCheckerboard(Odd, src_o, src); + result_o = zero; - SchurDiagMooeeOperator HermOpEO(Ddwf); - ConjugateGradient CG(1.0e-8,10000); - CG(HermOpEO,src_o,result_o); + GridStopWatch CGTimer; + + SchurDiagMooeeOperator HermOpEO(Ddwf); + ConjugateGradient CG(1.0e-8, 10000, 0);// switch off the assert + + CGTimer.Start(); + CG(HermOpEO, src_o, result_o); + CGTimer.Stop(); + + std::cout << GridLogMessage << "Total CG time : " << CGTimer.Elapsed() + << std::endl; + + std::cout << GridLogMessage << "######## Dhop calls summary" << std::endl; + Ddwf.Report(); Grid_finalize(); } diff --git a/tests/solver/Test_wilson_cg_prec.cc b/tests/solver/Test_wilson_cg_prec.cc index 5fd2ec8c..7cc9d574 100644 --- a/tests/solver/Test_wilson_cg_prec.cc +++ b/tests/solver/Test_wilson_cg_prec.cc @@ -83,6 +83,7 @@ int main (int argc, char ** argv) SchurDiagMooeeOperator HermOpEO(Dw); ConjugateGradient CG(1.0e-8,10000); CG(HermOpEO,src_o,result_o); + Grid_finalize(); }