portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-19 16:27:05 +01:00
Code Releases Activity

Merge branch 'sycl' of https://github.com/paboyle/Grid into sycl

Browse Source
This commit is contained in:
Peter Boyle
2020-06-10 12:58:13 -04:00
parent 8720aecb80 cdf0a04fc5
commit fd97f64612
88 changed files with 2663 additions and 1359 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

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

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

@ -41,8 +41,8 @@ public:
public:
// override multiply
virtual RealD M (const FermionField &in, FermionField &out);
virtual RealD Mdag (const FermionField &in, FermionField &out);
virtual void M (const FermionField &in, FermionField &out);
virtual void Mdag (const FermionField &in, FermionField &out);
// half checkerboard operaions
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);
// override multiply
virtual RealD M (const FermionField& in, FermionField& out);
virtual RealD Mdag (const FermionField& in, FermionField& out);
virtual void M (const FermionField& in, FermionField& out);
virtual void Mdag (const FermionField& in, FermionField& out);
// half checkerboard operations
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;
// override multiply
virtual RealD M (const FermionField &in, FermionField &out)=0;
virtual RealD Mdag (const FermionField &in, FermionField &out)=0;
virtual void M (const FermionField &in, FermionField &out)=0;
virtual void Mdag (const FermionField &in, FermionField &out)=0;
// half checkerboard operaions
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 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 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 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());

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

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

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

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

2
Grid/qcd/action/fermion/MobiusFermion.h
View File

