From df586a142d4a10ddffaad9bbb74caaa9378e9938 Mon Sep 17 00:00:00 2001 From: ferben Date: Mon, 11 Nov 2019 17:29:55 +0000 Subject: [PATCH] added DistilPar-module and cleaned up some code --- Hadrons/Modules.hpp | 1 + Hadrons/Modules/MDistil/DistilCommon.hpp | 36 +---- Hadrons/Modules/MDistil/DistilPar.cc | 7 + Hadrons/Modules/MDistil/DistilPar.hpp | 96 ++++++++++++++ Hadrons/Modules/MDistil/DistilVectors.hpp | 138 ++++++++++---------- Hadrons/Modules/MDistil/Noises.hpp | 23 ++-- Hadrons/Modules/MDistil/PerambFromSolve.hpp | 72 ++++++---- Hadrons/Modules/MDistil/Perambulator.hpp | 90 +++++++------ Hadrons/Modules/MIO/LoadDistilNoise.hpp | 9 +- Hadrons/Modules/MIO/LoadPerambulator.hpp | 17 ++- Hadrons/modules.inc | 3 + 11 files changed, 302 insertions(+), 190 deletions(-) create mode 100644 Hadrons/Modules/MDistil/DistilPar.cc create mode 100644 Hadrons/Modules/MDistil/DistilPar.hpp diff --git a/Hadrons/Modules.hpp b/Hadrons/Modules.hpp index c3c0cf86..e0e8832e 100644 --- a/Hadrons/Modules.hpp +++ b/Hadrons/Modules.hpp @@ -70,6 +70,7 @@ #include #include #include +#include #include #include #include diff --git a/Hadrons/Modules/MDistil/DistilCommon.hpp b/Hadrons/Modules/MDistil/DistilCommon.hpp index da46b918..72ee7449 100644 --- a/Hadrons/Modules/MDistil/DistilCommon.hpp +++ b/Hadrons/Modules/MDistil/DistilCommon.hpp @@ -63,42 +63,14 @@ BEGIN_MODULE_NAMESPACE(MDistil) struct DistilParameters: Serializable { GRID_SERIALIZABLE_CLASS_MEMBERS(DistilParameters, + int, nvec, int, nnoise, int, tsrc, - std::string, TI, - std::string, LI, - std::string, SI ) - DistilParameters() = default; - template DistilParameters(Reader& Reader){read(Reader,"Distil",*this);} - - // Numeric parameter is allowed to be empty (in which case it = Default), - // but assert during setup() if specified but not numeric - - static int ParameterDefault( const std::string & s, int Default, bool bCalledFromSetup ) - { - int i = Default; - if( s.length() > 0 ) { - std::istringstream ss( s ); - ss >> i; - if( bCalledFromSetup ) - assert( !ss.fail() && "Parameter should either be empty or integer" ); - } - return i; - } + int, TI, + int, LI, + int, SI ) }; -#define DISTIL_PARAMETERS_DEFINE( inSetup ) \ -const int Nt{env().getDim(Tdir)}; \ -const int nvec{par().nvec}; \ -const int nnoise{par().Distil.nnoise}; \ -const int tsrc{par().Distil.tsrc}; \ -const int TI{Hadrons::MDistil::DistilParameters::ParameterDefault(par().Distil.TI, Nt, inSetup)}; \ -const int LI{Hadrons::MDistil::DistilParameters::ParameterDefault(par().Distil.LI, nvec, inSetup)}; \ -const int SI{Hadrons::MDistil::DistilParameters::ParameterDefault(par().Distil.SI, Ns, inSetup)}; \ -const bool full_tdil{ TI == Nt }; \ -const bool exact_distillation{ full_tdil && LI == nvec }; \ -const int Nt_inv{ full_tdil ? 1 : TI } - /****************************************************************************** Make a lower dimensional grid in preparation for local slice operations ******************************************************************************/ diff --git a/Hadrons/Modules/MDistil/DistilPar.cc b/Hadrons/Modules/MDistil/DistilPar.cc new file mode 100644 index 00000000..f0c03ec9 --- /dev/null +++ b/Hadrons/Modules/MDistil/DistilPar.cc @@ -0,0 +1,7 @@ +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MDistil; + +template class Grid::Hadrons::MDistil::TDistilPar; diff --git a/Hadrons/Modules/MDistil/DistilPar.hpp b/Hadrons/Modules/MDistil/DistilPar.hpp new file mode 100644 index 00000000..328280bc --- /dev/null +++ b/Hadrons/Modules/MDistil/DistilPar.hpp @@ -0,0 +1,96 @@ +#ifndef Hadrons_MDistil_DistilPar_hpp_ +#define Hadrons_MDistil_DistilPar_hpp_ + +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * DistilPar * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MDistil) + + +class DistilParPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(DistilParPar, + int, nnoise, + int, tsrc, + int, TI, + int, LI, + int, SI ); +}; + +template +class TDistilPar: public Module +{ +public: + // constructor + TDistilPar(const std::string name); + // destructor + virtual ~TDistilPar(void) {}; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_TMP(DistilPar, TDistilPar, MDistil); + +/****************************************************************************** + * TDistilPar implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TDistilPar::TDistilPar(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TDistilPar::getInput(void) +{ + std::vector in; + + return in; +} + +template +std::vector TDistilPar::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TDistilPar::setup(void) +{ + +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TDistilPar::execute(void) +{ + Hadrons::MDistil::DistilParameters &out = envGet(Hadrons::MDistil::DistilParameters, getName()); + out.nnoise=par().nnoise; + out.tsrc=par().tsrc; + out.TI=par().TI; + out.LI=par().LI; + out.SI=par().SI; +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MDistil_DistilPar_hpp_ diff --git a/Hadrons/Modules/MDistil/DistilVectors.hpp b/Hadrons/Modules/MDistil/DistilVectors.hpp index 643e7d63..3f9909b7 100644 --- a/Hadrons/Modules/MDistil/DistilVectors.hpp +++ b/Hadrons/Modules/MDistil/DistilVectors.hpp @@ -47,11 +47,9 @@ public: std::string, noise, std::string, perambulator, std::string, lapevec, - std::string, source, - std::string, sink, - int, tsrc, - std::string, nvec, - std::string, TI) + std::string, rho, + std::string, phi, + MDistil::DistilParameters, DistilPar); }; template @@ -80,10 +78,10 @@ public: std::string PerambulatorName; std::string NoiseVectorName; std::string LapEvecName; - bool bMakeSource; - bool bMakeSink; - std::string SourceName; - std::string SinkName; + bool bMakeRho; + bool bMakePhi; + std::string RhoName; + std::string PhiName; }; MODULE_REGISTER_TMP(DistilVectors, TDistilVectors, MDistil); @@ -131,24 +129,24 @@ std::vector TDistilVectors::getInput(void) template std::vector TDistilVectors::getOutput(void) { - SourceName = par().source; - SinkName = par().sink; - bMakeSource = ( SourceName.size() > 0 ); - bMakeSink = ( SinkName.size() > 0 ); - if (!bMakeSource && !bMakeSink) + RhoName = par().rho; + PhiName = par().phi; + bMakeRho = ( RhoName.size() > 0 ); + bMakePhi = ( PhiName.size() > 0 ); + if (!bMakeRho && !bMakePhi) { - SourceName = getName(); - SinkName = SourceName; - SourceName.append("_rho"); - SinkName.append("_phi"); - bMakeSource = true; - bMakeSink = true; + RhoName = getName(); + PhiName = RhoName; + RhoName.append("_rho"); + PhiName.append("_phi"); + bMakeRho = true; + bMakePhi = true; } std::vector out; - if (bMakeSource) - out.push_back(SourceName); - if (bMakeSink) - out.push_back(SinkName); + if (bMakeRho) + out.push_back(RhoName); + if (bMakePhi) + out.push_back(PhiName); return out; } @@ -157,28 +155,25 @@ template void TDistilVectors::setup(void) { Cleanup(); - auto &noise = envGet(NoiseTensor, NoiseVectorName); + auto &noise = envGet(NoiseTensor, NoiseVectorName); auto &perambulator = envGet(PerambTensor, PerambulatorName); // We expect the perambulator to have been created with these indices assert( perambulator.ValidateIndexNames() && "Perambulator index names bad" ); - const int Nt{ env().getDim(Tdir) }; - assert( Nt == static_cast( perambulator.tensor.dimension(0) ) && "PerambTensor time dimensionality bad" ); - 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.tensor.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" ); + const int Nt{env().getDim(Tdir)}; + const int nvec{par().DistilPar.nvec}; + const int nnoise{par().DistilPar.nnoise}; + const int TI{par().DistilPar.TI}; + const int LI{par().DistilPar.LI}; + const int SI{par().DistilPar.SI}; + const bool full_tdil{ TI == Nt }; + const int Nt_inv{ full_tdil ? 1 : TI }; - if (bMakeSource) - envCreate(std::vector, SourceName, 1, nnoise*LI*SI*Nt_inv, envGetGrid(FermionField)); - if (bMakeSink) - envCreate(std::vector, SinkName, 1, nnoise*LI*SI*Nt_inv, envGetGrid(FermionField)); + if (bMakeRho) + envCreate(std::vector, RhoName, 1, nnoise*LI*SI*Nt_inv, envGetGrid(FermionField)); + if (bMakePhi) + envCreate(std::vector, PhiName, 1, nnoise*LI*SI*Nt_inv, envGetGrid(FermionField)); grid4d = env().getGrid(); Coordinate latt_size = GridDefaultLatt(); @@ -189,10 +184,10 @@ void TDistilVectors::setup(void) mpi_layout[Nd-1] = 1; grid3d = MakeLowerDimGrid(grid4d); - envTmp(LatticeSpinColourVector, "tmp2",1,LatticeSpinColourVector(grid4d)); - envTmp(LatticeSpinColourVector, "tmp3d",1,LatticeSpinColourVector(grid3d)); - envTmp(LatticeColourVector, "tmp3d_nospin",1,LatticeColourVector(grid3d)); - envTmp(LatticeSpinColourVector, "sink_tslice",1,LatticeSpinColourVector(grid3d)); + envTmp(LatticeSpinColourVector, "source4d",1,LatticeSpinColourVector(grid4d)); + envTmp(LatticeSpinColourVector, "source3d",1,LatticeSpinColourVector(grid3d)); + envTmp(LatticeColourVector, "source3d_nospin",1,LatticeColourVector(grid3d)); + envTmp(LatticeSpinColourVector, "sink3d",1,LatticeSpinColourVector(grid3d)); envTmp(LatticeColourVector, "evec3d",1,LatticeColourVector(grid3d)); } @@ -216,50 +211,49 @@ void TDistilVectors::execute(void) auto &perambulator = envGet(PerambTensor, PerambulatorName); auto &epack = envGet(Grid::Hadrons::EigenPack, LapEvecName); - envGetTmp(LatticeSpinColourVector, tmp2); - envGetTmp(LatticeSpinColourVector, tmp3d); - envGetTmp(LatticeColourVector, tmp3d_nospin); - envGetTmp(LatticeSpinColourVector, sink_tslice); + envGetTmp(LatticeSpinColourVector, source4d); + envGetTmp(LatticeSpinColourVector, source3d); + envGetTmp(LatticeColourVector, source3d_nospin); + envGetTmp(LatticeSpinColourVector, sink3d); envGetTmp(LatticeColourVector, evec3d); const int Ntlocal{ grid4d->LocalDimensions()[3] }; const int Ntfirst{ grid4d->LocalStarts()[3] }; - 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 }; + const int Nt{env().getDim(Tdir)}; + const int nvec{par().DistilPar.nvec}; + const int nnoise{par().DistilPar.nnoise}; + const int tsrc{par().DistilPar.tsrc}; + const int TI{par().DistilPar.TI}; + const int LI{par().DistilPar.LI}; + const int SI{par().DistilPar.SI}; + const bool full_tdil{ TI == Nt }; + const int Nt_inv{ full_tdil ? 1 : TI }; int vecindex; int t_inv; - if (bMakeSource) + if (bMakeRho) { - auto &rho = envGet(std::vector, SourceName); + auto &rho = envGet(std::vector, RhoName); 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] = 0; - tmp3d_nospin = 0; + source3d_nospin = 0; 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-Ntfirst,Tdir); - tmp3d_nospin = evec3d * noise.tensor(inoise, t_inv, ik, is); - tmp3d=0; - pokeSpin(tmp3d,tmp3d_nospin,is); - tmp2=0; - InsertSliceLocal(tmp3d,tmp2,0,t_inv-Ntfirst,Tdir); - rho[vecindex] += tmp2; + source3d_nospin = evec3d * noise.tensor(inoise, t_inv, ik, is); + source3d=0; + pokeSpin(source3d,source3d_nospin,is); + source4d=0; + InsertSliceLocal(source3d,source4d,0,t_inv-Ntfirst,Tdir); + rho[vecindex] += source4d; } } } @@ -269,8 +263,8 @@ void TDistilVectors::execute(void) } } } - if (bMakeSink) { - auto &phi = envGet(std::vector, SinkName); + if (bMakePhi) { + auto &phi = envGet(std::vector, PhiName); for (int inoise = 0; inoise < nnoise; inoise++) { for (int dk = 0; dk < LI; dk++) { for (int dt = 0; dt < Nt_inv; dt++) { @@ -278,12 +272,12 @@ void TDistilVectors::execute(void) vecindex = inoise + nnoise * dk + nnoise * LI * ds + nnoise *LI * SI*dt; phi[vecindex] = 0; for (int t = Ntfirst; t < Ntfirst + Ntlocal; t++) { - sink_tslice=0; + sink3d=0; for (int ivec = 0; ivec < nvec; ivec++) { ExtractSliceLocal(evec3d,epack.evec[ivec],0,t-Ntfirst,Tdir); - sink_tslice += evec3d * perambulator.tensor(t, ivec, dk, inoise,dt,ds); + sink3d += evec3d * perambulator.tensor(t, ivec, dk, inoise,dt,ds); } - InsertSliceLocal(sink_tslice,phi[vecindex],0,t-Ntfirst,Tdir); + InsertSliceLocal(sink3d,phi[vecindex],0,t-Ntfirst,Tdir); } } } diff --git a/Hadrons/Modules/MDistil/Noises.hpp b/Hadrons/Modules/MDistil/Noises.hpp index c17a29fa..1af023e4 100644 --- a/Hadrons/Modules/MDistil/Noises.hpp +++ b/Hadrons/Modules/MDistil/Noises.hpp @@ -43,10 +43,7 @@ class NoisesPar: Serializable { public: GRID_SERIALIZABLE_CLASS_MEMBERS(NoisesPar, - int, nnoise, - int, nvec, - std::string, TI, - std::string, LI, + MDistil::DistilParameters, DistilPar, std::string, NoiseFileName) }; @@ -96,9 +93,9 @@ std::vector TNoises::getOutput(void) template void TNoises::setup(void) { - const int Nt{env().getDim(Tdir)}; - const int nnoise{par().nnoise}; - const int nvec{par().nvec}; + const int Nt{env().getDim(Tdir)}; + const int nvec{par().DistilPar.nvec}; + const int nnoise{par().DistilPar.nnoise}; envCreate(NoiseTensor, getName(), 1, nnoise, Nt, nvec, Ns); } @@ -107,12 +104,12 @@ template void TNoises::execute(void) { const int Nt{env().getDim(Tdir)}; - const int nnoise{par().nnoise}; - const int nvec{par().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 }; \ + const int nnoise{par().DistilPar.nnoise}; + const int nvec{par().DistilPar.nvec}; + const int TI{par().DistilPar.TI}; + const int LI{par().DistilPar.LI}; + const bool full_tdil{ TI == Nt }; + const bool exact_distillation{ full_tdil && LI == nvec }; // We use our own seeds so we can specify different noises per quark Real rn; diff --git a/Hadrons/Modules/MDistil/PerambFromSolve.hpp b/Hadrons/Modules/MDistil/PerambFromSolve.hpp index 14340408..662055a0 100644 --- a/Hadrons/Modules/MDistil/PerambFromSolve.hpp +++ b/Hadrons/Modules/MDistil/PerambFromSolve.hpp @@ -51,10 +51,9 @@ public: std::string, eigenPack, std::string, PerambFileName, std::string, solve, - int, nvec, std::string, nvec_reduced, std::string, LI_reduced, - DistilParameters, Distil); + MDistil::DistilParameters, DistilPar); }; template @@ -116,16 +115,23 @@ template void TPerambFromSolve::setup(void) { Cleanup(); - DISTIL_PARAMETERS_DEFINE( true ); - const int nvec_reduced{ Hadrons::MDistil::DistilParameters::ParameterDefault( par().nvec_reduced, nvec, true) }; - const int LI_reduced{ Hadrons::MDistil::DistilParameters::ParameterDefault( par().LI_reduced, LI, true) }; + const int Nt{env().getDim(Tdir)}; + const int nvec{par().DistilPar.nvec}; + const int nnoise{par().DistilPar.nnoise}; + const int LI{par().DistilPar.LI}; + const int TI{par().DistilPar.TI}; + const int SI{par().DistilPar.SI}; + const bool full_tdil{ TI == Nt }; + const int Nt_inv{ full_tdil ? 1 : TI }; + const int nvec_reduced{ par().nvec_reduced.empty() ? nvec:std::stoi(par().nvec_reduced)}; + const int LI_reduced{ par().LI_reduced.empty() ? LI:std::stoi(par().LI_reduced)}; grid4d = env().getGrid(); grid3d = MakeLowerDimGrid(grid4d); envCreate(PerambTensor, getName(), 1, Nt,nvec_reduced,LI_reduced,nnoise,Nt_inv,SI); envCreate(NoiseTensor, getName() + "_noise", 1, nnoise, Nt, nvec, Ns ); - envTmp(LatticeColourVector, "result_3d",1,LatticeColourVector(grid3d)); + envTmp(LatticeColourVector, "result3d_nospin",1,LatticeColourVector(grid3d)); envTmp(LatticeColourVector, "evec3d",1,LatticeColourVector(grid3d)); - envTmpLat(LatticeColourVector, "result_nospin"); + envTmpLat(LatticeColourVector, "result4d_nospin"); } template @@ -146,28 +152,42 @@ void TPerambFromSolve::execute(void) GridCartesian * grid4d = env().getGrid(); const int Ntlocal{grid4d->LocalDimensions()[3]}; const int Ntfirst{grid4d->LocalStarts()[3]}; - DISTIL_PARAMETERS_DEFINE( false ); - const int nvec_reduced{ Hadrons::MDistil::DistilParameters::ParameterDefault( par().nvec_reduced, nvec, false) }; - const int LI_reduced{ Hadrons::MDistil::DistilParameters::ParameterDefault( par().LI_reduced, LI, false) }; - auto &perambulator = envGet(PerambTensor, getName()); - auto &solve = envGet(std::vector, par().solve); - auto &epack = envGet(Grid::Hadrons::EigenPack, par().eigenPack); + const int Nt{env().getDim(Tdir)}; + const int nvec{par().DistilPar.nvec}; + const int nnoise{par().DistilPar.nnoise}; + const int TI{par().DistilPar.TI}; + const int LI{par().DistilPar.LI}; + const int SI{par().DistilPar.SI}; + const bool full_tdil{ TI == Nt }; + const int Nt_inv{ full_tdil ? 1 : TI }; + const int nvec_reduced{ par().nvec_reduced.empty() ? nvec:std::stoi(par().nvec_reduced)}; + const int LI_reduced{ par().LI_reduced.empty() ? LI:std::stoi(par().LI_reduced)}; + auto &perambulator = envGet(PerambTensor, getName()); + auto &solve = envGet(std::vector, par().solve); + auto &epack = envGet(Grid::Hadrons::EigenPack, par().eigenPack); - envGetTmp(LatticeColourVector, result_nospin); - envGetTmp(LatticeColourVector, result_3d); + envGetTmp(LatticeColourVector, result4d_nospin); + envGetTmp(LatticeColourVector, result3d_nospin); envGetTmp(LatticeColourVector, evec3d); - for (int inoise = 0; inoise < nnoise; inoise++) { - for (int dk = 0; dk < LI_reduced; dk++) { - for (int dt = 0; dt < Nt_inv; dt++) { - for (int ds = 0; ds < SI; ds++) { - for (int is = 0; is < Ns; is++) { - result_nospin = peekSpin(solve[inoise+nnoise*(dk+LI*(dt+Nt_inv*ds))],is); - for (int t = Ntfirst; t < Ntfirst + Ntlocal; t++) { - ExtractSliceLocal(result_3d,result_nospin,0,t-Ntfirst,Tdir); - for (int ivec = 0; ivec < nvec_reduced; ivec++) { + for (int inoise = 0; inoise < nnoise; inoise++) + { + for (int dk = 0; dk < LI_reduced; dk++) + { + for (int dt = 0; dt < Nt_inv; dt++) + { + for (int ds = 0; ds < SI; ds++) + { + for (int is = 0; is < Ns; is++) + { + result4d_nospin = peekSpin(solve[inoise+nnoise*(dk+LI*(dt+Nt_inv*ds))],is); + for (int t = Ntfirst; t < Ntfirst + Ntlocal; t++) + { + ExtractSliceLocal(result3d_nospin,result4d_nospin,0,t-Ntfirst,Tdir); + for (int ivec = 0; ivec < nvec_reduced; ivec++) + { ExtractSliceLocal(evec3d,epack.evec[ivec],0,t-Ntfirst,Tdir); - pokeSpin(perambulator.tensor(t, ivec, dk, inoise,dt,ds),static_cast(innerProduct(evec3d, result_3d)),is); + pokeSpin(perambulator.tensor(t, ivec, dk, inoise,dt,ds),static_cast(innerProduct(evec3d, result3d_nospin)),is); LOG(Message) << "perambulator(t, ivec, dk, inoise,dt,ds)(is) = (" << t << "," << ivec << "," << dk << "," << inoise << "," << dt << "," << ds << ")(" << is << ") = " << perambulator.tensor(t, ivec, dk, inoise,dt,ds)()(is)() << std::endl; } } @@ -179,8 +199,6 @@ void TPerambFromSolve::execute(void) if(grid4d->IsBoss()) { std::string sPerambName{par().PerambFileName}; - if (sPerambName.empty()) - sPerambName = getName(); sPerambName.append( "." ); sPerambName.append( std::to_string(vm().getTrajectory())); perambulator.write(sPerambName.c_str()); diff --git a/Hadrons/Modules/MDistil/Perambulator.hpp b/Hadrons/Modules/MDistil/Perambulator.hpp index 17ea41a2..3cc444e3 100644 --- a/Hadrons/Modules/MDistil/Perambulator.hpp +++ b/Hadrons/Modules/MDistil/Perambulator.hpp @@ -49,8 +49,7 @@ public: std::string, PerambFileName, std::string, UnsmearedSinkFileName, std::string, UnsmearedSinkMultiFile, - int, nvec, - DistilParameters, Distil); + MDistil::DistilParameters, DistilPar); }; template @@ -106,8 +105,6 @@ std::vector TPerambulator::getInput(void) { sLapEvecName = par().lapevec; sNoiseName = par().noise; - if( sNoiseName.length() == 0 ) - sNoiseName = getName() + "_noise"; return {sLapEvecName, par().solver, sNoiseName }; } @@ -124,23 +121,32 @@ void TPerambulator::setup(void) Cleanup(); grid4d = env().getGrid(); grid3d = MakeLowerDimGrid(grid4d); - DISTIL_PARAMETERS_DEFINE( true ); + const int Nt{env().getDim(Tdir)}; + const int nvec{par().DistilPar.nvec}; + const int nnoise{par().DistilPar.nnoise}; + const int tsrc{par().DistilPar.tsrc}; + const int TI{par().DistilPar.TI}; + const int LI{par().DistilPar.LI}; + const int SI{par().DistilPar.SI}; + const bool full_tdil{ TI == Nt }; + const int Nt_inv{ full_tdil ? 1 : TI }; const std::string UnsmearedSinkFileName{ par().UnsmearedSinkFileName }; if (!UnsmearedSinkFileName.empty()) - bool bMulti = ( Hadrons::MDistil::DistilParameters::ParameterDefault( par().UnsmearedSinkMultiFile, 1, true ) != 0 ); + //bool bMulti = ( Hadrons::MDistil::DistilParameters::ParameterDefault( par().UnsmearedSinkMultiFile, 1, true ) != 0 ); + bool bMulti = 0; envCreate(PerambTensor, getName(), 1, Nt,nvec,LI,nnoise,Nt_inv,SI); envCreate(std::vector, getName() + "_unsmeared_sink", 1, nnoise*LI*Ns*Nt_inv, envGetGrid(FermionField)); - envTmpLat(LatticeSpinColourVector, "dist_source"); - envTmpLat(LatticeSpinColourVector, "tmp2"); - envTmpLat(LatticeSpinColourVector, "result"); - envTmpLat(LatticeColourVector, "result_nospin"); - envTmp(LatticeSpinColourVector, "tmp3d",1,LatticeSpinColourVector(grid3d)); - envTmp(LatticeColourVector, "tmp3d_nospin",1,LatticeColourVector(grid3d)); - envTmp(LatticeColourVector, "result_3d",1,LatticeColourVector(grid3d)); - envTmp(LatticeColourVector, "evec3d",1,LatticeColourVector(grid3d)); + envTmpLat(LatticeSpinColourVector, "dist_source"); + envTmpLat(LatticeSpinColourVector, "source4d"); + envTmp(LatticeSpinColourVector, "source3d",1,LatticeSpinColourVector(grid3d)); + envTmp(LatticeColourVector, "source3d_nospin",1,LatticeColourVector(grid3d)); + envTmpLat(LatticeSpinColourVector, "result4d"); + envTmpLat(LatticeColourVector, "result4d_nospin"); + envTmp(LatticeColourVector, "result3d_nospin",1,LatticeColourVector(grid3d)); + envTmp(LatticeColourVector, "evec3d",1,LatticeColourVector(grid3d)); Ls_ = env().getObjectLs(par().solver); envTmpLat(FermionField, "v4dtmp"); @@ -164,7 +170,16 @@ void TPerambulator::Cleanup(void) template void TPerambulator::execute(void) { - DISTIL_PARAMETERS_DEFINE( false ); + const int Nt{env().getDim(Tdir)}; + const int nvec{par().DistilPar.nvec}; + const int nnoise{par().DistilPar.nnoise}; + const int tsrc{par().DistilPar.tsrc}; + const int TI{par().DistilPar.TI}; + const int LI{par().DistilPar.LI}; + const int SI{par().DistilPar.SI}; + const bool full_tdil{ TI == Nt }; + const int Nt_inv{ full_tdil ? 1 : TI }; + auto &solver=envGet(Solver, par().solver); auto &mat = solver.getFMat(); envGetTmp(FermionField, v4dtmp); @@ -175,12 +190,12 @@ void TPerambulator::execute(void) auto &epack = envGet(LapEvecs, sLapEvecName); auto &unsmeared_sink = envGet(std::vector, getName() + "_unsmeared_sink"); envGetTmp(LatticeSpinColourVector, dist_source); - envGetTmp(LatticeSpinColourVector, tmp2); - envGetTmp(LatticeSpinColourVector, result); - envGetTmp(LatticeColourVector, result_nospin); - envGetTmp(LatticeSpinColourVector, tmp3d); - envGetTmp(LatticeColourVector, tmp3d_nospin); - envGetTmp(LatticeColourVector, result_3d); + envGetTmp(LatticeSpinColourVector, source4d); + envGetTmp(LatticeSpinColourVector, source3d); + envGetTmp(LatticeColourVector, source3d_nospin); + envGetTmp(LatticeSpinColourVector, result4d); + envGetTmp(LatticeColourVector, result4d_nospin); + envGetTmp(LatticeColourVector, result3d_nospin); envGetTmp(LatticeColourVector, evec3d); const int Ntlocal{grid4d->LocalDimensions()[3]}; const int Ntfirst{grid4d->LocalStarts()[3]}; @@ -196,7 +211,7 @@ void TPerambulator::execute(void) { LOG(Message) << "LapH source vector from noise " << inoise << " and dilution component (d_k,d_t,d_alpha) : (" << dk << ","<< dt << "," << ds << ")" << std::endl; dist_source = 0; - tmp3d_nospin = 0; + source3d_nospin = 0; evec3d = 0; for (int it = dt; it < Nt; it += TI) { @@ -208,39 +223,39 @@ void TPerambulator::execute(void) for (int is = ds; is < Ns; is += SI) { ExtractSliceLocal(evec3d,epack.evec[ik],0,t_inv-Ntfirst,Tdir); - tmp3d_nospin = evec3d * noise.tensor(inoise, t_inv, ik, is); - tmp3d=0; - pokeSpin(tmp3d,tmp3d_nospin,is); - tmp2=0; - InsertSliceLocal(tmp3d,tmp2,0,t_inv-Ntfirst,Tdir); - dist_source += tmp2; + source3d_nospin = evec3d * noise.tensor(inoise, t_inv, ik, is); + source3d=0; + pokeSpin(source3d,source3d_nospin,is); + source4d=0; + InsertSliceLocal(source3d,source4d,0,t_inv-Ntfirst,Tdir); + dist_source += source4d; } } } } - result=0; + result4d=0; v4dtmp = dist_source; if (Ls_ == 1) - solver(result, v4dtmp); + solver(result4d, v4dtmp); else { mat.ImportPhysicalFermionSource(v4dtmp, v5dtmp); solver(v5dtmp_sol, v5dtmp); mat.ExportPhysicalFermionSolution(v5dtmp_sol, v4dtmp); - result = v4dtmp; + result4d = v4dtmp; } if (!UnsmearedSinkFileName.empty()) - unsmeared_sink[inoise+nnoise*(dk+LI*(dt+Nt_inv*ds))] = result; + unsmeared_sink[inoise+nnoise*(dk+LI*(dt+Nt_inv*ds))] = result4d; for (int is = 0; is < Ns; is++) { - result_nospin = peekSpin(result,is); + result4d_nospin = peekSpin(result4d,is); for (int t = Ntfirst; t < Ntfirst + Ntlocal; t++) { - ExtractSliceLocal(result_3d,result_nospin,0,t-Ntfirst,Tdir); - for (int ivec = 0; ivec < nvec; ivec++) + ExtractSliceLocal(result3d_nospin,result4d_nospin,0,t-Ntfirst,Tdir); + for (int ivec = 0; ivec < nvec; ivec++) { ExtractSliceLocal(evec3d,epack.evec[ivec],0,t-Ntfirst,Tdir); - pokeSpin(perambulator.tensor(t, ivec, dk, inoise,dt,ds),static_cast(innerProduct(evec3d, result_3d)),is); + pokeSpin(perambulator.tensor(t, ivec, dk, inoise,dt,ds),static_cast(innerProduct(evec3d, result3d_nospin)),is); } } } @@ -284,7 +299,8 @@ void TPerambulator::execute(void) //Save the unsmeared sinks if filename specified if (!UnsmearedSinkFileName.empty()) { - bool bMulti = ( Hadrons::MDistil::DistilParameters::ParameterDefault( par().UnsmearedSinkMultiFile, 1, false ) != 0 ); + // bool bMulti = ( Hadrons::MDistil::DistilParameters::ParameterDefault( par().UnsmearedSinkMultiFile, 1, false ) != 0 ); + bool bMulti =0; LOG(Message) << "Writing unsmeared sink to " << UnsmearedSinkFileName << std::endl; A2AVectorsIo::write(UnsmearedSinkFileName, unsmeared_sink, bMulti, vm().getTrajectory()); } diff --git a/Hadrons/Modules/MIO/LoadDistilNoise.hpp b/Hadrons/Modules/MIO/LoadDistilNoise.hpp index 27b7e76a..a1a20ca1 100644 --- a/Hadrons/Modules/MIO/LoadDistilNoise.hpp +++ b/Hadrons/Modules/MIO/LoadDistilNoise.hpp @@ -44,8 +44,7 @@ class LoadDistilNoisePar: Serializable public: GRID_SERIALIZABLE_CLASS_MEMBERS(LoadDistilNoisePar, std::string, NoiseFileName, - int, nvec, - MDistil::DistilParameters, Distil); + MDistil::DistilParameters, DistilPar); }; template @@ -97,8 +96,10 @@ std::vector TLoadDistilNoise::getOutput(void) template void TLoadDistilNoise::setup(void) { - DISTIL_PARAMETERS_DEFINE( true ); - envCreate(MDistil::NoiseTensor, getName(), 1, nnoise, Nt, nvec, Ns); + const int Nt{env().getDim(Tdir)}; + const int nvec{par().DistilPar.nvec}; + const int nnoise{par().DistilPar.nnoise}; + envCreate(MDistil::NoiseTensor, getName(), 1, nnoise, Nt, nvec, Ns); } // execution /////////////////////////////////////////////////////////////////// diff --git a/Hadrons/Modules/MIO/LoadPerambulator.hpp b/Hadrons/Modules/MIO/LoadPerambulator.hpp index 23244512..45f31b2d 100644 --- a/Hadrons/Modules/MIO/LoadPerambulator.hpp +++ b/Hadrons/Modules/MIO/LoadPerambulator.hpp @@ -43,9 +43,8 @@ class LoadPerambulatorPar: Serializable { public: GRID_SERIALIZABLE_CLASS_MEMBERS(LoadPerambulatorPar, - std::string, PerambFileName, //stem!!! - int, nvec, - MDistil::DistilParameters, Distil); + std::string, PerambFileName, + MDistil::DistilParameters, DistilPar); }; template @@ -97,8 +96,16 @@ std::vector TLoadPerambulator::getOutput(void) template void TLoadPerambulator::setup(void) { - DISTIL_PARAMETERS_DEFINE( true ); - //std::array sIndexNames{"Nt", "nvec", "LI", "nnoise", "Nt_inv", "SI"}; + const int Nt{env().getDim(Tdir)}; + const int nvec{par().DistilPar.nvec}; + const int nnoise{par().DistilPar.nnoise}; + const int tsrc{par().DistilPar.tsrc}; + const int TI{par().DistilPar.TI}; + const int LI{par().DistilPar.LI}; + const int SI{par().DistilPar.SI}; + const bool full_tdil{ TI == Nt }; + const bool exact_distillation{ full_tdil && LI == nvec }; + const int Nt_inv{ full_tdil ? 1 : TI }; envCreate(MDistil::PerambTensor, getName(), 1, Nt,nvec,LI,nnoise,Nt_inv,SI); } diff --git a/Hadrons/modules.inc b/Hadrons/modules.inc index f3e43d28..30a7cc6d 100644 --- a/Hadrons/modules.inc +++ b/Hadrons/modules.inc @@ -65,7 +65,9 @@ modules_cc =\ Modules/MScalarSUN/TrMag.cc \ Modules/MDistil/PerambFromSolve.cc \ Modules/MDistil/DistilVectors.cc \ + Modules/MDistil/DistilPar.cc \ Modules/MDistil/LapEvec.cc \ + Modules/MDistil/DistilParameters.cc \ Modules/MDistil/Perambulator.cc \ Modules/MDistil/Noises.cc \ Modules/MAction/MobiusDWF.cc \ @@ -151,6 +153,7 @@ modules_hpp =\ Modules/MDistil/Noises.hpp \ Modules/MDistil/LapEvec.hpp \ Modules/MDistil/Perambulator.hpp \ + Modules/MDistil/DistilPar.hpp \ Modules/MDistil/DistilCommon.hpp \ Modules/MAction/ZMobiusDWF.hpp \ Modules/MAction/ScaledDWF.hpp \