Bounds checking

Grid/qcd/action/pseudofermion/OneFlavourEvenOddRational.h

@@ -157,6 +157,13 @@ class OneFlavourEvenOddRationalPseudoFermionAction
 
     msCG(Mpc, PhiOdd, Y);
 
+    if ( (rand()%param.BoundsCheckFreq)==0 ) { 
+      FermionField gauss(FermOp.FermionRedBlackGrid());
+      gauss = PhiOdd;
+      HighBoundCheck(Mpc,gauss,param.hi);
+      InverseSqrtBoundsCheck(param.MaxIter,param.tolerance*100,Mpc,gauss,PowerNegHalf);
+    }
+
     RealD action = norm2(Y);
     std::cout << GridLogMessage << "Pseudofermion action FIXME -- is -1/4 "
                                    "solve or -1/2 solve faster??? "

Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h

@@ -170,6 +170,14 @@ namespace Grid{
 	ConjugateGradientMultiShift<FermionField> msCG_M(param.MaxIter,PowerNegQuarter);
 	msCG_M(MdagM,X,Y);
 
+	// Randomly apply rational bounds checks.
+	if ( (rand()%param.BoundsCheckFreq)==0 ) { 
+	  FermionField gauss(NumOp.FermionRedBlackGrid());
+	  gauss = PhiOdd;
+	  HighBoundCheck(MdagM,gauss,param.hi);
+	  InverseSqrtBoundsCheck(param.MaxIter,param.tolerance*100,MdagM,gauss,PowerNegHalf);
+	}
+
 	//  Phidag VdagV^1/4 MdagM^-1/4  MdagM^-1/4 VdagV^1/4 Phi
 	RealD action = norm2(Y);
 

Grid/qcd/action/pseudofermion/OneFlavourRational.h

@@ -143,6 +143,14 @@ namespace Grid{
 
 	msCG(MdagMOp,Phi,Y);
 
+	if ( (rand()%param.BoundsCheckFreq)==0 ) { 
+	  FermionField gauss(FermOp.FermionGrid());
+	  gauss = Phi;
+	  HighBoundCheck(MdagMOp,gauss,param.hi);
+	  InverseSqrtBoundsCheck(param.MaxIter,param.tolerance*100,MdagMOp,gauss,PowerNegHalf);
+	}
+
+
 	RealD action = norm2(Y);
 	std::cout << GridLogMessage << "Pseudofermion action FIXME -- is -1/4 solve or -1/2 solve faster??? "<<action<<std::endl;
 	return action;

Grid/qcd/action/pseudofermion/OneFlavourRationalRatio.h

@@ -156,6 +156,14 @@ namespace Grid{
 	ConjugateGradientMultiShift<FermionField> msCG_M(param.MaxIter,PowerNegQuarter);
 	msCG_M(MdagM,X,Y);
 
+	// Randomly apply rational bounds checks.
+	if ( (rand()%param.BoundsCheckFreq)==0 ) { 
+	  FermionField gauss(NumOp.FermionGrid());
+	  gauss = Phi;
+	  HighBoundCheck(MdagM,gauss,param.hi);
+	  InverseSqrtBoundsCheck(param.MaxIter,param.tolerance*100,MdagM,gauss,PowerNegHalf);
+	}
+
 	//  Phidag VdagV^1/4 MdagM^-1/4  MdagM^-1/4 VdagV^1/4 Phi
 	RealD action = norm2(Y);
 

Grid/qcd/action/pseudofermion/TwoFlavour.h

@@ -85,21 +85,20 @@ class TwoFlavourPseudoFermionAction : public Action<typename Impl::GaugeField> {
     // and must multiply by 0.707....
     //
     // Chroma has this scale factor: two_flavor_monomial_w.h
+    // CPS uses this factor
     // IroIro: does not use this scale. It is absorbed by a change of vars
     //         in the Phi integral, and thus is only an irrelevant prefactor for
     //         the partition function.
     //
 
-    RealD scale = std::sqrt(0.5);
+    const RealD scale = std::sqrt(0.5);
 
     FermionField eta(FermOp.FermionGrid());
 
-    gaussian(pRNG, eta);
+    gaussian(pRNG, eta); eta = scale *eta;
 
     FermOp.ImportGauge(U);
     FermOp.Mdag(eta, Phi);
-
-    Phi = Phi * scale;
   };
 
   //////////////////////////////////////////////////////