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

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

lib/algorithms/approx/MultiShiftFunction.h

@@ -11,20 +11,27 @@ public:
   std::vector<RealD> 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()),
+  void Init(AlgRemez & remez,double tol,bool inverse) 
-      tolerances(remez.getDegree(),tol),
-      poles(remez.getDegree()),
-      residues(remez.getDegree())
   {
+    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

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<int> 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

lib/qcd/action/Actions.h

@@ -14,6 +14,7 @@
 // Abstract base interface
 ////////////////////////////////////////////
 #include <qcd/action/ActionBase.h>
+#include <qcd/action/ActionParams.h>
 
 ////////////////////////////////////////////
 // Gauge Actions
@@ -157,9 +158,9 @@ typedef DomainWallFermion<GparityWilsonImplD> GparityDomainWallFermionD;
 #include <qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h>
 
 //Todo: RHMC
-//#include <qcd/action/pseudofermion/OneFlavour.h>
+#include <qcd/action/pseudofermion/OneFlavourRational.h>
-//#include <qcd/action/pseudofermion/OneFlavourRatio.h>
+//#include <qcd/action/pseudofermion/OneFlavourRationalRatio.h>
-//#include <qcd/action/pseudofermion/OneFlavourEvenOdd.h>
+//#include <qcd/action/pseudofermion/OneFlavourEvenOddRational.h>
-//#include <qcd/action/pseudofermion/OneFlavourEvenOddRatio.h>
+//#include <qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h>
 
 #endif

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;

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<SiteDoubledGaugeField> DoubledGaugeField;
 
       typedef WilsonCompressor<SiteHalfSpinor,SiteSpinor> Compressor;
-
+      typedef WilsonImplParams ImplParams;
-      typedef struct 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 S,int Nrepresentation>
     class GparityWilsonImpl : public ImplGauge<S,Nrepresentation> { 
     public:
@@ -198,7 +199,7 @@ PARALLEL_FOR_LOOP
 
       typedef WilsonCompressor<SiteHalfSpinor,SiteSpinor> Compressor;
 
-      typedef struct GparityWilsonImplParams {std::vector<int> twists; } ImplParams;
+      typedef GparityWilsonImplParams ImplParams;
       ImplParams Params;
       GparityWilsonImpl(const ImplParams &p= ImplParams()) : Params(p) {}; 
       

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 Impl>
+    class OneFlavourRationalPseudoFermionAction : public Action<typename Impl::GaugeField> {
+    public:
+      INHERIT_IMPL_TYPES(Impl);
+
+      typedef OneFlavourRationalParams Params;
+      Params param;
+
+      MultiShiftFunction PowerHalf   ;
+      MultiShiftFunction PowerNegHalf;
+      MultiShiftFunction PowerQuarter;
+      MultiShiftFunction PowerNegQuarter;
+
+    private:
+     
+      FermionOperator<Impl> & 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<Impl>  &Op, 
+					    Params & p
+					    ) : FermOp(Op), Phi(Op.FermionGrid()), param(p) 
+      {
+	AlgRemez remez(param.lo,param.hi,param.precision);
+
+	// MdagM^(+- 1/2)
+	std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/2)"<<std::endl;
+	remez.generateApprox(param.degree,1,2);
+	PowerHalf.Init(remez,param.tolerance,false);
+	PowerNegHalf.Init(remez,param.tolerance,true);
+
+	// MdagM^(+- 1/4)
+	std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/4)"<<std::endl;
+	remez.generateApprox(param.degree,1,4);
+   	PowerQuarter.Init(remez,param.tolerance,false);
+	PowerNegQuarter.Init(remez,param.tolerance,true);
+      };
+      
+      virtual void init(const GaugeField &U, GridParallelRNG& pRNG) {
+
+	// P(phi) = e^{- phi^dag (MdagM)^-1/2 phi}
+	//        = e^{- phi^dag (MdagM)^-1/4 (MdagM)^-1/4 phi}
+	// Phi = Mdag^{1/4} eta 
+	// P(eta) = e^{- eta^dag eta}
+	//
+	// e^{x^2/2 sig^2} => 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<FermionOperator<Impl> ,FermionField> MdagMOp(FermOp);
+	ConjugateGradientMultiShift<FermionField> 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<FermionOperator<Impl> ,FermionField> MdagMOp(FermOp);
+
+	ConjugateGradientMultiShift<FermionField> 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??? "<<action<<std::endl;
+	return action;
+      };
+
+      //////////////////////////////////////////////////////
+      // Need
+      // dS_f/dU = chi^dag   d[N/D]  chi
+      //
+      // N/D is expressed as partial fraction expansion:
+      //
+      //           a0 + \sum_k ak/(M^dagM + bk)
+      //
+      // d[N/D] is then
+      //
+      //          \sum_k -ak [M^dagM+bk]^{-1}  [ dM^dag M + M^dag dM ] [M^dag M + bk]^{-1}
+      //
+      // Need
+      //       Mf Phi_k = [MdagM+bk]^{-1} Phi
+      //       Mf Phi   = \sum_k ak [MdagM+bk]^{-1} Phi
+      //
+      // With these building blocks
+      //
+      //       dS/dU =  \sum_k -ak Mf Phi_k^dag      [ dM^dag M + M^dag dM ] Mf Phi_k
+      //        S    = innerprodReal(Phi,Mf Phi);
+      //////////////////////////////////////////////////////
+      virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
+
+	const int Npole = PowerNegHalf.poles.size();
+
+	std::vector<FermionField> MPhi_k (Npole,FermOp.FermionGrid());
+
+	FermionField X(FermOp.FermionGrid());
+	FermionField Y(FermOp.FermionGrid());
+
+	GaugeField   tmp(FermOp.GaugeGrid());
+
+	FermOp.ImportGauge(U);
+
+	MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(FermOp);
+
+	ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,PowerNegHalf);
+
+	msCG(MdagMOp,Phi,MPhi_k);
+
+	dSdU = zero;
+	for(int k=0;k<Npole;k++){
+
+	  RealD ak = PowerNegHalf.residues[k];
+
+	  X  = MPhi_k[k];
+
+	  FermOp.M(X,Y);
+
+	  FermOp.MDeriv(tmp , Y, X,DaggerNo );  dSdU=dSdU+ak*tmp;
+	  FermOp.MDeriv(tmp , X, Y,DaggerYes);  dSdU=dSdU+ak*tmp;
+
+	}
+
+	dSdU = Ta(dSdU);
+
+      };
+    };
+  }
+}
+
+
+#endif

