Generalized GeneralEvenOddRatioRationalPseudoFermionAction such that the multi-shift CG algorithm can be overridden by derived classes

Added a mixed-precision variant of GeneralEvenOddRatioRationalPseudoFermionAction and a verification test against double prec class
Fixed non-const reference used in passing RHMC approx to multishift classes

Grid/algorithms/iterative/ConjugateGradientMultiShift.h

@@ -52,7 +52,7 @@ public:
   MultiShiftFunction shifts;
   std::vector<RealD> TrueResidualShift;
 
-  ConjugateGradientMultiShift(Integer maxit,MultiShiftFunction &_shifts) : 
+  ConjugateGradientMultiShift(Integer maxit, const MultiShiftFunction &_shifts) : 
     MaxIterations(maxit),
     shifts(_shifts)
   { 

Grid/algorithms/iterative/ConjugateGradientMultiShiftMixedPrec.h

@@ -93,7 +93,7 @@ public:
   GridBase* SinglePrecGrid; //Grid for single-precision fields
   LinearOperatorBase<FieldF> &Linop_f; //single precision
 
-  ConjugateGradientMultiShiftMixedPrec(Integer maxit, MultiShiftFunction &_shifts,
+  ConjugateGradientMultiShiftMixedPrec(Integer maxit, const MultiShiftFunction &_shifts,
 				       GridBase* _SinglePrecGrid, LinearOperatorBase<FieldF> &_Linop_f,
 				       int _ReliableUpdateFreq
 				       ) : 

Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatio.h

@@ -79,6 +79,27 @@ NAMESPACE_BEGIN(Grid);
       FermionField PhiEven; // the pseudo fermion field for this trajectory
       FermionField PhiOdd; // the pseudo fermion field for this trajectory
 
+    protected:
+      static constexpr bool Numerator = true;
+      static constexpr bool Denominator = false;
+
+      //Allow derived classes to override the multishift CG
+      virtual void multiShiftInverse(bool numerator, const MultiShiftFunction &approx, const Integer MaxIter, const FermionField &in, FermionField &out){
+	SchurDifferentiableOperator<Impl> schurOp(numerator ? NumOp : DenOp);
+	ConjugateGradientMultiShift<FermionField> msCG(MaxIter, approx);
+	msCG(schurOp,in, out);
+      }
+      virtual void multiShiftInverse(bool numerator, const MultiShiftFunction &approx, const Integer MaxIter, const FermionField &in, std::vector<FermionField> &out_elems, FermionField &out){
+	SchurDifferentiableOperator<Impl> schurOp(numerator ? NumOp : DenOp);
+	ConjugateGradientMultiShift<FermionField> msCG(MaxIter, approx);
+	msCG(schurOp,in, out_elems, out);
+      }
+      //Allow derived classes to override the gauge import
+      virtual void ImportGauge(const GaugeField &U){
+	NumOp.ImportGauge(U);
+	DenOp.ImportGauge(U);
+      }
+      
     public:
 
       GeneralEvenOddRatioRationalPseudoFermionAction(FermionOperator<Impl>  &_NumOp, 
@@ -188,22 +209,16 @@ NAMESPACE_BEGIN(Grid);
 	pickCheckerboard(Even,etaEven,eta);
 	pickCheckerboard(Odd,etaOdd,eta);
 
-	NumOp.ImportGauge(U);
+	ImportGauge(U);
-	DenOp.ImportGauge(U);
-
 
 	// MdagM^1/(2*inv_pow) eta
 	std::cout<<GridLogMessage << action_name() << " refresh: doing (M^dag M)^{1/" << 2*param.inv_pow << "} eta" << std::endl;
-	SchurDifferentiableOperator<Impl> MdagM(DenOp);
+	multiShiftInverse(Denominator, ApproxHalfPowerAction, param.MaxIter, etaOdd, tmp);
-	ConjugateGradientMultiShift<FermionField> msCG_M(param.MaxIter,ApproxHalfPowerAction);
-	msCG_M(MdagM,etaOdd,tmp);
 
 	// VdagV^-1/(2*inv_pow) MdagM^1/(2*inv_pow) eta
 	std::cout<<GridLogMessage << action_name() << " refresh: doing (V^dag V)^{-1/" << 2*param.inv_pow << "} ( (M^dag M)^{1/" << 2*param.inv_pow << "} eta)" << std::endl;
-	SchurDifferentiableOperator<Impl> VdagV(NumOp);
+	multiShiftInverse(Numerator, ApproxNegHalfPowerAction, param.MaxIter, tmp, PhiOdd);
-	ConjugateGradientMultiShift<FermionField> msCG_V(param.MaxIter,ApproxNegHalfPowerAction);
+		
-	msCG_V(VdagV,tmp,PhiOdd);
-
 	assert(NumOp.ConstEE() == 1);
 	assert(DenOp.ConstEE() == 1);
 	PhiEven = Zero();
@@ -215,29 +230,25 @@ NAMESPACE_BEGIN(Grid);
       //////////////////////////////////////////////////////
       virtual RealD S(const GaugeField &U) {
 	std::cout<<GridLogMessage << action_name() << " compute action: starting" << std::endl;
-	NumOp.ImportGauge(U);
+	ImportGauge(U);
-	DenOp.ImportGauge(U);
 
 	FermionField X(NumOp.FermionRedBlackGrid());
 	FermionField Y(NumOp.FermionRedBlackGrid());
 
 	// VdagV^1/(2*inv_pow) Phi
 	std::cout<<GridLogMessage << action_name() << " compute action: doing (V^dag V)^{1/" << 2*param.inv_pow << "} Phi" << std::endl;
-	SchurDifferentiableOperator<Impl> VdagV(NumOp);
+	multiShiftInverse(Numerator, ApproxHalfPowerAction, param.MaxIter, PhiOdd,X);
-	ConjugateGradientMultiShift<FermionField> msCG_V(param.MaxIter,ApproxHalfPowerAction);
-	msCG_V(VdagV,PhiOdd,X);
 
 	// MdagM^-1/(2*inv_pow) VdagV^1/(2*inv_pow) Phi
 	std::cout<<GridLogMessage << action_name() << " compute action: doing (M^dag M)^{-1/" << 2*param.inv_pow << "} ( (V^dag V)^{1/" << 2*param.inv_pow << "} Phi)" << std::endl;
-	SchurDifferentiableOperator<Impl> MdagM(DenOp);
+	multiShiftInverse(Denominator, ApproxNegHalfPowerAction, param.MaxIter, X,Y);
-	ConjugateGradientMultiShift<FermionField> msCG_M(param.MaxIter,ApproxNegHalfPowerAction);
-	msCG_M(MdagM,X,Y);
 
 	// Randomly apply rational bounds checks.
 	if ( param.BoundsCheckFreq != 0 && (rand()%param.BoundsCheckFreq)==0 ) { 
 	  std::cout<<GridLogMessage << action_name() << " compute action: doing bounds check" << std::endl;
 	  FermionField gauss(NumOp.FermionRedBlackGrid());
 	  gauss = PhiOdd;
+	  SchurDifferentiableOperator<Impl> MdagM(DenOp);
 	  HighBoundCheck(MdagM,gauss,param.hi);
 	  InversePowerBoundsCheck(param.inv_pow,param.MaxIter,param.action_tolerance*100,MdagM,gauss,ApproxNegPowerAction);
 	}
@@ -295,21 +306,21 @@ NAMESPACE_BEGIN(Grid);
 
 	GaugeField   tmp(NumOp.GaugeGrid());
 
-	NumOp.ImportGauge(U);
+	ImportGauge(U);
-	DenOp.ImportGauge(U);
-
-	SchurDifferentiableOperator<Impl> VdagV(NumOp);
-	SchurDifferentiableOperator<Impl> MdagM(DenOp);
-
-	ConjugateGradientMultiShift<FermionField> msCG_V(param.MaxIter,ApproxHalfPowerMD);
-	ConjugateGradientMultiShift<FermionField> msCG_M(param.MaxIter,ApproxNegPowerMD);
 
 	std::cout<<GridLogMessage << action_name() << " deriv: doing (V^dag V)^{1/" << 2*param.inv_pow << "} Phi" << std::endl;
-	msCG_V(VdagV,PhiOdd,MpvPhi_k,MpvPhi);
+	multiShiftInverse(Numerator, ApproxHalfPowerMD, param.MaxIter, PhiOdd,MpvPhi_k,MpvPhi);
+
 	std::cout<<GridLogMessage << action_name() << " deriv: doing (M^dag M)^{-1/" << param.inv_pow << "} ( (V^dag V)^{1/" << 2*param.inv_pow << "} Phi)" << std::endl;
-	msCG_M(MdagM,MpvPhi,MfMpvPhi_k,MfMpvPhi);
+	multiShiftInverse(Denominator, ApproxNegPowerMD, param.MaxIter, MpvPhi,MfMpvPhi_k,MfMpvPhi);
+
 	std::cout<<GridLogMessage << action_name() << " deriv: doing (V^dag V)^{1/" << 2*param.inv_pow << "} ( (M^dag M)^{-1/" << param.inv_pow << "} (V^dag V)^{1/" << 2*param.inv_pow << "} Phi)" << std::endl;
-	msCG_V(VdagV,MfMpvPhi,MpvMfMpvPhi_k,MpvMfMpvPhi);
+	multiShiftInverse(Numerator, ApproxHalfPowerMD, param.MaxIter, MfMpvPhi,MpvMfMpvPhi_k,MpvMfMpvPhi);
+		
+
+	SchurDifferentiableOperator<Impl> MdagM(DenOp);
+	SchurDifferentiableOperator<Impl> VdagV(NumOp);
+
 
 	RealD ak;
 

Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatioMixedPrec.h

@@ -0,0 +1,93 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/qcd/action/pseudofermion/GeneralEvenOddRationalRatioMixedPrec.h
+
+    Copyright (C) 2015
+
+    Author: Christopher Kelly <ckelly@bnl.gov>
+    Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+    /*  END LEGAL */
+#ifndef QCD_PSEUDOFERMION_GENERAL_EVEN_ODD_RATIONAL_RATIO_MIXED_PREC_H
+#define QCD_PSEUDOFERMION_GENERAL_EVEN_ODD_RATIONAL_RATIO_MIXED_PREC_H
+
+NAMESPACE_BEGIN(Grid);
+
+    /////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Generic rational approximation for ratios of operators utilizing the mixed precision multishift algorithm
+    // cf. GeneralEvenOddRational.h for details
+    /////////////////////////////////////////////////////////////////////////////////////////////////////////////
+      
+    template<class ImplD, class ImplF>
+    class GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction : public GeneralEvenOddRatioRationalPseudoFermionAction<ImplD> {
+    private:
+      typedef typename ImplD::FermionField FermionFieldD;
+      typedef typename ImplF::FermionField FermionFieldF;
+
+      FermionOperator<ImplD> & NumOpD;
+      FermionOperator<ImplD> & DenOpD;
+     
+      FermionOperator<ImplF> & NumOpF;
+      FermionOperator<ImplF> & DenOpF;
+
+      Integer ReliableUpdateFreq;
+    protected:
+
+      //Allow derived classes to override the multishift CG
+      virtual void multiShiftInverse(bool numerator, const MultiShiftFunction &approx, const Integer MaxIter, const FermionFieldD &in, FermionFieldD &out){
+	SchurDifferentiableOperator<ImplD> schurOpD(numerator ? NumOpD : DenOpD);
+	SchurDifferentiableOperator<ImplF> schurOpF(numerator ? NumOpF : DenOpF);
+
+	ConjugateGradientMultiShiftMixedPrec<FermionFieldD, FermionFieldF> msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq);
+	msCG(schurOpD, in, out);
+      }
+      virtual void multiShiftInverse(bool numerator, const MultiShiftFunction &approx, const Integer MaxIter, const FermionFieldD &in, std::vector<FermionFieldD> &out_elems, FermionFieldD &out){
+	SchurDifferentiableOperator<ImplD> schurOpD(numerator ? NumOpD : DenOpD);
+	SchurDifferentiableOperator<ImplF> schurOpF(numerator ? NumOpF : DenOpF);
+
+	ConjugateGradientMultiShiftMixedPrec<FermionFieldD, FermionFieldF> msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq);
+	msCG(schurOpD, in, out_elems, out);
+      }
+      //Allow derived classes to override the gauge import
+      virtual void ImportGauge(const typename ImplD::GaugeField &Ud){
+	typename ImplF::GaugeField Uf(NumOpF.GaugeGrid());
+	precisionChange(Uf, Ud);
+	
+	NumOpD.ImportGauge(Ud);
+	DenOpD.ImportGauge(Ud);
+
+	NumOpF.ImportGauge(Uf);
+	DenOpF.ImportGauge(Uf);
+      }
+      
+    public:
+      GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction(FermionOperator<ImplD>  &_NumOpD, FermionOperator<ImplD>  &_DenOpD, 
+							      FermionOperator<ImplF>  &_NumOpF, FermionOperator<ImplF>  &_DenOpF, 
+							      const RationalActionParams & p, Integer _ReliableUpdateFreq
+							      ) : GeneralEvenOddRatioRationalPseudoFermionAction<ImplD>(_NumOpD, _DenOpD, p),
+								  ReliableUpdateFreq(_ReliableUpdateFreq), NumOpD(_NumOpD), DenOpD(_DenOpD), NumOpF(_NumOpF), DenOpF(_DenOpF){}
+      
+      virtual std::string action_name(){return "GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction";}
+    };
+
+NAMESPACE_END(Grid);
+
+#endif

Grid/qcd/action/pseudofermion/PseudoFermion.h

@@ -41,6 +41,7 @@ directory
 #include <Grid/qcd/action/pseudofermion/OneFlavourRationalRatio.h>
 #include <Grid/qcd/action/pseudofermion/OneFlavourEvenOddRational.h>
 #include <Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatio.h>
+#include <Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatioMixedPrec.h>
 #include <Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h>
 #include <Grid/qcd/action/pseudofermion/ExactOneFlavourRatio.h>
 

tests/hmc/Test_rhmc_EOWilsonRatio_doubleVsMixedPrec.cc

@@ -0,0 +1,119 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/Test_rhmc_EOWilsonRatio_doubleVsMixedPrec.cc
+
+    Copyright (C) 2015
+
+Author: Christopher Kelly <ckelly@bnl.gov>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+    /*  END LEGAL */
+#include <Grid/Grid.h>
+
+//This test ensures the mixed precision RHMC gives the same result as the regular double precision
+int main(int argc, char **argv) {
+  using namespace Grid;
+
+  Grid_init(&argc, &argv);
+  int threads = GridThread::GetThreads();
+  std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
+
+  typedef GenericHMCRunner<MinimumNorm2> HMCWrapper;  // Uses the default minimum norm
+
+  typedef WilsonImplD FermionImplPolicyD;
+  typedef WilsonFermionD FermionActionD;
+  typedef typename FermionActionD::FermionField FermionFieldD;
+
+  typedef WilsonImplF FermionImplPolicyF;
+  typedef WilsonFermionF FermionActionF;
+  typedef typename FermionActionF::FermionField FermionFieldF;
+
+  //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+  HMCWrapper TheHMC;
+  TheHMC.Resources.AddFourDimGrid("gauge");
+
+  RNGModuleParameters RNGpar;
+  RNGpar.serial_seeds = "1 2 3 4 5";
+  RNGpar.parallel_seeds = "6 7 8 9 10";
+  TheHMC.Resources.SetRNGSeeds(RNGpar);
+
+  auto GridPtrD = TheHMC.Resources.GetCartesian();
+  auto GridRBPtrD = TheHMC.Resources.GetRBCartesian();
+
+  GridCartesian* GridPtrF = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexF::Nsimd()), GridDefaultMpi());
+  GridRedBlackCartesian* GridRBPtrF = SpaceTimeGrid::makeFourDimRedBlackGrid(GridPtrF);
+
+  // temporarily need a gauge field
+  LatticeGaugeFieldD Ud(GridPtrD);
+  LatticeGaugeFieldF Uf(GridPtrF);
+
+  Real mass = -0.77;
+  Real pv   = 0.0;
+
+  FermionActionD DenOpD(Ud, *GridPtrD, *GridRBPtrD, mass);
+  FermionActionD NumOpD(Ud, *GridPtrD, *GridRBPtrD, pv);
+
+  FermionActionF DenOpF(Uf, *GridPtrF, *GridRBPtrF, mass);
+  FermionActionF NumOpF(Uf, *GridPtrF, *GridRBPtrF, pv);
+
+  TheHMC.Resources.AddRNGs();
+  PeriodicGimplR::HotConfiguration(TheHMC.Resources.GetParallelRNG(), Ud);
+
+  std::string seed_string = "the_seed";
+
+  //Setup the pseudofermion actions
+  RationalActionParams GenParams;
+  GenParams.inv_pow = 2;
+  GenParams.lo = 1e-2;
+  GenParams.hi = 64.0;
+  GenParams.MaxIter = 1000;
+  GenParams.action_tolerance = GenParams.md_tolerance = 1e-6;
+  GenParams.action_degree = GenParams.md_degree = 6;
+  GenParams.precision = 64;
+  GenParams.BoundsCheckFreq = 20;
+
+  GeneralEvenOddRatioRationalPseudoFermionAction<FermionImplPolicyD> GenD(NumOpD,DenOpD,GenParams);
+  GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicyD, FermionImplPolicyF> GenFD(NumOpD, DenOpD, 
+													NumOpF, DenOpF, 
+													GenParams, 50);
+  TheHMC.Resources.GetParallelRNG().SeedUniqueString(seed_string);
+  GenD.refresh(Ud, TheHMC.Resources.GetParallelRNG());    
+  RealD Sd = GenD.S(Ud);
+  LatticeGaugeField derivD(Ud);
+  GenD.deriv(Ud,derivD);   
+
+  TheHMC.Resources.GetParallelRNG().SeedUniqueString(seed_string);
+  GenFD.refresh(Ud, TheHMC.Resources.GetParallelRNG());    
+  RealD Sfd = GenFD.S(Ud);
+  LatticeGaugeField derivFD(Ud);
+  GenFD.deriv(Ud,derivFD);   
+
+  //Compare
+  std::cout << "Action : " << Sd << " " << Sfd << " reldiff " << (Sd - Sfd)/Sd << std::endl;
+  
+  LatticeGaugeField diff(Ud);
+  axpy(diff, -1.0, derivD, derivFD);
+  std::cout << "Norm of difference in deriv " << sqrt(norm2(diff)) << std::endl;
+
+  Grid_finalize();
+  return 0;
+}
+