From d9d4c5916a3e77c7563175ed3042695643d216b4 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sun, 26 Jul 2015 10:54:38 +0900 Subject: [PATCH] Elemental force term for Wilson dslash added and tests thereof passing. Now need to construct pseudofermion two flavour, ratio, one flavour, ratio action fragments. --- lib/GridConfig.h.in | 3 - lib/lattice/Lattice_base.h | 5 + lib/qcd/action/fermion/FermionOperator.h | 19 ++- lib/qcd/action/fermion/WilsonFermion.cc | 88 +++++++++++- lib/qcd/action/fermion/WilsonFermion.h | 88 +++++++++++- lib/qcd/action/fermion/WilsonFermion5D.cc | 95 ++++++++++++- lib/qcd/action/fermion/WilsonFermion5D.h | 18 ++- lib/qcd/action/fermion/WilsonKernels.cc | 44 +++--- lib/qcd/action/fermion/WilsonKernels.h | 2 +- lib/qcd/hmc/HMC.h | 2 - lib/simd/Grid_vector_types.h | 2 +- lib/tensors/Tensor_outer.h | 6 +- lib/tensors/Tensor_trace.h | 10 ++ tests/Make.inc | 10 +- tests/Test_dwf_cg_prec.cc | 4 +- tests/Test_wilson_force.cc | 102 ++++++++++++++ tests/Test_wilson_force_phiMdagMphi.cc | 140 +++++++++++++++++++ tests/Test_wilson_force_phiMphi.cc | 162 ++++++++++++++++++++++ 18 files changed, 758 insertions(+), 42 deletions(-) create mode 100644 tests/Test_wilson_force.cc create mode 100644 tests/Test_wilson_force_phiMdagMphi.cc create mode 100644 tests/Test_wilson_force_phiMphi.cc diff --git a/lib/GridConfig.h.in b/lib/GridConfig.h.in index 1ec0ecf5..bcb8d811 100644 --- a/lib/GridConfig.h.in +++ b/lib/GridConfig.h.in @@ -33,9 +33,6 @@ /* Support AVX2 (Advanced Vector Extensions 2) instructions */ #undef HAVE_AVX2 -/* define if the compiler supports basic C++11 syntax */ -#undef HAVE_CXX11 - /* Define to 1 if you have the declaration of `be64toh', and to 0 if you don't. */ #undef HAVE_DECL_BE64TOH diff --git a/lib/lattice/Lattice_base.h b/lib/lattice/Lattice_base.h index 2a365cdd..9c9ef044 100644 --- a/lib/lattice/Lattice_base.h +++ b/lib/lattice/Lattice_base.h @@ -60,6 +60,11 @@ public: GridBase *_grid; int checkerboard; std::vector > _odata; + + // to pthread need a computable loop where loop induction is not required + int begin(void) { return 0;}; + int end(void) { return _odata.size(); } + vobj & operator[](int i) { return _odata[i]; }; public: typedef typename vobj::scalar_type scalar_type; diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h index 1aa2f8f0..a8f10380 100644 --- a/lib/qcd/action/fermion/FermionOperator.h +++ b/lib/qcd/action/fermion/FermionOperator.h @@ -39,10 +39,23 @@ namespace Grid { virtual void Dhop (const FermionField &in, FermionField &out,int dag)=0; virtual void DhopOE(const FermionField &in, FermionField &out,int dag)=0; virtual void DhopEO(const FermionField &in, FermionField &out,int dag)=0; + virtual void DhopDir(const FermionField &in, FermionField &out,int dir,int disp)=0; // implemented by WilsonFermion and WilsonFermion5D + + // force terms; five routines; default to Dhop on diagonal + virtual void MDeriv (LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag){DhopDeriv(mat,U,V,dag);}; + virtual void MoeDeriv(LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag){DhopDerivOE(mat,U,V,dag);}; + virtual void MeoDeriv(LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag){DhopDerivEO(mat,U,V,dag);}; + virtual void MooDeriv(LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag){mat=zero;}; + virtual void MeeDeriv(LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag){mat=zero;}; + + virtual void DhopDeriv (LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0; + virtual void DhopDerivEO(LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0; + virtual void DhopDerivOE(LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0; + + + virtual void Mdiag (const FermionField &in, FermionField &out) { Mooee(in,out);}; // Same as Mooee applied to both CB's + virtual void Mdir (const FermionField &in, FermionField &out,int dir,int disp)=0; // case by case Wilson, Clover, Cayley, ContFrac, PartFrac - virtual void Mdiag(const FermionField &in, FermionField &out) { Mooee(in,out);}; // Same as Mooee applied to both CB's - virtual void Mdir (const FermionField &in, FermionField &out,int dir,int disp)=0; // case by case Wilson, Clover, Cayley, ContFrac, PartFrac - virtual void DhopDir(const FermionField &in, FermionField &out,int dir,int disp)=0; // implemented by WilsonFermion and WilsonFermion5D }; diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index f5ff7cb0..afd7edd7 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -98,7 +98,9 @@ void WilsonFermion::Mdir (const LatticeFermion &in, LatticeFermion &out,int dir, DhopDir(in,out,dir,disp); } void WilsonFermion::DhopDir(const LatticeFermion &in, LatticeFermion &out,int dir,int disp){ + WilsonCompressor compressor(DaggerNo); + Stencil.HaloExchange(in,comm_buf,compressor); assert( (disp==1)||(disp==-1) ); @@ -109,9 +111,22 @@ void WilsonFermion::DhopDir(const LatticeFermion &in, LatticeFermion &out,int di PARALLEL_FOR_LOOP for(int sss=0;sssoSites();sss++){ - DiracOptDhopDir(Stencil,Umu,comm_buf,sss,sss,in,out,dirdisp); + DiracOptDhopDir(Stencil,Umu,comm_buf,sss,sss,in,out,dirdisp,dirdisp); } + +}; +void WilsonFermion::DhopDirDisp(const LatticeFermion &in, LatticeFermion &out,int dirdisp,int gamma,int dag) +{ + WilsonCompressor compressor(dag); + + Stencil.HaloExchange(in,comm_buf,compressor); + +PARALLEL_FOR_LOOP + for(int sss=0;sssoSites();sss++){ + DiracOptDhopDir(Stencil,Umu,comm_buf,sss,sss,in,out,dirdisp,gamma); + } + }; void WilsonFermion::DhopInternal(CartesianStencil & st,LatticeDoubledGaugeField & U, @@ -177,6 +192,77 @@ void WilsonFermion::Dhop(const LatticeFermion &in, LatticeFermion &out,int dag) DhopInternal(Stencil,Umu,in,out,dag); } +void WilsonFermion::DerivInternal(CartesianStencil & st,LatticeDoubledGaugeField & U, + LatticeGaugeField &mat,const LatticeFermion &A,const LatticeFermion &B,int dag) +{ + assert((dag==DaggerNo) ||(dag==DaggerYes)); + + WilsonCompressor compressor(dag); + + LatticeColourMatrix tmp(B._grid); + LatticeFermion Btilde(B._grid); + + st.HaloExchange(B,comm_buf,compressor); + + for(int mu=0;mu(1-g) if dag + //////////////////////////////////////////////////////////////////////// + int gamma = mu; + if ( dag ) gamma+= Nd; + + //////////////////////// + // Call the single hop + //////////////////////// +PARALLEL_FOR_LOOP + for(int sss=0;sssoSites();sss++){ + DiracOptDhopDir(st,U,comm_buf,sss,sss,B,Btilde,mu,gamma); + } + + ////////////////////////////////////////////////// + // spin trace outer product + ////////////////////////////////////////////////// + tmp = TraceIndex(outerProduct(Btilde,A)); + PokeIndex(mat,tmp,mu); + + } +} +void WilsonFermion::DhopDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag) +{ + conformable(U._grid,_grid); + conformable(U._grid,V._grid); + conformable(U._grid,mat._grid); + + mat.checkerboard = U.checkerboard; + + DerivInternal(Stencil,Umu,mat,U,V,dag); +} +void WilsonFermion::DhopDerivOE(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag) +{ + conformable(U._grid,_cbgrid); + conformable(U._grid,V._grid); + conformable(U._grid,mat._grid); + + assert(V.checkerboard==Even); + assert(U.checkerboard==Odd); + mat.checkerboard = Odd; + + DerivInternal(StencilEven,UmuOdd,mat,U,V,dag); +} +void WilsonFermion::DhopDerivEO(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag) +{ + conformable(U._grid,_cbgrid); + conformable(U._grid,V._grid); + conformable(U._grid,mat._grid); + + assert(V.checkerboard==Odd); + assert(U.checkerboard==Even); + mat.checkerboard = Even; + + DerivInternal(StencilOdd,UmuEven,mat,U,V,dag); +} + }} diff --git a/lib/qcd/action/fermion/WilsonFermion.h b/lib/qcd/action/fermion/WilsonFermion.h index 60726fdd..5422ca52 100644 --- a/lib/qcd/action/fermion/WilsonFermion.h +++ b/lib/qcd/action/fermion/WilsonFermion.h @@ -24,11 +24,95 @@ namespace Grid { // half checkerboard operaions void Meooe (const LatticeFermion &in, LatticeFermion &out); void MeooeDag (const LatticeFermion &in, LatticeFermion &out); + virtual void Mooee (const LatticeFermion &in, LatticeFermion &out); // remain virtual so we virtual void MooeeDag (const LatticeFermion &in, LatticeFermion &out); // can derive Clover virtual void MooeeInv (const LatticeFermion &in, LatticeFermion &out); // from Wilson base virtual void MooeeInvDag (const LatticeFermion &in, LatticeFermion &out); + //////////////////////// + // + // Force term: d/dtau S = 0 + // + // It is simplest to consider the two flavour force term + // + // S[U,phi] = phidag (MdagM)^-1 phi + // + // But simplify even this to + // + // S[U,phi] = phidag MdagM phi + // + // (other options exist depending on nature of action fragment.) + // + // Require momentum be traceless anti-hermitian to move within group manifold [ P = i P^a T^a ] + // + // Define the HMC hamiltonian + // + // H = 1/2 Tr P^2 + S(U,phi) + // + // . + // U = P U (lorentz & color indices multiplied) + // + // Hence + // + // .c c c c + // U = U P = - U P (c == dagger) + // + // So, taking some liberty with implicit indices + // . . .c c + // dH/dt = 0 = Tr P P +Tr[ U dS/dU + U dS/dU ] + // + // . c c + // = Tr P P + i Tr[ P U dS/dU - U P dS/dU ] + // + // . c c + // = Tr P (P + i ( U dS/dU - P dS/dU U ] + // + // . c c + // => P = -i [ U dS/dU - dS/dU U ] generates HMC EoM + // + // Simple case work this out using S = phi^dag MdagM phi for wilson: + // c c + // dSdt = dU_xdt dSdUx + dUxdt dSdUx + // + // = Tr i P U_x [ (\phi^\dag)_x (1+g) (M \phi)_x+\mu +(\phi^\dag M^\dag)_x (1-g) \phi_{x+\mu} ] + // c + // - i U_x P [ (\phi^\dag)_x+mu (1-g) (M \phi)_x +(\phi^\dag M^\dag)_(x+\mu) (1+g) \phi_{x} ] + // + // = i [(\phi^\dag)_x ]_j P_jk [U_x(1+g) (M \phi)_x+\mu]_k (1) + // + i [(\phi^\dagM^\dag)_x]_j P_jk [U_x(1-g) (\phi)_x+\mu]_k (2) + // - i [(\phi^\dag)_x+mu (1-g) U^dag_x]_j P_jk [(M \phi)_xk (3) + // - i [(\phi^\dagM^\dag)_x+mu (1+g) U^dag_x]_j P_jk [ \phi]_xk (4) + // + // Observe that (1)* = (4) + // (2)* = (3) + // + // Write as . + // P_{kj} = - i ( [U_x(1+g) (M \phi)_x+\mu] (x) [(\phi^\dag)_x] + [U_x(1-g) (\phi)_x+\mu] (x) [(\phi^\dagM^\dag)_x] - h.c ) + // + // where (x) denotes outer product in colour and spins are traced. + // + // Need only evaluate (1) and (2) [Chroma] or (2) and (4) [IroIro] and take the + // traceless anti hermitian part (of term in brackets w/o the "i") + // + // Generalisation to S=phi^dag (MdagM)^{-1} phi is simple: + // + // For more complicated DWF etc... apply product rule in differentiation + // + //////////////////////// + void DhopDeriv (LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag); + void DhopDerivEO(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag); + void DhopDerivOE(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag); + + // Extra support internal + void DerivInternal(CartesianStencil & st, + LatticeDoubledGaugeField & U, + LatticeGaugeField &mat, + const LatticeFermion &A, + const LatticeFermion &B, + int dag); + + // non-hermitian hopping term; half cb or both void Dhop (const LatticeFermion &in, LatticeFermion &out,int dag); void DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag); @@ -37,6 +121,7 @@ namespace Grid { // Multigrid assistance void Mdir (const LatticeFermion &in, LatticeFermion &out,int dir,int disp); void DhopDir(const LatticeFermion &in, LatticeFermion &out,int dir,int disp); + void DhopDirDisp(const LatticeFermion &in, LatticeFermion &out,int dirdisp,int gamma,int dag); /////////////////////////////////////////////////////////////// // Extra methods added by derived @@ -59,7 +144,8 @@ namespace Grid { static int HandOptDslash; // these are a temporary hack static int MortonOrder; - protected: + // protected: + public: RealD mass; diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 43829ea7..523934c4 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -105,14 +105,105 @@ PARALLEL_FOR_LOOP for(int s=0;s(B,comm_buf,compressor); + + for(int mu=0;muoSites();sss++){ + for(int s=0;s(outerProduct(Btilde,A)); // ordering here + PokeIndex(mat,tmp,mu); + + } +} + +void WilsonFermion5D::DhopDeriv( LatticeGaugeField &mat, + const LatticeFermion &A, + const LatticeFermion &B, + int dag) +{ + conformable(A._grid,FermionGrid()); + conformable(A._grid,B._grid); + conformable(GaugeGrid(),mat._grid); + + mat.checkerboard = A.checkerboard; + + DerivInternal(Stencil,Umu,mat,A,B,dag); +} + +void WilsonFermion5D::DhopDerivEO(LatticeGaugeField &mat, + const LatticeFermion &A, + const LatticeFermion &B, + int dag) +{ + conformable(A._grid,FermionRedBlackGrid()); + conformable(GaugeRedBlackGrid(),mat._grid); + conformable(A._grid,B._grid); + + assert(B.checkerboard==Odd); + assert(A.checkerboard==Even); + mat.checkerboard = Even; + + DerivInternal(StencilOdd,UmuEven,mat,A,B,dag); +} + +void WilsonFermion5D::DhopDerivOE(LatticeGaugeField &mat, + const LatticeFermion &A, + const LatticeFermion &B, + int dag) +{ + conformable(A._grid,FermionRedBlackGrid()); + conformable(GaugeRedBlackGrid(),mat._grid); + conformable(A._grid,B._grid); + + assert(B.checkerboard==Even); + assert(A.checkerboard==Odd); + mat.checkerboard = Odd; + + DerivInternal(StencilEven,UmuOdd,mat,A,B,dag); +} + void WilsonFermion5D::DhopInternal(CartesianStencil & st, LebesgueOrder &lo, LatticeDoubledGaugeField & U, - const LatticeFermion &in, LatticeFermion &out,int dag) + const LatticeFermion &in, LatticeFermion &out,int dag) { // assert((dag==DaggerNo) ||(dag==DaggerYes)); diff --git a/lib/qcd/action/fermion/WilsonFermion5D.h b/lib/qcd/action/fermion/WilsonFermion5D.h index e345d2a4..7ff5d1c3 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.h +++ b/lib/qcd/action/fermion/WilsonFermion5D.h @@ -44,12 +44,18 @@ namespace Grid { // half checkerboard operations; leave unimplemented as abstract for now virtual void Meooe (const LatticeFermion &in, LatticeFermion &out){assert(0);}; - virtual void MeooeDag (const LatticeFermion &in, LatticeFermion &out){assert(0);}; virtual void Mooee (const LatticeFermion &in, LatticeFermion &out){assert(0);}; - virtual void MooeeDag (const LatticeFermion &in, LatticeFermion &out){assert(0);}; virtual void MooeeInv (const LatticeFermion &in, LatticeFermion &out){assert(0);}; + + virtual void MeooeDag (const LatticeFermion &in, LatticeFermion &out){assert(0);}; + virtual void MooeeDag (const LatticeFermion &in, LatticeFermion &out){assert(0);}; virtual void MooeeInvDag (const LatticeFermion &in, LatticeFermion &out){assert(0);}; + // These can be overridden by fancy 5d chiral actions + virtual void DhopDeriv (LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag); + virtual void DhopDerivEO(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag); + virtual void DhopDerivOE(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag); + // Implement hopping term non-hermitian hopping term; half cb or both // Implement s-diagonal DW void DW (const LatticeFermion &in, LatticeFermion &out,int dag); @@ -64,6 +70,14 @@ namespace Grid { /////////////////////////////////////////////////////////////// // New methods added /////////////////////////////////////////////////////////////// + + void DerivInternal(CartesianStencil & st, + LatticeDoubledGaugeField & U, + LatticeGaugeField &mat, + const LatticeFermion &A, + const LatticeFermion &B, + int dag); + void DhopInternal(CartesianStencil & st, LebesgueOrder &lo, LatticeDoubledGaugeField &U, diff --git a/lib/qcd/action/fermion/WilsonKernels.cc b/lib/qcd/action/fermion/WilsonKernels.cc index bdc18eee..43ccff63 100644 --- a/lib/qcd/action/fermion/WilsonKernels.cc +++ b/lib/qcd/action/fermion/WilsonKernels.cc @@ -294,8 +294,8 @@ void DiracOptDhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U, } void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U, - std::vector > &buf, - int sF,int sU,const LatticeFermion &in, LatticeFermion &out,int dirdisp) + std::vector > &buf, + int sF,int sU,const LatticeFermion &in, LatticeFermion &out,int dir,int gamma) { vHalfSpinColourVector tmp; vHalfSpinColourVector chi; @@ -304,13 +304,13 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U, int offset,local,perm, ptype; int ss=sF; - offset = st._offsets [dirdisp][ss]; - local = st._is_local[dirdisp][ss]; - perm = st._permute[dirdisp][ss]; - ptype = st._permute_type[dirdisp]; + offset = st._offsets [dir][ss]; + local = st._is_local[dir][ss]; + perm = st._permute[dir][ss]; + ptype = st._permute_type[dir]; // Xp - if(dirdisp==Xp){ + if(gamma==Xp){ if ( local && perm ) { spProjXp(tmp,in._odata[offset]); permute(chi,tmp,ptype); @@ -319,12 +319,12 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[sU](Xp),&chi()); + mult(&Uchi(),&U._odata[sU](dir),&chi()); spReconXp(result,Uchi); } // Yp - if ( dirdisp==Yp ){ + if ( gamma==Yp ){ if ( local && perm ) { spProjYp(tmp,in._odata[offset]); permute(chi,tmp,ptype); @@ -333,12 +333,12 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[sU](Yp),&chi()); + mult(&Uchi(),&U._odata[sU](dir),&chi()); spReconYp(result,Uchi); } // Zp - if ( dirdisp ==Zp ){ + if ( gamma ==Zp ){ if ( local && perm ) { spProjZp(tmp,in._odata[offset]); permute(chi,tmp,ptype); @@ -347,12 +347,12 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[sU](Zp),&chi()); + mult(&Uchi(),&U._odata[sU](dir),&chi()); spReconZp(result,Uchi); } // Tp - if ( dirdisp ==Tp ){ + if ( gamma ==Tp ){ if ( local && perm ) { spProjTp(tmp,in._odata[offset]); permute(chi,tmp,ptype); @@ -361,12 +361,12 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[sU](Tp),&chi()); + mult(&Uchi(),&U._odata[sU](dir),&chi()); spReconTp(result,Uchi); } // Xm - if ( dirdisp==Xm ){ + if ( gamma==Xm ){ if ( local && perm ) { spProjXm(tmp,in._odata[offset]); permute(chi,tmp,ptype); @@ -375,12 +375,12 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[sU](Xm),&chi()); + mult(&Uchi(),&U._odata[sU](dir),&chi()); spReconXm(result,Uchi); } // Ym - if ( dirdisp == Ym ){ + if ( gamma == Ym ){ if ( local && perm ) { spProjYm(tmp,in._odata[offset]); permute(chi,tmp,ptype); @@ -389,12 +389,12 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[sU](Ym),&chi()); + mult(&Uchi(),&U._odata[sU](dir),&chi()); spReconYm(result,Uchi); } // Zm - if ( dirdisp == Zm ){ + if ( gamma == Zm ){ if ( local && perm ) { spProjZm(tmp,in._odata[offset]); permute(chi,tmp,ptype); @@ -403,12 +403,12 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[sU](Zm),&chi()); + mult(&Uchi(),&U._odata[sU](dir),&chi()); spReconZm(result,Uchi); } // Tm - if ( dirdisp==Tm ) { + if ( gamma==Tm ) { if ( local && perm ) { spProjTm(tmp,in._odata[offset]); permute(chi,tmp,ptype); @@ -417,7 +417,7 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[sU](Tm),&chi()); + mult(&Uchi(),&U._odata[sU](dir),&chi()); spReconTm(result,Uchi); } diff --git a/lib/qcd/action/fermion/WilsonKernels.h b/lib/qcd/action/fermion/WilsonKernels.h index ec807af5..cb0e07bb 100644 --- a/lib/qcd/action/fermion/WilsonKernels.h +++ b/lib/qcd/action/fermion/WilsonKernels.h @@ -22,7 +22,7 @@ namespace Grid { int sF,int sU,const LatticeFermion &in, LatticeFermion &out); void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U, std::vector > &buf, - int sF,int sU,const LatticeFermion &in, LatticeFermion &out,int dirdisp); + int sF,int sU,const LatticeFermion &in, LatticeFermion &out,int dirdisp,int gamma); // }; diff --git a/lib/qcd/hmc/HMC.h b/lib/qcd/hmc/HMC.h index e3fe57ed..21c6bb0d 100644 --- a/lib/qcd/hmc/HMC.h +++ b/lib/qcd/hmc/HMC.h @@ -68,14 +68,12 @@ namespace Grid{ Integrator& MolDyn): Params(Pms),MD(MolDyn){ //FIXME - // initialize RNGs also with seed sRNG.SeedRandomDevice(); } ~HybridMonteCarlo(){}; - void evolve(LatticeLorentzColourMatrix& Uin){ Real DeltaH; diff --git a/lib/simd/Grid_vector_types.h b/lib/simd/Grid_vector_types.h index 034d5314..71472407 100644 --- a/lib/simd/Grid_vector_types.h +++ b/lib/simd/Grid_vector_types.h @@ -414,7 +414,7 @@ namespace Grid { template inline Grid_simd< S, V> outerProduct(const Grid_simd< S, V> &l, const Grid_simd< S, V> & r) { - return l*r; + return l*conjugate(r); } template diff --git a/lib/tensors/Tensor_outer.h b/lib/tensors/Tensor_outer.h index 816f6ad8..d0150709 100644 --- a/lib/tensors/Tensor_outer.h +++ b/lib/tensors/Tensor_outer.h @@ -28,11 +28,13 @@ auto outerProduct (const iScalar& lhs,const iScalar& rhs) -> iScalar &arg) -> iScalar + inline auto trace(const iVector &arg) -> iVector +{ + iVector ret; + for(int i=0;i U(4,UGrid); for(int mu=0;mu + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +#define parallel_for PARALLEL_FOR_LOOP for + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + + GridCartesian Grid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); + + int threads = GridThread::GetThreads(); + std::cout< seeds({1,2,3,4}); + + GridParallelRNG pRNG(&Grid); + pRNG.SeedRandomDevice(); + + LatticeFermion phi (&Grid); gaussian(pRNG,phi); + LatticeFermion Mphi (&Grid); + LatticeFermion MphiPrime (&Grid); + + LatticeGaugeField U(&Grid); + + SU3::HotConfiguration(pRNG,U); + // SU3::ColdConfiguration(pRNG,U); + + //////////////////////////////////// + // Unmodified matrix element + //////////////////////////////////// + RealD mass=-4.0; //kills the diagonal term + WilsonFermion Dw (U, Grid,RBGrid,mass); + Dw.M (phi,Mphi); + + ComplexD S = innerProduct(Mphi,Mphi); // pdag MdagM p + + // get the deriv of phidag MdagM phi with respect to "U" + LatticeGaugeField UdSdU(&Grid); + LatticeGaugeField tmp(&Grid); + + Dw.MDeriv(tmp , Mphi, phi,DaggerNo ); UdSdU=tmp; + Dw.MDeriv(tmp , phi, Mphi,DaggerYes ); UdSdU=UdSdU+tmp; + + + LatticeFermion Ftmp (&Grid); + + //////////////////////////////////// + // Modify the gauge field a little + //////////////////////////////////// + RealD dt = 1.0e-6; + + LatticeColourMatrix mommu(&Grid); + LatticeGaugeField mom(&Grid); + LatticeGaugeField Uprime(&Grid); + + for(int mu=0;mu(mom,mommu,mu); + parallel_for(auto i=mom.begin();i + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +#define parallel_for PARALLEL_FOR_LOOP for + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + + GridCartesian Grid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); + + int threads = GridThread::GetThreads(); + std::cout< seeds({1,2,3,4}); + + GridParallelRNG pRNG(&Grid); + pRNG.SeedRandomDevice(); + + LatticeFermion phi (&Grid); gaussian(pRNG,phi); + LatticeFermion Mphi (&Grid); + LatticeFermion Mdagphi (&Grid); + LatticeFermion MphiPrime (&Grid); + LatticeFermion MdagphiPrime (&Grid); + LatticeFermion dMphi (&Grid); + + LatticeGaugeField U(&Grid); + + + SU3::HotConfiguration(pRNG,U); + // SU3::ColdConfiguration(pRNG,U); + + //////////////////////////////////// + // Unmodified matrix element + //////////////////////////////////// + RealD mass=-4.0; //kills the diagonal term + WilsonFermion Dw (U, Grid,RBGrid,mass); + Dw.M (phi,Mphi); + Dw.Mdag(phi,Mdagphi); + + ComplexD S = innerProduct(Mphi,Mphi); // pdag MdagM p + ComplexD Sdag = innerProduct(Mdagphi,Mdagphi); // pdag MMdag p + + // get the deriv of phidag MdagM phi with respect to "U" + LatticeGaugeField UdSdU(&Grid); + LatticeGaugeField UdSdUdag(&Grid); + LatticeGaugeField tmp(&Grid); + + Dw.MDeriv(tmp , Mphi, phi,DaggerNo ); UdSdU=tmp; + + Dw.MDeriv(tmp , Mdagphi, phi,DaggerYes ); UdSdUdag=tmp; + + + LatticeFermion dMdagphi (&Grid); dMdagphi=zero; + LatticeFermion Ftmp (&Grid); + + + // Dw.MDeriv(UdSdU,Mdagphi, phi,DaggerYes );// UdSdU =UdSdU +tmp; + + //////////////////////////////////// + // Modify the gauge field a little in one dir + //////////////////////////////////// + RealD dt = 1.0e-3; + + LatticeColourMatrix mommu(&Grid); + LatticeGaugeField mom(&Grid); + LatticeGaugeField Uprime(&Grid); + + for(int mu=0;mu(mom,mommu,mu); + + parallel_for(auto i=mom.begin();i + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +#define parallel_for PARALLEL_FOR_LOOP for + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + + GridCartesian Grid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); + + int threads = GridThread::GetThreads(); + std::cout< seeds({1,2,3,4}); + + GridParallelRNG pRNG(&Grid); + pRNG.SeedRandomDevice(); + + LatticeFermion phi (&Grid); gaussian(pRNG,phi); + LatticeFermion Mphi (&Grid); + LatticeFermion MphiPrime (&Grid); + LatticeFermion dMphi (&Grid); + + LatticeGaugeField U(&Grid); + + + SU3::HotConfiguration(pRNG,U); + + //////////////////////////////////// + // Unmodified matrix element + //////////////////////////////////// + RealD mass=-4.0; //kills the diagonal term + WilsonFermion Dw (U, Grid,RBGrid,mass); + Dw.M(phi,Mphi); + + ComplexD S = innerProduct(phi,Mphi); + + // get the deriv + LatticeGaugeField UdSdU(&Grid); + Dw.MDeriv(UdSdU,phi, phi,DaggerNo ); + + + //////////////////////////////////// + // Modify the gauge field a little in one dir + //////////////////////////////////// + RealD dt = 1.0e-3; + Complex Complex_i(0,1); + + LatticeColourMatrix Umu(&Grid); + LatticeColourMatrix Umu_save(&Grid); + LatticeColourMatrix dU (&Grid); + LatticeColourMatrix mom(&Grid); + SU3::GaussianLieAlgebraMatrix(pRNG, mom); // Traceless antihermitian momentum; gaussian in lie alg + + + // check mom is as i expect + LatticeColourMatrix tmpmom(&Grid); + tmpmom = mom+adj(mom); + std::cout << GridLogMessage << "mom anti-herm check "<< norm2(tmpmom)<(U,mu); + Umu_save=Umu; + dU = mom * Umu * dt; + Umu= Umu+dU; + PokeIndex(Dw.Umu,Umu,mu); + + Dw.M(phi,MphiPrime); + + ComplexD Sprime = innerProduct(phi,MphiPrime); + + std::cout << GridLogMessage << " S "<(Dw.Umu,dU,mu); + Dw.M(phi,dMphi); + + + ComplexD deltaS = innerProduct(phi,dMphi); + std::cout << GridLogMessage << "deltaS "<(Dw.Umu,Umu_save,mu); + Dw.Mdir(phi,dMphi,mu,1); + dMphi = dt*mom*dMphi; + + deltaS = innerProduct(phi,dMphi); + std::cout << GridLogMessage << "deltaS from inner prod of mom* M[u] "<