lib/qcd/action/pseudofermion/TwoFlavour.h

@@ -4,85 +4,6 @@
 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 
-    //
-    // Need
-    // dS_f/dU = chi^dag   P/Q d[N/D]  P/Q  chi
-    //
-    // Here N/D \sim R_{-1/2} ~ (M^dagM)^{-1/2}
-    //
-    // N/D is expressed as partial fraction expansion:
-    //
-    //           a0 + \sum_k ak/(M^dagM + bk)
-    //
-    // d[N/D] is then
-    //
-    //          \sum_k -ak [M^dagM+bk]^{-1}  [ dM^dag M + M^dag dM ] [M^dag M + bk]^{-1}
-    //
-    // Need
-    //
-    //       Mf Phi_k = [MdagM+bk]^{-1} Phi
-    //       Mf Phi   = \sum_k ak [MdagM+bk]^{-1} Phi
-    //
-    // With these building blocks
-    //
-    //       dS/dU =  \sum_k -ak Mf Phi_k^dag      [ dM^dag M + M^dag dM ] Mf Phi_k
-    //        S    = innerprodReal(Phi,Mf Phi);
-    
-    ///////////////////////////////////////
-    // One flavour rational ratio
-    ///////////////////////////////////////
-
-    // S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi
-    //
-    // Here, M is some 5D operator and V is the Pauli-Villars field
-    // N and D makeup the rat. poly of the M term and P and & makeup the rat.poly of the denom term
-    //
-    // Need
-    // dS_f/dU =  chi^dag d[P/Q]  N/D   P/Q  chi
-    //         +  chi^dag   P/Q d[N/D]  P/Q  chi
-    //         +  chi^dag   P/Q   N/D d[P/Q] chi
-    //
-    // Here P/Q \sim R_{1/4}  ~ (V^dagV)^{1/4}
-    // Here N/D \sim R_{-1/2} ~ (M^dagM)^{-1/2}
-    //
-    // P/Q is expressed as partial fraction expansion:
-    //
-    //           a0 + \sum_k ak/(V^dagV + bk)
-    //
-    // d[P/Q] is then
-    //
-    //          \sum_k -ak [V^dagV+bk]^{-1}  [ dV^dag V + V^dag dV ] [V^dag V + bk]^{-1}
-    //
-    // and similar for N/D.
-    // 
-    // Need
-    //       MpvPhi_k   = [Vdag V + bk]^{-1} chi
-    //
-    //       MpvPhi     = {a0 +  \sum_k ak [Vdag V + bk]^{-1} }chi
-    //
-    //       MfMpvPhi_k = [MdagM+bk]^{-1} MpvPhi
-    //      
-    //       MfMpvPhi   = {a0 +  \sum_k ak [Mdag M + bk]^{-1} } MpvPhi
-    //
-    //       MpvMfMpvPhi_k = [Vdag V + bk]^{-1} MfMpvchi
-    //
-    // With these building blocks
-    //
-    //       dS/dU =  
-    //                 \sum_k -ak MpvPhi_k^dag        [ dV^dag V + V^dag dV ] MpvMfMpvPhi_k           <- deriv on P left
-    //             +   \sum_k -ak MpvMfMpvPhi_k^\dag  [ dV^dag V + V^dag dV ] MpvPhi_k
-    //             +   \sum_k -ak MfMpvPhi_k^dag      [ dM^dag M + M^dag dM ] MfMpvPhi_k
-
-    
     ////////////////////////////////////////////////////////////////////////
     // Two flavour pseudofermion action for any dop
     ////////////////////////////////////////////////////////////////////////

lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h

@@ -95,8 +95,8 @@ namespace Grid{
 
 	// The EE factorised block; normally can replace with zero if det is constant (gauge field indept)
 	// Only really clover term that creates this.
-	//	FermOp.MooeeInvDag(PhiEven,Y);
+	FermOp.MooeeInvDag(PhiEven,Y);
-	//	action = action + norm2(Y);
+	action = action + norm2(Y);
 
 	std::cout << GridLogMessage << "Pseudofermion EO action "<<action<<std::endl;
 	return action;
@@ -135,6 +135,9 @@ namespace Grid{
 	//      FermOp.MooeeInv(Y,X);
 	//	FermOp.MeeDeriv(tmp , Y, X,DaggerNo );    dSdU=tmp;
 	//  FermOp.MeeDeriv(tmp , X, Y,DaggerYes);  dSdU=dSdU+tmp;
+
+	assert(FermOp.ConstEE() == 1);
+
 	/*
         FermOp.MooeeInvDag(PhiOdd,Y);
         FermOp.MooeeInv(Y,X);

lib/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h

@@ -109,9 +109,9 @@ namespace Grid{
 	// Only really clover term that creates this. Leave the EE portion as a future to do to make most
 	// rapid progresss on DWF for now.
 	//
-	// Vpc.MooeeDag(PhiEven,X);
+	NumOp.MooeeDag(PhiEven,X);
-	// Mpc.MooeeInvDag(X,Y);
+	DenOp.MooeeInvDag(X,Y);
-	// action = action + norm2(Y);
+	action = action + norm2(Y);
 
 	return action;
       };
@@ -154,6 +154,11 @@ namespace Grid{
 	Mpc.MpcDeriv(force,Y,X);   dSdU=dSdU-force;
 	Mpc.MpcDagDeriv(force,X,Y);  dSdU=dSdU-force;
 
+	// FIXME No force contribution from EvenEven assumed here
+	// Needs a fix for clover.
+	assert(NumOp.ConstEE() == 1);
+	assert(DenOp.ConstEE() == 1);
+
 	dSdU = -Ta(dSdU);
 
       };