1
0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 07:55:35 +00:00

Moving some things around for pretty

This commit is contained in:
Peter Boyle 2015-05-11 19:09:49 +01:00
parent a411b48a91
commit 8b765be2b1
6 changed files with 336 additions and 324 deletions

View File

@ -1,117 +1,6 @@
#ifndef GRID_COMMUNICATOR_H #ifndef GRID_COMMUNICATOR_H
#define GRID_COMMUNICATOR_H #define GRID_COMMUNICATOR_H
/////////////////////////////////// #include <communicator/Grid_communicator_base.h>
// Processor layout information
///////////////////////////////////
#ifdef GRID_COMMS_MPI
#include <mpi.h>
#endif
namespace Grid {
class CartesianCommunicator {
public:
// Communicator should know nothing of the physics grid, only processor grid.
int _Nprocessors; // How many in all
std::vector<int> _processors; // Which dimensions get relayed out over processors lanes.
int _processor; // linear processor rank
std::vector<int> _processor_coor; // linear processor coordinate
unsigned long _ndimension;
#ifdef GRID_COMMS_MPI
MPI_Comm communicator;
typedef MPI_Request CommsRequest_t;
#else
typedef int CommsRequest_t;
#endif
// Constructor
CartesianCommunicator(std::vector<int> &pdimensions_in);
// Wraps MPI_Cart routines
void ShiftedRanks(int dim,int shift,int & source, int & dest);
int RankFromProcessorCoor(std::vector<int> &coor);
void ProcessorCoorFromRank(int rank,std::vector<int> &coor);
/////////////////////////////////
// Grid information queries
/////////////////////////////////
int IsBoss(void) { return _processor==0; };
int BossRank(void) { return 0; };
int ThisRank(void) { return _processor; };
const std::vector<int> & ThisProcessorCoor(void) { return _processor_coor; };
const std::vector<int> & ProcessorGrid(void) { return _processors; };
int ProcessorCount(void) { return _Nprocessors; };
////////////////////////////////////////////////////////////
// Reduction
////////////////////////////////////////////////////////////
void GlobalSum(RealF &);
void GlobalSumVector(RealF *,int N);
void GlobalSum(RealD &);
void GlobalSumVector(RealD *,int N);
void GlobalSum(uint32_t &);
void GlobalSum(ComplexF &c)
{
GlobalSumVector((float *)&c,2);
}
void GlobalSumVector(ComplexF *c,int N)
{
GlobalSumVector((float *)c,2*N);
}
void GlobalSum(ComplexD &c)
{
GlobalSumVector((double *)&c,2);
}
void GlobalSumVector(ComplexD *c,int N)
{
GlobalSumVector((double *)c,2*N);
}
template<class obj> void GlobalSum(obj &o){
typedef typename obj::scalar_type scalar_type;
int words = sizeof(obj)/sizeof(scalar_type);
scalar_type * ptr = (scalar_type *)& o;
GlobalSumVector(ptr,words);
}
////////////////////////////////////////////////////////////
// Face exchange, buffer swap in translational invariant way
////////////////////////////////////////////////////////////
void SendToRecvFrom(void *xmit,
int xmit_to_rank,
void *recv,
int recv_from_rank,
int bytes);
void SendToRecvFromBegin(std::vector<CommsRequest_t> &list,
void *xmit,
int xmit_to_rank,
void *recv,
int recv_from_rank,
int bytes);
void SendToRecvFromComplete(std::vector<CommsRequest_t> &waitall);
////////////////////////////////////////////////////////////
// Barrier
////////////////////////////////////////////////////////////
void Barrier(void);
////////////////////////////////////////////////////////////
// Broadcast a buffer and composite larger
////////////////////////////////////////////////////////////
void Broadcast(int root,void* data, int bytes);
template<class obj> void Broadcast(int root,obj &data)
{
Broadcast(root,(void *)&data,sizeof(data));
};
static void BroadcastWorld(int root,void* data, int bytes);
};
}
#endif #endif

View File

