1
0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-04 19:25:56 +01:00

Typo fixes and check-in of G-parity action test for DWF

This commit is contained in:
David Murphy 2017-08-15 22:47:00 -04:00
parent b7f93aeb4d
commit 6d0786ff9d
3 changed files with 336 additions and 109 deletions

View File

@ -35,66 +35,66 @@ See the full license in the file "LICENSE" in the top level distribution directo
namespace Grid {
namespace QCD {
// DJM: Abstract base class for EOFA fermion types.
// Defines layout of additional EOFA-specific parameters and operators.
// Use to construct EOFA pseudofermion actions that are agnostic to Shamir / Mobius / etc.,
// and ensure that no one can construct EOFA pseudofermion action with non-EOFA fermion type.
template<class Impl>
class AbstractEOFAFermion : public CayleyFermion5D<Impl> {
public:
INHERIT_IMPL_TYPES(Impl);
// DJM: Abstract base class for EOFA fermion types.
// Defines layout of additional EOFA-specific parameters and operators.
// Use to construct EOFA pseudofermion actions that are agnostic to
// Shamir / Mobius / etc., and ensure that no one can construct EOFA
// pseudofermion action with non-EOFA fermion type.
template<class Impl>
class AbstractEOFAFermion : public CayleyFermion5D<Impl> {
public:
INHERIT_IMPL_TYPES(Impl);
public:
// Fermion operator: D(mq1) + shift*\gamma_{5}*R_{5}*\Delta_{\pm}(mq2,mq3)*P_{\pm}
RealD mq1;
RealD mq2;
RealD mq3;
RealD shift;
int pm;
public:
// Fermion operator: D(mq1) + shift*\gamma_{5}*R_{5}*\Delta_{\pm}(mq2,mq3)*P_{\pm}
RealD mq1;
RealD mq2;
RealD mq3;
RealD shift;
int pm;
RealD alpha; // Mobius scale
RealD k; // EOFA normalization constant
RealD alpha; // Mobius scale
RealD k; // EOFA normalization constant
virtual void Instantiatable(void) = 0;
virtual void Instantiatable(void) = 0;
// EOFA-specific operations
// Force user to implement in derived classes
virtual void Omega (const FermionField& in, FermionField& out, int sign, int dag) = 0;
virtual void Dtilde (const FermionField& in, FermionField& out) = 0;
virtual void DtildeInv(const FermionField& in, FermionField& out) = 0;
// EOFA-specific operations
// Force user to implement in derived classes
virtual void Omega (const FermionField& in, FermionField& out, int sign, int dag) = 0;
virtual void Dtilde (const FermionField& in, FermionField& out) = 0;
virtual void DtildeInv(const FermionField& in, FermionField& out) = 0;
// Implement derivatives in base clcass: for EOFA both DWF and Mobius just need d(Dw)/dU
virtual void MDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag){
this->DhopDeriv(mat, U, V, dag);
};
virtual void MoeDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag){
this->DhopDerivOE(mat, U, V, dag);
};
virtual void MeoDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag){
this->DhopDerivEO(mat, U, V, dag);
};
// Implement derivatives in base class:
// for EOFA both DWF and Mobius just need d(Dw)/dU
virtual void MDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag){
this->DhopDeriv(mat, U, V, dag);
};
virtual void MoeDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag){
this->DhopDerivOE(mat, U, V, dag);
};
virtual void MeoDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag){
this->DhopDerivEO(mat, U, V, dag);
};
// Recompute 5D coefficients for different value of shift constant (needed for heatbath loop over poles)
virtual void RefreshShiftCoefficients(RealD new_shift) = 0;
// Recompute 5D coefficients for different value of shift constant
// (needed for heatbath loop over poles)
virtual void RefreshShiftCoefficients(RealD new_shift) = 0;
// Constructors
AbstractEOFAFermion(GaugeField& _Umu,
GridCartesian& FiveDimGrid,
GridRedBlackCartesian& FiveDimRedBlackGrid,
GridCartesian& FourDimGrid,
GridRedBlackCartesian& FourDimRedBlackGrid,
RealD _mq1, RealD _mq2, RealD _mq3, RealD _shift, int _pm,
RealD _M5, RealD _b, RealD _c, const ImplParams &p= ImplParams()) :
CayleyFermion5D<Impl>(_Umu, FiveDimGrid, FiveDimRedBlackGrid, FourDimGrid, FourDimRedBlackGrid, _mq1, _M5, p),
mq1(_mq1), mq2(_mq2), mq3(_mq3), shift(_shift), pm(_pm)
{
int Ls = this->Ls;
this->alpha = _b + _c;
this->k = this->alpha * (_mq3 - _mq2) * std::pow(this->alpha+1.0,2*Ls) /
( std::pow(alpha+1.0,Ls) + _mq2*std::pow(alpha-1.0,Ls) ) /
( std::pow(alpha+1.0,Ls) + _mq3*std::pow(alpha-1.0,Ls) );
}
};
// Constructors
AbstractEOFAFermion(GaugeField& _Umu, GridCartesian& FiveDimGrid, GridRedBlackCartesian& FiveDimRedBlackGrid,
GridCartesian& FourDimGrid, GridRedBlackCartesian& FourDimRedBlackGrid,
RealD _mq1, RealD _mq2, RealD _mq3, RealD _shift, int _pm,
RealD _M5, RealD _b, RealD _c, const ImplParams& p=ImplParams())
: CayleyFermion5D<Impl>(_Umu, FiveDimGrid, FiveDimRedBlackGrid, FourDimGrid, FourDimRedBlackGrid,
_mq1, _M5, p), mq1(_mq1), mq2(_mq2), mq3(_mq3), shift(_shift), pm(_pm)
{
int Ls = this->Ls;
this->alpha = _b + _c;
this->k = this->alpha * (_mq3-_mq2) * std::pow(this->alpha+1.0,2*Ls) /
( std::pow(this->alpha+1.0,Ls) + _mq2*std::pow(this->alpha-1.0,Ls) ) /
( std::pow(this->alpha+1.0,Ls) + _mq3*std::pow(this->alpha-1.0,Ls) );
};
};
}}
#endif

