From fd6031b005ffd31fdf29409ac4d77de3ed26aba1 Mon Sep 17 00:00:00 2001
From: paboyle <paboyle@ph.ed.ac.uk>
Date: Sun, 14 Jan 2018 22:46:17 +0000
Subject: [PATCH] Namespace

---
 .../action/pseudofermion/OneFlavourRational.h | 263 +++++++++---------
 1 file changed, 131 insertions(+), 132 deletions(-)

diff --git a/lib/qcd/action/pseudofermion/OneFlavourRational.h b/lib/qcd/action/pseudofermion/OneFlavourRational.h
index 078aebf1..ca324e94 100644
--- a/lib/qcd/action/pseudofermion/OneFlavourRational.h
+++ b/lib/qcd/action/pseudofermion/OneFlavourRational.h
@@ -1,4 +1,4 @@
-    /*************************************************************************************
+/*************************************************************************************
 
     Grid physics library, www.github.com/paboyle/Grid 
 
@@ -23,191 +23,190 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     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 */
+*************************************************************************************/
+/*  END LEGAL */
 #ifndef QCD_PSEUDOFERMION_ONE_FLAVOUR_RATIONAL_H
 #define QCD_PSEUDOFERMION_ONE_FLAVOUR_RATIONAL_H
 
-namespace Grid{
-  namespace QCD{
+NAMESPACE_BEGIN(Grid);
 
-    ///////////////////////////////////////
-    // One flavour rational
-    ///////////////////////////////////////
+///////////////////////////////////////
+// One flavour rational
+///////////////////////////////////////
 
-    // S_f = chi^dag *  N(M^dag*M)/D(M^dag*M) * chi
-    //
-    // Here, M is some operator 
-    // N and D makeup the rat. poly 
-    //
+// S_f = chi^dag *  N(M^dag*M)/D(M^dag*M) * chi
+//
+// Here, M is some operator 
+// N and D makeup the rat. poly 
+//
   
-    template<class Impl>
-    class OneFlavourRationalPseudoFermionAction : public Action<typename Impl::GaugeField> {
-    public:
-      INHERIT_IMPL_TYPES(Impl);
+template<class Impl>
+class OneFlavourRationalPseudoFermionAction : public Action<typename Impl::GaugeField> {
+public:
+  INHERIT_IMPL_TYPES(Impl);
 
-      typedef OneFlavourRationalParams Params;
-      Params param;
+  typedef OneFlavourRationalParams Params;
+  Params param;
 
-      MultiShiftFunction PowerHalf   ;
-      MultiShiftFunction PowerNegHalf;
-      MultiShiftFunction PowerQuarter;
-      MultiShiftFunction PowerNegQuarter;
+  MultiShiftFunction PowerHalf   ;
+  MultiShiftFunction PowerNegHalf;
+  MultiShiftFunction PowerQuarter;
+  MultiShiftFunction PowerNegQuarter;
 
-    private:
+private:
      
-      FermionOperator<Impl> & FermOp;// the basic operator
+  FermionOperator<Impl> & FermOp;// the basic operator
 
-      // NOT using "Nroots"; IroIro is -- perhaps later, but this wasn't good for us historically
-      // and hasenbusch works better
+  // NOT using "Nroots"; IroIro is -- perhaps later, but this wasn't good for us historically
+  // and hasenbusch works better
 
-      FermionField Phi; // the pseudo fermion field for this trajectory
+  FermionField Phi; // the pseudo fermion field for this trajectory
 
-    public:
+public:
 
-      OneFlavourRationalPseudoFermionAction(FermionOperator<Impl>  &Op, 
-					    Params & p
-					    ) : FermOp(Op), Phi(Op.FermionGrid()), param(p) 
-      {
-	AlgRemez remez(param.lo,param.hi,param.precision);
+  OneFlavourRationalPseudoFermionAction(FermionOperator<Impl>  &Op, 
+					Params & p
+					) : FermOp(Op), Phi(Op.FermionGrid()), param(p) 
+  {
+    AlgRemez remez(param.lo,param.hi,param.precision);
 
-	// MdagM^(+- 1/2)
-	std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/2)"<<std::endl;
-	remez.generateApprox(param.degree,1,2);
-	PowerHalf.Init(remez,param.tolerance,false);
-	PowerNegHalf.Init(remez,param.tolerance,true);
+    // MdagM^(+- 1/2)
+    std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/2)"<<std::endl;
+    remez.generateApprox(param.degree,1,2);
+    PowerHalf.Init(remez,param.tolerance,false);
+    PowerNegHalf.Init(remez,param.tolerance,true);
 
-	// MdagM^(+- 1/4)
-	std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/4)"<<std::endl;
-	remez.generateApprox(param.degree,1,4);
-   	PowerQuarter.Init(remez,param.tolerance,false);
-	PowerNegQuarter.Init(remez,param.tolerance,true);
-      };
+    // MdagM^(+- 1/4)
+    std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/4)"<<std::endl;
+    remez.generateApprox(param.degree,1,4);
+    PowerQuarter.Init(remez,param.tolerance,false);
+    PowerNegQuarter.Init(remez,param.tolerance,true);
+  };
 
