mirror of
https://github.com/paboyle/Grid.git
synced 2025-04-04 19:25:56 +01:00
Push 5d tests
This commit is contained in:
Azusa Yamaguchi
parent
389e0a77bd
commit
77fb25fb29
344
lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc
Normal file
344
lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc
Normal file
@ -0,0 +1,344 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc
|
||||
|
||||
Copyright (C) 2015
|
||||
|
||||
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
|
||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
|
||||
This program is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License along
|
||||
with this program; if not, write to the Free Software Foundation, Inc.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
|
||||
See the full license in the file "LICENSE" in the top level distribution directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#include <Grid.h>
|
||||
#include <PerfCount.h>
|
||||
|
||||
namespace Grid {
|
||||
namespace QCD {
|
||||
|
||||
// S-direction is INNERMOST and takes no part in the parity.
|
||||
const std::vector<int>
|
||||
ImprovedStaggeredFermion5DStatic::directions({1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4});
|
||||
const std::vector<int>
|
||||
ImprovedStaggeredFermion5DStatic::displacements({1, 1, 1, 1, -1, -1, -1, -1, 3, 3, 3, 3, -3, -3, -3, -3});
|
||||
|
||||
// 5d lattice for DWF.
|
||||
template<class Impl>
|
||||
ImprovedStaggeredFermion5D<Impl>::ImprovedStaggeredFermion5D(GaugeField &_Uthin,GaugeField &_Ufat,
|
||||
GridCartesian &FiveDimGrid,
|
||||
GridRedBlackCartesian &FiveDimRedBlackGrid,
|
||||
GridCartesian &FourDimGrid,
|
||||
GridRedBlackCartesian &FourDimRedBlackGrid,
|
||||
RealD _mass,
|
||||
RealD _c1,RealD _c2, RealD _u0,
|
||||
const ImplParams &p) :
|
||||
Kernels(p),
|
||||
_FiveDimGrid (&FiveDimGrid),
|
||||
_FiveDimRedBlackGrid(&FiveDimRedBlackGrid),
|
||||
_FourDimGrid (&FourDimGrid),
|
||||
_FourDimRedBlackGrid(&FourDimRedBlackGrid),
|
||||
Stencil (_FiveDimGrid,npoint,Even,directions,displacements),
|
||||
StencilEven(_FiveDimRedBlackGrid,npoint,Even,directions,displacements), // source is Even
|
||||
StencilOdd (_FiveDimRedBlackGrid,npoint,Odd ,directions,displacements), // source is Odd
|
||||
mass(_mass),
|
||||
c1(_c1),
|
||||
c2(_c2),
|
||||
u0(_u0),
|
||||
Umu(_FourDimGrid),
|
||||
UmuEven(_FourDimRedBlackGrid),
|
||||
UmuOdd (_FourDimRedBlackGrid),
|
||||
UUUmu(_FourDimGrid),
|
||||
UUUmuEven(_FourDimRedBlackGrid),
|
||||
UUUmuOdd(_FourDimRedBlackGrid),
|
||||
Lebesgue(_FourDimGrid),
|
||||
LebesgueEvenOdd(_FourDimRedBlackGrid)
|
||||
{
|
||||
// some assertions
|
||||
assert(FiveDimGrid._ndimension==5);
|
||||
assert(FourDimGrid._ndimension==4);
|
||||
assert(FiveDimRedBlackGrid._ndimension==5);
|
||||
assert(FourDimRedBlackGrid._ndimension==4);
|
||||
assert(FiveDimRedBlackGrid._checker_dim==1);
|
||||
|
||||
// Dimension zero of the five-d is the Ls direction
|
||||
Ls=FiveDimGrid._fdimensions[0];
|
||||
assert(FiveDimRedBlackGrid._fdimensions[0]==Ls);
|
||||
assert(FiveDimRedBlackGrid._processors[0] ==1);
|
||||
assert(FiveDimRedBlackGrid._simd_layout[0]==1);
|
||||
assert(FiveDimGrid._processors[0] ==1);
|
||||
assert(FiveDimGrid._simd_layout[0] ==1);
|
||||
|
||||
// Other dimensions must match the decomposition of the four-D fields
|
||||
for(int d=0;d<4;d++){
|
||||
assert(FourDimRedBlackGrid._fdimensions[d] ==FourDimGrid._fdimensions[d]);
|
||||
assert(FiveDimRedBlackGrid._fdimensions[d+1]==FourDimGrid._fdimensions[d]);
|
||||
|
||||
assert(FourDimRedBlackGrid._processors[d] ==FourDimGrid._processors[d]);
|
||||
assert(FiveDimRedBlackGrid._processors[d+1] ==FourDimGrid._processors[d]);
|
||||
|
||||
assert(FourDimRedBlackGrid._simd_layout[d] ==FourDimGrid._simd_layout[d]);
|
||||
assert(FiveDimRedBlackGrid._simd_layout[d+1]==FourDimGrid._simd_layout[d]);
|
||||
|
||||
assert(FiveDimGrid._fdimensions[d+1] ==FourDimGrid._fdimensions[d]);
|
||||
assert(FiveDimGrid._processors[d+1] ==FourDimGrid._processors[d]);
|
||||
assert(FiveDimGrid._simd_layout[d+1] ==FourDimGrid._simd_layout[d]);
|
||||
}
|
||||
|
||||
// Allocate the required comms buffer
|
||||
ImportGauge(_Uthin,_Ufat);
|
||||
}
|
||||
|
||||
template <class Impl>
|
||||
void ImprovedStaggeredFermion5D<Impl>::ImportGauge(const GaugeField &_Uthin)
|
||||
{
|
||||
ImportGauge(_Uthin,_Uthin);
|
||||
};
|
||||
template<class Impl>
|
||||
void ImprovedStaggeredFermion5D<Impl>::ImportGauge(const GaugeField &_Uthin,const GaugeField &_Ufat)
|
||||
{
|
||||
GaugeLinkField U(GaugeGrid());
|
||||
|
||||
////////////////////////////////////////////////////////
|
||||
// Double Store should take two fields for Naik and one hop separately.
|
||||
////////////////////////////////////////////////////////
|
||||
Impl::DoubleStore(GaugeGrid(), UUUmu, Umu, _Uthin, _Ufat );
|
||||
|
||||
////////////////////////////////////////////////////////
|
||||
// Apply scale factors to get the right fermion Kinetic term
|
||||
// Could pass coeffs into the double store to save work.
|
||||
// 0.5 ( U p(x+mu) - Udag(x-mu) p(x-mu) )
|
||||
////////////////////////////////////////////////////////
|
||||
for (int mu = 0; mu < Nd; mu++) {
|
||||
|
||||
U = PeekIndex<LorentzIndex>(Umu, mu);
|
||||
PokeIndex<LorentzIndex>(Umu, U*( 0.5*c1/u0), mu );
|
||||
|
||||
U = PeekIndex<LorentzIndex>(Umu, mu+4);
|
||||
PokeIndex<LorentzIndex>(Umu, U*(-0.5*c1/u0), mu+4);
|
||||
|
||||
U = PeekIndex<LorentzIndex>(UUUmu, mu);
|
||||
PokeIndex<LorentzIndex>(UUUmu, U*( 0.5*c2/u0/u0/u0), mu );
|
||||
|
||||
U = PeekIndex<LorentzIndex>(UUUmu, mu+4);
|
||||
PokeIndex<LorentzIndex>(UUUmu, U*(-0.5*c2/u0/u0/u0), mu+4);
|
||||
}
|
||||
|
||||
pickCheckerboard(Even, UmuEven, Umu);
|
||||
pickCheckerboard(Odd, UmuOdd , Umu);
|
||||
pickCheckerboard(Even, UUUmuEven, UUUmu);
|
||||
pickCheckerboard(Odd, UUUmuOdd, UUUmu);
|
||||
}
|
||||
template<class Impl>
|
||||
void ImprovedStaggeredFermion5D<Impl>::DhopDir(const FermionField &in, FermionField &out,int dir5,int disp)
|
||||
{
|
||||
int dir = dir5-1; // Maps to the ordering above in "directions" that is passed to stencil
|
||||
// we drop off the innermost fifth dimension
|
||||
|
||||
Compressor compressor;
|
||||
Stencil.HaloExchange(in,compressor);
|
||||
|
||||
PARALLEL_FOR_LOOP
|
||||
for(int ss=0;ss<Umu._grid->oSites();ss++){
|
||||
for(int s=0;s<Ls;s++){
|
||||
int sU=ss;
|
||||
int sF = s+Ls*sU;
|
||||
Kernels::DhopDir(Stencil, Umu, UUUmu, Stencil.CommBuf(), sF, sU, in, out, dir, disp);
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
template<class Impl>
|
||||
void ImprovedStaggeredFermion5D<Impl>::DerivInternal(StencilImpl & st,
|
||||
DoubledGaugeField & U,
|
||||
DoubledGaugeField & UUU,
|
||||
GaugeField &mat,
|
||||
const FermionField &A,
|
||||
const FermionField &B,
|
||||
int dag)
|
||||
{
|
||||
// No force terms in multi-rhs solver staggered
|
||||
assert(0);
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void ImprovedStaggeredFermion5D<Impl>::DhopDeriv(GaugeField &mat,
|
||||
const FermionField &A,
|
||||
const FermionField &B,
|
||||
int dag)
|
||||
{
|
||||
assert(0);
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void ImprovedStaggeredFermion5D<Impl>::DhopDerivEO(GaugeField &mat,
|
||||
const FermionField &A,
|
||||
const FermionField &B,
|
||||
int dag)
|
||||
{
|
||||
assert(0);
|
||||
}
|
||||
|
||||
|
||||
template<class Impl>
|
||||
void ImprovedStaggeredFermion5D<Impl>::DhopDerivOE(GaugeField &mat,
|
||||
const FermionField &A,
|
||||
const FermionField &B,
|
||||
int dag)
|
||||
{
|
||||
assert(0);
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void ImprovedStaggeredFermion5D<Impl>::DhopInternal(StencilImpl & st, LebesgueOrder &lo,
|
||||
DoubledGaugeField & U,DoubledGaugeField & UUU,
|
||||
const FermionField &in, FermionField &out,int dag)
|
||||
{
|
||||
Compressor compressor;
|
||||
|
||||
int LLs = in._grid->_rdimensions[0];
|
||||
|
||||
st.HaloExchange(in,compressor);
|
||||
|
||||
// Dhop takes the 4d grid from U, and makes a 5d index for fermion
|
||||
if (dag == DaggerYes) {
|
||||
PARALLEL_FOR_LOOP
|
||||
for (int ss = 0; ss < U._grid->oSites(); ss++) {
|
||||
for(int s=0;s<LLs;s++){
|
||||
int sU=ss;
|
||||
int sF = s+LLs*sU;
|
||||
Kernels::DhopSiteDag(st, lo, U, UUU, st.CommBuf(), sF, sU, in, out);
|
||||
}}
|
||||
} else {
|
||||
PARALLEL_FOR_LOOP
|
||||
for (int ss = 0; ss < U._grid->oSites(); ss++) {
|
||||
for(int s=0;s<LLs;s++){
|
||||
int sU=ss;
|
||||
int sF = s+LLs*sU;
|
||||
Kernels::DhopSite(st,lo,U,UUU,st.CommBuf(),sF,sU,in,out);
|
||||
}}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class Impl>
|
||||
void ImprovedStaggeredFermion5D<Impl>::DhopOE(const FermionField &in, FermionField &out,int dag)
|
||||
{
|
||||
conformable(in._grid,FermionRedBlackGrid()); // verifies half grid
|
||||
conformable(in._grid,out._grid); // drops the cb check
|
||||
|
||||
assert(in.checkerboard==Even);
|
||||
out.checkerboard = Odd;
|
||||
|
||||
DhopInternal(StencilEven,LebesgueEvenOdd,UmuOdd,UUUmuOdd,in,out,dag);
|
||||
}
|
||||
template<class Impl>
|
||||
void ImprovedStaggeredFermion5D<Impl>::DhopEO(const FermionField &in, FermionField &out,int dag)
|
||||
{
|
||||
conformable(in._grid,FermionRedBlackGrid()); // verifies half grid
|
||||
conformable(in._grid,out._grid); // drops the cb check
|
||||
|
||||
assert(in.checkerboard==Odd);
|
||||
out.checkerboard = Even;
|
||||
|
||||
DhopInternal(StencilOdd,LebesgueEvenOdd,UmuEven,UUUmuEven,in,out,dag);
|
||||
}
|
||||
template<class Impl>
|
||||
void ImprovedStaggeredFermion5D<Impl>::Dhop(const FermionField &in, FermionField &out,int dag)
|
||||
{
|
||||
conformable(in._grid,FermionGrid()); // verifies full grid
|
||||
conformable(in._grid,out._grid);
|
||||
|
||||
out.checkerboard = in.checkerboard;
|
||||
|
||||
DhopInternal(Stencil,Lebesgue,Umu,UUUmu,in,out,dag);
|
||||
}
|
||||
|
||||
|
||||
/////////////////////////////////////////////////////////////////////////
|
||||
// Implement the general interface. Here we use SAME mass on all slices
|
||||
/////////////////////////////////////////////////////////////////////////
|
||||
template <class Impl>
|
||||
void ImprovedStaggeredFermion5D<Impl>::Mdir(const FermionField &in, FermionField &out, int dir, int disp) {
|
||||
DhopDir(in, out, dir, disp);
|
||||
}
|
||||
template <class Impl>
|
||||
RealD ImprovedStaggeredFermion5D<Impl>::M(const FermionField &in, FermionField &out) {
|
||||
out.checkerboard = in.checkerboard;
|
||||
Dhop(in, out, DaggerNo);
|
||||
return axpy_norm(out, mass, in, out);
|
||||
}
|
||||
|
||||
template <class Impl>
|
||||
RealD ImprovedStaggeredFermion5D<Impl>::Mdag(const FermionField &in, FermionField &out) {
|
||||
out.checkerboard = in.checkerboard;
|
||||
Dhop(in, out, DaggerYes);
|
||||
return axpy_norm(out, mass, in, out);
|
||||
}
|
||||
|
||||
template <class Impl>
|
||||
void ImprovedStaggeredFermion5D<Impl>::Meooe(const FermionField &in, FermionField &out) {
|
||||
if (in.checkerboard == Odd) {
|
||||
DhopEO(in, out, DaggerNo);
|
||||
} else {
|
||||
DhopOE(in, out, DaggerNo);
|
||||
}
|
||||
}
|
||||
template <class Impl>
|
||||
void ImprovedStaggeredFermion5D<Impl>::MeooeDag(const FermionField &in, FermionField &out) {
|
||||
if (in.checkerboard == Odd) {
|
||||
DhopEO(in, out, DaggerYes);
|
||||
} else {
|
||||
DhopOE(in, out, DaggerYes);
|
||||
}
|
||||
}
|
||||
|
||||
template <class Impl>
|
||||
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) {
|
||||
out.checkerboard = in.checkerboard;
|
||||
Mooee(in, out);
|
||||
}
|
||||
|
||||
template <class Impl>
|
||||
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) {
|
||||
out.checkerboard = in.checkerboard;
|
||||
MooeeInv(in, out);
|
||||
}
|
||||
|
||||
|
||||
|
||||
FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion5D);
|
||||
|
||||
}}
|
||||
|
||||
|
||||
|
||||
164
lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h
Normal file
164
lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h
Normal file
@ -0,0 +1,164 @@
|
||||
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h
|
||||
|
||||
Copyright (C) 2015
|
||||
|
||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
Author: AzusaYamaguchi <ayamaguc@staffmail.ed.ac.uk>
|
||||
|
||||
This program is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License along
|
||||
with this program; if not, write to the Free Software Foundation, Inc.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
|
||||
See the full license in the file "LICENSE" in the top level distribution directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#ifndef GRID_QCD_IMPROVED_STAGGERED_FERMION_5D_H
|
||||
#define GRID_QCD_IMPROVED_STAGGERED_FERMION_5D_H
|
||||
|
||||
namespace Grid {
|
||||
namespace QCD {
|
||||
|
||||
////////////////////////////////////////////////////////////////////////////////
|
||||
// This is the 4d red black case appropriate to support
|
||||
////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
class ImprovedStaggeredFermion5DStatic {
|
||||
public:
|
||||
// S-direction is INNERMOST and takes no part in the parity.
|
||||
static const std::vector<int> directions;
|
||||
static const std::vector<int> displacements;
|
||||
const int npoint = 16;
|
||||
};
|
||||
|
||||
template<class Impl>
|
||||
class ImprovedStaggeredFermion5D : public StaggeredKernels<Impl>, public ImprovedStaggeredFermion5DStatic
|
||||
{
|
||||
public:
|
||||
INHERIT_IMPL_TYPES(Impl);
|
||||
typedef StaggeredKernels<Impl> Kernels;
|
||||
|
||||
///////////////////////////////////////////////////////////////
|
||||
// Implement the abstract base
|
||||
///////////////////////////////////////////////////////////////
|
||||
GridBase *GaugeGrid(void) { return _FourDimGrid ;}
|
||||
GridBase *GaugeRedBlackGrid(void) { return _FourDimRedBlackGrid ;}
|
||||
GridBase *FermionGrid(void) { return _FiveDimGrid;}
|
||||
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);
|
||||
|
||||
// half checkerboard operations
|
||||
void Meooe (const FermionField &in, FermionField &out);
|
||||
void Mooee (const FermionField &in, FermionField &out);
|
||||
void MooeeInv (const FermionField &in, FermionField &out);
|
||||
|
||||
void MeooeDag (const FermionField &in, FermionField &out);
|
||||
void MooeeDag (const FermionField &in, FermionField &out);
|
||||
void MooeeInvDag (const FermionField &in, FermionField &out);
|
||||
|
||||
void Mdir (const FermionField &in, FermionField &out,int dir,int disp);
|
||||
void DhopDir(const FermionField &in, FermionField &out,int dir,int disp);
|
||||
|
||||
// These can be overridden by fancy 5d chiral action
|
||||
void DhopDeriv (GaugeField &mat,const FermionField &U,const FermionField &V,int dag);
|
||||
void DhopDerivEO(GaugeField &mat,const FermionField &U,const FermionField &V,int dag);
|
||||
void DhopDerivOE(GaugeField &mat,const FermionField &U,const FermionField &V,int dag);
|
||||
|
||||
// Implement hopping term non-hermitian hopping term; half cb or both
|
||||
void Dhop (const FermionField &in, FermionField &out,int dag);
|
||||
void DhopOE(const FermionField &in, FermionField &out,int dag);
|
||||
void DhopEO(const FermionField &in, FermionField &out,int dag);
|
||||
|
||||
|
||||
///////////////////////////////////////////////////////////////
|
||||
// New methods added
|
||||
///////////////////////////////////////////////////////////////
|
||||
void DerivInternal(StencilImpl & st,
|
||||
DoubledGaugeField & U,
|
||||
DoubledGaugeField & UUU,
|
||||
GaugeField &mat,
|
||||
const FermionField &A,
|
||||
const FermionField &B,
|
||||
int dag);
|
||||
|
||||
void DhopInternal(StencilImpl & st,
|
||||
LebesgueOrder &lo,
|
||||
DoubledGaugeField &U,
|
||||
DoubledGaugeField &UUU,
|
||||
const FermionField &in,
|
||||
FermionField &out,
|
||||
int dag);
|
||||
|
||||
// Constructors
|
||||
ImprovedStaggeredFermion5D(GaugeField &_Uthin,
|
||||
GaugeField &_Ufat,
|
||||
GridCartesian &FiveDimGrid,
|
||||
GridRedBlackCartesian &FiveDimRedBlackGrid,
|
||||
GridCartesian &FourDimGrid,
|
||||
GridRedBlackCartesian &FourDimRedBlackGrid,
|
||||
double _mass,
|
||||
RealD _c1=9.0/8.0, RealD _c2=-1.0/24.0,RealD _u0=1.0,
|
||||
const ImplParams &p= ImplParams());
|
||||
|
||||
// DoubleStore
|
||||
void ImportGauge(const GaugeField &_U);
|
||||
void ImportGauge(const GaugeField &_Uthin,const GaugeField &_Ufat);
|
||||
|
||||
///////////////////////////////////////////////////////////////
|
||||
// Data members require to support the functionality
|
||||
///////////////////////////////////////////////////////////////
|
||||
public:
|
||||
|
||||
GridBase *_FourDimGrid;
|
||||
GridBase *_FourDimRedBlackGrid;
|
||||
GridBase *_FiveDimGrid;
|
||||
GridBase *_FiveDimRedBlackGrid;
|
||||
|
||||
RealD mass;
|
||||
RealD c1;
|
||||
RealD c2;
|
||||
RealD u0;
|
||||
int Ls;
|
||||
|
||||
//Defines the stencils for even and odd
|
||||
StencilImpl Stencil;
|
||||
StencilImpl StencilEven;
|
||||
StencilImpl StencilOdd;
|
||||
|
||||
// Copy of the gauge field , with even and odd subsets
|
||||
DoubledGaugeField Umu;
|
||||
DoubledGaugeField UmuEven;
|
||||
DoubledGaugeField UmuOdd;
|
||||
|
||||
DoubledGaugeField UUUmu;
|
||||
DoubledGaugeField UUUmuEven;
|
||||
DoubledGaugeField UUUmuOdd;
|
||||
|
||||
LebesgueOrder Lebesgue;
|
||||
LebesgueOrder LebesgueEvenOdd;
|
||||
|
||||
// Comms buffer
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > comm_buf;
|
||||
|
||||
};
|
||||
|
||||
}}
|
||||
|
||||
#endif
|
||||
291
tests/core/Test_staggered.cc
Normal file
291
tests/core/Test_staggered.cc
Normal file
@ -0,0 +1,291 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./benchmarks/Benchmark_wilson.cc
|
||||
|
||||
Copyright (C) 2015
|
||||
|
||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
|
||||
This program is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License along
|
||||
with this program; if not, write to the Free Software Foundation, Inc.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
|
||||
See the full license in the file "LICENSE" in the top level distribution directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#include <Grid/Grid.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
std::vector<int> latt_size = GridDefaultLatt();
|
||||
std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
|
||||
std::vector<int> mpi_layout = GridDefaultMpi();
|
||||
GridCartesian Grid(latt_size,simd_layout,mpi_layout);
|
||||
GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout);
|
||||
|
||||
int threads = GridThread::GetThreads();
|
||||
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
|
||||
std::cout<<GridLogMessage << "Grid floating point word size is REALF"<< sizeof(RealF)<<std::endl;
|
||||
std::cout<<GridLogMessage << "Grid floating point word size is REALD"<< sizeof(RealD)<<std::endl;
|
||||
std::cout<<GridLogMessage << "Grid floating point word size is REAL"<< sizeof(Real)<<std::endl;
|
||||
|
||||
std::vector<int> seeds({1,2,3,4});
|
||||
GridParallelRNG pRNG(&Grid);
|
||||
pRNG.SeedFixedIntegers(seeds);
|
||||
// pRNG.SeedRandomDevice();
|
||||
|
||||
typedef typename ImprovedStaggeredFermionR::FermionField FermionField;
|
||||
typedef typename ImprovedStaggeredFermionR::ComplexField ComplexField;
|
||||
typename ImprovedStaggeredFermionR::ImplParams params;
|
||||
|
||||
FermionField src (&Grid); random(pRNG,src);
|
||||
FermionField result(&Grid); result=zero;
|
||||
FermionField ref(&Grid); ref=zero;
|
||||
FermionField tmp(&Grid); tmp=zero;
|
||||
FermionField err(&Grid); tmp=zero;
|
||||
FermionField phi (&Grid); random(pRNG,phi);
|
||||
FermionField chi (&Grid); random(pRNG,chi);
|
||||
LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu);
|
||||
std::vector<LatticeColourMatrix> U(4,&Grid);
|
||||
|
||||
|
||||
double volume=1;
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
volume=volume*latt_size[mu];
|
||||
}
|
||||
|
||||
// Only one non-zero (y)
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
|
||||
/* Debug force unit
|
||||
U[mu] = 1.0;
|
||||
PokeIndex<LorentzIndex>(Umu,U[mu],mu);
|
||||
*/
|
||||
}
|
||||
|
||||
ref = zero;
|
||||
|
||||
RealD mass=0.1;
|
||||
RealD c1=9.0/8.0;
|
||||
RealD c2=-1.0/24.0;
|
||||
RealD u0=1.0;
|
||||
|
||||
{ // Simple improved staggered implementation
|
||||
ref = zero;
|
||||
RealD c1tad = 0.5*c1/u0;
|
||||
RealD c2tad = 0.5*c2/u0/u0/u0;
|
||||
|
||||
Lattice<iScalar<vInteger> > coor(&Grid);
|
||||
|
||||
Lattice<iScalar<vInteger> > x(&Grid); LatticeCoordinate(x,0);
|
||||
Lattice<iScalar<vInteger> > y(&Grid); LatticeCoordinate(y,1);
|
||||
Lattice<iScalar<vInteger> > z(&Grid); LatticeCoordinate(z,2);
|
||||
Lattice<iScalar<vInteger> > t(&Grid); LatticeCoordinate(t,3);
|
||||
|
||||
Lattice<iScalar<vInteger> > lin_z(&Grid); lin_z=x+y;
|
||||
Lattice<iScalar<vInteger> > lin_t(&Grid); lin_t=x+y+z;
|
||||
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
|
||||
// Staggered Phase.
|
||||
ComplexField phases(&Grid); phases=1.0;
|
||||
|
||||
if ( mu == 1 ) phases = where( mod(x ,2)==(Integer)0, phases,-phases);
|
||||
if ( mu == 2 ) phases = where( mod(lin_z,2)==(Integer)0, phases,-phases);
|
||||
if ( mu == 3 ) phases = where( mod(lin_t,2)==(Integer)0, phases,-phases);
|
||||
|
||||
tmp = PeriodicBC::CovShiftForward(U[mu],mu,src);
|
||||
ref = ref +c1tad*tmp*phases; // Forward 1 hop
|
||||
|
||||
tmp = PeriodicBC::CovShiftForward(U[mu],mu,tmp); // 2 hop
|
||||
tmp = PeriodicBC::CovShiftForward(U[mu],mu,tmp); // 3 hop
|
||||
ref = ref +c2tad*tmp*phases; // Forward 3 hop
|
||||
|
||||
tmp = PeriodicBC::CovShiftBackward(U[mu],mu,src);
|
||||
ref = ref -c1tad*tmp*phases; // Forward 3 hop
|
||||
|
||||
tmp = PeriodicBC::CovShiftBackward(U[mu],mu,tmp);
|
||||
tmp = PeriodicBC::CovShiftBackward(U[mu],mu,tmp);
|
||||
ref = ref -c2tad*tmp*phases; // Forward 3 hop
|
||||
}
|
||||
// ref = ref + mass * src;
|
||||
}
|
||||
|
||||
ImprovedStaggeredFermionR Ds(Umu,Umu,Grid,RBGrid,mass,c1,c2,u0,params);
|
||||
|
||||
|
||||
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
|
||||
std::cout<<GridLogMessage<<"= Testing Dhop against cshift implementation "<<std::endl;
|
||||
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
|
||||
|
||||
std::cout<<GridLogMessage << "Calling Ds"<<std::endl;
|
||||
int ncall=1000;
|
||||
double t0=usecond();
|
||||
for(int i=0;i<ncall;i++){
|
||||
Ds.Dhop(src,result,0);
|
||||
}
|
||||
double t1=usecond();
|
||||
double t2;
|
||||
double flops=(16*(3*(6+8+8)) + 15*3*2)*volume*ncall; // == 66*16 + == 1146
|
||||
|
||||
std::cout<<GridLogMessage << "Called Ds"<<std::endl;
|
||||
std::cout<<GridLogMessage << "norm result "<< norm2(result)<<std::endl;
|
||||
std::cout<<GridLogMessage << "norm ref "<< norm2(ref)<<std::endl;
|
||||
std::cout<<GridLogMessage << "mflop/s = "<< flops/(t1-t0)<<std::endl;
|
||||
|
||||
err = ref-result;
|
||||
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<<std::endl;
|
||||
|
||||
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
|
||||
std::cout<<GridLogMessage<<"= Testing that Deo + Doe = Dunprec "<<std::endl;
|
||||
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
|
||||
|
||||
FermionField src_e (&RBGrid);
|
||||
FermionField src_o (&RBGrid);
|
||||
FermionField r_e (&RBGrid);
|
||||
FermionField r_o (&RBGrid);
|
||||
FermionField r_eo (&Grid);
|
||||
pickCheckerboard(Even,src_e,src);
|
||||
pickCheckerboard(Odd,src_o,src);
|
||||
|
||||
Ds.Meooe(src_e,r_o); std::cout<<GridLogMessage<<"Applied Meo"<<std::endl;
|
||||
Ds.Meooe(src_o,r_e); std::cout<<GridLogMessage<<"Applied Moe"<<std::endl;
|
||||
Ds.Dhop (src,ref,DaggerNo);
|
||||
|
||||
setCheckerboard(r_eo,r_o);
|
||||
setCheckerboard(r_eo,r_e);
|
||||
|
||||
err= ref - r_eo;
|
||||
std::cout<<GridLogMessage << "EO norm diff "<< norm2(err)<< " "<<norm2(ref)<< " " << norm2(r_eo) <<std::endl;
|
||||
|
||||
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||
std::cout<<GridLogMessage<<"= Test Ddagger is the dagger of D by requiring "<<std::endl;
|
||||
std::cout<<GridLogMessage<<"= < phi | Deo | chi > * = < chi | Deo^dag| phi> "<<std::endl;
|
||||
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||
|
||||
FermionField chi_e (&RBGrid);
|
||||
FermionField chi_o (&RBGrid);
|
||||
|
||||
FermionField dchi_e (&RBGrid);
|
||||
FermionField dchi_o (&RBGrid);
|
||||
|
||||
FermionField phi_e (&RBGrid);
|
||||
FermionField phi_o (&RBGrid);
|
||||
|
||||
FermionField dphi_e (&RBGrid);
|
||||
FermionField dphi_o (&RBGrid);
|
||||
|
||||
pickCheckerboard(Even,chi_e,chi);
|
||||
pickCheckerboard(Odd ,chi_o,chi);
|
||||
pickCheckerboard(Even,phi_e,phi);
|
||||
pickCheckerboard(Odd ,phi_o,phi);
|
||||
|
||||
Ds.Meooe(chi_e,dchi_o);
|
||||
Ds.Meooe(chi_o,dchi_e);
|
||||
Ds.MeooeDag(phi_e,dphi_o);
|
||||
Ds.MeooeDag(phi_o,dphi_e);
|
||||
|
||||
ComplexD pDce = innerProduct(phi_e,dchi_e);
|
||||
ComplexD pDco = innerProduct(phi_o,dchi_o);
|
||||
ComplexD cDpe = innerProduct(chi_e,dphi_e);
|
||||
ComplexD cDpo = innerProduct(chi_o,dphi_o);
|
||||
|
||||
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
|
||||
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
|
||||
|
||||
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
|
||||
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
|
||||
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
|
||||
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
|
||||
|
||||
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
|
||||
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
|
||||
|
||||
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||
std::cout<<GridLogMessage<<"= Test MeeInv Mee = 1 "<<std::endl;
|
||||
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||
|
||||
pickCheckerboard(Even,chi_e,chi);
|
||||
pickCheckerboard(Odd ,chi_o,chi);
|
||||
|
||||
Ds.Mooee(chi_e,src_e);
|
||||
Ds.MooeeInv(src_e,phi_e);
|
||||
|
||||
Ds.Mooee(chi_o,src_o);
|
||||
Ds.MooeeInv(src_o,phi_o);
|
||||
|
||||
setCheckerboard(phi,phi_e);
|
||||
setCheckerboard(phi,phi_o);
|
||||
|
||||
err = phi-chi;
|
||||
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
|
||||
|
||||
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||
std::cout<<GridLogMessage<<"= Test MeeInvDag MeeDag = 1 "<<std::endl;
|
||||
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||
|
||||
pickCheckerboard(Even,chi_e,chi);
|
||||
pickCheckerboard(Odd ,chi_o,chi);
|
||||
|
||||
Ds.MooeeDag(chi_e,src_e);
|
||||
Ds.MooeeInvDag(src_e,phi_e);
|
||||
|
||||
Ds.MooeeDag(chi_o,src_o);
|
||||
Ds.MooeeInvDag(src_o,phi_o);
|
||||
|
||||
setCheckerboard(phi,phi_e);
|
||||
setCheckerboard(phi,phi_o);
|
||||
|
||||
err = phi-chi;
|
||||
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
|
||||
|
||||
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||
std::cout<<GridLogMessage<<"= Test MpcDagMpc is Hermitian "<<std::endl;
|
||||
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||
|
||||
random(pRNG,phi);
|
||||
random(pRNG,chi);
|
||||
pickCheckerboard(Even,chi_e,chi);
|
||||
pickCheckerboard(Odd ,chi_o,chi);
|
||||
pickCheckerboard(Even,phi_e,phi);
|
||||
pickCheckerboard(Odd ,phi_o,phi);
|
||||
|
||||
SchurDiagMooeeOperator<ImprovedStaggeredFermionR,FermionField> HermOpEO(Ds);
|
||||
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2);
|
||||
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2);
|
||||
|
||||
HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2);
|
||||
HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2);
|
||||
|
||||
pDce = innerProduct(phi_e,dchi_e);
|
||||
pDco = innerProduct(phi_o,dchi_o);
|
||||
cDpe = innerProduct(chi_e,dphi_e);
|
||||
cDpo = innerProduct(chi_o,dphi_o);
|
||||
|
||||
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
|
||||
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
|
||||
|
||||
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDco-conj(cDpo) <<std::endl;
|
||||
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDce-conj(cDpe) <<std::endl;
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
314
tests/core/Test_staggered5D.cc
Normal file
314
tests/core/Test_staggered5D.cc
Normal file
@ -0,0 +1,314 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./benchmarks/Benchmark_wilson.cc
|
||||
|
||||
Copyright (C) 2015
|
||||
|
||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
|
||||
This program is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License along
|
||||
with this program; if not, write to the Free Software Foundation, Inc.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
|
||||
See the full license in the file "LICENSE" in the top level distribution directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#include <Grid/Grid.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
std::vector<int> latt_size = GridDefaultLatt();
|
||||
std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
|
||||
std::vector<int> mpi_layout = GridDefaultMpi();
|
||||
|
||||
std::cout << GridLogMessage << "Making s innermost grids"<<std::endl;
|
||||
|
||||
const int Ls=16;
|
||||
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
|
||||
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
|
||||
|
||||
int threads = GridThread::GetThreads();
|
||||
|
||||
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
|
||||
|
||||
std::cout<<GridLogMessage << "Grid floating point word size is REALF"<< sizeof(RealF)<<std::endl;
|
||||
std::cout<<GridLogMessage << "Grid floating point word size is REALD"<< sizeof(RealD)<<std::endl;
|
||||
std::cout<<GridLogMessage << "Grid floating point word size is REAL"<< sizeof(Real)<<std::endl;
|
||||
|
||||
std::vector<int> seeds({1,2,3,4});
|
||||
GridParallelRNG pRNG4(UGrid);
|
||||
GridParallelRNG pRNG5(FGrid);
|
||||
pRNG4.SeedFixedIntegers(seeds);
|
||||
pRNG5.SeedFixedIntegers(seeds);
|
||||
|
||||
typedef typename ImprovedStaggeredFermion5DR::FermionField FermionField;
|
||||
typedef typename ImprovedStaggeredFermion5DR::ComplexField ComplexField;
|
||||
typename ImprovedStaggeredFermion5DR::ImplParams params;
|
||||
|
||||
FermionField src (FGrid); random(pRNG5,src);
|
||||
FermionField result(FGrid); result=zero;
|
||||
FermionField ref(FGrid); ref=zero;
|
||||
FermionField tmp(FGrid); tmp=zero;
|
||||
FermionField err(FGrid); tmp=zero;
|
||||
FermionField phi (FGrid); random(pRNG5,phi);
|
||||
FermionField chi (FGrid); random(pRNG5,chi);
|
||||
|
||||
LatticeGaugeField Umu(UGrid); SU3::ColdConfiguration(pRNG4,Umu);
|
||||
|
||||
double volume=Ls;
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
volume=volume*latt_size[mu];
|
||||
}
|
||||
|
||||
////////////////////////////////////
|
||||
// Naive implementation needs to
|
||||
// replicate across fifth dimension
|
||||
////////////////////////////////////
|
||||
LatticeGaugeField Umu5d(FGrid);
|
||||
for(int ss=0;ss<Umu._grid->oSites();ss++){
|
||||
for(int s=0;s<Ls;s++){
|
||||
Umu5d._odata[Ls*ss+s] = Umu._odata[ss];
|
||||
}
|
||||
}
|
||||
|
||||
std::vector<LatticeColourMatrix> U(4,FGrid);
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
U[mu] = PeekIndex<LorentzIndex>(Umu5d,mu);
|
||||
if ( mu!=0 ) U[mu]=zero;
|
||||
PokeIndex<LorentzIndex>(Umu5d,U[mu],mu);
|
||||
}
|
||||
|
||||
std::vector<LatticeColourMatrix> Ua(4,UGrid);
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
Ua[mu] = PeekIndex<LorentzIndex>(Umu,mu);
|
||||
if ( mu!=0 ) {
|
||||
Ua[mu]=zero;
|
||||
}
|
||||
PokeIndex<LorentzIndex>(Umu,Ua[mu],mu);
|
||||
}
|
||||
|
||||
RealD mass=0.1;
|
||||
RealD c1=9.0/8.0;
|
||||
RealD c2=-1.0/24.0;
|
||||
RealD u0=1.0;
|
||||
|
||||
{ // Simple improved staggered implementation
|
||||
ref = zero;
|
||||
RealD c1tad = 0.5*c1/u0;
|
||||
RealD c2tad = 0.5*c2/u0/u0/u0;
|
||||
|
||||
Lattice<iScalar<vInteger> > coor(FGrid);
|
||||
|
||||
Lattice<iScalar<vInteger> > x(FGrid); LatticeCoordinate(x,1); // s innermost
|
||||
Lattice<iScalar<vInteger> > y(FGrid); LatticeCoordinate(y,2);
|
||||
Lattice<iScalar<vInteger> > z(FGrid); LatticeCoordinate(z,3);
|
||||
Lattice<iScalar<vInteger> > t(FGrid); LatticeCoordinate(t,4);
|
||||
|
||||
Lattice<iScalar<vInteger> > lin_z(FGrid); lin_z=x+y;
|
||||
Lattice<iScalar<vInteger> > lin_t(FGrid); lin_t=x+y+z;
|
||||
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
|
||||
// Staggered Phase.
|
||||
ComplexField phases(FGrid); phases=1.0;
|
||||
|
||||
if ( mu == 1 ) phases = where( mod(x ,2)==(Integer)0, phases,-phases);
|
||||
if ( mu == 2 ) phases = where( mod(lin_z,2)==(Integer)0, phases,-phases);
|
||||
if ( mu == 3 ) phases = where( mod(lin_t,2)==(Integer)0, phases,-phases);
|
||||
|
||||
tmp = PeriodicBC::CovShiftForward(U[mu],mu+1,src);
|
||||
ref = ref +c1tad*tmp*phases; // Forward 1 hop
|
||||
|
||||
tmp = PeriodicBC::CovShiftForward(U[mu],mu+1,tmp); // 2 hop
|
||||
tmp = PeriodicBC::CovShiftForward(U[mu],mu+1,tmp); // 3 hop
|
||||
ref = ref +c2tad*tmp*phases; // Forward 3 hop
|
||||
|
||||
tmp = PeriodicBC::CovShiftBackward(U[mu],mu+1,src);
|
||||
ref = ref -c1tad*tmp*phases; // Forward 3 hop
|
||||
|
||||
tmp = PeriodicBC::CovShiftBackward(U[mu],mu+1,tmp);
|
||||
tmp = PeriodicBC::CovShiftBackward(U[mu],mu+1,tmp);
|
||||
ref = ref -c2tad*tmp*phases; // Forward 3 hop
|
||||
}
|
||||
}
|
||||
|
||||
ImprovedStaggeredFermion5DR Ds(Umu,Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,c1,c2,u0,params);
|
||||
|
||||
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
|
||||
std::cout<<GridLogMessage<<"= Testing Dhop against cshift implementation "<<std::endl;
|
||||
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
|
||||
|
||||
std::cout<<GridLogMessage << "Calling Ds"<<std::endl;
|
||||
int ncall=1000;
|
||||
double t0=usecond();
|
||||
for(int i=0;i<ncall;i++){
|
||||
Ds.Dhop(src,result,0);
|
||||
}
|
||||
double t1=usecond();
|
||||
|
||||
double t2;
|
||||
double flops=(16*(3*(6+8+8)) + 15*3*2)*volume*ncall; // == 66*16 + == 1146
|
||||
|
||||
std::cout<<GridLogMessage << "Called Ds"<<std::endl;
|
||||
std::cout<<GridLogMessage << "norm result "<< norm2(result)<<std::endl;
|
||||
std::cout<<GridLogMessage << "norm ref "<< norm2(ref)<<std::endl;
|
||||
std::cout<<GridLogMessage << "mflop/s = "<< flops/(t1-t0)<<std::endl;
|
||||
|
||||
err = ref-result;
|
||||
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<<std::endl;
|
||||
|
||||
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
|
||||
std::cout<<GridLogMessage<<"= Testing that Deo + Doe = Dunprec "<<std::endl;
|
||||
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
|
||||
|
||||
FermionField src_e (FrbGrid);
|
||||
FermionField src_o (FrbGrid);
|
||||
FermionField r_e (FrbGrid);
|
||||
FermionField r_o (FrbGrid);
|
||||
FermionField r_eo (FGrid);
|
||||
pickCheckerboard(Even,src_e,src);
|
||||
pickCheckerboard(Odd,src_o,src);
|
||||
|
||||
Ds.Meooe(src_e,r_o); std::cout<<GridLogMessage<<"Applied Meo"<<std::endl;
|
||||
Ds.Meooe(src_o,r_e); std::cout<<GridLogMessage<<"Applied Moe"<<std::endl;
|
||||
Ds.Dhop (src,ref,DaggerNo);
|
||||
|
||||
setCheckerboard(r_eo,r_o);
|
||||
setCheckerboard(r_eo,r_e);
|
||||
|
||||
err= ref - r_eo;
|
||||
std::cout<<GridLogMessage << "EO norm diff "<< norm2(err)<< " "<<norm2(ref)<< " " << norm2(r_eo) <<std::endl;
|
||||
|
||||
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||
std::cout<<GridLogMessage<<"= Test Ddagger is the dagger of D by requiring "<<std::endl;
|
||||
std::cout<<GridLogMessage<<"= < phi | Deo | chi > * = < chi | Deo^dag| phi> "<<std::endl;
|
||||
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||
|
||||
FermionField chi_e (FrbGrid);
|
||||
FermionField chi_o (FrbGrid);
|
||||
|
||||
FermionField dchi_e (FrbGrid);
|
||||
FermionField dchi_o (FrbGrid);
|
||||
|
||||
FermionField phi_e (FrbGrid);
|
||||
FermionField phi_o (FrbGrid);
|
||||
|
||||
FermionField dphi_e (FrbGrid);
|
||||
FermionField dphi_o (FrbGrid);
|
||||
|
||||
pickCheckerboard(Even,chi_e,chi);
|
||||
pickCheckerboard(Odd ,chi_o,chi);
|
||||
pickCheckerboard(Even,phi_e,phi);
|
||||
pickCheckerboard(Odd ,phi_o,phi);
|
||||
|
||||
Ds.Meooe(chi_e,dchi_o);
|
||||
Ds.Meooe(chi_o,dchi_e);
|
||||
Ds.MeooeDag(phi_e,dphi_o);
|
||||
Ds.MeooeDag(phi_o,dphi_e);
|
||||
|
||||
ComplexD pDce = innerProduct(phi_e,dchi_e);
|
||||
ComplexD pDco = innerProduct(phi_o,dchi_o);
|
||||
ComplexD cDpe = innerProduct(chi_e,dphi_e);
|
||||
ComplexD cDpo = innerProduct(chi_o,dphi_o);
|
||||
|
||||
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
|
||||
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
|
||||
|
||||
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
|
||||
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
|
||||
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
|
||||
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
|
||||
|
||||
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
|
||||
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
|
||||
|
||||
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||
std::cout<<GridLogMessage<<"= Test MeeInv Mee = 1 "<<std::endl;
|
||||
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||
|
||||
pickCheckerboard(Even,chi_e,chi);
|
||||
pickCheckerboard(Odd ,chi_o,chi);
|
||||
|
||||
Ds.Mooee(chi_e,src_e);
|
||||
Ds.MooeeInv(src_e,phi_e);
|
||||
|
||||
Ds.Mooee(chi_o,src_o);
|
||||
Ds.MooeeInv(src_o,phi_o);
|
||||
|
||||
setCheckerboard(phi,phi_e);
|
||||
setCheckerboard(phi,phi_o);
|
||||
|
||||
err = phi-chi;
|
||||
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
|
||||
|
||||
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||
std::cout<<GridLogMessage<<"= Test MeeInvDag MeeDag = 1 "<<std::endl;
|
||||
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||
|
||||
pickCheckerboard(Even,chi_e,chi);
|
||||
pickCheckerboard(Odd ,chi_o,chi);
|
||||
|
||||
Ds.MooeeDag(chi_e,src_e);
|
||||
Ds.MooeeInvDag(src_e,phi_e);
|
||||
|
||||
Ds.MooeeDag(chi_o,src_o);
|
||||
Ds.MooeeInvDag(src_o,phi_o);
|
||||
|
||||
setCheckerboard(phi,phi_e);
|
||||
setCheckerboard(phi,phi_o);
|
||||
|
||||
err = phi-chi;
|
||||
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
|
||||
|
||||
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||
std::cout<<GridLogMessage<<"= Test MpcDagMpc is Hermitian "<<std::endl;
|
||||
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||
|
||||
random(pRNG5,phi);
|
||||
random(pRNG5,chi);
|
||||
pickCheckerboard(Even,chi_e,chi);
|
||||
pickCheckerboard(Odd ,chi_o,chi);
|
||||
pickCheckerboard(Even,phi_e,phi);
|
||||
pickCheckerboard(Odd ,phi_o,phi);
|
||||
|
||||
SchurDiagMooeeOperator<ImprovedStaggeredFermion5DR,FermionField> HermOpEO(Ds);
|
||||
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2);
|
||||
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2);
|
||||
|
||||
HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2);
|
||||
HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2);
|
||||
|
||||
pDce = innerProduct(phi_e,dchi_e);
|
||||
pDco = innerProduct(phi_o,dchi_o);
|
||||
cDpe = innerProduct(chi_e,dphi_e);
|
||||
cDpo = innerProduct(chi_o,dphi_o);
|
||||
|
||||
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
|
||||
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
|
||||
|
||||
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDco-conj(cDpo) <<std::endl;
|
||||
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDce-conj(cDpe) <<std::endl;
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
Loading…
x
Reference in New Issue
Block a user