@ -1,203 +1,6 @@
#ifndef GRID_LATTICE_H #ifndef GRID_LATTICE_H
#define GRID_LATTICE_H #define GRID_LATTICE_H
namespace Grid { #include <lattice/Grid_lattice_base.h>
// TODO:
// mac,real,imag
// Functionality:
// -=,+=,*=,()
// add,+,sub,-,mult,mac,*
// adj,conj
// real,imag
// transpose,transposeIndex
// trace,traceIndex
// peekIndex
// innerProduct,outerProduct,
// localNorm2
// localInnerProduct
extern int GridCshiftPermuteMap[4][16];
////////////////////////////////////////////////
// Basic expressions used in Expression Template
////////////////////////////////////////////////
class LatticeBase {};
class LatticeExpressionBase {};
template <typename Op, typename T1>
class LatticeUnaryExpression : public std::pair<Op,std::tuple<T1> > , public LatticeExpressionBase {
public:
LatticeUnaryExpression(const std::pair<Op,std::tuple<T1> > &arg): std::pair<Op,std::tuple<T1> >(arg) {};
};
template <typename Op, typename T1, typename T2>
class LatticeBinaryExpression : public std::pair<Op,std::tuple<T1,T2> > , public LatticeExpressionBase {
public:
LatticeBinaryExpression(const std::pair<Op,std::tuple<T1,T2> > &arg): std::pair<Op,std::tuple<T1,T2> >(arg) {};
};
template <typename Op, typename T1, typename T2, typename T3>
class LatticeTrinaryExpression :public std::pair<Op,std::tuple<T1,T2,T3> >, public LatticeExpressionBase {
public:
LatticeTrinaryExpression(const std::pair<Op,std::tuple<T1,T2,T3> > &arg): std::pair<Op,std::tuple<T1,T2,T3> >(arg) {};
};
template<class vobj>
class Lattice : public LatticeBase
{
public:
GridBase *_grid;
int checkerboard;
std::vector<vobj,alignedAllocator<vobj> > _odata;
public:
typedef typename vobj::scalar_type scalar_type;
typedef typename vobj::vector_type vector_type;
typedef vobj vector_object;
////////////////////////////////////////////////////////////////////////////////
// Expression Template closure support
////////////////////////////////////////////////////////////////////////////////
template <typename Op, typename T1> inline Lattice<vobj> & operator=(const LatticeUnaryExpression<Op,T1> &expr)
{
#pragma omp parallel for
for(int ss=0;ss<_grid->oSites();ss++){
vobj tmp= eval(ss,expr);
vstream(_odata[ss] ,tmp);
}
return *this;
}
template <typename Op, typename T1,typename T2> inline Lattice<vobj> & operator=(const LatticeBinaryExpression<Op,T1,T2> &expr)
{
#pragma omp parallel for
for(int ss=0;ss<_grid->oSites();ss++){
vobj tmp= eval(ss,expr);
vstream(_odata[ss] ,tmp);
}
return *this;
}
template <typename Op, typename T1,typename T2,typename T3> inline Lattice<vobj> & operator=(const LatticeTrinaryExpression<Op,T1,T2,T3> &expr)
{
#pragma omp parallel for
for(int ss=0;ss<_grid->oSites();ss++){
vobj tmp= eval(ss,expr);
vstream(_odata[ss] ,tmp);
}
return *this;
}
//GridFromExpression is tricky to do
template<class Op,class T1>
Lattice(const LatticeUnaryExpression<Op,T1> & expr): _grid(nullptr){
GridFromExpression(_grid,expr);
assert(_grid!=nullptr);
_odata.resize(_grid->oSites());
#pragma omp parallel for
for(int ss=0;ss<_grid->oSites();ss++){
_odata[ss] = eval(ss,expr);
}
};
template<class Op,class T1, class T2>
Lattice(const LatticeBinaryExpression<Op,T1,T2> & expr): _grid(nullptr){
GridFromExpression(_grid,expr);
assert(_grid!=nullptr);
_odata.resize(_grid->oSites());
#pragma omp parallel for
for(int ss=0;ss<_grid->oSites();ss++){
_odata[ss] = eval(ss,expr);
}
};
template<class Op,class T1, class T2, class T3>
Lattice(const LatticeTrinaryExpression<Op,T1,T2,T3> & expr): _grid(nullptr){
GridFromExpression(_grid,expr);
assert(_grid!=nullptr);
_odata.resize(_grid->oSites());
#pragma omp parallel for
for(int ss=0;ss<_grid->oSites();ss++){
_odata[ss] = eval(ss,expr);
}
};
//////////////////////////////////////////////////////////////////
// Constructor requires "grid" passed.
// what about a default grid?
//////////////////////////////////////////////////////////////////
Lattice(GridBase *grid) : _grid(grid), _odata(_grid->oSites()) {
// _odata.reserve(_grid->oSites());
// _odata.resize(_grid->oSites());
assert((((uint64_t)&_odata[0])&0xF) ==0);
checkerboard=0;
}
template<class sobj> inline Lattice<vobj> & operator = (const sobj & r){
#pragma omp parallel for
for(int ss=0;ss<_grid->oSites();ss++){
this->_odata[ss]=r;
}
return *this;
}
template<class robj> inline Lattice<vobj> & operator = (const Lattice<robj> & r){
conformable(*this,r);
std::cout<<"Lattice operator ="<<std::endl;
#pragma omp parallel for
for(int ss=0;ss<_grid->oSites();ss++){
this->_odata[ss]=r._odata[ss];
}
return *this;
}
// *=,+=,-= operators inherit behvour from correspond */+/- operation
template<class T> inline Lattice<vobj> &operator *=(const T &r) {
*this = (*this)*r;
return *this;
}
template<class T> inline Lattice<vobj> &operator -=(const T &r) {
*this = (*this)-r;
return *this;
}
template<class T> inline Lattice<vobj> &operator +=(const T &r) {
*this = (*this)+r;
return *this;
}
inline friend Lattice<vobj> operator / (const Lattice<vobj> &lhs,const Lattice<vobj> &rhs){
conformable(lhs,rhs);
Lattice<vobj> ret(lhs._grid);
#pragma omp parallel for
for(int ss=0;ss<lhs._grid->oSites();ss++){
ret._odata[ss] = lhs._odata[ss]/rhs._odata[ss];
}
return ret;
};
}; // class Lattice
}
#undef GRID_LATTICE_EXPRESSION_TEMPLATES
#include <lattice/Grid_lattice_conformable.h>
#ifdef GRID_LATTICE_EXPRESSION_TEMPLATES
#include <lattice/Grid_lattice_ET.h>
#else
#include <lattice/Grid_lattice_overload.h>
#endif
#include <lattice/Grid_lattice_arith.h>
#include <lattice/Grid_lattice_trace.h>
#include <lattice/Grid_lattice_transpose.h>
#include <lattice/Grid_lattice_local.h>
#include <lattice/Grid_lattice_reduction.h>
#include <lattice/Grid_lattice_peekpoke.h>
#include <lattice/Grid_lattice_reality.h>
#include <Grid_extract.h>
#include <lattice/Grid_lattice_coordinate.h>
#include <lattice/Grid_lattice_rng.h>
#include <lattice/Grid_lattice_transfer.h>
#endif #endif

