portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-09 21:50:45 +01:00
Code Releases Activity

Merge pull request #5 from paboyle/develop

Browse Source 
Sync upstream
This commit is contained in:
Christoph Lehner 2020-05-13 09:02:56 +02:00 committed by GitHub
commit 32fbdf4fb1
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
45 changed files with 552 additions and 541 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)); });
} }
} }

628
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);
} }
@ -208,339 +207,305 @@ public:
} }
}; };
////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////
// Even Odd Schur decomp operators; there are several // Even Odd Schur decomp operators; there are several
// ways to introduce the even odd checkerboarding // ways to introduce the even odd checkerboarding
////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////
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){
out.Checkerboard() = in.Checkerboard();
MpcDagMpc(in,out,n1,n2);
}
virtual void HermOp(const Field &in, Field &out){
RealD n1,n2;
HermOpAndNorm(in,out,n1,n2);
}
void Op (const Field &in, Field &out){
Mpc(in,out);
}
void AdjOp (const Field &in, Field &out){
MpcDag(in,out);
}
// Support for coarsening to a multigrid
void OpDiag (const Field &in, Field &out) {
assert(0); // must coarsen the unpreconditioned system
}
void OpDir (const Field &in, Field &out,int dir,int disp) {
assert(0);
}
void OpDirAll (const Field &in, std::vector<Field> &out){
assert(0);
};
};
template<class Matrix,class Field>
class SchurDiagMooeeOperator : public SchurOperatorBase<Field> {
public:
Matrix &_Mat;
SchurDiagMooeeOperator (Matrix &Mat): _Mat(Mat){};
virtual RealD Mpc (const Field &in, Field &out) {
Field tmp(in.Grid());
tmp.Checkerboard() = !in.Checkerboard();
_Mat.Meooe(in,tmp);
_Mat.MooeeInv(tmp,out);
_Mat.Meooe(out,tmp);
_Mat.Mooee(in,out);
return axpy_norm(out,-1.0,tmp,out);
}
virtual RealD MpcDag (const Field &in, Field &out){
Field tmp(in.Grid());
_Mat.MeooeDag(in,tmp);
_Mat.MooeeInvDag(tmp,out);
_Mat.MeooeDag(out,tmp);
_Mat.MooeeDag(in,out);
return axpy_norm(out,-1.0,tmp,out);
}
};
template<class Matrix,class Field>
class SchurDiagOneOperator : public SchurOperatorBase<Field> {
protected:
Matrix &_Mat;
public:
SchurDiagOneOperator (Matrix &Mat): _Mat(Mat){};
virtual RealD Mpc (const Field &in, Field &out) {
Field tmp(in.Grid());
_Mat.Meooe(in,out);
_Mat.MooeeInv(out,tmp);
_Mat.Meooe(tmp,out);
_Mat.MooeeInv(out,tmp);
return axpy_norm(out,-1.0,tmp,in);
}
virtual RealD MpcDag (const Field &in, Field &out){
Field tmp(in.Grid());
_Mat.MooeeInvDag(in,out);
_Mat.MeooeDag(out,tmp);
_Mat.MooeeInvDag(tmp,out);
_Mat.MeooeDag(out,tmp);
return axpy_norm(out,-1.0,tmp,in);
}
};
template<class Matrix,class Field>
class SchurDiagTwoOperator : public SchurOperatorBase<Field> {
protected:
Matrix &_Mat;
public:
SchurDiagTwoOperator (Matrix &Mat): _Mat(Mat){};
virtual RealD Mpc (const Field &in, Field &out) {
Field tmp(in.Grid());
_Mat.MooeeInv(in,out);
_Mat.Meooe(out,tmp);
_Mat.MooeeInv(tmp,out);
_Mat.Meooe(out,tmp);
return axpy_norm(out,-1.0,tmp,in);
}
virtual RealD MpcDag (const Field &in, Field &out){
Field tmp(in.Grid());
_Mat.MeooeDag(in,out);
_Mat.MooeeInvDag(out,tmp);
_Mat.MeooeDag(tmp,out);
_Mat.MooeeInvDag(out,tmp);
return axpy_norm(out,-1.0,tmp,in);
}
};
template<class Field>
class NonHermitianSchurOperatorBase : public LinearOperatorBase<Field>
{
public:
virtual RealD Mpc (const Field& in, Field& out) = 0;
virtual RealD MpcDag (const Field& in, Field& out) = 0;
virtual void MpcDagMpc(const Field& in, Field& out, RealD& ni, RealD& no) {
Field tmp(in.Grid());
tmp.Checkerboard() = in.Checkerboard();
ni = Mpc(in,tmp);
no = MpcDag(tmp,out);
}
virtual void HermOpAndNorm(const Field& in, Field& out, RealD& n1, RealD& n2) {
assert(0);
}
virtual void HermOp(const Field& in, Field& out) {
assert(0);
}
void Op(const Field& in, Field& out) {
Mpc(in, out);
}
void AdjOp(const Field& in, Field& out) {
MpcDag(in, out);
}
// Support for coarsening to a multigrid
void OpDiag(const Field& in, Field& out) {
assert(0); // must coarsen the unpreconditioned system
}
void OpDir(const Field& in, Field& out, int dir, int disp) {
assert(0);
}
void OpDirAll(const Field& in, std::vector<Field>& out){
assert(0);
};
};
template<class Matrix, class Field>
class NonHermitianSchurDiagMooeeOperator : public NonHermitianSchurOperatorBase<Field>
{
public:
Matrix& _Mat;
NonHermitianSchurDiagMooeeOperator(Matrix& Mat): _Mat(Mat){};
virtual RealD Mpc(const Field& in, Field& out) {
Field tmp(in.Grid());
tmp.Checkerboard() = !in.Checkerboard();
_Mat.Meooe(in, tmp);
_Mat.MooeeInv(tmp, out);
_Mat.Meooe(out, tmp);
_Mat.Mooee(in, out);
return axpy_norm(out, -1.0, tmp, out);
}
virtual RealD MpcDag(const Field& in, Field& out) {
Field tmp(in.Grid());
_Mat.MeooeDag(in, tmp);
_Mat.MooeeInvDag(tmp, out);
_Mat.MeooeDag(out, tmp);
_Mat.MooeeDag(in, out);
return axpy_norm(out, -1.0, tmp, out);
}
};
template<class Matrix,class Field>
class NonHermitianSchurDiagOneOperator : public NonHermitianSchurOperatorBase<Field>
{
protected:
Matrix &_Mat;
public:
NonHermitianSchurDiagOneOperator (Matrix& Mat): _Mat(Mat){};
virtual RealD Mpc(const Field& in, Field& out) {
Field tmp(in.Grid());
_Mat.Meooe(in, out);
_Mat.MooeeInv(out, tmp);
_Mat.Meooe(tmp, out);
_Mat.MooeeInv(out, tmp);
return axpy_norm(out, -1.0, tmp, in);
}
virtual RealD MpcDag(const Field& in, Field& out) {
Field tmp(in.Grid());
_Mat.MooeeInvDag(in, out);
_Mat.MeooeDag(out, tmp);
_Mat.MooeeInvDag(tmp, out);
_Mat.MeooeDag(out, tmp);
return axpy_norm(out, -1.0, tmp, in);
}
};
template<class Matrix, class Field>
class NonHermitianSchurDiagTwoOperator : public NonHermitianSchurOperatorBase<Field>
{
protected:
Matrix& _Mat;
public:
NonHermitianSchurDiagTwoOperator(Matrix& Mat): _Mat(Mat){};
virtual RealD Mpc(const Field& in, Field& out) {
Field tmp(in.Grid());
_Mat.MooeeInv(in, out);
_Mat.Meooe(out, tmp);
_Mat.MooeeInv(tmp, out);
_Mat.Meooe(out, tmp);
return axpy_norm(out, -1.0, tmp, in);
}
virtual RealD MpcDag(const Field& in, Field& out) {
Field tmp(in.Grid());
_Mat.MeooeDag(in, out);
_Mat.MooeeInvDag(out, tmp);
_Mat.MeooeDag(tmp, out);
_Mat.MooeeInvDag(out, tmp);
return axpy_norm(out, -1.0, tmp, in);
}
};
///////////////////////////////////////////////////////////////////////////////////////////////////
// Left handed Moo^-1 ; (Moo - Moe Mee^-1 Meo) psi = eta --> ( 1 - Moo^-1 Moe Mee^-1 Meo ) psi = Moo^-1 eta
// Right handed Moo^-1 ; (Moo - Moe Mee^-1 Meo) Moo^-1 Moo psi = eta --> ( 1 - Moe Mee^-1 Meo Moo^-1) phi=eta ; psi = Moo^-1 phi
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class Matrix,class Field> using SchurDiagOneRH = SchurDiagTwoOperator<Matrix,Field> ;
template<class Matrix,class Field> using SchurDiagOneLH = SchurDiagOneOperator<Matrix,Field> ;
///////////////////////////////////////////////////////////////////////////////////////////////////
// Staggered use
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class Matrix,class Field>
class SchurStaggeredOperator : public SchurOperatorBase<Field> {
protected:
Matrix &_Mat;
Field tmp;
RealD mass;
double tMpc;
double tIP;
double tMeo;
double taxpby_norm;
uint64_t ncall;
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())
{
assert( _Mat.isTrivialEE() );
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++; out.Checkerboard() = in.Checkerboard();
tMpc-=usecond(); MpcDagMpc(in,out);
n2 = Mpc(in,out); ComplexD dot= innerProduct(in,out);
tMpc+=usecond(); n1=real(dot);
tIP-=usecond(); n2=norm2(out);
ComplexD dot= innerProduct(in,out);
tIP+=usecond();
n1 = real(dot);
} }
virtual void HermOp(const Field &in, Field &out){ virtual void HermOp(const Field &in, Field &out){
ncall++; out.Checkerboard() = in.Checkerboard();
tMpc-=usecond(); MpcDagMpc(in,out);
_Mat.Meooe(in,out);
_Mat.Meooe(out,tmp);
tMpc+=usecond();
taxpby_norm-=usecond();
axpby(out,-1.0,mass*mass,tmp,in);
taxpby_norm+=usecond();
} }
virtual RealD Mpc (const Field &in, Field &out) void Op (const Field &in, Field &out){
{ Mpc(in,out);
}
void AdjOp (const Field &in, Field &out){
MpcDag(in,out);
}
// Support for coarsening to a multigrid
void OpDiag (const Field &in, Field &out) {
assert(0); // must coarsen the unpreconditioned system
}
void OpDir (const Field &in, Field &out,int dir,int disp) {
assert(0);
}
void OpDirAll (const Field &in, std::vector<Field> &out){
assert(0);
};
};
template<class Matrix,class Field>
class SchurDiagMooeeOperator : public SchurOperatorBase<Field> {
public:
Matrix &_Mat;
SchurDiagMooeeOperator (Matrix &Mat): _Mat(Mat){};
virtual void Mpc (const Field &in, Field &out) {
Field tmp(in.Grid());
tmp.Checkerboard() = !in.Checkerboard();
_Mat.Meooe(in,tmp);
_Mat.MooeeInv(tmp,out);
_Mat.Meooe(out,tmp);
_Mat.Mooee(in,out);
axpy(out,-1.0,tmp,out);
}
virtual void MpcDag (const Field &in, Field &out){
Field tmp(in.Grid());
_Mat.MeooeDag(in,tmp);
_Mat.MooeeInvDag(tmp,out);
_Mat.MeooeDag(out,tmp);
_Mat.MooeeDag(in,out);
axpy(out,-1.0,tmp,out);
}
};
template<class Matrix,class Field>
class SchurDiagOneOperator : public SchurOperatorBase<Field> {
protected:
Matrix &_Mat;
public:
SchurDiagOneOperator (Matrix &Mat): _Mat(Mat){};
virtual void Mpc (const Field &in, Field &out) {
Field tmp(in.Grid());
_Mat.Meooe(in,out);
_Mat.MooeeInv(out,tmp);
_Mat.Meooe(tmp,out);
_Mat.MooeeInv(out,tmp);
axpy(out,-1.0,tmp,in);
}
virtual void MpcDag (const Field &in, Field &out){
Field tmp(in.Grid());
_Mat.MooeeInvDag(in,out);
_Mat.MeooeDag(out,tmp);
_Mat.MooeeInvDag(tmp,out);
_Mat.MeooeDag(out,tmp);
axpy(out,-1.0,tmp,in);
}
};
template<class Matrix,class Field>
class SchurDiagTwoOperator : public SchurOperatorBase<Field> {
protected:
Matrix &_Mat;
public:
SchurDiagTwoOperator (Matrix &Mat): _Mat(Mat){};
virtual void Mpc (const Field &in, Field &out) {
Field tmp(in.Grid());
_Mat.MooeeInv(in,out);
_Mat.Meooe(out,tmp);
_Mat.MooeeInv(tmp,out);
_Mat.Meooe(out,tmp);
axpy(out,-1.0,tmp,in);
}
virtual void MpcDag (const Field &in, Field &out){
Field tmp(in.Grid());
_Mat.MeooeDag(in,out);
_Mat.MooeeInvDag(out,tmp);
_Mat.MeooeDag(tmp,out);
_Mat.MooeeInvDag(out,tmp);
axpy(out,-1.0,tmp,in);
}
};
template<class Field>
class NonHermitianSchurOperatorBase : public LinearOperatorBase<Field>
{
public:
virtual void Mpc (const Field& in, Field& out) = 0;
virtual void MpcDag (const Field& in, Field& out) = 0;
virtual void MpcDagMpc(const Field& in, Field& out) {
Field tmp(in.Grid());
tmp.Checkerboard() = in.Checkerboard();
Mpc(in,tmp);
MpcDag(tmp,out);
}
virtual void HermOpAndNorm(const Field& in, Field& out, RealD& n1, RealD& n2) {
assert(0);
}
virtual void HermOp(const Field& in, Field& out) {
assert(0);
}
void Op(const Field& in, Field& out) {
Mpc(in, out);
}
void AdjOp(const Field& in, Field& out) {
MpcDag(in, out);
}
// Support for coarsening to a multigrid
void OpDiag(const Field& in, Field& out) {
assert(0); // must coarsen the unpreconditioned system
}
void OpDir(const Field& in, Field& out, int dir, int disp) {
assert(0);
}
void OpDirAll(const Field& in, std::vector<Field>& out){
assert(0);
};
};
template<class Matrix, class Field>
class NonHermitianSchurDiagMooeeOperator : public NonHermitianSchurOperatorBase<Field>
{
public:
Matrix& _Mat;
NonHermitianSchurDiagMooeeOperator(Matrix& Mat): _Mat(Mat){};
virtual void Mpc(const Field& in, Field& out) {
Field tmp(in.Grid());
tmp.Checkerboard() = !in.Checkerboard();
_Mat.Meooe(in, tmp);
_Mat.MooeeInv(tmp, out);
_Mat.Meooe(out, tmp);
_Mat.Mooee(in, out);
axpy(out, -1.0, tmp, out);
}
virtual void MpcDag(const Field& in, Field& out) {
Field tmp(in.Grid());
_Mat.MeooeDag(in, tmp);
_Mat.MooeeInvDag(tmp, out);
_Mat.MeooeDag(out, tmp);
_Mat.MooeeDag(in, out);
axpy(out, -1.0, tmp, out);
}
};
template<class Matrix,class Field>
class NonHermitianSchurDiagOneOperator : public NonHermitianSchurOperatorBase<Field>
{
protected:
Matrix &_Mat;
public:
NonHermitianSchurDiagOneOperator (Matrix& Mat): _Mat(Mat){};
virtual void Mpc(const Field& in, Field& out) {
Field tmp(in.Grid());
_Mat.Meooe(in, out);
_Mat.MooeeInv(out, tmp);
_Mat.Meooe(tmp, out);
_Mat.MooeeInv(out, tmp);
axpy(out, -1.0, tmp, in);
}
virtual void MpcDag(const Field& in, Field& out) {
Field tmp(in.Grid());
_Mat.MooeeInvDag(in, out);
_Mat.MeooeDag(out, tmp);
_Mat.MooeeInvDag(tmp, out);
_Mat.MeooeDag(out, tmp);
axpy(out, -1.0, tmp, in);
}
};
template<class Matrix, class Field>
class NonHermitianSchurDiagTwoOperator : public NonHermitianSchurOperatorBase<Field>
{
protected:
Matrix& _Mat;
public:
NonHermitianSchurDiagTwoOperator(Matrix& Mat): _Mat(Mat){};
virtual void Mpc(const Field& in, Field& out) {
Field tmp(in.Grid());
_Mat.MooeeInv(in, out);
_Mat.Meooe(out, tmp);
_Mat.MooeeInv(tmp, out);
_Mat.Meooe(out, tmp);
axpy(out, -1.0, tmp, in);
}
virtual void MpcDag(const Field& in, Field& out) {
Field tmp(in.Grid());
_Mat.MeooeDag(in, out);
_Mat.MooeeInvDag(out, tmp);
_Mat.MeooeDag(tmp, out);
_Mat.MooeeInvDag(out, tmp);
axpy(out, -1.0, tmp, in);
}
};
///////////////////////////////////////////////////////////////////////////////////////////////////
// Left handed Moo^-1 ; (Moo - Moe Mee^-1 Meo) psi = eta --> ( 1 - Moo^-1 Moe Mee^-1 Meo ) psi = Moo^-1 eta
// Right handed Moo^-1 ; (Moo - Moe Mee^-1 Meo) Moo^-1 Moo psi = eta --> ( 1 - Moe Mee^-1 Meo Moo^-1) phi=eta ; psi = Moo^-1 phi
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class Matrix,class Field> using SchurDiagOneRH = SchurDiagTwoOperator<Matrix,Field> ;
template<class Matrix,class Field> using SchurDiagOneLH = SchurDiagOneOperator<Matrix,Field> ;
///////////////////////////////////////////////////////////////////////////////////////////////////
// Staggered use
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class Matrix,class Field>
class SchurStaggeredOperator : public SchurOperatorBase<Field> {
protected:
Matrix &_Mat;
Field tmp;
RealD mass;
public:
SchurStaggeredOperator (Matrix &Mat): _Mat(Mat), tmp(_Mat.RedBlackGrid())
{
assert( _Mat.isTrivialEE() );
mass = _Mat.Mass();
}
virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
Mpc(in,out);
ComplexD dot= innerProduct(in,out);
n1 = real(dot);
n2 =0.0;
}
virtual void HermOp(const Field &in, Field &out){
Mpc(in,out);
// _Mat.Meooe(in,out);
// _Mat.Meooe(out,tmp);
// axpby(out,-1.0,mass*mass,tmp,in);
}
virtual void Mpc (const Field &in, Field &out)
{
Field tmp(in.Grid()); Field tmp(in.Grid());
Field tmp2(in.Grid()); Field tmp2(in.Grid());
// _Mat.Mooee(in,out);
// _Mat.Mooee(out,tmp);
// std::cout << GridLogIterative << " HermOp.Mpc "<<std::endl;
_Mat.Mooee(in,out);
_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;

25
Grid/cshift/Cshift.h
View File

@ -49,4 +49,29 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
#ifdef GRID_COMMS_SHMEM #ifdef GRID_COMMS_SHMEM
#include <Grid/cshift/Cshift_mpi.h> // uses same implementation of communicator #include <Grid/cshift/Cshift_mpi.h> // uses same implementation of communicator
#endif #endif
NAMESPACE_BEGIN(Grid);
template<typename Op, typename T1>
auto Cshift(const LatticeUnaryExpression<Op,T1> &expr,int dim,int shift)
-> Lattice<decltype(expr.op.func(eval(0, expr.arg1)))>
{
return Cshift(closure(expr),dim,shift);
}
template <class Op, class T1, class T2>
auto Cshift(const LatticeBinaryExpression<Op,T1,T2> &expr,int dim,int shift)
-> Lattice<decltype(expr.op.func(eval(0, expr.arg1),eval(0, expr.arg2)))>
{
return Cshift(closure(expr),dim,shift);
}
template <class Op, class T1, class T2, class T3>
auto Cshift(const LatticeTrinaryExpression<Op,T1,T2,T3> &expr,int dim,int shift)
-> Lattice<decltype(expr.op.func(eval(0, expr.arg1),
eval(0, expr.arg2),
eval(0, expr.arg3)))>
{
return Cshift(closure(expr),dim,shift);
}
NAMESPACE_END(Grid);
#endif #endif

2
Grid/lattice/Lattice.h
View File

@ -31,7 +31,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
#include <Grid/lattice/Lattice_ET.h> #include <Grid/lattice/Lattice_ET.h>
#include <Grid/lattice/Lattice_arith.h> #include <Grid/lattice/Lattice_arith.h>
#include <Grid/lattice/Lattice_trace.h> #include <Grid/lattice/Lattice_trace.h>
//#include <Grid/lattice/Lattice_transpose.h> #include <Grid/lattice/Lattice_transpose.h>
#include <Grid/lattice/Lattice_local.h> #include <Grid/lattice/Lattice_local.h>
#include <Grid/lattice/Lattice_reduction.h> #include <Grid/lattice/Lattice_reduction.h>
#include <Grid/lattice/Lattice_peekpoke.h> #include <Grid/lattice/Lattice_peekpoke.h>

4
Grid/lattice/Lattice_transpose.h
View File

@ -38,6 +38,7 @@ NAMESPACE_BEGIN(Grid);
//////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////////
// Transpose // Transpose
//////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////////
/*
template<class vobj> template<class vobj>
inline Lattice<vobj> transpose(const Lattice<vobj> &lhs){ inline Lattice<vobj> transpose(const Lattice<vobj> &lhs){
Lattice<vobj> ret(lhs.Grid()); Lattice<vobj> ret(lhs.Grid());
@ -48,7 +49,8 @@ inline Lattice<vobj> transpose(const Lattice<vobj> &lhs){
}); });
return ret; return ret;
}; };
*/
//////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////////
// Index level dependent transpose // Index level dependent transpose
//////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////////

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);

