portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-09 23:45:36 +00:00
Code Releases Activity

Merge branch 'develop' of https://github.com/paboyle/Grid into develop

Browse Source
This commit is contained in:
u61464 2021-05-04 08:40:38 -07:00
commit 15ae317858
25 changed files with 358 additions and 159 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
Grid/cshift/Cshift_mpi.h
View File

@ -122,8 +122,8 @@ template<class vobj> void Cshift_comms(Lattice<vobj> &ret,const Lattice<vobj> &r
assert(shift<fd);
int buffer_size = rhs.Grid()->_slice_nblock[dimension]*rhs.Grid()->_slice_block[dimension];
cshiftVector<vobj> send_buf(buffer_size);
cshiftVector<vobj> recv_buf(buffer_size);
static cshiftVector<vobj> send_buf; send_buf.resize(buffer_size);
static cshiftVector<vobj> 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<class vobj> void Cshift_comms_simd(Lattice<vobj> &ret,const Lattice<vo
int buffer_size = grid->_slice_nblock[dimension]*grid->_slice_block[dimension];
// int words = sizeof(vobj)/sizeof(vector_type);
std::vector<cshiftVector<scalar_object> > send_buf_extract(Nsimd);
std::vector<cshiftVector<scalar_object> > recv_buf_extract(Nsimd);
static std::vector<cshiftVector<scalar_object> > send_buf_extract; send_buf_extract.resize(Nsimd);
static std::vector<cshiftVector<scalar_object> > 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<class vobj> void Cshift_comms(Lattice<vobj> &ret,const Lattice<vobj> &r
assert(shift<fd);
int buffer_size = rhs.Grid()->_slice_nblock[dimension]*rhs.Grid()->_slice_block[dimension];
cshiftVector<vobj> send_buf_v(buffer_size);
cshiftVector<vobj> recv_buf_v(buffer_size);
static cshiftVector<vobj> send_buf_v; send_buf_v.resize(buffer_size);
static cshiftVector<vobj> recv_buf_v; recv_buf_v.resize(buffer_size);
vobj *send_buf;
vobj *recv_buf;
{
@ -381,8 +381,8 @@ template<class vobj> void Cshift_comms_simd(Lattice<vobj> &ret,const Lattice<vo
int buffer_size = grid->_slice_nblock[dimension]*grid->_slice_block[dimension];
// int words = sizeof(vobj)/sizeof(vector_type);
std::vector<cshiftVector<scalar_object> > send_buf_extract(Nsimd);
std::vector<cshiftVector<scalar_object> > recv_buf_extract(Nsimd);
static std::vector<cshiftVector<scalar_object> > send_buf_extract; send_buf_extract.resize(Nsimd);
static std::vector<cshiftVector<scalar_object> > recv_buf_extract; recv_buf_extract.resize(Nsimd);
scalar_object * recv_buf_extract_mpi;
scalar_object * send_buf_extract_mpi;
{

2
Grid/parallelIO/MetaData.h
View File

@ -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;

17
Grid/parallelIO/NerscIO.h
View File

@ -205,11 +205,20 @@ public:
std::cout<<GridLogMessage <<"NERSC Configuration "<<file<< " and plaquette, link trace, and checksum agree"<<std::endl;
}
// Preferred interface
template<class GaugeStats=PeriodicGaugeStatistics>
static inline void writeConfiguration(Lattice<vLorentzColourMatrixD > &Umu,
std::string file,
std::string ens_label = std::string("DWF"))
{
writeConfiguration(Umu,file,0,1,ens_label);
}
template<class GaugeStats=PeriodicGaugeStatistics>
static inline void writeConfiguration(Lattice<vLorentzColourMatrixD > &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");

16
Grid/qcd/action/fermion/StaggeredImpl.h
View File

@ -72,19 +72,23 @@ public:
StaggeredImpl(const ImplParams &p = ImplParams()) : Params(p){};
static accelerator_inline void multLink(SiteSpinor &phi,
template<class _Spinor>
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<class _Spinor>
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 <class ref>

35
Grid/qcd/action/fermion/WilsonImpl.h
View File

@ -184,18 +184,22 @@ public:
mat = TraceIndex<SpinIndex>(P);
}
inline void extractLinkField(std::vector<GaugeLinkField> &mat, DoubledGaugeField &Uds){
inline void extractLinkField(std::vector<GaugeLinkField> &mat, DoubledGaugeField &Uds)
{
for (int mu = 0; mu < Nd; mu++)
mat[mu] = PeekIndex<LorentzIndex>(Uds, mu);
}
inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde,int mu){
inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde,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<LorentzIndex>(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<Ls;s++){
int sF = s+Ls*sU;
for(int spn=0;spn<Ns;spn++){ //sum over spin
auto bb = coalescedRead(Btilde_v[sF]()(spn) ); //color vector
auto aa = coalescedRead(Atilde_v[sF]()(spn) );
auto op = outerProduct(bb,aa);
sum = sum + op;
}
}
coalescedWrite(mat_v[sU](mu)(), sum);
});
}
#endif
}
};

14
Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h
View File

@ -880,11 +880,23 @@ void CayleyFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
}
std::vector<RealD> G_s(Ls,1.0);
Integer sign = 1; // sign flip for vector/tadpole
if ( curr_type == Current::Axial ) {
for(int s=0;s<Ls/2;s++){
G_s[s] = -1.0;
}
}
else if ( curr_type == Current::Tadpole ) {
auto b=this->_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<Ls;s++){
@ -907,7 +919,7 @@ void CayleyFermion5D<Impl>::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

76
Grid/qcd/action/fermion/implementation/StaggeredKernelsImplementation.h
View File

@ -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 = &chi; \
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 = &chi; \
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 <class Impl>
@ -84,12 +77,14 @@ void StaggeredKernels<Impl>::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<LLs;s++){
//
@ -118,7 +113,7 @@ void StaggeredKernels<Impl>::DhopSiteGeneric(StencilView &st,
if ( dag ) {
Uchi = - Uchi;
}
vstream(out[sF], Uchi);
coalescedWrite(out[sF], Uchi,lane);
}
};
@ -130,13 +125,16 @@ template <int Naik> accelerator_inline
void StaggeredKernels<Impl>::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<LLs;s++){
// int sF=LLs*sU+s;
@ -165,7 +163,7 @@ void StaggeredKernels<Impl>::DhopSiteGenericInt(StencilView &st,
if ( dag ) {
Uchi = - Uchi;
}
vstream(out[sF], Uchi);
coalescedWrite(out[sF], Uchi,lane);
}
};
@ -178,14 +176,17 @@ template <int Naik> accelerator_inline
void StaggeredKernels<Impl>::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<LLs;s++){
// int sF=LLs*sU+s;
@ -211,11 +212,12 @@ void StaggeredKernels<Impl>::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<Impl>::DhopImproved(StencilImpl &st, LebesgueOrder &lo,
GridBase *FGrid=in.Grid();
GridBase *UGrid=U.Grid();
typedef StaggeredKernels<Impl> 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<Impl>::DhopNaive(StencilImpl &st, LebesgueOrder &lo,
GridBase *FGrid=in.Grid();
GridBase *UGrid=U.Grid();
typedef StaggeredKernels<Impl> 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);

21
Grid/qcd/smearing/StoutSmearing.h
View File

@ -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";
};

8
Grid/qcd/utils/Metric.h
View File

@ -93,13 +93,13 @@ public:
GeneralisedMomenta(GridBase* grid, Metric<MomentaField>& 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

21
Grid/tensors/Tensor_outer.h
View File

@ -34,6 +34,16 @@ NAMESPACE_BEGIN(Grid);
// outerProduct Scalar x Scalar -> Scalar
// Vector x Vector -> Matrix
///////////////////////////////////////////////////////////////////////////////////////
template<class CC,IfComplex<CC> = 0>
accelerator_inline CC outerProduct(const CC &l, const CC& r)
{
return l*conj(r);
}
template<class RR,IfReal<RR> = 0>
accelerator_inline RR outerProduct(const RR &l, const RR& r)
{
return l*r;
}
template<class l,class r,int N> accelerator_inline
auto outerProduct (const iVector<l,N>& lhs,const iVector<r,N>& rhs) -> iMatrix<decltype(outerProduct(lhs._internal[0],rhs._internal[0])),N>
@ -57,17 +67,6 @@ auto outerProduct (const iScalar<l>& lhs,const iScalar<r>& rhs) -> iScalar<declt
return ret;
}
template<class CC,IfComplex<CC> = 0>
accelerator_inline CC outerProduct(const CC &l, const CC& r)
{
return l*conj(r);
}
template<class RR,IfReal<RR> = 0>
accelerator_inline RR outerProduct(const RR &l, const RR& r)
{
return l*r;
}
NAMESPACE_END(Grid);
#endif

3
Grid/util/Init.cc
View File

@ -153,6 +153,9 @@ void GridCmdOptionIntVector(const std::string &str,VectorInt & vec)
return;
}
template void GridCmdOptionIntVector(const std::string &str,std::vector<int> & vec);
template void GridCmdOptionIntVector(const std::string &str,Coordinate & vec);
void GridCmdOptionInt(std::string &str,int & val)
{
std::stringstream ss(str);

6
HMC/Mobius2p1fRHMC.cc
View File

@ -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);

