portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-25 13:15:55 +01:00
Code Releases Activity

Remove the norm in MdagM

Browse Source
This commit is contained in:
Peter Boyle 2020-05-12 17:55:53 -04:00
parent ea08f193e7
commit 82f71643a4
27 changed files with 469 additions and 476 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

31
Grid/algorithms/CoarsenedMatrix.h
View File

@ -541,17 +541,14 @@ public:
/////////////////////// ///////////////////////
GridBase * Grid(void) { return _grid; }; // this is all the linalg routines need to know GridBase * Grid(void) { return _grid; }; // this is all the linalg routines need to know
RealD M (const CoarseVector &in, CoarseVector &out){ void M (const CoarseVector &in, CoarseVector &out)
{
conformable(_grid,in.Grid()); conformable(_grid,in.Grid());
conformable(in.Grid(),out.Grid()); conformable(in.Grid(),out.Grid());
// RealD Nin = norm2(in);
SimpleCompressor<siteVector> compressor; SimpleCompressor<siteVector> compressor;
double comms_usec = -usecond();
Stencil.HaloExchange(in,compressor); Stencil.HaloExchange(in,compressor);
comms_usec += usecond();
auto in_v = in.View(); auto in_v = in.View();
auto out_v = out.View(); auto out_v = out.View();
@ -565,12 +562,7 @@ public:
typedef decltype(coalescedRead(in_v[0])) calcVector; typedef decltype(coalescedRead(in_v[0])) calcVector;
typedef decltype(coalescedRead(in_v[0](0))) calcComplex; typedef decltype(coalescedRead(in_v[0](0))) calcComplex;
GridStopWatch ArithmeticTimer;
int osites=Grid()->oSites(); int osites=Grid()->oSites();
// double flops = osites*Nsimd*nbasis*nbasis*8.0*geom.npoint;
// double bytes = osites*nbasis*nbasis*geom.npoint*sizeof(CComplex);
double usecs =-usecond();
// assert(geom.npoint==9);
accelerator_for(sss, Grid()->oSites()*nbasis, Nsimd, { accelerator_for(sss, Grid()->oSites()*nbasis, Nsimd, {
int ss = sss/nbasis; int ss = sss/nbasis;
@ -598,23 +590,9 @@ public:
} }
coalescedWrite(out_v[ss](b),res,lane); coalescedWrite(out_v[ss](b),res,lane);
}); });
usecs +=usecond();
double nrm_usec=-usecond();
RealD Nout= norm2(out);
nrm_usec+=usecond();
/*
std::cout << GridLogMessage << "\tNorm " << nrm_usec << " us" <<std::endl;
std::cout << GridLogMessage << "\tHalo " << comms_usec << " us" <<std::endl;
std::cout << GridLogMessage << "\tMatrix " << usecs << " us" <<std::endl;
std::cout << GridLogMessage << "\t mflop/s " << flops/usecs<<std::endl;
std::cout << GridLogMessage << "\t MB/s " << bytes/usecs<<std::endl;
*/
return Nout;
}; };
RealD Mdag (const CoarseVector &in, CoarseVector &out) void Mdag (const CoarseVector &in, CoarseVector &out)
{ {
if(hermitian) { if(hermitian) {
// corresponds to Petrov-Galerkin coarsening // corresponds to Petrov-Galerkin coarsening
@ -625,7 +603,6 @@ public:
G5C(tmp, in); G5C(tmp, in);
M(tmp, out); M(tmp, out);
G5C(out, out); G5C(out, out);
return norm2(out);
} }
}; };
void MdirComms(const CoarseVector &in) void MdirComms(const CoarseVector &in)
@ -870,8 +847,6 @@ public:
auto A_self = A[self_stencil].View(); auto A_self = A[self_stencil].View();
accelerator_for(ss, Grid()->oSites(), Fobj::Nsimd(),{ coalescedWrite(A_p[ss](j,i),oZProj_v(ss)); }); accelerator_for(ss, Grid()->oSites(), Fobj::Nsimd(),{ coalescedWrite(A_p[ss](j,i),oZProj_v(ss)); });
// if( disp!= 0 ) { accelerator_for(ss, Grid()->oSites(), Fobj::Nsimd(),{ coalescedWrite(A_p[ss](j,i),oZProj_v(ss)); });}
// accelerator_for(ss, Grid()->oSites(), Fobj::Nsimd(),{ coalescedWrite(A_self[ss](j,i),A_self(ss)(j,i)+iZProj_v(ss)); });
} }
} }

158
Grid/algorithms/LinearOperator.h
View File