9
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,15 +86,14 @@ 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
virtual void MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _m,std::vector<double> twist) { assert(0);}; virtual void MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _m,std::vector<double> twist) { assert(0);};
virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector<Complex> boundary,std::vector<double> twist) virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector<Complex> boundary,std::vector<double> twist)
{ {
FFT theFFT((GridCartesian *) in.Grid()); FFT theFFT((GridCartesian *) in.Grid());

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

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);
}
} }

6
tests/IO/Test_openqcd_io.cc
View File

@ -28,13 +28,12 @@ See the full license in the file "LICENSE" in the top level distribution directo
#include <Grid/Grid.h> #include <Grid/Grid.h>
#if defined(GRID_COMMS_NONE)
#error This test requires Grid compiled with MPI
#endif
using namespace Grid; using namespace Grid;
int main(int argc, char** argv) { int main(int argc, char** argv) {
#if !defined(GRID_COMMS_NONE)
Grid_init(&argc, &argv); Grid_init(&argc, &argv);
auto simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd()); auto simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
@ -81,4 +80,5 @@ int main(int argc, char** argv) {
// clang-format on // clang-format on
Grid_finalize(); Grid_finalize();
#endif
} }

8
tests/core/Test_contfrac_even_odd.cc
View File

@ -238,11 +238,11 @@ void TestWhat(What & Ddwf,
RealD t1,t2; RealD t1,t2;
SchurDiagMooeeOperator<What,LatticeFermion> HermOpEO(Ddwf); SchurDiagMooeeOperator<What,LatticeFermion> HermOpEO(Ddwf);
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2); HermOpEO.MpcDagMpc(chi_e,dchi_e);
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2); HermOpEO.MpcDagMpc(chi_o,dchi_o);
HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2); HermOpEO.MpcDagMpc(phi_e,dphi_e);
HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2); HermOpEO.MpcDagMpc(phi_o,dphi_o);
pDce = innerProduct(phi_e,dchi_e); pDce = innerProduct(phi_e,dchi_e);
pDco = innerProduct(phi_o,dchi_o); pDco = innerProduct(phi_o,dchi_o);

