mirror of
https://github.com/paboyle/Grid.git
synced 2024-09-20 01:05:38 +01:00
Rationalisation of NamedTensor (Perambulator)
This commit is contained in:
Michael Marshall
parent
52d8d576d0
commit
1c10933db1
@ -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 <ferben@ed.ac.uk>
|
||||
Author: Michael Marshall <Michael.Marshall@ed.ac.uk>
|
||||
|
||||
|
||||
This program is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
@ -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 <Hadrons/Global.hpp>
|
||||
#include <Hadrons/Module.hpp>
|
||||
#include <Hadrons/ModuleFactory.hpp>
|
||||
#include <Hadrons/Solver.hpp>
|
||||
#include <Hadrons/EigenPack.hpp>
|
||||
#include <Hadrons/A2AVectors.hpp>
|
||||
#include <Hadrons/DilutedNoise.hpp>
|
||||
|
||||
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<typename Scalar_, int NumIndices_, const std::string &Name_, const std::array<std::string,NumIndices_> &IndexNames_>
|
||||
class NamedTensor : Serializable
|
||||
{
|
||||
public:
|
||||
using Scalar = Scalar_;
|
||||
static constexpr int NumIndices = NumIndices_;
|
||||
using ET = Eigen::Tensor<Scalar_, NumIndices_, Eigen::RowMajor>;
|
||||
using Index = typename ET::Index;
|
||||
GRID_SERIALIZABLE_CLASS_MEMBERS(NamedTensor,
|
||||
ET, tensor,
|
||||
std::vector<std::string>, IndexNames );
|
||||
|
||||
// Get the default index names as std::vector
|
||||
std::vector<std::string> DefaultIndexNames()
|
||||
{
|
||||
std::vector<std::string> 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<typename... IndexTypes>
|
||||
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<typename Writer> 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<typename Reader> void read(Reader &r, bool bValidate = true, const std::string &Tag = Name_)
|
||||
{
|
||||
// Grab index names and dimensions
|
||||
std::vector<std::string> 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<LatticeColourVector>;
|
||||
|
||||
// Noise vector index order: nnoise, nt, nvec, ns
|
||||
// Noise vector (index order: nnoise, nt, nvec, ns)
|
||||
using NoiseTensor = Eigen::Tensor<Complex, 4, Eigen::RowMajor>;
|
||||
|
||||
struct DistilParameters: Serializable {
|
||||
GRID_SERIALIZABLE_CLASS_MEMBERS(DistilParameters,
|
||||
int, nnoise,
|
||||
int, tsrc,
|
||||
std::string, TI,
|
||||
std::string, LI,
|
||||
std::string, SI )
|
||||
DistilParameters() = default;
|
||||
template <class ReaderClass> DistilParameters(Reader<ReaderClass>& 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<typename Scalar_, int NumIndices_>
|
||||
class NamedTensor : Serializable
|
||||
{
|
||||
public:
|
||||
using Scalar = Scalar_;
|
||||
static constexpr int NumIndices = NumIndices_;
|
||||
using ET = Eigen::Tensor<Scalar_, NumIndices_, Eigen::RowMajor>;
|
||||
using Index = typename ET::Index;
|
||||
GRID_SERIALIZABLE_CLASS_MEMBERS(NamedTensor
|
||||
, ET, tensor
|
||||
, std::vector<std::string>, IndexNames
|
||||
);
|
||||
// Named tensors are intended to be a superset of Eigen tensor
|
||||
inline operator ET&() { return tensor; }
|
||||
template<typename... IndexTypes>
|
||||
inline const Scalar_& operator()(const std::array<Eigen::Index, NumIndices_> &Indices) const
|
||||
{ return tensor.operator()(Indices); }
|
||||
inline Scalar_& operator()(const std::array<Eigen::Index, NumIndices_> &Indices)
|
||||
{ return tensor.operator()(Indices); }
|
||||
template<typename... IndexTypes>
|
||||
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<Eigen::Index, NumIndices_>{{firstDimension, otherDimensions...}});
|
||||
}
|
||||
template<typename... IndexTypes>
|
||||
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<Eigen::Index, NumIndices_>{{firstDimension, otherDimensions...}});
|
||||
}
|
||||
|
||||
// Construct a named tensor explicitly specifying size of each dimension
|
||||
template<typename... IndexTypes>
|
||||
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NamedTensor(const std::array<std::string,NumIndices_> &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<std::string,NumIndices_> &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<typename Reader> inline void read (Reader &r, const char * pszTag = nullptr);
|
||||
template<typename Writer> 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<typename T, typename V = void> struct is_named_tensor : public std::false_type {};
|
||||
template<typename Scalar_, int NumIndices_> struct is_named_tensor<NamedTensor<Scalar_, NumIndices_>> : public std::true_type {};
|
||||
template<typename T> struct is_named_tensor<T, typename std::enable_if<std::is_base_of<NamedTensor<typename T::Scalar, T::NumIndices>, T>::value>::type> : public std::true_type {};
|
||||
|
||||
/******************************************************************************
|
||||
PerambTensor object
|
||||
******************************************************************************/
|
||||
|
||||
using PerambTensor = NamedTensor<SpinVector, 6>;
|
||||
static const std::array<std::string, 6> PerambIndexNames{"nT", "nVec", "LI", "nNoise", "nT_inv", "SI"};
|
||||
|
||||
/******************************************************************************
|
||||
Write NamedTensor
|
||||
******************************************************************************/
|
||||
|
||||
template<typename Scalar_, int NumIndices_>
|
||||
template<typename Writer>
|
||||
void NamedTensor<Scalar_, NumIndices_>::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<typename Scalar_, int NumIndices_>
|
||||
void NamedTensor<Scalar_, NumIndices_>::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<typename Scalar_, int NumIndices_>
|
||||
bool NamedTensor<Scalar_, NumIndices_>::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<typename Scalar_, int NumIndices_>
|
||||
template<typename Reader>
|
||||
void NamedTensor<Scalar_, NumIndices_>::read(Reader &r, const char * pszTag) {
|
||||
if( pszTag == nullptr )
|
||||
pszTag = "NamedTensor";
|
||||
// Grab index names and dimensions
|
||||
std::vector<std::string> 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<typename Scalar_, int NumIndices_>
|
||||
void NamedTensor<Scalar_, NumIndices_>::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<int>(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<LatticeColourVector> & 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<std::string, 6> PerambIndexNames;
|
||||
using PerambTensor = NamedTensor<SpinVector, 6, PerambTensorName, PerambIndexNames>;
|
||||
|
||||
END_MODULE_NAMESPACE
|
||||
END_HADRONS_NAMESPACE
|
||||
#endif // Hadrons_MDistil_Distil_hpp_
|
||||
#endif // Hadrons_Distil_hpp_
|
||||
|
||||
143
Hadrons/Modules/MDistil/DistilCommon.hpp
Normal file
143
Hadrons/Modules/MDistil/DistilCommon.hpp
Normal file
@ -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 <ferben@ed.ac.uk>
|
||||
Author: Michael Marshall <Michael.Marshall@ed.ac.uk>
|
||||
|
||||
This program is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License along
|
||||
with this program; if not, write to the Free Software Foundation, Inc.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
|
||||
See the full license in the file "LICENSE" in the top level distribution directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
|
||||
#ifndef Hadrons_MDistil_DistilCommon_hpp_
|
||||
#define Hadrons_MDistil_DistilCommon_hpp_
|
||||
|
||||
#include <Hadrons/Distil.hpp>
|
||||
#include <Hadrons/Module.hpp>
|
||||
#include <Hadrons/ModuleFactory.hpp>
|
||||
#include <Hadrons/Solver.hpp>
|
||||
#include <Hadrons/A2AVectors.hpp>
|
||||
#include <Hadrons/DilutedNoise.hpp>
|
||||
|
||||
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 <class ReaderClass> DistilParameters(Reader<ReaderClass>& 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<int>(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<LatticeColourVector> & 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_
|
||||
@ -30,24 +30,15 @@
|
||||
#ifndef Hadrons_MDistil_DistilVectors_hpp_
|
||||
#define Hadrons_MDistil_DistilVectors_hpp_
|
||||
|
||||
#include <Hadrons/Global.hpp>
|
||||
#include <Hadrons/Module.hpp>
|
||||
#include <Hadrons/ModuleFactory.hpp>
|
||||
#include <Hadrons/Solver.hpp>
|
||||
#include <Hadrons/EigenPack.hpp>
|
||||
#include <Hadrons/A2AVectors.hpp>
|
||||
#include <Hadrons/DilutedNoise.hpp>
|
||||
|
||||
// These are members of Distillation
|
||||
#include <Hadrons/Distil.hpp>
|
||||
#include <Hadrons/Modules/MDistil/DistilCommon.hpp>
|
||||
|
||||
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<FImpl>::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<int>( perambulator.tensor.dimension(0) ) && "PerambTensor time dimensionality bad" );
|
||||
@ -287,7 +278,7 @@ void TDistilVectors<FImpl>::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);
|
||||
}
|
||||
|
||||
@ -30,16 +30,9 @@
|
||||
#ifndef Hadrons_MDistil_LapEvec_hpp_
|
||||
#define Hadrons_MDistil_LapEvec_hpp_
|
||||
|
||||
#include <Hadrons/Global.hpp>
|
||||
#include <Hadrons/Module.hpp>
|
||||
#include <Hadrons/ModuleFactory.hpp>
|
||||
#include <Hadrons/EigenPack.hpp>
|
||||
|
||||
// These are members of Distillation
|
||||
#include <Hadrons/Distil.hpp>
|
||||
#include <Hadrons/Modules/MDistil/DistilCommon.hpp>
|
||||
|
||||
BEGIN_HADRONS_NAMESPACE
|
||||
|
||||
BEGIN_MODULE_NAMESPACE(MDistil)
|
||||
|
||||
/******************************************************************************
|
||||
|
||||
@ -30,23 +30,14 @@
|
||||
#ifndef Hadrons_MDistil_Noises_hpp_
|
||||
#define Hadrons_MDistil_Noises_hpp_
|
||||
|
||||
#include <Hadrons/Global.hpp>
|
||||
#include <Hadrons/Module.hpp>
|
||||
#include <Hadrons/ModuleFactory.hpp>
|
||||
#include <Hadrons/Solver.hpp>
|
||||
#include <Hadrons/EigenPack.hpp>
|
||||
#include <Hadrons/A2AVectors.hpp>
|
||||
#include <Hadrons/DilutedNoise.hpp>
|
||||
|
||||
// These are members of Distillation
|
||||
#include <Hadrons/Distil.hpp>
|
||||
#include <Hadrons/Modules/MDistil/DistilCommon.hpp>
|
||||
|
||||
BEGIN_HADRONS_NAMESPACE
|
||||
BEGIN_MODULE_NAMESPACE(MDistil)
|
||||
|
||||
/******************************************************************************
|
||||
* Noises *
|
||||
******************************************************************************/
|
||||
BEGIN_MODULE_NAMESPACE(MDistil)
|
||||
|
||||
class NoisesPar: Serializable
|
||||
{
|
||||
|
||||
@ -30,23 +30,14 @@
|
||||
#ifndef Hadrons_MDistil_PerambFromSolve_hpp_
|
||||
#define Hadrons_MDistil_PerambFromSolve_hpp_
|
||||
|
||||
#include <Hadrons/Global.hpp>
|
||||
#include <Hadrons/Module.hpp>
|
||||
#include <Hadrons/ModuleFactory.hpp>
|
||||
#include <Hadrons/Solver.hpp>
|
||||
#include <Hadrons/EigenPack.hpp>
|
||||
#include <Hadrons/A2AVectors.hpp>
|
||||
#include <Hadrons/DilutedNoise.hpp>
|
||||
#include <Hadrons/Modules/MDistil/DistilCommon.hpp>
|
||||
|
||||
// These are members of Distillation
|
||||
#include <Hadrons/Distil.hpp>
|
||||
|
||||
BEGIN_HADRONS_NAMESPACE
|
||||
BEGIN_MODULE_NAMESPACE(MDistil)
|
||||
|
||||
/******************************************************************************
|
||||
* PerambFromSolve *
|
||||
******************************************************************************/
|
||||
BEGIN_MODULE_NAMESPACE(MDistil)
|
||||
|
||||
class PerambFromSolvePar: Serializable
|
||||
{
|
||||
@ -126,7 +117,7 @@ void TPerambFromSolve<FImpl>::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<FImpl>::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<Complex>(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<Complex>(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;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -34,3 +34,14 @@ using namespace Hadrons;
|
||||
using namespace MDistil;
|
||||
|
||||
template class Grid::Hadrons::MDistil::TPerambulator<FIMPL>;
|
||||
|
||||
// Global constants for distillation
|
||||
|
||||
const std::string Grid::Hadrons::MDistil::PerambTensorName{ "Perambulator" };
|
||||
const std::array<std::string, 6> 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
|
||||
|
||||
@ -30,7 +30,7 @@
|
||||
#ifndef Hadrons_MDistil_Perambulator_hpp_
|
||||
#define Hadrons_MDistil_Perambulator_hpp_
|
||||
|
||||
#include <Hadrons/Distil.hpp>
|
||||
#include <Hadrons/Modules/MDistil/DistilCommon.hpp>
|
||||
|
||||
BEGIN_HADRONS_NAMESPACE
|
||||
BEGIN_MODULE_NAMESPACE(MDistil)
|
||||
@ -129,7 +129,7 @@ void TPerambulator<FImpl>::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<FermionField>, getName() + "_unsmeared_sink", 1,
|
||||
nnoise*LI*Ns*Nt_inv, envGetGrid(FermionField));
|
||||
|
||||
@ -244,7 +244,7 @@ void TPerambulator<FImpl>::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<Complex>(innerProduct(evec3d, result_3d)),is);
|
||||
pokeSpin(perambulator.tensor(t, ivec, dk, inoise,dt,ds),static_cast<Complex>(innerProduct(evec3d, result_3d)),is);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -30,18 +30,14 @@
|
||||
#ifndef Hadrons_MIO_LoadPerambulator_hpp_
|
||||
#define Hadrons_MIO_LoadPerambulator_hpp_
|
||||
|
||||
#include <Hadrons/Global.hpp>
|
||||
#include <Hadrons/Module.hpp>
|
||||
#include <Hadrons/ModuleFactory.hpp>
|
||||
#include <Hadrons/EigenPack.hpp>
|
||||
#include <Hadrons/Distil.hpp>
|
||||
#include <Hadrons/Modules/MDistil/DistilCommon.hpp>
|
||||
|
||||
BEGIN_HADRONS_NAMESPACE
|
||||
BEGIN_MODULE_NAMESPACE(MIO)
|
||||
|
||||
/******************************************************************************
|
||||
* LoadPerambulator *
|
||||
******************************************************************************/
|
||||
BEGIN_MODULE_NAMESPACE(MIO)
|
||||
|
||||
class LoadPerambulatorPar: Serializable
|
||||
{
|
||||
@ -103,7 +99,7 @@ void TLoadPerambulator<FImpl>::setup(void)
|
||||
{
|
||||
DISTIL_PARAMETERS_DEFINE( true );
|
||||
//std::array<std::string,6> 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 ///////////////////////////////////////////////////////////////////
|
||||
|
||||
@ -30,6 +30,7 @@ Author: Michael Marshall <michael.marshall@ed.ac.uk>
|
||||
/* END LEGAL */
|
||||
#include <Grid/Grid.h>
|
||||
#include <typeinfo>
|
||||
#include <Hadrons/Distil.hpp>
|
||||
|
||||
using namespace Grid;
|
||||
|
||||
@ -252,6 +253,14 @@ void tensorConvTestFn(GridSerialRNG &rng, const std::string label)
|
||||
}
|
||||
|
||||
#define tensorConvTest(rng, type) tensorConvTestFn<type>(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<std::string,2> GreatestShowOnEarth{ "Burnum", "Burnum" };
|
||||
|
||||
using MyPerambTensor = Grid::Hadrons::NamedTensor<SpinVector, 6, Grid::Hadrons::MDistil::PerambTensorName, Grid::Hadrons::MDistil::PerambIndexNames>;
|
||||
static const std::string MPTName{ "Perambulator" };
|
||||
static const std::array<std::string, 6> MPTIndexNames{"nT", "nVec", "LI", "nNoise", "nT_inv", "SI"};
|
||||
using MyPerambTensor2 = Grid::Hadrons::NamedTensor<SpinVector, 6, MPTName, MPTIndexNames>;
|
||||
|
||||
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<double,2,Coco,GreatestShowOnEarth>;
|
||||
using Tst2=Hadrons::NamedTensor<double,2,TheyCallMeLurve,GreatestShowOnEarth>;
|
||||
using Tst3=Hadrons::NamedTensor<double,2,TheyCallMeLurve,GreatestShowOnEarth>;
|
||||
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();
|
||||
}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user