mirror of
https://github.com/paboyle/Grid.git
synced 2025-04-25 13:15:55 +01:00
Namespace
This commit is contained in:
paboyle
parent
fd6031b005
commit
3e139b52d3
@ -1,4 +1,4 @@
|
|||||||
/*************************************************************************************
|
/*************************************************************************************
|
||||||
|
|
||||||
Grid physics library, www.github.com/paboyle/Grid
|
Grid physics library, www.github.com/paboyle/Grid
|
||||||
|
|
||||||
@ -23,259 +23,257 @@ 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_EVEN_ODD_RATIONAL_RATIO_H
|
#ifndef QCD_PSEUDOFERMION_ONE_FLAVOUR_EVEN_ODD_RATIONAL_RATIO_H
|
||||||
#define QCD_PSEUDOFERMION_ONE_FLAVOUR_EVEN_ODD_RATIONAL_RATIO_H
|
#define QCD_PSEUDOFERMION_ONE_FLAVOUR_EVEN_ODD_RATIONAL_RATIO_H
|
||||||
|
|
||||||
namespace Grid{
|
NAMESPACE_BEGIN(Grid);
|
||||||
namespace QCD{
|
|
||||||
|
|
||||||
///////////////////////////////////////
|
///////////////////////////////////////
|
||||||
// One flavour rational
|
// One flavour rational
|
||||||
///////////////////////////////////////
|
///////////////////////////////////////
|
||||||
|
|
||||||
|
// S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi
|
||||||
|
//
|
||||||
|
// Here P/Q \sim R_{1/4} ~ (V^dagV)^{1/4}
|
||||||
|
// Here N/D \sim R_{-1/2} ~ (M^dagM)^{-1/2}
|
||||||
|
|
||||||
|
template<class Impl>
|
||||||
|
class OneFlavourEvenOddRatioRationalPseudoFermionAction : public Action<typename Impl::GaugeField> {
|
||||||
|
public:
|
||||||
|
|
||||||
|
INHERIT_IMPL_TYPES(Impl);
|
||||||
|
|
||||||
|
typedef OneFlavourRationalParams Params;
|
||||||
|
Params param;
|
||||||
|
|
||||||
|
MultiShiftFunction PowerHalf ;
|
||||||
|
MultiShiftFunction PowerNegHalf;
|
||||||
|
MultiShiftFunction PowerQuarter;
|
||||||
|
MultiShiftFunction PowerNegQuarter;
|
||||||
|
|
||||||
|
private:
|
||||||
|
|
||||||
|
FermionOperator<Impl> & NumOp;// the basic operator
|
||||||
|
FermionOperator<Impl> & DenOp;// the basic operator
|
||||||
|
FermionField PhiEven; // the pseudo fermion field for this trajectory
|
||||||
|
FermionField PhiOdd; // the pseudo fermion field for this trajectory
|
||||||
|
|
||||||
|
public:
|
||||||
|
|
||||||
|
OneFlavourEvenOddRatioRationalPseudoFermionAction(FermionOperator<Impl> &_NumOp,
|
||||||
|
FermionOperator<Impl> &_DenOp,
|
||||||
|
Params & p
|
||||||
|
) :
|
||||||
|
NumOp(_NumOp),
|
||||||
|
DenOp(_DenOp),
|
||||||
|
PhiOdd (_NumOp.FermionRedBlackGrid()),
|
||||||
|
PhiEven(_NumOp.FermionRedBlackGrid()),
|
||||||
|
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/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 "OneFlavourEvenOddRatioRationalPseudoFermionAction";}
|
||||||
|
|
||||||
|
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) {
|
||||||
|
|
||||||
// S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi
|
// S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi
|
||||||
//
|
//
|
||||||
// Here P/Q \sim R_{1/4} ~ (V^dagV)^{1/4}
|
// P(phi) = e^{- phi^dag (VdagV)^1/4 (MdagM)^-1/2 (VdagV)^1/4 phi}
|
||||||
// Here N/D \sim R_{-1/2} ~ (M^dagM)^{-1/2}
|
// = e^{- phi^dag (VdagV)^1/4 (MdagM)^-1/4 (MdagM)^-1/4 (VdagV)^1/4 phi}
|
||||||
|
//
|
||||||
|
// Phi = (VdagV)^-1/4 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).
|
||||||
|
|
||||||
template<class Impl>
|
RealD scale = std::sqrt(0.5);
|
||||||
class OneFlavourEvenOddRatioRationalPseudoFermionAction : public Action<typename Impl::GaugeField> {
|
|
||||||
public:
|
|
||||||
|
|
||||||
INHERIT_IMPL_TYPES(Impl);
|
FermionField eta(NumOp.FermionGrid());
|
||||||
|
FermionField etaOdd (NumOp.FermionRedBlackGrid());
|
||||||
|
FermionField etaEven(NumOp.FermionRedBlackGrid());
|
||||||
|
FermionField tmp(NumOp.FermionRedBlackGrid());
|
||||||
|
|
||||||
typedef OneFlavourRationalParams Params;
|
gaussian(pRNG,eta); eta=eta*scale;
|
||||||
Params param;
|
|
||||||
|
|
||||||
MultiShiftFunction PowerHalf ;
|
pickCheckerboard(Even,etaEven,eta);
|
||||||
MultiShiftFunction PowerNegHalf;
|
pickCheckerboard(Odd,etaOdd,eta);
|
||||||
MultiShiftFunction PowerQuarter;
|
|
||||||
MultiShiftFunction PowerNegQuarter;
|
|
||||||
|
|
||||||
private:
|
NumOp.ImportGauge(U);
|
||||||
|
DenOp.ImportGauge(U);
|
||||||
FermionOperator<Impl> & NumOp;// the basic operator
|
|
||||||
FermionOperator<Impl> & DenOp;// the basic operator
|
|
||||||
FermionField PhiEven; // the pseudo fermion field for this trajectory
|
|
||||||
FermionField PhiOdd; // the pseudo fermion field for this trajectory
|
|
||||||
|
|
||||||
public:
|
|
||||||
|
|
||||||
OneFlavourEvenOddRatioRationalPseudoFermionAction(FermionOperator<Impl> &_NumOp,
|
|
||||||
FermionOperator<Impl> &_DenOp,
|
|
||||||
Params & p
|
|
||||||
) :
|
|
||||||
NumOp(_NumOp),
|
|
||||||
DenOp(_DenOp),
|
|
||||||
PhiOdd (_NumOp.FermionRedBlackGrid()),
|
|
||||||
PhiEven(_NumOp.FermionRedBlackGrid()),
|
|
||||||
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/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 "OneFlavourEvenOddRatioRationalPseudoFermionAction";}
|
|
||||||
|
|
||||||
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) {
|
// MdagM^1/4 eta
|
||||||
|
SchurDifferentiableOperator<Impl> MdagM(DenOp);
|
||||||
|
ConjugateGradientMultiShift<FermionField> msCG_M(param.MaxIter,PowerQuarter);
|
||||||
|
msCG_M(MdagM,etaOdd,tmp);
|
||||||
|
|
||||||
// S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi
|
// VdagV^-1/4 MdagM^1/4 eta
|
||||||
//
|
SchurDifferentiableOperator<Impl> VdagV(NumOp);
|
||||||
// P(phi) = e^{- phi^dag (VdagV)^1/4 (MdagM)^-1/2 (VdagV)^1/4 phi}
|
ConjugateGradientMultiShift<FermionField> msCG_V(param.MaxIter,PowerNegQuarter);
|
||||||
// = e^{- phi^dag (VdagV)^1/4 (MdagM)^-1/4 (MdagM)^-1/4 (VdagV)^1/4 phi}
|
msCG_V(VdagV,tmp,PhiOdd);
|
||||||
//
|
|
||||||
// Phi = (VdagV)^-1/4 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);
|
assert(NumOp.ConstEE() == 1);
|
||||||
|
assert(DenOp.ConstEE() == 1);
|
||||||
|
PhiEven = zero;
|
||||||
|
|
||||||
FermionField eta(NumOp.FermionGrid());
|
};
|
||||||
FermionField etaOdd (NumOp.FermionRedBlackGrid());
|
|
||||||
FermionField etaEven(NumOp.FermionRedBlackGrid());
|
|
||||||
FermionField tmp(NumOp.FermionRedBlackGrid());
|
|
||||||
|
|
||||||
gaussian(pRNG,eta); eta=eta*scale;
|
//////////////////////////////////////////////////////
|
||||||
|
// S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi
|
||||||
|
//////////////////////////////////////////////////////
|
||||||
|
virtual RealD S(const GaugeField &U) {
|
||||||
|
|
||||||
pickCheckerboard(Even,etaEven,eta);
|
NumOp.ImportGauge(U);
|
||||||
pickCheckerboard(Odd,etaOdd,eta);
|
DenOp.ImportGauge(U);
|
||||||
|
|
||||||
NumOp.ImportGauge(U);
|
FermionField X(NumOp.FermionRedBlackGrid());
|
||||||
DenOp.ImportGauge(U);
|
FermionField Y(NumOp.FermionRedBlackGrid());
|
||||||
|
|
||||||
|
// VdagV^1/4 Phi
|
||||||
|
SchurDifferentiableOperator<Impl> VdagV(NumOp);
|
||||||
|
ConjugateGradientMultiShift<FermionField> msCG_V(param.MaxIter,PowerQuarter);
|
||||||
|
msCG_V(VdagV,PhiOdd,X);
|
||||||
|
|
||||||
// MdagM^1/4 eta
|
// MdagM^-1/4 VdagV^1/4 Phi
|
||||||
SchurDifferentiableOperator<Impl> MdagM(DenOp);
|
SchurDifferentiableOperator<Impl> MdagM(DenOp);
|
||||||
ConjugateGradientMultiShift<FermionField> msCG_M(param.MaxIter,PowerQuarter);
|
ConjugateGradientMultiShift<FermionField> msCG_M(param.MaxIter,PowerNegQuarter);
|
||||||
msCG_M(MdagM,etaOdd,tmp);
|
msCG_M(MdagM,X,Y);
|
||||||
|
|
||||||
// VdagV^-1/4 MdagM^1/4 eta
|
// Phidag VdagV^1/4 MdagM^-1/4 MdagM^-1/4 VdagV^1/4 Phi
|
||||||
SchurDifferentiableOperator<Impl> VdagV(NumOp);
|
RealD action = norm2(Y);
|
||||||
ConjugateGradientMultiShift<FermionField> msCG_V(param.MaxIter,PowerNegQuarter);
|
|
||||||
msCG_V(VdagV,tmp,PhiOdd);
|
|
||||||
|
|
||||||
assert(NumOp.ConstEE() == 1);
|
return action;
|
||||||
assert(DenOp.ConstEE() == 1);
|
};
|
||||||
PhiEven = zero;
|
|
||||||
|
|
||||||
};
|
// S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi
|
||||||
|
//
|
||||||
|
// Here, M is some 5D operator and V is the Pauli-Villars field
|
||||||
|
// N and D makeup the rat. poly of the M term and P and & makeup the rat.poly of the denom term
|
||||||
|
//
|
||||||
|
// Need
|
||||||
|
// dS_f/dU = chi^dag d[P/Q] N/D P/Q chi
|
||||||
|
// + chi^dag P/Q d[N/D] P/Q chi
|
||||||
|
// + chi^dag P/Q N/D d[P/Q] chi
|
||||||
|
//
|
||||||
|
// P/Q is expressed as partial fraction expansion:
|
||||||
|
//
|
||||||
|
// a0 + \sum_k ak/(V^dagV + bk)
|
||||||
|
//
|
||||||
|
// d[P/Q] is then
|
||||||
|
//
|
||||||
|
// \sum_k -ak [V^dagV+bk]^{-1} [ dV^dag V + V^dag dV ] [V^dag V + bk]^{-1}
|
||||||
|
//
|
||||||
|
// and similar for N/D.
|
||||||
|
//
|
||||||
|
// Need
|
||||||
|
// MpvPhi_k = [Vdag V + bk]^{-1} chi
|
||||||
|
// MpvPhi = {a0 + \sum_k ak [Vdag V + bk]^{-1} }chi
|
||||||
|
//
|
||||||
|
// MfMpvPhi_k = [MdagM+bk]^{-1} MpvPhi
|
||||||
|
// MfMpvPhi = {a0 + \sum_k ak [Mdag M + bk]^{-1} } MpvPhi
|
||||||
|
//
|
||||||
|
// MpvMfMpvPhi_k = [Vdag V + bk]^{-1} MfMpvchi
|
||||||
|
//
|
||||||
|
|
||||||
//////////////////////////////////////////////////////
|
virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
|
||||||
// S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi
|
|
||||||
//////////////////////////////////////////////////////
|
|
||||||
virtual RealD S(const GaugeField &U) {
|
|
||||||
|
|
||||||
NumOp.ImportGauge(U);
|
const int n_f = PowerNegHalf.poles.size();
|
||||||
DenOp.ImportGauge(U);
|
const int n_pv = PowerQuarter.poles.size();
|
||||||
|
|
||||||
FermionField X(NumOp.FermionRedBlackGrid());
|
std::vector<FermionField> MpvPhi_k (n_pv,NumOp.FermionRedBlackGrid());
|
||||||
FermionField Y(NumOp.FermionRedBlackGrid());
|
std::vector<FermionField> MpvMfMpvPhi_k(n_pv,NumOp.FermionRedBlackGrid());
|
||||||
|
std::vector<FermionField> MfMpvPhi_k (n_f ,NumOp.FermionRedBlackGrid());
|
||||||
|
|
||||||
// VdagV^1/4 Phi
|
FermionField MpvPhi(NumOp.FermionRedBlackGrid());
|
||||||
SchurDifferentiableOperator<Impl> VdagV(NumOp);
|
FermionField MfMpvPhi(NumOp.FermionRedBlackGrid());
|
||||||
ConjugateGradientMultiShift<FermionField> msCG_V(param.MaxIter,PowerQuarter);
|
FermionField MpvMfMpvPhi(NumOp.FermionRedBlackGrid());
|
||||||
msCG_V(VdagV,PhiOdd,X);
|
FermionField Y(NumOp.FermionRedBlackGrid());
|
||||||
|
|
||||||
// MdagM^-1/4 VdagV^1/4 Phi
|
GaugeField tmp(NumOp.GaugeGrid());
|
||||||
SchurDifferentiableOperator<Impl> MdagM(DenOp);
|
|
||||||
ConjugateGradientMultiShift<FermionField> msCG_M(param.MaxIter,PowerNegQuarter);
|
|
||||||
msCG_M(MdagM,X,Y);
|
|
||||||
|
|
||||||
// Phidag VdagV^1/4 MdagM^-1/4 MdagM^-1/4 VdagV^1/4 Phi
|
NumOp.ImportGauge(U);
|
||||||
RealD action = norm2(Y);
|
DenOp.ImportGauge(U);
|
||||||
|
|
||||||
return action;
|
SchurDifferentiableOperator<Impl> VdagV(NumOp);
|
||||||
};
|
SchurDifferentiableOperator<Impl> MdagM(DenOp);
|
||||||
|
|
||||||
// S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi
|
ConjugateGradientMultiShift<FermionField> msCG_V(param.MaxIter,PowerQuarter);
|
||||||
//
|
ConjugateGradientMultiShift<FermionField> msCG_M(param.MaxIter,PowerNegHalf);
|
||||||
// Here, M is some 5D operator and V is the Pauli-Villars field
|
|
||||||
// N and D makeup the rat. poly of the M term and P and & makeup the rat.poly of the denom term
|
|
||||||
//
|
|
||||||
// Need
|
|
||||||
// dS_f/dU = chi^dag d[P/Q] N/D P/Q chi
|
|
||||||
// + chi^dag P/Q d[N/D] P/Q chi
|
|
||||||
// + chi^dag P/Q N/D d[P/Q] chi
|
|
||||||
//
|
|
||||||
// P/Q is expressed as partial fraction expansion:
|
|
||||||
//
|
|
||||||
// a0 + \sum_k ak/(V^dagV + bk)
|
|
||||||
//
|
|
||||||
// d[P/Q] is then
|
|
||||||
//
|
|
||||||
// \sum_k -ak [V^dagV+bk]^{-1} [ dV^dag V + V^dag dV ] [V^dag V + bk]^{-1}
|
|
||||||
//
|
|
||||||
// and similar for N/D.
|
|
||||||
//
|
|
||||||
// Need
|
|
||||||
// MpvPhi_k = [Vdag V + bk]^{-1} chi
|
|
||||||
// MpvPhi = {a0 + \sum_k ak [Vdag V + bk]^{-1} }chi
|
|
||||||
//
|
|
||||||
// MfMpvPhi_k = [MdagM+bk]^{-1} MpvPhi
|
|
||||||
// MfMpvPhi = {a0 + \sum_k ak [Mdag M + bk]^{-1} } MpvPhi
|
|
||||||
//
|
|
||||||
// MpvMfMpvPhi_k = [Vdag V + bk]^{-1} MfMpvchi
|
|
||||||
//
|
|
||||||
|
|
||||||
virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
|
msCG_V(VdagV,PhiOdd,MpvPhi_k,MpvPhi);
|
||||||
|
msCG_M(MdagM,MpvPhi,MfMpvPhi_k,MfMpvPhi);
|
||||||
|
msCG_V(VdagV,MfMpvPhi,MpvMfMpvPhi_k,MpvMfMpvPhi);
|
||||||
|
|
||||||
const int n_f = PowerNegHalf.poles.size();
|
RealD ak;
|
||||||
const int n_pv = PowerQuarter.poles.size();
|
|
||||||
|
|
||||||
std::vector<FermionField> MpvPhi_k (n_pv,NumOp.FermionRedBlackGrid());
|
dSdU = zero;
|
||||||
std::vector<FermionField> MpvMfMpvPhi_k(n_pv,NumOp.FermionRedBlackGrid());
|
|
||||||
std::vector<FermionField> MfMpvPhi_k (n_f ,NumOp.FermionRedBlackGrid());
|
|
||||||
|
|
||||||
FermionField MpvPhi(NumOp.FermionRedBlackGrid());
|
// With these building blocks
|
||||||
FermionField MfMpvPhi(NumOp.FermionRedBlackGrid());
|
//
|
||||||
FermionField MpvMfMpvPhi(NumOp.FermionRedBlackGrid());
|
// dS/dU =
|
||||||
FermionField Y(NumOp.FermionRedBlackGrid());
|
// \sum_k -ak MfMpvPhi_k^dag [ dM^dag M + M^dag dM ] MfMpvPhi_k (1)
|
||||||
|
// + \sum_k -ak MpvMfMpvPhi_k^\dag [ dV^dag V + V^dag dV ] MpvPhi_k (2)
|
||||||
|
// -ak MpvPhi_k^dag [ dV^dag V + V^dag dV ] MpvMfMpvPhi_k (3)
|
||||||
|
|
||||||
GaugeField tmp(NumOp.GaugeGrid());
|
//(1)
|
||||||
|
for(int k=0;k<n_f;k++){
|
||||||
|
ak = PowerNegHalf.residues[k];
|
||||||
|
MdagM.Mpc(MfMpvPhi_k[k],Y);
|
||||||
|
MdagM.MpcDagDeriv(tmp , MfMpvPhi_k[k], Y ); dSdU=dSdU+ak*tmp;
|
||||||
|
MdagM.MpcDeriv(tmp , Y, MfMpvPhi_k[k] ); dSdU=dSdU+ak*tmp;
|
||||||
|
}
|
||||||
|
|
||||||
NumOp.ImportGauge(U);
|
//(2)
|
||||||
DenOp.ImportGauge(U);
|
//(3)
|
||||||
|
for(int k=0;k<n_pv;k++){
|
||||||
|
|
||||||
SchurDifferentiableOperator<Impl> VdagV(NumOp);
|
ak = PowerQuarter.residues[k];
|
||||||
SchurDifferentiableOperator<Impl> MdagM(DenOp);
|
|
||||||
|
|
||||||
ConjugateGradientMultiShift<FermionField> msCG_V(param.MaxIter,PowerQuarter);
|
VdagV.Mpc(MpvPhi_k[k],Y);
|
||||||
ConjugateGradientMultiShift<FermionField> msCG_M(param.MaxIter,PowerNegHalf);
|
VdagV.MpcDagDeriv(tmp,MpvMfMpvPhi_k[k],Y); dSdU=dSdU+ak*tmp;
|
||||||
|
VdagV.MpcDeriv (tmp,Y,MpvMfMpvPhi_k[k]); dSdU=dSdU+ak*tmp;
|
||||||
|
|
||||||
msCG_V(VdagV,PhiOdd,MpvPhi_k,MpvPhi);
|
VdagV.Mpc(MpvMfMpvPhi_k[k],Y); // V as we take Ydag
|
||||||
msCG_M(MdagM,MpvPhi,MfMpvPhi_k,MfMpvPhi);
|
VdagV.MpcDeriv (tmp,Y, MpvPhi_k[k]); dSdU=dSdU+ak*tmp;
|
||||||
msCG_V(VdagV,MfMpvPhi,MpvMfMpvPhi_k,MpvMfMpvPhi);
|
VdagV.MpcDagDeriv(tmp,MpvPhi_k[k], Y); dSdU=dSdU+ak*tmp;
|
||||||
|
|
||||||
RealD ak;
|
}
|
||||||
|
|
||||||
dSdU = zero;
|
//dSdU = Ta(dSdU);
|
||||||
|
|
||||||
// With these building blocks
|
};
|
||||||
//
|
};
|
||||||
// dS/dU =
|
|
||||||
// \sum_k -ak MfMpvPhi_k^dag [ dM^dag M + M^dag dM ] MfMpvPhi_k (1)
|
|
||||||
// + \sum_k -ak MpvMfMpvPhi_k^\dag [ dV^dag V + V^dag dV ] MpvPhi_k (2)
|
|
||||||
// -ak MpvPhi_k^dag [ dV^dag V + V^dag dV ] MpvMfMpvPhi_k (3)
|
|
||||||
|
|
||||||
//(1)
|
|
||||||
for(int k=0;k<n_f;k++){
|
|
||||||
ak = PowerNegHalf.residues[k];
|
|
||||||
MdagM.Mpc(MfMpvPhi_k[k],Y);
|
|
||||||
MdagM.MpcDagDeriv(tmp , MfMpvPhi_k[k], Y ); dSdU=dSdU+ak*tmp;
|
|
||||||
MdagM.MpcDeriv(tmp , Y, MfMpvPhi_k[k] ); dSdU=dSdU+ak*tmp;
|
|
||||||
}
|
|
||||||
|
|
||||||
//(2)
|
|
||||||
//(3)
|
|
||||||
for(int k=0;k<n_pv;k++){
|
|
||||||
|
|
||||||
ak = PowerQuarter.residues[k];
|
|
||||||
|
|
||||||
VdagV.Mpc(MpvPhi_k[k],Y);
|
|
||||||
VdagV.MpcDagDeriv(tmp,MpvMfMpvPhi_k[k],Y); dSdU=dSdU+ak*tmp;
|
|
||||||
VdagV.MpcDeriv (tmp,Y,MpvMfMpvPhi_k[k]); dSdU=dSdU+ak*tmp;
|
|
||||||
|
|
||||||
VdagV.Mpc(MpvMfMpvPhi_k[k],Y); // V as we take Ydag
|
|
||||||
VdagV.MpcDeriv (tmp,Y, MpvPhi_k[k]); dSdU=dSdU+ak*tmp;
|
|
||||||
VdagV.MpcDagDeriv(tmp,MpvPhi_k[k], Y); dSdU=dSdU+ak*tmp;
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
//dSdU = Ta(dSdU);
|
|
||||||
|
|
||||||
};
|
|
||||||
};
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
|
NAMESPACE_END(Grid);
|
||||||
|
|
||||||
#endif
|
#endif
|
||||||
|
|||||||
Loading…
x
Reference in New Issue
Block a user