portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-11-04 14:04:32 +00:00
Code Releases Activity

Moving some things around for pretty

Browse Source
This commit is contained in:
Peter Boyle
2015-05-11 19:09:49 +01:00
parent 22d384b07d
commit b1d2c60d07
6 changed files with 336 additions and 324 deletions
Download Patch File Download Diff File Expand all files Collapse all files

113
lib/Grid_communicator.h
View File

@@ -1,117 +1,6 @@
#ifndef GRID_COMMUNICATOR_H
#define GRID_COMMUNICATOR_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);
};
}
#include <communicator/Grid_communicator_base.h>
#endif

199
lib/Grid_lattice.h
View File

@@ -1,203 +1,6 @@
#ifndef GRID_LATTICE_H
#define GRID_LATTICE_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>
#include <lattice/Grid_lattice_base.h>
#endif

117
lib/communicator/Grid_communicator_base.h Normal file
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

203
lib/lattice/Grid_lattice_base.h Normal file
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

2
lib/math/Grid_math_inner.h
View File

@@ -6,7 +6,7 @@ namespace Grid {
// innerProduct Vector x Vector -> 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;
decltype(innerProduct(arg,arg)) nrm;
nrm = innerProduct(arg,arg);

26
tests/Grid_gamma.cc
View File

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