Debug code

Grid/qcd/action/fermion/implementation/WilsonFermion5DImplementation.h

@@ -196,7 +196,6 @@ void WilsonFermion5D<Impl>::DhopDir(const FermionField &in, FermionField &out,in
   
   uint64_t Nsite = Umu.Grid()->oSites();
   Kernels::DhopDirKernel(Stencil,Umu,Stencil.CommBuf(),Ls,Nsite,in,out,dirdisp,gamma);
-
 };
 template<class Impl>
 void WilsonFermion5D<Impl>::DhopDirAll(const FermionField &in, std::vector<FermionField> &out)
@@ -247,10 +246,14 @@ void WilsonFermion5D<Impl>::DerivInternal(StencilImpl & st,
 
     Kernels::DhopDirKernel(st, U, st.CommBuf(), Ls, Usites, B, Btilde, mu,gamma);
 
+    std::cout << " InsertForce Btilde "<< norm2(Btilde)<<std::endl;
+
     ////////////////////////////
     // spin trace outer product
     ////////////////////////////
     Impl::InsertForce5D(mat, Btilde, Atilde, mu);
+
+    std::cout << " InsertForce "<< norm2(mat)<<std::endl;
   }
 }
 

Grid/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h

@@ -119,13 +119,19 @@ public:
     //  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
+    std::cout << " tmp 1" << norm2(tmp1)<<std::endl;
     this->_Mat.MooeeInvDag(tmp1,tmp2);   // even->even 
+    std::cout << " tmp 1" << norm2(tmp2)<<std::endl;
     this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerYes);
+    std::cout << " ForceO " << norm2(ForceO)<<std::endl;
           
     //  Accumulate X^dag M_oe MeeInv Der_eo Y
     this->_Mat.Meooe   (U,tmp1);    // even->odd -- implicit -0.5 factor to be applied
+    std::cout << " tmp 1" << norm2(tmp1)<<std::endl;
     this->_Mat.MooeeInv(tmp1,tmp2); // even->even 
+    std::cout << " tmp 2" << norm2(tmp2)<<std::endl;
     this->_Mat.MeoDeriv(ForceE,tmp2,V,DaggerYes);
+    std::cout << " ForceE " << norm2(ForceE)<<std::endl;
 
     assert(ForceE.Checkerboard()==Even);
     assert(ForceO.Checkerboard()==Odd);

Grid/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h

@@ -207,20 +207,27 @@ NAMESPACE_BEGIN(Grid);
         //X = (Mdag M)^-1 V^dag phi
         //Y = (Mdag)^-1 V^dag  phi
         Vpc.MpcDag(PhiOdd,Y);          // Y= Vdag phi
+	std::cout << GridLogMessage <<" Y "<<norm2(Y)<<std::endl;
         X=Zero();
         DerivativeSolver(Mpc,Y,X);     // X= (MdagM)^-1 Vdag phi
+	std::cout << GridLogMessage <<" X "<<norm2(X)<<std::endl;
         Mpc.Mpc(X,Y);                  // Y=  Mdag^-1 Vdag phi
+	std::cout << GridLogMessage <<" Y "<<norm2(Y)<<std::endl;
 
         // phi^dag V (Mdag M)^-1 dV^dag  phi
         Vpc.MpcDagDeriv(force , X, PhiOdd );   dSdU = force;
+	std::cout << GridLogMessage <<" deriv "<<norm2(force)<<std::endl;
   
         // phi^dag dV (Mdag M)^-1 V^dag  phi
         Vpc.MpcDeriv(force , PhiOdd, X );      dSdU = dSdU+force;
+	std::cout << GridLogMessage <<" deriv "<<norm2(force)<<std::endl;
 
         //    -    phi^dag V (Mdag M)^-1 Mdag dM   (Mdag M)^-1 V^dag  phi
         //    -    phi^dag V (Mdag M)^-1 dMdag M   (Mdag M)^-1 V^dag  phi
         Mpc.MpcDeriv(force,Y,X);              dSdU = dSdU-force;
+	std::cout << GridLogMessage <<" deriv "<<norm2(force)<<std::endl;
         Mpc.MpcDagDeriv(force,X,Y);           dSdU = dSdU-force;
+	std::cout << GridLogMessage <<" deriv "<<norm2(force)<<std::endl;
 
         // FIXME No force contribution from EvenEven assumed here
         // Needs a fix for clover.

HMC/Mobius2p1f_EOFA_96I_hmc.cc

@@ -299,11 +299,11 @@ int main(int argc, char **argv) {
   // Probably dominates the force - back to EOFA.
   OneFlavourRationalParams SFRp;
   SFRp.lo       = 0.1;
-  SFRp.hi       = 25.0;
+  SFRp.hi       = 30.0;
   SFRp.MaxIter  = 10000;
-  SFRp.tolerance= 1.0e-5;
+  SFRp.tolerance= 1.0e-8;
-  SFRp.mdtolerance= 2.0e-4;
+  SFRp.mdtolerance= 2.0e-6;
-  SFRp.degree   = 8;
+  SFRp.degree   = 10;
   SFRp.precision= 50;
   
   MobiusEOFAFermionD Strange_Op_L (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , strange_mass, strange_mass, pv_mass, 0.0, -1, M5, b, c);

tests/hmc/Test_hmc_Mobius2p1f.cc

@@ -85,7 +85,7 @@ int main(int argc, char **argv) {
   TheHMC.Resources.AddObservable<PlaqObs>();
   //////////////////////////////////////////////
 
-  const int Ls      = 4;
+  const int Ls      = 8;
   Real beta         = 2.13;
   Real light_mass   = 0.01;
   Real strange_mass = 0.04;