-      virtual std::string action_name(){return "OneFlavourRationalPseudoFermionAction";}
+  virtual std::string action_name(){return "OneFlavourRationalPseudoFermionAction";}
 
-      virtual std::string LogParameters(){
-	std::stringstream sstream;
-	sstream << GridLogMessage << "["<<action_name()<<"] Low            :" << param.lo <<  std::endl;
-	sstream << GridLogMessage << "["<<action_name()<<"] High           :" << param.hi <<  std::endl;
-	sstream << GridLogMessage << "["<<action_name()<<"] Max iterations :" << param.MaxIter <<  std::endl;
-	sstream << GridLogMessage << "["<<action_name()<<"] Tolerance      :" << param.tolerance <<  std::endl;
-	sstream << GridLogMessage << "["<<action_name()<<"] Degree         :" << param.degree <<  std::endl;
-	sstream << GridLogMessage << "["<<action_name()<<"] Precision      :" << param.precision <<  std::endl;
-	return sstream.str();
-      }  
+  virtual std::string LogParameters(){
+    std::stringstream sstream;
+    sstream << GridLogMessage << "["<<action_name()<<"] Low            :" << param.lo <<  std::endl;
+    sstream << GridLogMessage << "["<<action_name()<<"] High           :" << param.hi <<  std::endl;
+    sstream << GridLogMessage << "["<<action_name()<<"] Max iterations :" << param.MaxIter <<  std::endl;
+    sstream << GridLogMessage << "["<<action_name()<<"] Tolerance      :" << param.tolerance <<  std::endl;
+    sstream << GridLogMessage << "["<<action_name()<<"] Degree         :" << param.degree <<  std::endl;
+    sstream << GridLogMessage << "["<<action_name()<<"] Precision      :" << param.precision <<  std::endl;
+    return sstream.str();
+  }  
 
 
       
-      virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {
+  virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {
 
 	
-	// P(phi) = e^{- phi^dag (MdagM)^-1/2 phi}
-	//        = e^{- phi^dag (MdagM)^-1/4 (MdagM)^-1/4 phi}
-	// Phi = Mdag^{1/4} eta 
-	// P(eta) = e^{- eta^dag eta}
-	//
-	// e^{x^2/2 sig^2} => sig^2 = 0.5.
-	// 
-	// So eta should be of width sig = 1/sqrt(2).
+    // P(phi) = e^{- phi^dag (MdagM)^-1/2 phi}
+    //        = e^{- phi^dag (MdagM)^-1/4 (MdagM)^-1/4 phi}
+    // Phi = Mdag^{1/4} eta 
+    // P(eta) = e^{- eta^dag eta}
+    //
+    // e^{x^2/2 sig^2} => sig^2 = 0.5.
+    // 
+    // So eta should be of width sig = 1/sqrt(2).
 
-	RealD scale = std::sqrt(0.5);
+    RealD scale = std::sqrt(0.5);
 
-	FermionField eta(FermOp.FermionGrid());
+    FermionField eta(FermOp.FermionGrid());
 
-	gaussian(pRNG,eta);
+    gaussian(pRNG,eta);
 
-	FermOp.ImportGauge(U);
+    FermOp.ImportGauge(U);
 
-	// mutishift CG
-	MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(FermOp);
-	ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,PowerQuarter);
-	msCG(MdagMOp,eta,Phi);
+    // mutishift CG
+    MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(FermOp);
+    ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,PowerQuarter);
+    msCG(MdagMOp,eta,Phi);
 
-	Phi=Phi*scale;
+    Phi=Phi*scale;
 	
-      };
+  };
 
-      //////////////////////////////////////////////////////
-      // S = phi^dag (Mdag M)^-1/2 phi
-      //////////////////////////////////////////////////////
-      virtual RealD S(const GaugeField &U) {
+  //////////////////////////////////////////////////////
+  // S = phi^dag (Mdag M)^-1/2 phi
+  //////////////////////////////////////////////////////
+  virtual RealD S(const GaugeField &U) {
 
-	FermOp.ImportGauge(U);
+    FermOp.ImportGauge(U);
 
-	FermionField Y(FermOp.FermionGrid());
+    FermionField Y(FermOp.FermionGrid());
 	
-	MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(FermOp);
+    MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(FermOp);
 
-	ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,PowerNegQuarter);
+    ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,PowerNegQuarter);
 