@ -43,7 +43,6 @@ NAMESPACE_BEGIN(Grid);
///////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////////
template<class Field> class LinearOperatorBase { template<class Field> class LinearOperatorBase {
public: public:
// Support for coarsening to a multigrid // Support for coarsening to a multigrid
virtual void OpDiag (const Field &in, Field &out) = 0; // Abstract base virtual void OpDiag (const Field &in, Field &out) = 0; // Abstract base
virtual void OpDir (const Field &in, Field &out,int dir,int disp) = 0; // Abstract base virtual void OpDir (const Field &in, Field &out,int dir,int disp) = 0; // Abstract base
@ -94,7 +93,10 @@ public:
_Mat.Mdag(in,out); _Mat.Mdag(in,out);
} }
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
_Mat.MdagM(in,out,n1,n2); _Mat.MdagM(in,out);
ComplexD dot = innerProduct(in,out);
n1=real(dot);
n2=norm2(out);
} }
void HermOp(const Field &in, Field &out){ void HermOp(const Field &in, Field &out){
_Mat.MdagM(in,out); _Mat.MdagM(in,out);
@ -131,17 +133,14 @@ public:
assert(0); assert(0);
} }
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
_Mat.MdagM(in,out,n1,n2); HermOp(in,out);
out = out + _shift*in; ComplexD dot = innerProduct(in,out);
ComplexD dot;
dot= innerProduct(in,out);
n1=real(dot); n1=real(dot);
n2=norm2(out); n2=norm2(out);
} }
void HermOp(const Field &in, Field &out){ void HermOp(const Field &in, Field &out){
RealD n1,n2; _Mat.MdagM(in,out);
HermOpAndNorm(in,out,n1,n2); out = out + _shift*in;
} }
}; };
@ -170,7 +169,7 @@ public:
_Mat.M(in,out); _Mat.M(in,out);
} }
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
_Mat.M(in,out); HermOp(in,out);
ComplexD dot= innerProduct(in,out); n1=real(dot); ComplexD dot= innerProduct(in,out); n1=real(dot);
n2=norm2(out); n2=norm2(out);
} }
@ -216,21 +215,24 @@ public:
template<class Field> template<class Field>
class SchurOperatorBase : public LinearOperatorBase<Field> { class SchurOperatorBase : public LinearOperatorBase<Field> {
public: public:
virtual RealD Mpc (const Field &in, Field &out) =0; virtual void Mpc (const Field &in, Field &out) =0;
virtual RealD MpcDag (const Field &in, Field &out) =0; virtual void MpcDag (const Field &in, Field &out) =0;
virtual void MpcDagMpc(const Field &in, Field &out,RealD &ni,RealD &no) { virtual void MpcDagMpc(const Field &in, Field &out) {
Field tmp(in.Grid()); Field tmp(in.Grid());
tmp.Checkerboard() = in.Checkerboard(); tmp.Checkerboard() = in.Checkerboard();
ni=Mpc(in,tmp); Mpc(in,tmp);
no=MpcDag(tmp,out); MpcDag(tmp,out);
} }
virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
MpcDagMpc(in,out,n1,n2); MpcDagMpc(in,out);
ComplexD dot= innerProduct(in,out);
n1=real(dot);
n2=norm2(out);
} }
virtual void HermOp(const Field &in, Field &out){ virtual void HermOp(const Field &in, Field &out){
RealD n1,n2; out.Checkerboard() = in.Checkerboard();
HermOpAndNorm(in,out,n1,n2); MpcDagMpc(in,out);
} }
void Op (const Field &in, Field &out){ void Op (const Field &in, Field &out){
Mpc(in,out); Mpc(in,out);
@ -254,26 +256,24 @@ public:
public: public:
Matrix &_Mat; Matrix &_Mat;
SchurDiagMooeeOperator (Matrix &Mat): _Mat(Mat){}; SchurDiagMooeeOperator (Matrix &Mat): _Mat(Mat){};
virtual RealD Mpc (const Field &in, Field &out) { virtual void Mpc (const Field &in, Field &out) {
Field tmp(in.Grid()); Field tmp(in.Grid());
tmp.Checkerboard() = !in.Checkerboard(); tmp.Checkerboard() = !in.Checkerboard();
_Mat.Meooe(in,tmp); _Mat.Meooe(in,tmp);
_Mat.MooeeInv(tmp,out); _Mat.MooeeInv(tmp,out);
_Mat.Meooe(out,tmp); _Mat.Meooe(out,tmp);
_Mat.Mooee(in,out); _Mat.Mooee(in,out);
return axpy_norm(out,-1.0,tmp,out); axpy(out,-1.0,tmp,out);
} }
virtual RealD MpcDag (const Field &in, Field &out){ virtual void MpcDag (const Field &in, Field &out){
Field tmp(in.Grid()); Field tmp(in.Grid());
_Mat.MeooeDag(in,tmp); _Mat.MeooeDag(in,tmp);
_Mat.MooeeInvDag(tmp,out); _Mat.MooeeInvDag(tmp,out);
_Mat.MeooeDag(out,tmp); _Mat.MeooeDag(out,tmp);
_Mat.MooeeDag(in,out); _Mat.MooeeDag(in,out);
return axpy_norm(out,-1.0,tmp,out); axpy(out,-1.0,tmp,out);
} }
}; };
template<class Matrix,class Field> template<class Matrix,class Field>
@ -283,25 +283,23 @@ public:
public: public:
SchurDiagOneOperator (Matrix &Mat): _Mat(Mat){}; SchurDiagOneOperator (Matrix &Mat): _Mat(Mat){};
virtual RealD Mpc (const Field &in, Field &out) { virtual void Mpc (const Field &in, Field &out) {
Field tmp(in.Grid()); Field tmp(in.Grid());
_Mat.Meooe(in,out); _Mat.Meooe(in,out);
_Mat.MooeeInv(out,tmp); _Mat.MooeeInv(out,tmp);
_Mat.Meooe(tmp,out); _Mat.Meooe(tmp,out);
_Mat.MooeeInv(out,tmp); _Mat.MooeeInv(out,tmp);
axpy(out,-1.0,tmp,in);
return axpy_norm(out,-1.0,tmp,in);
} }
virtual RealD MpcDag (const Field &in, Field &out){ virtual void MpcDag (const Field &in, Field &out){
Field tmp(in.Grid()); Field tmp(in.Grid());
_Mat.MooeeInvDag(in,out); _Mat.MooeeInvDag(in,out);
_Mat.MeooeDag(out,tmp); _Mat.MeooeDag(out,tmp);
_Mat.MooeeInvDag(tmp,out); _Mat.MooeeInvDag(tmp,out);
_Mat.MeooeDag(out,tmp); _Mat.MeooeDag(out,tmp);
axpy(out,-1.0,tmp,in);
return axpy_norm(out,-1.0,tmp,in);
} }
}; };
template<class Matrix,class Field> template<class Matrix,class Field>
@ -311,7 +309,7 @@ public:
public: public:
SchurDiagTwoOperator (Matrix &Mat): _Mat(Mat){}; SchurDiagTwoOperator (Matrix &Mat): _Mat(Mat){};
virtual RealD Mpc (const Field &in, Field &out) { virtual void Mpc (const Field &in, Field &out) {
Field tmp(in.Grid()); Field tmp(in.Grid());
_Mat.MooeeInv(in,out); _Mat.MooeeInv(in,out);
@ -319,9 +317,9 @@ public:
_Mat.MooeeInv(tmp,out); _Mat.MooeeInv(tmp,out);
_Mat.Meooe(out,tmp); _Mat.Meooe(out,tmp);
return axpy_norm(out,-1.0,tmp,in); axpy(out,-1.0,tmp,in);
} }
virtual RealD MpcDag (const Field &in, Field &out){ virtual void MpcDag (const Field &in, Field &out){
Field tmp(in.Grid()); Field tmp(in.Grid());
_Mat.MeooeDag(in,out); _Mat.MeooeDag(in,out);
@ -329,7 +327,7 @@ public:
_Mat.MeooeDag(tmp,out); _Mat.MeooeDag(tmp,out);
_Mat.MooeeInvDag(out,tmp); _Mat.MooeeInvDag(out,tmp);
return axpy_norm(out,-1.0,tmp,in); axpy(out,-1.0,tmp,in);
} }
}; };
@ -337,13 +335,13 @@ public:
class NonHermitianSchurOperatorBase : public LinearOperatorBase<Field> class NonHermitianSchurOperatorBase : public LinearOperatorBase<Field>
{ {
public: public:
virtual RealD Mpc (const Field& in, Field& out) = 0; virtual void Mpc (const Field& in, Field& out) = 0;
virtual RealD MpcDag (const Field& in, Field& out) = 0; virtual void MpcDag (const Field& in, Field& out) = 0;
virtual void MpcDagMpc(const Field& in, Field& out, RealD& ni, RealD& no) { virtual void MpcDagMpc(const Field& in, Field& out) {
Field tmp(in.Grid()); Field tmp(in.Grid());
tmp.Checkerboard() = in.Checkerboard(); tmp.Checkerboard() = in.Checkerboard();
ni = Mpc(in,tmp); Mpc(in,tmp);
no = MpcDag(tmp,out); MpcDag(tmp,out);
} }
virtual void HermOpAndNorm(const Field& in, Field& out, RealD& n1, RealD& n2) { virtual void HermOpAndNorm(const Field& in, Field& out, RealD& n1, RealD& n2) {
assert(0); assert(0);
@ -375,7 +373,7 @@ public:
public: public:
Matrix& _Mat; Matrix& _Mat;
NonHermitianSchurDiagMooeeOperator(Matrix& Mat): _Mat(Mat){}; NonHermitianSchurDiagMooeeOperator(Matrix& Mat): _Mat(Mat){};
virtual RealD Mpc(const Field& in, Field& out) { virtual void Mpc(const Field& in, Field& out) {
Field tmp(in.Grid()); Field tmp(in.Grid());
tmp.Checkerboard() = !in.Checkerboard(); tmp.Checkerboard() = !in.Checkerboard();
@ -385,9 +383,9 @@ public:
_Mat.Mooee(in, out); _Mat.Mooee(in, out);
return axpy_norm(out, -1.0, tmp, out); axpy(out, -1.0, tmp, out);
} }
virtual RealD MpcDag(const Field& in, Field& out) { virtual void MpcDag(const Field& in, Field& out) {
Field tmp(in.Grid()); Field tmp(in.Grid());
_Mat.MeooeDag(in, tmp); _Mat.MeooeDag(in, tmp);
@ -396,7 +394,7 @@ public:
_Mat.MooeeDag(in, out); _Mat.MooeeDag(in, out);
return axpy_norm(out, -1.0, tmp, out); axpy(out, -1.0, tmp, out);
} }
}; };
@ -408,7 +406,7 @@ public:
public: public:
NonHermitianSchurDiagOneOperator (Matrix& Mat): _Mat(Mat){}; NonHermitianSchurDiagOneOperator (Matrix& Mat): _Mat(Mat){};
virtual RealD Mpc(const Field& in, Field& out) { virtual void Mpc(const Field& in, Field& out) {
Field tmp(in.Grid()); Field tmp(in.Grid());
_Mat.Meooe(in, out); _Mat.Meooe(in, out);
@ -416,9 +414,9 @@ public:
_Mat.Meooe(tmp, out); _Mat.Meooe(tmp, out);
_Mat.MooeeInv(out, tmp); _Mat.MooeeInv(out, tmp);
return axpy_norm(out, -1.0, tmp, in); axpy(out, -1.0, tmp, in);
} }
virtual RealD MpcDag(const Field& in, Field& out) { virtual void MpcDag(const Field& in, Field& out) {
Field tmp(in.Grid()); Field tmp(in.Grid());
_Mat.MooeeInvDag(in, out); _Mat.MooeeInvDag(in, out);
@ -426,7 +424,7 @@ public:
_Mat.MooeeInvDag(tmp, out); _Mat.MooeeInvDag(tmp, out);
_Mat.MeooeDag(out, tmp); _Mat.MeooeDag(out, tmp);
return axpy_norm(out, -1.0, tmp, in); axpy(out, -1.0, tmp, in);
} }
}; };
@ -439,7 +437,7 @@ public:
public: public:
NonHermitianSchurDiagTwoOperator(Matrix& Mat): _Mat(Mat){}; NonHermitianSchurDiagTwoOperator(Matrix& Mat): _Mat(Mat){};
virtual RealD Mpc(const Field& in, Field& out) { virtual void Mpc(const Field& in, Field& out) {
Field tmp(in.Grid()); Field tmp(in.Grid());
_Mat.MooeeInv(in, out); _Mat.MooeeInv(in, out);
@ -447,9 +445,9 @@ public:
_Mat.MooeeInv(tmp, out); _Mat.MooeeInv(tmp, out);
_Mat.Meooe(out, tmp); _Mat.Meooe(out, tmp);
return axpy_norm(out, -1.0, tmp, in); axpy(out, -1.0, tmp, in);
} }
virtual RealD MpcDag(const Field& in, Field& out) { virtual void MpcDag(const Field& in, Field& out) {
Field tmp(in.Grid()); Field tmp(in.Grid());
_Mat.MeooeDag(in, out); _Mat.MeooeDag(in, out);
@ -457,7 +455,7 @@ public:
_Mat.MeooeDag(tmp, out); _Mat.MeooeDag(tmp, out);
_Mat.MooeeInvDag(out, tmp); _Mat.MooeeInvDag(out, tmp);
return axpy_norm(out, -1.0, tmp, in); axpy(out, -1.0, tmp, in);
} }
}; };
@ -476,71 +474,38 @@ public:
Matrix &_Mat; Matrix &_Mat;
Field tmp; Field tmp;
RealD mass; RealD mass;
double tMpc;
double tIP;
double tMeo;
double taxpby_norm;
uint64_t ncall;
public: public:
void Report(void)
{
std::cout << GridLogMessage << " HermOpAndNorm.Mpc "<< tMpc/ncall<<" usec "<<std::endl;
std::cout << GridLogMessage << " HermOpAndNorm.IP "<< tIP /ncall<<" usec "<<std::endl;
std::cout << GridLogMessage << " Mpc.MeoMoe "<< tMeo/ncall<<" usec "<<std::endl;
std::cout << GridLogMessage << " Mpc.axpby_norm "<< taxpby_norm/ncall<<" usec "<<std::endl;
}
SchurStaggeredOperator (Matrix &Mat): _Mat(Mat), tmp(_Mat.RedBlackGrid()) SchurStaggeredOperator (Matrix &Mat): _Mat(Mat), tmp(_Mat.RedBlackGrid())
{ {
assert( _Mat.isTrivialEE() ); assert( _Mat.isTrivialEE() );
mass = _Mat.Mass(); mass = _Mat.Mass();
tMpc=0;
tIP =0;
tMeo=0;
taxpby_norm=0;
ncall=0;
} }
virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
ncall++; Mpc(in,out);
tMpc-=usecond();
n2 = Mpc(in,out);
tMpc+=usecond();
tIP-=usecond();
ComplexD dot= innerProduct(in,out); ComplexD dot= innerProduct(in,out);
tIP+=usecond();
n1 = real(dot); n1 = real(dot);
n2 =0.0;
} }
virtual void HermOp(const Field &in, Field &out){ virtual void HermOp(const Field &in, Field &out){
ncall++; Mpc(in,out);
tMpc-=usecond(); // _Mat.Meooe(in,out);
_Mat.Meooe(in,out); // _Mat.Meooe(out,tmp);
_Mat.Meooe(out,tmp); // axpby(out,-1.0,mass*mass,tmp,in);
tMpc+=usecond();
taxpby_norm-=usecond();
axpby(out,-1.0,mass*mass,tmp,in);
taxpby_norm+=usecond();
} }
virtual RealD Mpc (const Field &in, Field &out) virtual void Mpc (const Field &in, Field &out)
{ {
Field tmp(in.Grid()); Field tmp(in.Grid());
Field tmp2(in.Grid()); Field tmp2(in.Grid());
// std::cout << GridLogIterative << " HermOp.Mpc "<<std::endl; // _Mat.Mooee(in,out);
_Mat.Mooee(in,out); // _Mat.Mooee(out,tmp);
_Mat.Mooee(out,tmp);
// std::cout << GridLogIterative << " HermOp.MooeeMooee "<<std::endl;
tMeo-=usecond();
_Mat.Meooe(in,out); _Mat.Meooe(in,out);
_Mat.Meooe(out,tmp); _Mat.Meooe(out,tmp);
tMeo+=usecond(); axpby(out,-1.0,mass*mass,tmp,in);
taxpby_norm-=usecond();
RealD nn=axpby_norm(out,-1.0,mass*mass,tmp,in);
taxpby_norm+=usecond();
return nn;
} }
virtual RealD MpcDag (const Field &in, Field &out){ virtual void MpcDag (const Field &in, Field &out){
return Mpc(in,out); Mpc(in,out);
} }
virtual void MpcDagMpc(const Field &in, Field &out,RealD &ni,RealD &no) { virtual void MpcDagMpc(const Field &in, Field &out,RealD &ni,RealD &no) {
assert(0);// Never need with staggered assert(0);// Never need with staggered
@ -548,7 +513,6 @@ public:
}; };
template<class Matrix,class Field> using SchurStagOperator = SchurStaggeredOperator<Matrix,Field>; template<class Matrix,class Field> using SchurStagOperator = SchurStaggeredOperator<Matrix,Field>;
///////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////
// Base classes for functions of operators // Base classes for functions of operators
///////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////