8
tests/core/Test_dwf_eofa_even_odd.cc
View File

@ -218,11 +218,11 @@ int main (int argc, char ** argv)
RealD t1,t2; RealD t1,t2;
SchurDiagMooeeOperator<DomainWallEOFAFermionR,LatticeFermion> HermOpEO(Ddwf); SchurDiagMooeeOperator<DomainWallEOFAFermionR,LatticeFermion> HermOpEO(Ddwf);
HermOpEO.MpcDagMpc(chi_e, dchi_e, t1, t2); HermOpEO.MpcDagMpc(chi_e, dchi_e);
HermOpEO.MpcDagMpc(chi_o, dchi_o, t1, t2); HermOpEO.MpcDagMpc(chi_o, dchi_o);
HermOpEO.MpcDagMpc(phi_e, dphi_e, t1, t2); HermOpEO.MpcDagMpc(phi_e, dphi_e);
HermOpEO.MpcDagMpc(phi_o, dphi_o, t1, t2); HermOpEO.MpcDagMpc(phi_o, dphi_o);
pDce = innerProduct(phi_e, dchi_e); pDce = innerProduct(phi_e, dchi_e);
pDco = innerProduct(phi_o, dchi_o); pDco = innerProduct(phi_o, dchi_o);

8
tests/core/Test_dwf_even_odd.cc
View File