-	msCG(MdagMOp,Phi,Y);
+    msCG(MdagMOp,Phi,Y);
 
-	RealD action = norm2(Y);
-	std::cout << GridLogMessage << "Pseudofermion action FIXME -- is -1/4 solve or -1/2 solve faster??? "<<action<<std::endl;
-	return action;
-      };
+    RealD action = norm2(Y);
+    std::cout << GridLogMessage << "Pseudofermion action FIXME -- is -1/4 solve or -1/2 solve faster??? "<<action<<std::endl;
+    return action;
+  };
 
-      //////////////////////////////////////////////////////
-      // Need
-      // dS_f/dU = chi^dag   d[N/D]  chi
-      //
-      // N/D is expressed as partial fraction expansion:
-      //
-      //           a0 + \sum_k ak/(M^dagM + bk)
-      //
-      // d[N/D] is then
-      //
-      //          \sum_k -ak [M^dagM+bk]^{-1}  [ dM^dag M + M^dag dM ] [M^dag M + bk]^{-1}
-      //
-      // Need
-      //       Mf Phi_k = [MdagM+bk]^{-1} Phi
-      //       Mf Phi   = \sum_k ak [MdagM+bk]^{-1} Phi
-      //
-      // With these building blocks
-      //
-      //       dS/dU =  \sum_k -ak Mf Phi_k^dag      [ dM^dag M + M^dag dM ] Mf Phi_k
-      //        S    = innerprodReal(Phi,Mf Phi);
-      //////////////////////////////////////////////////////
-      virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
+  //////////////////////////////////////////////////////
+  // Need
+  // dS_f/dU = chi^dag   d[N/D]  chi
+  //
+  // N/D is expressed as partial fraction expansion:
+  //
+  //           a0 + \sum_k ak/(M^dagM + bk)
+  //
+  // d[N/D] is then
+  //
+  //          \sum_k -ak [M^dagM+bk]^{-1}  [ dM^dag M + M^dag dM ] [M^dag M + bk]^{-1}
+  //
+  // Need
+  //       Mf Phi_k = [MdagM+bk]^{-1} Phi
+  //       Mf Phi   = \sum_k ak [MdagM+bk]^{-1} Phi
+  //
+  // With these building blocks
+  //
+  //       dS/dU =  \sum_k -ak Mf Phi_k^dag      [ dM^dag M + M^dag dM ] Mf Phi_k
+  //        S    = innerprodReal(Phi,Mf Phi);
+  //////////////////////////////////////////////////////
+  virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
 
-	const int Npole = PowerNegHalf.poles.size();
+    const int Npole = PowerNegHalf.poles.size();
 
-	std::vector<FermionField> MPhi_k (Npole,FermOp.FermionGrid());
+    std::vector<FermionField> MPhi_k (Npole,FermOp.FermionGrid());
 
-	FermionField X(FermOp.FermionGrid());
-	FermionField Y(FermOp.FermionGrid());
+    FermionField X(FermOp.FermionGrid());
+    FermionField Y(FermOp.FermionGrid());
 
-	GaugeField   tmp(FermOp.GaugeGrid());
+    GaugeField   tmp(FermOp.GaugeGrid());
 
-	FermOp.ImportGauge(U);
+    FermOp.ImportGauge(U);
 
-	MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(FermOp);
+    MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(FermOp);
 
-	ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,PowerNegHalf);
+    ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,PowerNegHalf);
 
-	msCG(MdagMOp,Phi,MPhi_k);
+    msCG(MdagMOp,Phi,MPhi_k);
 
-	dSdU = zero;
-	for(int k=0;k<Npole;k++){
+    dSdU = zero;
+    for(int k=0;k<Npole;k++){
 
-	  RealD ak = PowerNegHalf.residues[k];
+      RealD ak = PowerNegHalf.residues[k];
 
-	  X  = MPhi_k[k];
+      X  = MPhi_k[k];
 
-	  FermOp.M(X,Y);
+      FermOp.M(X,Y);
 
-	  FermOp.MDeriv(tmp , Y, X,DaggerNo );  dSdU=dSdU+ak*tmp;
-	  FermOp.MDeriv(tmp , X, Y,DaggerYes);  dSdU=dSdU+ak*tmp;
+      FermOp.MDeriv(tmp , Y, X,DaggerNo );  dSdU=dSdU+ak*tmp;
+      FermOp.MDeriv(tmp , X, Y,DaggerYes);  dSdU=dSdU+ak*tmp;
 
-	}
+    }
 
-	//dSdU = Ta(dSdU);
+    //dSdU = Ta(dSdU);
 
-      };
-    };
-  }
-}
+  };
+};
+
+NAMESPACE_END(Grid);
 
 
 #endif