Completed implementation of Meofa method of ExactOneFlavourRatio pseudofermion action

Added tests to tests/forces/Test_mobius_force_eofa.cc testing that the EOFA heatbath results in Phi = M^{-1/2} eta

Grid/qcd/action/pseudofermion/ExactOneFlavourRatio.h

@@ -224,14 +224,12 @@ NAMESPACE_BEGIN(Grid);
 
       void Meofa(const GaugeField& U,const FermionField &phi, FermionField & Mphi) 
       {
-#if 0
         Lop.ImportGauge(U);
         Rop.ImportGauge(U);
 
         FermionField spProj_Phi(Lop.FermionGrid());
-	FermionField mPhi(Lop.FermionGrid());
         std::vector<FermionField> tmp(2, Lop.FermionGrid());
-	mPhi = phi;
+	Mphi = phi;
 	
         // LH term: S = S - k <\Phi| P_{-} \Omega_{-}^{\dagger} H(mf)^{-1} \Omega_{-} P_{-} |\Phi>
         spProj(Phi, spProj_Phi, -1, Lop.Ls);
@@ -241,10 +239,12 @@ NAMESPACE_BEGIN(Grid);
         SolverL(Lop, tmp[1], tmp[0]);
         Lop.Dtilde(tmp[0], tmp[1]); // We actually solved Cayley preconditioned system: transform back
         Lop.Omega(tmp[1], tmp[0], -1, 1);
-	mPhi = mPhi -  Lop.k * innerProduct(spProj_Phi, tmp[0]).real();
+	spProj(tmp[0], tmp[1], -1, Lop.Ls);
+
+	Mphi = Mphi -  Lop.k * tmp[1];
 
         // RH term: S = S + k <\Phi| P_{+} \Omega_{+}^{\dagger} ( H(mb)
-        //               - \Delta_{+}(mf,mb) P_{+} )^{-1} \Omega_{-} P_{-} |\Phi>
+        //               - \Delta_{+}(mf,mb) P_{+} )^{-1} \Omega_{+} P_{+} |\Phi>
         spProj(Phi, spProj_Phi, 1, Rop.Ls);
         Rop.Omega(spProj_Phi, tmp[0], 1, 0);
         G5R5(tmp[1], tmp[0]);
@@ -252,8 +252,9 @@ NAMESPACE_BEGIN(Grid);
         SolverR(Rop, tmp[1], tmp[0]);
         Rop.Dtilde(tmp[0], tmp[1]);
         Rop.Omega(tmp[1], tmp[0], 1, 1);
-        action += Rop.k * innerProduct(spProj_Phi, tmp[0]).real();
+	spProj(tmp[0], tmp[1], 1, Rop.Ls);
-#endif
+
+        Mphi = Mphi + Rop.k * tmp[1];
       }
 
       // EOFA action: see Eqn. (10) of arXiv:1706.05843
@@ -279,7 +280,7 @@ NAMESPACE_BEGIN(Grid);
         action -= Lop.k * innerProduct(spProj_Phi, tmp[0]).real();
 
         // RH term: S = S + k <\Phi| P_{+} \Omega_{+}^{\dagger} ( H(mb)
-        //               - \Delta_{+}(mf,mb) P_{+} )^{-1} \Omega_{-} P_{-} |\Phi>
+        //               - \Delta_{+}(mf,mb) P_{+} )^{-1} \Omega_{+} P_{+} |\Phi>
         spProj(Phi, spProj_Phi, 1, Rop.Ls);
         Rop.Omega(spProj_Phi, tmp[0], 1, 0);
         G5R5(tmp[1], tmp[0]);

tests/forces/Test_mobius_force_eofa.cc

@@ -89,7 +89,64 @@ int main (int argc, char** argv)
   ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> Meofa(Lop, Rop, CG, CG, CG, CG, CG, Params, false);
 
   GridSerialRNG  sRNG; sRNG.SeedFixedIntegers(seeds4);
+
+  //Check the rational approximation
+  {    
+    RealD scale = std::sqrt(0.5);
+    LatticeFermion eta    (Lop.FermionGrid());
+    gaussian(RNG5,eta); eta = eta * scale;
+
+    Meofa.refresh(U, eta);
+    
+    //Phi = M^{-1/2} eta
+    //M is Hermitian    
+    //(Phi, M Phi) = eta^\dagger  M^{-1/2} M M^{-1/2} eta = eta^\dagger eta
+    LatticeFermion phi = Meofa.getPhi();
+    LatticeFermion Mphi(FGrid);
+    
+    Meofa.Meofa(U, phi, Mphi);
+    std::cout << "Computing inner product" << std::endl;
+    ComplexD inner = innerProduct(phi, Mphi);
+    ComplexD test = inner - norm2(eta);
+    
+    std::cout << "(phi, Mphi) - (eta,eta): " << test << "  expect 0" << std::endl;
+
+    assert(test.real() < 1e-8);
+    assert(test.imag() < 1e-8);
+
+    //Another test is to use heatbath twice to apply M^{-1/2} to Phi then apply M
+    // M  Phi' 
+    //= M M^{-1/2} Phi 
+    //= M M^{-1/2} M^{-1/2} eta 
+    //= eta
+    Meofa.refresh(U, phi);
+    LatticeFermion phi2 = Meofa.getPhi();
+    LatticeFermion test2(FGrid);
+    Meofa.Meofa(U, phi2, test2);
+    test2  = test2 - eta;
+    RealD test2_norm = norm2(test2);
+    std::cout << "|M M^{-1/2} M^{-1/2} eta - eta|^2 = " << test2_norm << " expect 0" << std::endl;
+    assert( test2_norm < 1e-8 );
+  }
+
+
   Meofa.refresh(U, sRNG, RNG5 );
+
+
+    
+
+
+
+
+
+
+
+
+
+
+
+
+
   RealD S = Meofa.S(U); // pdag M p
 
   // get the deriv of phidag M phi with respect to "U"