mirror of
https://github.com/paboyle/Grid.git
synced 2025-04-25 05:05:56 +01:00
RdagR inverse
This commit is contained in:
Peter Boyle
parent
34ca4dd53a
commit
09b233b82e
@ -38,7 +38,6 @@ public:
|
||||
INHERIT_IMPL_TYPES(Impl);
|
||||
|
||||
private:
|
||||
SchurFactoredFermionOperator<Impl> & NumOp;// the basic operator
|
||||
SchurFactoredFermionOperator<Impl> & DenOp;// the basic operator
|
||||
|
||||
OperatorFunction<FermionField> &DerivativeSolver;
|
||||
@ -46,16 +45,16 @@ private:
|
||||
|
||||
FermionField Phi; // the pseudo fermion field for this trajectory
|
||||
|
||||
RealD refresh_action;
|
||||
public:
|
||||
DomainBoundaryPseudoFermionAction(SchurFactoredFermionOperator<Impl> &_NumOp,
|
||||
SchurFactoredFermionOperator<Impl> &_DenOp,
|
||||
DomainBoundaryPseudoFermionAction(SchurFactoredFermionOperator<Impl> &_DenOp,
|
||||
OperatorFunction<FermionField> & DS,
|
||||
OperatorFunction<FermionField> & AS
|
||||
) : NumOp(_NumOp), DenOp(_DenOp),
|
||||
) : DenOp(_DenOp),
|
||||
DerivativeSolver(DS), ActionSolver(AS),
|
||||
Phi(_NumOp.FermionGrid()) {};
|
||||
Phi(_DenOp.FermOp.FermionGrid()) {};
|
||||
|
||||
virtual std::string action_name(){return "DomainBoundaryPseudoFermionRatioAction";}
|
||||
virtual std::string action_name(){return "DomainBoundaryPseudoFermionAction";}
|
||||
|
||||
|
||||
virtual std::string LogParameters(){
|
||||
@ -63,14 +62,13 @@ public:
|
||||
return sstream.str();
|
||||
}
|
||||
|
||||
virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG)
|
||||
virtual void refresh(const GaugeField &U, GridSerialRNG& sRNG, GridParallelRNG& pRNG)
|
||||
{
|
||||
// P(phi) = e^{- phi^dag P^dag Rdag^-1 R^-1 P phi}
|
||||
// P(phi) = e^{- phi^dag Rdag^-1 R^-1 phi}
|
||||
//
|
||||
// NumOp == P
|
||||
// DenOp == R
|
||||
//
|
||||
// Take phi = P^{-1} R eta ; eta = R^-1 P Phi
|
||||
// Take phi = R eta ; eta = R^-1 Phi
|
||||
//
|
||||
// P(eta) = e^{- eta^dag eta}
|
||||
//
|
||||
@ -80,55 +78,40 @@ public:
|
||||
//
|
||||
RealD scale = std::sqrt(0.5);
|
||||
|
||||
NumOp.ImportGauge(U);
|
||||
DenOp.ImportGauge(U);
|
||||
|
||||
FermionField eta(NumOp.FermOp.FermionGrid());
|
||||
FermionField tmp(NumOp.FermOp.FermionGrid());
|
||||
FermionField eta(DenOp.FermOp.FermionGrid());
|
||||
|
||||
gaussian(pRNG,eta);
|
||||
gaussian(pRNG,eta); eta=eta*scale;
|
||||
|
||||
DenOp.ProjectBoundaryBar(eta);
|
||||
DenOp.R(eta,Phi);
|
||||
|
||||
NumOp.ProjectBoundaryBar(eta);
|
||||
|
||||
NumOp.R(eta,tmp);
|
||||
DenOp.RInv(tmp,Phi);
|
||||
|
||||
Phi=Phi*scale;
|
||||
|
||||
NumOp.ProjectBoundaryBar(Phi);
|
||||
refresh_action = norm2(eta);
|
||||
};
|
||||
|
||||
//////////////////////////////////////////////////////
|
||||
// S = phi^dag V (Mdag M)^-1 Vdag phi
|
||||
// S = phi^dag Rdag^-1 R^-1 phi
|
||||
//////////////////////////////////////////////////////
|
||||
virtual RealD S(const GaugeField &U) {
|
||||
|
||||
NumOp.ImportGauge(U);
|
||||
DenOp.ImportGauge(U);
|
||||
|
||||
FermionField X(NumOp.FermOp.FermionGrid());
|
||||
FermionField Y(NumOp.FermOp.FermionGrid());
|
||||
FermionField X(DenOp.FermOp.FermionGrid());
|
||||
|
||||
NumOp.R(Phi,Y);
|
||||
DenOp.RInv(Y,X);
|
||||
DenOp.RInv(Phi,X);
|
||||
|
||||
RealD action = norm2(X);
|
||||
|
||||
return action;
|
||||
};
|
||||
|
||||
//////////////////////////////////////////////////////
|
||||
// dS/du = phi^dag dV (Mdag M)^-1 V^dag phi
|
||||
// - phi^dag V (Mdag M)^-1 [ Mdag dM + dMdag M ] (Mdag M)^-1 V^dag phi
|
||||
// + phi^dag V (Mdag M)^-1 dV^dag phi
|
||||
//////////////////////////////////////////////////////
|
||||
virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
|
||||
|
||||
NumOp.ImportGauge(U);
|
||||
virtual void deriv(const GaugeField &U,GaugeField & dSdU)
|
||||
{
|
||||
DenOp.ImportGauge(U);
|
||||
|
||||
GridBase *fgrid = NumOp.FermOp.FermionGrid();
|
||||
GridBase *ugrid = NumOp.FermOp.GaugeGrid();
|
||||
GridBase *fgrid = DenOp.FermOp.FermionGrid();
|
||||
GridBase *ugrid = DenOp.FermOp.GaugeGrid();
|
||||
|
||||
FermionField X(fgrid);
|
||||
FermionField Y(fgrid);
|
||||
@ -136,66 +119,34 @@ public:
|
||||
|
||||
GaugeField force(ugrid);
|
||||
|
||||
FermionField DobiDdbPhi(fgrid); // Vector A in my notes
|
||||
FermionField DoiDdDobiDdbPhi(fgrid); // Vector B in my notes
|
||||
FermionField DiDdbP_Phi(fgrid); // Vector C in my notes
|
||||
FermionField DidRinvP_Phi(fgrid); // Vector D in my notes
|
||||
FermionField DoidRinvDagRinvP_Phi(fgrid); // Vector E in my notes
|
||||
FermionField DobidDddDoidRinvDagRinvP_Phi(fgrid); // Vector F in my notes
|
||||
FermionField DiDdb_Phi(fgrid); // Vector C in my notes
|
||||
FermionField DidRinv_Phi(fgrid); // Vector D in my notes
|
||||
FermionField DdbdDidRinv_Phi(fgrid);
|
||||
|
||||
FermionField P_Phi(fgrid);
|
||||
FermionField RinvP_Phi(fgrid);
|
||||
FermionField RinvDagRinvP_Phi(fgrid);
|
||||
FermionField Rinv_Phi(fgrid);
|
||||
FermionField RinvDagRinv_Phi(fgrid);
|
||||
|
||||
// P term
|
||||
NumOp.dBoundaryBar(Phi,tmp);
|
||||
NumOp.dOmegaBarInv(tmp,DobiDdbPhi); // Vector A
|
||||
NumOp.dBoundary(DobiDdbPhi,tmp);
|
||||
NumOp.dOmegaInv(tmp,DoiDdDobiDdbPhi); // Vector B
|
||||
P_Phi = Phi - DoiDdDobiDdbPhi;
|
||||
NumOp.ProjectBoundaryBar(P_Phi);
|
||||
// R^-1 term
|
||||
DenOp.dBoundaryBar(Phi,tmp);
|
||||
DenOp.Dinverse(tmp,DiDdb_Phi); // Vector C
|
||||
Rinv_Phi = Phi - DiDdb_Phi;
|
||||
DenOp.ProjectBoundaryBar(Rinv_Phi);
|
||||
|
||||
// R^-dagger R^-1 term
|
||||
DenOp.DinverseDag(Rinv_Phi,DidRinv_Phi); // Vector D
|
||||
DenOp.dBoundaryBarDag(DidRinv_Phi,DdbdDidRinv_Phi);
|
||||
RinvDagRinv_Phi = Rinv_Phi - DdbdDidRinv_Phi;
|
||||
DenOp.ProjectBoundaryBar(RinvDagRinv_Phi);
|
||||
|
||||
// R^-1 P term
|
||||
DenOp.dBoundaryBar(P_Phi,tmp);
|
||||
DenOp.Dinverse(tmp,DiDdbP_Phi); // Vector C
|
||||
RinvP_Phi = P_Phi - DiDdbP_Phi;
|
||||
DenOp.ProjectBoundaryBar(RinvP_Phi);
|
||||
|
||||
// R^-dagger R^-1 P term
|
||||
DenOp.DinverseDag(RinvP_Phi,DidRinvP_Phi); // Vector D
|
||||
RinvDagRinvP_Phi = RinvP_Phi - DidRinvP_Phi;
|
||||
DenOp.ProjectBoundaryBar(RinvDagRinvP_Phi);
|
||||
|
||||
// P^dag R^-dagger R^-1 P term
|
||||
NumOp.dOmegaDagInv(RinvDagRinvP_Phi,DoidRinvDagRinvP_Phi); // Vector E
|
||||
NumOp.dBoundaryDag(DoidRinvDagRinvP_Phi,tmp);
|
||||
NumOp.dOmegaBarDagInv(tmp,DobidDddDoidRinvDagRinvP_Phi); // Vector F
|
||||
|
||||
|
||||
dSdU=Zero();
|
||||
|
||||
// phi^dag V (Mdag M)^-1 dV^dag phi
|
||||
X = DobiDdbPhi;
|
||||
Y = DobidDddDoidRinvDagRinvP_Phi;
|
||||
NumOp.DirichletFermOp.MDeriv(force,X,Y,DaggerYes); dSdU=dsdU+force;
|
||||
NumOp.DirichletFermOp.MDeriv(force,Y,X,DaggerNo); dSdU=dSdU+force;
|
||||
|
||||
X = DoiDdDobiDdbPhi;
|
||||
Y = DoidRinvDagRinvP_Phi;
|
||||
NumOp.DirichletFermOp.MDeriv(force,X,Y,DaggerYes); dSdU=dsdU+force;
|
||||
NumOp.DirichletFermOp.MDeriv(force,Y,X,DaggerNo); dSdU=dSdU+force;
|
||||
|
||||
X = DiDdbP_Phi;
|
||||
Y = DidRinvP_Phi;
|
||||
NumOp.DirichletFermOp.MDeriv(force,X,Y,DaggerYes); dSdU=dsdU+force;
|
||||
NumOp.DirichletFermOp.MDeriv(force,Y,X,DaggerNo); dSdU=dSdU+force;
|
||||
X = DiDdb_Phi;
|
||||
Y = DidRinv_Phi;
|
||||
DenOp.FermOp.MDeriv(force,Y,X,DaggerNo); dSdU=force;
|
||||
DenOp.FermOp.MDeriv(force,X,Y,DaggerYes); dSdU=dSdU+force;
|
||||
|
||||
dSdU *= -1.0;
|
||||
//dSdU = - Ta(dSdU);
|
||||
|
||||
};
|
||||
};
|
||||
|
||||
NAMESPACE_END(Grid);
|
||||
|
||||
#endif
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user