@ -216,11 +216,11 @@ int main (int argc, char ** argv)
SchurDiagMooeeOperator<DomainWallFermionR,LatticeFermion> HermOpEO(Ddwf); SchurDiagMooeeOperator<DomainWallFermionR,LatticeFermion> HermOpEO(Ddwf);
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2); HermOpEO.MpcDagMpc(chi_e,dchi_e);
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2); HermOpEO.MpcDagMpc(chi_o,dchi_o);
HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2); HermOpEO.MpcDagMpc(phi_e,dphi_e);
HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2); HermOpEO.MpcDagMpc(phi_o,dphi_o);
pDce = innerProduct(phi_e,dchi_e); pDce = innerProduct(phi_e,dchi_e);
pDco = innerProduct(phi_o,dchi_o); pDco = innerProduct(phi_o,dchi_o);

8
tests/core/Test_gpwilson_even_odd.cc
View File

@ -201,11 +201,11 @@ int main (int argc, char ** argv)
RealD t1,t2; RealD t1,t2;
SchurDiagMooeeOperator<GparityWilsonFermionR,FermionField> HermOpEO(Dw); SchurDiagMooeeOperator<GparityWilsonFermionR,FermionField> HermOpEO(Dw);
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2); HermOpEO.MpcDagMpc(chi_e,dchi_e);
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2); HermOpEO.MpcDagMpc(chi_o,dchi_o);
HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2); HermOpEO.MpcDagMpc(phi_e,dphi_e);
HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2); HermOpEO.MpcDagMpc(phi_o,dphi_o);
pDce = innerProduct(phi_e,dchi_e); pDce = innerProduct(phi_e,dchi_e);
pDco = innerProduct(phi_o,dchi_o); pDco = innerProduct(phi_o,dchi_o);