View File

@ -0,0 +1,117 @@
#ifndef GRID_COMMUNICATOR_BASE_H
#define GRID_COMMUNICATOR_BASE_H
///////////////////////////////////
// Processor layout information
///////////////////////////////////
#ifdef GRID_COMMS_MPI
#include <mpi.h>
#endif
namespace Grid {
class CartesianCommunicator {
public:
// Communicator should know nothing of the physics grid, only processor grid.
int _Nprocessors; // How many in all
std::vector<int> _processors; // Which dimensions get relayed out over processors lanes.
int _processor; // linear processor rank
std::vector<int> _processor_coor; // linear processor coordinate
unsigned long _ndimension;
#ifdef GRID_COMMS_MPI
MPI_Comm communicator;
typedef MPI_Request CommsRequest_t;
#else
typedef int CommsRequest_t;
#endif
// Constructor
CartesianCommunicator(std::vector<int> &pdimensions_in);
// Wraps MPI_Cart routines
void ShiftedRanks(int dim,int shift,int & source, int & dest);
int RankFromProcessorCoor(std::vector<int> &coor);
void ProcessorCoorFromRank(int rank,std::vector<int> &coor);
/////////////////////////////////
// Grid information queries
/////////////////////////////////
int IsBoss(void) { return _processor==0; };
int BossRank(void) { return 0; };
int ThisRank(void) { return _processor; };
const std::vector<int> & ThisProcessorCoor(void) { return _processor_coor; };
const std::vector<int> & ProcessorGrid(void) { return _processors; };
int ProcessorCount(void) { return _Nprocessors; };
////////////////////////////////////////////////////////////
// Reduction
////////////////////////////////////////////////////////////
void GlobalSum(RealF &);
void GlobalSumVector(RealF *,int N);
void GlobalSum(RealD &);
void GlobalSumVector(RealD *,int N);
void GlobalSum(uint32_t &);
void GlobalSum(ComplexF &c)
{
GlobalSumVector((float *)&c,2);
}
void GlobalSumVector(ComplexF *c,int N)
{
GlobalSumVector((float *)c,2*N);
}
void GlobalSum(ComplexD &c)
{
GlobalSumVector((double *)&c,2);
}
void GlobalSumVector(ComplexD *c,int N)
{
GlobalSumVector((double *)c,2*N);
}
template<class obj> void GlobalSum(obj &o){
typedef typename obj::scalar_type scalar_type;
int words = sizeof(obj)/sizeof(scalar_type);
scalar_type * ptr = (scalar_type *)& o;
GlobalSumVector(ptr,words);
}
////////////////////////////////////////////////////////////
// Face exchange, buffer swap in translational invariant way
////////////////////////////////////////////////////////////
void SendToRecvFrom(void *xmit,
int xmit_to_rank,
void *recv,
int recv_from_rank,
int bytes);
void SendToRecvFromBegin(std::vector<CommsRequest_t> &list,
void *xmit,
int xmit_to_rank,
void *recv,
int recv_from_rank,
int bytes);
void SendToRecvFromComplete(std::vector<CommsRequest_t> &waitall);
////////////////////////////////////////////////////////////
// Barrier
////////////////////////////////////////////////////////////
void Barrier(void);
////////////////////////////////////////////////////////////
// Broadcast a buffer and composite larger
////////////////////////////////////////////////////////////
void Broadcast(int root,void* data, int bytes);
template<class obj> void Broadcast(int root,obj &data)
{
Broadcast(root,(void *)&data,sizeof(data));
};
static void BroadcastWorld(int root,void* data, int bytes);
};
}
#endif

