diff --git a/Hadrons/Distil.hpp b/Hadrons/Distil.hpp index 9d4423ad..ef2b3cbf 100644 --- a/Hadrons/Distil.hpp +++ b/Hadrons/Distil.hpp @@ -2,13 +2,13 @@ Grid physics library, www.github.com/paboyle/Grid - Source file: Hadrons/Modules/MDistil/Distil.hpp - + Source file: Hadrons/Distil.hpp + Copyright (C) 2015-2019 Author: Felix Erben Author: Michael Marshall - + This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or @@ -27,295 +27,132 @@ *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_MDistil_Distil_hpp_ -#define Hadrons_MDistil_Distil_hpp_ +#ifndef Hadrons_Distil_hpp_ +#define Hadrons_Distil_hpp_ #include -#include -#include -#include #include -#include -#include BEGIN_HADRONS_NAMESPACE -BEGIN_MODULE_NAMESPACE(MDistil) + +/****************************************************************************** + NamedTensor + This is an Eigen::Tensor of type Scalar_ and rank NumIndices_ (row-major order) + They can be persisted to disk in tag Name_, and IndexNames are validated on load. + TODO: WHAT TO SAVE / VALIDATE ON LOAD (Override to warn instead of assert on load) + Ensemble string + Configuration number + Noise unique string + Distillation parameters + ******************************************************************************/ + +extern const std::string NamedTensorFileExtension; + +template &IndexNames_> +class NamedTensor : Serializable +{ +public: + using Scalar = Scalar_; + static constexpr int NumIndices = NumIndices_; + using ET = Eigen::Tensor; + using Index = typename ET::Index; + GRID_SERIALIZABLE_CLASS_MEMBERS(NamedTensor, + ET, tensor, + std::vector, IndexNames ); + + // Get the default index names as std::vector + std::vector DefaultIndexNames() + { + std::vector names{NumIndices_}; + for (std::size_t i = 0; i < NumIndices_; i++) + names[i] = IndexNames_[i]; + return names; + } + + // Default constructor (assumes tensor will be loaded from file) + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NamedTensor() : IndexNames{DefaultIndexNames()} {} + + // Construct a named tensor explicitly specifying size of each dimension + template + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NamedTensor(Eigen::Index firstDimension, IndexTypes... otherDimensions) + : tensor(firstDimension, otherDimensions...), IndexNames{DefaultIndexNames()} + { + assert(sizeof...(otherDimensions) + 1 == NumIndices_ && "NamedTensor: dimensions != tensor rank"); + } + + // Do my index names match the default for my type? + bool ValidateIndexNames() const + { + bool bSame{ IndexNames.size() == NumIndices_ }; + for( std::size_t i = 0; bSame && i < NumIndices_; i++ ) + { + bSame = IndexNames[i].size() == IndexNames_[i].size() + && std::equal( IndexNames[i].begin(), IndexNames[i].end(), IndexNames_[i].begin(), + [](const char & c1, const char & c2){ return c1 == c2 || std::toupper(c1) == std::toupper(c2); }); + } + return bSame; + } + +#ifdef HAVE_HDF5 + using Default_Reader = Grid::Hdf5Reader; + using Default_Writer = Grid::Hdf5Writer; +#else + using Default_Reader = Grid::BinaryReader; + using Default_Writer = Grid::BinaryWriter; +#endif + + template void write(Writer &w, const std::string &Tag = Name_) const + { write(w, Tag, *this); } + + inline void write(const std::string &filename, const std::string &Tag = Name_) const + { + std::string sFileName{filename}; + sFileName.append( NamedTensorFileExtension ); + LOG(Message) << "Writing " << Name_ << " to file " << sFileName << " tag " << Tag << std::endl; + Default_Writer w( sFileName ); + write( w, Tag ); + } + + // Read and validate index names + template void read(Reader &r, bool bValidate = true, const std::string &Tag = Name_) + { + // Grab index names and dimensions + std::vector OldIndexNames{std::move(IndexNames)}; + typename ET::Dimensions OldDimensions{tensor.dimensions()}; + read(r, Tag, *this); + const typename ET::Dimensions & NewDimensions{tensor.dimensions()}; + for (int i = 0; i < NumIndices_; i++) + assert(OldDimensions[i] == 0 || OldDimensions[i] == NewDimensions[i] && "NamedTensor::read dimension size"); + if (bValidate) + assert(ValidateIndexNames() && "NamedTensor::read dimension name"); + } + + inline void read (const std::string &filename, bool bValidate = true, const std::string &Tag = Name_) + { + std::string sFileName{filename}; + sFileName.append( NamedTensorFileExtension ); + LOG(Message) << "Reading " << Name_ << " from file " << sFileName << " tag " << Tag << std::endl; + Default_Reader r(sFileName); + read(r, bValidate, Tag); + } +}; /****************************************************************************** Common elements for distillation ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MDistil) + +//Eigenvectors of the Laplacian using LapEvecs = Grid::Hadrons::EigenPack; -// Noise vector index order: nnoise, nt, nvec, ns +// Noise vector (index order: nnoise, nt, nvec, ns) using NoiseTensor = Eigen::Tensor; -struct DistilParameters: Serializable { - GRID_SERIALIZABLE_CLASS_MEMBERS(DistilParameters, - int, nnoise, - int, tsrc, - std::string, TI, - std::string, LI, - std::string, SI ) - DistilParameters() = default; - template DistilParameters(Reader& Reader){read(Reader,"Distil",*this);} - - // Numeric parameter is allowed to be empty (in which case it = Default), - // but assert during setup() if specified but not numeric - - static int ParameterDefault( const std::string & s, int Default, bool bCalledFromSetup ) - { - int i = Default; - if( s.length() > 0 ) { - std::istringstream ss( s ); - ss >> i; - if( bCalledFromSetup ) - assert( !ss.fail() && "Parameter should either be empty or integer" ); - } - return i; - } -}; - -#define DISTIL_PARAMETERS_DEFINE( inSetup ) \ -const int Nt{env().getDim(Tdir)}; \ -const int nvec{par().nvec}; \ -const int nnoise{par().Distil.nnoise}; \ -const int tsrc{par().Distil.tsrc}; \ -const int TI{Hadrons::MDistil::DistilParameters::ParameterDefault(par().Distil.TI, Nt, inSetup)}; \ -const int LI{Hadrons::MDistil::DistilParameters::ParameterDefault(par().Distil.LI, nvec, inSetup)}; \ -const int SI{Hadrons::MDistil::DistilParameters::ParameterDefault(par().Distil.SI, Ns, inSetup)}; \ -const bool full_tdil{ TI == Nt }; \ -const bool exact_distillation{ full_tdil && LI == nvec }; \ -const int Nt_inv{ full_tdil ? 1 : TI } - -/****************************************************************************** - Default for distillation file operations. For now only used by NamedTensor - ******************************************************************************/ - -#ifdef HAVE_HDF5 -using Default_Reader = Grid::Hdf5Reader; -using Default_Writer = Grid::Hdf5Writer; -static const char * FileExtension = ".h5"; -#else -using Default_Reader = Grid::BinaryReader; -using Default_Writer = Grid::BinaryWriter; -static const char * FileExtension = ".dat"; -#endif - -/****************************************************************************** - NamedTensor object - This is an Eigen::Tensor of type Scalar_ and rank NumIndices_ (row-major order) - They can be persisted to disk - IndexNames contains one name for each index, and IndexNames are validated on load. - WHAT TO SAVE / VALIDATE ON LOAD (Override to warn instead of assert on load) - Ensemble string - Configuration number - Noise unique string - Distillation parameters - - ******************************************************************************/ - -template -class NamedTensor : Serializable -{ -public: - using Scalar = Scalar_; - static constexpr int NumIndices = NumIndices_; - using ET = Eigen::Tensor; - using Index = typename ET::Index; - GRID_SERIALIZABLE_CLASS_MEMBERS(NamedTensor - , ET, tensor - , std::vector, IndexNames - ); - // Named tensors are intended to be a superset of Eigen tensor - inline operator ET&() { return tensor; } - template - inline const Scalar_& operator()(const std::array &Indices) const - { return tensor.operator()(Indices); } - inline Scalar_& operator()(const std::array &Indices) - { return tensor.operator()(Indices); } - template - inline const Scalar_& operator()(Eigen::Index firstDimension, IndexTypes... otherDimensions) const - { - // The number of indices used to access a tensor coefficient must be equal to the rank of the tensor. - assert(sizeof...(otherDimensions) + 1 == NumIndices_ && "NamedTensor: dimensions != tensor rank"); - return tensor.operator()(std::array{{firstDimension, otherDimensions...}}); - } - template - inline Scalar_& operator()(Eigen::Index firstDimension, IndexTypes... otherDimensions) - { - // The number of indices used to access a tensor coefficient must be equal to the rank of the tensor. - assert(sizeof...(otherDimensions) + 1 == NumIndices_ && "NamedTensor: dimensions != tensor rank"); - return tensor.operator()(std::array{{firstDimension, otherDimensions...}}); - } - - // Construct a named tensor explicitly specifying size of each dimension - template - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NamedTensor(const std::array &IndexNames_, Eigen::Index firstDimension, IndexTypes... otherDimensions) - : tensor(firstDimension, otherDimensions...), IndexNames{NumIndices} - { - // The number of dimensions used to construct a tensor must be equal to the rank of the tensor. - assert(sizeof...(otherDimensions) + 1 == NumIndices_ && "NamedTensor: dimensions != tensor rank"); - for( int i = 0; i < NumIndices_; i++ ) - IndexNames[i] = IndexNames_[i]; - } - - // Default constructor (assumes tensor will be loaded from file) - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NamedTensor() : IndexNames{NumIndices_} {} - - // Construct a named tensor without specifying size of each dimension (because it will be loaded from file) - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NamedTensor(const std::array &IndexNames_) - : IndexNames{NumIndices_} - { - for( int i = 0; i < NumIndices_; i++ ) - IndexNames[i] = IndexNames_[i]; - } - - bool ValidateIndexNames( std::size_t NumNames, const std::string * MatchNames ) const; - - // Read/Write in any format - template inline void read (Reader &r, const char * pszTag = nullptr); - template inline void write(Writer &w, const char * pszTag = nullptr) const; - // Read/Write in default format, i.e. HDF5 if present, else binary - inline void read (const char * filename, const char * pszTag = nullptr); - inline void write(const char * filename, const char * pszTag = nullptr) const; -}; - -// Is this a named tensor -template struct is_named_tensor : public std::false_type {}; -template struct is_named_tensor> : public std::true_type {}; -template struct is_named_tensor, T>::value>::type> : public std::true_type {}; - -/****************************************************************************** - PerambTensor object - ******************************************************************************/ - -using PerambTensor = NamedTensor; -static const std::array PerambIndexNames{"nT", "nVec", "LI", "nNoise", "nT_inv", "SI"}; - -/****************************************************************************** - Write NamedTensor - ******************************************************************************/ - -template -template -void NamedTensor::write(Writer &w, const char * pszTag)const{ - if( pszTag == nullptr ) - pszTag = "NamedTensor"; - LOG(Message) << "Writing NamedTensor to tag " << pszTag << std::endl; - write(w, pszTag, *this); -} - -template -void NamedTensor::write(const char * filename, const char * pszTag)const{ - std::string sFileName{filename}; - sFileName.append( MDistil::FileExtension ); - LOG(Message) << "Writing NamedTensor to file " << sFileName << std::endl; - MDistil::Default_Writer w( sFileName ); - write( w, pszTag ); -} - -/****************************************************************************** - Validate named tensor index names - ******************************************************************************/ - -template -bool NamedTensor::ValidateIndexNames( std::size_t NumNames, const std::string * MatchNames ) const { - bool bSame{ NumNames == NumIndices_ && IndexNames.size() == NumIndices_ }; - for( std::size_t i = 0; bSame && i < NumIndices_; i++ ) - { - bSame = MatchNames[i].size() == IndexNames[i].size() - && std::equal( MatchNames[i].begin(), MatchNames[i].end(), IndexNames[i].begin(), - [](const char & c1, const char & c2){ return c1 == c2 || std::toupper(c1) == std::toupper(c2); }); - } - return bSame; -} - -/****************************************************************************** - Read NamedTensor - ******************************************************************************/ - -template -template -void NamedTensor::read(Reader &r, const char * pszTag) { - if( pszTag == nullptr ) - pszTag = "NamedTensor"; - // Grab index names and dimensions - std::vector OldIndexNames{std::move(IndexNames)}; - typename ET::Dimensions OldDimensions{tensor.dimensions()}; - LOG(Message) << "Reading NamedTensor from tag " << pszTag << std::endl; - read(r, pszTag, *this); - const typename ET::Dimensions & NewDimensions{tensor.dimensions()}; - for( int i = 0; i < NumIndices_; i++ ) - assert(OldDimensions[i] == 0 || OldDimensions[i] == NewDimensions[i] && "NamedTensor::load dimension size"); - assert( ValidateIndexNames( OldIndexNames.size(), &OldIndexNames[0] ) && "NamedTensor::load dimension name" ); -} - -template -void NamedTensor::read(const char * filename, const char * pszTag) { - std::string sFileName{filename}; - sFileName.append( MDistil::FileExtension ); - LOG(Message) << "Reading NamedTensor from file " << sFileName << std::endl; - MDistil::Default_Reader r( sFileName ); - read( r, pszTag ); -} - -/****************************************************************************** - Make a lower dimensional grid in preparation for local slice operations - ******************************************************************************/ - -inline GridCartesian * MakeLowerDimGrid( GridCartesian * gridHD ) -{ - int nd{static_cast(gridHD->_ndimension)}; - Coordinate latt_size = gridHD->_gdimensions; - latt_size[nd-1] = 1; - Coordinate simd_layout = GridDefaultSimd(nd-1, vComplex::Nsimd()); - simd_layout.push_back( 1 ); - Coordinate mpi_layout = gridHD->_processors; - mpi_layout[nd-1] = 1; - GridCartesian * gridLD = new GridCartesian(latt_size,simd_layout,mpi_layout,*gridHD); - return gridLD; -} - -/************************************************************************************* - Rotate eigenvectors into our phase convention - First component of first eigenvector is real and positive - *************************************************************************************/ - -inline void RotateEigen(std::vector & evec) -{ - ColourVector cv0; - auto grid = evec[0].Grid(); - Coordinate siteFirst(grid->Nd(),0); - peekSite(cv0, evec[0], siteFirst); - Grid::Complex cplx0 = cv0()()(0); - if( cplx0.imag() == 0 ) - std::cout << GridLogMessage << "RotateEigen() : Site 0 : " << cplx0 << " => already meets phase convention" << std::endl; - else { - const Real cplx0_mag = Grid::abs(cplx0); -#ifdef GRID_NVCC - const Grid::Complex phase = thrust::conj(cplx0 / cplx0_mag); - const Real argphase = thrust::arg(phase); -#else - const Grid::Complex phase = std::conj(cplx0 / cplx0_mag); - const Real argphase = std::arg(phase); -#endif - std::cout << GridLogMessage << "RotateEigen() : Site 0 : |" << cplx0 << "|=" << cplx0_mag << " => phase=" << (argphase / 3.14159265) << " pi" << std::endl; - { - // TODO: Only really needed on the master slice - for( int k = 0 ; k < evec.size() ; k++ ) - evec[k] *= phase; - if(grid->IsBoss()){ - for( int c = 0 ; c < Nc ; c++ ) - cv0()()(c) *= phase; - cplx0.imag(0); // This assumes phase convention is real, positive (so I get rid of rounding error) - //pokeSite(cv0, evec[0], siteFirst); - pokeLocalSite(cv0, evec[0], siteFirst); - } - } - } -} +extern const std::string PerambTensorName; +extern const std::array PerambIndexNames; +using PerambTensor = NamedTensor; END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_MDistil_Distil_hpp_ +#endif // Hadrons_Distil_hpp_ diff --git a/Hadrons/Modules/MDistil/DistilCommon.hpp b/Hadrons/Modules/MDistil/DistilCommon.hpp new file mode 100644 index 00000000..bb224447 --- /dev/null +++ b/Hadrons/Modules/MDistil/DistilCommon.hpp @@ -0,0 +1,143 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: Hadrons/Modules/MDistil/DistilCommon.hpp + + Copyright (C) 2015-2019 + + Author: Felix Erben + Author: Michael Marshall + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_MDistil_DistilCommon_hpp_ +#define Hadrons_MDistil_DistilCommon_hpp_ + +#include +#include +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE +BEGIN_MODULE_NAMESPACE(MDistil) + +/****************************************************************************** + Common elements for distillation + ******************************************************************************/ + +struct DistilParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(DistilParameters, + int, nnoise, + int, tsrc, + std::string, TI, + std::string, LI, + std::string, SI ) + DistilParameters() = default; + template DistilParameters(Reader& Reader){read(Reader,"Distil",*this);} + + // Numeric parameter is allowed to be empty (in which case it = Default), + // but assert during setup() if specified but not numeric + + static int ParameterDefault( const std::string & s, int Default, bool bCalledFromSetup ) + { + int i = Default; + if( s.length() > 0 ) { + std::istringstream ss( s ); + ss >> i; + if( bCalledFromSetup ) + assert( !ss.fail() && "Parameter should either be empty or integer" ); + } + return i; + } +}; + +#define DISTIL_PARAMETERS_DEFINE( inSetup ) \ +const int Nt{env().getDim(Tdir)}; \ +const int nvec{par().nvec}; \ +const int nnoise{par().Distil.nnoise}; \ +const int tsrc{par().Distil.tsrc}; \ +const int TI{Hadrons::MDistil::DistilParameters::ParameterDefault(par().Distil.TI, Nt, inSetup)}; \ +const int LI{Hadrons::MDistil::DistilParameters::ParameterDefault(par().Distil.LI, nvec, inSetup)}; \ +const int SI{Hadrons::MDistil::DistilParameters::ParameterDefault(par().Distil.SI, Ns, inSetup)}; \ +const bool full_tdil{ TI == Nt }; \ +const bool exact_distillation{ full_tdil && LI == nvec }; \ +const int Nt_inv{ full_tdil ? 1 : TI } + +/****************************************************************************** + Make a lower dimensional grid in preparation for local slice operations + ******************************************************************************/ + +inline GridCartesian * MakeLowerDimGrid( GridCartesian * gridHD ) +{ + int nd{static_cast(gridHD->_ndimension)}; + Coordinate latt_size = gridHD->_gdimensions; + latt_size[nd-1] = 1; + Coordinate simd_layout = GridDefaultSimd(nd-1, vComplex::Nsimd()); + simd_layout.push_back( 1 ); + Coordinate mpi_layout = gridHD->_processors; + mpi_layout[nd-1] = 1; + GridCartesian * gridLD = new GridCartesian(latt_size,simd_layout,mpi_layout,*gridHD); + return gridLD; +} + +/************************************************************************************* + Rotate eigenvectors into our phase convention + First component of first eigenvector is real and positive + *************************************************************************************/ + +inline void RotateEigen(std::vector & evec) +{ + ColourVector cv0; + auto grid = evec[0].Grid(); + Coordinate siteFirst(grid->Nd(),0); + peekSite(cv0, evec[0], siteFirst); + Grid::Complex cplx0 = cv0()()(0); + if( cplx0.imag() == 0 ) + std::cout << GridLogMessage << "RotateEigen() : Site 0 : " << cplx0 << " => already meets phase convention" << std::endl; + else { + const Real cplx0_mag = Grid::abs(cplx0); +#ifdef GRID_NVCC + const Grid::Complex phase = thrust::conj(cplx0 / cplx0_mag); + const Real argphase = thrust::arg(phase); +#else + const Grid::Complex phase = std::conj(cplx0 / cplx0_mag); + const Real argphase = std::arg(phase); +#endif + std::cout << GridLogMessage << "RotateEigen() : Site 0 : |" << cplx0 << "|=" << cplx0_mag << " => phase=" << (argphase / 3.14159265) << " pi" << std::endl; + { + // TODO: Only really needed on the master slice + for( int k = 0 ; k < evec.size() ; k++ ) + evec[k] *= phase; + if(grid->IsBoss()){ + for( int c = 0 ; c < Nc ; c++ ) + cv0()()(c) *= phase; + cplx0.imag(0); // This assumes phase convention is real, positive (so I get rid of rounding error) + //pokeSite(cv0, evec[0], siteFirst); + pokeLocalSite(cv0, evec[0], siteFirst); + } + } + } +} + +END_MODULE_NAMESPACE +END_HADRONS_NAMESPACE +#endif // Hadrons_MDistil_DistilCommon_hpp_ diff --git a/Hadrons/Modules/MDistil/DistilVectors.hpp b/Hadrons/Modules/MDistil/DistilVectors.hpp index d7097a06..9ac3750d 100644 --- a/Hadrons/Modules/MDistil/DistilVectors.hpp +++ b/Hadrons/Modules/MDistil/DistilVectors.hpp @@ -30,24 +30,15 @@ #ifndef Hadrons_MDistil_DistilVectors_hpp_ #define Hadrons_MDistil_DistilVectors_hpp_ -#include -#include -#include -#include -#include -#include -#include - -// These are members of Distillation -#include +#include BEGIN_HADRONS_NAMESPACE +BEGIN_MODULE_NAMESPACE(MDistil) /****************************************************************************** * DistilVectors * * (Create rho and/or phi vectors) * ******************************************************************************/ -BEGIN_MODULE_NAMESPACE(MDistil) class DistilVectorsPar: Serializable { @@ -166,7 +157,7 @@ void TDistilVectors::setup(void) auto &perambulator = envGet(PerambTensor, PerambulatorName); // We expect the perambulator to have been created with these indices - assert( perambulator.ValidateIndexNames( PerambIndexNames.size(), &PerambIndexNames[0] ) && "Perambulator index names bad" ); + assert( perambulator.ValidateIndexNames() && "Perambulator index names bad" ); const int Nt{ env().getDim(Tdir) }; assert( Nt == static_cast( perambulator.tensor.dimension(0) ) && "PerambTensor time dimensionality bad" ); @@ -287,7 +278,7 @@ void TDistilVectors::execute(void) sink_tslice=0; for (int ivec = 0; ivec < nvec; ivec++) { ExtractSliceLocal(evec3d,epack.evec[ivec],0,t-Ntfirst,Tdir); - sink_tslice += evec3d * perambulator(t, ivec, dk, inoise,dt,ds); + sink_tslice += evec3d * perambulator.tensor(t, ivec, dk, inoise,dt,ds); } InsertSliceLocal(sink_tslice,phi[vecindex],0,t-Ntfirst,Tdir); } diff --git a/Hadrons/Modules/MDistil/LapEvec.hpp b/Hadrons/Modules/MDistil/LapEvec.hpp index 5ac8e258..022e88d1 100644 --- a/Hadrons/Modules/MDistil/LapEvec.hpp +++ b/Hadrons/Modules/MDistil/LapEvec.hpp @@ -30,16 +30,9 @@ #ifndef Hadrons_MDistil_LapEvec_hpp_ #define Hadrons_MDistil_LapEvec_hpp_ -#include -#include -#include -#include - -// These are members of Distillation -#include +#include BEGIN_HADRONS_NAMESPACE - BEGIN_MODULE_NAMESPACE(MDistil) /****************************************************************************** diff --git a/Hadrons/Modules/MDistil/Noises.hpp b/Hadrons/Modules/MDistil/Noises.hpp index a8c612cd..58c658dc 100644 --- a/Hadrons/Modules/MDistil/Noises.hpp +++ b/Hadrons/Modules/MDistil/Noises.hpp @@ -30,23 +30,14 @@ #ifndef Hadrons_MDistil_Noises_hpp_ #define Hadrons_MDistil_Noises_hpp_ -#include -#include -#include -#include -#include -#include -#include - -// These are members of Distillation - #include +#include BEGIN_HADRONS_NAMESPACE +BEGIN_MODULE_NAMESPACE(MDistil) /****************************************************************************** * Noises * ******************************************************************************/ -BEGIN_MODULE_NAMESPACE(MDistil) class NoisesPar: Serializable { diff --git a/Hadrons/Modules/MDistil/PerambFromSolve.hpp b/Hadrons/Modules/MDistil/PerambFromSolve.hpp index fc580b55..3562c4d2 100644 --- a/Hadrons/Modules/MDistil/PerambFromSolve.hpp +++ b/Hadrons/Modules/MDistil/PerambFromSolve.hpp @@ -30,23 +30,14 @@ #ifndef Hadrons_MDistil_PerambFromSolve_hpp_ #define Hadrons_MDistil_PerambFromSolve_hpp_ -#include -#include -#include -#include -#include -#include -#include +#include -// These are members of Distillation -#include - BEGIN_HADRONS_NAMESPACE +BEGIN_MODULE_NAMESPACE(MDistil) /****************************************************************************** * PerambFromSolve * ******************************************************************************/ -BEGIN_MODULE_NAMESPACE(MDistil) class PerambFromSolvePar: Serializable { @@ -126,7 +117,7 @@ void TPerambFromSolve::setup(void) const int LI_reduced{ Hadrons::MDistil::DistilParameters::ParameterDefault( par().LI_reduced, LI, true) }; grid4d = env().getGrid(); grid3d = MakeLowerDimGrid(grid4d); - envCreate(PerambTensor, getName(), 1, PerambIndexNames,Nt,nvec_reduced,LI_reduced,nnoise,Nt_inv,SI); + envCreate(PerambTensor, getName(), 1, Nt,nvec_reduced,LI_reduced,nnoise,Nt_inv,SI); envCreate(NoiseTensor, getName() + "_noise", 1, nnoise, Nt, nvec, Ns ); envTmp(LatticeColourVector, "result_3d",1,LatticeColourVector(grid3d)); envTmp(LatticeColourVector, "evec3d",1,LatticeColourVector(grid3d)); @@ -171,8 +162,8 @@ void TPerambFromSolve::execute(void) ExtractSliceLocal(result_3d,result_nospin,0,t-Ntfirst,Tdir); for (int ivec = 0; ivec < nvec_reduced; ivec++) { ExtractSliceLocal(evec3d,epack.evec[ivec],0,t-Ntfirst,Tdir); - pokeSpin(perambulator(t, ivec, dk, inoise,dt,ds),static_cast(innerProduct(evec3d, result_3d)),is); - LOG(Message) << "perambulator(t, ivec, dk, inoise,dt,ds)(is) = (" << t << "," << ivec << "," << dk << "," << inoise << "," << dt << "," << ds << ")(" << is << ") = " << perambulator(t, ivec, dk, inoise,dt,ds)()(is)() << std::endl; + pokeSpin(perambulator.tensor(t, ivec, dk, inoise,dt,ds),static_cast(innerProduct(evec3d, result_3d)),is); + LOG(Message) << "perambulator(t, ivec, dk, inoise,dt,ds)(is) = (" << t << "," << ivec << "," << dk << "," << inoise << "," << dt << "," << ds << ")(" << is << ") = " << perambulator.tensor(t, ivec, dk, inoise,dt,ds)()(is)() << std::endl; } } } diff --git a/Hadrons/Modules/MDistil/Perambulator.cc b/Hadrons/Modules/MDistil/Perambulator.cc index 055b8814..b141eb95 100644 --- a/Hadrons/Modules/MDistil/Perambulator.cc +++ b/Hadrons/Modules/MDistil/Perambulator.cc @@ -34,3 +34,14 @@ using namespace Hadrons; using namespace MDistil; template class Grid::Hadrons::MDistil::TPerambulator; + +// Global constants for distillation + +const std::string Grid::Hadrons::MDistil::PerambTensorName{ "Perambulator" }; +const std::array Grid::Hadrons::MDistil::PerambIndexNames{"nT", "nVec", "LI", "nNoise", "nT_inv", "SI"}; + +#ifdef HAVE_HDF5 +extern const std::string Grid::Hadrons::NamedTensorFileExtension{".h5"}; +#else +extern const std::string Grid::Hadrons::NamedTensorFileExtension{".dat"}; +#endif diff --git a/Hadrons/Modules/MDistil/Perambulator.hpp b/Hadrons/Modules/MDistil/Perambulator.hpp index 86740283..52ec7501 100644 --- a/Hadrons/Modules/MDistil/Perambulator.hpp +++ b/Hadrons/Modules/MDistil/Perambulator.hpp @@ -30,7 +30,7 @@ #ifndef Hadrons_MDistil_Perambulator_hpp_ #define Hadrons_MDistil_Perambulator_hpp_ -#include +#include BEGIN_HADRONS_NAMESPACE BEGIN_MODULE_NAMESPACE(MDistil) @@ -129,7 +129,7 @@ void TPerambulator::setup(void) if( !UnsmearedSinkFileName.empty() ) bool bMulti = ( Hadrons::MDistil::DistilParameters::ParameterDefault( par().UnsmearedSinkMultiFile, 1, true ) != 0 ); - envCreate(PerambTensor, getName(), 1, PerambIndexNames,Nt,nvec,LI,nnoise,Nt_inv,SI); + envCreate(PerambTensor, getName(), 1, Nt,nvec,LI,nnoise,Nt_inv,SI); envCreate(std::vector, getName() + "_unsmeared_sink", 1, nnoise*LI*Ns*Nt_inv, envGetGrid(FermionField)); @@ -244,7 +244,7 @@ void TPerambulator::execute(void) for (int ivec = 0; ivec < nvec; ivec++) { ExtractSliceLocal(evec3d,epack.evec[ivec],0,t-Ntfirst,Tdir); - pokeSpin(perambulator(t, ivec, dk, inoise,dt,ds),static_cast(innerProduct(evec3d, result_3d)),is); + pokeSpin(perambulator.tensor(t, ivec, dk, inoise,dt,ds),static_cast(innerProduct(evec3d, result_3d)),is); } } } diff --git a/Hadrons/Modules/MIO/LoadPerambulator.hpp b/Hadrons/Modules/MIO/LoadPerambulator.hpp index 558bf145..9d30f9fa 100644 --- a/Hadrons/Modules/MIO/LoadPerambulator.hpp +++ b/Hadrons/Modules/MIO/LoadPerambulator.hpp @@ -30,18 +30,14 @@ #ifndef Hadrons_MIO_LoadPerambulator_hpp_ #define Hadrons_MIO_LoadPerambulator_hpp_ -#include -#include -#include -#include -#include +#include BEGIN_HADRONS_NAMESPACE +BEGIN_MODULE_NAMESPACE(MIO) /****************************************************************************** * LoadPerambulator * ******************************************************************************/ -BEGIN_MODULE_NAMESPACE(MIO) class LoadPerambulatorPar: Serializable { @@ -103,7 +99,7 @@ void TLoadPerambulator::setup(void) { DISTIL_PARAMETERS_DEFINE( true ); //std::array sIndexNames{"Nt", "nvec", "LI", "nnoise", "Nt_inv", "SI"}; - envCreate(MDistil::PerambTensor, getName(), 1, MDistil::PerambIndexNames,Nt,nvec,LI,nnoise,Nt_inv,SI); + envCreate(MDistil::PerambTensor, getName(), 1, Nt,nvec,LI,nnoise,Nt_inv,SI); } // execution /////////////////////////////////////////////////////////////////// diff --git a/tests/IO/Test_serialisation.cc b/tests/IO/Test_serialisation.cc index 7f3e5bca..fd0b8b52 100644 --- a/tests/IO/Test_serialisation.cc +++ b/tests/IO/Test_serialisation.cc @@ -30,6 +30,7 @@ Author: Michael Marshall /* END LEGAL */ #include #include +#include using namespace Grid; @@ -252,6 +253,14 @@ void tensorConvTestFn(GridSerialRNG &rng, const std::string label) } #define tensorConvTest(rng, type) tensorConvTestFn(rng, #type) +static const std::string TheyCallMeLurve{"That's not my name"}; +static const std::string Coco{"Coco the dancing seal"}; +static const std::array GreatestShowOnEarth{ "Burnum", "Burnum" }; + +using MyPerambTensor = Grid::Hadrons::NamedTensor; +static const std::string MPTName{ "Perambulator" }; +static const std::array MPTIndexNames{"nT", "nVec", "LI", "nNoise", "nT_inv", "SI"}; +using MyPerambTensor2 = Grid::Hadrons::NamedTensor; int main(int argc,char **argv) { @@ -383,5 +392,33 @@ int main(int argc,char **argv) XmlWriter HMCwr("HMCparameters.xml"); write(HMCwr,"HMCparameters",HMCparams); } + + using Tst1=Hadrons::NamedTensor; + using Tst2=Hadrons::NamedTensor; + using Tst3=Hadrons::NamedTensor; + std::cout << "typeid(Tst1)=" << typeid(Tst1).name() << std::endl; + std::cout << "typeid(Tst2)=" << typeid(Tst2).name() << std::endl; + std::cout << "typeid(Tst3)=" << typeid(Tst3).name() << std::endl; + std::cout << "Tst1 is " << ( typeid(Tst1) == typeid(Tst2) ? "" : "not " ) + << "the same as Tst2" << std::endl; + std::cout << "Tst3 is " << ( typeid(Tst3) == typeid(Tst2) ? "" : "not " ) + << "the same as Tst2" << std::endl; + std::cout << "typeid(PerambTensor)=" << typeid(Grid::Hadrons::MDistil::PerambTensor).name() << std::endl; + std::cout << "typeid(MyPerambTensor)=" << typeid(MyPerambTensor).name() << std::endl; + std::cout << "Grid::Hadrons::MDistil::PerambTensor is " + << ( typeid(Grid::Hadrons::MDistil::PerambTensor) == typeid(MyPerambTensor) ? "" : "not " ) + << "the same as MyPerambTensor" << std::endl; + std::cout << "Grid::Hadrons::MDistil::PerambTensor is " + << ( typeid(Grid::Hadrons::MDistil::PerambTensor) == typeid(MyPerambTensor2) ? "" : "not " ) + << "the same as MyPerambTensor2" << std::endl; + + const std::string FileOriginal{"Peramb.3000"}; + const std::string FileCopy{"Peramb_deleteme.3000"}; + MyPerambTensor p; + p.read(FileOriginal,true,std::string("NamedTensor")); + //p.IndexNames[3] = "Lemons " + p.IndexNames[3]; + p.write(FileCopy); + p.read(FileCopy,false); + p.read(FileCopy); Grid_finalize(); }