8
tests/core/Test_mobius_eofa_even_odd.cc
View File

@ -220,11 +220,11 @@ int main (int argc, char ** argv)
RealD t1,t2; RealD t1,t2;
SchurDiagMooeeOperator<MobiusEOFAFermionR,LatticeFermion> HermOpEO(Ddwf); SchurDiagMooeeOperator<MobiusEOFAFermionR,LatticeFermion> HermOpEO(Ddwf);
HermOpEO.MpcDagMpc(chi_e, dchi_e, t1, t2); HermOpEO.MpcDagMpc(chi_e, dchi_e);
HermOpEO.MpcDagMpc(chi_o, dchi_o, t1, t2); HermOpEO.MpcDagMpc(chi_o, dchi_o);
HermOpEO.MpcDagMpc(phi_e, dphi_e, t1, t2); HermOpEO.MpcDagMpc(phi_e, dphi_e);
HermOpEO.MpcDagMpc(phi_o, dphi_o, t1, t2); HermOpEO.MpcDagMpc(phi_o, dphi_o);
pDce = innerProduct(phi_e, dchi_e); pDce = innerProduct(phi_e, dchi_e);
pDco = innerProduct(phi_o, dchi_o); pDco = innerProduct(phi_o, dchi_o);

8
tests/core/Test_mobius_even_odd.cc
View File