6
tests/debug/Test_heatbath_dwf_eofa.cc
View File

@ -66,7 +66,9 @@ int main(int argc, char** argv)
// Set up RNGs
std::vector<int> seeds4({1, 2, 3, 4});
std::vector<int> 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<LatticeFermion> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> Meofa(Lop, Rop, CG, Params, false);
Meofa.refresh(Umu, RNG5);
Meofa.refresh(Umu,sRNG, RNG5);
printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
}
@ -94,7 +96,7 @@ int main(int argc, char** argv)
ConjugateGradient<LatticeFermion> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> Meofa(Lop, Rop, CG, Params, true);
Meofa.refresh(Umu, RNG5);
Meofa.refresh(Umu,sRNG, RNG5);
printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
}

7
tests/debug/Test_heatbath_dwf_eofa_gparity.cc
View File

@ -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<FermionField> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, Params, false);
Meofa.refresh(Umu, RNG5);
Meofa.refresh(Umu,sRNG, RNG5);
printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
}
@ -100,7 +103,7 @@ int main(int argc, char** argv)
ConjugateGradient<FermionField> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, Params, true);
Meofa.refresh(Umu, RNG5);
Meofa.refresh(Umu,sRNG, RNG5);
printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
}

6
tests/debug/Test_heatbath_mobius_eofa.cc
View File