View File

@ -35,75 +35,93 @@ See the full license in the file "LICENSE" in the top level distribution directo
namespace Grid {
namespace QCD {
// DJM: EOFA with (Shamir) domain wall fermions.
// We overload and re-implement only the routines which have a different operator
// structure than the CayleyFermion5D base class.
template<class Impl>
class DomainWallEOFAFermion : public AbstractEOFAFermion<Impl>
{
public:
INHERIT_IMPL_TYPES(Impl);
/*template<typename T> struct switcheroo {
static inline int iscomplex() { return 0; }
template<class vec>
static inline vec mult(vec a, vec b){ return real_mult(a,b); }
};
public:
// Modified (0,Ls-1) and (Ls-1,0) elements of Mooee for red-black preconditioned Shamir EOFA
Coeff_t dm;
Coeff_t dp;
template<> struct switcheroo<ComplexD> {
static inline int iscomplex() { return 1; }
template<class vec>
static inline vec mult(vec a, vec b){ return a*b; }
};
virtual void Instantiatable(void) {};
template<> struct switcheroo<ComplexF> {
static inline int iscomplex() { return 1; }
template<class vec>
static inline vec mult(vec a, vec b) { return a*b; }
};*/
// EOFA specific operators
virtual void Omega (const FermionField& in, FermionField &out, int sign, int dag);
virtual void Dtilde (const FermionField& in, FermionField &out);
virtual void DtildeInv (const FermionField& in, FermionField &out);
template<class Impl>
class DomainWallEOFAFermion : public AbstractEOFAFermion<Impl>
{
public:
INHERIT_IMPL_TYPES(Impl);
// override multiply
virtual RealD M (const FermionField& in, FermionField& out);
virtual RealD Mdag (const FermionField& in, FermionField& out);
public:
// Modified (0,Ls-1) and (Ls-1,0) elements of Mooee
// for red-black preconditioned Shamir EOFA
Coeff_t dm;
Coeff_t dp;
// half checkerboard operations
virtual void Mooee (const FermionField &in, FermionField &out);
virtual void MooeeDag (const FermionField &in, FermionField &out);
virtual void MooeeInv (const FermionField &in, FermionField &out);
virtual void MooeeInvDag(const FermionField &in, FermionField &out);
virtual void Instantiatable(void) {};
virtual void M5D (const FermionField &psi, FermionField &chi);
virtual void M5Ddag (const FermionField &psi, FermionField &chi);
// EOFA-specific operations
virtual void Omega (const FermionField& in, FermionField& out, int sign, int dag);
virtual void Dtilde (const FermionField& in, FermionField& out);
virtual void DtildeInv (const FermionField& in, FermionField& out);
/////////////////////////////////////////////////////
// Instantiate different versions depending on Impl
/////////////////////////////////////////////////////
void M5D (const FermionField &psi, const FermionField &phi, FermionField &chi,
std::vector<Coeff_t> &lower, std::vector<Coeff_t> &diag, std::vector<Coeff_t> &upper);
void M5Ddag (const FermionField &psi, const FermionField &phi, FermionField &chi,
std::vector<Coeff_t> &lower, std::vector<Coeff_t> &diag, std::vector<Coeff_t> &upper);
void MooeeInternal (const FermionField &in, FermionField &out, int dag, int inv);
void MooeeInternalCompute(int dag, int inv, Vector<iSinglet<Simd> >& Matp, Vector<iSinglet<Simd> >& Matm);
void MooeeInternalAsm (const FermionField &in, FermionField &out, int LLs, int site,
Vector<iSinglet<Simd> >& Matp, Vector<iSinglet<Simd> >& Matm);
void MooeeInternalZAsm (const FermionField &in, FermionField &out, int LLs, int site,
Vector<iSinglet<Simd> >& Matp, Vector<iSinglet<Simd> >& Matm);
// override multiply
virtual RealD M (const FermionField& in, FermionField& out);
virtual RealD Mdag (const FermionField& in, FermionField& out);
virtual void RefreshShiftCoefficients(RealD new_shift);
// half checkerboard operations
virtual void Mooee (const FermionField& in, FermionField& out);
virtual void MooeeDag (const FermionField& in, FermionField& out);
virtual void MooeeInv (const FermionField& in, FermionField& out);
virtual void MooeeInvDag(const FermionField& in, FermionField& out);
// Constructors
DomainWallEOFAFermion(GaugeField& _Umu,
GridCartesian& FiveDimGrid,
GridRedBlackCartesian& FiveDimRedBlackGrid,
GridCartesian& FourDimGrid,
GridRedBlackCartesian& FourDimRedBlackGrid,
RealD _mq1, RealD _mq2, RealD _mq3, RealD _shift,
int pm, RealD _M5, const ImplParams& p=ImplParams());
virtual void M5D (const FermionField& psi, FermionField& chi);
virtual void M5Ddag (const FermionField& psi, FermionField& chi);
protected:
virtual void SetCoefficientsInternal(RealD zolo_hi, std::vector<Coeff_t> &gamma, RealD b, RealD c);
};
/////////////////////////////////////////////////////
// Instantiate different versions depending on Impl
/////////////////////////////////////////////////////
void M5D(const FermionField& psi, const FermionField& phi, FermionField& chi,
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper);
void M5Ddag(const FermionField& psi, const FermionField& phi, FermionField& chi,
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper);
void MooeeInternal(const FermionField& in, FermionField& out, int dag, int inv);
void MooeeInternalCompute(int dag, int inv, Vector<iSinglet<Simd>>& Matp, Vector<iSinglet<Simd>>& Matm);
void MooeeInternalAsm(const FermionField& in, FermionField& out, int LLs, int site,
Vector<iSinglet<Simd>>& Matp, Vector<iSinglet<Simd>>& Matm);
void MooeeInternalZAsm(const FermionField& in, FermionField& out, int LLs, int site,
Vector<iSinglet<Simd>>& Matp, Vector<iSinglet<Simd>>& Matm);
virtual void RefreshShiftCoefficients(RealD new_shift);
// Constructors
DomainWallEOFAFermion(GaugeField& _Umu, GridCartesian& FiveDimGrid, GridRedBlackCartesian& FiveDimRedBlackGrid,
GridCartesian& FourDimGrid, GridRedBlackCartesian& FourDimRedBlackGrid,
RealD _mq1, RealD _mq2, RealD _mq3, RealD _shift, int pm,
RealD _M5, const ImplParams& p=ImplParams());
protected:
void SetCoefficientsInternal(RealD zolo_hi, std::vector<Coeff_t>& gamma, RealD b, RealD c);
};
}}
#define INSTANTIATE_DPERP_DWF_EOFA(A)\
template void DomainWallEOFAFermion<A>::M5D(const FermionField& psi, const FermionField& phi, FermionField& chi,\
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper); \
template void DomainWallEOFAFermion<A>::M5Ddag(const FermionField& psi, const FermionField& phi, FermionField& chi,\
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper); \
template void DomainWallEOFAFermion<A>::M5D(const FermionField& psi, const FermionField& phi, FermionField& chi, \
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper); \
template void DomainWallEOFAFermion<A>::M5Ddag(const FermionField& psi, const FermionField& phi, FermionField& chi, \
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper); \
template void DomainWallEOFAFermion<A>::MooeeInv(const FermionField& psi, FermionField& chi); \
template void DomainWallEOFAFermion<A>::MooeeInvDag(const FermionField& psi, FermionField& chi);