14
Grid/algorithms/SparseMatrix.h
View File

@ -38,16 +38,12 @@ template<class Field> class SparseMatrixBase {
public: public:
virtual GridBase *Grid(void) =0; virtual GridBase *Grid(void) =0;
// Full checkerboar operations // Full checkerboar operations
virtual RealD M (const Field &in, Field &out)=0; virtual void M (const Field &in, Field &out)=0;
virtual RealD Mdag (const Field &in, Field &out)=0; virtual void Mdag (const Field &in, Field &out)=0;
virtual void MdagM(const Field &in, Field &out,RealD &ni,RealD &no) {
Field tmp (in.Grid());
ni=M(in,tmp);
no=Mdag(tmp,out);
}
virtual void MdagM(const Field &in, Field &out) { virtual void MdagM(const Field &in, Field &out) {
RealD ni, no; Field tmp (in.Grid());
MdagM(in,out,ni,no); M(in,tmp);
Mdag(tmp,out);
} }
virtual void Mdiag (const Field &in, Field &out)=0; virtual void Mdiag (const Field &in, Field &out)=0;
virtual void Mdir (const Field &in, Field &out,int dir, int disp)=0; virtual void Mdir (const Field &in, Field &out,int dir, int disp)=0;

4
Grid/qcd/action/fermion/CayleyFermion5D.h
View File

@ -40,8 +40,8 @@ public:
public: public:
// override multiply // override multiply
virtual RealD M (const FermionField &in, FermionField &out); virtual void M (const FermionField &in, FermionField &out);
virtual RealD Mdag (const FermionField &in, FermionField &out); virtual void Mdag (const FermionField &in, FermionField &out);
// half checkerboard operations // half checkerboard operations
virtual void Meooe (const FermionField &in, FermionField &out); virtual void Meooe (const FermionField &in, FermionField &out);

4
Grid/qcd/action/fermion/ContinuedFractionFermion5D.h
View File

@ -41,8 +41,8 @@ public:
public: public:
// override multiply // override multiply
virtual RealD M (const FermionField &in, FermionField &out); virtual void M (const FermionField &in, FermionField &out);
virtual RealD Mdag (const FermionField &in, FermionField &out); virtual void Mdag (const FermionField &in, FermionField &out);
// half checkerboard operaions // half checkerboard operaions
virtual void Meooe (const FermionField &in, FermionField &out); virtual void Meooe (const FermionField &in, FermionField &out);

4
Grid/qcd/action/fermion/DomainWallEOFAFermion.h
View File

@ -53,8 +53,8 @@ public:
virtual void DtildeInv (const FermionField& in, FermionField& out); virtual void DtildeInv (const FermionField& in, FermionField& out);
// override multiply // override multiply
virtual RealD M (const FermionField& in, FermionField& out); virtual void M (const FermionField& in, FermionField& out);
virtual RealD Mdag (const FermionField& in, FermionField& out); virtual void Mdag (const FermionField& in, FermionField& out);
// half checkerboard operations // half checkerboard operations
virtual void Mooee (const FermionField& in, FermionField& out); virtual void Mooee (const FermionField& in, FermionField& out);

5
Grid/qcd/action/fermion/FermionOperator.h
View File

@ -58,8 +58,8 @@ public:
virtual GridBase *GaugeRedBlackGrid(void) =0; virtual GridBase *GaugeRedBlackGrid(void) =0;
// override multiply // override multiply
virtual RealD M (const FermionField &in, FermionField &out)=0; virtual void M (const FermionField &in, FermionField &out)=0;
virtual RealD Mdag (const FermionField &in, FermionField &out)=0; virtual void Mdag (const FermionField &in, FermionField &out)=0;
// half checkerboard operaions // half checkerboard operaions
virtual void Meooe (const FermionField &in, FermionField &out)=0; virtual void Meooe (const FermionField &in, FermionField &out)=0;
@ -86,7 +86,6 @@ public:
virtual void DhopDerivEO(GaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0; virtual void DhopDerivEO(GaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0;
virtual void DhopDerivOE(GaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0; virtual void DhopDerivOE(GaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0;
virtual void Mdiag (const FermionField &in, FermionField &out) { Mooee(in,out);}; // Same as Mooee applied to both CB's virtual void Mdiag (const FermionField &in, FermionField &out) { Mooee(in,out);}; // Same as Mooee applied to both CB's
virtual void Mdir (const FermionField &in, FermionField &out,int dir,int disp)=0; // case by case Wilson, Clover, Cayley, ContFrac, PartFrac virtual void Mdir (const FermionField &in, FermionField &out,int dir,int disp)=0; // case by case Wilson, Clover, Cayley, ContFrac, PartFrac
virtual void MdirAll(const FermionField &in, std::vector<FermionField> &out)=0; // case by case Wilson, Clover, Cayley, ContFrac, PartFrac virtual void MdirAll(const FermionField &in, std::vector<FermionField> &out)=0; // case by case Wilson, Clover, Cayley, ContFrac, PartFrac

4
Grid/qcd/action/fermion/ImprovedStaggeredFermion.h
View File

@ -71,8 +71,8 @@ public:
// override multiply; cut number routines if pass dagger argument // override multiply; cut number routines if pass dagger argument
// and also make interface more uniformly consistent // and also make interface more uniformly consistent
////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////
RealD M(const FermionField &in, FermionField &out); void M(const FermionField &in, FermionField &out);
RealD Mdag(const FermionField &in, FermionField &out); void Mdag(const FermionField &in, FermionField &out);
///////////////////////////////////////////////////////// /////////////////////////////////////////////////////////
// half checkerboard operations // half checkerboard operations

5
Grid/qcd/action/fermion/ImprovedStaggeredFermion5D.h
View File

@ -1,4 +1,3 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
@ -74,8 +73,8 @@ public:
GridBase *FermionRedBlackGrid(void) { return _FiveDimRedBlackGrid;} GridBase *FermionRedBlackGrid(void) { return _FiveDimRedBlackGrid;}
// full checkerboard operations; leave unimplemented as abstract for now // full checkerboard operations; leave unimplemented as abstract for now
RealD M (const FermionField &in, FermionField &out); void M (const FermionField &in, FermionField &out);
RealD Mdag (const FermionField &in, FermionField &out); void Mdag (const FermionField &in, FermionField &out);
// half checkerboard operations // half checkerboard operations
void Meooe (const FermionField &in, FermionField &out); void Meooe (const FermionField &in, FermionField &out);

4
Grid/qcd/action/fermion/MobiusEOFAFermion.h
View File

@ -56,8 +56,8 @@ public:
virtual void DtildeInv (const FermionField& in, FermionField& out); virtual void DtildeInv (const FermionField& in, FermionField& out);
// override multiply // override multiply
virtual RealD M (const FermionField& in, FermionField& out); virtual void M (const FermionField& in, FermionField& out);
virtual RealD Mdag (const FermionField& in, FermionField& out); virtual void Mdag (const FermionField& in, FermionField& out);
// half checkerboard operations // half checkerboard operations
virtual void Mooee (const FermionField& in, FermionField& out); virtual void Mooee (const FermionField& in, FermionField& out);

4
Grid/qcd/action/fermion/PartialFractionFermion5D.h
View File

@ -47,8 +47,8 @@ public:
void M_internal(const FermionField &in, FermionField &out,int dag); void M_internal(const FermionField &in, FermionField &out,int dag);
// override multiply // override multiply
virtual RealD M (const FermionField &in, FermionField &out); virtual void M (const FermionField &in, FermionField &out);
virtual RealD Mdag (const FermionField &in, FermionField &out); virtual void Mdag (const FermionField &in, FermionField &out);
// half checkerboard operaions // half checkerboard operaions
virtual void Meooe (const FermionField &in, FermionField &out); virtual void Meooe (const FermionField &in, FermionField &out);

5
Grid/qcd/action/fermion/WilsonCloverFermion.h
View File

@ -109,9 +109,8 @@ public:
ImportGauge(_Umu); ImportGauge(_Umu);
} }
virtual RealD M(const FermionField &in, FermionField &out); virtual void M(const FermionField &in, FermionField &out);
virtual RealD Mdag(const FermionField &in, FermionField &out); virtual void Mdag(const FermionField &in, FermionField &out);
virtual void Mooee(const FermionField &in, FermionField &out); virtual void Mooee(const FermionField &in, FermionField &out);
virtual void MooeeDag(const FermionField &in, FermionField &out); virtual void MooeeDag(const FermionField &in, FermionField &out);
virtual void MooeeInv(const FermionField &in, FermionField &out); virtual void MooeeInv(const FermionField &in, FermionField &out);

4
Grid/qcd/action/fermion/WilsonFermion.h
View File

@ -78,8 +78,8 @@ public:
// override multiply; cut number routines if pass dagger argument // override multiply; cut number routines if pass dagger argument
// and also make interface more uniformly consistent // and also make interface more uniformly consistent
////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////
virtual RealD M(const FermionField &in, FermionField &out); virtual void M(const FermionField &in, FermionField &out);
virtual RealD Mdag(const FermionField &in, FermionField &out); virtual void Mdag(const FermionField &in, FermionField &out);
///////////////////////////////////////////////////////// /////////////////////////////////////////////////////////
// half checkerboard operations // half checkerboard operations

5
Grid/qcd/action/fermion/WilsonFermion5D.h
View File

@ -1,4 +1,3 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
@ -99,8 +98,8 @@ public:
GridBase *FermionRedBlackGrid(void) { return _FiveDimRedBlackGrid;} GridBase *FermionRedBlackGrid(void) { return _FiveDimRedBlackGrid;}
// full checkerboard operations; leave unimplemented as abstract for now // full checkerboard operations; leave unimplemented as abstract for now
virtual RealD M (const FermionField &in, FermionField &out){assert(0); return 0.0;}; virtual void M (const FermionField &in, FermionField &out){assert(0);};
virtual RealD Mdag (const FermionField &in, FermionField &out){assert(0); return 0.0;}; virtual void Mdag (const FermionField &in, FermionField &out){assert(0);};
// half checkerboard operations; leave unimplemented as abstract for now // half checkerboard operations; leave unimplemented as abstract for now
virtual void Meooe (const FermionField &in, FermionField &out){assert(0);}; virtual void Meooe (const FermionField &in, FermionField &out){assert(0);};

8
Grid/qcd/action/fermion/WilsonTMFermion5D.h
View File

@ -120,7 +120,8 @@ class WilsonTMFermion5D : public WilsonFermion5D<Impl>
} }
} }
virtual RealD M(const FermionField &in, FermionField &out) { virtual void M(const FermionField &in, FermionField &out)
{
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
this->Dhop(in, out, DaggerNo); this->Dhop(in, out, DaggerNo);
FermionField tmp(out.Grid()); FermionField tmp(out.Grid());
@ -129,11 +130,12 @@ class WilsonTMFermion5D : public WilsonFermion5D<Impl>
ComplexD b(0.0,this->mu[s]); ComplexD b(0.0,this->mu[s]);
axpbg5y_ssp(tmp,a,in,b,in,s,s); axpbg5y_ssp(tmp,a,in,b,in,s,s);
} }
return axpy_norm(out, 1.0, tmp, out); axpy(out, 1.0, tmp, out);
} }
// needed for fast PV // needed for fast PV
void update(const std::vector<RealD>& _mass, const std::vector<RealD>& _mu) { void update(const std::vector<RealD>& _mass, const std::vector<RealD>& _mu)
{
assert(_mass.size() == _mu.size()); assert(_mass.size() == _mu.size());
assert(_mass.size() == this->FermionGrid()->_fdimensions[0]); assert(_mass.size() == this->FermionGrid()->_fdimensions[0]);
this->mass = _mass; this->mass = _mass;

6
Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h
View File

@ -323,7 +323,7 @@ void CayleyFermion5D<Impl>::MeooeDag5D (const FermionField &psi, FermionField
} }
template<class Impl> template<class Impl>
RealD CayleyFermion5D<Impl>::M (const FermionField &psi, FermionField &chi) void CayleyFermion5D<Impl>::M (const FermionField &psi, FermionField &chi)
{ {
FermionField Din(psi.Grid()); FermionField Din(psi.Grid());
@ -335,11 +335,10 @@ RealD CayleyFermion5D<Impl>::M (const FermionField &psi, FermionField &chi)
axpby(chi,1.0,1.0,chi,psi); axpby(chi,1.0,1.0,chi,psi);
M5D(psi,chi); M5D(psi,chi);
return(norm2(chi));
} }
template<class Impl> template<class Impl>
RealD CayleyFermion5D<Impl>::Mdag (const FermionField &psi, FermionField &chi) void CayleyFermion5D<Impl>::Mdag (const FermionField &psi, FermionField &chi)
{ {
// Under adjoint // Under adjoint
//D1+ D1- P- -> D1+^dag P+ D2-^dag //D1+ D1- P- -> D1+^dag P+ D2-^dag
@ -354,7 +353,6 @@ RealD CayleyFermion5D<Impl>::Mdag (const FermionField &psi, FermionField &chi)
M5Ddag(psi,chi); M5Ddag(psi,chi);
// ((b D_W + D_w hop terms +1) on s-diag // ((b D_W + D_w hop terms +1) on s-diag
axpby (chi,1.0,1.0,chi,psi); axpby (chi,1.0,1.0,chi,psi);
return norm2(chi);
} }
// half checkerboard operations // half checkerboard operations

7
Grid/qcd/action/fermion/implementation/ContinuedFractionFermion5DImplementation.h
View File

@ -94,7 +94,7 @@ void ContinuedFractionFermion5D<Impl>::SetCoefficientsZolotarev(RealD zolo_hi,Ap
template<class Impl> template<class Impl>
RealD ContinuedFractionFermion5D<Impl>::M (const FermionField &psi, FermionField &chi) void ContinuedFractionFermion5D<Impl>::M (const FermionField &psi, FermionField &chi)
{ {
int Ls = this->Ls; int Ls = this->Ls;
@ -116,15 +116,14 @@ RealD ContinuedFractionFermion5D<Impl>::M (const FermionField &psi, F
} }
sign=-sign; sign=-sign;
} }
return norm2(chi);
} }
template<class Impl> template<class Impl>
RealD ContinuedFractionFermion5D<Impl>::Mdag (const FermionField &psi, FermionField &chi) void ContinuedFractionFermion5D<Impl>::Mdag (const FermionField &psi, FermionField &chi)
{ {
// This matrix is already hermitian. (g5 Dw) = Dw dag g5 = (g5 Dw)dag // This matrix is already hermitian. (g5 Dw) = Dw dag g5 = (g5 Dw)dag
// The rest of matrix is symmetric. // The rest of matrix is symmetric.
// Can ignore "dag" // Can ignore "dag"
return M(psi,chi); M(psi,chi);
} }
template<class Impl> template<class Impl>
void ContinuedFractionFermion5D<Impl>::Mdir (const FermionField &psi, FermionField &chi,int dir,int disp){ void ContinuedFractionFermion5D<Impl>::Mdir (const FermionField &psi, FermionField &chi,int dir,int disp){

6
Grid/qcd/action/fermion/implementation/DomainWallEOFAFermionImplementation.h
View File

@ -89,7 +89,7 @@ void DomainWallEOFAFermion<Impl>::DtildeInv(const FermionField& psi, FermionFiel
/*****************************************************************************************************/ /*****************************************************************************************************/
template<class Impl> template<class Impl>
RealD DomainWallEOFAFermion<Impl>::M(const FermionField& psi, FermionField& chi) void DomainWallEOFAFermion<Impl>::M(const FermionField& psi, FermionField& chi)
{ {
FermionField Din(psi.Grid()); FermionField Din(psi.Grid());
@ -97,11 +97,10 @@ RealD DomainWallEOFAFermion<Impl>::M(const FermionField& psi, FermionField& chi)
this->DW(Din, chi, DaggerNo); this->DW(Din, chi, DaggerNo);
axpby(chi, 1.0, 1.0, chi, psi); axpby(chi, 1.0, 1.0, chi, psi);
this->M5D(psi, chi); this->M5D(psi, chi);
return(norm2(chi));
} }
template<class Impl> template<class Impl>
RealD DomainWallEOFAFermion<Impl>::Mdag(const FermionField& psi, FermionField& chi) void DomainWallEOFAFermion<Impl>::Mdag(const FermionField& psi, FermionField& chi)
{ {
FermionField Din(psi.Grid()); FermionField Din(psi.Grid());
@ -109,7 +108,6 @@ RealD DomainWallEOFAFermion<Impl>::Mdag(const FermionField& psi, FermionField& c
this->MeooeDag5D(Din, chi); this->MeooeDag5D(Din, chi);
this->M5Ddag(psi, chi); this->M5Ddag(psi, chi);
axpby(chi, 1.0, 1.0, chi, psi); axpby(chi, 1.0, 1.0, chi, psi);
return(norm2(chi));
} }
/******************************************************************** /********************************************************************

29
Grid/qcd/action/fermion/implementation/ImprovedStaggeredFermion5DImplementation.h
View File

@ -548,21 +548,24 @@ void ImprovedStaggeredFermion5D<Impl>::MdirAll(const FermionField &in, std::vect
assert(0); assert(0);
} }
template <class Impl> template <class Impl>
RealD ImprovedStaggeredFermion5D<Impl>::M(const FermionField &in, FermionField &out) { void ImprovedStaggeredFermion5D<Impl>::M(const FermionField &in, FermionField &out)
{
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
Dhop(in, out, DaggerNo); Dhop(in, out, DaggerNo);
return axpy_norm(out, mass, in, out); axpy(out, mass, in, out);
} }
template <class Impl> template <class Impl>
RealD ImprovedStaggeredFermion5D<Impl>::Mdag(const FermionField &in, FermionField &out) { void ImprovedStaggeredFermion5D<Impl>::Mdag(const FermionField &in, FermionField &out)
{
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
Dhop(in, out, DaggerYes); Dhop(in, out, DaggerYes);
return axpy_norm(out, mass, in, out); axpy(out, mass, in, out);
} }
template <class Impl> template <class Impl>
void ImprovedStaggeredFermion5D<Impl>::Meooe(const FermionField &in, FermionField &out) { void ImprovedStaggeredFermion5D<Impl>::Meooe(const FermionField &in, FermionField &out)
{
if (in.Checkerboard() == Odd) { if (in.Checkerboard() == Odd) {
DhopEO(in, out, DaggerNo); DhopEO(in, out, DaggerNo);
} else { } else {
@ -570,7 +573,8 @@ void ImprovedStaggeredFermion5D<Impl>::Meooe(const FermionField &in, FermionFiel
} }
} }
template <class Impl> template <class Impl>
void ImprovedStaggeredFermion5D<Impl>::MeooeDag(const FermionField &in, FermionField &out) { void ImprovedStaggeredFermion5D<Impl>::MeooeDag(const FermionField &in, FermionField &out)
{
if (in.Checkerboard() == Odd) { if (in.Checkerboard() == Odd) {
DhopEO(in, out, DaggerYes); DhopEO(in, out, DaggerYes);
} else { } else {
@ -579,27 +583,30 @@ void ImprovedStaggeredFermion5D<Impl>::MeooeDag(const FermionField &in, FermionF
} }
template <class Impl> template <class Impl>
void ImprovedStaggeredFermion5D<Impl>::Mooee(const FermionField &in, FermionField &out) { void ImprovedStaggeredFermion5D<Impl>::Mooee(const FermionField &in, FermionField &out)
{
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
typename FermionField::scalar_type scal(mass); typename FermionField::scalar_type scal(mass);
out = scal * in; out = scal * in;
} }
template <class Impl> template <class Impl>
void ImprovedStaggeredFermion5D<Impl>::MooeeDag(const FermionField &in, FermionField &out) { void ImprovedStaggeredFermion5D<Impl>::MooeeDag(const FermionField &in, FermionField &out)
{
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
Mooee(in, out); Mooee(in, out);
} }
template <class Impl> template <class Impl>
void ImprovedStaggeredFermion5D<Impl>::MooeeInv(const FermionField &in, FermionField &out) { void ImprovedStaggeredFermion5D<Impl>::MooeeInv(const FermionField &in, FermionField &out)
{
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
out = (1.0 / (mass)) * in; out = (1.0 / (mass)) * in;
} }
template <class Impl> template <class Impl>
void ImprovedStaggeredFermion5D<Impl>::MooeeInvDag(const FermionField &in, void ImprovedStaggeredFermion5D<Impl>::MooeeInvDag(const FermionField &in,FermionField &out)
FermionField &out) { {
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
MooeeInv(in, out); MooeeInv(in, out);
} }

44
Grid/qcd/action/fermion/implementation/ImprovedStaggeredFermionImplementation.h
View File

@ -171,21 +171,24 @@ void ImprovedStaggeredFermion<Impl>::ImportGauge(const GaugeField &_Uthin,const
///////////////////////////// /////////////////////////////
template <class Impl> template <class Impl>
RealD ImprovedStaggeredFermion<Impl>::M(const FermionField &in, FermionField &out) { void ImprovedStaggeredFermion<Impl>::M(const FermionField &in, FermionField &out)
{
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
Dhop(in, out, DaggerNo); Dhop(in, out, DaggerNo);
return axpy_norm(out, mass, in, out); axpy(out, mass, in, out);
} }
template <class Impl> template <class Impl>
RealD ImprovedStaggeredFermion<Impl>::Mdag(const FermionField &in, FermionField &out) { void ImprovedStaggeredFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
{
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
Dhop(in, out, DaggerYes); Dhop(in, out, DaggerYes);
return axpy_norm(out, mass, in, out); axpy(out, mass, in, out);
} }
template <class Impl> template <class Impl>
void ImprovedStaggeredFermion<Impl>::Meooe(const FermionField &in, FermionField &out) { void ImprovedStaggeredFermion<Impl>::Meooe(const FermionField &in, FermionField &out)
{
if (in.Checkerboard() == Odd) { if (in.Checkerboard() == Odd) {
DhopEO(in, out, DaggerNo); DhopEO(in, out, DaggerNo);
} else { } else {
@ -193,7 +196,8 @@ void ImprovedStaggeredFermion<Impl>::Meooe(const FermionField &in, FermionField
} }
} }
template <class Impl> template <class Impl>
void ImprovedStaggeredFermion<Impl>::MeooeDag(const FermionField &in, FermionField &out) { void ImprovedStaggeredFermion<Impl>::MeooeDag(const FermionField &in, FermionField &out)
{
if (in.Checkerboard() == Odd) { if (in.Checkerboard() == Odd) {
DhopEO(in, out, DaggerYes); DhopEO(in, out, DaggerYes);
} else { } else {
@ -202,27 +206,30 @@ void ImprovedStaggeredFermion<Impl>::MeooeDag(const FermionField &in, FermionFie
} }
template <class Impl> template <class Impl>
void ImprovedStaggeredFermion<Impl>::Mooee(const FermionField &in, FermionField &out) { void ImprovedStaggeredFermion<Impl>::Mooee(const FermionField &in, FermionField &out)
{
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
typename FermionField::scalar_type scal(mass); typename FermionField::scalar_type scal(mass);
out = scal * in; out = scal * in;
} }
template <class Impl> template <class Impl>
void ImprovedStaggeredFermion<Impl>::MooeeDag(const FermionField &in, FermionField &out) { void ImprovedStaggeredFermion<Impl>::MooeeDag(const FermionField &in, FermionField &out)
{
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
Mooee(in, out); Mooee(in, out);
} }
template <class Impl> template <class Impl>
void ImprovedStaggeredFermion<Impl>::MooeeInv(const FermionField &in, FermionField &out) { void ImprovedStaggeredFermion<Impl>::MooeeInv(const FermionField &in, FermionField &out)
{
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
out = (1.0 / (mass)) * in; out = (1.0 / (mass)) * in;
} }
template <class Impl> template <class Impl>
void ImprovedStaggeredFermion<Impl>::MooeeInvDag(const FermionField &in, void ImprovedStaggeredFermion<Impl>::MooeeInvDag(const FermionField &in,FermionField &out)
FermionField &out) { {
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
MooeeInv(in, out); MooeeInv(in, out);
} }
@ -234,7 +241,8 @@ void ImprovedStaggeredFermion<Impl>::MooeeInvDag(const FermionField &in,
template <class Impl> template <class Impl>
void ImprovedStaggeredFermion<Impl>::DerivInternal(StencilImpl &st, DoubledGaugeField &U, DoubledGaugeField &UUU, void ImprovedStaggeredFermion<Impl>::DerivInternal(StencilImpl &st, DoubledGaugeField &U, DoubledGaugeField &UUU,
GaugeField & mat, GaugeField & mat,
const FermionField &A, const FermionField &B, int dag) { const FermionField &A, const FermionField &B, int dag)
{
assert((dag == DaggerNo) || (dag == DaggerYes)); assert((dag == DaggerNo) || (dag == DaggerYes));
Compressor compressor; Compressor compressor;
@ -284,8 +292,8 @@ void ImprovedStaggeredFermion<Impl>::DerivInternal(StencilImpl &st, DoubledGauge
} }
template <class Impl> template <class Impl>
void ImprovedStaggeredFermion<Impl>::DhopDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) { void ImprovedStaggeredFermion<Impl>::DhopDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag)
{
conformable(U.Grid(), _grid); conformable(U.Grid(), _grid);
conformable(U.Grid(), V.Grid()); conformable(U.Grid(), V.Grid());
conformable(U.Grid(), mat.Grid()); conformable(U.Grid(), mat.Grid());
@ -296,8 +304,8 @@ void ImprovedStaggeredFermion<Impl>::DhopDeriv(GaugeField &mat, const FermionFie
} }
template <class Impl> template <class Impl>
void ImprovedStaggeredFermion<Impl>::DhopDerivOE(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) { void ImprovedStaggeredFermion<Impl>::DhopDerivOE(GaugeField &mat, const FermionField &U, const FermionField &V, int dag)
{
conformable(U.Grid(), _cbgrid); conformable(U.Grid(), _cbgrid);
conformable(U.Grid(), V.Grid()); conformable(U.Grid(), V.Grid());
conformable(U.Grid(), mat.Grid()); conformable(U.Grid(), mat.Grid());
@ -310,8 +318,8 @@ void ImprovedStaggeredFermion<Impl>::DhopDerivOE(GaugeField &mat, const FermionF
} }
template <class Impl> template <class Impl>
void ImprovedStaggeredFermion<Impl>::DhopDerivEO(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) { void ImprovedStaggeredFermion<Impl>::DhopDerivEO(GaugeField &mat, const FermionField &U, const FermionField &V, int dag)
{
conformable(U.Grid(), _cbgrid); conformable(U.Grid(), _cbgrid);
conformable(U.Grid(), V.Grid()); conformable(U.Grid(), V.Grid());
conformable(U.Grid(), mat.Grid()); conformable(U.Grid(), mat.Grid());

6
Grid/qcd/action/fermion/implementation/MobiusEOFAFermionImplementation.h
View File

@ -166,7 +166,7 @@ void MobiusEOFAFermion<Impl>::DtildeInv(const FermionField& psi, FermionField& c
/*****************************************************************************************************/ /*****************************************************************************************************/
template<class Impl> template<class Impl>
RealD MobiusEOFAFermion<Impl>::M(const FermionField& psi, FermionField& chi) void MobiusEOFAFermion<Impl>::M(const FermionField& psi, FermionField& chi)
{ {
FermionField Din(psi.Grid()); FermionField Din(psi.Grid());
@ -174,11 +174,10 @@ RealD MobiusEOFAFermion<Impl>::M(const FermionField& psi, FermionField& chi)
this->DW(Din, chi, DaggerNo); this->DW(Din, chi, DaggerNo);
axpby(chi, 1.0, 1.0, chi, psi); axpby(chi, 1.0, 1.0, chi, psi);
this->M5D(psi, chi); this->M5D(psi, chi);
return(norm2(chi));
} }
template<class Impl> template<class Impl>
RealD MobiusEOFAFermion<Impl>::Mdag(const FermionField& psi, FermionField& chi) void MobiusEOFAFermion<Impl>::Mdag(const FermionField& psi, FermionField& chi)
{ {
FermionField Din(psi.Grid()); FermionField Din(psi.Grid());
@ -186,7 +185,6 @@ RealD MobiusEOFAFermion<Impl>::Mdag(const FermionField& psi, FermionField& chi)
this->MeooeDag5D(Din, chi); this->MeooeDag5D(Din, chi);
this->M5Ddag(psi, chi); this->M5Ddag(psi, chi);
axpby(chi, 1.0, 1.0, chi, psi); axpby(chi, 1.0, 1.0, chi, psi);
return(norm2(chi));
} }
/******************************************************************** /********************************************************************

6
Grid/qcd/action/fermion/implementation/PartialFractionFermion5DImplementation.h
View File

@ -269,16 +269,14 @@ void PartialFractionFermion5D<Impl>::M_internal(const FermionField &psi, Fermi
} }
template<class Impl> template<class Impl>
RealD PartialFractionFermion5D<Impl>::M (const FermionField &in, FermionField &out) void PartialFractionFermion5D<Impl>::M (const FermionField &in, FermionField &out)
{ {
M_internal(in,out,DaggerNo); M_internal(in,out,DaggerNo);
return norm2(out);
} }
template<class Impl> template<class Impl>
RealD PartialFractionFermion5D<Impl>::Mdag (const FermionField &in, FermionField &out) void PartialFractionFermion5D<Impl>::Mdag (const FermionField &in, FermionField &out)
{ {
M_internal(in,out,DaggerYes); M_internal(in,out,DaggerYes);
return norm2(out);
} }
template<class Impl> template<class Impl>

6
Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h
View File

@ -35,7 +35,7 @@ NAMESPACE_BEGIN(Grid);
// *NOT* EO // *NOT* EO
template <class Impl> template <class Impl>
RealD WilsonCloverFermion<Impl>::M(const FermionField &in, FermionField &out) void WilsonCloverFermion<Impl>::M(const FermionField &in, FermionField &out)
{ {
FermionField temp(out.Grid()); FermionField temp(out.Grid());
@ -47,11 +47,10 @@ RealD WilsonCloverFermion<Impl>::M(const FermionField &in, FermionField &out)
Mooee(in, temp); Mooee(in, temp);
out += temp; out += temp;
return norm2(out);
} }
template <class Impl> template <class Impl>
RealD WilsonCloverFermion<Impl>::Mdag(const FermionField &in, FermionField &out) void WilsonCloverFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
{ {
FermionField temp(out.Grid()); FermionField temp(out.Grid());
@ -63,7 +62,6 @@ RealD WilsonCloverFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
MooeeDag(in, temp); MooeeDag(in, temp);
out += temp; out += temp;
return norm2(out);
} }
template <class Impl> template <class Impl>

52
Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h
View File

@ -102,21 +102,24 @@ void WilsonFermion<Impl>::ImportGauge(const GaugeField &_Umu)
///////////////////////////// /////////////////////////////
template <class Impl> template <class Impl>
RealD WilsonFermion<Impl>::M(const FermionField &in, FermionField &out) { void WilsonFermion<Impl>::M(const FermionField &in, FermionField &out)
{
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
Dhop(in, out, DaggerNo); Dhop(in, out, DaggerNo);
return axpy_norm(out, diag_mass, in, out); axpy(out, diag_mass, in, out);
} }
template <class Impl> template <class Impl>
RealD WilsonFermion<Impl>::Mdag(const FermionField &in, FermionField &out) { void WilsonFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
{
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
Dhop(in, out, DaggerYes); Dhop(in, out, DaggerYes);
return axpy_norm(out, diag_mass, in, out); axpy(out, diag_mass, in, out);
} }
template <class Impl> template <class Impl>
void WilsonFermion<Impl>::Meooe(const FermionField &in, FermionField &out) { void WilsonFermion<Impl>::Meooe(const FermionField &in, FermionField &out)
{
if (in.Checkerboard() == Odd) { if (in.Checkerboard() == Odd) {
DhopEO(in, out, DaggerNo); DhopEO(in, out, DaggerNo);
} else { } else {
@ -125,7 +128,8 @@ void WilsonFermion<Impl>::Meooe(const FermionField &in, FermionField &out) {
} }
template <class Impl> template <class Impl>
void WilsonFermion<Impl>::MeooeDag(const FermionField &in, FermionField &out) { void WilsonFermion<Impl>::MeooeDag(const FermionField &in, FermionField &out)
{
if (in.Checkerboard() == Odd) { if (in.Checkerboard() == Odd) {
DhopEO(in, out, DaggerYes); DhopEO(in, out, DaggerYes);
} else { } else {
@ -134,26 +138,30 @@ void WilsonFermion<Impl>::MeooeDag(const FermionField &in, FermionField &out) {
} }
template <class Impl> template <class Impl>
void WilsonFermion<Impl>::Mooee(const FermionField &in, FermionField &out) { void WilsonFermion<Impl>::Mooee(const FermionField &in, FermionField &out)
{
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
typename FermionField::scalar_type scal(diag_mass); typename FermionField::scalar_type scal(diag_mass);
out = scal * in; out = scal * in;
} }
template <class Impl> template <class Impl>
void WilsonFermion<Impl>::MooeeDag(const FermionField &in, FermionField &out) { void WilsonFermion<Impl>::MooeeDag(const FermionField &in, FermionField &out)
{
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
Mooee(in, out); Mooee(in, out);
} }
template<class Impl> template<class Impl>
void WilsonFermion<Impl>::MooeeInv(const FermionField &in, FermionField &out) { void WilsonFermion<Impl>::MooeeInv(const FermionField &in, FermionField &out)
{
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
out = (1.0/(diag_mass))*in; out = (1.0/(diag_mass))*in;
} }
template<class Impl> template<class Impl>
void WilsonFermion<Impl>::MooeeInvDag(const FermionField &in, FermionField &out) { void WilsonFermion<Impl>::MooeeInvDag(const FermionField &in, FermionField &out)
{
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
MooeeInv(in,out); MooeeInv(in,out);
} }
@ -249,7 +257,8 @@ void WilsonFermion<Impl>::DerivInternal(StencilImpl &st, DoubledGaugeField &U,
} }
template <class Impl> template <class Impl>
void WilsonFermion<Impl>::DhopDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) { void WilsonFermion<Impl>::DhopDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag)
{
conformable(U.Grid(), _grid); conformable(U.Grid(), _grid);
conformable(U.Grid(), V.Grid()); conformable(U.Grid(), V.Grid());
conformable(U.Grid(), mat.Grid()); conformable(U.Grid(), mat.Grid());
@ -260,7 +269,8 @@ void WilsonFermion<Impl>::DhopDeriv(GaugeField &mat, const FermionField &U, cons
} }
template <class Impl> template <class Impl>
void WilsonFermion<Impl>::DhopDerivOE(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) { void WilsonFermion<Impl>::DhopDerivOE(GaugeField &mat, const FermionField &U, const FermionField &V, int dag)
{
conformable(U.Grid(), _cbgrid); conformable(U.Grid(), _cbgrid);
conformable(U.Grid(), V.Grid()); conformable(U.Grid(), V.Grid());
//conformable(U.Grid(), mat.Grid()); not general, leaving as a comment (Guido) //conformable(U.Grid(), mat.Grid()); not general, leaving as a comment (Guido)
@ -274,7 +284,8 @@ void WilsonFermion<Impl>::DhopDerivOE(GaugeField &mat, const FermionField &U, co
} }
template <class Impl> template <class Impl>
void WilsonFermion<Impl>::DhopDerivEO(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) { void WilsonFermion<Impl>::DhopDerivEO(GaugeField &mat, const FermionField &U, const FermionField &V, int dag)
{
conformable(U.Grid(), _cbgrid); conformable(U.Grid(), _cbgrid);
conformable(U.Grid(), V.Grid()); conformable(U.Grid(), V.Grid());
//conformable(U.Grid(), mat.Grid()); //conformable(U.Grid(), mat.Grid());
@ -287,7 +298,8 @@ void WilsonFermion<Impl>::DhopDerivEO(GaugeField &mat, const FermionField &U, co
} }
template <class Impl> template <class Impl>
void WilsonFermion<Impl>::Dhop(const FermionField &in, FermionField &out, int dag) { void WilsonFermion<Impl>::Dhop(const FermionField &in, FermionField &out, int dag)
{
conformable(in.Grid(), _grid); // verifies full grid conformable(in.Grid(), _grid); // verifies full grid
conformable(in.Grid(), out.Grid()); conformable(in.Grid(), out.Grid());
@ -297,7 +309,8 @@ void WilsonFermion<Impl>::Dhop(const FermionField &in, FermionField &out, int da
} }
template <class Impl> template <class Impl>
void WilsonFermion<Impl>::DhopOE(const FermionField &in, FermionField &out, int dag) { void WilsonFermion<Impl>::DhopOE(const FermionField &in, FermionField &out, int dag)
{
conformable(in.Grid(), _cbgrid); // verifies half grid conformable(in.Grid(), _cbgrid); // verifies half grid
conformable(in.Grid(), out.Grid()); // drops the cb check conformable(in.Grid(), out.Grid()); // drops the cb check
@ -308,7 +321,8 @@ void WilsonFermion<Impl>::DhopOE(const FermionField &in, FermionField &out, int
} }
template <class Impl> template <class Impl>
void WilsonFermion<Impl>::DhopEO(const FermionField &in, FermionField &out,int dag) { void WilsonFermion<Impl>::DhopEO(const FermionField &in, FermionField &out,int dag)
{
conformable(in.Grid(), _cbgrid); // verifies half grid conformable(in.Grid(), _cbgrid); // verifies half grid
conformable(in.Grid(), out.Grid()); // drops the cb check conformable(in.Grid(), out.Grid()); // drops the cb check
@ -386,7 +400,8 @@ template <class Impl>
void WilsonFermion<Impl>::DhopInternalOverlappedComms(StencilImpl &st, LebesgueOrder &lo, void WilsonFermion<Impl>::DhopInternalOverlappedComms(StencilImpl &st, LebesgueOrder &lo,
DoubledGaugeField &U, DoubledGaugeField &U,
const FermionField &in, const FermionField &in,
FermionField &out, int dag) { FermionField &out, int dag)
{
assert((dag == DaggerNo) || (dag == DaggerYes)); assert((dag == DaggerNo) || (dag == DaggerYes));
Compressor compressor(dag); Compressor compressor(dag);
@ -436,7 +451,8 @@ template <class Impl>
void WilsonFermion<Impl>::DhopInternalSerial(StencilImpl &st, LebesgueOrder &lo, void WilsonFermion<Impl>::DhopInternalSerial(StencilImpl &st, LebesgueOrder &lo,
DoubledGaugeField &U, DoubledGaugeField &U,
const FermionField &in, const FermionField &in,
FermionField &out, int dag) { FermionField &out, int dag)
{
assert((dag == DaggerNo) || (dag == DaggerYes)); assert((dag == DaggerNo) || (dag == DaggerYes));
Compressor compressor(dag); Compressor compressor(dag);
st.HaloExchange(in, compressor); st.HaloExchange(in, compressor);

8
Grid/qcd/modules/Registration.h
View File

@ -80,10 +80,10 @@ static Registrar<OneFlavourRatioEOFModule<FermionImplementationPolicy>,
static Registrar< ConjugateGradientModule<WilsonFermionR::FermionField>, static Registrar< ConjugateGradientModule<WilsonFermionR::FermionField>,
HMC_SolverModuleFactory<solver_string, WilsonFermionR::FermionField, Serialiser> > __CGWFmodXMLInit("ConjugateGradient"); HMC_SolverModuleFactory<solver_string, WilsonFermionR::FermionField, Serialiser> > __CGWFmodXMLInit("ConjugateGradient");
static Registrar< BiCGSTABModule<WilsonFermionR::FermionField>, //static Registrar< BiCGSTABModule<WilsonFermionR::FermionField>,
HMC_SolverModuleFactory<solver_string, WilsonFermionR::FermionField, Serialiser> > __CGWFmodXMLInit("BiCGSTAB"); // HMC_SolverModuleFactory<solver_string, WilsonFermionR::FermionField, Serialiser> > __CGWFmodXMLInit("BiCGSTAB");
static Registrar< ConjugateResidualModule<WilsonFermionR::FermionField>, //static Registrar< ConjugateResidualModule<WilsonFermionR::FermionField>,
HMC_SolverModuleFactory<solver_string, WilsonFermionR::FermionField, Serialiser> > __CRWFmodXMLInit("ConjugateResidual"); // HMC_SolverModuleFactory<solver_string, WilsonFermionR::FermionField, Serialiser> > __CRWFmodXMLInit("ConjugateResidual");
// add the staggered, scalar versions here // add the staggered, scalar versions here

2
Grid/qcd/smearing/GaugeConfiguration.h
View File

@ -49,7 +49,7 @@ public:
private: private:
const unsigned int smearingLevels; const unsigned int smearingLevels;
Smear_Stout<Gimpl> StoutSmearing; Smear_Stout<Gimpl> &StoutSmearing;
std::vector<GaugeField> SmearedSet; std::vector<GaugeField> SmearedSet;
// Member functions // Member functions

40
Grid/qcd/utils/CovariantCshift.h
View File

@ -52,6 +52,26 @@ namespace PeriodicBC {
tmp = adj(Link)*field; tmp = adj(Link)*field;
return Cshift(tmp,mu,-1);// moves towards positive mu return Cshift(tmp,mu,-1);// moves towards positive mu
} }
template<class gauge,typename Op, typename T1> auto
CovShiftForward(const Lattice<gauge> &Link,
int mu,
const LatticeUnaryExpression<Op,T1> &expr)
-> Lattice<decltype(expr.op.func(eval(0, expr.arg1)))>
{
Lattice<decltype(expr.op.func(eval(0, expr.arg1)))> arg(expr);
return CovShiftForward(Link,mu,arg);
}
template<class gauge,typename Op, typename T1> auto
CovShiftBackward(const Lattice<gauge> &Link,
int mu,
const LatticeUnaryExpression<Op,T1> &expr)
-> Lattice<decltype(expr.op.func(eval(0, expr.arg1)))>
{
Lattice<decltype(expr.op.func(eval(0, expr.arg1)))> arg(expr);
return CovShiftForward(Link,mu,arg);
}
} }
@ -122,6 +142,26 @@ namespace ConjugateBC {
return Cshift(tmp,mu,-1);// moves towards positive mu return Cshift(tmp,mu,-1);// moves towards positive mu
} }
template<class gauge,typename Op, typename T1> auto
CovShiftForward(const Lattice<gauge> &Link,
int mu,
const LatticeUnaryExpression<Op,T1> &expr)
-> Lattice<decltype(expr.op.func(eval(0, expr.arg1)))>
{
Lattice<decltype(expr.op.func(eval(0, expr.arg1)))> arg(expr);
return CovShiftForward(Link,mu,arg);
}
template<class gauge,typename Op, typename T1> auto
CovShiftBackward(const Lattice<gauge> &Link,
int mu,
const LatticeUnaryExpression<Op,T1> &expr)
-> Lattice<decltype(expr.op.func(eval(0, expr.arg1)))>
{
Lattice<decltype(expr.op.func(eval(0, expr.arg1)))> arg(expr);
return CovShiftForward(Link,mu,arg);
}
} }