portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-23 18:22:02 +01:00
Code Releases Activity

Merge develop

Browse Source
This commit is contained in:
Peter Boyle
2019-07-16 11:59:56 +01:00
parent fa9cd50c5b
commit 08904f830e
164 changed files with 5901 additions and 228 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
Grid/qcd/action/pseudofermion/Bounds.h
View File

@ -1,7 +1,6 @@
#pragma once
namespace Grid{
namespace QCD{
NAMESPACE_BEGIN(Grid);
template<class Field>
void HighBoundCheck(LinearOperatorBase<Field> &HermOp,
@ -20,7 +19,7 @@ namespace Grid{
Field &GaussNoise,
MultiShiftFunction &PowerNegHalf)
{
GridBase *FermionGrid = GaussNoise._grid;
GridBase *FermionGrid = GaussNoise.Grid();
Field X(FermionGrid);
Field Y(FermionGrid);
@ -49,5 +48,5 @@ namespace Grid{
assert( (std::sqrt(Nd/Nx)<tol) && " InverseSqrtBoundsCheck ");
}
}
}
NAMESPACE_END(Grid);

33
Grid/qcd/action/pseudofermion/ExactOneFlavourRatio.h
View File

@ -38,8 +38,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
#ifndef QCD_PSEUDOFERMION_EXACT_ONE_FLAVOUR_RATIO_H
#define QCD_PSEUDOFERMION_EXACT_ONE_FLAVOUR_RATIO_H
namespace Grid{
namespace QCD{
NAMESPACE_BEGIN(Grid);
///////////////////////////////////////////////////////////////
// Exact one flavour implementation of DWF determinant ratio //
@ -66,6 +65,13 @@ namespace QCD{
FermionField Phi; // the pseudofermion field for this trajectory
public:
ExactOneFlavourRatioPseudoFermionAction(AbstractEOFAFermion<Impl>& _Lop,
AbstractEOFAFermion<Impl>& _Rop,
OperatorFunction<FermionField>& CG,
Params& p,
bool use_fc=false)
: ExactOneFlavourRatioPseudoFermionAction(_Lop,_Rop,CG,CG,CG,CG,CG,p,use_fc) {};
ExactOneFlavourRatioPseudoFermionAction(AbstractEOFAFermion<Impl>& _Lop,
AbstractEOFAFermion<Impl>& _Rop,
OperatorFunction<FermionField>& HeatbathCG,
@ -150,7 +156,7 @@ namespace QCD{
spProj(eta, tmp[0], -1, Lop.Ls);
Lop.Omega(tmp[0], tmp[1], -1, 0);
G5R5(CG_src, tmp[1]);
tmp[1] = zero;
tmp[1] = Zero();
for(int k=0; k<param.degree; ++k){
gamma_l = 1.0 / ( 1.0 + PowerNegHalf.poles[k] );
Lop.RefreshShiftCoefficients(-gamma_l);
@ -160,7 +166,7 @@ namespace QCD{
SolverHB(Lop, CG_src, CG_soln);
prev_solns.push_back(CG_soln);
} else {
CG_soln = zero; // Just use zero as the initial guess
CG_soln = Zero(); // Just use zero as the initial guess
SolverHB(Lop, CG_src, CG_soln);
}
Lop.Dtilde(CG_soln, tmp[0]); // We actually solved Cayley preconditioned system: transform back
@ -176,7 +182,7 @@ namespace QCD{
spProj(eta, tmp[0], 1, Rop.Ls);
Rop.Omega(tmp[0], tmp[1], 1, 0);
G5R5(CG_src, tmp[1]);
tmp[1] = zero;
tmp[1] = Zero();
if(use_heatbath_forecasting){ prev_solns.clear(); } // empirically, LH solns don't help for RH solves
for(int k=0; k<param.degree; ++k){
gamma_l = 1.0 / ( 1.0 + PowerNegHalf.poles[k] );
@ -187,7 +193,7 @@ namespace QCD{
SolverHB(Rop, CG_src, CG_soln);
prev_solns.push_back(CG_soln);
} else {
CG_soln = zero;
CG_soln = Zero();
SolverHB(Rop, CG_src, CG_soln);
}
Rop.Dtilde(CG_soln, tmp[0]); // We actually solved Cayley preconditioned system: transform back
@ -222,7 +228,7 @@ namespace QCD{
spProj(Phi, spProj_Phi, -1, Lop.Ls);
Lop.Omega(spProj_Phi, tmp[0], -1, 0);
G5R5(tmp[1], tmp[0]);
tmp[0] = zero;
tmp[0] = Zero();
SolverL(Lop, tmp[1], tmp[0]);
Lop.Dtilde(tmp[0], tmp[1]); // We actually solved Cayley preconditioned system: transform back
Lop.Omega(tmp[1], tmp[0], -1, 1);
@ -233,7 +239,7 @@ namespace QCD{
spProj(Phi, spProj_Phi, 1, Rop.Ls);
Rop.Omega(spProj_Phi, tmp[0], 1, 0);
G5R5(tmp[1], tmp[0]);
tmp[0] = zero;
tmp[0] = Zero();
SolverR(Rop, tmp[1], tmp[0]);
Rop.Dtilde(tmp[0], tmp[1]);
Rop.Omega(tmp[1], tmp[0], 1, 1);
@ -257,7 +263,7 @@ namespace QCD{
spProj(Phi, spProj_Phi, -1, Lop.Ls);
Lop.Omega(spProj_Phi, tmp[0], -1, 0);
G5R5(tmp[1], tmp[0]);
tmp[0] = zero;
tmp[0] = Zero();
SolverL(Lop, tmp[1], tmp[0]);
Lop.Dtilde(tmp[0], tmp[1]); // We actually solved Cayley preconditioned system: transform back
Lop.Omega(tmp[1], tmp[0], -1, 1);
@ -268,7 +274,7 @@ namespace QCD{
spProj(Phi, spProj_Phi, 1, Rop.Ls);
Rop.Omega(spProj_Phi, tmp[0], 1, 0);
G5R5(tmp[1], tmp[0]);
tmp[0] = zero;
tmp[0] = Zero();
SolverR(Rop, tmp[1], tmp[0]);
Rop.Dtilde(tmp[0], tmp[1]);
Rop.Omega(tmp[1], tmp[0], 1, 1);
@ -301,7 +307,7 @@ namespace QCD{
spProj(Phi, spProj_Phi, -1, Lop.Ls);
Lop.Omega(spProj_Phi, Omega_spProj_Phi, -1, 0);
G5R5(CG_src, Omega_spProj_Phi);
spProj_Phi = zero;
spProj_Phi = Zero();
DerivativeSolverL(Lop, CG_src, spProj_Phi);
Lop.Dtilde(spProj_Phi, Chi);
G5R5(g5_R5_Chi, Chi);
@ -313,7 +319,7 @@ namespace QCD{
spProj(Phi, spProj_Phi, 1, Rop.Ls);
Rop.Omega(spProj_Phi, Omega_spProj_Phi, 1, 0);
G5R5(CG_src, Omega_spProj_Phi);
spProj_Phi = zero;
spProj_Phi = Zero();
DerivativeSolverR(Rop, CG_src, spProj_Phi);
Rop.Dtilde(spProj_Phi, Chi);
G5R5(g5_R5_Chi, Chi);
@ -321,6 +327,7 @@ namespace QCD{
dSdU = dSdU + Rop.k * force;
};
};
}}
NAMESPACE_END(Grid);
#endif

10
Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h
View File

@ -28,8 +28,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
#ifndef QCD_PSEUDOFERMION_ONE_FLAVOUR_EVEN_ODD_RATIONAL_RATIO_H
#define QCD_PSEUDOFERMION_ONE_FLAVOUR_EVEN_ODD_RATIONAL_RATIO_H
namespace Grid{
namespace QCD{
NAMESPACE_BEGIN(Grid);
///////////////////////////////////////
// One flavour rational
@ -145,7 +144,7 @@ namespace Grid{
assert(NumOp.ConstEE() == 1);
assert(DenOp.ConstEE() == 1);
PhiEven = zero;
PhiEven = Zero();
};
@ -245,7 +244,7 @@ namespace Grid{
RealD ak;
dSdU = zero;
dSdU = Zero();
// With these building blocks
//
@ -282,8 +281,7 @@ namespace Grid{
};
};
}
}
NAMESPACE_END(Grid);
#endif

8
Grid/qcd/action/pseudofermion/OneFlavourRational.h
View File

@ -28,8 +28,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
#ifndef QCD_PSEUDOFERMION_ONE_FLAVOUR_RATIONAL_H
#define QCD_PSEUDOFERMION_ONE_FLAVOUR_RATIONAL_H
namespace Grid{
namespace QCD{
NAMESPACE_BEGIN(Grid);
///////////////////////////////////////
// One flavour rational
@ -196,7 +195,7 @@ namespace Grid{
msCG(MdagMOp,Phi,MPhi_k);
dSdU = zero;
dSdU = Zero();
for(int k=0;k<Npole;k++){
RealD ak = PowerNegHalf.residues[k];
@ -214,8 +213,7 @@ namespace Grid{
};
};
}
}
NAMESPACE_END(Grid);
#endif

8
Grid/qcd/action/pseudofermion/OneFlavourRationalRatio.h
View File

@ -28,8 +28,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
#ifndef QCD_PSEUDOFERMION_ONE_FLAVOUR_RATIONAL_RATIO_H
#define QCD_PSEUDOFERMION_ONE_FLAVOUR_RATIONAL_RATIO_H
namespace Grid{
namespace QCD{
NAMESPACE_BEGIN(Grid);
///////////////////////////////////////
// One flavour rational
@ -231,7 +230,7 @@ namespace Grid{
RealD ak;
dSdU = zero;
dSdU = Zero();
// With these building blocks
//
@ -268,8 +267,7 @@ namespace Grid{
};
};
}
}
NAMESPACE_END(Grid);
#endif

14
Grid/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h
View File

@ -29,8 +29,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
#ifndef QCD_PSEUDOFERMION_TWO_FLAVOUR_EVEN_ODD_RATIO_H
#define QCD_PSEUDOFERMION_TWO_FLAVOUR_EVEN_ODD_RATIO_H
namespace Grid{
namespace QCD{
NAMESPACE_BEGIN(Grid);
///////////////////////////////////////
// Two flavour ratio
@ -118,7 +117,7 @@ namespace Grid{
// Odd det factors
Mpc.MpcDag(etaOdd,PhiOdd);
tmp=zero;
tmp=Zero();
HeatbathSolver(Vpc,PhiOdd,tmp);
Vpc.Mpc(tmp,PhiOdd);
@ -146,7 +145,7 @@ namespace Grid{
FermionField Y(NumOp.FermionRedBlackGrid());
Vpc.MpcDag(PhiOdd,Y); // Y= Vdag phi
X=zero;
X=Zero();
ActionSolver(Mpc,Y,X); // X= (MdagM)^-1 Vdag phi
//Mpc.Mpc(X,Y); // Y= Mdag^-1 Vdag phi
// Multiply by Ydag
@ -182,13 +181,13 @@ namespace Grid{
FermionField Y(NumOp.FermionRedBlackGrid());
// This assignment is necessary to be compliant with the HMC grids
GaugeField force(dSdU._grid);
GaugeField force(dSdU.Grid());
//Y=Vdag phi
//X = (Mdag M)^-1 V^dag phi
//Y = (Mdag)^-1 V^dag phi
Vpc.MpcDag(PhiOdd,Y); // Y= Vdag phi
X=zero;
X=Zero();
DerivativeSolver(Mpc,Y,X); // X= (MdagM)^-1 Vdag phi
Mpc.Mpc(X,Y); // Y= Mdag^-1 Vdag phi
@ -212,6 +211,5 @@ namespace Grid{
};
};
}
}
NAMESPACE_END(Grid);
#endif