@ -68,8 +68,10 @@ int main(int argc, char** argv)
// Set up RNGs
std::vector<int> seeds4({1, 2, 3, 4});
std::vector<int> 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<LatticeFermion> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> Meofa(Lop, Rop, CG, Params, false);
Meofa.refresh(Umu, RNG5);
Meofa.refresh(Umu, sRNG,RNG5);
printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
}
@ -96,7 +98,7 @@ int main(int argc, char** argv)
ConjugateGradient<LatticeFermion> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> Meofa(Lop, Rop, CG, Params, true);
Meofa.refresh(Umu, RNG5);
Meofa.refresh(Umu, sRNG,RNG5);
printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
}

6
tests/debug/Test_heatbath_mobius_eofa_gparity.cc
View File

@ -73,7 +73,9 @@ int main(int argc, char** argv)
std::vector<int> seeds4({1, 2, 3, 4});
std::vector<int> 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<FermionField> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, Params, false);
Meofa.refresh(Umu, RNG5);
Meofa.refresh(Umu, sRNG, RNG5);
printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
}
@ -101,7 +103,7 @@ int main(int argc, char** argv)
ConjugateGradient<FermionField> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, Params, true);
Meofa.refresh(Umu, RNG5);
Meofa.refresh(Umu, sRNG, RNG5);
printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
}

4
tests/forces/Test_dwf_force_eofa.cc
View File

@ -86,7 +86,9 @@ int main (int argc, char** argv)
ConjugateGradient<LatticeFermion> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> 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"

123
tests/forces/Test_dwf_gpforce.cc
View File

