diff --git a/Grid/cshift/Cshift_mpi.h b/Grid/cshift/Cshift_mpi.h index 375d004e..7e93e260 100644 --- a/Grid/cshift/Cshift_mpi.h +++ b/Grid/cshift/Cshift_mpi.h @@ -122,8 +122,8 @@ template void Cshift_comms(Lattice &ret,const Lattice &r assert(shift_slice_nblock[dimension]*rhs.Grid()->_slice_block[dimension]; - cshiftVector send_buf(buffer_size); - cshiftVector recv_buf(buffer_size); + static cshiftVector send_buf; send_buf.resize(buffer_size); + static cshiftVector recv_buf; recv_buf.resize(buffer_size); int cb= (cbmask==0x2)? Odd : Even; int sshift= rhs.Grid()->CheckerBoardShiftForCB(rhs.Checkerboard(),dimension,shift,cb); @@ -198,8 +198,8 @@ template void Cshift_comms_simd(Lattice &ret,const Lattice_slice_nblock[dimension]*grid->_slice_block[dimension]; // int words = sizeof(vobj)/sizeof(vector_type); - std::vector > send_buf_extract(Nsimd); - std::vector > recv_buf_extract(Nsimd); + static std::vector > send_buf_extract; send_buf_extract.resize(Nsimd); + static std::vector > recv_buf_extract; recv_buf_extract.resize(Nsimd); scalar_object * recv_buf_extract_mpi; scalar_object * send_buf_extract_mpi; @@ -294,8 +294,8 @@ template void Cshift_comms(Lattice &ret,const Lattice &r assert(shift_slice_nblock[dimension]*rhs.Grid()->_slice_block[dimension]; - cshiftVector send_buf_v(buffer_size); - cshiftVector recv_buf_v(buffer_size); + static cshiftVector send_buf_v; send_buf_v.resize(buffer_size); + static cshiftVector recv_buf_v; recv_buf_v.resize(buffer_size); vobj *send_buf; vobj *recv_buf; { @@ -381,8 +381,8 @@ template void Cshift_comms_simd(Lattice &ret,const Lattice_slice_nblock[dimension]*grid->_slice_block[dimension]; // int words = sizeof(vobj)/sizeof(vector_type); - std::vector > send_buf_extract(Nsimd); - std::vector > recv_buf_extract(Nsimd); + static std::vector > send_buf_extract; send_buf_extract.resize(Nsimd); + static std::vector > recv_buf_extract; recv_buf_extract.resize(Nsimd); scalar_object * recv_buf_extract_mpi; scalar_object * send_buf_extract_mpi; { diff --git a/Grid/parallelIO/MetaData.h b/Grid/parallelIO/MetaData.h index d30ba523..af8b3f76 100644 --- a/Grid/parallelIO/MetaData.h +++ b/Grid/parallelIO/MetaData.h @@ -128,7 +128,7 @@ inline void MachineCharacteristics(FieldMetaData &header) std::time_t t = std::time(nullptr); std::tm tm_ = *std::localtime(&t); std::ostringstream oss; - // oss << std::put_time(&tm_, "%c %Z"); + oss << std::put_time(&tm_, "%c %Z"); header.creation_date = oss.str(); header.archive_date = header.creation_date; diff --git a/Grid/parallelIO/NerscIO.h b/Grid/parallelIO/NerscIO.h index 3ebdf0cc..99011e25 100644 --- a/Grid/parallelIO/NerscIO.h +++ b/Grid/parallelIO/NerscIO.h @@ -205,11 +205,20 @@ public: std::cout< + static inline void writeConfiguration(Lattice &Umu, + std::string file, + std::string ens_label = std::string("DWF")) + { + writeConfiguration(Umu,file,0,1,ens_label); + } template static inline void writeConfiguration(Lattice &Umu, std::string file, int two_row, - int bits32) + int bits32, + std::string ens_label = std::string("DWF")) { typedef vLorentzColourMatrixD vobj; typedef typename vobj::scalar_object sobj; @@ -219,8 +228,8 @@ public: // Following should become arguments /////////////////////////////////////////// header.sequence_number = 1; - header.ensemble_id = "UKQCD"; - header.ensemble_label = "DWF"; + header.ensemble_id = std::string("UKQCD"); + header.ensemble_label = ens_label; typedef LorentzColourMatrixD fobj3D; typedef LorentzColour2x3D fobj2D; @@ -232,7 +241,7 @@ public: GaugeStats Stats; Stats(Umu,header); MachineCharacteristics(header); - uint64_t offset; + uint64_t offset; // Sod it -- always write 3x3 double header.floating_point = std::string("IEEE64BIG"); diff --git a/Grid/qcd/action/fermion/StaggeredImpl.h b/Grid/qcd/action/fermion/StaggeredImpl.h index 8adf45a4..f44d12f4 100644 --- a/Grid/qcd/action/fermion/StaggeredImpl.h +++ b/Grid/qcd/action/fermion/StaggeredImpl.h @@ -72,19 +72,23 @@ public: StaggeredImpl(const ImplParams &p = ImplParams()) : Params(p){}; - static accelerator_inline void multLink(SiteSpinor &phi, + template + static accelerator_inline void multLink(_Spinor &phi, const SiteDoubledGaugeField &U, - const SiteSpinor &chi, + const _Spinor &chi, int mu) { - mult(&phi(), &U(mu), &chi()); + auto UU = coalescedRead(U(mu)); + mult(&phi(), &UU, &chi()); } - static accelerator_inline void multLinkAdd(SiteSpinor &phi, + template + static accelerator_inline void multLinkAdd(_Spinor &phi, const SiteDoubledGaugeField &U, - const SiteSpinor &chi, + const _Spinor &chi, int mu) { - mac(&phi(), &U(mu), &chi()); + auto UU = coalescedRead(U(mu)); + mac(&phi(), &UU, &chi()); } template diff --git a/Grid/qcd/action/fermion/WilsonImpl.h b/Grid/qcd/action/fermion/WilsonImpl.h index 94676b6b..2ff6feba 100644 --- a/Grid/qcd/action/fermion/WilsonImpl.h +++ b/Grid/qcd/action/fermion/WilsonImpl.h @@ -184,18 +184,22 @@ public: mat = TraceIndex(P); } - inline void extractLinkField(std::vector &mat, DoubledGaugeField &Uds){ + inline void extractLinkField(std::vector &mat, DoubledGaugeField &Uds) + { for (int mu = 0; mu < Nd; mu++) mat[mu] = PeekIndex(Uds, mu); } - - inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField Ã,int mu){ - + inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField Ã,int mu) + { +#undef USE_OLD_INSERT_FORCE int Ls=Btilde.Grid()->_fdimensions[0]; + autoView( mat_v , mat, AcceleratorWrite); +#ifdef USE_OLD_INSERT_FORCE GaugeLinkField tmp(mat.Grid()); tmp = Zero(); { + const int Nsimd = SiteSpinor::Nsimd(); autoView( tmp_v , tmp, AcceleratorWrite); autoView( Btilde_v , Btilde, AcceleratorRead); autoView( Atilde_v , Atilde, AcceleratorRead); @@ -208,6 +212,29 @@ public: }); } PokeIndex(mat,tmp,mu); +#else + { + const int Nsimd = SiteSpinor::Nsimd(); + autoView( Btilde_v , Btilde, AcceleratorRead); + autoView( Atilde_v , Atilde, AcceleratorRead); + accelerator_for(sss,mat.Grid()->oSites(),Nsimd,{ + int sU=sss; + typedef decltype(coalescedRead(mat_v[sU](mu)() )) ColorMatrixType; + ColorMatrixType sum; + zeroit(sum); + for(int s=0;s::SeqConservedCurrent(PropagatorField &q_in, } std::vector G_s(Ls,1.0); + Integer sign = 1; // sign flip for vector/tadpole if ( curr_type == Current::Axial ) { for(int s=0;s_b; + auto c=this->_c; + if ( b == 1 && c == 0 ) { + sign = -1; + } + else { + std::cerr << "Error: Tadpole implementation currently unavailable for non-Shamir actions." << std::endl; + assert(b==1 && c==0); + } + } for(int s=0;s::SeqConservedCurrent(PropagatorField &q_in, tmp = Cshift(tmp,mu,1); Impl::multLinkField(Utmp,this->Umu,tmp,mu); - tmp = G_s[s]*( Utmp*ph - gmu*Utmp*ph ); // Forward hop + tmp = sign*G_s[s]*( Utmp*ph - gmu*Utmp*ph ); // Forward hop tmp = where((lcoor>=tmin),tmp,zz); // Mask the time L_Q = where((lcoor<=tmax),tmp,zz); // Position of current complicated diff --git a/Grid/qcd/action/fermion/implementation/StaggeredKernelsImplementation.h b/Grid/qcd/action/fermion/implementation/StaggeredKernelsImplementation.h index 0b6f9fb0..dd62e109 100644 --- a/Grid/qcd/action/fermion/implementation/StaggeredKernelsImplementation.h +++ b/Grid/qcd/action/fermion/implementation/StaggeredKernelsImplementation.h @@ -35,39 +35,32 @@ NAMESPACE_BEGIN(Grid); #define GENERIC_STENCIL_LEG(U,Dir,skew,multLink) \ SE = st.GetEntry(ptype, Dir+skew, sF); \ if (SE->_is_local ) { \ - if (SE->_permute) { \ - chi_p = χ \ - permute(chi, in[SE->_offset], ptype); \ - } else { \ - chi_p = &in[SE->_offset]; \ - } \ + int perm= SE->_permute; \ + chi = coalescedReadPermute(in[SE->_offset],ptype,perm,lane);\ } else { \ - chi_p = &buf[SE->_offset]; \ + chi = coalescedRead(buf[SE->_offset],lane); \ } \ - multLink(Uchi, U[sU], *chi_p, Dir); + acceleratorSynchronise(); \ + multLink(Uchi, U[sU], chi, Dir); #define GENERIC_STENCIL_LEG_INT(U,Dir,skew,multLink) \ SE = st.GetEntry(ptype, Dir+skew, sF); \ if (SE->_is_local ) { \ - if (SE->_permute) { \ - chi_p = χ \ - permute(chi, in[SE->_offset], ptype); \ - } else { \ - chi_p = &in[SE->_offset]; \ - } \ + int perm= SE->_permute; \ + chi = coalescedReadPermute(in[SE->_offset],ptype,perm,lane);\ } else if ( st.same_node[Dir] ) { \ - chi_p = &buf[SE->_offset]; \ + chi = coalescedRead(buf[SE->_offset],lane); \ } \ if (SE->_is_local || st.same_node[Dir] ) { \ - multLink(Uchi, U[sU], *chi_p, Dir); \ + multLink(Uchi, U[sU], chi, Dir); \ } #define GENERIC_STENCIL_LEG_EXT(U,Dir,skew,multLink) \ SE = st.GetEntry(ptype, Dir+skew, sF); \ if ((!SE->_is_local) && (!st.same_node[Dir]) ) { \ nmu++; \ - chi_p = &buf[SE->_offset]; \ - multLink(Uchi, U[sU], *chi_p, Dir); \ + chi = coalescedRead(buf[SE->_offset],lane); \ + multLink(Uchi, U[sU], chi, Dir); \ } template @@ -84,12 +77,14 @@ void StaggeredKernels::DhopSiteGeneric(StencilView &st, SiteSpinor *buf, int sF, int sU, const FermionFieldView &in, FermionFieldView &out, int dag) { - const SiteSpinor *chi_p; - SiteSpinor chi; - SiteSpinor Uchi; + typedef decltype(coalescedRead(in[0])) calcSpinor; + calcSpinor chi; + calcSpinor Uchi; StencilEntry *SE; int ptype; int skew; + const int Nsimd = SiteHalfSpinor::Nsimd(); + const int lane=acceleratorSIMTlane(Nsimd); // for(int s=0;s::DhopSiteGeneric(StencilView &st, if ( dag ) { Uchi = - Uchi; } - vstream(out[sF], Uchi); + coalescedWrite(out[sF], Uchi,lane); } }; @@ -130,13 +125,16 @@ template accelerator_inline void StaggeredKernels::DhopSiteGenericInt(StencilView &st, DoubledGaugeFieldView &U, DoubledGaugeFieldView &UUU, SiteSpinor *buf, int sF, int sU, - const FermionFieldView &in, FermionFieldView &out,int dag) { - const SiteSpinor *chi_p; - SiteSpinor chi; - SiteSpinor Uchi; + const FermionFieldView &in, FermionFieldView &out,int dag) +{ + typedef decltype(coalescedRead(in[0])) calcSpinor; + calcSpinor chi; + calcSpinor Uchi; StencilEntry *SE; int ptype; int skew ; + const int Nsimd = SiteHalfSpinor::Nsimd(); + const int lane=acceleratorSIMTlane(Nsimd); // for(int s=0;s::DhopSiteGenericInt(StencilView &st, if ( dag ) { Uchi = - Uchi; } - vstream(out[sF], Uchi); + coalescedWrite(out[sF], Uchi,lane); } }; @@ -178,14 +176,17 @@ template accelerator_inline void StaggeredKernels::DhopSiteGenericExt(StencilView &st, DoubledGaugeFieldView &U, DoubledGaugeFieldView &UUU, SiteSpinor *buf, int sF, int sU, - const FermionFieldView &in, FermionFieldView &out,int dag) { - const SiteSpinor *chi_p; - // SiteSpinor chi; - SiteSpinor Uchi; + const FermionFieldView &in, FermionFieldView &out,int dag) +{ + typedef decltype(coalescedRead(in[0])) calcSpinor; + calcSpinor chi; + calcSpinor Uchi; StencilEntry *SE; int ptype; int nmu=0; int skew ; + const int Nsimd = SiteHalfSpinor::Nsimd(); + const int lane=acceleratorSIMTlane(Nsimd); // for(int s=0;s::DhopSiteGenericExt(StencilView &st, GENERIC_STENCIL_LEG_EXT(UUU,Zm,skew,Impl::multLinkAdd); GENERIC_STENCIL_LEG_EXT(UUU,Tm,skew,Impl::multLinkAdd); } - if ( nmu ) { - if ( dag ) { - out[sF] = out[sF] - Uchi; + if ( nmu ) { + auto _out = coalescedRead(out[sF],lane); + if ( dag ) { + coalescedWrite(out[sF], _out-Uchi,lane); } else { - out[sF] = out[sF] + Uchi; + coalescedWrite(out[sF], _out+Uchi,lane); } } } @@ -261,6 +263,8 @@ void StaggeredKernels::DhopImproved(StencilImpl &st, LebesgueOrder &lo, GridBase *FGrid=in.Grid(); GridBase *UGrid=U.Grid(); typedef StaggeredKernels ThisKernel; + const int Nsimd = SiteHalfSpinor::Nsimd(); + const int lane=acceleratorSIMTlane(Nsimd); autoView( UUU_v , UUU, AcceleratorRead); autoView( U_v , U, AcceleratorRead); autoView( in_v , in, AcceleratorRead); @@ -301,6 +305,8 @@ void StaggeredKernels::DhopNaive(StencilImpl &st, LebesgueOrder &lo, GridBase *FGrid=in.Grid(); GridBase *UGrid=U.Grid(); typedef StaggeredKernels ThisKernel; + const int Nsimd = SiteHalfSpinor::Nsimd(); + const int lane=acceleratorSIMTlane(Nsimd); autoView( UUU_v , U, AcceleratorRead); autoView( U_v , U, AcceleratorRead); autoView( in_v , in, AcceleratorRead); diff --git a/Grid/qcd/smearing/StoutSmearing.h b/Grid/qcd/smearing/StoutSmearing.h index ed2ccdb6..6ee78e8c 100644 --- a/Grid/qcd/smearing/StoutSmearing.h +++ b/Grid/qcd/smearing/StoutSmearing.h @@ -85,21 +85,18 @@ public: std::cout << GridLogDebug << "Stout smearing started\n"; - // Smear the configurations + // C contains the staples multiplied by some rho + u_smr = U ; // set the smeared field to the current gauge field SmearBase->smear(C, U); for (int mu = 0; mu < Nd; mu++) { - if( mu == OrthogDim ) - tmp = 1.0; // Don't smear in the orthogonal direction - else { - tmp = peekLorentz(C, mu); - Umu = peekLorentz(U, mu); - iq_mu = Ta( - tmp * - adj(Umu)); // iq_mu = Ta(Omega_mu) to match the signs with the paper - exponentiate_iQ(tmp, iq_mu); - } - pokeLorentz(u_smr, tmp * Umu, mu); // u_smr = exp(iQ_mu)*U_mu + if( mu == OrthogDim ) continue ; + // u_smr = exp(iQ_mu)*U_mu apart from Orthogdim + Umu = peekLorentz(U, mu); + tmp = peekLorentz(C, mu); + iq_mu = Ta( tmp * adj(Umu)); + exponentiate_iQ(tmp, iq_mu); + pokeLorentz(u_smr, tmp * Umu, mu); } std::cout << GridLogDebug << "Stout smearing completed\n"; }; diff --git a/Grid/qcd/utils/Metric.h b/Grid/qcd/utils/Metric.h index 10d06de8..d8ae27dc 100644 --- a/Grid/qcd/utils/Metric.h +++ b/Grid/qcd/utils/Metric.h @@ -93,13 +93,13 @@ public: GeneralisedMomenta(GridBase* grid, Metric& M): M(M), Mom(grid), AuxMom(grid), AuxField(grid){} // Correct - void MomentaDistribution(GridParallelRNG& pRNG){ + void MomentaDistribution(GridSerialRNG & sRNG, GridParallelRNG& pRNG){ // Generate a distribution for // P^dag G P // where G = M^-1 // Generate gaussian momenta - Implementation::generate_momenta(Mom, pRNG); + Implementation::generate_momenta(Mom, sRNG, pRNG); // Modify the distribution with the metric M.MSquareRoot(Mom); @@ -107,8 +107,8 @@ public: // Auxiliary momenta // do nothing if trivial, so hide in the metric MomentaField AuxMomTemp(Mom.Grid()); - Implementation::generate_momenta(AuxMom, pRNG); - Implementation::generate_momenta(AuxField, pRNG); + Implementation::generate_momenta(AuxMom, sRNG, pRNG); + Implementation::generate_momenta(AuxField, sRNG, pRNG); // Modify the distribution with the metric // Aux^dag M Aux M.MInvSquareRoot(AuxMom); // AuxMom = M^{-1/2} AuxMomTemp diff --git a/Grid/tensors/Tensor_outer.h b/Grid/tensors/Tensor_outer.h index 4902c22f..a32a2a91 100644 --- a/Grid/tensors/Tensor_outer.h +++ b/Grid/tensors/Tensor_outer.h @@ -34,6 +34,16 @@ NAMESPACE_BEGIN(Grid); // outerProduct Scalar x Scalar -> Scalar // Vector x Vector -> Matrix /////////////////////////////////////////////////////////////////////////////////////// +template = 0> +accelerator_inline CC outerProduct(const CC &l, const CC& r) +{ + return l*conj(r); +} +template = 0> +accelerator_inline RR outerProduct(const RR &l, const RR& r) +{ + return l*r; +} template accelerator_inline auto outerProduct (const iVector& lhs,const iVector& rhs) -> iMatrix @@ -57,17 +67,6 @@ auto outerProduct (const iScalar& lhs,const iScalar& rhs) -> iScalar = 0> -accelerator_inline CC outerProduct(const CC &l, const CC& r) -{ - return l*conj(r); -} -template = 0> -accelerator_inline RR outerProduct(const RR &l, const RR& r) -{ - return l*r; -} - NAMESPACE_END(Grid); #endif diff --git a/Grid/util/Init.cc b/Grid/util/Init.cc index 55d8c5bf..bfbc464d 100644 --- a/Grid/util/Init.cc +++ b/Grid/util/Init.cc @@ -153,6 +153,9 @@ void GridCmdOptionIntVector(const std::string &str,VectorInt & vec) return; } +template void GridCmdOptionIntVector(const std::string &str,std::vector & vec); +template void GridCmdOptionIntVector(const std::string &str,Coordinate & vec); + void GridCmdOptionInt(std::string &str,int & val) { std::stringstream ss(str); diff --git a/HMC/Mobius2p1fRHMC.cc b/HMC/Mobius2p1fRHMC.cc index 82ca4d37..b958d548 100644 --- a/HMC/Mobius2p1fRHMC.cc +++ b/HMC/Mobius2p1fRHMC.cc @@ -56,12 +56,12 @@ int main(int argc, char **argv) { MD.trajL = 1.0; HMCparameters HMCparams; - HMCparams.StartTrajectory = 30; + HMCparams.StartTrajectory = 0; HMCparams.Trajectories = 200; HMCparams.NoMetropolisUntil= 0; // "[HotStart, ColdStart, TepidStart, CheckpointStart]\n"; - // HMCparams.StartingType =std::string("ColdStart"); - HMCparams.StartingType =std::string("CheckpointStart"); + HMCparams.StartingType =std::string("ColdStart"); + // HMCparams.StartingType =std::string("CheckpointStart"); HMCparams.MD = MD; HMCWrapper TheHMC(HMCparams); diff --git a/tests/debug/Test_heatbath_dwf_eofa.cc b/tests/debug/Test_heatbath_dwf_eofa.cc index 9d453a96..e1c18021 100644 --- a/tests/debug/Test_heatbath_dwf_eofa.cc +++ b/tests/debug/Test_heatbath_dwf_eofa.cc @@ -66,7 +66,9 @@ int main(int argc, char** argv) // Set up RNGs std::vector seeds4({1, 2, 3, 4}); std::vector seeds5({5, 6, 7, 8}); + GridSerialRNG sRNG; GridParallelRNG RNG5(FGrid); + sRNG.SeedFixedIntegers(seeds5); RNG5.SeedFixedIntegers(seeds5); GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); @@ -84,7 +86,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, false); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu,sRNG, RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } @@ -94,7 +96,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, true); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu,sRNG, RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } diff --git a/tests/debug/Test_heatbath_dwf_eofa_gparity.cc b/tests/debug/Test_heatbath_dwf_eofa_gparity.cc index 22cc1e90..7eabfc65 100644 --- a/tests/debug/Test_heatbath_dwf_eofa_gparity.cc +++ b/tests/debug/Test_heatbath_dwf_eofa_gparity.cc @@ -74,6 +74,9 @@ int main(int argc, char** argv) RNG5.SeedFixedIntegers(seeds5); GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + GridSerialRNG sRNG; + RNG4.SeedFixedIntegers(seeds4); + sRNG.SeedFixedIntegers(seeds5); // Random gauge field LatticeGaugeField Umu(UGrid); @@ -90,7 +93,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, false); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu,sRNG, RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } @@ -100,7 +103,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, true); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu,sRNG, RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } diff --git a/tests/debug/Test_heatbath_mobius_eofa.cc b/tests/debug/Test_heatbath_mobius_eofa.cc index 4cf4bf53..48806642 100644 --- a/tests/debug/Test_heatbath_mobius_eofa.cc +++ b/tests/debug/Test_heatbath_mobius_eofa.cc @@ -68,8 +68,10 @@ int main(int argc, char** argv) // Set up RNGs std::vector seeds4({1, 2, 3, 4}); std::vector seeds5({5, 6, 7, 8}); + GridSerialRNG sRNG; GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + sRNG.SeedFixedIntegers(seeds5); GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); @@ -86,7 +88,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, false); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu, sRNG,RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } @@ -96,7 +98,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, true); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu, sRNG,RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } diff --git a/tests/debug/Test_heatbath_mobius_eofa_gparity.cc b/tests/debug/Test_heatbath_mobius_eofa_gparity.cc index 2fcb4b9f..52447e5e 100644 --- a/tests/debug/Test_heatbath_mobius_eofa_gparity.cc +++ b/tests/debug/Test_heatbath_mobius_eofa_gparity.cc @@ -73,7 +73,9 @@ int main(int argc, char** argv) std::vector seeds4({1, 2, 3, 4}); std::vector seeds5({5, 6, 7, 8}); GridParallelRNG RNG5(FGrid); + GridSerialRNG sRNG; RNG5.SeedFixedIntegers(seeds5); + sRNG.SeedFixedIntegers(seeds5); GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); @@ -91,7 +93,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, false); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu, sRNG, RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } @@ -101,7 +103,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, true); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu, sRNG, RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } diff --git a/tests/forces/Test_dwf_force_eofa.cc b/tests/forces/Test_dwf_force_eofa.cc index 80d36934..525178d0 100644 --- a/tests/forces/Test_dwf_force_eofa.cc +++ b/tests/forces/Test_dwf_force_eofa.cc @@ -86,7 +86,9 @@ int main (int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, CG, CG, CG, CG, Params, true); - Meofa.refresh(U, RNG5); + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds4); + Meofa.refresh(U, sRNG, RNG5 ); + RealD S = Meofa.S(U); // pdag M p // get the deriv of phidag M phi with respect to "U" diff --git a/tests/forces/Test_dwf_gpforce.cc b/tests/forces/Test_dwf_gpforce.cc index 28133cc6..1fa1c6e4 100644 --- a/tests/forces/Test_dwf_gpforce.cc +++ b/tests/forces/Test_dwf_gpforce.cc @@ -84,6 +84,13 @@ int main (int argc, char ** argv) GparityDomainWallFermionR::ImplParams params; params.twists = twists; + /* + params.boundary_phases[0] = 1.0; + params.boundary_phases[1] = 1.0; + params.boundary_phases[2] = 1.0; + params.boundary_phases[3] =- 1.0; + */ + GparityDomainWallFermionR Dw(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,params); Dw.M (phi,Mphi); @@ -96,6 +103,16 @@ int main (int argc, char ** argv) Dw.MDeriv(tmp , Mphi, phi,DaggerNo ); UdSdU=tmp; Dw.MDeriv(tmp , phi, Mphi,DaggerYes ); UdSdU=(UdSdU+tmp); + + // ***************************************************************************************** + // *** There is a funny negative sign in all derivatives. This is - UdSdU. *** + // *** *** + // *** Deriv in both Wilson gauge action and the TwoFlavour.h seems to miss a minus sign *** + // *** UdSdU is negated relative to what I think - call what is returned mUdSdU, *** + // *** and insert minus sign *** + // ***************************************************************************************** + + UdSdU = - UdSdU ; // Follow sign convention of actions in Grid. Seems crazy. FermionField Ftmp (FGrid); @@ -106,7 +123,7 @@ int main (int argc, char ** argv) RealD Hmom = 0.0; RealD Hmomprime = 0.0; LatticeColourMatrix mommu(UGrid); - LatticeColourMatrix forcemu(UGrid); + LatticeColourMatrix mUdSdUmu(UGrid); LatticeGaugeField mom(UGrid); LatticeGaugeField Uprime(UGrid); @@ -114,10 +131,20 @@ int main (int argc, char ** argv) SU::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg - Hmom -= real(sum(trace(mommu*mommu))); + // Momentum Hamiltonian is - trace(p^2)/HMC_MOM_DENOMINATOR + // + // Integrator.h: RealD H = - FieldImplementation::FieldSquareNorm(P)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom // GaugeImplTypes.h: Hloc += trace(Pmu * Pmu); + // Sign comes from a sneaky multiply by "i" in GaussianFundemantalLie algebra + // P is i P^a_\mu T^a, not Pa Ta + // + // Integrator.h: H = Hmom + sum S(action) + Hmom -= real(sum(trace(mommu*mommu)))/ HMC_MOMENTUM_DENOMINATOR; PokeIndex(mom,mommu,mu); + // -- Drops factor of "i" in the U update: U' = e^{P dt} U [ _not_ e^{iPdt}U ]. P is anti hermitian already + // -- Udot = p U + // fourth order exponential approx autoView( mom_v, mom, CpuRead); autoView( U_v , U, CpuRead); @@ -134,8 +161,8 @@ int main (int argc, char ** argv) ; }); } - std::cout << GridLogMessage <<"Initial mom hamiltonian is "<< Hmom <(mom,mu); - std::cout << GridLogMessage<< " Mommu " << norm2(mommu)<(UdSdU,mu); - std::cout << GridLogMessage<< " dsdumu " << norm2(mommu)<(UdSdU,mu); - mommu=Ta(mommu)*2.0; + mommu=Ta(mommu); // projectForce , GaugeImplTypes.h PokeIndex(UdSdU,mommu,mu); } for(int mu=0;mu(mom,mu); - std::cout << GridLogMessage<< " Mommu " << norm2(mommu)<(UdSdU,mu); - std::cout << GridLogMessage<< " dsdumu " << norm2(mommu)<(UdSdU,mu); + mUdSdUmu= PeekIndex(UdSdU,mu); mommu = PeekIndex(mom,mu); - // Update PF action density - dS = dS+trace(mommu*forcemu)*dt; + // + // Derive HMC eom: + // + // Sdot = - 2 trace( p p^dot ) / D - trace( p [ mUdSdU - h.c. ] ) = 0 + // + // + // Sdot = 0 = - 2 trace( p p^dot ) / D - 2 trace( p Ta( mUdSdU ) = 0 + // + // EOM: + // + // pdot = - D Ta( mUdSdU ) -- source of sign is the "funny sign" above + // + // dSqcd_dt = - 2.0*trace(mommu* Ta(mUdSdU) )*dt -- i.e. mUdSdU with adjoint term -> force has a 2x implicit + // + // dH_mom/dt = - 2 trace (p pdot)/Denom + // + // dH_tot / dt = 0 <= pdot = - Denom * mUdSdU + // + // dH_mom/dt = 2 trace (p mUdSdU ) + // + // True Momentum delta H has a dt^2 piece + // + // dSmom = [ trace mom*mom - trace ( (mom-Denom*f*dt)(mom-Denom*f*dt) ) ] / Denom + // = 2*trace(mom*f) dt - Denom*dt*dt * trace(f*f). + // = dSmom + dSmom2 + // - dSmom = dSmom - trace(mommu*forcemu) * dt; - dSmom2 = dSmom2 - trace(forcemu*forcemu) *(0.25* dt*dt); + dS = dS - 2.0*trace(mommu*mUdSdUmu)*dt; // U and Udagger derivs hence 2x. - // Update mom action density - mommu = mommu + forcemu*(dt*0.5); + dSmom = dSmom + 2.0*trace(mommu*mUdSdUmu) * dt; // this 2.0 coms from derivative of p^2 + + dSmom2 = dSmom2 - trace(mUdSdUmu*mUdSdUmu) * dt*dt* HMC_MOMENTUM_DENOMINATOR; // Remnant - Hmomprime -= real(sum(trace(mommu*mommu))); + // Update mom action density . Verbatim update_P in Integrator.h + mommu = mommu - mUdSdUmu * dt* HMC_MOMENTUM_DENOMINATOR;; + + Hmomprime -= real(sum(trace(mommu*mommu))) / HMC_MOMENTUM_DENOMINATOR; } @@ -199,20 +233,25 @@ int main (int argc, char ** argv) ComplexD dSm = sum(dSmom); ComplexD dSm2 = sum(dSmom2); + std::cout << GridLogMessage <<"dSm "<< dSm< CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, CG, CG, CG, CG, Params, true); - Meofa.refresh(U, RNG5); + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds4); + Meofa.refresh(U, sRNG, RNG5); RealD S = Meofa.S(U); // pdag M p // get the deriv of phidag M phi with respect to "U" diff --git a/tests/forces/Test_gp_rect_force.cc b/tests/forces/Test_gp_rect_force.cc index 98ebb2fa..e277ea6b 100644 --- a/tests/forces/Test_gp_rect_force.cc +++ b/tests/forces/Test_gp_rect_force.cc @@ -29,7 +29,6 @@ Author: paboyle using namespace std; using namespace Grid; - ; @@ -59,6 +58,10 @@ int main (int argc, char ** argv) double beta = 1.0; double c1 = 0.331; + const int nu = 1; + std::vector twists(Nd,0); + twists[nu] = 1; + ConjugateGimplD::setDirections(twists); ConjugatePlaqPlusRectangleActionR Action(beta,c1); //ConjugateWilsonGaugeActionR Action(beta); //WilsonGaugeActionR Action(beta); diff --git a/tests/forces/Test_laplacian_force.cc b/tests/forces/Test_laplacian_force.cc index 18508860..dbaf1cbd 100644 --- a/tests/forces/Test_laplacian_force.cc +++ b/tests/forces/Test_laplacian_force.cc @@ -46,6 +46,7 @@ int main (int argc, char ** argv) std::vector seeds({1,2,3,4}); + GridSerialRNG sRNG; sRNG.SeedFixedIntegers({4,5,6,7}); GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(std::vector({15,91,21,3})); @@ -67,7 +68,7 @@ int main (int argc, char ** argv) LaplacianAdjointField Laplacian(&Grid, CG, LapPar, Kappa); GeneralisedMomenta LaplacianMomenta(&Grid, Laplacian); LaplacianMomenta.M.ImportGauge(U); - LaplacianMomenta.MomentaDistribution(pRNG);// fills the Momenta with the correct distr + LaplacianMomenta.MomentaDistribution(sRNG,pRNG);// fills the Momenta with the correct distr std::cout << std::setprecision(15); diff --git a/tests/forces/Test_mobius_force.cc b/tests/forces/Test_mobius_force.cc index ba7bc363..d2326a81 100644 --- a/tests/forces/Test_mobius_force.cc +++ b/tests/forces/Test_mobius_force.cc @@ -69,7 +69,14 @@ int main (int argc, char ** argv) RealD M5=1.8; RealD b=0.5; RealD c=0.5; - MobiusFermionR Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c); + + WilsonImplParams p; + p.boundary_phases[0] = 1.0; + p.boundary_phases[1] = 1.0; + p.boundary_phases[2] = 1.0; + p.boundary_phases[3] =- 1.0; + + MobiusFermionR Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c,p); Ddwf.M (phi,Mphi); ComplexD S = innerProduct(Mphi,Mphi); // pdag MdagM p @@ -82,24 +89,44 @@ int main (int argc, char ** argv) Ddwf.MDeriv(tmp , Mphi, phi,DaggerNo ); UdSdU=tmp; Ddwf.MDeriv(tmp , phi, Mphi,DaggerYes ); UdSdU=(UdSdU+tmp); + // ***************************************************************************************** + // *** There is a funny negative sign in all derivatives. This is - UdSdU. *** + // *** *** + // *** Deriv in both Wilson gauge action and the TwoFlavour.h seems to miss a minus sign *** + // *** UdSdU is negated relative to what I think - call what is returned mUdSdU, *** + // *** and insert minus sign *** + // ***************************************************************************************** + + UdSdU = - UdSdU ; // Follow sign convention of actions in Grid. Seems crazy. + LatticeFermion Ftmp (FGrid); //////////////////////////////////// // Modify the gauge field a little //////////////////////////////////// - RealD dt = 0.0001; + RealD dt = 0.001; + RealD Hmom = 0.0; + RealD Hmomprime = 0.0; LatticeColourMatrix mommu(UGrid); - LatticeColourMatrix forcemu(UGrid); + LatticeColourMatrix mUdSdUmu(UGrid); LatticeGaugeField mom(UGrid); LatticeGaugeField Uprime(UGrid); for(int mu=0;mu::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg - PokeIndex(mom,mommu,mu); + // Momentum Hamiltonian is - trace(p^2)/HMC_MOM_DENOMINATOR + // + // Integrator.h: RealD H = - FieldImplementation::FieldSquareNorm(P)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom // GaugeImplTypes.h: Hloc += trace(Pmu * Pmu); + // Sign comes from a sneaky multiply by "i" in GaussianFundemantalLie algebra + // P is i P^a_\mu T^a, not Pa Ta + // + // Integrator.h: H = Hmom + sum S(action) + Hmom -= real(sum(trace(mommu*mommu)))/ HMC_MOMENTUM_DENOMINATOR; + // fourth order exponential approx autoView( U_v , U, CpuRead); autoView( mom_v, mom, CpuRead); @@ -115,6 +142,7 @@ int main (int argc, char ** argv) ; }); } + std::cout << GridLogMessage <<"Initial mom hamiltonian is "<< Hmom <(UdSdU,mu); - mommu=Ta(mommu)*2.0; + mommu=Ta(mommu); PokeIndex(UdSdU,mommu,mu); } for(int mu=0;mu(UdSdU,mu); + + mUdSdUmu= PeekIndex(UdSdU,mu); mommu = PeekIndex(mom,mu); - // Update PF action density - dS = dS+trace(mommu*forcemu)*dt; + // + // Derive HMC eom: + // + // Sdot = - 2 trace( p p^dot ) / D - trace( p [ mUdSdU - h.c. ] ) = 0 + // + // + // Sdot = 0 = - 2 trace( p p^dot ) / D - 2 trace( p Ta( mUdSdU ) = 0 + // + // EOM: + // + // pdot = - D Ta( mUdSdU ) -- source of sign is the "funny sign" above + // + // dSqcd_dt = - 2.0*trace(mommu* Ta(mUdSdU) )*dt -- i.e. mUdSdU with adjoint term -> force has a 2x implicit + // + // dH_mom/dt = - 2 trace (p pdot)/Denom + // + // dH_tot / dt = 0 <= pdot = - Denom * mUdSdU + // + // dH_mom/dt = 2 trace (p mUdSdU ) + // + // True Momentum delta H has a dt^2 piece + // + // dSmom = [ trace mom*mom - trace ( (mom-Denom*f*dt)(mom-Denom*f*dt) ) ] / Denom + // = 2*trace(mom*f) dt - Denom*dt*dt * trace(f*f). + // = dSmom + dSmom2 + // + + dS = dS - 2.0*trace(mommu*mUdSdUmu)*dt; // U and Udagger derivs hence 2x. + + dSmom = dSmom + 2.0*trace(mommu*mUdSdUmu) * dt; // this 2.0 coms from derivative of p^2 + + dSmom2 = dSmom2 - trace(mUdSdUmu*mUdSdUmu) * dt*dt* HMC_MOMENTUM_DENOMINATOR; // Remnant + + mommu = mommu - mUdSdUmu * dt* HMC_MOMENTUM_DENOMINATOR;; + + Hmomprime -= real(sum(trace(mommu*mommu))) / HMC_MOMENTUM_DENOMINATOR; } ComplexD dSpred = sum(dS); + ComplexD dSm = sum(dSmom); + ComplexD dSm2 = sum(dSmom2); - std::cout << GridLogMessage << " -- S "< CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, CG, CG, CG, CG, Params, false); - Meofa.refresh(U, RNG5); + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds4); + Meofa.refresh(U, sRNG, RNG5 ); RealD S = Meofa.S(U); // pdag M p // get the deriv of phidag M phi with respect to "U" diff --git a/tests/forces/Test_mobius_gpforce_eofa.cc b/tests/forces/Test_mobius_gpforce_eofa.cc index 9c80b2aa..7f114615 100644 --- a/tests/forces/Test_mobius_gpforce_eofa.cc +++ b/tests/forces/Test_mobius_gpforce_eofa.cc @@ -93,7 +93,8 @@ int main (int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, CG, CG, CG, CG, Params, false); - Meofa.refresh(U, RNG5); + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds4); + Meofa.refresh(U, sRNG, RNG5 ); RealD S = Meofa.S(U); // pdag M p // get the deriv of phidag M phi with respect to "U" diff --git a/tests/forces/Test_momentum_filter.cc b/tests/forces/Test_momentum_filter.cc index 856ea0f2..794b5fa0 100644 --- a/tests/forces/Test_momentum_filter.cc +++ b/tests/forces/Test_momentum_filter.cc @@ -61,7 +61,9 @@ int main (int argc, char ** argv) std::vector seeds({1,2,3,4}); GridParallelRNG pRNG(&Grid); + GridSerialRNG sRNG; pRNG.SeedFixedIntegers(seeds); + sRNG.SeedFixedIntegers(seeds); typedef PeriodicGimplR Gimpl; typedef WilsonGaugeAction GaugeAction; @@ -115,7 +117,7 @@ int main (int argc, char ** argv) integrator.setMomentumFilter(filter); - integrator.refresh(U, pRNG); //doesn't actually change the gauge field + integrator.refresh(U, sRNG, pRNG); //doesn't actually change the gauge field //Check the momentum is zero on the boundary const auto &P = integrator.getMomentum();