From 5c364f80823fe58d2d9edb78c2cb7a9a0776adc6 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Tue, 18 Aug 2015 14:40:08 +0100 Subject: [PATCH] One flavour rational unprec added; untested but does compile. Moving param structs into a single header for later connection to file I/O using macromagic.h --- TODO | 5 +- lib/algorithms/approx/MultiShiftFunction.h | 19 +- lib/qcd/action/ActionParams.h | 26 +++ lib/qcd/action/Actions.h | 9 +- lib/qcd/action/fermion/FermionOperator.h | 4 +- lib/qcd/action/fermion/FermionOperatorImpl.h | 7 +- .../action/pseudofermion/OneFlavourRational.h | 170 ++++++++++++++++++ lib/qcd/action/pseudofermion/TwoFlavour.h | 79 -------- .../action/pseudofermion/TwoFlavourEvenOdd.h | 7 +- .../pseudofermion/TwoFlavourEvenOddRatio.h | 11 +- 10 files changed, 238 insertions(+), 99 deletions(-) create mode 100644 lib/qcd/action/ActionParams.h create mode 100644 lib/qcd/action/pseudofermion/OneFlavourRational.h diff --git a/TODO b/TODO index 432b9028..dc7065d1 100644 --- a/TODO +++ b/TODO @@ -21,8 +21,11 @@ fill in: - Force Gradient - Multi-timescale looks broken and operating on single timescale for now. Fix/debug/rewrite this + - Sign of force term. + - Prefer "RefreshInternal" or such like to "init" in naming + - Rename "Ta" as too unclear -- MacroMagic -> readers +- MacroMagic -> virtual reader class. - Link smearing/boundary conds; Policy class based implementation diff --git a/lib/algorithms/approx/MultiShiftFunction.h b/lib/algorithms/approx/MultiShiftFunction.h index fceed2eb..c9ed39e3 100644 --- a/lib/algorithms/approx/MultiShiftFunction.h +++ b/lib/algorithms/approx/MultiShiftFunction.h @@ -11,20 +11,27 @@ public: std::vector tolerances; RealD norm; RealD lo,hi; + MultiShiftFunction(int n,RealD _lo,RealD _hi): poles(n), residues(n), lo(_lo), hi(_hi) {;}; RealD approx(RealD x); void csv(std::ostream &out); void gnuplot(std::ostream &out); - MultiShiftFunction(AlgRemez & remez,double tol,bool inverse) : - order(remez.getDegree()), - tolerances(remez.getDegree(),tol), - poles(remez.getDegree()), - residues(remez.getDegree()) + + void Init(AlgRemez & remez,double tol,bool inverse) { + order=remez.getDegree(); + tolerances.resize(remez.getDegree(),tol); + poles.resize(remez.getDegree()); + residues.resize(remez.getDegree()); remez.getBounds(lo,hi); if ( inverse ) remez.getIPFE (&residues[0],&poles[0],&norm); - else remez.getPFE (&residues[0],&poles[0],&norm); + else remez.getPFE (&residues[0],&poles[0],&norm); } + MultiShiftFunction(AlgRemez & remez,double tol,bool inverse) + { + Init(remez,tol,inverse); + } + }; } #endif diff --git a/lib/qcd/action/ActionParams.h b/lib/qcd/action/ActionParams.h new file mode 100644 index 00000000..d866339c --- /dev/null +++ b/lib/qcd/action/ActionParams.h @@ -0,0 +1,26 @@ +#ifndef GRID_QCD_ACTION_PARAMS_H +#define GRID_QCD_ACTION_PARAMS_H + +namespace Grid { +namespace QCD { + + // These can move into a params header and be given MacroMagic serialisation + struct GparityWilsonImplParams { + std::vector twists; + }; + + struct WilsonImplParams { }; + + struct OneFlavourRationalParams { + RealD lo; + RealD hi; + int precision=64; + int degree=10; + RealD tolerance; // Vector? + RealD MaxIter; // Vector? + OneFlavourRationalParams (RealD lo,RealD hi,int precision=64,int degree = 10); + }; + +}} + +#endif diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h index bc9526ca..3e62a509 100644 --- a/lib/qcd/action/Actions.h +++ b/lib/qcd/action/Actions.h @@ -14,6 +14,7 @@ // Abstract base interface //////////////////////////////////////////// #include +#include //////////////////////////////////////////// // Gauge Actions @@ -157,9 +158,9 @@ typedef DomainWallFermion GparityDomainWallFermionD; #include //Todo: RHMC -//#include -//#include -//#include -//#include +#include +//#include +//#include +//#include #endif diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h index e728909e..7cdfc0d1 100644 --- a/lib/qcd/action/fermion/FermionOperator.h +++ b/lib/qcd/action/fermion/FermionOperator.h @@ -32,6 +32,8 @@ namespace Grid { virtual RealD Mdag (const FermionField &in, FermionField &out)=0; // half checkerboard operaions + virtual int ConstEE(void) { return 1; }; // clover returns zero as EE depends on gauge field + virtual void Meooe (const FermionField &in, FermionField &out)=0; virtual void MeooeDag (const FermionField &in, FermionField &out)=0; virtual void Mooee (const FermionField &in, FermionField &out)=0; @@ -49,7 +51,7 @@ namespace Grid { virtual void MDeriv (GaugeField &mat,const FermionField &U,const FermionField &V,int dag){DhopDeriv(mat,U,V,dag);}; virtual void MoeDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag){DhopDerivOE(mat,U,V,dag);}; virtual void MeoDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag){DhopDerivEO(mat,U,V,dag);}; - virtual void MooDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag){mat=zero;}; + virtual void MooDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag){mat=zero;}; // Clover can override these virtual void MeeDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag){mat=zero;}; virtual void DhopDeriv (GaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0; diff --git a/lib/qcd/action/fermion/FermionOperatorImpl.h b/lib/qcd/action/fermion/FermionOperatorImpl.h index 24c64149..1d90cace 100644 --- a/lib/qcd/action/fermion/FermionOperatorImpl.h +++ b/lib/qcd/action/fermion/FermionOperatorImpl.h @@ -5,6 +5,7 @@ namespace Grid { namespace QCD { + ////////////////////////////////////////////// // Template parameter class constructs to package // externally control Fermion implementations @@ -126,8 +127,7 @@ namespace Grid { typedef Lattice DoubledGaugeField; typedef WilsonCompressor Compressor; - - typedef struct WilsonImplParams { } ImplParams; + typedef WilsonImplParams ImplParams; ImplParams Params; WilsonImpl(const ImplParams &p= ImplParams()) : Params(p) {}; @@ -177,6 +177,7 @@ PARALLEL_FOR_LOOP //////////////////////////////////////////////////////////////////////////////////////// // Flavour doubled spinors; is Gparity the only? what about C*? //////////////////////////////////////////////////////////////////////////////////////// + template class GparityWilsonImpl : public ImplGauge { public: @@ -198,7 +199,7 @@ PARALLEL_FOR_LOOP typedef WilsonCompressor Compressor; - typedef struct GparityWilsonImplParams {std::vector twists; } ImplParams; + typedef GparityWilsonImplParams ImplParams; ImplParams Params; GparityWilsonImpl(const ImplParams &p= ImplParams()) : Params(p) {}; diff --git a/lib/qcd/action/pseudofermion/OneFlavourRational.h b/lib/qcd/action/pseudofermion/OneFlavourRational.h new file mode 100644 index 00000000..e7596065 --- /dev/null +++ b/lib/qcd/action/pseudofermion/OneFlavourRational.h @@ -0,0 +1,170 @@ +#ifndef QCD_PSEUDOFERMION_ONE_FLAVOUR_RATIONAL_H +#define QCD_PSEUDOFERMION_ONE_FLAVOUR_RATIONAL_H + +namespace Grid{ + namespace QCD{ + + /////////////////////////////////////// + // One flavour rational + /////////////////////////////////////// + + // S_f = chi^dag * N(M^dag*M)/D(M^dag*M) * chi + // + // Here, M is some operator + // N and D makeup the rat. poly + // + + template + class OneFlavourRationalPseudoFermionAction : public Action { + public: + INHERIT_IMPL_TYPES(Impl); + + typedef OneFlavourRationalParams Params; + Params param; + + MultiShiftFunction PowerHalf ; + MultiShiftFunction PowerNegHalf; + MultiShiftFunction PowerQuarter; + MultiShiftFunction PowerNegQuarter; + + private: + + FermionOperator & FermOp;// the basic operator + + // NOT using "Nroots"; IroIro is -- perhaps later, but this wasn't good for us historically + // and hasenbusch works better + + FermionField Phi; // the pseudo fermion field for this trajectory + + public: + + OneFlavourRationalPseudoFermionAction(FermionOperator &Op, + Params & p + ) : FermOp(Op), Phi(Op.FermionGrid()), param(p) + { + AlgRemez remez(param.lo,param.hi,param.precision); + + // MdagM^(+- 1/2) + std::cout< sig^2 = 0.5. + // + // So eta should be of width sig = 1/sqrt(2). + + RealD scale = std::sqrt(0.5); + + FermionField eta(FermOp.FermionGrid()); + + gaussian(pRNG,eta); + + FermOp.ImportGauge(U); + + // mutishift CG + MdagMLinearOperator ,FermionField> MdagMOp(FermOp); + ConjugateGradientMultiShift msCG(param.MaxIter,PowerQuarter); + msCG(MdagMOp,eta,Phi); + + Phi=Phi*scale; + + }; + + ////////////////////////////////////////////////////// + // S = phi^dag (Mdag M)^-1/2 phi + ////////////////////////////////////////////////////// + virtual RealD S(const GaugeField &U) { + + FermOp.ImportGauge(U); + + FermionField Y(FermOp.FermionGrid()); + + MdagMLinearOperator ,FermionField> MdagMOp(FermOp); + + ConjugateGradientMultiShift msCG(param.MaxIter,PowerNegQuarter); + + msCG(MdagMOp,Phi,Y); + + RealD action = norm2(Y); + std::cout << GridLogMessage << "Pseudofermion action FIXME -- is -1/4 solve or -1/2 solve faster??? "< MPhi_k (Npole,FermOp.FermionGrid()); + + FermionField X(FermOp.FermionGrid()); + FermionField Y(FermOp.FermionGrid()); + + GaugeField tmp(FermOp.GaugeGrid()); + + FermOp.ImportGauge(U); + + MdagMLinearOperator ,FermionField> MdagMOp(FermOp); + + ConjugateGradientMultiShift msCG(param.MaxIter,PowerNegHalf); + + msCG(MdagMOp,Phi,MPhi_k); + + dSdU = zero; + for(int k=0;k