View File

@ -0,0 +1,203 @@
#ifndef GRID_LATTICE_BASE_H
#define GRID_LATTICE_BASE_H
namespace Grid {
// TODO:
// mac,real,imag
// Functionality:
// -=,+=,*=,()
// add,+,sub,-,mult,mac,*
// adj,conj
// real,imag
// transpose,transposeIndex
// trace,traceIndex
// peekIndex
// innerProduct,outerProduct,
// localNorm2
// localInnerProduct
extern int GridCshiftPermuteMap[4][16];
////////////////////////////////////////////////
// Basic expressions used in Expression Template
////////////////////////////////////////////////
class LatticeBase {};
class LatticeExpressionBase {};
template <typename Op, typename T1>
class LatticeUnaryExpression : public std::pair<Op,std::tuple<T1> > , public LatticeExpressionBase {
public:
LatticeUnaryExpression(const std::pair<Op,std::tuple<T1> > &arg): std::pair<Op,std::tuple<T1> >(arg) {};
};
template <typename Op, typename T1, typename T2>
class LatticeBinaryExpression : public std::pair<Op,std::tuple<T1,T2> > , public LatticeExpressionBase {
public:
LatticeBinaryExpression(const std::pair<Op,std::tuple<T1,T2> > &arg): std::pair<Op,std::tuple<T1,T2> >(arg) {};
};
template <typename Op, typename T1, typename T2, typename T3>
class LatticeTrinaryExpression :public std::pair<Op,std::tuple<T1,T2,T3> >, public LatticeExpressionBase {
public:
LatticeTrinaryExpression(const std::pair<Op,std::tuple<T1,T2,T3> > &arg): std::pair<Op,std::tuple<T1,T2,T3> >(arg) {};
};
template<class vobj>
class Lattice : public LatticeBase
{
public:
GridBase *_grid;
int checkerboard;
std::vector<vobj,alignedAllocator<vobj> > _odata;
public:
typedef typename vobj::scalar_type scalar_type;
typedef typename vobj::vector_type vector_type;
typedef vobj vector_object;
////////////////////////////////////////////////////////////////////////////////
// Expression Template closure support
////////////////////////////////////////////////////////////////////////////////
template <typename Op, typename T1> inline Lattice<vobj> & operator=(const LatticeUnaryExpression<Op,T1> &expr)
{
#pragma omp parallel for
for(int ss=0;ss<_grid->oSites();ss++){
vobj tmp= eval(ss,expr);
vstream(_odata[ss] ,tmp);
}
return *this;
}
template <typename Op, typename T1,typename T2> inline Lattice<vobj> & operator=(const LatticeBinaryExpression<Op,T1,T2> &expr)
{
#pragma omp parallel for
for(int ss=0;ss<_grid->oSites();ss++){
vobj tmp= eval(ss,expr);
vstream(_odata[ss] ,tmp);
}
return *this;
}
template <typename Op, typename T1,typename T2,typename T3> inline Lattice<vobj> & operator=(const LatticeTrinaryExpression<Op,T1,T2,T3> &expr)
{
#pragma omp parallel for
for(int ss=0;ss<_grid->oSites();ss++){
vobj tmp= eval(ss,expr);
vstream(_odata[ss] ,tmp);
}
return *this;
}
//GridFromExpression is tricky to do
template<class Op,class T1>
Lattice(const LatticeUnaryExpression<Op,T1> & expr): _grid(nullptr){
GridFromExpression(_grid,expr);
assert(_grid!=nullptr);
_odata.resize(_grid->oSites());
#pragma omp parallel for
for(int ss=0;ss<_grid->oSites();ss++){
_odata[ss] = eval(ss,expr);
}
};
template<class Op,class T1, class T2>
Lattice(const LatticeBinaryExpression<Op,T1,T2> & expr): _grid(nullptr){
GridFromExpression(_grid,expr);
assert(_grid!=nullptr);
_odata.resize(_grid->oSites());
#pragma omp parallel for
for(int ss=0;ss<_grid->oSites();ss++){
_odata[ss] = eval(ss,expr);
}
};
template<class Op,class T1, class T2, class T3>
Lattice(const LatticeTrinaryExpression<Op,T1,T2,T3> & expr): _grid(nullptr){
GridFromExpression(_grid,expr);
assert(_grid!=nullptr);
_odata.resize(_grid->oSites());
#pragma omp parallel for
for(int ss=0;ss<_grid->oSites();ss++){
_odata[ss] = eval(ss,expr);
}
};
//////////////////////////////////////////////////////////////////
// Constructor requires "grid" passed.
// what about a default grid?
//////////////////////////////////////////////////////////////////
Lattice(GridBase *grid) : _grid(grid), _odata(_grid->oSites()) {
// _odata.reserve(_grid->oSites());
// _odata.resize(_grid->oSites());
assert((((uint64_t)&_odata[0])&0xF) ==0);
checkerboard=0;
}
template<class sobj> inline Lattice<vobj> & operator = (const sobj & r){
#pragma omp parallel for
for(int ss=0;ss<_grid->oSites();ss++){
this->_odata[ss]=r;
}
return *this;
}
template<class robj> inline Lattice<vobj> & operator = (const Lattice<robj> & r){
conformable(*this,r);
std::cout<<"Lattice operator ="<<std::endl;
#pragma omp parallel for
for(int ss=0;ss<_grid->oSites();ss++){
this->_odata[ss]=r._odata[ss];
}
return *this;
}
// *=,+=,-= operators inherit behvour from correspond */+/- operation
template<class T> inline Lattice<vobj> &operator *=(const T &r) {
*this = (*this)*r;
return *this;
}
template<class T> inline Lattice<vobj> &operator -=(const T &r) {
*this = (*this)-r;
return *this;
}
template<class T> inline Lattice<vobj> &operator +=(const T &r) {
*this = (*this)+r;
return *this;
}
inline friend Lattice<vobj> operator / (const Lattice<vobj> &lhs,const Lattice<vobj> &rhs){
conformable(lhs,rhs);
Lattice<vobj> ret(lhs._grid);
#pragma omp parallel for
for(int ss=0;ss<lhs._grid->oSites();ss++){
ret._odata[ss] = lhs._odata[ss]/rhs._odata[ss];
}
return ret;
};
}; // class Lattice
}
#undef GRID_LATTICE_EXPRESSION_TEMPLATES
#include <lattice/Grid_lattice_conformable.h>
#ifdef GRID_LATTICE_EXPRESSION_TEMPLATES
#include <lattice/Grid_lattice_ET.h>
#else
#include <lattice/Grid_lattice_overload.h>
#endif
#include <lattice/Grid_lattice_arith.h>
#include <lattice/Grid_lattice_trace.h>
#include <lattice/Grid_lattice_transpose.h>
#include <lattice/Grid_lattice_local.h>
#include <lattice/Grid_lattice_reduction.h>
#include <lattice/Grid_lattice_peekpoke.h>
#include <lattice/Grid_lattice_reality.h>
#include <Grid_extract.h>
#include <lattice/Grid_lattice_coordinate.h>
#include <lattice/Grid_lattice_rng.h>
#include <lattice/Grid_lattice_transfer.h>
#endif