@ -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<Nc>::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<LorentzIndex>(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 <<std::endl;
Dw.ImportGauge(Uprime);
Dw.M (phi,MphiPrime);
@ -145,53 +172,60 @@ int main (int argc, char ** argv)
// Use derivative to estimate dS
//////////////////////////////////////////////
for(int mu=0;mu<Nd;mu++){
std::cout << "" <<std::endl;
mommu = PeekIndex<LorentzIndex>(mom,mu);
std::cout << GridLogMessage<< " Mommu " << norm2(mommu)<<std::endl;
mommu = mommu+adj(mommu);
std::cout << GridLogMessage<< " Mommu + Mommudag " << norm2(mommu)<<std::endl;
mommu = PeekIndex<LorentzIndex>(UdSdU,mu);
std::cout << GridLogMessage<< " dsdumu " << norm2(mommu)<<std::endl;
mommu = mommu+adj(mommu);
std::cout << GridLogMessage<< " dsdumu + dag " << norm2(mommu)<<std::endl;
}
//
// Ta has 1/2([ F - adj(F) ])_traceless and want the UdSdU _and_ UdagdSdUdag terms so 2x.
//
LatticeComplex dS(UGrid); dS = Zero();
LatticeComplex dSmom(UGrid); dSmom = Zero();
LatticeComplex dSmom2(UGrid); dSmom2 = Zero();
for(int mu=0;mu<Nd;mu++){
mommu = PeekIndex<LorentzIndex>(UdSdU,mu);
mommu=Ta(mommu)*2.0;
mommu=Ta(mommu); // projectForce , GaugeImplTypes.h
PokeIndex<LorentzIndex>(UdSdU,mommu,mu);
}
for(int mu=0;mu<Nd;mu++){
mommu = PeekIndex<LorentzIndex>(mom,mu);
std::cout << GridLogMessage<< " Mommu " << norm2(mommu)<<std::endl;
mommu = mommu+adj(mommu);
std::cout << GridLogMessage<< " Mommu + Mommudag " << norm2(mommu)<<std::endl;
mommu = PeekIndex<LorentzIndex>(UdSdU,mu);
std::cout << GridLogMessage<< " dsdumu " << norm2(mommu)<<std::endl;
mommu = mommu+adj(mommu);
std::cout << GridLogMessage<< " dsdumu + dag " << norm2(mommu)<<std::endl;
}
for(int mu=0;mu<Nd;mu++){
forcemu = PeekIndex<LorentzIndex>(UdSdU,mu);
mUdSdUmu= PeekIndex<LorentzIndex>(UdSdU,mu);
mommu = PeekIndex<LorentzIndex>(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<<std::endl;
std::cout << GridLogMessage <<"dSm2 "<< dSm2<<std::endl;
std::cout << GridLogMessage <<"Initial mom hamiltonian is "<< Hmom <<std::endl;
std::cout << GridLogMessage <<"Final mom hamiltonian is "<< Hmomprime <<std::endl;
std::cout << GridLogMessage <<"Delta mom hamiltonian is "<< Hmomprime-Hmom <<std::endl;
std::cout << GridLogMessage << " S "<<S<<std::endl;
std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;
std::cout << GridLogMessage << "dS "<<Sprime-S<<std::endl;
std::cout << GridLogMessage << "predict dS "<< dSpred <<std::endl;
std::cout << GridLogMessage <<"dSm "<< dSm<<std::endl;
std::cout << GridLogMessage <<"dSm2"<< dSm2<<std::endl;
std::cout << GridLogMessage <<"Delta mom hamiltonian is "<< Hmomprime-Hmom <<std::endl;
std::cout << GridLogMessage <<"predict Delta mom hamiltonian is "<< dSm+dSm2 <<std::endl;
std::cout << GridLogMessage << "Initial S "<<S<<std::endl;
std::cout << GridLogMessage << "Final S "<<Sprime<<std::endl;
std::cout << GridLogMessage << "Delta S "<<Sprime-S<<std::endl;
std::cout << GridLogMessage << "predict delta S"<< dSpred <<std::endl;
std::cout << GridLogMessage << "defect "<< Sprime-S-dSpred <<std::endl;
std::cout << GridLogMessage << "Total dS "<< Hmomprime - Hmom + Sprime - S <<std::endl;
std::cout << GridLogMessage << "dS - dt^2 term "<< Hmomprime - Hmom + Sprime - S - dSm2 <<std::endl;
assert( fabs(real(Sprime-S-dSpred)) < 5.0 ) ;
std::cout<< GridLogMessage << "Done" <<std::endl;

3
tests/forces/Test_dwf_gpforce_eofa.cc
View File

@ -90,7 +90,8 @@ int main (int argc, char** argv)
ConjugateGradient<FermionField> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> 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"

5
tests/forces/Test_gp_rect_force.cc
View File

@ -29,7 +29,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
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<int> twists(Nd,0);
twists[nu] = 1;
ConjugateGimplD::setDirections(twists);
ConjugatePlaqPlusRectangleActionR Action(beta,c1);
//ConjugateWilsonGaugeActionR Action(beta);
//WilsonGaugeActionR Action(beta);

3
tests/forces/Test_laplacian_force.cc
View File

@ -46,6 +46,7 @@ int main (int argc, char ** argv)
std::vector<int> seeds({1,2,3,4});
GridSerialRNG sRNG; sRNG.SeedFixedIntegers({4,5,6,7});
GridParallelRNG pRNG(&Grid);
pRNG.SeedFixedIntegers(std::vector<int>({15,91,21,3}));
@ -67,7 +68,7 @@ int main (int argc, char ** argv)
LaplacianAdjointField<PeriodicGimplR> Laplacian(&Grid, CG, LapPar, Kappa);
GeneralisedMomenta<PeriodicGimplR> 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);

109
tests/forces/Test_mobius_force.cc
View File

@ -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<Nd;mu++){
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(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 <<std::endl;
Ddwf.ImportGauge(Uprime);
Ddwf.M (phi,MphiPrime);
@ -125,32 +153,87 @@ int main (int argc, char ** argv)
// Use derivative to estimate dS
//////////////////////////////////////////////
LatticeComplex dS(UGrid); dS = Zero();
LatticeComplex dSmom(UGrid); dSmom = Zero();
LatticeComplex dSmom2(UGrid); dSmom2 = Zero();
for(int mu=0;mu<Nd;mu++){
mommu = PeekIndex<LorentzIndex>(UdSdU,mu);
mommu=Ta(mommu)*2.0;
mommu=Ta(mommu);
PokeIndex<LorentzIndex>(UdSdU,mommu,mu);
}
for(int mu=0;mu<Nd;mu++){
forcemu = PeekIndex<LorentzIndex>(UdSdU,mu);
mUdSdUmu= PeekIndex<LorentzIndex>(UdSdU,mu);
mommu = PeekIndex<LorentzIndex>(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 "<<S<<std::endl;
std::cout << GridLogMessage << " -- Sprime "<<Sprime<<std::endl;
std::cout << GridLogMessage << "dS "<<Sprime-S<<std::endl;
std::cout << GridLogMessage << "predict dS "<< dSpred <<std::endl;
std::cout << GridLogMessage <<"dSm "<< dSm<<std::endl;
std::cout << GridLogMessage <<"dSm2 "<< dSm2<<std::endl;
std::cout << GridLogMessage <<"Initial mom hamiltonian is "<< Hmom <<std::endl;
std::cout << GridLogMessage <<"Final mom hamiltonian is "<< Hmomprime <<std::endl;
std::cout << GridLogMessage <<"Delta mom hamiltonian is "<< Hmomprime-Hmom <<std::endl;
std::cout << GridLogMessage <<"predict Delta mom hamiltonian is "<< dSm+dSm2 <<std::endl;
std::cout << GridLogMessage << "Initial S "<<S<<std::endl;
std::cout << GridLogMessage << "Final S "<<Sprime<<std::endl;
std::cout << GridLogMessage << "Delta S "<<Sprime-S<<std::endl;
std::cout << GridLogMessage << "predict delta S"<< dSpred <<std::endl;
std::cout << GridLogMessage << "defect "<< Sprime-S-dSpred <<std::endl;
std::cout << GridLogMessage << "Total dS "<< Hmomprime - Hmom + Sprime - S <<std::endl;
std::cout << GridLogMessage << "dS - dt^2 term "<< Hmomprime - Hmom + Sprime - S - dSm2 <<std::endl;
assert( fabs(real(Sprime-S-dSpred)) < 1.0 ) ;
std::cout<< GridLogMessage << "Done" <<std::endl;
Grid_finalize();
}

3
tests/forces/Test_mobius_force_eofa.cc
View File

@ -88,7 +88,8 @@ int main (int argc, char** argv)
ConjugateGradient<LatticeFermion> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> 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"

3
tests/forces/Test_mobius_gpforce_eofa.cc
View File

@ -93,7 +93,8 @@ int main (int argc, char** argv)
ConjugateGradient<FermionField> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> 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"

4
tests/forces/Test_momentum_filter.cc
View File

@ -61,7 +61,9 @@ int main (int argc, char ** argv)
std::vector<int> seeds({1,2,3,4});
GridParallelRNG pRNG(&Grid);
GridSerialRNG sRNG;
pRNG.SeedFixedIntegers(seeds);
sRNG.SeedFixedIntegers(seeds);
typedef PeriodicGimplR Gimpl;
typedef WilsonGaugeAction<Gimpl> 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();