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

Fairly close to ready for release. Felix and I to review, then submit for release

Browse Source
This commit is contained in:
Michael Marshall 2019-04-30 23:53:57 +01:00
parent e56ead55ef
commit 334f29becb
14 changed files with 290 additions and 1036 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

62
Hadrons/Distil.hpp
View File

@ -33,16 +33,10 @@
#include <Hadrons/Global.hpp>
#include <Hadrons/Module.hpp>
#include <Hadrons/ModuleFactory.hpp>
/******************************************************************************
Needed to make sure envCreate() (see Hadrons) work with specialisations
with more than one parameter, eg obj<T1 COMMA T2>
I imagine this exists already?
******************************************************************************/
#ifndef COMMA
#define COMMA ,
#endif
#include <Hadrons/Solver.hpp>
#include <Hadrons/EigenPack.hpp>
#include <Hadrons/A2AVectors.hpp>
#include <Hadrons/DilutedNoise.hpp>
/******************************************************************************
A consistent set of cross-platform methods for big endian <-> host byte ordering
@ -207,24 +201,53 @@ BEGIN_HADRONS_NAMESPACE
BEGIN_MODULE_NAMESPACE(MDistil)
using DistilEP = Grid::Hadrons::EigenPack<LatticeColourVector>;
using LapEvecs = Grid::Hadrons::EigenPack<LatticeColourVector>;
// 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, TI,
int, LI,
int, nnoise,
int, tsrc,
int, SI,
int, Ns,
int, Nt_inv)
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;
}
int getTI( const Environment & env, bool bCalledFromSetup ) const {
return ParameterDefault( TI, env.getDim(Tdir), bCalledFromSetup ); }
};
#define DISTIL_PARAMETERS_DEFINE( inSetup ) \
const int Nt{env().getDim(Tdir)}; \
const int nvec{par().nvec}; \
const int Ns{Grid::QCD::Ns}; \
const int nnoise{par().Distil.nnoise}; \
const int tsrc{par().Distil.tsrc}; \
const int TI{par().Distil.getTI(env(), 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 }
class BFieldIO: Serializable{
public:
using BaryonTensorSet = Eigen::Tensor<Complex, 7>;
@ -299,7 +322,7 @@ public:
// Construct a named tensor explicitly specifying size of each dimension
template<typename... IndexTypes>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NamedTensor(std::array<std::string,NumIndices_> &IndexNames_, Eigen::Index firstDimension, IndexTypes... otherDimensions)
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.
@ -343,8 +366,9 @@ template<typename T> struct is_named_tensor<T, typename std::enable_if<std::is_b
Endian_Scalar_Size can be removed once (deprecated) NamedTensor::ReadBinary & WriteBinary methods deleted
******************************************************************************/
template<typename Scalar_, int NumIndices_, uint16_t Endian_Scalar_Size = sizeof(Scalar_)>
using Perambulator = NamedTensor<Scalar_, NumIndices_, Endian_Scalar_Size>;
//template<typename Scalar_, int NumIndices_, uint16_t Endian_Scalar_Size = sizeof(Scalar_)>
using Perambulator = NamedTensor<SpinVector, 6, sizeof(Real)>;
static const std::array<std::string, 6> PerambIndexNames{"Nt", "nvec", "LI", "nnoise", "Nt_inv", "SI"};
/******************************************************************************
Save NamedTensor binary format (NB: On-disk format is Big Endian)

2
Hadrons/Modules.hpp
View File

@ -39,10 +39,8 @@
#include <Hadrons/Modules/MScalar/FreeProp.hpp>
#include <Hadrons/Modules/MScalar/Scalar.hpp>
#include <Hadrons/Modules/MScalar/ChargedProp.hpp>
#include <Hadrons/Modules/MDistil/DistilSource.hpp>
#include <Hadrons/Modules/MDistil/LapEvec.hpp>
#include <Hadrons/Modules/MDistil/DistilVectors.hpp>
#include <Hadrons/Modules/MDistil/DistilSink.hpp>
#include <Hadrons/Modules/MDistil/Noises.hpp>
#include <Hadrons/Modules/MDistil/BC2.hpp>
#include <Hadrons/Modules/MDistil/Perambulator.hpp>

36
Hadrons/Modules/MDistil/DistilSink.cc
View File

@ -1,36 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MDistil/DistilSink.cc
Copyright (C) 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 */
#include <Hadrons/Modules/MDistil/DistilSink.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MDistil;
template class Grid::Hadrons::MDistil::TDistilSink<FIMPL>;

236
Hadrons/Modules/MDistil/DistilSink.hpp
View File

@ -1,236 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MDistil/DistilSink.hpp
Copyright (C) 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_DistilSink_hpp_
#define Hadrons_MDistil_DistilSink_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>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* DistilSink *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MDistil)
class DistilSinkPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(DistilSinkPar,
std::string, perambulator,
std::string, eigenPack,
bool, multiFile,
int, tsrc,
int, nnoise,
int, LI,
int, SI,
int, TI,
int, nvec,
int, Ns,
int, Nt,
int, Nt_inv);
};
template <typename FImpl>
class TDistilSink: public Module<DistilSinkPar>
{
public:
FERM_TYPE_ALIASES(FImpl,);
public:
// constructor
TDistilSink(const std::string name);
// destructor
virtual ~TDistilSink(void);
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
protected:
// These variables are created in setup() and freed in Cleanup()
GridCartesian * grid3d; // Owned by me, so I must delete it
GridCartesian * grid4d; // Owned by environment (so I won't delete it)
protected:
virtual void Cleanup(void);
};
MODULE_REGISTER_TMP(DistilSink, TDistilSink<FIMPL>, MDistil);
/******************************************************************************
* TDistilSink implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TDistilSink<FImpl>::TDistilSink(const std::string name)
: grid3d{nullptr}, grid4d{nullptr}, Module<DistilSinkPar>(name)
{}
// destructor
template <typename FImpl>
TDistilSink<FImpl>::~TDistilSink(void)
{
Cleanup();
};
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TDistilSink<FImpl>::getInput(void)
{
std::vector<std::string> in;
in.push_back(par().perambulator);
in.push_back(par().eigenPack);
return in;
}
template <typename FImpl>
std::vector<std::string> TDistilSink<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TDistilSink<FImpl>::setup(void)
{
Cleanup();
int nnoise=par().nnoise;
int LI=par().LI;
int Ns=par().Ns;
int SI=par().SI;
int Nt_inv=par().Nt_inv;
envCreate(std::vector<FermionField>, getName(), 1,
nnoise*LI*SI*Nt_inv, envGetGrid(FermionField));
grid4d = env().getGrid();
std::vector<int> latt_size = GridDefaultLatt();
std::vector<int> simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
std::vector<int> mpi_layout = GridDefaultMpi();
std::vector<int> simd_layout_3 = GridDefaultSimd(Nd-1, vComplex::Nsimd());
latt_size[Nd-1] = 1;
simd_layout_3.push_back( 1 );
mpi_layout[Nd-1] = 1;
grid3d = MakeLowerDimGrid(grid4d);
envTmp(LatticeSpinColourVector, "tmp2",1,LatticeSpinColourVector(grid4d));
envTmp(LatticeColourVector, "tmp_nospin",1,LatticeColourVector(grid4d));
envTmp(LatticeSpinColourVector, "tmp3d",1,LatticeSpinColourVector(grid3d));
envTmp(LatticeColourVector, "tmp3d_nospin",1,LatticeColourVector(grid3d));
envTmp(LatticeSpinColourVector, "sink_tslice",1,LatticeSpinColourVector(grid3d));
envTmp(LatticeColourVector, "evec3d",1,LatticeColourVector(grid3d));
}
// clean up any temporaries created by setup (that aren't stored in the environment)
template <typename FImpl>
void TDistilSink<FImpl>::Cleanup(void)
{
if( grid3d != nullptr ) {
delete grid3d;
grid3d = nullptr;
}
grid4d = nullptr;
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TDistilSink<FImpl>::execute(void)
{
auto &perambulator = envGet(Perambulator<SpinVector COMMA 6 COMMA sizeof(Real)>, par().perambulator);
auto &epack = envGet(Grid::Hadrons::EigenPack<LatticeColourVector>, par().eigenPack);
auto &phi = envGet(std::vector<FermionField>, getName());
envGetTmp(LatticeSpinColourVector, tmp2);
envGetTmp(LatticeColourVector, tmp_nospin);
envGetTmp(LatticeSpinColourVector, tmp3d);
envGetTmp(LatticeColourVector, tmp3d_nospin);
envGetTmp(LatticeSpinColourVector, sink_tslice);
envGetTmp(LatticeColourVector, evec3d);
int Ntlocal = grid4d->LocalDimensions()[3];
int Ntfirst = grid4d->LocalStarts()[3];
int tsrc=par().tsrc;
int nnoise=par().nnoise;
int LI=par().LI;
int SI=par().SI;
int Ns=par().Ns;
int Nt_inv=par().Nt_inv; // TODO: No input, but define through Nt, TI
int Nt=par().Nt;
int TI=par().TI;
int nvec=par().nvec;
bool full_tdil=(TI==Nt);
int vecindex;
int t_inv;
for (int inoise = 0; inoise < nnoise; inoise++) {
for (int dk = 0; dk < LI; dk++) {
for (int dt = 0; dt < Nt_inv; dt++) {
for (int ds = 0; ds < SI; ds++) {
vecindex = inoise + nnoise * dk + nnoise * LI * ds + nnoise *LI * SI*dt;
phi[vecindex] = zero;
for (int t = Ntfirst; t < Ntfirst + Ntlocal; t++) {
sink_tslice=zero;
for (int ivec = 0; ivec < nvec; ivec++) {
ExtractSliceLocal(evec3d,epack.evec[ivec],0,t,3);
sink_tslice += evec3d * perambulator(t, ivec, dk, inoise,dt,ds);
}
InsertSliceLocal(sink_tslice,phi[vecindex],0,t-Ntfirst,Grid::QCD::Tdir);
}
}
}
}
}
// TEST TO SEE WHETHER THIS MIGHT BE THE MEMORY LEAK
Cleanup();
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MDistil_DistilSink_hpp_

36
Hadrons/Modules/MDistil/DistilSource.cc
View File

@ -1,36 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MDistil/DistilSource.cc
Copyright (C) 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 */
#include <Hadrons/Modules/MDistil/DistilSource.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MDistil;
template class Grid::Hadrons::MDistil::TDistilSource<FIMPL>;