@ -266,11 +266,11 @@ int main (int argc, char ** argv)
SchurDiagMooeeOperator<MobiusFermionR,LatticeFermion> HermOpEO(Ddwf); SchurDiagMooeeOperator<MobiusFermionR,LatticeFermion> HermOpEO(Ddwf);
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2); HermOpEO.MpcDagMpc(chi_e,dchi_e);
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2); HermOpEO.MpcDagMpc(chi_o,dchi_o);
HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2); HermOpEO.MpcDagMpc(phi_e,dphi_e);
HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2); HermOpEO.MpcDagMpc(phi_o,dphi_o);
pDce = innerProduct(phi_e,dchi_e); pDce = innerProduct(phi_e,dchi_e);
pDco = innerProduct(phi_o,dchi_o); pDco = innerProduct(phi_o,dchi_o);

8
tests/core/Test_staggered.cc
View File

@ -270,11 +270,11 @@ int main (int argc, char ** argv)
pickCheckerboard(Odd ,phi_o,phi); pickCheckerboard(Odd ,phi_o,phi);
SchurDiagMooeeOperator<ImprovedStaggeredFermionR,FermionField> HermOpEO(Ds); SchurDiagMooeeOperator<ImprovedStaggeredFermionR,FermionField> HermOpEO(Ds);
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2); HermOpEO.MpcDagMpc(chi_e,dchi_e);
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2); HermOpEO.MpcDagMpc(chi_o,dchi_o);
HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2); HermOpEO.MpcDagMpc(phi_e,dphi_e);
HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2); HermOpEO.MpcDagMpc(phi_o,dphi_o);
pDce = innerProduct(phi_e,dchi_e); pDce = innerProduct(phi_e,dchi_e);
pDco = innerProduct(phi_o,dchi_o); pDco = innerProduct(phi_o,dchi_o);

8
tests/core/Test_staggered5D.cc
View File

@ -289,11 +289,11 @@ int main (int argc, char ** argv)
pickCheckerboard(Odd ,phi_o,phi); pickCheckerboard(Odd ,phi_o,phi);
SchurDiagMooeeOperator<ImprovedStaggeredFermion5DR,FermionField> HermOpEO(Ds); SchurDiagMooeeOperator<ImprovedStaggeredFermion5DR,FermionField> HermOpEO(Ds);
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2); HermOpEO.MpcDagMpc(chi_e,dchi_e);
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2); HermOpEO.MpcDagMpc(chi_o,dchi_o);
HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2); HermOpEO.MpcDagMpc(phi_e,dphi_e);
HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2); HermOpEO.MpcDagMpc(phi_o,dphi_o);
pDce = innerProduct(phi_e,dchi_e); pDce = innerProduct(phi_e,dchi_e);
pDco = innerProduct(phi_o,dchi_o); pDco = innerProduct(phi_o,dchi_o);

8
tests/core/Test_wilson_even_odd.cc
View File

@ -207,11 +207,11 @@ int main (int argc, char ** argv)
RealD t1,t2; RealD t1,t2;
SchurDiagMooeeOperator<WilsonFermionR,LatticeFermion> HermOpEO(Dw); SchurDiagMooeeOperator<WilsonFermionR,LatticeFermion> HermOpEO(Dw);
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2); HermOpEO.MpcDagMpc(chi_e,dchi_e);
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2); HermOpEO.MpcDagMpc(chi_o,dchi_o);
HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2); HermOpEO.MpcDagMpc(phi_e,dphi_e);
HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2); HermOpEO.MpcDagMpc(phi_o,dphi_o);
pDce = innerProduct(phi_e,dchi_e); pDce = innerProduct(phi_e,dchi_e);
pDco = innerProduct(phi_o,dchi_o); pDco = innerProduct(phi_o,dchi_o);

