From 334f29becbb7ac67af3af91c19e584ade1c6a2b3 Mon Sep 17 00:00:00 2001 From: Michael Marshall <43034299+mmphys@users.noreply.github.com> Date: Tue, 30 Apr 2019 23:53:57 +0100 Subject: [PATCH] Fairly close to ready for release. Felix and I to review, then submit for release --- Hadrons/Distil.hpp | 62 +++-- Hadrons/Modules.hpp | 2 - Hadrons/Modules/MDistil/DistilSink.cc | 36 --- Hadrons/Modules/MDistil/DistilSink.hpp | 236 ---------------- Hadrons/Modules/MDistil/DistilSource.cc | 36 --- Hadrons/Modules/MDistil/DistilSource.hpp | 250 ----------------- Hadrons/Modules/MDistil/DistilVectors.hpp | 75 +++-- Hadrons/Modules/MDistil/LapEvec.hpp | 77 ++--- Hadrons/Modules/MDistil/Noises.hpp | 43 +-- Hadrons/Modules/MDistil/PerambFromSolve.hpp | 94 ++----- Hadrons/Modules/MDistil/Perambulator.hpp | 105 +++---- Hadrons/Modules/MIO/LoadPerambulator.hpp | 12 +- Hadrons/modules.inc | 4 - tests/hadrons/Test_distil.cc | 294 +++++++------------- 14 files changed, 290 insertions(+), 1036 deletions(-) delete mode 100644 Hadrons/Modules/MDistil/DistilSink.cc delete mode 100644 Hadrons/Modules/MDistil/DistilSink.hpp delete mode 100644 Hadrons/Modules/MDistil/DistilSource.cc delete mode 100644 Hadrons/Modules/MDistil/DistilSource.hpp diff --git a/Hadrons/Distil.hpp b/Hadrons/Distil.hpp index 3728f0d5..58bfc532 100644 --- a/Hadrons/Distil.hpp +++ b/Hadrons/Distil.hpp @@ -33,16 +33,10 @@ #include #include #include - -/****************************************************************************** - Needed to make sure envCreate() (see Hadrons) work with specialisations - with more than one parameter, eg obj - I imagine this exists already? - ******************************************************************************/ - -#ifndef COMMA -#define COMMA , -#endif +#include +#include +#include +#include /****************************************************************************** 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; +using LapEvecs = Grid::Hadrons::EigenPack; // Noise vector index order: nnoise, nt, nvec, ns using NoiseTensor = Eigen::Tensor; 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 DistilParameters(Reader& Reader){read(Reader,"Distil",*this);} + + // Numeric parameter is allowed to be empty (in which case it = Default), + // but assert during setup() if specified but not numeric + + static int ParameterDefault( const std::string & s, int Default, bool bCalledFromSetup ) + { + int i = Default; + if( s.length() > 0 ) { + std::istringstream ss( s ); + ss >> i; + if( bCalledFromSetup ) + assert( !ss.fail() && "Parameter should either be empty or integer" ); + } + return i; + } + + 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; @@ -299,7 +322,7 @@ public: // Construct a named tensor explicitly specifying size of each dimension template - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NamedTensor(std::array &IndexNames_, Eigen::Index firstDimension, IndexTypes... otherDimensions) + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NamedTensor(const std::array &IndexNames_, Eigen::Index firstDimension, IndexTypes... otherDimensions) : tensor(firstDimension, otherDimensions...), IndexNames{NumIndices} { // The number of dimensions used to construct a tensor must be equal to the rank of the tensor. @@ -343,8 +366,9 @@ template struct is_named_tensor -using Perambulator = NamedTensor; +//template +using Perambulator = NamedTensor; +static const std::array PerambIndexNames{"Nt", "nvec", "LI", "nnoise", "Nt_inv", "SI"}; /****************************************************************************** Save NamedTensor binary format (NB: On-disk format is Big Endian) diff --git a/Hadrons/Modules.hpp b/Hadrons/Modules.hpp index 0ffe6b6d..774eafa4 100644 --- a/Hadrons/Modules.hpp +++ b/Hadrons/Modules.hpp @@ -39,10 +39,8 @@ #include #include #include -#include #include #include -#include #include #include #include diff --git a/Hadrons/Modules/MDistil/DistilSink.cc b/Hadrons/Modules/MDistil/DistilSink.cc deleted file mode 100644 index c438dd8b..00000000 --- a/Hadrons/Modules/MDistil/DistilSink.cc +++ /dev/null @@ -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 - Author: Michael Marshall - - This program is free software; you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation; either version 2 of the License, or - (at your option) any later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License along - with this program; if not, write to the Free Software Foundation, Inc., - 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - - See the full license in the file "LICENSE" in the top level distribution directory - *************************************************************************************/ -/* END LEGAL */ - -#include - -using namespace Grid; -using namespace Hadrons; -using namespace MDistil; - -template class Grid::Hadrons::MDistil::TDistilSink; diff --git a/Hadrons/Modules/MDistil/DistilSink.hpp b/Hadrons/Modules/MDistil/DistilSink.hpp deleted file mode 100644 index 9435bcc4..00000000 --- a/Hadrons/Modules/MDistil/DistilSink.hpp +++ /dev/null @@ -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 - Author: Michael Marshall - - This program is free software; you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation; either version 2 of the License, or - (at your option) any later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License along - with this program; if not, write to the Free Software Foundation, Inc., - 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - - See the full license in the file "LICENSE" in the top level distribution directory - *************************************************************************************/ -/* END LEGAL */ - -#ifndef Hadrons_MDistil_DistilSink_hpp_ -#define Hadrons_MDistil_DistilSink_hpp_ - -#include -#include -#include -#include -#include -#include -#include - -// These are members of Distillation -#include - -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 -class TDistilSink: public Module -{ -public: - FERM_TYPE_ALIASES(FImpl,); -public: - // constructor - TDistilSink(const std::string name); - // destructor - virtual ~TDistilSink(void); - // dependency relation - virtual std::vector getInput(void); - virtual std::vector 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, MDistil); - -/****************************************************************************** - * TDistilSink implementation * - ******************************************************************************/ -// constructor ///////////////////////////////////////////////////////////////// -template -TDistilSink::TDistilSink(const std::string name) -: grid3d{nullptr}, grid4d{nullptr}, Module(name) -{} -// destructor -template -TDistilSink::~TDistilSink(void) -{ - Cleanup(); -}; - -// dependencies/products /////////////////////////////////////////////////////// -template -std::vector TDistilSink::getInput(void) -{ - std::vector in; - - in.push_back(par().perambulator); - in.push_back(par().eigenPack); - - return in; -} - -template -std::vector TDistilSink::getOutput(void) -{ - std::vector out = {getName()}; - - return out; -} - -// setup /////////////////////////////////////////////////////////////////////// -template -void TDistilSink::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, getName(), 1, - nnoise*LI*SI*Nt_inv, envGetGrid(FermionField)); - - grid4d = env().getGrid(); - std::vector latt_size = GridDefaultLatt(); - std::vector simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd()); - std::vector mpi_layout = GridDefaultMpi(); - std::vector 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 -void TDistilSink::Cleanup(void) -{ - if( grid3d != nullptr ) { - delete grid3d; - grid3d = nullptr; - } - grid4d = nullptr; -} - -// execution /////////////////////////////////////////////////////////////////// -template -void TDistilSink::execute(void) -{ - - auto &perambulator = envGet(Perambulator, par().perambulator); - auto &epack = envGet(Grid::Hadrons::EigenPack, par().eigenPack); - auto &phi = envGet(std::vector, 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_ diff --git a/Hadrons/Modules/MDistil/DistilSource.cc b/Hadrons/Modules/MDistil/DistilSource.cc deleted file mode 100644 index 3438b8ea..00000000 --- a/Hadrons/Modules/MDistil/DistilSource.cc +++ /dev/null @@ -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 - Author: Michael Marshall - - This program is free software; you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation; either version 2 of the License, or - (at your option) any later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License along - with this program; if not, write to the Free Software Foundation, Inc., - 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - - See the full license in the file "LICENSE" in the top level distribution directory - *************************************************************************************/ -/* END LEGAL */ - -#include - -using namespace Grid; -using namespace Hadrons; -using namespace MDistil; - -template class Grid::Hadrons::MDistil::TDistilSource; diff --git a/Hadrons/Modules/MDistil/DistilSource.hpp b/Hadrons/Modules/MDistil/DistilSource.hpp deleted file mode 100644 index da0c5028..00000000 --- a/Hadrons/Modules/MDistil/DistilSource.hpp +++ /dev/null @@ -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 - Author: Michael Marshall - - This program is free software; you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation; either version 2 of the License, or - (at your option) any later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License along - with this program; if not, write to the Free Software Foundation, Inc., - 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - - See the full license in the file "LICENSE" in the top level distribution directory - *************************************************************************************/ -/* END LEGAL */ - -#ifndef Hadrons_MDistil_DistilSource_hpp_ -#define Hadrons_MDistil_DistilSource_hpp_ - -#include -#include -#include -#include -#include -#include -#include - -// These are members of Distillation -#include - -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 -class TDistilSource: public Module -{ -public: - FERM_TYPE_ALIASES(FImpl,); -public: - // constructor - TDistilSource(const std::string name); - // destructor - virtual ~TDistilSource(void); - // dependency relation - virtual std::vector getInput(void); - virtual std::vector 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, MDistil); - -/****************************************************************************** - * TDistilSource implementation * - ******************************************************************************/ -// constructor ///////////////////////////////////////////////////////////////// -template -TDistilSource::TDistilSource(const std::string name) -: grid3d{nullptr}, grid4d{nullptr}, Module(name) -{} -// destructor -template -TDistilSource::~TDistilSource(void) -{ - Cleanup(); -}; - -// dependencies/products /////////////////////////////////////////////////////// -template -std::vector TDistilSource::getInput(void) -{ - std::vector in; - - in.push_back(par().noise); - in.push_back(par().eigenPack); - - return in; -} - -template -std::vector TDistilSource::getOutput(void) -{ - std::vector out = {getName()}; - - return out; -} - -// setup /////////////////////////////////////////////////////////////////////// -template -void TDistilSource::setup(void) -{ - Cleanup(); - auto &noise = envGet(std::vector, 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, getName(), 1, - nnoise*LI*SI*Nt_inv, envGetGrid(FermionField)); - - grid4d = env().getGrid(); - std::vector latt_size = GridDefaultLatt(); - std::vector simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd()); - std::vector mpi_layout = GridDefaultMpi(); - std::vector 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 -void TDistilSource::Cleanup(void) -{ - if( grid3d != nullptr ) { - delete grid3d; - grid3d = nullptr; - } - grid4d = nullptr; -} - -// execution /////////////////////////////////////////////////////////////////// -template -void TDistilSource::execute(void) -{ - - auto &noise = envGet(std::vector, par().noise); - auto &epack = envGet(Grid::Hadrons::EigenPack, par().eigenPack); - auto &rho = envGet(std::vector, 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_ diff --git a/Hadrons/Modules/MDistil/DistilVectors.hpp b/Hadrons/Modules/MDistil/DistilVectors.hpp index 4a76c68c..a9f6ad14 100644 --- a/Hadrons/Modules/MDistil/DistilVectors.hpp +++ b/Hadrons/Modules/MDistil/DistilVectors.hpp @@ -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 @@ -73,8 +69,6 @@ class TDistilVectors: public Module { public: FERM_TYPE_ALIASES(FImpl,); - // This is the type of perambulator I expect - using DistilPeramb = Perambulator; // 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, MDistil); @@ -171,19 +164,24 @@ void TDistilVectors::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 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( 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( perambulator.tensor.dimension(0) ) && "Perambulator time dimensionality bad" ); + const int TI{ Hadrons::MDistil::DistilParameters::ParameterDefault( par().TI, Nt, true) }; + const int LI{ static_cast( perambulator.tensor.dimension(2) ) }; + const int SI{ static_cast( perambulator.tensor.dimension(5) ) }; + const int Nt_inv{ static_cast( perambulator.tensor.dimension(4) ) }; + const int nnoise{ static_cast( perambulator.tensor.dimension(3) ) }; + assert( nnoise >= static_cast( 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( 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, SourceName, 1, nnoise*LI*SI*Nt_inv, envGetGrid(FermionField)); @@ -225,31 +223,30 @@ template void TDistilVectors::execute(void) { auto &noise = envGet(NoiseTensor, NoiseVectorName); - auto &perambulator = envGet(DistilPeramb, PerambulatorName); + auto &perambulator = envGet(Perambulator, PerambulatorName); auto &epack = envGet(Grid::Hadrons::EigenPack, 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( perambulator.tensor.dimension(2) ) }; + const int SI{ static_cast( perambulator.tensor.dimension(5) ) }; + const int Nt_inv{ static_cast( perambulator.tensor.dimension(4) ) }; + const int nnoise{ static_cast( 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( perambulator.tensor.dimension(1) ), false)}; + const int tsrc{ par().tsrc }; + const bool full_tdil{ TI==Nt }; + int vecindex; int t_inv; if( bMakeSource ) { diff --git a/Hadrons/Modules/MDistil/LapEvec.hpp b/Hadrons/Modules/MDistil/LapEvec.hpp index c2139b0e..1b603ef6 100644 --- a/Hadrons/Modules/MDistil/LapEvec.hpp +++ b/Hadrons/Modules/MDistil/LapEvec.hpp @@ -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 StoutParameters(Reader& 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::~TLapEvec() template std::vector TLapEvec::getInput(void) { - std::vector in = {par().gauge}; - return in; + sGaugeName = par().gauge; + if( sGaugeName.size() == 0 ) { + sGaugeName = getName(); + sGaugeName.append( "_gauge" ); + } + return std::vector{ sGaugeName }; } template @@ -170,19 +164,16 @@ void TLapEvec::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, "eig",1,std::vector(Nt)); + envTmp(std::vector, "eig",1,std::vector(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::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::avgPlaquette(Umu) << std::endl; + envGetTmp(GaugeField, Umu_smear); { + auto &Umu = envGet(GaugeField, sGaugeName); + LOG(Message) << "Initial plaquette: " << WilsonLoops::avgPlaquette(Umu) << std::endl; + Umu_smear = Umu; const StoutParameters &Stout{par().Stout}; envGetTmp(GaugeField, Umu_stout); - Smear_Stout LS(Stout.parm, Tdir); // spatial smearing only + Smear_Stout 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::execute(void) // Invert Peardon Nabla operator separately on each time-slice //////////////////////////////////////////////////////////////////////// - auto & eig4d = envGet(DistilEP, getName() ); - envGetTmp(std::vector, eig); // Eigenpack for each timeslice + auto & eig4d = envGet(LapEvecs, getName() ); + envGetTmp(std::vector, 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::execute(void) << " with parameters (alpha,beta) = (" << ChebPar.alpha << "," << ChebPar.beta << ")" << std::endl; Chebyshev 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); diff --git a/Hadrons/Modules/MDistil/Noises.hpp b/Hadrons/Modules/MDistil/Noises.hpp index fcb102f5..324318ec 100644 --- a/Hadrons/Modules/MDistil/Noises.hpp +++ b/Hadrons/Modules/MDistil/Noises.hpp @@ -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 @@ -103,36 +105,41 @@ std::vector TNoises::getOutput(void) } // setup /////////////////////////////////////////////////////////////////////// + template void TNoises::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, 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 void TNoises::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++ ) { diff --git a/Hadrons/Modules/MDistil/PerambFromSolve.hpp b/Hadrons/Modules/MDistil/PerambFromSolve.hpp index e3160ee5..7c3f2729 100644 --- a/Hadrons/Modules/MDistil/PerambFromSolve.hpp +++ b/Hadrons/Modules/MDistil/PerambFromSolve.hpp @@ -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 DistilParameters(Reader& Reader){read(Reader,"Distil",*this);} -}; -*/ class PerambFromSolvePar: Serializable { @@ -122,91 +107,54 @@ TPerambFromSolve::~TPerambFromSolve(void) template std::vector TPerambFromSolve::getInput(void) { - std::vector in; - - in.push_back(par().solve); - in.push_back(par().eigenPack); - - return in; + return std::vector{ par().solve, par().eigenPack }; } template std::vector TPerambFromSolve::getOutput(void) { - std::vector out = {getName()}; - - return out; + return std::vector{ getName() }; } // setup /////////////////////////////////////////////////////////////////////// template void TPerambFromSolve::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 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, getName(), 1, - // sIndexNames,Distil.Nt,nvec,Distil.LI,Distil.nnoise,Distil.Nt_inv,Distil.SI); - envCreate(Perambulator, getName(), 1, - sIndexNames,Nt,nvec_reduced,LI_reduced,Distil.nnoise,Distil.Nt_inv,Distil.SI); - envCreate(std::vector, 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 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 void TPerambFromSolve::Cleanup(void) { - if( grid3d != nullptr ) { - delete grid3d; - grid3d = nullptr; - } - grid4d = nullptr; + if( grid3d != nullptr ) { + delete grid3d; + grid3d = nullptr; + } + grid4d = nullptr; } - // execution /////////////////////////////////////////////////////////////////// template void TPerambFromSolve::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, getName()); + DISTIL_PARAMETERS_DEFINE( false ); + auto &perambulator = envGet(Perambulator, getName()); auto &solve = envGet(std::vector, par().solve); auto &epack = envGet(Grid::Hadrons::EigenPack, par().eigenPack); diff --git a/Hadrons/Modules/MDistil/Perambulator.hpp b/Hadrons/Modules/MDistil/Perambulator.hpp index 202efae0..325a2386 100644 --- a/Hadrons/Modules/MDistil/Perambulator.hpp +++ b/Hadrons/Modules/MDistil/Perambulator.hpp @@ -30,14 +30,6 @@ #ifndef Hadrons_MDistil_Perambulator_hpp_ #define Hadrons_MDistil_Perambulator_hpp_ -#include -#include -#include -#include -#include -#include -#include - // These are members of Distillation #include @@ -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 class TPerambulator: public Module { public: - FERM_TYPE_ALIASES(FImpl,); - SOLVER_TYPE_ALIASES(FImpl,); - // constructor - TPerambulator(const std::string name); - // destructor - virtual ~TPerambulator(void); - // dependency relation - virtual std::vector getInput(void); - virtual std::vector 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 getInput(void); + virtual std::vector 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::~TPerambulator(void) template std::vector TPerambulator::getInput(void) { - std::vector 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 std::vector TPerambulator::getOutput(void) { - std::vector out = {getName(),getName() + "_unsmeared_sink"}; - - return out; + return {getName(), getName() + "_unsmeared_sink"}; } // setup /////////////////////////////////////////////////////////////////////// @@ -137,17 +126,9 @@ void TPerambulator::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 sIndexNames{"Nt", "nvec", "LI", "nnoise", "Nt_inv", "SI"}; + DISTIL_PARAMETERS_DEFINE( true ); - envCreate(Perambulator, 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, getName() + "_unsmeared_sink", 1, nnoise*LI*Ns*Nt_inv, envGetGrid(FermionField)); @@ -181,18 +162,7 @@ void TPerambulator::Cleanup(void) template void TPerambulator::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::execute(void) envGetTmp(FermionField, v5dtmp); envGetTmp(FermionField, v5dtmp_sol); - const std::string &UniqueIdentifier{par().UniqueIdentifier}; - - //auto &noise = envGet(std::vector, par().noise); - auto &noise = envGet(NoiseTensor, par().noise); - auto &perambulator = envGet(Perambulator, getName()); - auto &epack = envGet(Grid::Hadrons::EigenPack, par().eigenPack); + auto &noise = envGet(NoiseTensor, sNoiseName); + auto &perambulator = envGet(Perambulator, getName()); + auto &epack = envGet(LapEvecs, sLapEvecName); auto &unsmeared_sink = envGet(std::vector, getName() + "_unsmeared_sink"); diff --git a/Hadrons/Modules/MIO/LoadPerambulator.hpp b/Hadrons/Modules/MIO/LoadPerambulator.hpp index 6f7c5283..156fb5dc 100644 --- a/Hadrons/Modules/MIO/LoadPerambulator.hpp +++ b/Hadrons/Modules/MIO/LoadPerambulator.hpp @@ -101,20 +101,16 @@ std::vector TLoadPerambulator::getOutput(void) template void TLoadPerambulator::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 sIndexNames{"Nt", "nvec", "LI", "nnoise", "Nt_inv", "SI"}; - envCreate(MDistil::Perambulator, getName(), 1, - sIndexNames,Nt,nvec,Distil.LI,Distil.nnoise,Distil.Nt_inv,Distil.SI); + DISTIL_PARAMETERS_DEFINE( true ); + //std::array sIndexNames{"Nt", "nvec", "LI", "nnoise", "Nt_inv", "SI"}; + envCreate(MDistil::Perambulator, getName(), 1, MDistil::PerambIndexNames,Nt,nvec,LI,nnoise,Nt_inv,SI); } // execution /////////////////////////////////////////////////////////////////// template void TLoadPerambulator::execute(void) { - auto &perambulator = envGet(MDistil::Perambulator, getName()); + auto &perambulator = envGet(MDistil::Perambulator, getName()); const std::string sPerambName{par().PerambFileName + "." + std::to_string(vm().getTrajectory())}; perambulator.read(sPerambName.c_str()); } diff --git a/Hadrons/modules.inc b/Hadrons/modules.inc index d87aad20..c13c18a2 100644 --- a/Hadrons/modules.inc +++ b/Hadrons/modules.inc @@ -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 \ diff --git a/tests/hadrons/Test_distil.cc b/tests/hadrons/Test_distil.cc index a791858e..bc4578f1 100644 --- a/tests/hadrons/Test_distil.cc +++ b/tests/hadrons/Test_distil.cc @@ -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(szGaugeName); + std::string sGaugeName{szModuleName}; + sGaugeName.append( "_gauge" ); + application.createModule(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("LapEvec",p); + application.createModule(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(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("Peramb",PerambPar); + application.createModule( 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(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("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("Peramb2",SolvePar); SolvePar.PerambFileName="Peramb3"; SolvePar.nvec_reduced=3; SolvePar.LI_reduced=3; application.createModule("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("DistilVecs2",DistilVecPar); - DistilVecPar.perambulator="Peramb3"; - DistilVecPar.LI=3; - DistilVecPar.nvec=3; - application.createModule("DistilVecs3",DistilVecPar); - //DistilVecPar.perambulator="Peramb5_perambulator_light"; - DistilVecPar.perambulator="Peramb5"; - DistilVecPar.LI=5; - DistilVecPar.nvec=5; - application.createModule("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("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("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("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("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("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("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("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 p{as,2,7,2}; + Grid::Hadrons::MDistil::NamedTensor 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 sIndexNames{"Nt", "nvec", "LI", "nnoise", "Nt_inv", "SI"}; - Grid::Hadrons::MDistil::Perambulator 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 &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;