250
Hadrons/Modules/MDistil/DistilSource.hpp
View File

@ -1,250 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MDistil/DistilSource.hpp
Copyright (C) 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_DistilSource_hpp_
#define Hadrons_MDistil_DistilSource_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>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* DistilSource *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MDistil)
class DistilSourcePar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(DistilSourcePar,
std::string, noise,
std::string, eigenPack,
bool, multiFile,
int, tsrc,
int, nnoise,
int, LI,
int, SI,
int, TI,
int, nvec,
int, Ns,
int, Nt,
int, Nt_inv);
};
template <typename FImpl>
class TDistilSource: public Module<DistilSourcePar>
{
public:
FERM_TYPE_ALIASES(FImpl,);
public:
// constructor
TDistilSource(const std::string name);
// destructor
virtual ~TDistilSource(void);
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
protected:
// These variables are created in setup() and freed in Cleanup()
GridCartesian * grid3d; // Owned by me, so I must delete it
GridCartesian * grid4d; // Owned by environment (so I won't delete it)
protected:
virtual void Cleanup(void);
};
MODULE_REGISTER_TMP(DistilSource, TDistilSource<FIMPL>, MDistil);
/******************************************************************************
* TDistilSource implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TDistilSource<FImpl>::TDistilSource(const std::string name)
: grid3d{nullptr}, grid4d{nullptr}, Module<DistilSourcePar>(name)
{}
// destructor
template <typename FImpl>
TDistilSource<FImpl>::~TDistilSource(void)
{
Cleanup();
};
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TDistilSource<FImpl>::getInput(void)
{
std::vector<std::string> in;
in.push_back(par().noise);
in.push_back(par().eigenPack);
return in;
}
template <typename FImpl>
std::vector<std::string> TDistilSource<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TDistilSource<FImpl>::setup(void)
{
Cleanup();
auto &noise = envGet(std::vector<Complex>, par().noise);
int nnoise=par().nnoise;
int LI=par().LI;
int Ns=par().Ns;
int SI=par().SI;
int Nt_inv=par().Nt_inv;
envCreate(std::vector<FermionField>, getName(), 1,
nnoise*LI*SI*Nt_inv, envGetGrid(FermionField));
grid4d = env().getGrid();
std::vector<int> latt_size = GridDefaultLatt();
std::vector<int> simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
std::vector<int> mpi_layout = GridDefaultMpi();
std::vector<int> simd_layout_3 = GridDefaultSimd(Nd-1, vComplex::Nsimd());
latt_size[Nd-1] = 1;
simd_layout_3.push_back( 1 );
mpi_layout[Nd-1] = 1;
grid3d = MakeLowerDimGrid(grid4d);
envTmp(LatticeSpinColourVector, "tmp2",1,LatticeSpinColourVector(grid4d));
envTmp(LatticeColourVector, "tmp_nospin",1,LatticeColourVector(grid4d));
envTmp(LatticeSpinColourVector, "tmp3d",1,LatticeSpinColourVector(grid3d));
envTmp(LatticeColourVector, "tmp3d_nospin",1,LatticeColourVector(grid3d));
envTmp(LatticeSpinColourVector, "sink_tslice",1,LatticeSpinColourVector(grid3d));
envTmp(LatticeColourVector, "evec3d",1,LatticeColourVector(grid3d));
}
// clean up any temporaries created by setup (that aren't stored in the environment)
template <typename FImpl>
void TDistilSource<FImpl>::Cleanup(void)
{
if( grid3d != nullptr ) {
delete grid3d;
grid3d = nullptr;
}
grid4d = nullptr;
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TDistilSource<FImpl>::execute(void)
{
auto &noise = envGet(std::vector<Complex>, par().noise);
auto &epack = envGet(Grid::Hadrons::EigenPack<LatticeColourVector>, par().eigenPack);
auto &rho = envGet(std::vector<FermionField>, getName());
envGetTmp(LatticeSpinColourVector, tmp2);
envGetTmp(LatticeColourVector, tmp_nospin);
envGetTmp(LatticeSpinColourVector, tmp3d);
envGetTmp(LatticeColourVector, tmp3d_nospin);
envGetTmp(LatticeSpinColourVector, sink_tslice);
envGetTmp(LatticeColourVector, evec3d);
int Ntlocal = grid4d->LocalDimensions()[3];
int Ntfirst = grid4d->LocalStarts()[3];
int tsrc=par().tsrc;
int nnoise=par().nnoise;
int LI=par().LI;
int Ns=par().Ns;
int Nt_inv=par().Nt_inv; // TODO: No input, but define through Nt, TI
int Nt=par().Nt;
int TI=par().TI;
int nvec=par().nvec;
int SI=par().SI;
bool full_tdil=(TI==Nt);
int vecindex;
int t_inv;
for (int inoise = 0; inoise < nnoise; inoise++) {
for (int dk = 0; dk < LI; dk++) {
for (int dt = 0; dt < Nt_inv; dt++) {
for (int ds = 0; ds < SI; ds++) {
vecindex = inoise + nnoise * dk + nnoise * LI * ds + nnoise *LI * SI*dt;
rho[vecindex] = zero;
tmp3d_nospin = zero;
for (int it = dt; it < Nt; it += TI){
if (full_tdil) t_inv = tsrc; else t_inv = it;
if( t_inv >= Ntfirst && t_inv < Ntfirst + Ntlocal ) {
for (int ik = dk; ik < nvec; ik += LI){
for (int is = ds; is < Ns; is += SI){
ExtractSliceLocal(evec3d,epack.evec[ik],0,t_inv,3);
tmp3d_nospin = evec3d * noise[inoise + nnoise*(t_inv + Nt*(ik+nvec*is))];
tmp3d=zero;
pokeSpin(tmp3d,tmp3d_nospin,is);
tmp2=zero;
InsertSliceLocal(tmp3d,tmp2,0,t_inv-Ntfirst,Grid::QCD::Tdir);
rho[vecindex] += tmp2;
}
}
}
}
}
}
}
}
// TEST TO SEE WHETHER THIS MIGHT BE THE MEMORY LEAK
Cleanup();
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MDistil_DistilSource_hpp_

75
Hadrons/Modules/MDistil/DistilVectors.hpp
View File

@ -60,12 +60,8 @@ public:
std::string, sink,
bool, multiFile,
int, tsrc,
int, LI,
int, SI,
int, TI,
int, nvec,
int, Ns,
int, Nt_inv);
std::string, nvec,
std::string, TI)
};
template <typename FImpl>
@ -73,8 +69,6 @@ class TDistilVectors: public Module<DistilVectorsPar>
{
public:
FERM_TYPE_ALIASES(FImpl,);
// This is the type of perambulator I expect
using DistilPeramb = Perambulator<SpinVector, 6, sizeof(Real)>;
// constructor
TDistilVectors(const std::string name);
// destructor
@ -100,7 +94,6 @@ public:
bool bMakeSink;
std::string SourceName;
std::string SinkName;
int nnoise;
};
MODULE_REGISTER_TMP(DistilVectors, TDistilVectors<FIMPL>, MDistil);
@ -171,19 +164,24 @@ void TDistilVectors<FImpl>::setup(void)
{
Cleanup();
auto &noise = envGet(NoiseTensor, NoiseVectorName);
auto &perambulator = envGet(DistilPeramb, PerambulatorName);
auto &perambulator = envGet(Perambulator, PerambulatorName);
// We expect the perambulator to have been created with these indices
std::array<std::string,6> sIndexNames{"Nt", "nvec", "LI", "nnoise", "Nt_inv", "SI"};
for(int i = 0; i < DistilPeramb::NumIndices; i++ )
assert( sIndexNames[i] == perambulator.IndexNames[i] && "Perambulator indices bad" );
for(int i = 0; i < Perambulator::NumIndices; i++ )
assert( PerambIndexNames[i] == perambulator.IndexNames[i] && "Perambulator indices bad" );
nnoise = static_cast<int>( noise.dimension(0) );
LOG(Message) << "NoiseTensor has " << nnoise << " noise vectors" << std::endl;
int LI=par().LI;
int Ns=par().Ns;
int SI=par().SI;
int Nt_inv=par().Nt_inv;
const int Nt{ env().getDim(Tdir) };
assert( Nt == static_cast<int>( perambulator.tensor.dimension(0) ) && "Perambulator time dimensionality bad" );
const int TI{ Hadrons::MDistil::DistilParameters::ParameterDefault( par().TI, Nt, true) };
const int LI{ static_cast<int>( perambulator.tensor.dimension(2) ) };
const int SI{ static_cast<int>( perambulator.tensor.dimension(5) ) };
const int Nt_inv{ static_cast<int>( perambulator.tensor.dimension(4) ) };
const int nnoise{ static_cast<int>( perambulator.tensor.dimension(3) ) };
assert( nnoise >= static_cast<int>( noise.dimension(0) ) && "Not enough noise vectors for perambulator" );
// Nvec defaults to what's in the perambulator unless overriden
const int nvec_per{ static_cast<int>( perambulator.tensor.dimension(1) ) };
const int nvec{Hadrons::MDistil::DistilParameters::ParameterDefault(par().nvec, nvec_per, true) };
assert( nvec <= nvec_per && "Not enough distillation sub-space vectors" );
if( bMakeSource )
envCreate(std::vector<FermionField>, SourceName, 1, nnoise*LI*SI*Nt_inv, envGetGrid(FermionField));
@ -225,31 +223,30 @@ template <typename FImpl>
void TDistilVectors<FImpl>::execute(void)
{
auto &noise = envGet(NoiseTensor, NoiseVectorName);
auto &perambulator = envGet(DistilPeramb, PerambulatorName);
auto &perambulator = envGet(Perambulator, PerambulatorName);
auto &epack = envGet(Grid::Hadrons::EigenPack<LatticeColourVector>, LapEvecName);
envGetTmp(LatticeSpinColourVector, tmp2);
//envGetTmp(LatticeColourVector, tmp_nospin);
envGetTmp(LatticeSpinColourVector, tmp3d);
envGetTmp(LatticeColourVector, tmp3d_nospin);
envGetTmp(LatticeColourVector, tmp3d_nospin);
envGetTmp(LatticeSpinColourVector, sink_tslice);
envGetTmp(LatticeColourVector, evec3d);
int Ntlocal = grid4d->LocalDimensions()[3];
int Ntfirst = grid4d->LocalStarts()[3];
int tsrc=par().tsrc;
int LI=par().LI;
int Ns=par().Ns;
int Nt_inv=par().Nt_inv; // TODO: No input, but define through Nt, TI
const int Nt{grid4d->GlobalDimensions()[Tdir]};
int TI=par().TI;
int nvec=par().nvec;
int SI=par().SI;
bool full_tdil=(TI==Nt);
envGetTmp(LatticeColourVector, evec3d);
const int Ntlocal{ grid4d->LocalDimensions()[3] };
const int Ntfirst{ grid4d->LocalStarts()[3] };
const int Ns{ Grid::QCD::Ns };
const int Nt{ env().getDim(Tdir) };
const int TI{ Hadrons::MDistil::DistilParameters::ParameterDefault( par().TI, Nt, false ) };
const int LI{ static_cast<int>( perambulator.tensor.dimension(2) ) };
const int SI{ static_cast<int>( perambulator.tensor.dimension(5) ) };
const int Nt_inv{ static_cast<int>( perambulator.tensor.dimension(4) ) };
const int nnoise{ static_cast<int>( perambulator.tensor.dimension(3) ) };
// Nvec defaults to what's in the perambulator unless overriden
const int nvec{Hadrons::MDistil::DistilParameters::ParameterDefault(par().nvec, static_cast<int>( perambulator.tensor.dimension(1) ), false)};
const int tsrc{ par().tsrc };
const bool full_tdil{ TI==Nt };
int vecindex;
int t_inv;
if( bMakeSource ) {

77
Hadrons/Modules/MDistil/LapEvec.hpp
View File

@ -51,7 +51,7 @@ BEGIN_MODULE_NAMESPACE(MDistil)
struct StoutParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(StoutParameters,
int, steps,
double, parm) // TODO: change name of this to rho
double, rho) // TODO: change name of this to rho
StoutParameters() = default;
template <class ReaderClass> StoutParameters(Reader<ReaderClass>& Reader){read(Reader,"StoutSmearing",*this);}
};
@ -67,13 +67,10 @@ struct ChebyshevParameters: Serializable {
struct LanczosParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(LanczosParameters,
//int, Nstart,
int, Nvec,
int, Nk,
//int, Nm, // Not currently used
int, Np,
int, MaxIt,
//int, MinRes,
double, resid,
int, IRLLog)
LanczosParameters() = default;
@ -82,19 +79,12 @@ struct LanczosParameters: Serializable {
// These are the actual parameters passed to the module during construction
class LapEvecPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(LapEvecPar,
std::string, gauge,
// std::string, ConfigFileDir,
// std::string, ConfigFileName,
//,std::string, EigenPackName
StoutParameters, Stout
struct LapEvecPar: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(LapEvecPar
,std::string, gauge
,StoutParameters, Stout
,ChebyshevParameters, Cheby
,LanczosParameters, Lanczos
//,DistilParameters, Distil
)//,SolverParameters, Solver)
,LanczosParameters, Lanczos)
};
/******************************************************************************
@ -123,7 +113,7 @@ protected:
// These variables are created in setup() and freed in Cleanup()
GridCartesian * gridLD; // Owned by me, so I must delete it
GridCartesian * gridHD; // Owned by environment (so I won't delete it)
int Nx, Ny, Nz, Nt;
std::string sGaugeName;
protected:
virtual void Cleanup(void);
};
@ -151,8 +141,12 @@ TLapEvec<GImpl>::~TLapEvec()
template <typename GImpl>
std::vector<std::string> TLapEvec<GImpl>::getInput(void)
{
std::vector<std::string> in = {par().gauge};
return in;
sGaugeName = par().gauge;
if( sGaugeName.size() == 0 ) {
sGaugeName = getName();
sGaugeName.append( "_gauge" );
}
return std::vector<std::string>{ sGaugeName };
}
template <typename GImpl>
@ -170,19 +164,16 @@ void TLapEvec<GImpl>::setup(void)
Environment & e{env()};
gridHD = e.getGrid();
gridLD = MakeLowerDimGrid( gridHD );
Nx = gridHD->_fdimensions[Xdir];
Ny = gridHD->_fdimensions[Ydir];
Nz = gridHD->_fdimensions[Zdir];
Nt = gridHD->_fdimensions[Tdir];
const int Nt{e.getDim(Tdir)};
// Temporaries
//envTmpLat(GaugeField, "Umu");
envTmpLat(GaugeField, "Umu");
envTmpLat(GaugeField, "Umu_stout");
envTmpLat(GaugeField, "Umu_smear");
envTmp(LatticeGaugeField, "UmuNoTime",1,LatticeGaugeField(gridLD));
envTmp(LatticeColourVector, "src",1,LatticeColourVector(gridLD));
envTmp(std::vector<DistilEP>, "eig",1,std::vector<DistilEP>(Nt));
envTmp(std::vector<LapEvecs>, "eig",1,std::vector<LapEvecs>(Nt));
// Output objects
envCreate(DistilEP, getName(), 1, par().Lanczos.Nvec, gridHD );
envCreate(LapEvecs, getName(), 1, par().Lanczos.Nvec, gridHD );
}
// clean up any temporaries created by setup (that aren't stored in the environment)
@ -206,35 +197,21 @@ void TLapEvec<GImpl>::execute(void)
{
const ChebyshevParameters &ChebPar{par().Cheby};
const LanczosParameters &LPar{par().Lanczos};
const int &nvec{LPar.Nvec};
//const bool exact_distillation{TI==Nt && LI==nvec};
//const bool full_tdil{TI==Nt};
//const int &Nt_inv{full_tdil ? 1 : TI};
// Assertions on the parameters we read
//assert(TI>1);
//assert(LI>1);
//if(exact_distillation)
//assert(nnoise==1);
//else
//assert(nnoise>1);
// Disable IRL logging if requested
LOG(Message) << "IRLLog=" << LPar.IRLLog << std::endl;
const int PreviousIRLLogState{GridLogIRL.isActive()};
GridLogIRL.Active( LPar.IRLLog == 0 ? 0 : 1 );
auto &Umu = envGet(GaugeField, par().gauge);
envGetTmp(GaugeField, Umu_smear);
// Stout smearing
Umu_smear = Umu;
LOG(Message) << "Initial plaquette: " << WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu) << std::endl;
envGetTmp(GaugeField, Umu_smear);
{
auto &Umu = envGet(GaugeField, sGaugeName);
LOG(Message) << "Initial plaquette: " << WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu) << std::endl;
Umu_smear = Umu;
const StoutParameters &Stout{par().Stout};
envGetTmp(GaugeField, Umu_stout);
Smear_Stout<PeriodicGimplR> LS(Stout.parm, Tdir); // spatial smearing only
Smear_Stout<PeriodicGimplR> LS(Stout.rho, Tdir); // spatial smearing only
for (int i = 0; i < Stout.steps; i++) {
LS.smear(Umu_stout, Umu_smear);
Umu_smear = Umu_stout;
@ -246,8 +223,8 @@ void TLapEvec<GImpl>::execute(void)
// Invert Peardon Nabla operator separately on each time-slice
////////////////////////////////////////////////////////////////////////
auto & eig4d = envGet(DistilEP, getName() );
envGetTmp(std::vector<DistilEP>, eig); // Eigenpack for each timeslice
auto & eig4d = envGet(LapEvecs, getName() );
envGetTmp(std::vector<LapEvecs>, eig); // Eigenpack for each timeslice
envGetTmp(LatticeGaugeField, UmuNoTime); // Gauge field without time dimension
envGetTmp(LatticeColourVector, src);
const int Ntlocal{gridHD->LocalDimensions()[Tdir]};
@ -267,12 +244,6 @@ void TLapEvec<GImpl>::execute(void)
<< " with parameters (alpha,beta) = (" << ChebPar.alpha << "," << ChebPar.beta << ")" << std::endl;
Chebyshev<LatticeColourVector> Cheb(ChebPar.alpha,ChebPar.beta,ChebPar.PolyOrder);
//from Test_Cheby.cc
//if( Ntfirst == 0 && t==0) {
//std::ofstream of("cheby_" + std::to_string(ChebPar.alpha) + "_" + std::to_string(ChebPar.beta) + "_" + std::to_string(ChebPar.PolyOrder));
//Cheb.csv(of);
//}
// Construct source vector according to Test_dwf_compressed_lanczos.cc
src = 11.0;
RealD nn = norm2(src);

43
Hadrons/Modules/MDistil/Noises.hpp
View File

@ -52,9 +52,11 @@ class NoisesPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(NoisesPar,
std::string, UniqueIdentifier,
int, nnoise,
int, nvec,
DistilParameters, Distil);
std::string, UniqueIdentifier,
std::string, TI,
std::string, LI)
};
template <typename FImpl>
@ -103,36 +105,41 @@ std::vector<std::string> TNoises<FImpl>::getOutput(void)
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TNoises<FImpl>::setup(void)
{
const int Nt{env().getGrid()->GlobalDimensions()[Tdir]};
const int Nt{env().getDim(Tdir)};
const int Ns{Grid::QCD::Ns};
const int nnoise{par().nnoise};
const int nvec{par().nvec};
const DistilParameters & Distil{par().Distil};
//envCreate(std::vector<Complex>, getName(), 1, nvec*Distil.Ns*Distil.Nt*Distil.nnoise);
envCreate(NoiseTensor, getName(), 1, Distil.nnoise, Nt, nvec, Distil.Ns);
const int TI{ Hadrons::MDistil::DistilParameters::ParameterDefault( par().TI, Nt, true) };
const int LI{ Hadrons::MDistil::DistilParameters::ParameterDefault( par().LI, nvec, true) };
envCreate(NoiseTensor, getName(), 1, nnoise, Nt, nvec, Ns);
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TNoises<FImpl>::execute(void)
{
const std::string &UniqueIdentifier{par().UniqueIdentifier};
auto &noise = envGet(NoiseTensor, getName());
const int Nt{env().getDim(Tdir)};
const int Ns{Grid::QCD::Ns};
const int nnoise{par().nnoise};
const int nvec{par().nvec};
const DistilParameters & Distil{par().Distil};
const int nnoise{Distil.nnoise};
const int Nt{env().getGrid()->GlobalDimensions()[Tdir]};
const int Ns{Distil.Ns};
const int TI{Distil.TI};
const int LI{Distil.LI};
const bool full_tdil{TI==Nt};
const bool exact_distillation{full_tdil && LI==nvec};
const int TI{ Hadrons::MDistil::DistilParameters::ParameterDefault( par().TI, Nt, false) };
const int LI{ Hadrons::MDistil::DistilParameters::ParameterDefault( par().LI, nvec, false) };
const bool full_tdil{ TI == Nt }; \
const bool exact_distillation{ full_tdil && LI == nvec }; \
std::string UniqueIdentifier{par().UniqueIdentifier};
if( UniqueIdentifier.length() == 0 ) {
UniqueIdentifier = getName();
}
UniqueIdentifier.append( std::to_string( vm().getTrajectory() ) ); //maybe add more??
GridSerialRNG sRNG;
sRNG.SeedUniqueString(UniqueIdentifier + std::to_string(vm().getTrajectory())); //maybe add more??
sRNG.SeedUniqueString(UniqueIdentifier);
Real rn;
auto &noise = envGet(NoiseTensor, getName());
for( int inoise = 0; inoise < nnoise; inoise++ ) {
for( int t = 0; t < Nt; t++ ) {
for( int ivec = 0; ivec < nvec; ivec++ ) {

94
Hadrons/Modules/MDistil/PerambFromSolve.hpp
View File

@ -47,21 +47,6 @@ BEGIN_HADRONS_NAMESPACE
* PerambFromSolve *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MDistil)
/*
struct DistilParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(DistilParameters,
int, TI,
int, LI,
int, nnoise,
int, tsrc,
int, SI,
int, Ns,
int, Nt,
int, Nt_inv)
DistilParameters() = default;
template <class ReaderClass> DistilParameters(Reader<ReaderClass>& Reader){read(Reader,"Distil",*this);}
};
*/
class PerambFromSolvePar: Serializable
{
@ -122,91 +107,54 @@ TPerambFromSolve<FImpl>::~TPerambFromSolve(void)
template <typename FImpl>
std::vector<std::string> TPerambFromSolve<FImpl>::getInput(void)
{
std::vector<std::string> in;
in.push_back(par().solve);
in.push_back(par().eigenPack);
return in;
return std::vector<std::string>{ par().solve, par().eigenPack };
}
template <typename FImpl>
std::vector<std::string> TPerambFromSolve<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
return std::vector<std::string>{ getName() };
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TPerambFromSolve<FImpl>::setup(void)
{
Cleanup();
const int nvec{par().nvec};
const DistilParameters & Distil{par().Distil};
const int LI{Distil.LI};
const int nnoise{Distil.nnoise};
const int Nt_inv{Distil.Nt_inv};
const int Ns{Distil.Ns};
std::array<std::string,6> sIndexNames{"Nt", "nvec", "LI", "nnoise", "Nt_inv", "SI"};
grid4d = env().getGrid();
grid3d = MakeLowerDimGrid(grid4d);
const int Nt{grid4d->GlobalDimensions()[Tdir]};
const int nvec_reduced{par().nvec_reduced};
const int LI_reduced{par().LI_reduced};
//envCreate(Perambulator<SpinVector COMMA 6 COMMA sizeof(Real)>, getName(), 1,
// sIndexNames,Distil.Nt,nvec,Distil.LI,Distil.nnoise,Distil.Nt_inv,Distil.SI);
envCreate(Perambulator<SpinVector COMMA 6 COMMA sizeof(Real)>, getName(), 1,
sIndexNames,Nt,nvec_reduced,LI_reduced,Distil.nnoise,Distil.Nt_inv,Distil.SI);
envCreate(std::vector<Complex>, getName() + "_noise", 1,
nvec*Distil.Ns*Nt*Distil.nnoise);
envTmp(LatticeColourVector, "result_3d",1,LatticeColourVector(grid3d));
envTmp(LatticeColourVector, "evec3d",1,LatticeColourVector(grid3d));
envTmpLat(LatticeColourVector, "result_nospin");
Cleanup();
DISTIL_PARAMETERS_DEFINE( true );
grid4d = env().getGrid();
grid3d = MakeLowerDimGrid(grid4d);
const int nvec_reduced{par().nvec_reduced};
const int LI_reduced{par().LI_reduced};
//std::array<std::string,6> sIndexNames{"Nt", "nvec", "LI", "nnoise", "Nt_inv", "SI"};
envCreate(Perambulator, getName(), 1, PerambIndexNames,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));
envTmpLat(LatticeColourVector, "result_nospin");
}
template <typename FImpl>
void TPerambFromSolve<FImpl>::Cleanup(void)
{
if( grid3d != nullptr ) {
delete grid3d;
grid3d = nullptr;
}
grid4d = nullptr;
if( grid3d != nullptr ) {
delete grid3d;
grid3d = nullptr;
}
grid4d = nullptr;
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TPerambFromSolve<FImpl>::execute(void)
{
GridCartesian * grid4d = env().getGrid();
const int Nt{grid4d->GlobalDimensions()[Tdir]};
const int Ntlocal{grid4d->LocalDimensions()[3]};
const int Ntfirst{grid4d->LocalStarts()[3]};
const int nvec_reduced{par().nvec_reduced};
const int LI_reduced{par().LI_reduced};
const int nvec{par().nvec};
const DistilParameters & Distil{par().Distil};
const int LI{Distil.LI};
const int TI{Distil.TI};
const int SI{Distil.SI};
const int nnoise{Distil.nnoise};
const int Nt_inv{Distil.Nt_inv};
const int tsrc{Distil.tsrc};
const int Ns{Distil.Ns};
const bool full_tdil{TI==Nt};
const bool exact_distillation{full_tdil && LI==nvec};
auto &perambulator = envGet(Perambulator<SpinVector COMMA 6 COMMA sizeof(Real)>, getName());
DISTIL_PARAMETERS_DEFINE( false );
auto &perambulator = envGet(Perambulator, getName());
auto &solve = envGet(std::vector<FermionField>, par().solve);
auto &epack = envGet(Grid::Hadrons::EigenPack<LatticeColourVector>, par().eigenPack);

105
Hadrons/Modules/MDistil/Perambulator.hpp
View File

@ -30,14 +30,6 @@
#ifndef Hadrons_MDistil_Perambulator_hpp_
#define Hadrons_MDistil_Perambulator_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>
@ -53,41 +45,42 @@ class PerambulatorPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(PerambulatorPar,
std::string, eigenPack,
std::string, lapevec,
std::string, solver,
std::string, noise,
std::string, PerambFileName, //stem!!!
std::string, UniqueIdentifier,
bool, multiFile,
int, nvec,
DistilParameters, Distil,
std::string, solver);
DistilParameters, Distil);
};
template <typename FImpl>
class TPerambulator: public Module<PerambulatorPar>
{
public:
FERM_TYPE_ALIASES(FImpl,);
SOLVER_TYPE_ALIASES(FImpl,);
// constructor
TPerambulator(const std::string name);
// destructor
virtual ~TPerambulator(void);
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
FERM_TYPE_ALIASES(FImpl,);
SOLVER_TYPE_ALIASES(FImpl,);
// constructor
TPerambulator(const std::string name);
// destructor
virtual ~TPerambulator(void);
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
protected:
// These variables are created in setup() and freed in Cleanup()
GridCartesian * grid3d; // Owned by me, so I must delete it
GridCartesian * grid4d; // Owned by environment (so I won't delete it)
virtual void Cleanup(void);
protected:
virtual void Cleanup(void);
private:
unsigned int Ls_;
// These variables are created in setup() and freed in Cleanup()
GridCartesian * grid3d; // Owned by me, so I must delete it
GridCartesian * grid4d; // Owned by environment (so I won't delete it)
// Other members
unsigned int Ls_;
std::string sLapEvecName;
std::string sNoiseName;
};
// Can't name the module Perambulator, because that's what we've called the object
@ -113,21 +106,17 @@ TPerambulator<FImpl>::~TPerambulator(void)
template <typename FImpl>
std::vector<std::string> TPerambulator<FImpl>::getInput(void)
{
std::vector<std::string> in;
in.push_back(par().eigenPack);
in.push_back(par().solver);
in.push_back(par().noise);
return in;
sLapEvecName = par().lapevec;
sNoiseName = par().noise;
if( sNoiseName.length() == 0 )
sNoiseName = getName() + "_noise";
return {sLapEvecName, par().solver, sNoiseName };
}
template <typename FImpl>
std::vector<std::string> TPerambulator<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName(),getName() + "_unsmeared_sink"};
return out;
return {getName(), getName() + "_unsmeared_sink"};
}
// setup ///////////////////////////////////////////////////////////////////////
@ -137,17 +126,9 @@ void TPerambulator<FImpl>::setup(void)
Cleanup();
grid4d = env().getGrid();
grid3d = MakeLowerDimGrid(grid4d);
const int Nt{grid4d->GlobalDimensions()[Tdir]};
const int nvec{par().nvec};
const DistilParameters & Distil{par().Distil};
const int LI{Distil.LI};
const int nnoise{Distil.nnoise};
const int Nt_inv{Distil.Nt_inv}; // TODO: PROBABLY BETTER: if (full_tdil) Nt_inv=1; else Nt_inv = TI;
const int Ns{Distil.Ns};
std::array<std::string,6> sIndexNames{"Nt", "nvec", "LI", "nnoise", "Nt_inv", "SI"};
DISTIL_PARAMETERS_DEFINE( true );
envCreate(Perambulator<SpinVector COMMA 6 COMMA sizeof(Real)>, getName(), 1,
sIndexNames,Nt,nvec,Distil.LI,Distil.nnoise,Distil.Nt_inv,Distil.SI);
envCreate(Perambulator, getName(), 1, PerambIndexNames,Nt,nvec,LI,nnoise,Nt_inv,SI);
envCreate(std::vector<FermionField>, getName() + "_unsmeared_sink", 1,
nnoise*LI*Ns*Nt_inv, envGetGrid(FermionField));
@ -181,18 +162,7 @@ void TPerambulator<FImpl>::Cleanup(void)
template <typename FImpl>
void TPerambulator<FImpl>::execute(void)
{
const int Nt{grid4d->GlobalDimensions()[Tdir]};
const int nvec{par().nvec};
const DistilParameters & Distil{par().Distil};
const int LI{Distil.LI};
const int SI{Distil.SI};
const int TI{Distil.TI};
const int nnoise{Distil.nnoise};
const int tsrc{Distil.tsrc};
const int Ns{Distil.Ns};
const bool full_tdil{TI==Nt};
const bool exact_distillation{full_tdil && LI==nvec};
const int Nt_inv{full_tdil ? 1 : TI}; // TODO: PROBABLY BETTER: if (full_tdil) Nt_inv=1; else Nt_inv = TI;
DISTIL_PARAMETERS_DEFINE( false );
auto &solver=envGet(Solver, par().solver);
auto &mat = solver.getFMat();
@ -200,12 +170,9 @@ void TPerambulator<FImpl>::execute(void)
envGetTmp(FermionField, v5dtmp);
envGetTmp(FermionField, v5dtmp_sol);
const std::string &UniqueIdentifier{par().UniqueIdentifier};
//auto &noise = envGet(std::vector<Complex>, par().noise);
auto &noise = envGet(NoiseTensor, par().noise);
auto &perambulator = envGet(Perambulator<SpinVector COMMA 6 COMMA sizeof(Real)>, getName());
auto &epack = envGet(Grid::Hadrons::EigenPack<LatticeColourVector>, par().eigenPack);
auto &noise = envGet(NoiseTensor, sNoiseName);
auto &perambulator = envGet(Perambulator, getName());
auto &epack = envGet(LapEvecs, sLapEvecName);
auto &unsmeared_sink = envGet(std::vector<FermionField>, getName() + "_unsmeared_sink");

12
Hadrons/Modules/MIO/LoadPerambulator.hpp
View File

@ -101,20 +101,16 @@ std::vector<std::string> TLoadPerambulator<FImpl>::getOutput(void)
template <typename FImpl>
void TLoadPerambulator<FImpl>::setup(void)
{
GridCartesian * grid4d = env().getGrid();
const int Nt{grid4d->GlobalDimensions()[Tdir]};
const int nvec{par().nvec};
const MDistil::DistilParameters & Distil{par().Distil};
std::array<std::string,6> sIndexNames{"Nt", "nvec", "LI", "nnoise", "Nt_inv", "SI"};
envCreate(MDistil::Perambulator<SpinVector COMMA 6 COMMA sizeof(Real)>, getName(), 1,
sIndexNames,Nt,nvec,Distil.LI,Distil.nnoise,Distil.Nt_inv,Distil.SI);
DISTIL_PARAMETERS_DEFINE( true );
//std::array<std::string,6> sIndexNames{"Nt", "nvec", "LI", "nnoise", "Nt_inv", "SI"};
envCreate(MDistil::Perambulator, getName(), 1, MDistil::PerambIndexNames,Nt,nvec,LI,nnoise,Nt_inv,SI);
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TLoadPerambulator<FImpl>::execute(void)
{
auto &perambulator = envGet(MDistil::Perambulator<SpinVector COMMA 6 COMMA sizeof(Real)>, getName());
auto &perambulator = envGet(MDistil::Perambulator, getName());
const std::string sPerambName{par().PerambFileName + "." + std::to_string(vm().getTrajectory())};
perambulator.read(sPerambName.c_str());
}

4
Hadrons/modules.inc
View File

@ -43,9 +43,7 @@ modules_cc =\
Modules/MDistil/BC2.cc \
Modules/MDistil/Noises.cc \
Modules/MDistil/BContraction.cc \
Modules/MDistil/DistilSource.cc \
Modules/MDistil/LapEvec.cc \
Modules/MDistil/DistilSink.cc \
Modules/MDistil/Baryon2pt.cc \
Modules/MDistil/Perambulator.cc \
Modules/MDistil/PerambFromSolve.cc \
@ -118,10 +116,8 @@ modules_hpp =\
Modules/MScalar/FreeProp.hpp \
Modules/MScalar/Scalar.hpp \
Modules/MScalar/ChargedProp.hpp \
Modules/MDistil/DistilSource.hpp \
Modules/MDistil/LapEvec.hpp \
Modules/MDistil/DistilVectors.hpp \
Modules/MDistil/DistilSink.hpp \
Modules/MDistil/Noises.hpp \
Modules/MDistil/BC2.hpp \
Modules/MDistil/Perambulator.hpp \

294
tests/hadrons/Test_distil.cc
View File

@ -70,7 +70,7 @@ void test_SolverS(Application &application)
{
std::string boundary = "1 1 1 -1";
MAction::DWF::Par actionPar;
actionPar.gauge = "gauge";
actionPar.gauge = "LapEvec_gauge";
actionPar.Ls = 16;
actionPar.M5 = 1.8;
actionPar.mass = 0.005;
@ -90,19 +90,16 @@ void test_SolverS(Application &application)
void test_LapEvec(Application &application)
{
const char szGaugeName[] = "gauge";
const char szModuleName[] = "LapEvec";
// gauge field
application.createModule<MGauge::Random>(szGaugeName);
std::string sGaugeName{szModuleName};
sGaugeName.append( "_gauge" );
application.createModule<MGauge::Random>(sGaugeName);
// Now make an instance of the LapEvec object
MDistil::LapEvecPar p;
p.gauge = szGaugeName;
//p.EigenPackName = "ePack";
//p.Distil.TI = 8;
//p.Distil.LI = 3;
//p.Distil.Nnoise = 2;
//p.Distil.tSrc = 0;
//p.gauge = sGaugeName; // This is the default for LapEvec, so no need to be explicit
p.Stout.steps = 3;
p.Stout.parm = 0.2;
p.Stout.rho = 0.2;
p.Cheby.PolyOrder = 11;
p.Cheby.alpha = 0.3;
p.Cheby.beta = 12.5;
@ -112,58 +109,90 @@ void test_LapEvec(Application &application)
p.Lanczos.MaxIt = 1000;
p.Lanczos.resid = 1e-2;
p.Lanczos.IRLLog = 0;
application.createModule<MDistil::LapEvec>("LapEvec",p);
application.createModule<MDistil::LapEvec>(szModuleName,p);
}
/////////////////////////////////////////////////////////////
// Noises
/////////////////////////////////////////////////////////////
std::string test_Noises(Application &application, const std::string &sNoiseBaseName ) {
// DistilVectors parameters
MDistil::NoisesPar NoisePar;
NoisePar.nnoise = 1;
NoisePar.nvec = 5;
std::string sNoiseName{sNoiseBaseName + "_noise"};
application.createModule<MDistil::Noises>(sNoiseName,NoisePar);
return sNoiseName;
}
/////////////////////////////////////////////////////////////
// Perambulators
/////////////////////////////////////////////////////////////
void test_Perambulators(Application &application)
void test_Perambulators( Application &application, const char * pszSuffix = nullptr )
{
std::string sModuleName{ "Peramb" };
if( pszSuffix )
sModuleName.append( pszSuffix );
test_Noises(application, sModuleName);
// Perambulator parameters
MDistil::Peramb::Par PerambPar;
PerambPar.eigenPack="LapEvec";
PerambPar.noise="Peramb_noise";
PerambPar.PerambFileName="peramb.bin";
PerambPar.UniqueIdentifier="full_dilution";
PerambPar.lapevec = "LapEvec";
PerambPar.PerambFileName = sModuleName + ".bin";
PerambPar.solver="CG_s";
PerambPar.Distil.tsrc = 0;
PerambPar.Distil.nnoise = 1;
PerambPar.Distil.LI=5;
PerambPar.Distil.SI=4;
PerambPar.Distil.TI=8;
PerambPar.nvec=5;
PerambPar.Distil.Ns=4;
PerambPar.Distil.Nt_inv=1;
//PerambPar.Solver.mass=0.005;
//PerambPar.Solver.M5=1.8;
//PerambPar.Ls=16;
//PerambPar.Solver.CGPrecision=1e-8;
//PerambPar.Solver.MaxIterations=10000;
application.createModule<MDistil::Peramb>("Peramb",PerambPar);
application.createModule<MDistil::Peramb>( sModuleName, PerambPar );
}
/////////////////////////////////////////////////////////////
// DistilVectors
/////////////////////////////////////////////////////////////
#define TEST_DISTIL_VECTORS_COMMON \
std::string sModuleName{"DistilVecs"}; \
if( pszSuffix ) \
sModuleName.append( pszSuffix ); \
std::string sPerambName{"Peramb"}; \
if( pszSuffix ) \
sPerambName.append( pszSuffix ); \
MDistil::DistilVectors::Par DistilVecPar; \
DistilVecPar.noise = sPerambName + "_noise"; \
DistilVecPar.perambulator = sPerambName; \
DistilVecPar.lapevec = "LapEvec"; \
DistilVecPar.tsrc = 0; \
if( pszNvec ) \
DistilVecPar.nvec = pszNvec
#define TEST_DISTIL_VECTORS_COMMON_END \
application.createModule<MDistil::DistilVectors>(sModuleName,DistilVecPar)
void test_DistilVectors(Application &application, const char * pszSuffix = nullptr, const char * pszNvec = nullptr )
{
TEST_DISTIL_VECTORS_COMMON;
TEST_DISTIL_VECTORS_COMMON_END;
}
void test_DistilVectorsSS(Application &application, const char * pszSink, const char * pszSource,
const char * pszSuffix = nullptr, const char * pszNvec = nullptr )
{
TEST_DISTIL_VECTORS_COMMON;
if( pszSink )
DistilVecPar.sink = pszSink;
if( pszSource )
DistilVecPar.source = pszSource;
TEST_DISTIL_VECTORS_COMMON_END;
}
/////////////////////////////////////////////////////////////
// Multiple Perambulators
/////////////////////////////////////////////////////////////
void test_MultiPerambulators(Application &application)
{
// Perambulator parameters
MDistil::Peramb::Par PerambPar;
PerambPar.eigenPack="LapEvec";
PerambPar.UniqueIdentifier="full_dilution";
PerambPar.PerambFileName="Peramb5";
PerambPar.solver="CG_s";
PerambPar.Distil.tsrc = 0;
PerambPar.Distil.nnoise = 1;
PerambPar.Distil.LI=5;
PerambPar.Distil.SI=4;
PerambPar.Distil.TI=8;
PerambPar.nvec=5;
PerambPar.Distil.Ns=4;
PerambPar.Distil.Nt_inv=1;
application.createModule<MDistil::Peramb>("Peramb5",PerambPar);
test_Perambulators( application, "5" );
MDistil::PerambFromSolve::Par SolvePar;
SolvePar.eigenPack="LapEvec";
SolvePar.PerambFileName="Peramb2";
@ -175,35 +204,16 @@ void test_MultiPerambulators(Application &application)
SolvePar.nvec=5;
SolvePar.nvec_reduced=2;
SolvePar.LI_reduced=2;
SolvePar.Distil.Ns=4;
SolvePar.Distil.Nt_inv=1;
application.createModule<MDistil::PerambFromSolve>("Peramb2",SolvePar);
SolvePar.PerambFileName="Peramb3";
SolvePar.nvec_reduced=3;
SolvePar.LI_reduced=3;
application.createModule<MDistil::PerambFromSolve>("Peramb3",SolvePar);
MDistil::DistilVectors::Par DistilVecPar;
DistilVecPar.noise="Peramb5_noise";
DistilVecPar.perambulator="Peramb2";
DistilVecPar.lapevec ="LapEvec";
DistilVecPar.tsrc = 0;
//DistilVecPar.nnoise = 1;
DistilVecPar.LI=2;
DistilVecPar.SI=4;
DistilVecPar.TI=8;
DistilVecPar.nvec=2;
DistilVecPar.Ns=4;
DistilVecPar.Nt_inv=1;
application.createModule<MDistil::DistilVectors>("DistilVecs2",DistilVecPar);
DistilVecPar.perambulator="Peramb3";
DistilVecPar.LI=3;
DistilVecPar.nvec=3;
application.createModule<MDistil::DistilVectors>("DistilVecs3",DistilVecPar);
//DistilVecPar.perambulator="Peramb5_perambulator_light";
DistilVecPar.perambulator="Peramb5";
DistilVecPar.LI=5;
DistilVecPar.nvec=5;
application.createModule<MDistil::DistilVectors>("DistilVecs5",DistilVecPar);
test_DistilVectors( application, "2", "2" );
test_DistilVectors( application, "3", "3" );
test_DistilVectors( application, "5", "5" );
MContraction::A2AMesonField::Par A2AMesonFieldPar;
A2AMesonFieldPar.left="DistilVecs2_rho";
A2AMesonFieldPar.right="DistilVecs2_rho";
@ -235,88 +245,6 @@ void test_MultiPerambulators(Application &application)
application.createModule<MContraction::A2AMesonField>("DistilMesonFieldPhi5",A2AMesonFieldPar);
}
void test_Noises(Application &application) {
// DistilVectors parameters
MDistil::NoisesPar NoisePar;
NoisePar.UniqueIdentifier = "full_dilution";
NoisePar.nvec = 5;
NoisePar.Distil.TI = 8;
NoisePar.Distil.SI = 4;
NoisePar.Distil.LI = 5;
NoisePar.Distil.nnoise = 1;
NoisePar.Distil.Ns = 4;
NoisePar.Distil.Nt_inv = 1;
NoisePar.Distil.tsrc = 0;
application.createModule<MDistil::Noises>("Peramb_noise",NoisePar);
}
/////////////////////////////////////////////////////////////
// DistilVectors
/////////////////////////////////////////////////////////////
void test_DistilVectors(Application &application)
{
test_Noises(application);
// DistilVectors parameters
MDistil::DistilVectors::Par DistilVecPar;
DistilVecPar.noise="Peramb_noise";
//DistilVecPar.perambulator="Peramb_perambulator_light";
DistilVecPar.perambulator="Peramb";
DistilVecPar.lapevec="LapEvec";
DistilVecPar.tsrc = 0;
DistilVecPar.LI=5;
DistilVecPar.SI=4;
DistilVecPar.TI=8;
DistilVecPar.nvec=5;
DistilVecPar.Ns=4;
DistilVecPar.Nt_inv=1;
application.createModule<MDistil::DistilVectors>("DistilVecs",DistilVecPar);
}
void test_PerambulatorsS(Application &application)
{
// Perambulator parameters
MDistil::Peramb::Par PerambPar;
PerambPar.eigenPack="LapEvec";
PerambPar.PerambFileName="perambS.bin";
PerambPar.UniqueIdentifier="full_dilution";
PerambPar.solver="CG_s";
PerambPar.Distil.tsrc = 0;
PerambPar.Distil.nnoise = 1;
PerambPar.Distil.LI=5;
PerambPar.Distil.SI=4;
PerambPar.Distil.TI=8;
PerambPar.nvec=5;
PerambPar.Distil.Ns=4;
PerambPar.Distil.Nt_inv=1;
//PerambPar.Solver.mass=0.005; //strange mass???
//PerambPar.Solver.M5=1.8;
//PerambPar.Ls=16;
//PerambPar.Solver.CGPrecision=1e-8;
//PerambPar.Solver.MaxIterations=10000;
application.createModule<MDistil::Peramb>("PerambS",PerambPar);
}
/////////////////////////////////////////////////////////////
// DistilVectors
/////////////////////////////////////////////////////////////
void test_DistilVectorsS(Application &application)
{
// DistilVectors parameters
MDistil::DistilVectors::Par DistilVecPar;
DistilVecPar.noise="PerambS_noise";
//DistilVecPar.perambulator="PerambS_perambulator_light";
DistilVecPar.perambulator="PerambS";
DistilVecPar.lapevec="LapEvec";
DistilVecPar.tsrc = 0;
//DistilVecPar.nnoise = 1;
DistilVecPar.LI=5;
DistilVecPar.SI=4;
DistilVecPar.TI=32;
DistilVecPar.nvec=5;
DistilVecPar.Ns=4;
DistilVecPar.Nt_inv=1;
application.createModule<MDistil::DistilVectors>("DistilVecsS",DistilVecPar);
}
/////////////////////////////////////////////////////////////
// MesonSink
/////////////////////////////////////////////////////////////
@ -516,35 +444,10 @@ void test_PerambulatorsSolve(Application &application)
PerambFromSolvePar.PerambFileName="perambAslashS.bin";
PerambFromSolvePar.Distil.tsrc = 0;
PerambFromSolvePar.Distil.nnoise = 1;
PerambFromSolvePar.Distil.LI=5;
PerambFromSolvePar.Distil.SI=4;
PerambFromSolvePar.Distil.TI=8;
PerambFromSolvePar.nvec=5;
PerambFromSolvePar.Distil.Ns=4;
PerambFromSolvePar.Distil.Nt_inv=1;
application.createModule<MDistil::PerambFromSolve>("PerambAslashS",PerambFromSolvePar);
}
/////////////////////////////////////////////////////////////
// DistilVectors
/////////////////////////////////////////////////////////////
void test_DistilVectorsAslashSeq(Application &application)
{
// DistilVectors parameters
MDistil::DistilSink::Par DistilSinkPar;
DistilSinkPar.perambulator="PerambAslashS";
DistilSinkPar.eigenPack="LapEvec";
DistilSinkPar.tsrc = 0;
DistilSinkPar.nnoise = 1;
DistilSinkPar.LI=5;
DistilSinkPar.SI=4;
DistilSinkPar.TI=8;
DistilSinkPar.nvec=5;
DistilSinkPar.Ns=4;
DistilSinkPar.Nt=8;
DistilSinkPar.Nt_inv=1;
application.createModule<MDistil::DistilSink>("DistilVecsAslashSeq",DistilSinkPar);
}
bool bNumber( int &ri, const char * & pstr, bool bGobbleWhiteSpace = true )
{
if( bGobbleWhiteSpace )
@ -702,7 +605,7 @@ bool DebugEigenTest()
// Test initialisation of an array of strings
for( auto a : as )
std::cout << a << std::endl;
Grid::Hadrons::MDistil::Perambulator<Complex,3,sizeof(Real)> p{as,2,7,2};
Grid::Hadrons::MDistil::NamedTensor<Complex,3,sizeof(Real)> p{as,2,7,2};
DebugShowTensor(p, "p");
std::cout << "p.IndexNames follow" << std::endl;
for( auto a : p.IndexNames )
@ -950,7 +853,7 @@ bool DebugGridTensorTest( void )
bool ConvertPeramb(const char * pszSource, const char * pszDest) {
std::array<std::string,6> sIndexNames{"Nt", "nvec", "LI", "nnoise", "Nt_inv", "SI"};
Grid::Hadrons::MDistil::Perambulator<SpinVector, 6, sizeof(Real)> p(sIndexNames);
Grid::Hadrons::MDistil::Perambulator p(sIndexNames);
p.ReadBinary( pszSource );
p.write(pszDest);
return true;
@ -1016,21 +919,21 @@ int main(int argc, char *argv[])
break;
case 2:
test_Global( application );
test_SolverS( application );
test_LapEvec( application );
test_SolverS( application );
test_Perambulators( application );
break;
default: // 3
test_Global( application );
test_SolverS( application );
test_LapEvec( application );
test_SolverS( application );
test_Perambulators( application );
test_DistilVectors( application );
break;
case 4:
test_Global( application );
test_SolverS( application );
test_LapEvec( application );
test_SolverS( application );
test_Perambulators( application );
test_DistilVectors( application );
test_MesonField( application, "Phi", "_phi" );
@ -1038,28 +941,28 @@ int main(int argc, char *argv[])
break;
case 5: // 3
test_Global( application );
test_SolverS( application );
test_LapEvec( application );
test_SolverS( application );
test_Perambulators( application );
test_DistilVectors( application );
test_PerambulatorsS( application );
test_DistilVectorsS( application );
test_Perambulators( application, "S" );
test_DistilVectors( application, "S" );
test_MesonField( application, "SPhi", "S_phi" );
test_MesonField( application, "SRho", "S_rho" );
break;
#ifdef DISTIL_PRE_RELEASE
case 6: // 3
test_Global( application );
test_SolverS( application );
test_LapEvec( application );
test_SolverS( application );
test_Perambulators( application );
test_g5_sinks( application );
test_MesonSink( application );
break;
case 7: // 3
test_Global( application );
test_SolverS( application );
test_LapEvec( application );
test_SolverS( application );
test_Perambulators( application );
test_DistilVectors( application );
test_BaryonFieldPhi( application );
@ -1068,8 +971,8 @@ int main(int argc, char *argv[])
#endif
case 8: // 3
test_Global( application );
test_SolverS( application );
test_LapEvec( application );
test_SolverS( application );
test_Perambulators( application );
test_DistilVectors( application );
test_MesonField( application, "Phi", "_phi" );
@ -1083,8 +986,8 @@ int main(int argc, char *argv[])
break;
case 10: // 3
test_Global( application );
test_SolverS( application );
test_LapEvec( application );
test_SolverS( application );
test_Perambulators( application );
test_g5_sinks( application );
test_em( application );
@ -1092,8 +995,8 @@ int main(int argc, char *argv[])
break;
case 11: // 3
test_Global( application );
test_SolverS( application );
test_LapEvec( application );
test_SolverS( application );
test_Perambulators( application );
test_DistilVectors( application );
test_BaryonFieldPhi2( application );
@ -1102,26 +1005,31 @@ int main(int argc, char *argv[])
#endif
case 12: // 3
test_Global( application );
test_SolverS( application );
test_LapEvec( application );
test_PerambulatorsS( application );
test_SolverS( application );
test_Perambulators( application, "S" );
test_em( application );
test_AslashSeq( application );
test_PerambulatorsSolve( application );
test_DistilVectorsAslashSeq( application );
test_DistilVectorsSS( application, "AslashSeq", nullptr, "S" );
test_MesonField( application, "AslashSeq" );
break;
case 13:
test_Global( application );
test_SolverS( application );
test_LapEvec( application );
test_SolverS( application );
test_MultiPerambulators( application );
break;
}
// execution
application.saveParameterFile("test_distil.xml");
LOG(Warning) << "These parameters are designed to run on a laptop usid --grid 4.4.4.8" << std::endl;
application.run();
static const char XmlFileName[] = "test_distil.xml";
application.saveParameterFile( XmlFileName );
const std::vector<int> &lat{GridDefaultLatt()};
if( lat.size() == 4 && lat[0] == 4 && lat[1] == 4 && lat[2] == 4 && lat[3] == 8 )
application.run();
else
LOG(Warning) << "The parameters in " << XmlFileName << " are designed to run on a laptop usid --grid 4.4.4.8" << std::endl;
// epilogue
LOG(Message) << "Grid is finalizing now" << std::endl;