8
tests/core/Test_wilson_twisted_mass_even_odd.cc
View File

@ -208,11 +208,11 @@ int main (int argc, char ** argv)
RealD t1,t2; RealD t1,t2;
SchurDiagMooeeOperator<WilsonTMFermionR,LatticeFermion> HermOpEO(Dw); SchurDiagMooeeOperator<WilsonTMFermionR,LatticeFermion> HermOpEO(Dw);
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2); HermOpEO.MpcDagMpc(chi_e,dchi_e);
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2); HermOpEO.MpcDagMpc(chi_o,dchi_o);
HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2); HermOpEO.MpcDagMpc(phi_e,dphi_e);
HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2); HermOpEO.MpcDagMpc(phi_o,dphi_o);
pDce = innerProduct(phi_e,dchi_e); pDce = innerProduct(phi_e,dchi_e);
pDco = innerProduct(phi_o,dchi_o); pDco = innerProduct(phi_o,dchi_o);

8
tests/core/Test_zmobius_even_odd.cc
View File

@ -280,11 +280,11 @@ int main (int argc, char ** argv)
SchurDiagMooeeOperator<ZMobiusFermionR,LatticeFermion> HermOpEO(Ddwf); SchurDiagMooeeOperator<ZMobiusFermionR,LatticeFermion> HermOpEO(Ddwf);
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2); HermOpEO.MpcDagMpc(chi_e,dchi_e);
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2); HermOpEO.MpcDagMpc(chi_o,dchi_o);
HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2); HermOpEO.MpcDagMpc(phi_e,dphi_e);
HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2); HermOpEO.MpcDagMpc(phi_o,dphi_o);
pDce = innerProduct(phi_e,dchi_e); pDce = innerProduct(phi_e,dchi_e);
pDco = innerProduct(phi_o,dchi_o); pDco = innerProduct(phi_o,dchi_o);

1
tests/hmc/Test_multishift_sqrt.cc
View File

