portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 07:55:35 +00:00
Code Releases Activity

Namespace

Browse Source
This commit is contained in:
paboyle 2018-01-14 22:46:17 +00:00
parent fe44fc50d9
commit fd6031b005
1 changed files with 131 additions and 132 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

333
lib/qcd/action/pseudofermion/OneFlavourRational.h
View File

@ -1,4 +1,4 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid 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. 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory 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 #ifndef QCD_PSEUDOFERMION_ONE_FLAVOUR_RATIONAL_H
#define QCD_PSEUDOFERMION_ONE_FLAVOUR_RATIONAL_H #define QCD_PSEUDOFERMION_ONE_FLAVOUR_RATIONAL_H
namespace Grid{ NAMESPACE_BEGIN(Grid);
namespace QCD{
/////////////////////////////////////// ///////////////////////////////////////
// One flavour rational // One flavour rational
/////////////////////////////////////// ///////////////////////////////////////
// S_f = chi^dag * N(M^dag*M)/D(M^dag*M) * chi // S_f = chi^dag * N(M^dag*M)/D(M^dag*M) * chi
// //
// Here, M is some operator // Here, M is some operator
// N and D makeup the rat. poly // N and D makeup the rat. poly
// //
template<class Impl> template<class Impl>
class OneFlavourRationalPseudoFermionAction : public Action<typename Impl::GaugeField> { class OneFlavourRationalPseudoFermionAction : public Action<typename Impl::GaugeField> {
public: public:
INHERIT_IMPL_TYPES(Impl); INHERIT_IMPL_TYPES(Impl);
typedef OneFlavourRationalParams Params; typedef OneFlavourRationalParams Params;
Params param; Params param;
MultiShiftFunction PowerHalf ; MultiShiftFunction PowerHalf ;
MultiShiftFunction PowerNegHalf; MultiShiftFunction PowerNegHalf;
MultiShiftFunction PowerQuarter; MultiShiftFunction PowerQuarter;
MultiShiftFunction PowerNegQuarter; 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 // NOT using "Nroots"; IroIro is -- perhaps later, but this wasn't good for us historically
// and hasenbusch works better // 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, OneFlavourRationalPseudoFermionAction(FermionOperator<Impl> &Op,
Params & p Params & p
) : FermOp(Op), Phi(Op.FermionGrid()), param(p) ) : FermOp(Op), Phi(Op.FermionGrid()), param(p)
{ {
AlgRemez remez(param.lo,param.hi,param.precision); AlgRemez remez(param.lo,param.hi,param.precision);
// MdagM^(+- 1/2) // MdagM^(+- 1/2)
std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/2)"<<std::endl; std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/2)"<<std::endl;
remez.generateApprox(param.degree,1,2); remez.generateApprox(param.degree,1,2);
PowerHalf.Init(remez,param.tolerance,false); PowerHalf.Init(remez,param.tolerance,false);
PowerNegHalf.Init(remez,param.tolerance,true); PowerNegHalf.Init(remez,param.tolerance,true);
// MdagM^(+- 1/4) // MdagM^(+- 1/4)
std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/4)"<<std::endl; std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/4)"<<std::endl;
remez.generateApprox(param.degree,1,4); remez.generateApprox(param.degree,1,4);
PowerQuarter.Init(remez,param.tolerance,false); PowerQuarter.Init(remez,param.tolerance,false);
PowerNegQuarter.Init(remez,param.tolerance,true); PowerNegQuarter.Init(remez,param.tolerance,true);
}; };
virtual std::string action_name(){return "OneFlavourRationalPseudoFermionAction";} virtual std::string action_name(){return "OneFlavourRationalPseudoFermionAction";}
virtual std::string LogParameters(){ virtual std::string LogParameters(){
std::stringstream sstream; std::stringstream sstream;
sstream << GridLogMessage << "["<<action_name()<<"] Low :" << param.lo << std::endl; sstream << GridLogMessage << "["<<action_name()<<"] Low :" << param.lo << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] High :" << param.hi << 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()<<"] Max iterations :" << param.MaxIter << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] Tolerance :" << param.tolerance << std::endl; sstream << GridLogMessage << "["<<action_name()<<"] Tolerance :" << param.tolerance << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] Degree :" << param.degree << std::endl; sstream << GridLogMessage << "["<<action_name()<<"] Degree :" << param.degree << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] Precision :" << param.precision << std::endl; sstream << GridLogMessage << "["<<action_name()<<"] Precision :" << param.precision << std::endl;
return sstream.str(); return sstream.str();
}
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).
RealD scale = std::sqrt(0.5);
FermionField eta(FermOp.FermionGrid());
gaussian(pRNG,eta);
FermOp.ImportGauge(U);
// mutishift CG
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(FermOp);
ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,PowerQuarter);
msCG(MdagMOp,eta,Phi);
Phi=Phi*scale;
};
//////////////////////////////////////////////////////
// S = phi^dag (Mdag M)^-1/2 phi
//////////////////////////////////////////////////////
virtual RealD S(const GaugeField &U) {
FermOp.ImportGauge(U);
FermionField Y(FermOp.FermionGrid());
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(FermOp);
ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,PowerNegQuarter);
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;
};
//////////////////////////////////////////////////////
// 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();
std::vector<FermionField> MPhi_k (Npole,FermOp.FermionGrid());
FermionField X(FermOp.FermionGrid());
FermionField Y(FermOp.FermionGrid());
GaugeField tmp(FermOp.GaugeGrid());
FermOp.ImportGauge(U);
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(FermOp);
ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,PowerNegHalf);
msCG(MdagMOp,Phi,MPhi_k);
dSdU = zero;
for(int k=0;k<Npole;k++){
RealD ak = PowerNegHalf.residues[k];
X = MPhi_k[k];
FermOp.M(X,Y);
FermOp.MDeriv(tmp , Y, X,DaggerNo ); dSdU=dSdU+ak*tmp;
FermOp.MDeriv(tmp , X, Y,DaggerYes); dSdU=dSdU+ak*tmp;
}
//dSdU = Ta(dSdU);
};
};
} }
}
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).
RealD scale = std::sqrt(0.5);
FermionField eta(FermOp.FermionGrid());
gaussian(pRNG,eta);
FermOp.ImportGauge(U);
// mutishift CG
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(FermOp);
ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,PowerQuarter);
msCG(MdagMOp,eta,Phi);
Phi=Phi*scale;
};
//////////////////////////////////////////////////////
// S = phi^dag (Mdag M)^-1/2 phi
//////////////////////////////////////////////////////
virtual RealD S(const GaugeField &U) {
FermOp.ImportGauge(U);
FermionField Y(FermOp.FermionGrid());
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(FermOp);
ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,PowerNegQuarter);
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;
};
//////////////////////////////////////////////////////
// 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();
std::vector<FermionField> MPhi_k (Npole,FermOp.FermionGrid());
FermionField X(FermOp.FermionGrid());
FermionField Y(FermOp.FermionGrid());
GaugeField tmp(FermOp.GaugeGrid());
FermOp.ImportGauge(U);
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(FermOp);
ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,PowerNegHalf);
msCG(MdagMOp,Phi,MPhi_k);
dSdU = zero;
for(int k=0;k<Npole;k++){
RealD ak = PowerNegHalf.residues[k];
X = MPhi_k[k];
FermOp.M(X,Y);
FermOp.MDeriv(tmp , Y, X,DaggerNo ); dSdU=dSdU+ak*tmp;
FermOp.MDeriv(tmp , X, Y,DaggerYes); dSdU=dSdU+ak*tmp;
}
//dSdU = Ta(dSdU);
};
};
NAMESPACE_END(Grid);
#endif #endif