From 48db72259e66d6486d35e1e7ff2df6da74057648 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Tue, 18 Aug 2015 18:37:39 +0100
Subject: [PATCH] EvenOdd schur decomposed mpcdagmpc version of rhmc
 determinant. dH is also small and plaquette looks right.

---
 .../pseudofermion/OneFlavourEvenOddRational.h | 185 ++++++++++++++++++
 tests/Test_rhmc_EOWilson1p1.cc                |  55 ++++++
 2 files changed, 240 insertions(+)
 create mode 100644 lib/qcd/action/pseudofermion/OneFlavourEvenOddRational.h
 create mode 100644 tests/Test_rhmc_EOWilson1p1.cc

diff --git a/lib/qcd/action/pseudofermion/OneFlavourEvenOddRational.h b/lib/qcd/action/pseudofermion/OneFlavourEvenOddRational.h
new file mode 100644
index 00000000..beb3abc1
--- /dev/null
+++ b/lib/qcd/action/pseudofermion/OneFlavourEvenOddRational.h
@@ -0,0 +1,185 @@
+#ifndef QCD_PSEUDOFERMION_ONE_FLAVOUR_EVEN_ODD_RATIONAL_H
+#define QCD_PSEUDOFERMION_ONE_FLAVOUR_EVEN_ODD_RATIONAL_H
+
+namespace Grid{
+  namespace QCD{
+
+    ///////////////////////////////////////
+    // One flavour rational
+    ///////////////////////////////////////
+
+    // S_f = chi^dag *  N(Mpc^dag*Mpc)/D(Mpc^dag*Mpc) * chi
+    //
+    // Here, M is some operator 
+    // N and D makeup the rat. poly 
+    //
+  
+    template<class Impl>
+    class OneFlavourEvenOddRationalPseudoFermionAction : 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 PhiEven; // the pseudo fermion field for this trajectory
+      FermionField PhiOdd; // the pseudo fermion field for this trajectory
+                        
+
+    public:
+
+      OneFlavourEvenOddRationalPseudoFermionAction(FermionOperator<Impl>  &Op, 
+						   Params & p ) : FermOp(Op), 
+	PhiEven(Op.FermionRedBlackGrid()), 
+	PhiOdd (Op.FermionRedBlackGrid()), 
+	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 (MpcdagMpc)^-1/2 phi}
+	//        = e^{- phi^dag (MpcdagMpc)^-1/4 (MpcdagMpc)^-1/4 phi}
+	// Phi = MpcdagMpc^{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());
+	FermionField etaOdd (FermOp.FermionRedBlackGrid());
+	FermionField etaEven(FermOp.FermionRedBlackGrid());
+
+	gaussian(pRNG,eta);	eta=eta*scale;
+
+	pickCheckerboard(Even,etaEven,eta);
+	pickCheckerboard(Odd,etaOdd,eta);
+
+	FermOp.ImportGauge(U);
+
+	// mutishift CG
+	SchurDifferentiableOperator<Impl> Mpc(FermOp);
+	ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,PowerQuarter);
+	msCG(Mpc,etaOdd,PhiOdd);
+
+	//////////////////////////////////////////////////////
+	// FIXME : Clover term not yet..
+	//////////////////////////////////////////////////////
+
+	assert(FermOp.ConstEE() == 1);
+	PhiEven = zero;
+	
+      };
+
+      //////////////////////////////////////////////////////
+      // S = phi^dag (Mdag M)^-1/2 phi
+      //////////////////////////////////////////////////////
+      virtual RealD S(const GaugeField &U) {
+
+	FermOp.ImportGauge(U);
+
+	FermionField Y(FermOp.FermionRedBlackGrid());
+	
+	SchurDifferentiableOperator<Impl> Mpc(FermOp);
+
+	ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,PowerNegQuarter);
+
+	msCG(Mpc,PhiOdd,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.FermionRedBlackGrid());
+
+	FermionField X(FermOp.FermionRedBlackGrid());
+	FermionField Y(FermOp.FermionRedBlackGrid());
+
+	GaugeField   tmp(FermOp.GaugeGrid());
+
+	FermOp.ImportGauge(U);
+
+	SchurDifferentiableOperator<Impl> Mpc(FermOp);
+
+	ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,PowerNegHalf);
+
+	msCG(Mpc,PhiOdd,MPhi_k);
+
+	dSdU = zero;
+	for(int k=0;k<Npole;k++){
+
+	  RealD ak = PowerNegHalf.residues[k];
+
+	  X  = MPhi_k[k];
+
+	  Mpc.Mpc(X,Y);
+	  Mpc.MpcDeriv   (tmp , Y, X );  dSdU=dSdU+ak*tmp;
+	  Mpc.MpcDagDeriv(tmp , X, Y );  dSdU=dSdU+ak*tmp;
+
+	}
+
+	dSdU = Ta(dSdU);
+
+      };
+    };
+  }
+}
+
+
+#endif
diff --git a/tests/Test_rhmc_EOWilson1p1.cc b/tests/Test_rhmc_EOWilson1p1.cc
new file mode 100644
index 00000000..96be8b20
--- /dev/null
+++ b/tests/Test_rhmc_EOWilson1p1.cc
@@ -0,0 +1,55 @@
+#include "Grid.h"
+
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  std::vector<int> latt_size   = GridDefaultLatt();
+  std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
+  std::vector<int> mpi_layout  = GridDefaultMpi();
+  
+  GridCartesian            Fine(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian  RBFine(latt_size,simd_layout,mpi_layout);
+  GridParallelRNG  pRNG(&Fine);
+  pRNG.SeedRandomDevice();
+  LatticeLorentzColourMatrix     U(&Fine);
+
+  SU3::HotConfiguration(pRNG, U);
+
+  // simplify template declaration? Strip the lorentz from the second template
+  WilsonGaugeAction<LatticeLorentzColourMatrix, LatticeColourMatrix> Waction(5.6);
+
+  Real mass=-0.77;
+  WilsonFermionR FermOp(U,Fine,RBFine,mass);
+
+  // 1+1 flavour
+  OneFlavourRationalParams Params(1.0e-4,64.0,1000,1.0e-6);
+  OneFlavourEvenOddRationalPseudoFermionAction<WilsonImplR> WilsonNf1a(FermOp,Params);
+  OneFlavourEvenOddRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(FermOp,Params);
+  
+  //Collect actions
+  ActionLevel Level1;
+  Level1.push_back(&WilsonNf1a);
+  Level1.push_back(&WilsonNf1b);
+  Level1.push_back(&Waction);
+
+  ActionSet FullSet;
+  FullSet.push_back(Level1);
+
+  // Create integrator
+  typedef MinimumNorm2  IntegratorAlgorithm;// change here to change the algorithm
+
+  IntegratorParameters MDpar(12,20,1.0);
+  Integrator<IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
+
+  // Create HMC
+  HMCparameters HMCpar;
+  HybridMonteCarlo<IntegratorAlgorithm>  HMC(HMCpar, MDynamics);
+  HMC.evolve(U);
+
+}