mirror of
https://github.com/paboyle/Grid.git
synced 2024-11-10 07:55:35 +00:00
Got the NERSC IO working and fixed a bug in cshift.
This commit is contained in:
Peter Boyle
parent
42f167ea37
commit
b32c14b433
22
configure
vendored
22
configure
vendored
@ -4999,6 +4999,28 @@ _ACEOF
|
||||
fi
|
||||
done
|
||||
|
||||
for ac_func in ntohll
|
||||
do :
|
||||
ac_fn_c_check_func "$LINENO" "ntohll" "ac_cv_func_ntohll"
|
||||
if test "x$ac_cv_func_ntohll" = xyes; then :
|
||||
cat >>confdefs.h <<_ACEOF
|
||||
#define HAVE_NTOHLL 1
|
||||
_ACEOF
|
||||
|
||||
fi
|
||||
done
|
||||
|
||||
for ac_func in be64toh
|
||||
do :
|
||||
ac_fn_c_check_func "$LINENO" "be64toh" "ac_cv_func_be64toh"
|
||||
if test "x$ac_cv_func_be64toh" = xyes; then :
|
||||
cat >>confdefs.h <<_ACEOF
|
||||
#define HAVE_BE64TOH 1
|
||||
_ACEOF
|
||||
|
||||
fi
|
||||
done
|
||||
|
||||
|
||||
# Check whether --enable-simd was given.
|
||||
if test "${enable_simd+set}" = set; then :
|
||||
|
||||
@ -25,6 +25,8 @@ AC_TYPE_UINT64_T
|
||||
|
||||
# Checks for library functions.
|
||||
AC_CHECK_FUNCS([gettimeofday])
|
||||
AC_CHECK_FUNCS([ntohll])
|
||||
AC_CHECK_FUNCS([be64toh])
|
||||
|
||||
AC_ARG_ENABLE([simd],[AC_HELP_STRING([--enable-simd=SSE|AVX|AVX2|AVX512],[Select instructions])],[ac_SIMD=${enable_simd}],[ac_SIMD=AVX2])
|
||||
|
||||
|
||||
@ -51,6 +51,7 @@
|
||||
#include <Grid_where.h>
|
||||
#include <Grid_stencil.h>
|
||||
#include <qcd/Grid_QCD.h>
|
||||
#include <parallelIO/GridNerscIO.h>
|
||||
|
||||
namespace Grid {
|
||||
|
||||
|
||||
@ -35,6 +35,7 @@ class CartesianCommunicator {
|
||||
// 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; };
|
||||
@ -49,6 +50,8 @@ class CartesianCommunicator {
|
||||
void GlobalSum(RealD &);
|
||||
void GlobalSumVector(RealD *,int N);
|
||||
|
||||
void GlobalSum(uint32_t &);
|
||||
|
||||
void GlobalSum(ComplexF &c)
|
||||
{
|
||||
GlobalSumVector((float *)&c,2);
|
||||
@ -68,12 +71,10 @@ class CartesianCommunicator {
|
||||
}
|
||||
|
||||
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;
|
||||
GlobalSum(ptr,words);
|
||||
GlobalSumVector(ptr,words);
|
||||
}
|
||||
////////////////////////////////////////////////////////////
|
||||
// Face exchange
|
||||
|
||||
@ -8,42 +8,42 @@ namespace Grid {
|
||||
public:
|
||||
vInteger operator()(const lobj &lhs, const robj &rhs)
|
||||
{
|
||||
return lhs == rhs;
|
||||
return TensorRemove(lhs) == TensorRemove(rhs);
|
||||
}
|
||||
};
|
||||
template<class lobj,class robj> class vne {
|
||||
public:
|
||||
vInteger operator()(const lobj &lhs, const robj &rhs)
|
||||
{
|
||||
return lhs != rhs;
|
||||
return TensorRemove(lhs) != TensorRemove(rhs);
|
||||
}
|
||||
};
|
||||
template<class lobj,class robj> class vlt {
|
||||
public:
|
||||
vInteger operator()(const lobj &lhs, const robj &rhs)
|
||||
{
|
||||
return lhs < rhs;
|
||||
return TensorRemove(lhs) < TensorRemove(rhs);
|
||||
}
|
||||
};
|
||||
template<class lobj,class robj> class vle {
|
||||
public:
|
||||
vInteger operator()(const lobj &lhs, const robj &rhs)
|
||||
{
|
||||
return lhs <= rhs;
|
||||
return TensorRemove(lhs) <= TensorRemove(rhs);
|
||||
}
|
||||
};
|
||||
template<class lobj,class robj> class vgt {
|
||||
public:
|
||||
vInteger operator()(const lobj &lhs, const robj &rhs)
|
||||
{
|
||||
return lhs > rhs;
|
||||
return TensorRemove(lhs) > TensorRemove(rhs);
|
||||
}
|
||||
};
|
||||
template<class lobj,class robj> class vge {
|
||||
public:
|
||||
vInteger operator()(const lobj &lhs, const robj &rhs)
|
||||
{
|
||||
return lhs >= rhs;
|
||||
return TensorRemove(lhs) >= TensorRemove(rhs);
|
||||
}
|
||||
};
|
||||
|
||||
@ -52,42 +52,42 @@ namespace Grid {
|
||||
public:
|
||||
Integer operator()(const lobj &lhs, const robj &rhs)
|
||||
{
|
||||
return lhs == rhs;
|
||||
return TensorRemove(lhs) == TensorRemove(rhs);
|
||||
}
|
||||
};
|
||||
template<class lobj,class robj> class sne {
|
||||
public:
|
||||
Integer operator()(const lobj &lhs, const robj &rhs)
|
||||
{
|
||||
return lhs != rhs;
|
||||
return TensorRemove(lhs) != TensorRemove(rhs);
|
||||
}
|
||||
};
|
||||
template<class lobj,class robj> class slt {
|
||||
public:
|
||||
Integer operator()(const lobj &lhs, const robj &rhs)
|
||||
{
|
||||
return lhs < rhs;
|
||||
return TensorRemove(lhs) < TensorRemove(rhs);
|
||||
}
|
||||
};
|
||||
template<class lobj,class robj> class sle {
|
||||
public:
|
||||
Integer operator()(const lobj &lhs, const robj &rhs)
|
||||
{
|
||||
return lhs <= rhs;
|
||||
return TensorRemove(lhs) <= TensorRemove(rhs);
|
||||
}
|
||||
};
|
||||
template<class lobj,class robj> class sgt {
|
||||
public:
|
||||
Integer operator()(const lobj &lhs, const robj &rhs)
|
||||
{
|
||||
return lhs > rhs;
|
||||
return TensorRemove(lhs) > TensorRemove(rhs);
|
||||
}
|
||||
};
|
||||
template<class lobj,class robj> class sge {
|
||||
public:
|
||||
Integer operator()(const lobj &lhs, const robj &rhs)
|
||||
{
|
||||
return lhs >= rhs;
|
||||
return TensorRemove(lhs) >= TensorRemove(rhs);
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
@ -16,6 +16,9 @@
|
||||
/* GRID_COMMS_NONE */
|
||||
/* #undef GRID_COMMS_NONE */
|
||||
|
||||
/* Define to 1 if you have the `be64toh' function. */
|
||||
/* #undef HAVE_BE64TOH */
|
||||
|
||||
/* Define to 1 if you have the `gettimeofday' function. */
|
||||
#define HAVE_GETTIMEOFDAY 1
|
||||
|
||||
@ -31,6 +34,9 @@
|
||||
/* Define to 1 if you have the <memory.h> header file. */
|
||||
#define HAVE_MEMORY_H 1
|
||||
|
||||
/* Define to 1 if you have the `ntohll' function. */
|
||||
/* #undef HAVE_NTOHLL */
|
||||
|
||||
/* Define to 1 if you have the <stdint.h> header file. */
|
||||
#define HAVE_STDINT_H 1
|
||||
|
||||
|
||||
@ -15,6 +15,9 @@
|
||||
/* GRID_COMMS_NONE */
|
||||
#undef GRID_COMMS_NONE
|
||||
|
||||
/* Define to 1 if you have the `be64toh' function. */
|
||||
#undef HAVE_BE64TOH
|
||||
|
||||
/* Define to 1 if you have the `gettimeofday' function. */
|
||||
#undef HAVE_GETTIMEOFDAY
|
||||
|
||||
@ -30,6 +33,9 @@
|
||||
/* Define to 1 if you have the <memory.h> header file. */
|
||||
#undef HAVE_MEMORY_H
|
||||
|
||||
/* Define to 1 if you have the `ntohll' function. */
|
||||
#undef HAVE_NTOHLL
|
||||
|
||||
/* Define to 1 if you have the <stdint.h> header file. */
|
||||
#undef HAVE_STDINT_H
|
||||
|
||||
|
||||
@ -37,7 +37,8 @@ public:
|
||||
// what about a default grid?
|
||||
//////////////////////////////////////////////////////////////////
|
||||
Lattice(GridBase *grid) : _grid(grid) {
|
||||
_odata.reserve(_grid->oSites());
|
||||
// _odata.reserve(_grid->oSites());
|
||||
_odata.resize(_grid->oSites());
|
||||
assert((((uint64_t)&_odata[0])&0xF) ==0);
|
||||
checkerboard=0;
|
||||
}
|
||||
@ -49,6 +50,15 @@ public:
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
template<class robj> inline Lattice<vobj> & operator = (const Lattice<robj> & r){
|
||||
conformable(*this,r);
|
||||
std::cout<<"Lattice operator ="<<std::endl;
|
||||
|
||||
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) {
|
||||
|
||||
@ -6,34 +6,23 @@
|
||||
|
||||
namespace Grid{
|
||||
|
||||
class GridRNG ; // Forward declaration;
|
||||
|
||||
//////////////////////////////////////////////////////////////////////
|
||||
// Commicator provides information on the processor grid
|
||||
//////////////////////////////////////////////////////////////////////
|
||||
// unsigned long _ndimension;
|
||||
// std::vector<int> _processors; // processor grid
|
||||
// int _processor; // linear processor rank
|
||||
// std::vector<int> _processor_coor; // linear processor rank
|
||||
//////////////////////////////////////////////////////////////////////
|
||||
class GridBase : public CartesianCommunicator {
|
||||
|
||||
public:
|
||||
|
||||
// Give Lattice access
|
||||
template<class object> friend class Lattice;
|
||||
|
||||
GridRNG *_rng;
|
||||
|
||||
GridBase(std::vector<int> & processor_grid) : CartesianCommunicator(processor_grid) {};
|
||||
|
||||
//FIXME
|
||||
// protected:
|
||||
// Lattice wide random support. not yet fully implemented. Need seed strategy
|
||||
// and one generator per site.
|
||||
// std::default_random_engine generator;
|
||||
// static std::mt19937 generator( 9 );
|
||||
|
||||
//////////////////////////////////////////////////////////////////////
|
||||
// Commicator provides information on the processor grid
|
||||
//////////////////////////////////////////////////////////////////////
|
||||
// unsigned long _ndimension;
|
||||
// std::vector<int> _processors; // processor grid
|
||||
// int _processor; // linear processor rank
|
||||
// std::vector<int> _processor_coor; // linear processor rank
|
||||
//////////////////////////////////////////////////////////////////////
|
||||
|
||||
// Physics Grid information.
|
||||
std::vector<int> _simd_layout;// Which dimensions get relayed out over simd lanes.
|
||||
std::vector<int> _fdimensions;// Global dimensions of array prior to cb removal
|
||||
@ -101,6 +90,13 @@ public:
|
||||
for(int d=0;d<_ndimension;d++) idx+=_ostride[d]*ocoor[d];
|
||||
return idx;
|
||||
}
|
||||
inline void CoorFromIndex (std::vector<int>& coor,int index,std::vector<int> &dims){
|
||||
coor.resize(_ndimension);
|
||||
for(int d=0;d<_ndimension;d++){
|
||||
coor[d] = index % dims[d];
|
||||
index = index / dims[d];
|
||||
}
|
||||
}
|
||||
inline void oCoorFromOindex (std::vector<int>& coor,int Oindex){
|
||||
coor.resize(_ndimension);
|
||||
for(int d=0;d<_ndimension;d++){
|
||||
@ -146,6 +142,7 @@ public:
|
||||
inline int lSites(void) { return _isites*_osites; };
|
||||
inline int gSites(void) { return _isites*_osites*_Nprocessors; };
|
||||
inline int Nd (void) { return _ndimension;};
|
||||
|
||||
inline const std::vector<int> &FullDimensions(void) { return _fdimensions;};
|
||||
inline const std::vector<int> &GlobalDimensions(void) { return _gdimensions;};
|
||||
inline const std::vector<int> &LocalDimensions(void) { return _ldimensions;};
|
||||
@ -175,10 +172,10 @@ public:
|
||||
std::vector<int> coor(_ndimension);
|
||||
|
||||
ProcessorCoorFromRank(rank,coor);
|
||||
for(int mu=0;mu<_ndimension;mu++) gcoor[mu] = _ldimensions[mu]&coor[mu];
|
||||
for(int mu=0;mu<_ndimension;mu++) gcoor[mu] = _ldimensions[mu]*coor[mu];
|
||||
|
||||
iCoorFromIindex(coor,i_idx);
|
||||
for(int mu=0;mu<_ndimension;mu++) gcoor[mu] += _rdimensions[mu]&coor[mu];
|
||||
for(int mu=0;mu<_ndimension;mu++) gcoor[mu] += _rdimensions[mu]*coor[mu];
|
||||
|
||||
oCoorFromOindex (coor,o_idx);
|
||||
for(int mu=0;mu<_ndimension;mu++) gcoor[mu] += coor[mu];
|
||||
|
||||
@ -17,6 +17,7 @@ CartesianCommunicator::CartesianCommunicator(std::vector<int> &processors)
|
||||
void CartesianCommunicator::GlobalSum(float &){}
|
||||
void CartesianCommunicator::GlobalSumVector(float *,int N){}
|
||||
void CartesianCommunicator::GlobalSum(double &){}
|
||||
void CartesianCommunicator::GlobalSum(uint32_t &){}
|
||||
void CartesianCommunicator::GlobalSumVector(double *,int N){}
|
||||
|
||||
// Basic Halo comms primitive
|
||||
|
||||
@ -28,6 +28,9 @@ CartesianCommunicator::CartesianCommunicator(std::vector<int> &processors)
|
||||
assert(Size==_Nprocessors);
|
||||
}
|
||||
|
||||
void CartesianCommunicator::GlobalSum(uint32_t &u){
|
||||
MPI_Allreduce(MPI_IN_PLACE,&u,1,MPI_UINT32_T,MPI_SUM,communicator);
|
||||
}
|
||||
void CartesianCommunicator::GlobalSum(float &f){
|
||||
MPI_Allreduce(MPI_IN_PLACE,&f,1,MPI_FLOAT,MPI_SUM,communicator);
|
||||
}
|
||||
|
||||
@ -188,8 +188,8 @@ template<class vobj> void Cshift_comms_simd(Lattice<vobj> &ret,Lattice<vobj> &r
|
||||
|
||||
if ( comm_any ) {
|
||||
|
||||
for(int i=0;i<Nsimd;i++){
|
||||
pointers[i] = (scalar_type *)&send_buf_extract[i][0];
|
||||
for(int i=0;i<Nsimd;i++){ // there is a reversal in extract merge
|
||||
pointers[i] = (scalar_type *)&send_buf_extract[Nsimd-1-i][0];
|
||||
}
|
||||
Gather_plane_extract(rhs,pointers,dimension,sx,cbmask);
|
||||
|
||||
@ -238,9 +238,9 @@ template<class vobj> void Cshift_comms_simd(Lattice<vobj> &ret,Lattice<vobj> &r
|
||||
for(int i=0;i<Nsimd;i++){
|
||||
if ( permute_slice ) {
|
||||
PermuteMap=i^toggle_bit;
|
||||
pointers[i] = rpointers[PermuteMap];
|
||||
pointers[Nsimd-1-i] = rpointers[PermuteMap];
|
||||
} else {
|
||||
pointers[i] = rpointers[i];
|
||||
pointers[Nsimd-1-i] = rpointers[i];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -5,16 +5,17 @@ namespace Grid {
|
||||
|
||||
template<class iobj> inline void LatticeCoordinate(Lattice<iobj> &l,int mu)
|
||||
{
|
||||
typedef typename iobj::scalar_type scalar_type;
|
||||
typedef typename iobj::vector_type vector_type;
|
||||
GridBase *grid = l._grid;
|
||||
int Nsimd = grid->iSites();
|
||||
std::vector<int> gcoor;
|
||||
std::vector<Integer> mergebuf(Nsimd);
|
||||
std::vector<Integer *> mergeptr(Nsimd);
|
||||
vInteger vI;
|
||||
std::vector<scalar_type> mergebuf(Nsimd);
|
||||
std::vector<scalar_type *> mergeptr(Nsimd);
|
||||
vector_type vI;
|
||||
for(int o=0;o<grid->oSites();o++){
|
||||
for(int i=0;i<grid->iSites();i++){
|
||||
grid->RankIndexToGlobalCoor(grid->ThisRank(),o,i,gcoor);
|
||||
// grid->RankIndexToGlobalCoor(0,o,i,gcoor);
|
||||
mergebuf[i]=gcoor[mu];
|
||||
mergeptr[i]=&mergebuf[i];
|
||||
}
|
||||
|
||||
@ -43,6 +43,7 @@ namespace Grid {
|
||||
}
|
||||
return ret;
|
||||
};
|
||||
|
||||
////////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
// Poke internal indices of a Lattice object
|
||||
////////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
@ -88,26 +89,26 @@ namespace Grid {
|
||||
assert( sizeof(sobj)*Nsimd == sizeof(vobj));
|
||||
|
||||
int rank,odx,idx;
|
||||
grid->GlobalCoorToRankIndex(rank,odx,idx,site);
|
||||
|
||||
// Optional to broadcast from node 0.
|
||||
grid->Broadcast(0,s);
|
||||
grid->GlobalCoorToRankIndex(rank,odx,idx,site);
|
||||
grid->Broadcast(grid->BossRank(),s);
|
||||
|
||||
std::vector<sobj> buf(Nsimd);
|
||||
std::vector<scalar_type *> pointers(Nsimd);
|
||||
for(int i=0;i<Nsimd;i++) pointers[i] = (scalar_type *)&buf[i];
|
||||
|
||||
// extract-modify-merge cycle is easiest way and this is not perf critical
|
||||
extract(l._odata[odx],pointers);
|
||||
|
||||
buf[idx] = s;
|
||||
|
||||
for(int i=0;i<Nsimd;i++) pointers[i] = (scalar_type *)&buf[i];
|
||||
merge(l._odata[odx],pointers);
|
||||
if ( rank == grid->ThisRank() ) {
|
||||
for(int i=0;i<Nsimd;i++) pointers[i] = (scalar_type *)&buf[i];
|
||||
extract(l._odata[odx],pointers);
|
||||
buf[idx] = s;
|
||||
for(int i=0;i<Nsimd;i++) pointers[i] = (scalar_type *)&buf[i];
|
||||
merge(l._odata[odx],pointers);
|
||||
}
|
||||
|
||||
return;
|
||||
};
|
||||
|
||||
|
||||
//////////////////////////////////////////////////////////
|
||||
// Peek a scalar object from the SIMD array
|
||||
//////////////////////////////////////////////////////////
|
||||
@ -138,6 +139,69 @@ namespace Grid {
|
||||
return;
|
||||
};
|
||||
|
||||
//////////////////////////////////////////////////////////
|
||||
// Peek a scalar object from the SIMD array
|
||||
//////////////////////////////////////////////////////////
|
||||
template<class vobj,class sobj>
|
||||
void peekLocalSite(sobj &s,Lattice<vobj> &l,std::vector<int> &site){
|
||||
|
||||
GridBase *grid=l._grid;
|
||||
|
||||
typedef typename vobj::scalar_type scalar_type;
|
||||
typedef typename vobj::vector_type vector_type;
|
||||
|
||||
int Nsimd = grid->Nsimd();
|
||||
|
||||
assert( l.checkerboard== l._grid->CheckerBoard(site));
|
||||
assert( sizeof(sobj)*Nsimd == sizeof(vobj));
|
||||
|
||||
int odx,idx;
|
||||
idx= grid->iIndex(site);
|
||||
odx= grid->oIndex(site);
|
||||
|
||||
std::vector<sobj> buf(Nsimd);
|
||||
std::vector<scalar_type *> pointers(Nsimd);
|
||||
for(int i=0;i<Nsimd;i++) pointers[i] = (scalar_type *)&buf[i];
|
||||
|
||||
extract(l._odata[odx],pointers);
|
||||
|
||||
s = buf[idx];
|
||||
|
||||
return;
|
||||
};
|
||||
|
||||
template<class vobj,class sobj>
|
||||
void pokeLocalSite(const sobj &s,Lattice<vobj> &l,std::vector<int> &site){
|
||||
|
||||
GridBase *grid=l._grid;
|
||||
|
||||
typedef typename vobj::scalar_type scalar_type;
|
||||
typedef typename vobj::vector_type vector_type;
|
||||
|
||||
int Nsimd = grid->Nsimd();
|
||||
|
||||
assert( l.checkerboard== l._grid->CheckerBoard(site));
|
||||
assert( sizeof(sobj)*Nsimd == sizeof(vobj));
|
||||
|
||||
int odx,idx;
|
||||
idx= grid->iIndex(site);
|
||||
odx= grid->oIndex(site);
|
||||
|
||||
std::vector<sobj> buf(Nsimd);
|
||||
std::vector<scalar_type *> pointers(Nsimd);
|
||||
for(int i=0;i<Nsimd;i++) pointers[i] = (scalar_type *)&buf[i];
|
||||
|
||||
// extract-modify-merge cycle is easiest way and this is not perf critical
|
||||
extract(l._odata[odx],pointers);
|
||||
|
||||
buf[idx] = s;
|
||||
|
||||
for(int i=0;i<Nsimd;i++) pointers[i] = (scalar_type *)&buf[i];
|
||||
merge(l._odata[odx],pointers);
|
||||
|
||||
return;
|
||||
};
|
||||
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
@ -40,6 +40,36 @@ namespace Grid {
|
||||
return nrm;
|
||||
}
|
||||
|
||||
template<class vobj>
|
||||
inline typename vobj::scalar_object sum(const Lattice<vobj> &arg){
|
||||
|
||||
GridBase *grid=arg._grid;
|
||||
int Nsimd = grid->Nsimd();
|
||||
|
||||
typedef typename vobj::scalar_object sobj;
|
||||
typedef typename vobj::scalar_type scalar_type;
|
||||
|
||||
vobj vsum;
|
||||
sobj ssum;
|
||||
|
||||
vsum=zero;
|
||||
ssum=zero;
|
||||
for(int ss=0;ss<arg._grid->oSites(); ss++){
|
||||
vsum = vsum + arg._odata[ss];
|
||||
}
|
||||
|
||||
std::vector<sobj> buf(Nsimd);
|
||||
std::vector<scalar_type *> pointers(Nsimd);
|
||||
for(int i=0;i<Nsimd;i++) pointers[i] = (scalar_type *)&buf[i];
|
||||
extract(vsum,pointers);
|
||||
|
||||
for(int i=0;i<Nsimd;i++) ssum = ssum + buf[i];
|
||||
|
||||
arg._grid->GlobalSum(ssum);
|
||||
|
||||
return ssum;
|
||||
}
|
||||
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
@ -13,10 +13,12 @@ template<class vtype> class iScalar
|
||||
public:
|
||||
vtype _internal;
|
||||
|
||||
typedef typename GridTypeMapper<vtype>::scalar_type scalar_type;
|
||||
typedef typename GridTypeMapper<vtype>::scalar_type scalar_type;
|
||||
typedef typename GridTypeMapper<vtype>::vector_type vector_type;
|
||||
typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v;
|
||||
typedef iScalar<tensor_reduced_v> tensor_reduced;
|
||||
typedef typename GridTypeMapper<vtype>::scalar_object recurse_scalar_object;
|
||||
typedef iScalar<recurse_scalar_object> scalar_object;
|
||||
|
||||
enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1};
|
||||
|
||||
@ -27,7 +29,7 @@ public:
|
||||
|
||||
iScalar(scalar_type s) : _internal(s) {};// recurse down and hit the constructor for vector_type
|
||||
|
||||
iScalar(Zero &z){ *this = zero; };
|
||||
iScalar(const Zero &z){ *this = zero; };
|
||||
|
||||
iScalar<vtype> & operator= (const Zero &hero){
|
||||
zeroit(*this);
|
||||
@ -73,11 +75,18 @@ public:
|
||||
operator ComplexD () const { return(TensorRemove(_internal)); };
|
||||
operator RealD () const { return(real(TensorRemove(_internal))); }
|
||||
|
||||
|
||||
template<class T,typename std::enable_if<isGridTensor<T>::notvalue, T>::type* = nullptr > inline auto operator = (T arg) -> iScalar<vtype>
|
||||
{
|
||||
_internal = arg;
|
||||
return *this;
|
||||
}
|
||||
|
||||
};
|
||||
///////////////////////////////////////////////////////////
|
||||
// Allows to turn scalar<scalar<scalar<double>>>> back to double.
|
||||
///////////////////////////////////////////////////////////
|
||||
template<class T> inline typename std::enable_if<isGridTensor<T>::notvalue, T>::type TensorRemove(T arg) { return arg;}
|
||||
template<class T> inline typename std::enable_if<isGridTensor<T>::notvalue, T>::type TensorRemove(T arg) { return arg;}
|
||||
template<class vtype> inline auto TensorRemove(iScalar<vtype> arg) -> decltype(TensorRemove(arg._internal))
|
||||
{
|
||||
return TensorRemove(arg._internal);
|
||||
@ -91,14 +100,17 @@ public:
|
||||
typedef typename GridTypeMapper<vtype>::scalar_type scalar_type;
|
||||
typedef typename GridTypeMapper<vtype>::vector_type vector_type;
|
||||
typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v;
|
||||
typedef typename GridTypeMapper<vtype>::scalar_object recurse_scalar_object;
|
||||
typedef iScalar<tensor_reduced_v> tensor_reduced;
|
||||
typedef iVector<recurse_scalar_object,N> scalar_object;
|
||||
|
||||
|
||||
enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1};
|
||||
typedef iScalar<tensor_reduced_v> tensor_reduced;
|
||||
|
||||
iVector(Zero &z){ *this = zero; };
|
||||
iVector(const Zero &z){ *this = zero; };
|
||||
iVector() {};// Empty constructure
|
||||
|
||||
iVector<vtype,N> & operator= (Zero &hero){
|
||||
iVector<vtype,N> & operator= (const Zero &hero){
|
||||
zeroit(*this);
|
||||
return *this;
|
||||
}
|
||||
@ -153,18 +165,27 @@ public:
|
||||
|
||||
typedef typename GridTypeMapper<vtype>::scalar_type scalar_type;
|
||||
typedef typename GridTypeMapper<vtype>::vector_type vector_type;
|
||||
|
||||
typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v;
|
||||
typedef typename GridTypeMapper<vtype>::scalar_object recurse_scalar_object;
|
||||
typedef iScalar<tensor_reduced_v> tensor_reduced;
|
||||
typedef iMatrix<recurse_scalar_object,N> scalar_object;
|
||||
|
||||
enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1};
|
||||
typedef iScalar<tensor_reduced_v> tensor_reduced;
|
||||
|
||||
iMatrix(Zero &z){ *this = zero; };
|
||||
iMatrix(const Zero &z){ *this = zero; };
|
||||
iMatrix() {};
|
||||
iMatrix<vtype,N> & operator= (Zero &hero){
|
||||
iMatrix<vtype,N> & operator= (const Zero &hero){
|
||||
zeroit(*this);
|
||||
return *this;
|
||||
}
|
||||
template<class T,typename std::enable_if<isGridTensor<T>::notvalue, T>::type* = nullptr > inline auto operator = (T arg) -> iMatrix<vtype,N>
|
||||
{
|
||||
zeroit(*this);
|
||||
for(int i=0;i<N;i++)
|
||||
_internal[i][i] = arg;
|
||||
return *this;
|
||||
}
|
||||
|
||||
friend void zeroit(iMatrix<vtype,N> &that){
|
||||
for(int i=0;i<N;i++){
|
||||
for(int j=0;j<N;j++){
|
||||
|
||||
@ -25,6 +25,7 @@ namespace Grid {
|
||||
typedef typename T::scalar_type scalar_type;
|
||||
typedef typename T::vector_type vector_type;
|
||||
typedef typename T::tensor_reduced tensor_reduced;
|
||||
typedef typename T::scalar_object scalar_object;
|
||||
enum { TensorLevel = T::TensorLevel };
|
||||
};
|
||||
|
||||
@ -36,6 +37,7 @@ namespace Grid {
|
||||
typedef RealF scalar_type;
|
||||
typedef RealF vector_type;
|
||||
typedef RealF tensor_reduced ;
|
||||
typedef RealF scalar_object;
|
||||
enum { TensorLevel = 0 };
|
||||
};
|
||||
template<> class GridTypeMapper<RealD> {
|
||||
@ -43,6 +45,7 @@ namespace Grid {
|
||||
typedef RealD scalar_type;
|
||||
typedef RealD vector_type;
|
||||
typedef RealD tensor_reduced;
|
||||
typedef RealD scalar_object;
|
||||
enum { TensorLevel = 0 };
|
||||
};
|
||||
template<> class GridTypeMapper<ComplexF> {
|
||||
@ -50,6 +53,7 @@ namespace Grid {
|
||||
typedef ComplexF scalar_type;
|
||||
typedef ComplexF vector_type;
|
||||
typedef ComplexF tensor_reduced;
|
||||
typedef ComplexF scalar_object;
|
||||
enum { TensorLevel = 0 };
|
||||
};
|
||||
template<> class GridTypeMapper<ComplexD> {
|
||||
@ -57,6 +61,7 @@ namespace Grid {
|
||||
typedef ComplexD scalar_type;
|
||||
typedef ComplexD vector_type;
|
||||
typedef ComplexD tensor_reduced;
|
||||
typedef ComplexD scalar_object;
|
||||
enum { TensorLevel = 0 };
|
||||
};
|
||||
|
||||
@ -65,6 +70,7 @@ namespace Grid {
|
||||
typedef RealF scalar_type;
|
||||
typedef vRealF vector_type;
|
||||
typedef vRealF tensor_reduced;
|
||||
typedef RealF scalar_object;
|
||||
enum { TensorLevel = 0 };
|
||||
};
|
||||
template<> class GridTypeMapper<vRealD> {
|
||||
@ -72,6 +78,7 @@ namespace Grid {
|
||||
typedef RealD scalar_type;
|
||||
typedef vRealD vector_type;
|
||||
typedef vRealD tensor_reduced;
|
||||
typedef RealD scalar_object;
|
||||
enum { TensorLevel = 0 };
|
||||
};
|
||||
template<> class GridTypeMapper<vComplexF> {
|
||||
@ -79,6 +86,7 @@ namespace Grid {
|
||||
typedef ComplexF scalar_type;
|
||||
typedef vComplexF vector_type;
|
||||
typedef vComplexF tensor_reduced;
|
||||
typedef ComplexF scalar_object;
|
||||
enum { TensorLevel = 0 };
|
||||
};
|
||||
template<> class GridTypeMapper<vComplexD> {
|
||||
@ -86,6 +94,7 @@ namespace Grid {
|
||||
typedef ComplexD scalar_type;
|
||||
typedef vComplexD vector_type;
|
||||
typedef vComplexD tensor_reduced;
|
||||
typedef ComplexD scalar_object;
|
||||
enum { TensorLevel = 0 };
|
||||
};
|
||||
template<> class GridTypeMapper<vInteger> {
|
||||
@ -93,6 +102,7 @@ namespace Grid {
|
||||
typedef Integer scalar_type;
|
||||
typedef vInteger vector_type;
|
||||
typedef vInteger tensor_reduced;
|
||||
typedef Integer scalar_object;
|
||||
enum { TensorLevel = 0 };
|
||||
};
|
||||
|
||||
@ -102,6 +112,10 @@ namespace Grid {
|
||||
static const bool notvalue = false;
|
||||
};
|
||||
|
||||
template<> struct isGridTensor<int > {
|
||||
static const bool value = false;
|
||||
static const bool notvalue = true;
|
||||
};
|
||||
template<> struct isGridTensor<RealD > {
|
||||
static const bool value = false;
|
||||
static const bool notvalue = true;
|
||||
|
||||
524
lib/parallelIO/GridNerscIO.h
Normal file
524
lib/parallelIO/GridNerscIO.h
Normal file
@ -0,0 +1,524 @@
|
||||
#ifndef GRID_NERSC_IO_H
|
||||
#define GRID_NERSC_IO_H
|
||||
|
||||
#include <algorithm>
|
||||
#include <iostream>
|
||||
#include <iomanip>
|
||||
#include <fstream>
|
||||
|
||||
#include <arpa/inet.h>
|
||||
|
||||
namespace Grid {
|
||||
|
||||
using namespace QCD;
|
||||
|
||||
////////////////////////////////////////////////////////////////////////////////
|
||||
// Some data types for intermediate storage
|
||||
////////////////////////////////////////////////////////////////////////////////
|
||||
template<typename vtype> using iLorentzColour2x3 = iVector<iVector<iVector<vtype, Nc>, 2>, 4 >;
|
||||
typedef iLorentzColour2x3<Complex> LorentzColour2x3;
|
||||
typedef iLorentzColour2x3<ComplexF> LorentzColour2x3F;
|
||||
typedef iLorentzColour2x3<ComplexD> LorentzColour2x3D;
|
||||
|
||||
////////////////////////////////////////////////////////////////////////////////
|
||||
// header specification/interpretation
|
||||
////////////////////////////////////////////////////////////////////////////////
|
||||
class NerscField {
|
||||
public:
|
||||
// header strings (not in order)
|
||||
int dimension[4];
|
||||
std::string boundary[4];
|
||||
int data_start;
|
||||
std::string hdr_version;
|
||||
std::string storage_format;
|
||||
// Checks on data
|
||||
double link_trace;
|
||||
double plaquette;
|
||||
uint32_t checksum;
|
||||
unsigned int sequence_number;
|
||||
std::string data_type;
|
||||
std::string ensemble_id ;
|
||||
std::string ensemble_label ;
|
||||
std::string creator ;
|
||||
std::string creator_hardware ;
|
||||
std::string creation_date ;
|
||||
std::string archive_date ;
|
||||
std::string floating_point;
|
||||
};
|
||||
|
||||
|
||||
////////////////////////////////////////////////////////////////////////////////
|
||||
// Write and read from fstream; comput header offset for payload
|
||||
////////////////////////////////////////////////////////////////////////////////
|
||||
inline unsigned int writeNerscHeader(NerscField &field,std::string file)
|
||||
{
|
||||
std::ofstream fout(file,std::ios::out);
|
||||
|
||||
fout.seekp(0,std::ios::beg);
|
||||
fout << "BEGIN_HEADER" << std::endl;
|
||||
fout << "HDR_VERSION = " << field.hdr_version << std::endl;
|
||||
fout << "DATATYPE = " << field.data_type << std::endl;
|
||||
fout << "STORAGE_FORMAT = " << field.storage_format << std::endl;
|
||||
|
||||
for(int i=0;i<4;i++){
|
||||
fout << "DIMENSION_" << i+1 << " = " << field.dimension[i] << std::endl ;
|
||||
}
|
||||
// just to keep the space and write it later
|
||||
fout << "LINK_TRACE = " << std::setprecision(10) << field.link_trace << std::endl;
|
||||
fout << "PLAQUETTE = " << std::setprecision(10) << field.plaquette << std::endl;
|
||||
for(int i=0;i<4;i++){
|
||||
fout << "BOUNDARY_"<<i+1<<" = " << field.boundary[i] << std::endl;
|
||||
}
|
||||
fout << "CHECKSUM = "<< std::hex << std::setw(16) << 0 << field.checksum << std::endl;
|
||||
|
||||
fout << "ENSEMBLE_ID = " << field.ensemble_id << std::endl;
|
||||
fout << "ENSEMBLE_LABEL = " << field.ensemble_label << std::endl;
|
||||
fout << "SEQUENCE_NUMBER = " << field.sequence_number << std::endl;
|
||||
fout << "CREATOR = " << field.creator << std::endl;
|
||||
fout << "CREATOR_HARDWARE = "<< field.creator_hardware << std::endl;
|
||||
fout << "CREATION_DATE = " << field.creation_date << std::endl;
|
||||
fout << "ARCHIVE_DATE = " << field.archive_date << std::endl;
|
||||
fout << "FLOATING_POINT = " << field.floating_point << std::endl;
|
||||
fout << "END_HEADER" << std::endl;
|
||||
field.data_start = fout.tellp();
|
||||
return field.data_start;
|
||||
}
|
||||
|
||||
// A little helper
|
||||
inline void removeWhitespace(std::string &key)
|
||||
{
|
||||
key.erase(std::remove_if(key.begin(), key.end(), ::isspace),key.end());
|
||||
}
|
||||
// for the header-reader
|
||||
inline int readNerscHeader(std::string file,GridBase *grid, NerscField &field)
|
||||
{
|
||||
int offset=0;
|
||||
std::map<std::string,std::string> header;
|
||||
std::string line;
|
||||
|
||||
//////////////////////////////////////////////////
|
||||
// read the header
|
||||
//////////////////////////////////////////////////
|
||||
std::ifstream fin(file);
|
||||
|
||||
getline(fin,line); // read one line and insist is
|
||||
|
||||
removeWhitespace(line);
|
||||
assert(line==std::string("BEGIN_HEADER"));
|
||||
|
||||
do {
|
||||
getline(fin,line); // read one line
|
||||
int eq = line.find("=");
|
||||
if(eq >0) {
|
||||
std::string key=line.substr(0,eq);
|
||||
std::string val=line.substr(eq+1);
|
||||
removeWhitespace(key);
|
||||
removeWhitespace(val);
|
||||
|
||||
header[key] = val;
|
||||
}
|
||||
} while( line.find("END_HEADER") == std::string::npos );
|
||||
|
||||
field.data_start = fin.tellg();
|
||||
|
||||
//////////////////////////////////////////////////
|
||||
// chomp the values
|
||||
//////////////////////////////////////////////////
|
||||
|
||||
field.hdr_version = header[std::string("HDR_VERSION")];
|
||||
field.data_type = header[std::string("DATATYPE")];
|
||||
field.storage_format = header[std::string("STORAGE_FORMAT")];
|
||||
|
||||
field.dimension[0] = std::stol(header[std::string("DIMENSION_1")]);
|
||||
field.dimension[1] = std::stol(header[std::string("DIMENSION_2")]);
|
||||
field.dimension[2] = std::stol(header[std::string("DIMENSION_3")]);
|
||||
field.dimension[3] = std::stol(header[std::string("DIMENSION_4")]);
|
||||
|
||||
assert(grid->_ndimension == 4);
|
||||
for(int d=0;d<4;d++){
|
||||
assert(grid->_fdimensions[d]==field.dimension[d]);
|
||||
}
|
||||
|
||||
field.link_trace = std::stod(header[std::string("LINK_TRACE")]);
|
||||
field.plaquette = std::stod(header[std::string("PLAQUETTE")]);
|
||||
|
||||
field.boundary[0] = header[std::string("BOUNDARY_1")];
|
||||
field.boundary[1] = header[std::string("BOUNDARY_2")];
|
||||
field.boundary[2] = header[std::string("BOUNDARY_3")];
|
||||
field.boundary[3] = header[std::string("BOUNDARY_4")];
|
||||
|
||||
field.checksum = std::stoul(header[std::string("CHECKSUM")],0,16);
|
||||
field.ensemble_id = header[std::string("ENSEMBLE_ID")];
|
||||
field.ensemble_label = header[std::string("ENSEMBLE_LABEL")];
|
||||
field.sequence_number = std::stol(header[std::string("SEQUENCE_NUMBER")]);
|
||||
field.creator = header[std::string("CREATOR")];
|
||||
field.creator_hardware = header[std::string("CREATOR_HARDWARE")];
|
||||
field.creation_date = header[std::string("CREATION_DATE")];
|
||||
field.archive_date = header[std::string("ARCHIVE_DATE")];
|
||||
field.floating_point = header[std::string("FLOATING_POINT")];
|
||||
|
||||
return field.data_start;
|
||||
}
|
||||
|
||||
|
||||
//////////////////////////////////////////////////////////////////////
|
||||
// Utilities
|
||||
//////////////////////////////////////////////////////////////////////
|
||||
inline void reconstruct3(LorentzColourMatrix & cm)
|
||||
{
|
||||
const int x=0;
|
||||
const int y=1;
|
||||
const int z=2;
|
||||
for(int mu=0;mu<4;mu++){
|
||||
cm(mu)()(2,x) = adj(cm(mu)()(0,y)*cm(mu)()(1,z)-cm(mu)()(0,z)*cm(mu)()(1,y)); //x= yz-zy
|
||||
cm(mu)()(2,y) = adj(cm(mu)()(0,z)*cm(mu)()(1,x)-cm(mu)()(0,x)*cm(mu)()(1,z)); //y= zx-xz
|
||||
cm(mu)()(2,z) = adj(cm(mu)()(0,x)*cm(mu)()(1,y)-cm(mu)()(0,y)*cm(mu)()(1,x)); //z= xy-yx
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void inline be32toh(void *file_object,uint32_t bytes)
|
||||
{
|
||||
uint32_t * f = (uint32_t *)file_object;
|
||||
for(int i=0;i*sizeof(uint32_t)<bytes;i++){
|
||||
f[i] = ntohl(f[i]);
|
||||
}
|
||||
}
|
||||
void inline le32toh(void *file_object,uint32_t bytes)
|
||||
{
|
||||
uint32_t *fp = (uint32_t *)file_object;
|
||||
|
||||
uint32_t f;
|
||||
|
||||
for(int i=0;i*sizeof(uint32_t)<bytes;i++){
|
||||
f = fp[i];
|
||||
// got network order and the network to host
|
||||
f = ((f&0xFF)<<24) | ((f&0xFF00)<<8) | ((f&0xFF0000)>>8) | ((f&0xFF000000UL)>>24) ;
|
||||
fp[i] = ntohl(f);
|
||||
}
|
||||
}
|
||||
void inline be64toh(void *file_object,uint32_t bytes)
|
||||
{
|
||||
uint64_t * f = (uint64_t *)file_object;
|
||||
for(int i=0;i*sizeof(uint64_t)<bytes;i++){
|
||||
f[i] = ntohll(f[i]);
|
||||
}
|
||||
}
|
||||
void inline le64toh(void *file_object,uint32_t bytes)
|
||||
{
|
||||
uint64_t *fp = (uint64_t *)file_object;
|
||||
uint64_t f,g;
|
||||
|
||||
for(int i=0;i*sizeof(uint64_t)<bytes;i++){
|
||||
f = fp[i];
|
||||
// got network order and the network to host
|
||||
g = ((f&0xFF)<<24) | ((f&0xFF00)<<8) | ((f&0xFF0000)>>8) | ((f&0xFF000000UL)>>24) ;
|
||||
g = g << 32;
|
||||
f = f >> 32;
|
||||
g|= ((f&0xFF)<<24) | ((f&0xFF00)<<8) | ((f&0xFF0000)>>8) | ((f&0xFF000000UL)>>24) ;
|
||||
fp[i] = ntohl(g);
|
||||
}
|
||||
}
|
||||
|
||||
inline void NerscChecksum(uint32_t *buf,uint32_t buf_size,uint32_t &csum)
|
||||
{
|
||||
for(int i=0;i*sizeof(uint32_t)<buf_size;i++){
|
||||
csum=csum+buf[i];
|
||||
}
|
||||
}
|
||||
|
||||
template<class fobj,class sobj>
|
||||
struct NerscSimpleMunger{
|
||||
void operator() (fobj &in,sobj &out,uint32_t &csum){
|
||||
|
||||
for(int mu=0;mu<4;mu++){
|
||||
for(int i=0;i<3;i++){
|
||||
for(int j=0;j<3;j++){
|
||||
out(mu)()(i,j) = in(mu)()(i,j);
|
||||
}}}
|
||||
|
||||
NerscChecksum((uint32_t *)&in,sizeof(in),csum);
|
||||
};
|
||||
};
|
||||
template<class fobj,class sobj>
|
||||
struct NerscSimpleUnmunger{
|
||||
void operator() (sobj &in,fobj &out,uint32_t &csum){
|
||||
for(int mu=0;mu<4;mu++){
|
||||
for(int i=0;i<3;i++){
|
||||
for(int j=0;j<3;j++){
|
||||
out(mu)()(i,j) = in(mu)()(i,j);
|
||||
}}}
|
||||
NerscChecksum((uint32_t *)&out,sizeof(out),csum);
|
||||
};
|
||||
};
|
||||
|
||||
template<class fobj,class sobj>
|
||||
struct Nersc3x2munger{
|
||||
void operator() (fobj &in,sobj &out,uint32_t &csum){
|
||||
|
||||
NerscChecksum((uint32_t *)&in,sizeof(in),csum);
|
||||
|
||||
for(int mu=0;mu<4;mu++){
|
||||
for(int i=0;i<2;i++){
|
||||
for(int j=0;j<3;j++){
|
||||
out(mu)()(i,j) = in(mu)(i)(j);
|
||||
}}
|
||||
}
|
||||
reconstruct3(out);
|
||||
}
|
||||
};
|
||||
|
||||
template<class fobj,class sobj>
|
||||
struct Nersc3x2unmunger{
|
||||
|
||||
void operator() (sobj &in,fobj &out,uint32_t &csum){
|
||||
|
||||
NerscChecksum((uint32_t *)&out,sizeof(out),csum);
|
||||
|
||||
for(int mu=0;mu<4;mu++){
|
||||
for(int i=0;i<2;i++){
|
||||
for(int j=0;j<3;j++){
|
||||
out(mu)(i)(j) = in(mu)()(i,j);
|
||||
}}
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
////////////////////////////////////////////////////////////////////////////
|
||||
// Template wizardry to map types to strings for NERSC in an extensible way
|
||||
////////////////////////////////////////////////////////////////////////////
|
||||
template<class vobj> struct NerscDataType {
|
||||
static void DataType (std::string &str) { str = std::string("4D_BINARY_UNKNOWN"); };
|
||||
static void FloatingPoint(std::string &str) { str = std::string("IEEE64BIG"); };
|
||||
};
|
||||
|
||||
template<> struct NerscDataType<iColourMatrix<ComplexD> > {
|
||||
static void DataType (std::string &str) { str = std::string("4D_SU3_GAUGE_3X3"); };
|
||||
static void FloatingPoint(std::string &str) { str = std::string("IEEE64BIG");};
|
||||
};
|
||||
|
||||
template<> struct NerscDataType<iColourMatrix<ComplexF> > {
|
||||
static void DataType (std::string &str) { str = std::string("4D_SU3_GAUGE_3X3"); };
|
||||
static void FloatingPoint(std::string &str) { str = std::string("IEEE32BIG");};
|
||||
};
|
||||
|
||||
//////////////////////////////////////////////////////////////////////
|
||||
// Bit and Physical Checksumming and QA of data
|
||||
//////////////////////////////////////////////////////////////////////
|
||||
/*
|
||||
template<class vobj> inline uint32_t NerscChecksum(Lattice<vobj> & data)
|
||||
{
|
||||
uint32_t sum;
|
||||
for(int ss=0;ss<data._grid->Osites();ss++){
|
||||
uint32_t *iptr = (uint32_t *)& data._odata[0] ;
|
||||
for(int i=0;i<sizeof(vobj);i+=sizeof(uint32_t)){
|
||||
sum=sum+iptr[i];
|
||||
}
|
||||
}
|
||||
data._grid->globalSum(sum);
|
||||
return sum;
|
||||
}
|
||||
*/
|
||||
template<class vobj> inline void NerscPhysicalCharacteristics(Lattice<vobj> & data,NerscField &header)
|
||||
{
|
||||
header.data_type = NerscDataType<vobj>::DataType;
|
||||
header.floating_point = NerscDataType<vobj>::FloatingPoint;
|
||||
return;
|
||||
}
|
||||
|
||||
template<> inline void NerscPhysicalCharacteristics(LatticeGaugeField & data,NerscField &header)
|
||||
{
|
||||
NerscDataType<decltype(data._odata[0])>::DataType(header.data_type);
|
||||
NerscDataType<decltype(data._odata[0])>::FloatingPoint(header.floating_point);
|
||||
header.link_trace=1.0;
|
||||
header.plaquette =1.0;
|
||||
}
|
||||
|
||||
template<class vobj> inline void NerscStatisics(Lattice<vobj> & data,NerscField &header)
|
||||
{
|
||||
assert(data._grid->_ndimension==4);
|
||||
|
||||
for(int d=0;d<4;d++)
|
||||
header.dimension[d] = data._grid->_fdimensions[d];
|
||||
|
||||
// compute checksum and any physical properties contained for this type
|
||||
// header.checksum = NerscChecksum(data);
|
||||
|
||||
NerscPhysicalCharacteristics(data,header);
|
||||
}
|
||||
|
||||
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
// Now the meat: the object readers
|
||||
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
template<class vobj,class sobj,class fobj,class munger>
|
||||
inline void readNerscObject(Lattice<vobj> &Umu,std::string file,munger munge,int offset,std::string &format)
|
||||
{
|
||||
GridBase *grid = Umu._grid;
|
||||
|
||||
int ieee32big = (format == std::string("IEEE32BIG"));
|
||||
int ieee32 = (format == std::string("IEEE32"));
|
||||
int ieee64big = (format == std::string("IEEE64BIG"));
|
||||
int ieee64 = (format == std::string("IEEE64"));
|
||||
|
||||
// Find the location of each site and send to primary node
|
||||
// for(int site=0; site < Layout::vol(); ++site){
|
||||
// multi1d<int> coord = crtesn(site, Layout::lattSize());
|
||||
// for(int dd=0; dd<Nd; dd++){ /* dir */
|
||||
// cfg_in.readArray(su3_buffer, float_size, mat_size);
|
||||
//
|
||||
// Above from Chroma; defines loop order now that NERSC doc no longer
|
||||
// available (how short sighted is that?)
|
||||
{
|
||||
std::ifstream fin(file,std::ios::binary|std::ios::in);
|
||||
fin.seekg(offset);
|
||||
|
||||
Umu = zero;
|
||||
uint32_t csum=0;
|
||||
fobj file_object;
|
||||
sobj munged;
|
||||
|
||||
for(int t=0;t<grid->_fdimensions[3];t++){
|
||||
for(int z=0;z<grid->_fdimensions[2];z++){
|
||||
for(int y=0;y<grid->_fdimensions[1];y++){
|
||||
for(int x=0;x<grid->_fdimensions[0];x++){
|
||||
|
||||
std::vector<int> site({x,y,z,t});
|
||||
|
||||
if ( grid->IsBoss() ) {
|
||||
fin.read((char *)&file_object,sizeof(file_object));
|
||||
|
||||
if(ieee32big) be32toh((void *)&file_object,sizeof(file_object));
|
||||
if(ieee32) le32toh((void *)&file_object,sizeof(file_object));
|
||||
if(ieee64big) be64toh((void *)&file_object,sizeof(file_object));
|
||||
if(ieee64) le64toh((void *)&file_object,sizeof(file_object));
|
||||
|
||||
munge(file_object,munged,csum);
|
||||
}
|
||||
// The boss who read the file has their value poked
|
||||
pokeSite(munged,Umu,site);
|
||||
}}}}
|
||||
}
|
||||
}
|
||||
|
||||
template<class vobj,class sobj,class fobj,class munger>
|
||||
inline void writeNerscObject(Lattice<vobj> &Umu,std::string file,munger munge,int offset,
|
||||
int sequence,double lt,double pl)
|
||||
{
|
||||
GridBase *grid = Umu._grid;
|
||||
NerscField header;
|
||||
|
||||
//////////////////////////////////////////////////
|
||||
// First write the header; this is in wrong place
|
||||
//////////////////////////////////////////////////
|
||||
assert(grid->_ndimension == 4);
|
||||
for(int d=0;d<4;d++){
|
||||
header.dimension[d]=grid->_fdimensions[d];
|
||||
header.boundary [d]=std::string("PERIODIC");;
|
||||
}
|
||||
header.hdr_version=std::string("WHATDAHECK");
|
||||
// header.storage_format=storage_format<vobj>::string; // use template specialisation
|
||||
// header.data_type=data_type<vobj>::string;
|
||||
header.storage_format=std::string("debug");
|
||||
header.data_type =std::string("debug");
|
||||
|
||||
//FIXME; use template specialisation to fill these out
|
||||
header.link_trace =lt;
|
||||
header.plaquette =pl;
|
||||
header.checksum =0;
|
||||
|
||||
//
|
||||
header.sequence_number =sequence;
|
||||
header.ensemble_id =std::string("UKQCD");
|
||||
header.ensemble_label =std::string("UKQCD");
|
||||
header.creator =std::string("Tadahito");
|
||||
header.creator_hardware=std::string("BlueGene/Q");
|
||||
header.creation_date =std::string("AnnoDomini");
|
||||
header.archive_date =std::string("AnnoDomini");
|
||||
header.floating_point =std::string("IEEE64BIG");
|
||||
// header.data_start=;
|
||||
// unsigned int checksum;
|
||||
|
||||
//////////////////////////////////////////////////
|
||||
// Now write the body
|
||||
//////////////////////////////////////////////////
|
||||
{
|
||||
std::ofstream fout(file,std::ios::binary|std::ios::in);
|
||||
fout.seekp(offset);
|
||||
|
||||
Umu = zero;
|
||||
uint32_t csum=0;
|
||||
fobj file_object;
|
||||
sobj unmunged;
|
||||
for(int t=0;t<grid->_fdimensions[3];t++){
|
||||
for(int z=0;z<grid->_fdimensions[2];z++){
|
||||
for(int y=0;y<grid->_fdimensions[1];y++){
|
||||
for(int x=0;x<grid->_fdimensions[0];x++){
|
||||
std::vector<int> site({x,y,z,t});
|
||||
peekSite(unmunged,Umu,site);
|
||||
munge(unmunged,file_object,csum);
|
||||
// broadcast & insert
|
||||
fout.write((char *)&file_object,sizeof(file_object));
|
||||
}}}}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
inline void readNerscConfiguration(LatticeGaugeField &Umu,NerscField& header,std::string file)
|
||||
{
|
||||
GridBase *grid = Umu._grid;
|
||||
|
||||
int offset = readNerscHeader(file,Umu._grid,header);
|
||||
|
||||
std::string format(header.floating_point);
|
||||
|
||||
int ieee32big = (format == std::string("IEEE32BIG"));
|
||||
int ieee32 = (format == std::string("IEEE32"));
|
||||
int ieee64big = (format == std::string("IEEE64BIG"));
|
||||
int ieee64 = (format == std::string("IEEE64"));
|
||||
|
||||
// depending on datatype, set up munger;
|
||||
// munger is a function of <floating point, Real, data_type>
|
||||
if ( header.data_type == std::string("4D_SU3_GAUGE") ) {
|
||||
if ( ieee32 || ieee32big ) {
|
||||
readNerscObject<vLorentzColourMatrix, LorentzColourMatrix, LorentzColour2x3F>
|
||||
(Umu,file,
|
||||
Nersc3x2munger<LorentzColour2x3F,LorentzColourMatrix>(),
|
||||
offset,format);
|
||||
}
|
||||
if ( ieee64 || ieee64big ) {
|
||||
readNerscObject<vLorentzColourMatrix, LorentzColourMatrix, LorentzColour2x3D>
|
||||
(Umu,file,
|
||||
Nersc3x2munger<LorentzColour2x3D,LorentzColourMatrix>(),
|
||||
offset,format);
|
||||
}
|
||||
} else if ( header.data_type == std::string("4D_SU3_GAUGE_3X3") ) {
|
||||
if ( ieee32 || ieee32big ) {
|
||||
readNerscObject<vLorentzColourMatrix,LorentzColourMatrix,LorentzColourMatrixF>
|
||||
(Umu,file,NerscSimpleMunger<LorentzColourMatrixF,LorentzColourMatrix>(),offset,format);
|
||||
}
|
||||
if ( ieee64 || ieee64big ) {
|
||||
readNerscObject<vLorentzColourMatrix,LorentzColourMatrix,LorentzColourMatrixD>
|
||||
(Umu,file,NerscSimpleMunger<LorentzColourMatrixD,LorentzColourMatrix>(),offset,format);
|
||||
}
|
||||
} else {
|
||||
assert(0);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
template<class vobj>
|
||||
inline void writeNerscConfiguration(Lattice<vobj> &Umu,NerscField &header,std::string file)
|
||||
{
|
||||
GridBase &grid = Umu._grid;
|
||||
|
||||
NerscStatisics(Umu,header);
|
||||
|
||||
int offset = writeNerscHeader(header,file);
|
||||
|
||||
writeNerscObject(Umu,NerscSimpleMunger<vobj,vobj>(),offset);
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
#endif
|
||||
@ -42,6 +42,8 @@ namespace QCD {
|
||||
typedef iColourMatrix<Complex > ColourMatrix;
|
||||
typedef iSpinColourMatrix<Complex > SpinColourMatrix;
|
||||
typedef iLorentzColourMatrix<Complex > LorentzColourMatrix;
|
||||
typedef iLorentzColourMatrix<ComplexF > LorentzColourMatrixF;
|
||||
typedef iLorentzColourMatrix<ComplexD > LorentzColourMatrixD;
|
||||
|
||||
typedef iSpinVector<Complex > SpinVector;
|
||||
typedef iColourVector<Complex > ColourVector;
|
||||
@ -66,7 +68,7 @@ namespace QCD {
|
||||
|
||||
typedef Lattice<vTReal> LatticeReal;
|
||||
typedef Lattice<vTComplex> LatticeComplex;
|
||||
typedef Lattice<vInteger> LatticeInteger; // Predicates for "where"
|
||||
typedef Lattice<vTInteger> LatticeInteger; // Predicates for "where"
|
||||
|
||||
typedef Lattice<vColourMatrix> LatticeColourMatrix;
|
||||
typedef Lattice<vSpinMatrix> LatticeSpinMatrix;
|
||||
|
||||
84
tests/Grid_cshift.cc
Normal file
84
tests/Grid_cshift.cc
Normal file
@ -0,0 +1,84 @@
|
||||
#include <Grid.h>
|
||||
#include <parallelIO/GridNerscIO.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
std::vector<int> simd_layout({1,1,2,2});
|
||||
std::vector<int> mpi_layout ({2,2,2,2});
|
||||
std::vector<int> latt_size ({8,8,8,16});
|
||||
|
||||
GridCartesian Fine(latt_size,simd_layout,mpi_layout);
|
||||
GridRNG FineRNG(&Fine);
|
||||
|
||||
LatticeComplex U(&Fine);
|
||||
LatticeComplex ShiftU(&Fine);
|
||||
|
||||
LatticeComplex lex(&Fine);
|
||||
lex=zero;
|
||||
Integer stride =1;
|
||||
{
|
||||
double nrm;
|
||||
LatticeComplex coor(&Fine);
|
||||
|
||||
for(int d=0;d<4;d++){
|
||||
LatticeCoordinate(coor,d);
|
||||
lex = lex + coor*stride;
|
||||
stride=stride*latt_size[d];
|
||||
}
|
||||
U=lex;
|
||||
}
|
||||
|
||||
TComplex cm;
|
||||
|
||||
for(int dir=0;dir<4;dir++){
|
||||
for(int shift=0;shift<latt_size[dir];shift++){
|
||||
if ( Fine.IsBoss() )
|
||||
std::cout<<"Shifting by "<<shift<<" in direction"<<dir<<std::endl;
|
||||
|
||||
ShiftU = Cshift(U,dir,shift); // Shift everything
|
||||
|
||||
std::vector<int> coor(4);
|
||||
|
||||
for(coor[3]=0;coor[3]<latt_size[3];coor[3]++){
|
||||
for(coor[2]=0;coor[2]<latt_size[2];coor[2]++){
|
||||
for(coor[1]=0;coor[1]<latt_size[1];coor[1]++){
|
||||
for(coor[0]=0;coor[0]<latt_size[0];coor[0]++){
|
||||
|
||||
peekSite(cm,ShiftU,coor);
|
||||
|
||||
double nrm=norm2(U);
|
||||
|
||||
std::vector<int> scoor(coor);
|
||||
scoor[dir] = (scoor[dir]+shift)%latt_size[dir];
|
||||
|
||||
Integer slex = scoor[0]
|
||||
+ latt_size[0]*scoor[1]
|
||||
+ latt_size[0]*latt_size[1]*scoor[2]
|
||||
+ latt_size[0]*latt_size[1]*latt_size[2]*scoor[3];
|
||||
|
||||
Complex scm(slex);
|
||||
|
||||
nrm = abs(scm-cm()()());
|
||||
std::vector<int> peer(4);
|
||||
int index=real(cm);
|
||||
Fine.CoorFromIndex(peer,index,latt_size);
|
||||
|
||||
if (nrm > 0){
|
||||
cout<<"FAIL shift "<< shift<<" in dir "<< dir<<" ["<<coor[0]<<","<<coor[1]<<","<<coor[2]<<","<<coor[3]<<"] = "<< cm()()()<<" expect "<<scm<<" "<<nrm<<endl;
|
||||
cout<<"Got "<<index<<" " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
|
||||
index=real(scm);
|
||||
Fine.CoorFromIndex(peer,index,latt_size);
|
||||
cout<<"Expect "<<index<<" " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
|
||||
}
|
||||
}}}}
|
||||
}
|
||||
}
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
62
tests/Grid_nersc_io.cc
Normal file
62
tests/Grid_nersc_io.cc
Normal file
@ -0,0 +1,62 @@
|
||||
#include <Grid.h>
|
||||
#include <parallelIO/GridNerscIO.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
std::vector<int> simd_layout({1,1,2,2});
|
||||
std::vector<int> mpi_layout ({2,1,1,2});
|
||||
std::vector<int> latt_size ({16,16,16,32});
|
||||
|
||||
GridCartesian Fine(latt_size,simd_layout,mpi_layout);
|
||||
GridRNG FineRNG(&Fine);
|
||||
LatticeGaugeField Umu(&Fine);
|
||||
|
||||
std::vector<LatticeColourMatrix> U(4,&Fine);
|
||||
|
||||
NerscField header;
|
||||
|
||||
std::string file("./ckpoint_lat.4000");
|
||||
readNerscConfiguration(Umu,header,file);
|
||||
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
U[mu] = peekIndex<3>(Umu,mu);
|
||||
}
|
||||
|
||||
// Painful ; fix syntactical niceness
|
||||
LatticeComplex LinkTrace(&Fine);
|
||||
LinkTrace=zero;
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
LinkTrace = LinkTrace + trace(U[mu]);
|
||||
}
|
||||
|
||||
// (1+2+3)=6 = N(N-1)/2 terms
|
||||
LatticeComplex Plaq(&Fine);
|
||||
Plaq = zero;
|
||||
for(int mu=1;mu<Nd;mu++){
|
||||
for(int nu=0;nu<mu;nu++){
|
||||
Plaq = Plaq + trace(U[mu]*Cshift(U[nu],mu,1)*adj(Cshift(U[mu],nu,1))*adj(U[nu]));
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
double vol = Fine.gSites();
|
||||
|
||||
Complex PlaqScale(1.0/vol/6.0/3.0);
|
||||
TComplex Tp = sum(Plaq);
|
||||
Complex p = TensorRemove(Tp);
|
||||
std::cout << "calculated plaquettes " <<p*PlaqScale<<std::endl;
|
||||
|
||||
Complex LinkTraceScale(1.0/vol/4.0/3.0);
|
||||
TComplex Tl = sum(LinkTrace);
|
||||
Complex l = TensorRemove(Tl);
|
||||
std::cout << "calculated link trace " <<l*LinkTraceScale<<std::endl;
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
@ -5,10 +5,16 @@ AM_LDFLAGS = -L$(top_srcdir)/lib
|
||||
#
|
||||
# Test code
|
||||
#
|
||||
bin_PROGRAMS = Grid_main test_Grid_stencil
|
||||
bin_PROGRAMS = Grid_main test_Grid_stencil Grid_nersc_io Grid_cshift
|
||||
|
||||
Grid_main_SOURCES = Grid_main.cc
|
||||
Grid_main_LDADD = -lGrid
|
||||
|
||||
Grid_nersc_io_SOURCES = Grid_nersc_io.cc
|
||||
Grid_nersc_io_LDADD = -lGrid
|
||||
|
||||
Grid_cshift_SOURCES = Grid_cshift.cc
|
||||
Grid_cshift_LDADD = -lGrid
|
||||
|
||||
test_Grid_stencil_SOURCES = test_Grid_stencil.cc
|
||||
test_Grid_stencil_LDADD = -lGrid
|
||||
|
||||
Loading…
Reference in New Issue
Block a user