@ -59,7 +59,7 @@ public:
{
RealD eps = 1.0;
std::cout<<GridLogMessage << "MobiusFermion (b="<<b<<",c="<<c<<") with Ls= "<<this->Ls<<" Tanh approx"<<std::endl;
// std::cout<<GridLogMessage << "MobiusFermion (b="<<b<<",c="<<c<<") with Ls= "<<this->Ls<<" Tanh approx"<<std::endl;
Approx::zolotarev_data *zdata = Approx::higham(eps,this->Ls);// eps is ignored for higham
assert(zdata->n==this->Ls);

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

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

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

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

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

@ -109,9 +109,8 @@ public:
ImportGauge(_Umu);
}
virtual RealD M(const FermionField &in, FermionField &out);
virtual RealD Mdag(const FermionField &in, FermionField &out);
virtual void M(const FermionField &in, FermionField &out);
virtual void Mdag(const FermionField &in, FermionField &out);
virtual void Mooee(const FermionField &in, FermionField &out);
virtual void MooeeDag(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
// and also make interface more uniformly consistent
//////////////////////////////////////////////////////////////////
virtual RealD M(const FermionField &in, FermionField &out);
virtual RealD Mdag(const FermionField &in, FermionField &out);
virtual void M(const FermionField &in, FermionField &out);
virtual void Mdag(const FermionField &in, FermionField &out);
/////////////////////////////////////////////////////////
// half checkerboard operations

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

@ -1,4 +1,3 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
@ -99,8 +98,8 @@ public:
GridBase *FermionRedBlackGrid(void) { return _FiveDimRedBlackGrid;}
// full checkerboard operations; leave unimplemented as abstract for now
virtual RealD M (const FermionField &in, FermionField &out){assert(0); return 0.0;};
virtual RealD Mdag (const FermionField &in, FermionField &out){assert(0); return 0.0;};
virtual void M (const FermionField &in, FermionField &out){assert(0);};
virtual void Mdag (const FermionField &in, FermionField &out){assert(0);};
// half checkerboard operations; leave unimplemented as abstract for now
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();
this->Dhop(in, out, DaggerNo);
FermionField tmp(out.Grid());
@ -129,11 +130,12 @@ class WilsonTMFermion5D : public WilsonFermion5D<Impl>
ComplexD b(0.0,this->mu[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
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() == this->FermionGrid()->_fdimensions[0]);
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>
RealD CayleyFermion5D<Impl>::M (const FermionField &psi, FermionField &chi)
void CayleyFermion5D<Impl>::M (const FermionField &psi, FermionField &chi)
{
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);
M5D(psi,chi);
return(norm2(chi));
}
template<class Impl>
RealD CayleyFermion5D<Impl>::Mdag (const FermionField &psi, FermionField &chi)
void CayleyFermion5D<Impl>::Mdag (const FermionField &psi, FermionField &chi)
{
// Under adjoint
//D1+ D1- P- -> D1+^dag P+ D2-^dag
@ -354,7 +353,6 @@ RealD CayleyFermion5D<Impl>::Mdag (const FermionField &psi, FermionField &chi)
M5Ddag(psi,chi);
// ((b D_W + D_w hop terms +1) on s-diag
axpby (chi,1.0,1.0,chi,psi);
return norm2(chi);
}
// 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>
RealD ContinuedFractionFermion5D<Impl>::M (const FermionField &psi, FermionField &chi)
void ContinuedFractionFermion5D<Impl>::M (const FermionField &psi, FermionField &chi)
{
int Ls = this->Ls;
@ -116,15 +116,14 @@ RealD ContinuedFractionFermion5D<Impl>::M (const FermionField &psi, F
}
sign=-sign;
}
return norm2(chi);
}
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
// The rest of matrix is symmetric.
// Can ignore "dag"
return M(psi,chi);
M(psi,chi);
}
template<class Impl>
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>
RealD DomainWallEOFAFermion<Impl>::M(const FermionField& psi, FermionField& chi)
void DomainWallEOFAFermion<Impl>::M(const FermionField& psi, FermionField& chi)
{
FermionField Din(psi.Grid());
@ -97,11 +97,10 @@ RealD DomainWallEOFAFermion<Impl>::M(const FermionField& psi, FermionField& chi)
this->DW(Din, chi, DaggerNo);
axpby(chi, 1.0, 1.0, chi, psi);
this->M5D(psi, chi);
return(norm2(chi));
}
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());
@ -109,7 +108,6 @@ RealD DomainWallEOFAFermion<Impl>::Mdag(const FermionField& psi, FermionField& c
this->MeooeDag5D(Din, chi);
this->M5Ddag(psi, chi);
axpby(chi, 1.0, 1.0, chi, psi);
return(norm2(chi));
}
/********************************************************************

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

@ -470,21 +470,24 @@ void ImprovedStaggeredFermion5D<Impl>::MdirAll(const FermionField &in, std::vect
assert(0);
}
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();
Dhop(in, out, DaggerNo);
return axpy_norm(out, mass, in, out);
axpy(out, mass, in, out);
}
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();
Dhop(in, out, DaggerYes);
return axpy_norm(out, mass, in, out);
axpy(out, mass, in, out);
}
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) {
DhopEO(in, out, DaggerNo);
} else {
@ -492,7 +495,8 @@ void ImprovedStaggeredFermion5D<Impl>::Meooe(const FermionField &in, FermionFiel
}
}
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) {
DhopEO(in, out, DaggerYes);
} else {
@ -501,27 +505,30 @@ void ImprovedStaggeredFermion5D<Impl>::MeooeDag(const FermionField &in, FermionF
}
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();
typename FermionField::scalar_type scal(mass);
out = scal * in;
}
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();
Mooee(in, out);
}
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 = (1.0 / (mass)) * in;
}
template <class Impl>
void ImprovedStaggeredFermion5D<Impl>::MooeeInvDag(const FermionField &in,
FermionField &out) {
void ImprovedStaggeredFermion5D<Impl>::MooeeInvDag(const FermionField &in,FermionField &out)
{
out.Checkerboard() = in.Checkerboard();
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>
RealD ImprovedStaggeredFermion<Impl>::M(const FermionField &in, FermionField &out) {
void ImprovedStaggeredFermion<Impl>::M(const FermionField &in, FermionField &out)
{
out.Checkerboard() = in.Checkerboard();
Dhop(in, out, DaggerNo);
return axpy_norm(out, mass, in, out);
axpy(out, mass, in, out);
}
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();
Dhop(in, out, DaggerYes);
return axpy_norm(out, mass, in, out);
axpy(out, mass, in, out);
}
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) {
DhopEO(in, out, DaggerNo);
} else {
@ -193,7 +196,8 @@ void ImprovedStaggeredFermion<Impl>::Meooe(const FermionField &in, FermionField
}
}
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) {
DhopEO(in, out, DaggerYes);
} else {
@ -202,27 +206,30 @@ void ImprovedStaggeredFermion<Impl>::MeooeDag(const FermionField &in, FermionFie
}
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();
typename FermionField::scalar_type scal(mass);
out = scal * in;
}
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();
Mooee(in, out);
}
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 = (1.0 / (mass)) * in;
}
template <class Impl>
void ImprovedStaggeredFermion<Impl>::MooeeInvDag(const FermionField &in,
FermionField &out) {
void ImprovedStaggeredFermion<Impl>::MooeeInvDag(const FermionField &in,FermionField &out)
{
out.Checkerboard() = in.Checkerboard();
MooeeInv(in, out);
}
@ -234,7 +241,8 @@ void ImprovedStaggeredFermion<Impl>::MooeeInvDag(const FermionField &in,
template <class Impl>
void ImprovedStaggeredFermion<Impl>::DerivInternal(StencilImpl &st, DoubledGaugeField &U, DoubledGaugeField &UUU,
GaugeField & mat,
const FermionField &A, const FermionField &B, int dag) {
const FermionField &A, const FermionField &B, int dag)
{
assert((dag == DaggerNo) || (dag == DaggerYes));
Compressor compressor;
@ -284,8 +292,8 @@ void ImprovedStaggeredFermion<Impl>::DerivInternal(StencilImpl &st, DoubledGauge
}
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(), V.Grid());
conformable(U.Grid(), mat.Grid());
@ -296,8 +304,8 @@ void ImprovedStaggeredFermion<Impl>::DhopDeriv(GaugeField &mat, const FermionFie
}
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(), V.Grid());
conformable(U.Grid(), mat.Grid());
@ -310,8 +318,8 @@ void ImprovedStaggeredFermion<Impl>::DhopDerivOE(GaugeField &mat, const FermionF
}
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(), V.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>
RealD MobiusEOFAFermion<Impl>::M(const FermionField& psi, FermionField& chi)
void MobiusEOFAFermion<Impl>::M(const FermionField& psi, FermionField& chi)
{
FermionField Din(psi.Grid());
@ -174,11 +174,10 @@ RealD MobiusEOFAFermion<Impl>::M(const FermionField& psi, FermionField& chi)
this->DW(Din, chi, DaggerNo);
axpby(chi, 1.0, 1.0, chi, psi);
this->M5D(psi, chi);
return(norm2(chi));
}
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());
@ -186,7 +185,6 @@ RealD MobiusEOFAFermion<Impl>::Mdag(const FermionField& psi, FermionField& chi)
this->MeooeDag5D(Din, chi);
this->M5Ddag(psi, chi);
axpby(chi, 1.0, 1.0, chi, psi);
return(norm2(chi));
}
/********************************************************************

8
Grid/qcd/action/fermion/implementation/NaiveStaggeredFermionImplementation.h
View File

@ -128,17 +128,17 @@ void NaiveStaggeredFermion<Impl>::ImportGauge(const GaugeField &_U)
/////////////////////////////
template <class Impl>
RealD NaiveStaggeredFermion<Impl>::M(const FermionField &in, FermionField &out) {
void NaiveStaggeredFermion<Impl>::M(const FermionField &in, FermionField &out) {
out.Checkerboard() = in.Checkerboard();
Dhop(in, out, DaggerNo);
return axpy_norm(out, mass, in, out);
axpy(out, mass, in, out);
}
template <class Impl>
RealD NaiveStaggeredFermion<Impl>::Mdag(const FermionField &in, FermionField &out) {
void NaiveStaggeredFermion<Impl>::Mdag(const FermionField &in, FermionField &out) {
out.Checkerboard() = in.Checkerboard();
Dhop(in, out, DaggerYes);
return axpy_norm(out, mass, in, out);
axpy(out, mass, in, out);
}
template <class Impl>

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>
RealD PartialFractionFermion5D<Impl>::M (const FermionField &in, FermionField &out)
void PartialFractionFermion5D<Impl>::M (const FermionField &in, FermionField &out)
{
M_internal(in,out,DaggerNo);
return norm2(out);
}
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);
return norm2(out);
}
template<class Impl>

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

@ -35,7 +35,7 @@ NAMESPACE_BEGIN(Grid);
// *NOT* EO
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());
@ -47,11 +47,10 @@ RealD WilsonCloverFermion<Impl>::M(const FermionField &in, FermionField &out)
Mooee(in, temp);
out += temp;
return norm2(out);
}
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());
@ -63,7 +62,6 @@ RealD WilsonCloverFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
MooeeDag(in, temp);
out += temp;
return norm2(out);
}
template <class Impl>
@ -132,14 +130,14 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
pickCheckerboard(Even, CloverTermEven, CloverTerm);
pickCheckerboard(Odd, CloverTermOdd, CloverTerm);
pickCheckerboard(Even, CloverTermDagEven, adj(CloverTerm));
pickCheckerboard(Odd, CloverTermDagOdd, adj(CloverTerm));
pickCheckerboard(Even, CloverTermDagEven, closure(adj(CloverTerm)));
pickCheckerboard(Odd, CloverTermDagOdd, closure(adj(CloverTerm)));
pickCheckerboard(Even, CloverTermInvEven, CloverTermInv);
pickCheckerboard(Odd, CloverTermInvOdd, CloverTermInv);
pickCheckerboard(Even, CloverTermInvDagEven, adj(CloverTermInv));
pickCheckerboard(Odd, CloverTermInvDagOdd, adj(CloverTermInv));
pickCheckerboard(Even, CloverTermInvDagEven, closure(adj(CloverTermInv)));
pickCheckerboard(Odd, CloverTermInvDagOdd, closure(adj(CloverTermInv)));
}
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>
RealD WilsonFermion<Impl>::M(const FermionField &in, FermionField &out) {
void WilsonFermion<Impl>::M(const FermionField &in, FermionField &out)
{
out.Checkerboard() = in.Checkerboard();
Dhop(in, out, DaggerNo);
return axpy_norm(out, diag_mass, in, out);
axpy(out, diag_mass, in, out);
}
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();
Dhop(in, out, DaggerYes);
return axpy_norm(out, diag_mass, in, out);
axpy(out, diag_mass, in, out);
}
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) {
DhopEO(in, out, DaggerNo);
} else {
@ -125,7 +128,8 @@ void WilsonFermion<Impl>::Meooe(const FermionField &in, FermionField &out) {
}
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) {
DhopEO(in, out, DaggerYes);
} else {
@ -134,26 +138,30 @@ void WilsonFermion<Impl>::MeooeDag(const FermionField &in, FermionField &out) {
}
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();
typename FermionField::scalar_type scal(diag_mass);
out = scal * in;
}
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();
Mooee(in, out);
}
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 = (1.0/(diag_mass))*in;
}
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();
MooeeInv(in,out);
}
@ -249,7 +257,8 @@ void WilsonFermion<Impl>::DerivInternal(StencilImpl &st, DoubledGaugeField &U,
}
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(), V.Grid());
conformable(U.Grid(), mat.Grid());
@ -260,7 +269,8 @@ void WilsonFermion<Impl>::DhopDeriv(GaugeField &mat, const FermionField &U, cons
}
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(), V.Grid());
//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>
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(), V.Grid());
//conformable(U.Grid(), mat.Grid());
@ -287,7 +298,8 @@ void WilsonFermion<Impl>::DhopDerivEO(GaugeField &mat, const FermionField &U, co
}
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(), out.Grid());
@ -297,7 +309,8 @@ void WilsonFermion<Impl>::Dhop(const FermionField &in, FermionField &out, int da
}
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(), out.Grid()); // drops the cb check
@ -308,7 +321,8 @@ void WilsonFermion<Impl>::DhopOE(const FermionField &in, FermionField &out, int
}
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(), out.Grid()); // drops the cb check
@ -386,7 +400,8 @@ template <class Impl>
void WilsonFermion<Impl>::DhopInternalOverlappedComms(StencilImpl &st, LebesgueOrder &lo,
DoubledGaugeField &U,
const FermionField &in,
FermionField &out, int dag) {
FermionField &out, int dag)
{
assert((dag == DaggerNo) || (dag == DaggerYes));
Compressor compressor(dag);
@ -436,7 +451,8 @@ template <class Impl>
void WilsonFermion<Impl>::DhopInternalSerial(StencilImpl &st, LebesgueOrder &lo,
DoubledGaugeField &U,
const FermionField &in,
FermionField &out, int dag) {
FermionField &out, int dag)
{
assert((dag == DaggerNo) || (dag == DaggerYes));
Compressor compressor(dag);
st.HaloExchange(in, compressor);

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

@ -447,19 +447,19 @@ void WilsonKernels<Impl>::DhopKernel(int Opt,StencilImpl &st, DoubledGaugeField
if (Opt == WilsonKernelsStatic::OptGeneric ) { KERNEL_CALL(GenericDhopSite); return;}
#ifndef GRID_CUDA
if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL(HandDhopSite); return;}
if (Opt == WilsonKernelsStatic::OptInlineAsm ) { ASM_CALL(AsmDhopSite); printf("."); return;}
if (Opt == WilsonKernelsStatic::OptInlineAsm ) { ASM_CALL(AsmDhopSite); return;}
#endif
} else if( interior ) {
if (Opt == WilsonKernelsStatic::OptGeneric ) { KERNEL_CALLNB(GenericDhopSiteInt); return;}
#ifndef GRID_CUDA
if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALLNB(HandDhopSiteInt); return;}
if (Opt == WilsonKernelsStatic::OptInlineAsm ) { ASM_CALL(AsmDhopSiteInt); printf("-"); return;}
if (Opt == WilsonKernelsStatic::OptInlineAsm ) { ASM_CALL(AsmDhopSiteInt); return;}
#endif
} else if( exterior ) {
if (Opt == WilsonKernelsStatic::OptGeneric ) { KERNEL_CALL(GenericDhopSiteExt); return;}
#ifndef GRID_CUDA
if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL(HandDhopSiteExt); return;}
if (Opt == WilsonKernelsStatic::OptInlineAsm ) { ASM_CALL(AsmDhopSiteExt); printf("+"); return;}
if (Opt == WilsonKernelsStatic::OptInlineAsm ) { ASM_CALL(AsmDhopSiteExt); return;}
#endif
}
assert(0 && " Kernel optimisation case not covered ");

2
Grid/qcd/action/gauge/GaugeImplementations.h
View File

@ -59,7 +59,7 @@ public:
}
static inline GaugeLinkField
CovShiftIdentityBackward(const GaugeLinkField &Link, int mu) {
return Cshift(adj(Link), mu, -1);
return Cshift(closure(adj(Link)), mu, -1);
}
static inline GaugeLinkField
CovShiftIdentityForward(const GaugeLinkField &Link, int mu) {

4
Grid/qcd/hmc/HMC_aggregate.h
View File

@ -39,6 +39,10 @@ directory
#include <Grid/parallelIO/IldgIOtypes.h>
#include <Grid/parallelIO/IldgIO.h>
#include <Grid/parallelIO/NerscIO.h>
#include <Grid/parallelIO/OpenQcdIO.h>
#if !defined(GRID_COMMS_NONE)
#include <Grid/parallelIO/OpenQcdIOChromaReference.h>
#endif
NAMESPACE_CHECK(Ildg);
#include <Grid/qcd/hmc/checkpointers/CheckPointers.h>

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

@ -80,6 +80,7 @@ static Registrar<OneFlavourRatioEOFModule<FermionImplementationPolicy>,
static Registrar< ConjugateGradientModule<WilsonFermionR::FermionField>,
HMC_SolverModuleFactory<solver_string, WilsonFermionR::FermionField, Serialiser> > __CGWFmodXMLInit("ConjugateGradient");
static Registrar< BiCGSTABModule<WilsonFermionR::FermionField>,
HMC_SolverModuleFactory<solver_string, WilsonFermionR::FermionField, Serialiser> > __BiCGWFmodXMLInit("BiCGSTAB");
static Registrar< ConjugateResidualModule<WilsonFermionR::FermionField>,

156
Grid/qcd/utils/BaryonUtils.h
View File

@ -46,7 +46,7 @@ public:
typedef typename SpinMatrixField::vector_object sobj;
static const int epsilon[6][3] ;
static const Complex epsilon_sgn[6];
static const Real epsilon_sgn[6];
private:
template <class mobj, class robj>
@ -151,14 +151,18 @@ public:
template <class FImpl>
const int BaryonUtils<FImpl>::epsilon[6][3] = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}};
template <class FImpl>
/*template <class FImpl>
const Complex BaryonUtils<FImpl>::epsilon_sgn[6] = {Complex(1),
Complex(1),
Complex(1),
Complex(-1),
Complex(-1),
Complex(-1)};
*/
template <class FImpl>
const Real BaryonUtils<FImpl>::epsilon_sgn[6] = {1.,1.,1.,-1.,-1.,-1.};
//This is the old version
template <class FImpl>
template <class mobj, class robj>
void BaryonUtils<FImpl>::baryon_site(const mobj &D1,
@ -173,12 +177,14 @@ void BaryonUtils<FImpl>::baryon_site(const mobj &D1,
robj &result)
{
Gamma g4(Gamma::Algebra::GammaT); //needed for parity P_\pm = 0.5*(1 \pm \gamma_4)
Gamma g4(Gamma::Algebra::GammaT); //needed for parity P_\pm = 0.5*(1 \pm \gamma_4)
auto gD1a = GammaA_left * GammaA_right * D1;
auto gD1b = GammaA_left * g4 * GammaA_right * D1;
auto pD1 = 0.5* (gD1a + (double)parity * gD1b);
auto pD1 = 0.5* (gD1a + (Real)parity * gD1b);
auto gD3 = GammaB_right * D3;
auto D2g = D2 * GammaB_left;
auto pD1g = pD1 * GammaB_left;
auto gD3g = gD3 * GammaB_left;
for (int ie_left=0; ie_left < 6 ; ie_left++){
int a_left = epsilon[ie_left][0]; //a
@ -188,59 +194,78 @@ void BaryonUtils<FImpl>::baryon_site(const mobj &D1,
int a_right = epsilon[ie_right][0]; //a'
int b_right = epsilon[ie_right][1]; //b'
int c_right = epsilon[ie_right][2]; //c'
Real ee = epsilon_sgn[ie_left] * epsilon_sgn[ie_right];
//This is the \delta_{456}^{123} part
if (wick_contraction[0]){
auto D2g = D2 * GammaB_left;
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
for (int beta_left=0; beta_left<Ns; beta_left++){
for (int gamma_left=0; gamma_left<Ns; gamma_left++){
result()()() += epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1()(gamma_left,gamma_left)(c_right,c_left)*D2g()(alpha_right,beta_left)(a_right,a_left)*gD3()(alpha_right,beta_left)(b_right,b_left);
}}}
auto eepD1 = ee * pD1()(gamma_left,gamma_left)(c_right,c_left);
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
for (int beta_left=0; beta_left<Ns; beta_left++){
auto D2g_ab = D2g()(alpha_right,beta_left)(a_right,a_left);
auto gD3_ab = gD3()(alpha_right,beta_left)(b_right,b_left);
result()()() += eepD1*D2g_ab*gD3_ab;
}}
}
}
//This is the \delta_{456}^{231} part
if (wick_contraction[1]){
auto pD1g = pD1 * GammaB_left;
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
for (int beta_left=0; beta_left<Ns; beta_left++){
for (int gamma_left=0; gamma_left<Ns; gamma_left++){
result()()() += epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1g()(gamma_left,beta_left)(c_right,a_left)*D2()(alpha_right,beta_left)(a_right,b_left)*gD3()(alpha_right,gamma_left)(b_right,c_left);
}}}
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
auto gD3_ag = gD3()(alpha_right,gamma_left)(b_right,c_left);
for (int beta_left=0; beta_left<Ns; beta_left++){
auto eepD1g_gb = ee * pD1g()(gamma_left,beta_left)(c_right,a_left);
auto D2_ab = D2()(alpha_right,beta_left)(a_right,b_left);
result()()() += eepD1g_gb*D2_ab*gD3_ag;
}
}}
}
//This is the \delta_{456}^{312} part
if (wick_contraction[2]){
auto gD3g = gD3 * GammaB_left;
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
for (int beta_left=0; beta_left<Ns; beta_left++){
for (int gamma_left=0; gamma_left<Ns; gamma_left++){
result()()() += epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1()(gamma_left,beta_left)(c_right,b_left)*D2()(alpha_right,gamma_left)(a_right,c_left)*gD3g()(alpha_right,beta_left)(b_right,a_left);
}}}
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
auto D2_ag = D2()(alpha_right,gamma_left)(a_right,c_left);
for (int beta_left=0; beta_left<Ns; beta_left++){
auto eepD1_gb = ee * pD1()(gamma_left,beta_left)(c_right,b_left);
auto gD3g_ab = gD3g()(alpha_right,beta_left)(b_right,a_left);
result()()() += eepD1_gb*D2_ag*gD3g_ab;
}
}}
}
//This is the \delta_{456}^{132} part
if (wick_contraction[3]){
auto gD3g = gD3 * GammaB_left;
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
for (int beta_left=0; beta_left<Ns; beta_left++){
for (int gamma_left=0; gamma_left<Ns; gamma_left++){
result()()() -= epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1()(gamma_left,gamma_left)(c_right,c_left)*D2()(alpha_right,beta_left)(a_right,b_left)*gD3g()(alpha_right,beta_left)(b_right,a_left);
}}}
auto eepD1 = ee * pD1()(gamma_left,gamma_left)(c_right,c_left);
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
for (int beta_left=0; beta_left<Ns; beta_left++){
auto D2_ab = D2()(alpha_right,beta_left)(a_right,b_left);
auto gD3g_ab = gD3g()(alpha_right,beta_left)(b_right,a_left);
result()()() -= eepD1*D2_ab*gD3g_ab;
}}
}
}
//This is the \delta_{456}^{321} part
if (wick_contraction[4]){
auto D2g = D2 * GammaB_left;
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
for (int beta_left=0; beta_left<Ns; beta_left++){
for (int gamma_left=0; gamma_left<Ns; gamma_left++){
result()()() -= epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1()(gamma_left,beta_left)(c_right,b_left)*D2g()(alpha_right,beta_left)(a_right,a_left)*gD3()(alpha_right,gamma_left)(b_right,c_left);
}}}
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
auto gD3_ag = gD3()(alpha_right,gamma_left)(b_right,c_left);
for (int beta_left=0; beta_left<Ns; beta_left++){
auto eepD1_gb = ee * pD1()(gamma_left,beta_left)(c_right,b_left);
auto D2g_ab = D2g()(alpha_right,beta_left)(a_right,a_left);
result()()() -= eepD1_gb*D2g_ab*gD3_ag;
}
}}
}
//This is the \delta_{456}^{213} part
if (wick_contraction[5]){
auto pD1g = pD1 * GammaB_left;
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
for (int beta_left=0; beta_left<Ns; beta_left++){
for (int gamma_left=0; gamma_left<Ns; gamma_left++){
result()()() -= epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1g()(gamma_left,beta_left)(c_right,a_left)*D2()(alpha_right,gamma_left)(a_right,c_left)*gD3()(alpha_right,beta_left)(b_right,b_left);
}}}
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
auto D2_ag = D2()(alpha_right,gamma_left)(a_right,c_left);
for (int beta_left=0; beta_left<Ns; beta_left++){
auto eepD1g_gb = ee * pD1g()(gamma_left,beta_left)(c_right,a_left);
auto gD3_ab = gD3()(alpha_right,beta_left)(b_right,b_left);
result()()() -= eepD1g_gb*D2_ag*gD3_ab;
}
}}
}
}
}
@ -259,11 +284,15 @@ void BaryonUtils<FImpl>::ContractBaryons(const PropagatorField &q1_left,
const int parity,
ComplexField &baryon_corr)
{
assert(Ns==4 && "Baryon code only implemented for N_spin = 4");
assert(Nc==3 && "Baryon code only implemented for N_colour = 3");
std::cout << "Contraction <" << quarks_right[0] << quarks_right[1] << quarks_right[2] << "|" << quarks_left[0] << quarks_left[1] << quarks_left[2] << ">" << std::endl;
std::cout << "GammaA (left) " << (GammaA_left.g) << std::endl;
std::cout << "GammaB (left) " << (GammaB_left.g) << std::endl;
std::cout << "GammaA (right) " << (GammaA_right.g) << std::endl;
std::cout << "GammaB (right) " << (GammaB_right.g) << std::endl;
std::cout << "GammaA (left) " << (GammaA_left.g) << std::endl;
std::cout << "GammaB (left) " << (GammaB_left.g) << std::endl;
std::cout << "GammaA (right) " << (GammaA_right.g) << std::endl;
std::cout << "GammaB (right) " << (GammaB_right.g) << std::endl;
assert(parity==1 || parity == -1 && "Parity must be +1 or -1");
@ -278,18 +307,32 @@ void BaryonUtils<FImpl>::ContractBaryons(const PropagatorField &q1_left,
autoView( v2 , q2_left, CpuRead);
autoView( v3 , q3_left, CpuRead);
// accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
thread_for(ss,grid->oSites(),{
//for(int ss=0; ss < grid->oSites(); ss++){
Real bytes =0.;
bytes += grid->oSites() * (432.*sizeof(vComplex) + 126.*sizeof(int) + 36.*sizeof(Real));
for (int ie=0; ie < 6 ; ie++){
if(ie==0 or ie==3){
bytes += grid->oSites() * (4.*sizeof(int) + 4752.*sizeof(vComplex)) * wick_contraction[ie];
}
else{
bytes += grid->oSites() * (64.*sizeof(int) + 5184.*sizeof(vComplex)) * wick_contraction[ie];
}
}
Real t=0.;
t =-usecond();
accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
auto D1 = v1[ss];
auto D2 = v2[ss];
auto D3 = v3[ss];
vobj result=Zero();
baryon_site(D1,D2,D3,GammaA_left,GammaB_left,GammaA_right,GammaB_right,parity,wick_contraction,result);
vbaryon_corr[ss] = result;
} );//end loop over lattice sites
t += usecond();
std::cout << std::setw(10) << bytes/t*1.0e6/1024/1024/1024 << " GB/s " << std::endl;
}
template <class FImpl>
template <class mobj, class robj>
@ -305,11 +348,15 @@ void BaryonUtils<FImpl>::ContractBaryons_Sliced(const mobj &D1,
const int parity,
robj &result)
{
assert(Ns==4 && "Baryon code only implemented for N_spin = 4");
assert(Nc==3 && "Baryon code only implemented for N_colour = 3");
std::cout << "Contraction <" << quarks_right[0] << quarks_right[1] << quarks_right[2] << "|" << quarks_left[0] << quarks_left[1] << quarks_left[2] << ">" << std::endl;
std::cout << "GammaA (left) " << (GammaA_left.g) << std::endl;
std::cout << "GammaB (left) " << (GammaB_left.g) << std::endl;
std::cout << "GammaA (right) " << (GammaA_right.g) << std::endl;
std::cout << "GammaB (right) " << (GammaB_right.g) << std::endl;
std::cout << "GammaA (left) " << (GammaA_left.g) << std::endl;
std::cout << "GammaB (left) " << (GammaB_left.g) << std::endl;
std::cout << "GammaA (right) " << (GammaA_right.g) << std::endl;
std::cout << "GammaB (right) " << (GammaB_right.g) << std::endl;
assert(parity==1 || parity == -1 && "Parity must be +1 or -1");
@ -317,8 +364,8 @@ void BaryonUtils<FImpl>::ContractBaryons_Sliced(const mobj &D1,
for (int ie=0; ie < 6 ; ie++)
wick_contraction[ie] = (quarks_left[0] == quarks_right[epsilon[ie][0]] && quarks_left[1] == quarks_right[epsilon[ie][1]] && quarks_left[2] == quarks_right[epsilon[ie][2]]) ? 1 : 0;
result=Zero();
baryon_site(D1,D2,D3,GammaA_left,GammaB_left,GammaA_right,GammaB_right,parity,wick_contraction,result);
result=Zero();
baryon_site<decltype(D1),decltype(result)>(D1,D2,D3,GammaA_left,GammaB_left,GammaA_right,GammaB_right,parity,wick_contraction,result);
}
/***********************************************************************
@ -558,6 +605,10 @@ void BaryonUtils<FImpl>::Sigma_to_Nucleon_Eye(const PropagatorField &qq_loop,
const std::string op,
SpinMatrixField &stn_corr)
{
assert(Ns==4 && "Baryon code only implemented for N_spin = 4");
assert(Nc==3 && "Baryon code only implemented for N_colour = 3");
GridBase *grid = qs_ti.Grid();
autoView( vcorr, stn_corr, CpuWrite);
@ -565,8 +616,7 @@ void BaryonUtils<FImpl>::Sigma_to_Nucleon_Eye(const PropagatorField &qq_loop,
autoView( vd_tf , qd_tf, CpuRead);
autoView( vs_ti , qs_ti, CpuRead);
// accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
thread_for(ss,grid->oSites(),{
accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
auto Dq_loop = vq_loop[ss];
auto Dd_tf = vd_tf[ss];
auto Ds_ti = vs_ti[ss];
@ -595,6 +645,10 @@ void BaryonUtils<FImpl>::Sigma_to_Nucleon_NonEye(const PropagatorField &qq_ti,
const std::string op,
SpinMatrixField &stn_corr)
{
assert(Ns==4 && "Baryon code only implemented for N_spin = 4");
assert(Nc==3 && "Baryon code only implemented for N_colour = 3");
GridBase *grid = qs_ti.Grid();
autoView( vcorr , stn_corr, CpuWrite);

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

@ -52,6 +52,26 @@ namespace PeriodicBC {
tmp = adj(Link)*field;
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
}
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);
}
}

4
Grid/qcd/utils/WilsonLoops.h
View File

@ -485,7 +485,7 @@ public:
// Up staple ___ ___
// | |
tmp = Cshift(adj(U[nu]), nu, -1);
tmp = Cshift(closure(adj(U[nu])), nu, -1);
tmp = adj(U2[mu]) * tmp;
tmp = Cshift(tmp, mu, -2);
@ -519,7 +519,7 @@ public:
//
// | |
tmp = Cshift(adj(U2[nu]), nu, -2);
tmp = Cshift(closure(adj(U2[nu])), nu, -2);
tmp = Gimpl::CovShiftBackward(U[mu], mu, tmp);
tmp = U2[nu] * Cshift(tmp, nu, 2);
Stap += Cshift(tmp, mu, 1);