@ -57,6 +57,7 @@ public:
// Support for coarsening to a multigrid // Support for coarsening to a multigrid
void OpDiag (const Field &in, Field &out) {}; void OpDiag (const Field &in, Field &out) {};
void OpDir (const Field &in, Field &out,int dir,int disp){}; void OpDir (const Field &in, Field &out,int dir,int disp){};
void OpDirAll (const Field &in, std::vector<Field> &out) {}; // Abstract base
void Op (const Field &in, Field &out){ void Op (const Field &in, Field &out){
out = scale * in; out = scale * in;

11
tests/solver/Test_dwf_fpgcr.cc
View File

@ -70,9 +70,6 @@ int main (int argc, char ** argv)
SU3::HotConfiguration(RNG4,Umu); SU3::HotConfiguration(RNG4,Umu);
TrivialPrecon<LatticeFermion> simple;
PrecGeneralisedConjugateResidual<LatticeFermion> PGCR(1.0e-6,10000,simple,4,160);
ConjugateResidual<LatticeFermion> CR(1.0e-6,10000); ConjugateResidual<LatticeFermion> CR(1.0e-6,10000);
@ -86,15 +83,19 @@ int main (int argc, char ** argv)
std::cout<<GridLogMessage<<"* Solving with MdagM VPGCR "<<std::endl; std::cout<<GridLogMessage<<"* Solving with MdagM VPGCR "<<std::endl;
std::cout<<GridLogMessage<<"*********************************************************"<<std::endl; std::cout<<GridLogMessage<<"*********************************************************"<<std::endl;
MdagMLinearOperator<DomainWallFermionR,LatticeFermion> HermOp(Ddwf); MdagMLinearOperator<DomainWallFermionR,LatticeFermion> HermOp(Ddwf);
TrivialPrecon<LatticeFermion> simple;
PrecGeneralisedConjugateResidual<LatticeFermion> PGCR(1.0e-6,10000,HermOp,simple,4,160);
result=Zero(); result=Zero();
PGCR(HermOp,src,result); PGCR(src,result);
std::cout<<GridLogMessage<<"*********************************************************"<<std::endl; std::cout<<GridLogMessage<<"*********************************************************"<<std::endl;
std::cout<<GridLogMessage<<"* Solving with g5-VPGCR "<<std::endl; std::cout<<GridLogMessage<<"* Solving with g5-VPGCR "<<std::endl;
std::cout<<GridLogMessage<<"*********************************************************"<<std::endl; std::cout<<GridLogMessage<<"*********************************************************"<<std::endl;
Gamma5R5HermitianLinearOperator<DomainWallFermionR,LatticeFermion> g5HermOp(Ddwf); Gamma5R5HermitianLinearOperator<DomainWallFermionR,LatticeFermion> g5HermOp(Ddwf);
PrecGeneralisedConjugateResidual<LatticeFermion> PGCR5(1.0e-6,10000,g5HermOp,simple,4,160);
result=Zero(); result=Zero();
PGCR(g5HermOp,src,result); PGCR5(src,result);
std::cout<<GridLogMessage<<"*********************************************************"<<std::endl; std::cout<<GridLogMessage<<"*********************************************************"<<std::endl;
std::cout<<GridLogMessage<<"* Solving with MdagM-CR "<<std::endl; std::cout<<GridLogMessage<<"* Solving with MdagM-CR "<<std::endl;

11
tests/solver/Test_staggered_block_cg_prec.cc
View File

@ -128,9 +128,7 @@ int main (int argc, char ** argv)
std::cout<<GridLogMessage << "usec = "<< (t2-t1)<<std::endl; std::cout<<GridLogMessage << "usec = "<< (t2-t1)<<std::endl;
std::cout<<GridLogMessage << "flops = "<< flops<<std::endl; std::cout<<GridLogMessage << "flops = "<< flops<<std::endl;
std::cout<<GridLogMessage << "mflop/s = "<< flops/(t2-t1)<<std::endl; std::cout<<GridLogMessage << "mflop/s = "<< flops/(t2-t1)<<std::endl;
HermOp4d.Report();
} }
Ds4d.Report();
std::cout << GridLogMessage << "************************************************************************ "<<std::endl; std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
@ -148,9 +146,7 @@ int main (int argc, char ** argv)
std::cout<<GridLogMessage << "usec = "<< (t2-t1)<<std::endl; std::cout<<GridLogMessage << "usec = "<< (t2-t1)<<std::endl;
std::cout<<GridLogMessage << "flops = "<< flops<<std::endl; std::cout<<GridLogMessage << "flops = "<< flops<<std::endl;
std::cout<<GridLogMessage << "mflop/s = "<< flops/(t2-t1)<<std::endl; std::cout<<GridLogMessage << "mflop/s = "<< flops/(t2-t1)<<std::endl;
HermOp.Report();
} }
Ds.Report();
std::cout << GridLogMessage << "************************************************************************ "<<std::endl; std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
std::cout << GridLogMessage << "************************************************************************ "<<std::endl; std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
@ -167,10 +163,8 @@ int main (int argc, char ** argv)
std::cout<<GridLogMessage << "usec = "<< (t2-t1)<<std::endl; std::cout<<GridLogMessage << "usec = "<< (t2-t1)<<std::endl;
std::cout<<GridLogMessage << "flops = "<< flops<<std::endl; std::cout<<GridLogMessage << "flops = "<< flops<<std::endl;
std::cout<<GridLogMessage << "mflop/s = "<< flops/(t2-t1)<<std::endl; std::cout<<GridLogMessage << "mflop/s = "<< flops/(t2-t1)<<std::endl;
HermOp.Report();
} }
Ds.Report();
std::cout << GridLogMessage << "************************************************************************ "<<std::endl; std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
std::cout << GridLogMessage << "************************************************************************ "<<std::endl; std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
@ -187,9 +181,7 @@ int main (int argc, char ** argv)
std::cout<<GridLogMessage << "usec = "<< (t2-t1)<<std::endl; std::cout<<GridLogMessage << "usec = "<< (t2-t1)<<std::endl;
std::cout<<GridLogMessage << "flops = "<< flops<<std::endl; std::cout<<GridLogMessage << "flops = "<< flops<<std::endl;
std::cout<<GridLogMessage << "mflop/s = "<< flops/(t2-t1)<<std::endl; std::cout<<GridLogMessage << "mflop/s = "<< flops/(t2-t1)<<std::endl;
HermOp.Report();
} }
Ds.Report();
std::cout << GridLogMessage << "************************************************************************ "<<std::endl; std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
std::cout << GridLogMessage << "************************************************************************ "<<std::endl; std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
@ -206,9 +198,7 @@ int main (int argc, char ** argv)
std::cout<<GridLogMessage << "usec = "<< (t2-t1)<<std::endl; std::cout<<GridLogMessage << "usec = "<< (t2-t1)<<std::endl;
std::cout<<GridLogMessage << "flops = "<< flops<<std::endl; std::cout<<GridLogMessage << "flops = "<< flops<<std::endl;
std::cout<<GridLogMessage << "mflop/s = "<< flops/(t2-t1)<<std::endl; std::cout<<GridLogMessage << "mflop/s = "<< flops/(t2-t1)<<std::endl;
HermOp.Report();
} }
Ds.Report();
std::cout << GridLogMessage << "************************************************************************ "<<std::endl; std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
std::cout << GridLogMessage << "****************************************************************** "<<std::endl; std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
@ -232,7 +222,6 @@ int main (int argc, char ** argv)
std::cout<<GridLogMessage << "usec = "<< (t2-t1)<<std::endl; std::cout<<GridLogMessage << "usec = "<< (t2-t1)<<std::endl;
std::cout<<GridLogMessage << "flops = "<< flops<<std::endl; std::cout<<GridLogMessage << "flops = "<< flops<<std::endl;
std::cout<<GridLogMessage << "mflop/s = "<< flops/(t2-t1)<<std::endl; std::cout<<GridLogMessage << "mflop/s = "<< flops/(t2-t1)<<std::endl;
// HermOp4d.Report();
} }

2
tests/solver/Test_zMADWF_prec.cc
View File

@ -220,7 +220,7 @@ void run(const TestParams &params){
GridStopWatch CGTimer; GridStopWatch CGTimer;
typename RunParamsOuter::HermOpType<MobiusFermionD> HermOpEO_outer(D_outer); typename RunParamsOuter::template HermOpType<MobiusFermionD> HermOpEO_outer(D_outer);
CGTimer.Start(); CGTimer.Start();
CG_outer(HermOpEO_outer, src_o_outer, result_o_outer); CG_outer(HermOpEO_outer, src_o_outer, result_o_outer);