Adding components for even odd decomposed determinant in HMC.

dH not yet conserved, so something wrong in the eo force code still

lib/algorithms/LinearOperator.h

@@ -188,6 +188,7 @@ namespace Grid {
     };
     template<class Matrix,class Field>
       class SchurDiagMooeeOperator :  public SchurOperatorBase<Field> {
+    protected:
       Matrix &_Mat;
     public:
       SchurDiagMooeeOperator (Matrix &Mat): _Mat(Mat){};
@@ -214,6 +215,7 @@ namespace Grid {
     };
     template<class Matrix,class Field>
       class SchurDiagOneOperator :  public SchurOperatorBase<Field> {
+    protected:
       Matrix &_Mat;
     public:
       SchurDiagOneOperator (Matrix &Mat): _Mat(Mat){};
@@ -240,6 +242,7 @@ namespace Grid {
       }
     };
 
+
     /////////////////////////////////////////////////////////////
     // Base classes for functions of operators
     /////////////////////////////////////////////////////////////

lib/qcd/action/Actions.h

@@ -83,6 +83,7 @@
 // Pseudo fermion combinations
 ////////////////////////////////////////
 #include <qcd/action/pseudofermion/TwoFlavour.h>
+#include <qcd/action/pseudofermion/TwoFlavourEvenOdd.h>
 
 
 #endif

lib/qcd/action/pseudofermion/TwoFlavour.h

@@ -106,8 +106,6 @@ namespace Grid{
 
       OperatorFunction<FermionField> &ActionSolver;
 
-      GridBase &Grid;
-
       FermionField Phi; // the pseudo fermion field for this trajectory
 
     public:
@@ -116,9 +114,8 @@ namespace Grid{
       /////////////////////////////////////////////////
     TwoFlavourPseudoFermionAction(FermionOperator<FermionField,GaugeField>  &Op, 
 				  OperatorFunction<FermionField> & DS,
-				  OperatorFunction<FermionField> & AS,
+				  OperatorFunction<FermionField> & AS
-				  GridBase &_Grid
+				  ) : FermOp(Op), DerivativeSolver(DS), ActionSolver(AS), Phi(Op.FermionGrid()) {
-				  ) : FermOp(Op), DerivativeSolver(DS), ActionSolver(AS), Phi(&_Grid), Grid(_Grid) {
       };
       
       //////////////////////////////////////////////////////////////////////////////////////
@@ -140,7 +137,7 @@ namespace Grid{
 	//         in the Phi integral, and thus is only an irrelevant prefactor for the partition function.
 	//
 	RealD scale = std::sqrt(0.5);
-	FermionField eta(&Grid);
+	FermionField eta(FermOp.FermionGrid());
 
 	gaussian(pRNG,eta);
 
@@ -158,8 +155,8 @@ namespace Grid{
 
 	FermOp.ImportGauge(U);
 
-	FermionField X(&Grid);
+	FermionField X(FermOp.FermionGrid());
-	FermionField Y(&Grid);
+	FermionField Y(FermOp.FermionGrid());
 	
 	MdagMLinearOperator<FermionOperator<FermionField,GaugeField> ,FermionField> MdagMOp(FermOp);
 	X=zero;
@@ -182,9 +179,9 @@ namespace Grid{
 
 	FermOp.ImportGauge(U);
 
-	FermionField X(&Grid);
+	FermionField X(FermOp.FermionGrid());
-	FermionField Y(&Grid);
+	FermionField Y(FermOp.FermionGrid());
-	GaugeField   tmp(&Grid);
+	GaugeField   tmp(FermOp.GaugeGrid());
 
 	MdagMLinearOperator<FermionOperator<FermionField,GaugeField> ,FermionField> MdagMOp(FermOp);
 

lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h

@@ -0,0 +1,226 @@
+#ifndef QCD_PSEUDOFERMION_TWO_FLAVOUR_EVEN_ODD_H
+#define QCD_PSEUDOFERMION_TWO_FLAVOUR_EVEN_ODD_H
+
+namespace Grid{
+  namespace QCD{
+
+    template<class Matrix,class FermionField>
+      class SchurDifferentiableOperator :  public SchurDiagMooeeOperator<Matrix,FermionField> 
+      {
+    public:
+	SchurDifferentiableOperator (Matrix &Mat) : SchurDiagMooeeOperator<Matrix,FermionField>(Mat) {};
+
+	void MpcDeriv(LatticeGaugeField &Force,const FermionField &U,const FermionField &V) {
+	
+	  GridBase *fgrid   = this->_Mat.FermionGrid();
+	  GridBase *fcbgrid = this->_Mat.FermionRedBlackGrid();
+	  GridBase *ugrid = this->_Mat.GaugeGrid();
+	  GridBase *ucbgrid = this->_Mat.GaugeRedBlackGrid();
+
+	  Real coeff = 1.0;
+	  FermionField tmp1(fcbgrid);
+	  FermionField tmp2(fcbgrid);
+
+	  conformable(fcbgrid,U._grid);
+	  conformable(fcbgrid,V._grid);
+
+	  // Assert the checkerboard?? or code for either
+	  assert(U.checkerboard==Odd);
+	  assert(V.checkerboard==V.checkerboard);
+
+	  LatticeGaugeField ForceO(ucbgrid);
+	  LatticeGaugeField ForceE(ucbgrid);
+
+	  //  X^dag Der_oe MeeInv Meo Y
+	  // Use Mooee as nontrivial but gauge field indept
+	  this->_Mat.Meooe   (V,tmp1);      // odd->even -- implicit -0.5 factor to be applied
+	  this->_Mat.MooeeInv(tmp1,tmp2);   // even->even 
+	  this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerNo);
+	  
+	  //  Accumulate X^dag M_oe MeeInv Der_eo Y
+	  this->_Mat.MeooeDag   (U,tmp1);    // even->odd -- implicit -0.5 factor to be applied
+	  this->_Mat.MooeeInvDag(tmp1,tmp2); // even->even 
+	  this->_Mat.MeoDeriv(ForceE,tmp2,V,DaggerNo);
+	  
+	  setCheckerboard(Force,ForceE); 
+	  setCheckerboard(Force,ForceO);
+	  Force=-Force;
+	}
+
+
+	void MpcDagDeriv(LatticeGaugeField &Force,const FermionField &U,const FermionField &V) {
+	
+	  GridBase *fgrid   = this->_Mat.FermionGrid();
+	  GridBase *fcbgrid = this->_Mat.FermionRedBlackGrid();
+	  GridBase *ugrid   = this->_Mat.GaugeGrid();
+	  GridBase *ucbgrid = this->_Mat.GaugeRedBlackGrid();
+
+	  Real coeff = 1.0;
+	  FermionField tmp1(fcbgrid);
+	  FermionField tmp2(fcbgrid);
+
+	  conformable(fcbgrid,U._grid);
+	  conformable(fcbgrid,V._grid);
+
+	  // Assert the checkerboard?? or code for either
+	  assert(V.checkerboard==Odd);
+	  assert(V.checkerboard==V.checkerboard);
+
+	  LatticeGaugeField ForceO(ucbgrid);
+	  LatticeGaugeField ForceE(ucbgrid);
+
+	  //  X^dag Der_oe MeeInv Meo Y
+	  // Use Mooee as nontrivial but gauge field indept
+	  this->_Mat.MeooeDag   (V,tmp1);      // odd->even -- implicit -0.5 factor to be applied
+	  this->_Mat.MooeeInvDag(tmp1,tmp2);   // even->even 
+	  this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerYes);
+	  
+	  //  Accumulate X^dag M_oe MeeInv Der_eo Y
+	  this->_Mat.Meooe   (U,tmp1);    // even->odd -- implicit -0.5 factor to be applied
+	  this->_Mat.MooeeInv(tmp1,tmp2); // even->even 
+	  this->_Mat.MeoDeriv(ForceE,tmp2,V,DaggerYes);
+
+	  setCheckerboard(Force,ForceE); 
+	  setCheckerboard(Force,ForceO);
+	  Force=-Force;
+	}
+
+    };
+
+
+    ////////////////////////////////////////////////////////////////////////
+    // Two flavour pseudofermion action for any EO prec dop
+    ////////////////////////////////////////////////////////////////////////
+    template<class GaugeField,class MatrixField,class FermionField>
+      class TwoFlavourEvenOddPseudoFermionAction : public Action<GaugeField> {
+
+    private:
+      
+      FermionOperator<FermionField,GaugeField> & FermOp;// the basic operator
+
+      OperatorFunction<FermionField> &DerivativeSolver;
+
+      OperatorFunction<FermionField> &ActionSolver;
+
+      FermionField PhiOdd;   // the pseudo fermion field for this trajectory
+      FermionField PhiEven;  // the pseudo fermion field for this trajectory
+
+    public:
+      /////////////////////////////////////////////////
+      // Pass in required objects.
+      /////////////////////////////////////////////////
+      TwoFlavourEvenOddPseudoFermionAction(FermionOperator<FermionField,GaugeField>  &Op, 
+					 OperatorFunction<FermionField> & DS,
+					 OperatorFunction<FermionField> & AS
+					   ) : 
+        FermOp(Op), 
+	DerivativeSolver(DS), 
+	ActionSolver(AS), 
+        PhiEven(Op.FermionRedBlackGrid()),
+	PhiOdd(Op.FermionRedBlackGrid())
+		  {};
+      
+      //////////////////////////////////////////////////////////////////////////////////////
+      // Push the gauge field in to the dops. Assume any BC's and smearing already applied
+      //////////////////////////////////////////////////////////////////////////////////////
+      virtual void init(const GaugeField &U, GridParallelRNG& pRNG) {
+
+	// P(phi) = e^{- phi^dag (MpcdagMpc)^-1 phi}
+	// Phi = McpDag eta 
+	// P(eta) = e^{- eta^dag eta}
+	//
+	// e^{x^2/2 sig^2} => sig^2 = 0.5.
+	RealD scale = std::sqrt(0.5);
+
+	FermionField eta    (FermOp.FermionGrid());
+	FermionField etaOdd (FermOp.FermionRedBlackGrid());
+	FermionField etaEven(FermOp.FermionRedBlackGrid());
+
+	gaussian(pRNG,eta);
+	pickCheckerboard(Even,etaEven,eta);
+	pickCheckerboard(Odd,etaOdd,eta);
+
+	SchurDifferentiableOperator<FermionOperator<FermionField,GaugeField>,FermionField> PCop(FermOp);
+	
+	FermOp.ImportGauge(U);
+
+	PCop.MpcDag(etaOdd,PhiOdd);
+	FermOp.MooeeDag(etaEven,PhiEven);
+
+	PhiOdd =PhiOdd*scale;
+	PhiEven=PhiEven*scale;
+	
+      };
+
+      //////////////////////////////////////////////////////
+      // S = phi^dag (Mdag M)^-1 phi  (odd)
+      //   + phi^dag (Mdag M)^-1 phi  (even)
+      //////////////////////////////////////////////////////
+      virtual RealD S(const GaugeField &U) {
+
+	FermOp.ImportGauge(U);
+
+	FermionField X(FermOp.FermionRedBlackGrid());
+	FermionField Y(FermOp.FermionRedBlackGrid());
+	
+	SchurDifferentiableOperator<FermionOperator<FermionField,GaugeField>,FermionField> PCop(FermOp);
+
+	X=zero;
+	ActionSolver(PCop,PhiOdd,X);
+	PCop.Op(X,Y);
+	RealD action = norm2(Y);
+
+	// 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);
+	action = action + norm2(Y);
+
+	std::cout << GridLogMessage << "Pseudofermion EO action "<<action<<std::endl;
+	return action;
+      };
+
+      //////////////////////////////////////////////////////
+      //
+      // dS/du = - phi^dag  (Mdag M)^-1 [ Mdag dM + dMdag M ]  (Mdag M)^-1 phi
+      //       = - phi^dag M^-1 dM (MdagM)^-1 phi -  phi^dag (MdagM)^-1 dMdag dM (Mdag)^-1 phi 
+      //
+      //       = - Ydag dM X  - Xdag dMdag Y
+      //
+      //////////////////////////////////////////////////////
+      virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
+
+	FermOp.ImportGauge(U);
+
+	FermionField X(FermOp.FermionRedBlackGrid());
+	FermionField Y(FermOp.FermionRedBlackGrid());
+	GaugeField tmp(FermOp.GaugeGrid());
+
+	SchurDifferentiableOperator<FermionOperator<FermionField,GaugeField>,FermionField> PCop(FermOp);
+
+	X=zero;
+	DerivativeSolver(PCop,PhiOdd,X);
+	PCop.Op(X,Y);
+
+	// Our conventions really make this UdSdU; We do not differentiate wrt Udag here.
+	// So must take dSdU - adj(dSdU) and left multiply by mom to get dS/dt.
+
+  	PCop.MpcDeriv(tmp , Y, X );    dSdU=tmp;
+	PCop.MpcDagDeriv(tmp , X, Y);  dSdU=dSdU+tmp;
+
+	// Treat the EE case. (MdagM)^-1 = Minv Minvdag
+	// Deriv defaults to zero.
+        FermOp.MooeeInvDag(PhiOdd,Y);
+        FermOp.MooeeInv(Y,X);
+  	FermOp.MeeDeriv(tmp , Y, X,DaggerNo );    dSdU=tmp;
+	FermOp.MeeDeriv(tmp , X, Y,DaggerYes);  dSdU=dSdU+tmp;
+	
+	dSdU = Ta(dSdU);
+
+      };
+
+    };
+    
+  }
+}
+
+#endif

tests/Make.inc

@@ -1,5 +1,5 @@
 
-bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_cheby Test_contfrac_cg Test_contfrac_even_odd Test_contfrac_force Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_force Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_partfrac_force Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi
+bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_cheby Test_contfrac_cg Test_contfrac_even_odd Test_contfrac_force Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_force Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_hmc_EOWilsonFermionGauge Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_partfrac_force Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi
 
 
 Test_GaugeAction_SOURCES=Test_GaugeAction.cc
@@ -90,6 +90,10 @@ Test_gamma_SOURCES=Test_gamma.cc
 Test_gamma_LDADD=-lGrid
 
 
+Test_hmc_EOWilsonFermionGauge_SOURCES=Test_hmc_EOWilsonFermionGauge.cc
+Test_hmc_EOWilsonFermionGauge_LDADD=-lGrid
+
+
 Test_hmc_WilsonFermionGauge_SOURCES=Test_hmc_WilsonFermionGauge.cc
 Test_hmc_WilsonFermionGauge_LDADD=-lGrid
 

tests/Test_hmc_EOWilsonFermionGauge.cc

@@ -0,0 +1,58 @@
+#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;
+  WilsonFermion FermOp(U,Fine,RBFine,mass);
+  
+  ConjugateGradient<LatticeFermion>  CG(1.0e-8,10000);
+  ConjugateGradient<LatticeFermion>  CGmd(1.0e-6,10000);
+  
+  TwoFlavourEvenOddPseudoFermionAction<LatticeLorentzColourMatrix, LatticeColourMatrix,LatticeFermion> 
+    WilsonNf2(FermOp,CG,CG);
+  
+  //Collect actions
+  ActionLevel Level1;
+  Level1.push_back(&WilsonNf2);
+  ActionLevel Level2(3);
+  Level2.push_back(&Waction);
+  ActionSet FullSet;
+  FullSet.push_back(Level1);
+  FullSet.push_back(Level2);
+
+  // Create integrator
+  typedef MinimumNorm2  IntegratorAlgorithm;// change here to change the algorithm
+  //  typedef LeapFrog  IntegratorAlgorithm;// change here to change the algorithm
+  IntegratorParameters MDpar(12,20,1.0);
+  std::vector<int> rel ={1};
+  Integrator<IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
+
+  // Create HMC
+  HMCparameters HMCpar;
+  HybridMonteCarlo<IntegratorAlgorithm>  HMC(HMCpar, MDynamics);
+
+  HMC.evolve(U);
+
+}

tests/Test_hmc_WilsonFermionGauge.cc

@@ -31,7 +31,7 @@ int main (int argc, char ** argv)
   ConjugateGradient<LatticeFermion>  CGmd(1.0e-6,10000);
   
   TwoFlavourPseudoFermionAction<LatticeLorentzColourMatrix, LatticeColourMatrix,LatticeFermion> 
-    WilsonNf2(FermOp,CG,CG,Fine);
+    WilsonNf2(FermOp,CG,CG);
   
   //Collect actions
   ActionLevel Level1;