diff --git a/configure b/configure index caf00467..2f14088e 100755 --- a/configure +++ b/configure @@ -626,10 +626,14 @@ ac_subst_vars='am__EXEEXT_FALSE am__EXEEXT_TRUE LTLIBOBJS LIBOBJS +BUILD_CHROMA_REGRESSION_FALSE +BUILD_CHROMA_REGRESSION_TRUE BUILD_COMMS_NONE_FALSE BUILD_COMMS_NONE_TRUE BUILD_COMMS_MPI_FALSE BUILD_COMMS_MPI_TRUE +BUILD_ZMM_FALSE +BUILD_ZMM_TRUE EGREP GREP CXXCPP @@ -745,6 +749,7 @@ enable_openmp enable_simd enable_precision enable_comms +enable_chroma ' ac_precious_vars='build_alias host_alias @@ -1390,6 +1395,7 @@ Optional Features: --enable-precision=single|double Select default word size of Real --enable-comms=none|mpi Select communications + --enable-chroma Expect chroma compiled under c++11 Some influential environment variables: CXX C++ compiler command @@ -6392,6 +6398,8 @@ fi supported=no +ac_ZMM=no; + case ${ac_SIMD} in SSE4) echo Configuring for SSE4 @@ -6443,6 +6451,7 @@ $as_echo "$as_me: WARNING: Your processor does not support AVX2 instructions" >& $as_echo "#define AVX512 1" >>confdefs.h supported="cross compilation" + ac_ZMM=yes; ;; IMCI) echo Configuring for IMCI @@ -6450,6 +6459,7 @@ $as_echo "#define AVX512 1" >>confdefs.h $as_echo "#define IMCI 1" >>confdefs.h supported="cross compilation" + ac_ZMM=yes; ;; NEONv8) echo Configuring for experimental ARMv8a support @@ -6469,6 +6479,24 @@ $as_echo "#define EMPTY_SIMD 1" >>confdefs.h ;; esac +case ${ac_ZMM} in +yes) + echo Enabling ZMM source code +;; +no) + echo Disabling ZMM source code +;; +esac + + if test "X${ac_ZMM}X" == "XyesX" ; then + BUILD_ZMM_TRUE= + BUILD_ZMM_FALSE='#' +else + BUILD_ZMM_TRUE='#' + BUILD_ZMM_FALSE= +fi + + # Check whether --enable-precision was given. if test "${enable_precision+set}" = set; then : enableval=$enable_precision; ac_PRECISION=${enable_precision} @@ -6534,6 +6562,34 @@ else fi +# Check whether --enable-chroma was given. +if test "${enable_chroma+set}" = set; then : + enableval=$enable_chroma; ac_CHROMA=yes +else + ac_CHROMA=no +fi + + +case ${ac_CHROMA} in + yes) + echo Enabling tests regressing to Chroma + ;; + no) + echo Disabling tests regressing to Chroma + ;; + *) + as_fn_error $? "${ac_CHROMA} unsupported --enable-chroma option" "$LINENO" 5; + ;; +esac + + if test "X${ac_CHROMA}X" == "XyesX" ; then + BUILD_CHROMA_REGRESSION_TRUE= + BUILD_CHROMA_REGRESSION_FALSE='#' +else + BUILD_CHROMA_REGRESSION_TRUE='#' + BUILD_CHROMA_REGRESSION_FALSE= +fi + ################################################################### # Checks for doxygen support @@ -6557,6 +6613,8 @@ ac_config_files="$ac_config_files lib/Makefile" ac_config_files="$ac_config_files tests/Makefile" +ac_config_files="$ac_config_files tests/qdpxx/Makefile" + ac_config_files="$ac_config_files benchmarks/Makefile" cat >confcache <<\_ACEOF @@ -6696,6 +6754,10 @@ if test -z "${am__fastdepCC_TRUE}" && test -z "${am__fastdepCC_FALSE}"; then as_fn_error $? "conditional \"am__fastdepCC\" was never defined. Usually this means the macro was only invoked conditionally." "$LINENO" 5 fi +if test -z "${BUILD_ZMM_TRUE}" && test -z "${BUILD_ZMM_FALSE}"; then + as_fn_error $? "conditional \"BUILD_ZMM\" was never defined. +Usually this means the macro was only invoked conditionally." "$LINENO" 5 +fi if test -z "${BUILD_COMMS_MPI_TRUE}" && test -z "${BUILD_COMMS_MPI_FALSE}"; then as_fn_error $? "conditional \"BUILD_COMMS_MPI\" was never defined. Usually this means the macro was only invoked conditionally." "$LINENO" 5 @@ -6704,6 +6766,10 @@ if test -z "${BUILD_COMMS_NONE_TRUE}" && test -z "${BUILD_COMMS_NONE_FALSE}"; th as_fn_error $? "conditional \"BUILD_COMMS_NONE\" was never defined. Usually this means the macro was only invoked conditionally." "$LINENO" 5 fi +if test -z "${BUILD_CHROMA_REGRESSION_TRUE}" && test -z "${BUILD_CHROMA_REGRESSION_FALSE}"; then + as_fn_error $? "conditional \"BUILD_CHROMA_REGRESSION\" was never defined. +Usually this means the macro was only invoked conditionally." "$LINENO" 5 +fi : "${CONFIG_STATUS=./config.status}" ac_write_fail=0 @@ -7301,6 +7367,7 @@ do "Makefile") CONFIG_FILES="$CONFIG_FILES Makefile" ;; "lib/Makefile") CONFIG_FILES="$CONFIG_FILES lib/Makefile" ;; "tests/Makefile") CONFIG_FILES="$CONFIG_FILES tests/Makefile" ;; + "tests/qdpxx/Makefile") CONFIG_FILES="$CONFIG_FILES tests/qdpxx/Makefile" ;; "benchmarks/Makefile") CONFIG_FILES="$CONFIG_FILES benchmarks/Makefile" ;; *) as_fn_error $? "invalid argument: \`$ac_config_target'" "$LINENO" 5;; diff --git a/configure.ac b/configure.ac index c8fd8b1a..094a70bd 100644 --- a/configure.ac +++ b/configure.ac @@ -72,6 +72,8 @@ AC_ARG_ENABLE([simd],[AC_HELP_STRING([--enable-simd=SSE4|AVX|AVXFMA4|AVX2|AVX512 supported=no +ac_ZMM=no; + case ${ac_SIMD} in SSE4) echo Configuring for SSE4 @@ -113,11 +115,13 @@ case ${ac_SIMD} in echo Configuring for AVX512 AC_DEFINE([AVX512],[1],[AVX512 Intrinsics for Knights Landing] ) supported="cross compilation" + ac_ZMM=yes; ;; IMCI) echo Configuring for IMCI AC_DEFINE([IMCI],[1],[IMCI Intrinsics for Knights Corner] ) supported="cross compilation" + ac_ZMM=yes; ;; NEONv8) echo Configuring for experimental ARMv8a support @@ -133,6 +137,17 @@ case ${ac_SIMD} in ;; esac +case ${ac_ZMM} in +yes) + echo Enabling ZMM source code +;; +no) + echo Disabling ZMM source code +;; +esac + +AM_CONDITIONAL(BUILD_ZMM,[ test "X${ac_ZMM}X" == "XyesX" ]) + AC_ARG_ENABLE([precision],[AC_HELP_STRING([--enable-precision=single|double],[Select default word size of Real])],[ac_PRECISION=${enable_precision}],[ac_PRECISION=double]) case ${ac_PRECISION} in single) diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h index 236f0ba0..242e8b73 100644 --- a/lib/qcd/action/Actions.h +++ b/lib/qcd/action/Actions.h @@ -16,8 +16,6 @@ #include #include - - //////////////////////////////////////////// // Utility functions //////////////////////////////////////////// @@ -26,7 +24,6 @@ #include #include //used by all wilson type fermions - //////////////////////////////////////////// // Gauge Actions //////////////////////////////////////////// @@ -55,9 +52,9 @@ typedef WilsonGaugeAction WilsonGaugeActionD; #define FermOpTemplateInstantiate(A) \ template class A; \ - template class A; - // template class A; \ - // template class A; + template class A; \ + template class A; \ + template class A; //////////////////////////////////////////// // Fermion operators / actions diff --git a/lib/qcd/action/fermion/FermionOperatorImpl.h b/lib/qcd/action/fermion/FermionOperatorImpl.h index 580f5e30..21187870 100644 --- a/lib/qcd/action/fermion/FermionOperatorImpl.h +++ b/lib/qcd/action/fermion/FermionOperatorImpl.h @@ -335,13 +335,33 @@ PARALLEL_FOR_LOOP } inline void InsertForce4D(GaugeField &mat, FermionField &Btilde, FermionField &A,int mu){ - assert(0); - // Fixme + + // DhopDir provides U or Uconj depending on coor/flavour. + GaugeLinkField link(mat._grid); + // use lorentz for flavour as hack. + auto tmp = TraceIndex(outerProduct(Btilde,A)); +PARALLEL_FOR_LOOP + for(auto ss=tmp.begin();ss(mat,link,mu); return; } - inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &A,int mu){ - assert(0); - // Fixme + inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField Ã,int mu){ + + int Ls=Btilde._grid->_fdimensions[0]; + + GaugeLinkField tmp(mat._grid); + tmp = zero; +PARALLEL_FOR_LOOP + for(int ss=0;ssoSites();ss++){ + for(int s=0;s(outerProduct(Btilde[sF],Atilde[sF])); + tmp[ss]() = tmp[ss]()+ ttmp(0,0) + conjugate(ttmp(1,1)); + } + } + PokeIndex(mat,tmp,mu); return; } }; diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index 5e37238c..4189c608 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -131,7 +131,7 @@ namespace QCD { // Flip gamma (1+g)<->(1-g) if dag //////////////////////////////////////////////////////////////////////// int gamma = mu; - if ( dag ) gamma+= Nd; + if ( !dag ) gamma+= Nd; //////////////////////// // Call the single hop @@ -227,13 +227,15 @@ PARALLEL_FOR_LOOP DhopDir(in,out,dir,disp); } + template void WilsonFermion::DhopDir(const FermionField &in, FermionField &out,int dir,int disp){ int skip = (disp==1) ? 0 : 1; - int dirdisp = dir+skip*4; + int dirdisp = dir+skip*4; + int gamma = dir+(1-skip)*4; - DhopDirDisp(in,out,dirdisp,dirdisp,DaggerNo); + DhopDirDisp(in,out,dirdisp,gamma,DaggerNo); }; diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 410aa629..54e6691e 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -94,6 +94,7 @@ void WilsonFermion5D::DhopDir(const FermionField &in, FermionField &out,in int skip = (disp==1) ? 0 : 1; int dirdisp = dir+skip*4; + int gamma = dir+(1-skip)*4; assert(dirdisp<=7); assert(dirdisp>=0); @@ -103,7 +104,7 @@ PARALLEL_FOR_LOOP for(int s=0;s::DerivInternal(StencilImpl & st, // Flip gamma if dag //////////////////////////////////////////////////////////////////////// int gamma = mu; - if ( dag ) gamma+= Nd; + if ( !dag ) gamma+= Nd; //////////////////////// // Call the single hop diff --git a/lib/qcd/action/fermion/WilsonKernels.h b/lib/qcd/action/fermion/WilsonKernels.h index 065ce678..80d36948 100644 --- a/lib/qcd/action/fermion/WilsonKernels.h +++ b/lib/qcd/action/fermion/WilsonKernels.h @@ -33,7 +33,6 @@ namespace Grid { std::vector > &buf, int sF,int sU,const FermionField &in, FermionField &out,uint64_t *); -#define HANDOPT void DiracOptHandDhopSite(StencilImpl &st,DoubledGaugeField &U, std::vector > &buf, int sF,int sU,const FermionField &in, FermionField &out); diff --git a/lib/qcd/action/fermion/WilsonKernelsHand.cc b/lib/qcd/action/fermion/WilsonKernelsHand.cc index c1d01f51..51a9f19a 100644 --- a/lib/qcd/action/fermion/WilsonKernelsHand.cc +++ b/lib/qcd/action/fermion/WilsonKernelsHand.cc @@ -280,7 +280,7 @@ namespace Grid { namespace QCD { -#ifdef HANDOPT + template void WilsonKernels::DiracOptHandDhopSiteDag(StencilImpl &st,DoubledGaugeField &U, std::vector > &buf, @@ -767,24 +767,43 @@ void WilsonKernels::DiracOptHandDhopSite(StencilImpl &st,DoubledGaugeFiel vstream(ref()(3)(2),result_32*(-0.5)); } } -#else -template -void WilsonKernels::DiracOptHandDhopSite(StencilImpl &st,DoubledGaugeField &U, + + //////////////////////////////////////////////// + // Specialise Gparity to simple implementation + //////////////////////////////////////////////// +template<> +void WilsonKernels::DiracOptHandDhopSite(StencilImpl &st,DoubledGaugeField &U, std::vector > &buf, - int ss,int sU,const FermionField &in, FermionField &out) + int sF,int sU,const FermionField &in, FermionField &out) { DiracOptDhopSite(st,U,buf,sF,sU,in,out); // will template override for Wilson Nc=3 } -template -void WilsonKernels::DiracOptHandDhopSiteDag(StencilImpl &st,DoubledGaugeField &U, +template<> +void WilsonKernels::DiracOptHandDhopSiteDag(StencilImpl &st,DoubledGaugeField &U, std::vector > &buf, - int ss,int sU,const FermionField &in, FermionField &out) + int sF,int sU,const FermionField &in, FermionField &out) { DiracOptDhopSiteDag(st,U,buf,sF,sU,in,out); // will template override for Wilson Nc=3 } -#endif +template<> +void WilsonKernels::DiracOptHandDhopSite(StencilImpl &st,DoubledGaugeField &U, + std::vector > &buf, + int sF,int sU,const FermionField &in, FermionField &out) +{ + DiracOptDhopSite(st,U,buf,sF,sU,in,out); // will template override for Wilson Nc=3 +} + +template<> +void WilsonKernels::DiracOptHandDhopSiteDag(StencilImpl &st,DoubledGaugeField &U, + std::vector > &buf, + int sF,int sU,const FermionField &in, FermionField &out) +{ + DiracOptDhopSiteDag(st,U,buf,sF,sU,in,out); // will template override for Wilson Nc=3 +} + + template void WilsonKernels::DiracOptHandDhopSite(StencilImpl &st,DoubledGaugeField &U, std::vector > &buf, @@ -799,4 +818,18 @@ template void WilsonKernels::DiracOptHandDhopSiteDag(StencilImpl &s std::vector > &buf, int ss,int sU,const FermionField &in, FermionField &out); + +template void WilsonKernels::DiracOptHandDhopSite(StencilImpl &st,DoubledGaugeField &U, + std::vector > &buf, + int ss,int sU,const FermionField &in, FermionField &out); +template void WilsonKernels::DiracOptHandDhopSite(StencilImpl &st,DoubledGaugeField &U, + std::vector > &buf, + int ss,int sU,const FermionField &in, FermionField &out); +template void WilsonKernels::DiracOptHandDhopSiteDag(StencilImpl &st,DoubledGaugeField &U, + std::vector > &buf, + int ss,int sU,const FermionField &in, FermionField &out); +template void WilsonKernels::DiracOptHandDhopSiteDag(StencilImpl &st,DoubledGaugeField &U, + std::vector > &buf, + int ss,int sU,const FermionField &in, FermionField &out); + }} diff --git a/scripts/filelist b/scripts/filelist index fd0e74ad..9d376e42 100755 --- a/scripts/filelist +++ b/scripts/filelist @@ -18,7 +18,7 @@ TESTS=`ls T*.cc` TESTLIST=`echo ${TESTS} | sed s/.cc//g ` echo > Make.inc -echo bin_PROGRAMS = ${TESTLIST} >> Make.inc +echo bin_PROGRAMS = ${TESTLIST} | sed s/Test_zmm//g >> Make.inc echo >> Make.inc for f in $TESTS diff --git a/tests/Make.inc b/tests/Make.inc index 4ce52d10..a302effb 100644 --- a/tests/Make.inc +++ b/tests/Make.inc @@ -1,5 +1,5 @@ -bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_cheby Test_contfrac_cg Test_contfrac_even_odd Test_contfrac_force Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_force Test_dwf_fpgcr Test_dwf_hdcr Test_dwf_lanczos Test_gamma Test_hmc_EODWFRatio Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonRatio Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_partfrac_force Test_quenched_update Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio Test_rng Test_rng_fixed Test_serialisation Test_simd Test_stencil Test_synthetic_lanczos Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi +bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_cheby Test_contfrac_cg Test_contfrac_even_odd Test_contfrac_force Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_force Test_dwf_fpgcr Test_dwf_gpforce Test_dwf_hdcr Test_dwf_lanczos Test_gamma Test_gparity Test_gpdwf_force Test_gpwilson_even_odd Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonRatio Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_partfrac_force Test_quenched_update Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio Test_rng Test_rng_fixed Test_serialisation Test_simd Test_stencil Test_synthetic_lanczos Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi Test_GaugeAction_SOURCES=Test_GaugeAction.cc @@ -82,6 +82,10 @@ Test_dwf_fpgcr_SOURCES=Test_dwf_fpgcr.cc Test_dwf_fpgcr_LDADD=-lGrid +Test_dwf_gpforce_SOURCES=Test_dwf_gpforce.cc +Test_dwf_gpforce_LDADD=-lGrid + + Test_dwf_hdcr_SOURCES=Test_dwf_hdcr.cc Test_dwf_hdcr_LDADD=-lGrid @@ -98,6 +102,10 @@ Test_gparity_SOURCES=Test_gparity.cc Test_gparity_LDADD=-lGrid +Test_gpdwf_force_SOURCES=Test_gpdwf_force.cc +Test_gpdwf_force_LDADD=-lGrid + + Test_gpwilson_even_odd_SOURCES=Test_gpwilson_even_odd.cc Test_gpwilson_even_odd_LDADD=-lGrid @@ -106,6 +114,10 @@ Test_hmc_EODWFRatio_SOURCES=Test_hmc_EODWFRatio.cc Test_hmc_EODWFRatio_LDADD=-lGrid +Test_hmc_EODWFRatio_Gparity_SOURCES=Test_hmc_EODWFRatio_Gparity.cc +Test_hmc_EODWFRatio_Gparity_LDADD=-lGrid + + Test_hmc_EOWilsonFermionGauge_SOURCES=Test_hmc_EOWilsonFermionGauge.cc Test_hmc_EOWilsonFermionGauge_LDADD=-lGrid @@ -225,10 +237,7 @@ Test_wilson_force_phiMdagMphi_LDADD=-lGrid Test_wilson_force_phiMphi_SOURCES=Test_wilson_force_phiMphi.cc Test_wilson_force_phiMphi_LDADD=-lGrid + Test_zmm_SOURCES=Test_zmm.cc Test_zmm_LDADD=-lGrid -Test_RectPlaq_SOURCES=Test_RectPlaq.cc -Test_RectPlaq_LDADD=-lGrid - - diff --git a/tests/Makefile.am b/tests/Makefile.am index 30d22d84..336cb6e7 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -8,5 +8,8 @@ endif AM_CXXFLAGS = -I$(top_srcdir)/lib AM_LDFLAGS = -L$(top_builddir)/lib +if BUILD_ZMM + bin_PROGRAMS=Test_zmm +endif include Make.inc diff --git a/tests/Test_dwf_gpforce.cc b/tests/Test_dwf_gpforce.cc new file mode 100644 index 00000000..c72f40d5 --- /dev/null +++ b/tests/Test_dwf_gpforce.cc @@ -0,0 +1,191 @@ +#include + +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(); + + const int Ls=8; + 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); + + typedef typename GparityDomainWallFermionR::FermionField FermionField; + + int threads = GridThread::GetThreads(); + std::cout< seeds({1,2,3,4}); + + GridParallelRNG RNG5(FGrid); RNG5.SeedRandomDevice(); + GridParallelRNG RNG4(UGrid); RNG4.SeedRandomDevice(); + + FermionField phi (FGrid); gaussian(RNG5,phi); + FermionField Mphi (FGrid); + FermionField MphiPrime (FGrid); + + LatticeGaugeField U(UGrid); + + SU3::HotConfiguration(RNG4,U); + // SU3::ColdConfiguration(pRNG,U); + + //////////////////////////////////// + // Unmodified matrix element + //////////////////////////////////// + RealD mass=0.2; //kills the diagonal term + RealD M5=1.8; + // const int nu = 3; + // std::vector twists(Nd,0); // twists[nu] = 1; + // GparityDomainWallFermionR::ImplParams params; params.twists = twists; + // GparityDomainWallFermionR Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,params); + + // DomainWallFermionR Dw (U, Grid,RBGrid,mass,M5); + + const int nu = 3; + std::vector twists(Nd,0); + twists[nu] = 1; + GparityDomainWallFermionR::ImplParams params; + params.twists = twists; + + GparityDomainWallFermionR Dw(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,params); + + 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(UGrid); + LatticeGaugeField tmp(UGrid); + + Dw.MDeriv(tmp , Mphi, phi,DaggerNo ); UdSdU=tmp; + Dw.MDeriv(tmp , phi, Mphi,DaggerYes ); UdSdU=(UdSdU+tmp); + + FermionField Ftmp (FGrid); + + //////////////////////////////////// + // Modify the gauge field a little + //////////////////////////////////// + RealD dt = 0.0001; + RealD Hmom = 0.0; + RealD Hmomprime = 0.0; + RealD Hmompp = 0.0; + LatticeColourMatrix mommu(UGrid); + LatticeColourMatrix forcemu(UGrid); + LatticeGaugeField mom(UGrid); + LatticeGaugeField Uprime(UGrid); + + for(int mu=0;mu(mom,mommu,mu); + + // fourth order exponential approx + parallel_for(auto i=mom.begin();i(mom,mu); + std::cout << GridLogMessage<< " Mommu " << norm2(mommu)<(UdSdU,mu); + std::cout << GridLogMessage<< " dsdumu " << norm2(mommu)<(UdSdU,mu); + mommu=Ta(mommu)*2.0; + 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); + mommu = PeekIndex(mom,mu); + + // Update PF action density + dS = dS+trace(mommu*forcemu)*dt; + + dSmom = dSmom - trace(mommu*forcemu) * dt; + dSmom2 = dSmom2 - trace(forcemu*forcemu) *(0.25* dt*dt); + + // Update mom action density + mommu = mommu + forcemu*(dt*0.5); + + Hmomprime -= real(sum(trace(mommu*mommu))); + + } + + Complex dSpred = sum(dS); + Complex dSm = sum(dSmom); + Complex dSm2 = sum(dSmom2); + + + std::cout << GridLogMessage <<"Initial mom hamiltonian is "<< Hmom < twists(Nd,0); - twists[nu] = 1; const int Ls=4; const int L =4; @@ -99,6 +97,9 @@ int main (int argc, char ** argv) ConjugateGradient CG(1.0e-8,10000); CG(HermOpEO,src_o_1f,result_o_1f); + // const int nu = 3; + std::vector twists(Nd,0); + twists[nu] = 1; GparityDomainWallFermionR::ImplParams params; params.twists = twists; GparityDomainWallFermionR GPDdwf(Umu_2f,*FGrid_2f,*FrbGrid_2f,*UGrid_2f,*UrbGrid_2f,mass,M5,params); diff --git a/tests/Test_gpdwf_force.cc b/tests/Test_gpdwf_force.cc new file mode 100644 index 00000000..d4a85d77 --- /dev/null +++ b/tests/Test_gpdwf_force.cc @@ -0,0 +1,179 @@ +#include + +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(); + + const int Ls=8; + + 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); + + int threads = GridThread::GetThreads(); + std::cout< twists(Nd,0); twists[nu] = 1; + GparityDomainWallFermionR::ImplParams params; params.twists = twists; + GparityDomainWallFermionR Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,params); + Ddwf.M (phi,Mphi); + + ComplexD S = innerProduct(Mphi,Mphi); // pdag MdagM p + + // get the deriv of phidag MdagM phi with respect to "U" + LatticeGaugeField UdSdU(UGrid); + LatticeGaugeField tmp(UGrid); + + Ddwf.MDeriv(tmp , Mphi, phi,DaggerNo ); UdSdU=tmp; + Ddwf.MDeriv(tmp , phi, Mphi,DaggerYes ); UdSdU=(UdSdU+tmp); + + FermionField Ftmp (FGrid); + + //////////////////////////////////// + // Modify the gauge field a little + //////////////////////////////////// + RealD dt = 0.001; + + LatticeColourMatrix mommu(UGrid); + LatticeColourMatrix forcemu(UGrid); + LatticeGaugeField mom(UGrid); + LatticeGaugeField Uprime(UGrid); + + for(int mu=0;mu(mom,mommu,mu); + + // fourth order exponential approx + parallel_for(auto i=mom.begin();i(UdSdU,mu); + mommu=Ta(mommu)*2.0; + PokeIndex(UdSdU,mommu,mu); + } + + for(int mu=0;mu(UdSdU,mu); + mommu = PeekIndex(mom,mu); + + // Update PF action density + dS = dS+trace(mommu*forcemu)*dt; + + } + + Complex dSpred = sum(dS); + + // From TwoFlavourPseudoFermion: + ////////////////////////////////////////////////////// + // dS/du = - phi^dag (Mdag M)^-1 [ Mdag dM + dMdag M ] (Mdag M)^-1 phi + // = - phi^dag M^-1 dM (MdagM)^-1 phi - phi^dag (MdagM)^-1 dMdag dM (Mdag)^-1 phi + // + // = - Ydag dM X - Xdag dMdag Y + // + ////////////////////////////////////////////////////// + // Our conventions really make this UdSdU; We do not differentiate wrt Udag here. + // So must take dSdU - adj(dSdU) and left multiply by mom to get dS/dt. + // + // When we have Gparity -- U and Uconj enter. + // + // dU/dt = p U + // dUc/dt = p* Uc // Is p real, traceless, etc.. + // + // dS/dt = dUdt dSdU = p U dSdU + // + // Gparity --- deriv is pc Uc dSdUc + p U dSdU + // + // Pmu = zero; + // for(int mu=0;mu::GaussianLieAlgebraMatrix(pRNG, Pmu); + // PokeIndex(P, Pmu, mu); + // } + // + // + // GridBase *grid = out._grid; + // LatticeReal ca (grid); + // LatticeMatrix la (grid); + // Complex ci(0.0,scale); + // Matrix ta; + // out=zero; + // for(int a=0;a twists(Nd,0); twists[1] = 1; + params.twists = twists; + GparityWilsonFermionR Dw(Umu,Grid,RBGrid,mass,params); FermionField src_e (&RBGrid); FermionField src_o (&RBGrid); diff --git a/tests/Test_hmc_EODWFRatio_Gparity.cc b/tests/Test_hmc_EODWFRatio_Gparity.cc new file mode 100644 index 00000000..76343317 --- /dev/null +++ b/tests/Test_hmc_EODWFRatio_Gparity.cc @@ -0,0 +1,70 @@ +#include "Grid.h" + + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + const int Ls = 8; + + 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); + + GridSerialRNG sRNG; + GridParallelRNG pRNG(UGrid); + sRNG.SeedRandomDevice(); + pRNG.SeedRandomDevice(); + + LatticeLorentzColourMatrix U(UGrid); + + SU3::HotConfiguration(pRNG, U); + + // simplify template declaration? Strip the lorentz from the second template + WilsonGaugeActionR Waction(5.6); + + Real mass=0.04; + Real pv =1.0; + RealD M5=1.5; + + typedef typename GparityDomainWallFermionR::FermionField FermionField; + + const int nu = 3; + std::vector twists(Nd,0); + twists[nu] = 1; + GparityDomainWallFermionR::ImplParams params; + params.twists = twists; + + GparityDomainWallFermionR DenOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,params); + GparityDomainWallFermionR NumOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,pv,M5,params); + + ConjugateGradient CG(1.0e-8,10000); + + TwoFlavourEvenOddRatioPseudoFermionAction Nf2(NumOp, DenOp,CG,CG); + + //Collect actions + ActionLevel Level1; + Level1.push_back(&Nf2); + Level1.push_back(&Waction); + ActionSet FullSet; + FullSet.push_back(Level1); + + // Create integrator + typedef MinimumNorm2 IntegratorType;// change here to change the algorithm + IntegratorParameters MDpar(20); + IntegratorType MDynamics(UGrid,MDpar, FullSet); + + // Create HMC + HMCparameters HMCpar; + HybridMonteCarlo HMC(HMCpar, MDynamics,sRNG,pRNG); + + // Run it + HMC.evolve(U); + +}