View File

@ -6,7 +6,7 @@ namespace Grid {
// innerProduct Vector x Vector -> Scalar // innerProduct Vector x Vector -> Scalar
// innerProduct Matrix x Matrix -> Scalar // innerProduct Matrix x Matrix -> Scalar
/////////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////////
template<class sobj> inline RealD norm2l(const sobj &arg){ template<class sobj> inline RealD norm2(const sobj &arg){
typedef typename sobj::scalar_type scalar; typedef typename sobj::scalar_type scalar;
decltype(innerProduct(arg,arg)) nrm; decltype(innerProduct(arg,arg)) nrm;
nrm = innerProduct(arg,arg); nrm = innerProduct(arg,arg);

View File

@ -97,7 +97,7 @@ int main (int argc, char ** argv)
for(int mu=0;mu<6;mu++){ for(int mu=0;mu<6;mu++){
result = Gamma(g[mu])* ident * Gamma(g[mu]); result = Gamma(g[mu])* ident * Gamma(g[mu]);
result = result - ident; result = result - ident;
double mag = TensorRemove(norm2l(result)); double mag = TensorRemove(norm2(result));
std::cout << list[mu]<<" " << mag<<std::endl; std::cout << list[mu]<<" " << mag<<std::endl;
} }
@ -105,7 +105,7 @@ int main (int argc, char ** argv)
for(int mu=0;mu<6;mu++){ for(int mu=0;mu<6;mu++){
result = rnd * Gamma(g[mu]); result = rnd * Gamma(g[mu]);
result = result + rnd * Gamma(g[mu+6]); result = result + rnd * Gamma(g[mu+6]);
double mag = TensorRemove(norm2l(result)); double mag = TensorRemove(norm2(result));
std::cout << list[mu]<<" " << mag<<std::endl; std::cout << list[mu]<<" " << mag<<std::endl;
} }
@ -113,7 +113,7 @@ int main (int argc, char ** argv)
for(int mu=0;mu<6;mu++){ for(int mu=0;mu<6;mu++){
result = Gamma(g[mu]) *rnd; result = Gamma(g[mu]) *rnd;
result = result + Gamma(g[mu+6])*rnd; result = result + Gamma(g[mu+6])*rnd;
double mag = TensorRemove(norm2l(result)); double mag = TensorRemove(norm2(result));
std::cout << list[mu]<<" " << mag<<std::endl; std::cout << list[mu]<<" " << mag<<std::endl;
} }
@ -127,70 +127,70 @@ int main (int argc, char ** argv)
spProjXp(hsm,rv); spProjXp(hsm,rv);
spReconXp(recon,hsm); spReconXp(recon,hsm);
full = rv + Gamma(Gamma::GammaX) *rv; full = rv + Gamma(Gamma::GammaX) *rv;
mag = TensorRemove(norm2l(full-recon)); mag = TensorRemove(norm2(full-recon));
std::cout << "Xp "<< mag<<std::endl; std::cout << "Xp "<< mag<<std::endl;
// Xm // Xm
spProjXm(hsm,rv); spProjXm(hsm,rv);
spReconXm(recon,hsm); spReconXm(recon,hsm);
full = rv - Gamma(Gamma::GammaX) *rv; full = rv - Gamma(Gamma::GammaX) *rv;
mag = TensorRemove(norm2l(full-recon)); mag = TensorRemove(norm2(full-recon));
std::cout << "Xm "<< mag<<std::endl; std::cout << "Xm "<< mag<<std::endl;
// Yp // Yp
spProjYp(hsm,rv); spProjYp(hsm,rv);
spReconYp(recon,hsm); spReconYp(recon,hsm);
full = rv + Gamma(Gamma::GammaY) *rv; full = rv + Gamma(Gamma::GammaY) *rv;
mag = TensorRemove(norm2l(full-recon)); mag = TensorRemove(norm2(full-recon));
std::cout << "Yp "<< mag<<std::endl; std::cout << "Yp "<< mag<<std::endl;
// Ym // Ym
spProjYm(hsm,rv); spProjYm(hsm,rv);
spReconYm(recon,hsm); spReconYm(recon,hsm);
full = rv - Gamma(Gamma::GammaY) *rv; full = rv - Gamma(Gamma::GammaY) *rv;
mag = TensorRemove(norm2l(full-recon)); mag = TensorRemove(norm2(full-recon));
std::cout << "Ym "<< mag<<std::endl; std::cout << "Ym "<< mag<<std::endl;
// Zp // Zp
spProjZp(hsm,rv); spProjZp(hsm,rv);
spReconZp(recon,hsm); spReconZp(recon,hsm);
full = rv + Gamma(Gamma::GammaZ) *rv; full = rv + Gamma(Gamma::GammaZ) *rv;
mag = TensorRemove(norm2l(full-recon)); mag = TensorRemove(norm2(full-recon));
std::cout << "Zp "<< mag<<std::endl; std::cout << "Zp "<< mag<<std::endl;
// Zm // Zm
spProjZm(hsm,rv); spProjZm(hsm,rv);
spReconZm(recon,hsm); spReconZm(recon,hsm);
full = rv - Gamma(Gamma::GammaZ) *rv; full = rv - Gamma(Gamma::GammaZ) *rv;
mag = TensorRemove(norm2l(full-recon)); mag = TensorRemove(norm2(full-recon));
std::cout << "Zm "<< mag<<std::endl; std::cout << "Zm "<< mag<<std::endl;
// Tp // Tp
spProjTp(hsm,rv); spProjTp(hsm,rv);
spReconTp(recon,hsm); spReconTp(recon,hsm);
full = rv + Gamma(Gamma::GammaT) *rv; full = rv + Gamma(Gamma::GammaT) *rv;
mag = TensorRemove(norm2l(full-recon)); mag = TensorRemove(norm2(full-recon));
std::cout << "Tp "<< mag<<std::endl; std::cout << "Tp "<< mag<<std::endl;
// Tm // Tm
spProjTm(hsm,rv); spProjTm(hsm,rv);
spReconTm(recon,hsm); spReconTm(recon,hsm);
full = rv - Gamma(Gamma::GammaT) *rv; full = rv - Gamma(Gamma::GammaT) *rv;
mag = TensorRemove(norm2l(full-recon)); mag = TensorRemove(norm2(full-recon));
std::cout << "Tm "<< mag<<std::endl; std::cout << "Tm "<< mag<<std::endl;
// 5p // 5p
spProj5p(hsm,rv); spProj5p(hsm,rv);
spRecon5p(recon,hsm); spRecon5p(recon,hsm);
full = rv + Gamma(Gamma::Gamma5) *rv; full = rv + Gamma(Gamma::Gamma5) *rv;
mag = TensorRemove(norm2l(full-recon)); mag = TensorRemove(norm2(full-recon));
std::cout << "5p "<< mag<<std::endl; std::cout << "5p "<< mag<<std::endl;
// 5m // 5m
spProj5m(hsm,rv); spProj5m(hsm,rv);
spRecon5m(recon,hsm); spRecon5m(recon,hsm);
full = rv - Gamma(Gamma::Gamma5) *rv; full = rv - Gamma(Gamma::Gamma5) *rv;
mag = TensorRemove(norm2l(full-recon)); mag = TensorRemove(norm2(full-recon));
std::cout << "5m "<< mag<<std::endl; std::cout << "5m "<< mag<<std::endl;
Grid_finalize(); Grid_finalize();