View File

@ -0,0 +1,209 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/debug/Test_reweight_dwf_eofa.cc
Copyright (C) 2017
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: paboyle <paboyle@ph.ed.ac.uk>
Author: David Murphy <dmurphy@phys.columbia.edu>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
typedef typename GparityDomainWallFermionR::FermionField FermionField;
// parameters for test
const std::vector<int> grid_dim = { 8, 8, 8, 8 };
const int Ls = 8;
const int Nhits = 10;
const int max_iter = 5000;
const RealD mf = 0.1;
const RealD mb = 0.11;
const RealD M5 = 1.8;
const RealD stop_tol = 1.0e-12;
RealD mean(const std::vector<RealD>& data)
{
int N = data.size();
RealD mean(0.0);
for(int i=0; i<N; ++i){ mean += data[i]; }
return mean/RealD(N);
}
RealD jack_mean(const std::vector<RealD>& data, int sample)
{
int N = data.size();
RealD mean(0.0);
for(int i=0; i<N; ++i){ if(i != sample){ mean += data[i]; } }
return mean/RealD(N-1);
}
RealD jack_std(const std::vector<RealD>& jacks, RealD mean)
{
int N = jacks.size();
RealD std(0.0);
for(int i=0; i<N; ++i){ std += std::pow(jacks[i]-mean, 2.0); }
return std::sqrt(RealD(N-1)/RealD(N)*std);
}
std::vector<RealD> jack_stats(const std::vector<RealD>& data)
{
int N = data.size();
std::vector<RealD> jack_samples(N);
std::vector<RealD> jack_stats(2);
jack_stats[0] = mean(data);
for(int i=0; i<N; i++){ jack_samples[i] = jack_mean(data,i); }
jack_stats[1] = jack_std(jack_samples, jack_stats[0]);
return jack_stats;
}
int main(int argc, char **argv)
{
Grid_init(&argc, &argv);
// Initialize spacetime grid
std::cout << GridLogMessage << "Lattice dimensions: "
<< grid_dim << " Ls: " << Ls << std::endl;
GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid(grid_dim,
GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi());
GridRedBlackCartesian* UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian* FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid);
GridRedBlackCartesian* FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid);
// Set up RNGs
std::vector<int> seeds4({1, 2, 3, 4});
std::vector<int> seeds5({5, 6, 7, 8});
GridParallelRNG RNG5(FGrid);
RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid);
RNG4.SeedFixedIntegers(seeds4);
// Random gauge field
LatticeGaugeField Umu(UGrid);
SU3::HotConfiguration(RNG4, Umu);
// Initialize RHMC fermion operators
GparityDomainWallFermionR::ImplParams params;
GparityDomainWallFermionR Ddwf_f(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, M5, params);
GparityDomainWallFermionR Ddwf_b(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mb, M5, params);
SchurDiagMooeeOperator<GparityDomainWallFermionR, FermionField> MdagM(Ddwf_f);
SchurDiagMooeeOperator<GparityDomainWallFermionR, FermionField> VdagV(Ddwf_b);
// Degree 12 rational approximations to x^(1/4) and x^(-1/4)
double lo = 0.0001;
double hi = 95.0;
int precision = 64;
int degree = 12;
AlgRemez remez(lo, hi, precision);
std::cout << GridLogMessage << "Generating degree " << degree << " for x^(1/4)" << std::endl;
remez.generateApprox(degree, 1, 4);
MultiShiftFunction PowerQuarter(remez, stop_tol, false);
MultiShiftFunction PowerNegQuarter(remez, stop_tol, true);
// Stochastically estimate reweighting factor via RHMC
RealD scale = std::sqrt(0.5);
std::vector<RealD> rw_rhmc(Nhits);
ConjugateGradientMultiShift<FermionField> msCG_V(max_iter, PowerQuarter);
ConjugateGradientMultiShift<FermionField> msCG_M(max_iter, PowerNegQuarter);
std::cout.precision(12);
for(int hit=0; hit<Nhits; hit++){
// Gaussian source
FermionField Phi (Ddwf_f.FermionGrid());
FermionField PhiOdd (Ddwf_f.FermionRedBlackGrid());
std::vector<FermionField> tmp(2, Ddwf_f.FermionRedBlackGrid());
gaussian(RNG5, Phi);
Phi = Phi*scale;
pickCheckerboard(Odd, PhiOdd, Phi);
// evaluate -log(rw)
msCG_V(VdagV, PhiOdd, tmp[0]);
msCG_M(MdagM, tmp[0], tmp[1]);
rw_rhmc[hit] = norm2(tmp[1]) - norm2(PhiOdd);
std::cout << std::endl << "==================================================" << std::endl;
std::cout << " --- RHMC: Hit " << hit << ": rw = " << rw_rhmc[hit];
std::cout << std::endl << "==================================================" << std::endl << std::endl;
}
// Initialize EOFA fermion operators
RealD shift_L = 0.0;
RealD shift_R = -1.0;
int pm = 1;
GparityDomainWallEOFAFermionR Deofa_L(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, mf, mb, shift_L, pm, M5, params);
GparityDomainWallEOFAFermionR Deofa_R(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mb, mf, mb, shift_R, pm, M5, params);
MdagMLinearOperator<GparityDomainWallEOFAFermionR, FermionField> LdagL(Deofa_L);
MdagMLinearOperator<GparityDomainWallEOFAFermionR, FermionField> RdagR(Deofa_R);
// Stochastically estimate reweighting factor via EOFA
RealD k = Deofa_L.k;
std::vector<RealD> rw_eofa(Nhits);
ConjugateGradient<FermionField> CG(stop_tol, max_iter);
SchurRedBlackDiagMooeeSolve<FermionField> SchurSolver(CG);
for(int hit=0; hit<Nhits; hit++){
// Gaussian source
FermionField Phi (Deofa_L.FermionGrid());
FermionField spProj_Phi(Deofa_L.FermionGrid());
std::vector<FermionField> tmp(2, Deofa_L.FermionGrid());
gaussian(RNG5, Phi);
Phi = Phi*scale;
// evaluate -log(rw)
// LH term
for(int s=0; s<Ls; ++s){ axpby_ssp_pminus(spProj_Phi, 0.0, Phi, 1.0, Phi, s, s); }
Deofa_L.Omega(spProj_Phi, tmp[0], -1, 0);
G5R5(tmp[1], tmp[0]);
tmp[0] = zero;
SchurSolver(Deofa_L, tmp[1], tmp[0]);
Deofa_L.Omega(tmp[0], tmp[1], -1, 1);
rw_eofa[hit] = -k*innerProduct(spProj_Phi,tmp[1]).real();
// RH term
for(int s=0; s<Ls; ++s){ axpby_ssp_pplus(spProj_Phi, 0.0, Phi, 1.0, Phi, s, s); }
Deofa_R.Omega(spProj_Phi, tmp[0], 1, 0);
G5R5(tmp[1], tmp[0]);
tmp[0] = zero;
SchurSolver(Deofa_R, tmp[1], tmp[0]);
Deofa_R.Omega(tmp[0], tmp[1], 1, 1);
rw_eofa[hit] += k*innerProduct(spProj_Phi,tmp[1]).real();
std::cout << std::endl << "==================================================" << std::endl;
std::cout << " --- EOFA: Hit " << hit << ": rw = " << rw_eofa[hit];
std::cout << std::endl << "==================================================" << std::endl << std::endl;
}
std::vector<RealD> rhmc_result = jack_stats(rw_rhmc);
std::vector<RealD> eofa_result = jack_stats(rw_eofa);
std::cout << std::endl << "RHMC: rw = " << rhmc_result[0] << " +/- " << rhmc_result[1] << std::endl;
std::cout << std::endl << "EOFA: rw = " << eofa_result[0] << " +/- " << eofa_result[1] << std::endl;
Grid_finalize();
}