diff --git a/Grid/communicator/SharedMemoryMPI.cc b/Grid/communicator/SharedMemoryMPI.cc index 4af7c402..ed465252 100644 --- a/Grid/communicator/SharedMemoryMPI.cc +++ b/Grid/communicator/SharedMemoryMPI.cc @@ -162,11 +162,8 @@ static inline int divides(int a,int b) void GlobalSharedMemory::GetShmDims(const Coordinate &WorldDims,Coordinate &ShmDims) { //////////////////////////////////////////////////////////////// - // Assert power of two shm_size. + // Powers of 2,3,5 only in prime decomposition for now //////////////////////////////////////////////////////////////// - int log2size = Log2Size(WorldShmSize,MAXLOG2RANKSPERNODE); - assert(log2size != -1); - int ndimension = WorldDims.size(); ShmDims=Coordinate(ndimension,1); @@ -177,7 +174,8 @@ void GlobalSharedMemory::GetShmDims(const Coordinate &WorldDims,Coordinate &ShmD while(AutoShmSize != WorldShmSize) { for(int p=0;p ranks(size); for(int r=0;r #include #endif #ifdef __x86_64__ +#ifdef GRID_NVCC +accelerator_inline uint64_t __rdtsc(void) { return 0; } +accelerator_inline uint64_t __rdpmc(int ) { return 0; } +#else #include #endif +#endif NAMESPACE_BEGIN(Grid); @@ -89,13 +94,8 @@ inline uint64_t cyclecount(void){ return tmp; } #elif defined __x86_64__ -#ifdef GRID_NVCC -accelerator_inline uint64_t __rdtsc(void) { return 0; } -#endif inline uint64_t cyclecount(void){ return __rdtsc(); - // unsigned int dummy; - // return __rdtscp(&dummy); } #else diff --git a/Grid/qcd/action/fermion/instantiation/ImprovedStaggeredFermionInstantiation.cc b/Grid/qcd/action/fermion/instantiation/ImprovedStaggeredFermionInstantiation.cc index a617f6cb..b1afc20e 100644 --- a/Grid/qcd/action/fermion/instantiation/ImprovedStaggeredFermionInstantiation.cc +++ b/Grid/qcd/action/fermion/instantiation/ImprovedStaggeredFermionInstantiation.cc @@ -26,7 +26,7 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#include +#include NAMESPACE_BEGIN(Grid); diff --git a/Grid/qcd/action/fermion/instantiation/ImprovedStaggeredFermionInstantiation.cc.master b/Grid/qcd/action/fermion/instantiation/ImprovedStaggeredFermionInstantiation.cc.master index 2023adc2..d35b7349 100644 --- a/Grid/qcd/action/fermion/instantiation/ImprovedStaggeredFermionInstantiation.cc.master +++ b/Grid/qcd/action/fermion/instantiation/ImprovedStaggeredFermionInstantiation.cc.master @@ -26,7 +26,7 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#include +#include #include NAMESPACE_BEGIN(Grid); diff --git a/Grid/qcd/action/fermion/instantiation/StaggeredImplD/ImprovedStaggeredFermionInstantiationStaggeredImplD.cc b/Grid/qcd/action/fermion/instantiation/StaggeredImplD/ImprovedStaggeredFermionInstantiationStaggeredImplD.cc index 2023adc2..d35b7349 100644 --- a/Grid/qcd/action/fermion/instantiation/StaggeredImplD/ImprovedStaggeredFermionInstantiationStaggeredImplD.cc +++ b/Grid/qcd/action/fermion/instantiation/StaggeredImplD/ImprovedStaggeredFermionInstantiationStaggeredImplD.cc @@ -26,7 +26,7 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#include +#include #include NAMESPACE_BEGIN(Grid); diff --git a/Grid/qcd/action/fermion/instantiation/StaggeredImplF/ImprovedStaggeredFermionInstantiationStaggeredImplF.cc b/Grid/qcd/action/fermion/instantiation/StaggeredImplF/ImprovedStaggeredFermionInstantiationStaggeredImplF.cc index 2023adc2..d35b7349 100644 --- a/Grid/qcd/action/fermion/instantiation/StaggeredImplF/ImprovedStaggeredFermionInstantiationStaggeredImplF.cc +++ b/Grid/qcd/action/fermion/instantiation/StaggeredImplF/ImprovedStaggeredFermionInstantiationStaggeredImplF.cc @@ -26,7 +26,7 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#include +#include #include NAMESPACE_BEGIN(Grid); diff --git a/Grid/qcd/smearing/StoutSmearing.h b/Grid/qcd/smearing/StoutSmearing.h index f463dc81..ed2ccdb6 100644 --- a/Grid/qcd/smearing/StoutSmearing.h +++ b/Grid/qcd/smearing/StoutSmearing.h @@ -1,5 +1,34 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/qcd/smearing/StoutSmearing.h + + Copyright (C) 2019 + + Author: unknown + 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 + *************************************************************************************/ /* - @file stoutSmear.hpp + @file StoutSmearing.h @brief Declares Stout smearing class */ #pragma once @@ -9,19 +38,43 @@ NAMESPACE_BEGIN(Grid); /*! @brief Stout smearing of link variable. */ template class Smear_Stout : public Smear { -private: - const Smear* SmearBase; + private: + int OrthogDim = -1; + const std::vector SmearRho; + // Smear* ownership semantics: + // Smear* passed in to constructor are owned by caller, so we don't delete them here + // Smear* created within constructor need to be deleted as part of the destructor + const std::unique_ptr> OwnedBase; // deleted at destruction + const Smear* SmearBase; // Not owned by this object, so not deleted at destruction + // only anticipated to be used from default constructor + inline static std::vector rho3D(double rho, int orthogdim){ + std::vector rho3d(Nd*Nd); + for (int mu=0; mu* base) : SmearBase(base) { - assert(Nc == 3);// "Stout smearing currently implemented only for Nc==3"); + /*! Stout smearing with base explicitly specified */ + Smear_Stout(Smear* base) : SmearBase{base} { + assert(Nc == 3 && "Stout smearing currently implemented only for Nc==3"); } - /*! Default constructor */ - Smear_Stout(double rho = 1.0) : SmearBase(new Smear_APE(rho)) { - assert(Nc == 3);// "Stout smearing currently implemented only for Nc==3"); + /*! Construct stout smearing object from explicitly specified rho matrix */ + Smear_Stout(const std::vector& rho_) + : OwnedBase{new Smear_APE(rho_)}, SmearBase{OwnedBase.get()} { + std::cout << GridLogDebug << "Stout smearing constructor : Smear_Stout(const std::vector& " << rho_ << " )" << std::endl + assert(Nc == 3 && "Stout smearing currently implemented only for Nc==3"); + } + + /*! Default constructor. rho is constant in all directions, optionally except for orthogonal dimension */ + Smear_Stout(double rho = 1.0, int orthogdim = -1) + : OrthogDim{orthogdim}, SmearRho{ rho3D(rho,orthogdim) }, OwnedBase{ new Smear_APE(SmearRho) }, SmearBase{OwnedBase.get()} { + assert(Nc == 3 && "Stout smearing currently implemented only for Nc==3"); } ~Smear_Stout() {} // delete SmearBase... @@ -36,12 +89,16 @@ public: SmearBase->smear(C, U); for (int mu = 0; mu < Nd; mu++) { - tmp = peekLorentz(C, mu); - Umu = peekLorentz(U, mu); - iq_mu = Ta( - tmp * - adj(Umu)); // iq_mu = Ta(Omega_mu) to match the signs with the paper - exponentiate_iQ(tmp, iq_mu); + if( mu == OrthogDim ) + tmp = 1.0; // Don't smear in the orthogonal direction + else { + tmp = peekLorentz(C, mu); + Umu = peekLorentz(U, mu); + iq_mu = Ta( + tmp * + adj(Umu)); // iq_mu = Ta(Omega_mu) to match the signs with the paper + exponentiate_iQ(tmp, iq_mu); + } pokeLorentz(u_smr, tmp * Umu, mu); // u_smr = exp(iQ_mu)*U_mu } std::cout << GridLogDebug << "Stout smearing completed\n"; @@ -80,6 +137,7 @@ public: iQ2 = iQ * iQ; iQ3 = iQ * iQ2; + //We should check sgn(c0) here already and then apply eq (34) from 0311018 set_uw(u, w, iQ2, iQ3); set_fj(f0, f1, f2, u, w); @@ -139,9 +197,8 @@ public: } LatticeComplex func_xi0(const LatticeComplex& w) const { - // Define a function to do the check - // if( w < 1e-4 ) std::cout << GridLogWarning<< "[Smear_stout] w too small: - // "<< w <<"\n"; + // Definition from arxiv 0311018 + //if (abs(w) < 0.05) {w2 = w*w; return 1.0 - w2/6.0 * (1.0-w2/20.0 * (1.0-w2/42.0));} return sin(w) / w; } @@ -154,4 +211,3 @@ public: }; NAMESPACE_END(Grid); - diff --git a/Grid/qcd/utils/A2Autils.h b/Grid/qcd/utils/A2Autils.h index b603eede..c7c7d329 100644 --- a/Grid/qcd/utils/A2Autils.h +++ b/Grid/qcd/utils/A2Autils.h @@ -260,7 +260,7 @@ void A2Autils::MesonField(TensorType &mat, int ij_dx = m+Nmom*i + Nmom*Lblock * j + Nmom*Lblock * Rblock * lt; for(int mu=0;mu> SpinMatrixField; typedef typename SpinMatrixField::vector_object sobj; - static constexpr int epsilon[6][3] = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}}; - static constexpr int epsilon_sgn[6]= {1,1,1,-1,-1,-1}; + static const int epsilon[6][3] ; + static const Complex epsilon_sgn[6]; private: template @@ -149,10 +149,15 @@ public: SpinMatrixField &stn_corr); }; -template -constexpr int BaryonUtils::epsilon[6][3]; -template -constexpr int BaryonUtils::epsilon_sgn[6]; +template +const int BaryonUtils::epsilon[6][3] = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}}; +template +const Complex BaryonUtils::epsilon_sgn[6] = {Complex(1), + Complex(1), + Complex(1), + Complex(-1), + Complex(-1), + Complex(-1)}; template template @@ -189,7 +194,7 @@ void BaryonUtils::baryon_site(const mobj &D1, for (int alpha_right=0; alpha_right(epsilon_sgn[ie_left] * epsilon_sgn[ie_right]) * pD1()(gamma_left,gamma_left)(c_right,c_left)*D2g()(alpha_right,beta_left)(a_right,a_left)*gD3()(alpha_right,beta_left)(b_right,b_left); + result()()() += epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1()(gamma_left,gamma_left)(c_right,c_left)*D2g()(alpha_right,beta_left)(a_right,a_left)*gD3()(alpha_right,beta_left)(b_right,b_left); }}} } //This is the \delta_{456}^{231} part @@ -198,7 +203,7 @@ void BaryonUtils::baryon_site(const mobj &D1, for (int alpha_right=0; alpha_right(epsilon_sgn[ie_left] * epsilon_sgn[ie_right]) * pD1g()(gamma_left,beta_left)(c_right,a_left)*D2()(alpha_right,beta_left)(a_right,b_left)*gD3()(alpha_right,gamma_left)(b_right,c_left); + result()()() += epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1g()(gamma_left,beta_left)(c_right,a_left)*D2()(alpha_right,beta_left)(a_right,b_left)*gD3()(alpha_right,gamma_left)(b_right,c_left); }}} } //This is the \delta_{456}^{312} part @@ -207,7 +212,7 @@ void BaryonUtils::baryon_site(const mobj &D1, for (int alpha_right=0; alpha_right(epsilon_sgn[ie_left] * epsilon_sgn[ie_right]) * pD1()(gamma_left,beta_left)(c_right,b_left)*D2()(alpha_right,gamma_left)(a_right,c_left)*gD3g()(alpha_right,beta_left)(b_right,a_left); + result()()() += epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1()(gamma_left,beta_left)(c_right,b_left)*D2()(alpha_right,gamma_left)(a_right,c_left)*gD3g()(alpha_right,beta_left)(b_right,a_left); }}} } //This is the \delta_{456}^{132} part @@ -216,7 +221,7 @@ void BaryonUtils::baryon_site(const mobj &D1, for (int alpha_right=0; alpha_right(epsilon_sgn[ie_left] * epsilon_sgn[ie_right]) * pD1()(gamma_left,gamma_left)(c_right,c_left)*D2()(alpha_right,beta_left)(a_right,b_left)*gD3g()(alpha_right,beta_left)(b_right,a_left); + result()()() -= epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1()(gamma_left,gamma_left)(c_right,c_left)*D2()(alpha_right,beta_left)(a_right,b_left)*gD3g()(alpha_right,beta_left)(b_right,a_left); }}} } //This is the \delta_{456}^{321} part @@ -225,7 +230,7 @@ void BaryonUtils::baryon_site(const mobj &D1, for (int alpha_right=0; alpha_right(epsilon_sgn[ie_left] * epsilon_sgn[ie_right]) * pD1()(gamma_left,beta_left)(c_right,b_left)*D2g()(alpha_right,beta_left)(a_right,a_left)*gD3()(alpha_right,gamma_left)(b_right,c_left); + result()()() -= epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1()(gamma_left,beta_left)(c_right,b_left)*D2g()(alpha_right,beta_left)(a_right,a_left)*gD3()(alpha_right,gamma_left)(b_right,c_left); }}} } //This is the \delta_{456}^{213} part @@ -234,7 +239,7 @@ void BaryonUtils::baryon_site(const mobj &D1, for (int alpha_right=0; alpha_right(epsilon_sgn[ie_left] * epsilon_sgn[ie_right]) * pD1g()(gamma_left,beta_left)(c_right,a_left)*D2()(alpha_right,gamma_left)(a_right,c_left)*gD3()(alpha_right,beta_left)(b_right,b_left); + result()()() -= epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1g()(gamma_left,beta_left)(c_right,a_left)*D2()(alpha_right,gamma_left)(a_right,c_left)*gD3()(alpha_right,beta_left)(b_right,b_left); }}} } } @@ -414,10 +419,10 @@ void BaryonUtils::Sigma_to_Nucleon_Q1_NonEye_site(const mobj &Du_xi, for (int tau2=0; tau2(epsilon_sgn[ie_n] * epsilon_sgn[ie_s]) * GDsGDd_ab_bb * DuGH_at_aj; + auto ee_GDGDDG_a = epsilon_sgn[ie_n] * epsilon_sgn[ie_s] * GDsGDd_ab_bb * DuGH_at_aj; for (int gamma_s=0; gamma_s(epsilon_sgn[ie_n] * epsilon_sgn[ie_s]) * GDsGDd_ab_bb * DuGH_gt_cj; + auto ee_GDGDDG_c = epsilon_sgn[ie_n] * epsilon_sgn[ie_s] * GDsGDd_ab_bb * DuGH_gt_cj; for (int gamma_n=0; gamma_n::Sigma_to_Nucleon_Q2_Eye_site(const mobj &Dq_loop, auto GDsG_at_bi = GDsG()(alpha_s,tau)(b_s,i); for (int beta_n=0; beta_n(epsilon_sgn[ie_n] * epsilon_sgn[ie_s]) * GDsG_at_bi * DqGDd_tb_ib; + auto ee_GDGDGD = epsilon_sgn[ie_n] * epsilon_sgn[ie_s] * GDsG_at_bi * DqGDd_tb_ib; for (int gamma_s=0; gamma_s::Sigma_to_Nucleon_Q2_NonEye_site(const mobj &Du_xi, auto GDsG_at_bi = GDsG()(alpha_s,tau)(b_s,i); for (int beta_n=0; beta_n(epsilon_sgn[ie_n] * epsilon_sgn[ie_s]) * GDsG_at_bi * DuGDd_ab_ab; + auto ee_GDGDGD_a = epsilon_sgn[ie_n] * epsilon_sgn[ie_s] * GDsG_at_bi * DuGDd_ab_ab; for (int gamma_s=0; gamma_s(epsilon_sgn[ie_n] * epsilon_sgn[ie_s]) * GDsG_at_bi * DuGDd_gb_cb; + auto ee_GDGDGD_c = epsilon_sgn[ie_n] * epsilon_sgn[ie_s] * GDsG_at_bi * DuGDd_gb_cb; for (int gamma_n=0; gamma_nGlobalSum(A);A/=NP; +#define AVERAGE(A) #define PRINTIT(A) AVERAGE(A); std::cout << GridLogMessage << " Stencil " << #A << " "<< A/calls<_Nprocessors; RealD NN = _grid->NodeCount(); @@ -1281,11 +1281,13 @@ public: std::cout << GridLogMessage << " Stencil SHM mem " << (membytes)/gatheralltime/1000. << " GB/s per rank"< GridPt; typedef std::unique_ptr GridRbPt; typedef std::unique_ptr RngPt; + typedef std::unique_ptr SerialRngPt; enum class Storage {object, cache, temporary}; private: struct ObjInfo @@ -114,6 +115,7 @@ public: double getVolume(void) const; // random number generator GridParallelRNG * get4dRng(void); + GridSerialRNG * getSerialRng(void); // general memory management void addObject(const std::string name, const int moduleAddress = -1); @@ -183,6 +185,7 @@ private: unsigned int nd_; // random number generator RngPt rng4d_{nullptr}; + SerialRngPt rngSerial_{nullptr}; // object store std::vector object_; std::map objectAddress_; diff --git a/Hadrons/Module.cc b/Hadrons/Module.cc index 8e332d5a..17017513 100644 --- a/Hadrons/Module.cc +++ b/Hadrons/Module.cc @@ -93,3 +93,18 @@ GridParallelRNG & ModuleBase::rng4d(void) return r; } + +GridSerialRNG & ModuleBase::rngSerial(void) +{ + auto &r = *env().getSerialRng(); + + if (makeSeedString() != seed_) + { + seed_ = makeSeedString(); + LOG(Message) << "Seeding Serial RNG " << &r << " with string '" + << seed_ << "'" << std::endl; + r.SeedUniqueString(seed_); + } + + return r; +} diff --git a/Hadrons/Module.hpp b/Hadrons/Module.hpp index 95db9511..3bea7ff2 100644 --- a/Hadrons/Module.hpp +++ b/Hadrons/Module.hpp @@ -1,7 +1,6 @@ /************************************************************************************* Grid physics library, www.github.com/paboyle/Grid - Source file: Hadrons/Module.hpp Copyright (C) 2015-2019 @@ -196,6 +195,7 @@ protected: DEFINE_VM_ALIAS; // RNG seeded from module string GridParallelRNG &rng4d(void); + GridSerialRNG &rngSerial(void); private: std::string makeSeedString(void); private: diff --git a/Hadrons/Modules.hpp b/Hadrons/Modules.hpp index ebe82910..51cb147c 100644 --- a/Hadrons/Modules.hpp +++ b/Hadrons/Modules.hpp @@ -50,8 +50,10 @@ #include #include #include +#include #include #include +#include #include #include #include @@ -65,6 +67,13 @@ #include #include #include +#include +#include +#include +#include +#include +#include +#include #include #include #include diff --git a/Hadrons/Modules/MContraction/Meson.hpp b/Hadrons/Modules/MContraction/Meson.hpp index ee923341..34b46b73 100644 --- a/Hadrons/Modules/MContraction/Meson.hpp +++ b/Hadrons/Modules/MContraction/Meson.hpp @@ -199,7 +199,7 @@ void TMeson::execute(void) Gamma gSnk(gammaList[i].first); Gamma gSrc(gammaList[i].second); - for (unsigned int t = 0; t < buf.size(); ++t) + for (unsigned int t = 0; t < nt; ++t) { result[i].corr[t] = TensorRemove(trace(mesonConnected(q1[t], q2[t], gSnk, gSrc))); } diff --git a/Hadrons/Modules/MDistil/Distil.hpp b/Hadrons/Modules/MDistil/Distil.hpp new file mode 100644 index 00000000..01a6b9da --- /dev/null +++ b/Hadrons/Modules/MDistil/Distil.hpp @@ -0,0 +1,124 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: Hadrons/Modules/MDistil/Distil.hpp + + Copyright (C) 2015-2019 + + Author: Felix Erben + Author: Michael Marshall + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (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_Distil_hpp_ +#define Hadrons_MDistil_Distil_hpp_ + +#include +#include +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE +BEGIN_MODULE_NAMESPACE(MDistil) + +/****************************************************************************** + Distillation code that is common across modules + + Documentation on how to use this code available at + + * https://aportelli.github.io/Hadrons-doc/#/mdistil * + + Notation for (stochastic) DistilParameters taken from 1104.3870: + + TI is interlaced dilution in time (corresponding to Nt = time-dimension of the lattice) + LI is interlaced dilution in laplacian-eigenvector space (corresponding to nvec) + SI is interlaced dilution in spin (corresponding to Ns, taken from Grid, usually Ns=4) + + This code automatically computes perambulators using exact distillation if + * (TI,LI,SI) = (Nt,nvec,Ns) * + In this case, nnoise=1 and Noises is set to an array of values =1 as well. + tsrc then specifies the only timeslice on which the sources are supported. + (( for stochastic distillation, the vaue of tsrc has no meaning in this code )) + + ******************************************************************************/ + +struct DistilParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(DistilParameters, + int, nvec, + int, nnoise, + int, tsrc, + int, TI, + int, LI, + int, SI ) +}; + +/****************************************************************************** + Make a lower dimensional grid in preparation for local slice operations + ******************************************************************************/ + +inline void MakeLowerDimGrid( std::unique_ptr &up, GridCartesian * gridHD ) +{ + int nd{static_cast(gridHD->_ndimension)}; + Coordinate latt_size = gridHD->_gdimensions; + latt_size[nd-1] = 1; + Coordinate simd_layout = GridDefaultSimd(nd-1, vComplex::Nsimd()); + simd_layout.push_back( 1 ); + Coordinate mpi_layout = gridHD->_processors; + mpi_layout[nd-1] = 1; + up.reset( new GridCartesian(latt_size,simd_layout,mpi_layout,*gridHD) ); +} + +/************************************************************************************* + Rotate eigenvectors into our phase convention + First component of first eigenvector is real and positive + *************************************************************************************/ + +inline void RotateEigen(std::vector & evec) +{ + ColourVector cv0; + auto grid = evec[0].Grid(); + Coordinate siteFirst(grid->Nd(),0); + peekSite(cv0, evec[0], siteFirst); + const std::complex cplx0{cv0()()(0).real(), cv0()()(0).imag()}; + if( cplx0.imag() == 0 ) + LOG(Message) << "RotateEigen() : Site 0 : " << cplx0 << " => already meets phase convention" << std::endl; + else + { + const Real cplx0_mag{ std::abs(cplx0) }; + const std::complex std_phase{std::conj(cplx0/cplx0_mag)}; + LOG(Message) << "RotateEigen() : Site 0 : |" << cplx0 << "|=" << cplx0_mag + << " => phase=" << (std::arg(std_phase) / M_PI) << " pi" << std::endl; + { + const Grid::Complex phase{std_phase.real(),std_phase.imag()}; + for( int k = 0 ; k < evec.size() ; k++ ) + evec[k] *= phase; + // Get rid of the rounding error in imaginary phase on the very first site + peekSite(cv0, evec[0], siteFirst); + cv0()()(0).imag(0); // this should be zero after the phase multiply - force it to be so + pokeSite(cv0, evec[0], siteFirst); + } + } +} + +END_MODULE_NAMESPACE +END_HADRONS_NAMESPACE +#endif 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..f38e78d1 --- /dev/null +++ b/Hadrons/Modules/MDistil/DistilPar.hpp @@ -0,0 +1,97 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: Hadrons/Modules/MDistil/DistilPar.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_DistilPar_hpp_ +#define Hadrons_MDistil_DistilPar_hpp_ + +#include + +BEGIN_HADRONS_NAMESPACE +BEGIN_MODULE_NAMESPACE(MDistil) + +/****************************************************************************** + * DistilPar * + ******************************************************************************/ + +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) +{ + return {}; +} + +template +std::vector TDistilPar::getOutput(void) +{ + return {getName()}; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TDistilPar::setup(void) +{ + envCreate(DistilParameters, getName(), 1, par() ); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TDistilPar::execute(void) +{ + // Nothing to do. setup() created and initialised the output object +} + +END_MODULE_NAMESPACE +END_HADRONS_NAMESPACE +#endif diff --git a/Hadrons/Modules/MDistil/DistilVectors.cc b/Hadrons/Modules/MDistil/DistilVectors.cc new file mode 100644 index 00000000..69b20853 --- /dev/null +++ b/Hadrons/Modules/MDistil/DistilVectors.cc @@ -0,0 +1,36 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: Hadrons/Modules/MDistil/DistilVectors.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::TDistilVectors; diff --git a/Hadrons/Modules/MDistil/DistilVectors.hpp b/Hadrons/Modules/MDistil/DistilVectors.hpp new file mode 100644 index 00000000..de871a64 --- /dev/null +++ b/Hadrons/Modules/MDistil/DistilVectors.hpp @@ -0,0 +1,243 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: Hadrons/Modules/MDistil/DistilVectors.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_DistilVectors_hpp_ +#define Hadrons_MDistil_DistilVectors_hpp_ + +#include + +BEGIN_HADRONS_NAMESPACE +BEGIN_MODULE_NAMESPACE(MDistil) + +/****************************************************************************** + * DistilVectors * + * (Create rho and/or phi vectors) * + ******************************************************************************/ + +class DistilVectorsPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(DistilVectorsPar, + std::string, noise, + std::string, perambulator, + std::string, lapevec, + std::string, rho, + std::string, phi, + std::string, DistilParams); +}; + +template +class TDistilVectors: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl,); + // constructor + TDistilVectors(const std::string name); + // destructor + virtual ~TDistilVectors(void) {}; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +protected: + std::unique_ptr grid3d; // Owned by me, so I must delete it +public: + // These variables contain parameters + std::string RhoName; + std::string PhiName; +}; + +MODULE_REGISTER_TMP(DistilVectors, TDistilVectors, MDistil); + +/****************************************************************************** + * TDistilVectors implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TDistilVectors::TDistilVectors(const std::string name) : Module(name) {} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TDistilVectors::getInput(void) +{ + return {par().noise,par().perambulator,par().lapevec,par().DistilParams}; +} + +template +std::vector TDistilVectors::getOutput(void) +{ + RhoName = par().rho; + PhiName = par().phi; + if (RhoName.empty() && PhiName.empty()) + { + HADRONS_ERROR(Argument,"No output specified"); + } + std::vector out; + if (!RhoName.empty()) + out.push_back(RhoName); + if (!PhiName.empty()) + out.push_back(PhiName); + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TDistilVectors::setup(void) +{ + // We expect the perambulator to have been created with these indices + auto &perambulator = envGet(PerambTensor, par().perambulator); + if (!perambulator.ValidateIndexNames()) + { + HADRONS_ERROR(Range,"Perambulator index names bad"); + } + + const DistilParameters &dp{envGet(DistilParameters, par().DistilParams)}; + const int Nt{env().getDim(Tdir)}; + const bool full_tdil{ dp.TI == Nt }; + const int Nt_inv{ full_tdil ? 1 : dp.TI }; + + if (!RhoName.empty()) + envCreate(std::vector, RhoName, 1, dp.nnoise*dp.LI*dp.SI*Nt_inv, envGetGrid(FermionField)); + if (!PhiName.empty()) + envCreate(std::vector, PhiName, 1, dp.nnoise*dp.LI*dp.SI*Nt_inv, envGetGrid(FermionField)); + + Coordinate latt_size = GridDefaultLatt(); + Coordinate mpi_layout = GridDefaultMpi(); + Coordinate simd_layout_3 = GridDefaultSimd(Nd-1, vComplex::Nsimd()); + latt_size[Nd-1] = 1; + simd_layout_3.push_back( 1 ); + mpi_layout[Nd-1] = 1; + GridCartesian * const grid4d{env().getGrid()}; + MakeLowerDimGrid(grid3d, grid4d); + + envTmp(LatticeSpinColourVector, "source4d",1,LatticeSpinColourVector(grid4d)); + envTmp(LatticeSpinColourVector, "source3d",1,LatticeSpinColourVector(grid3d.get())); + envTmp(LatticeColourVector, "source3d_nospin",1,LatticeColourVector(grid3d.get())); + envTmp(LatticeSpinColourVector, "sink3d",1,LatticeSpinColourVector(grid3d.get())); + envTmp(LatticeColourVector, "evec3d",1,LatticeColourVector(grid3d.get())); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TDistilVectors::execute(void) +{ + auto &noise = envGet(NoiseTensor, par().noise); + auto &perambulator = envGet(PerambTensor, par().perambulator); + auto &epack = envGet(Grid::Hadrons::EigenPack, par().lapevec); + const DistilParameters &dp{envGet(DistilParameters, par().DistilParams)}; + + envGetTmp(LatticeSpinColourVector, source4d); + envGetTmp(LatticeSpinColourVector, source3d); + envGetTmp(LatticeColourVector, source3d_nospin); + envGetTmp(LatticeSpinColourVector, sink3d); + envGetTmp(LatticeColourVector, evec3d); + + GridCartesian * const grid4d{env().getGrid()}; + const int Ntlocal{ grid4d->LocalDimensions()[3] }; + const int Ntfirst{ grid4d->LocalStarts()[3] }; + + const int Nt{env().getDim(Tdir)}; + const bool full_tdil{ dp.TI == Nt }; + const int Nt_inv{ full_tdil ? 1 : dp.TI }; + + int vecindex; + if (!RhoName.empty()) + { + auto &rho = envGet(std::vector, RhoName); + for (int inoise = 0; inoise < dp.nnoise; inoise++) + { + for (int dk = 0; dk < dp.LI; dk++) + { + for (int dt = 0; dt < Nt_inv; dt++) + { + for (int ds = 0; ds < dp.SI; ds++) + { + vecindex = inoise + dp.nnoise * (dk + dp.LI * (ds + dp.SI * dt)); + rho[vecindex] = 0; + for (int it = dt; it < Nt; it += dp.TI) + { + const int t_inv{full_tdil ? dp.tsrc : it}; + if (t_inv >= Ntfirst && t_inv < Ntfirst + Ntlocal) + { + for (int ik = dk; ik < dp.nvec; ik += dp.LI) + { + for (int is = ds; is < Ns; is += dp.SI) + { + ExtractSliceLocal(evec3d,epack.evec[ik],0,t_inv-Ntfirst,Tdir); + 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; + } + } + } + } + } + } + } + } + } + if (!PhiName.empty()) + { + auto &phi = envGet(std::vector, PhiName); + for (int inoise = 0; inoise < dp.nnoise; inoise++) + { + for (int dk = 0; dk < dp.LI; dk++) + { + for (int dt = 0; dt < Nt_inv; dt++) + { + for (int ds = 0; ds < dp.SI; ds++) + { + vecindex = inoise + dp.nnoise * (dk + dp.LI * (ds + dp.SI * dt)); + phi[vecindex] = 0; + for (int t = Ntfirst; t < Ntfirst + Ntlocal; t++) + { + sink3d=0; + for (int ivec = 0; ivec < dp.nvec; ivec++) + { + ExtractSliceLocal(evec3d,epack.evec[ivec],0,t-Ntfirst,Tdir); + sink3d += evec3d * perambulator.tensor(t, ivec, dk, inoise,dt,ds); + } + InsertSliceLocal(sink3d,phi[vecindex],0,t-Ntfirst,Tdir); + } + } + } + } + } + } +} + +END_MODULE_NAMESPACE +END_HADRONS_NAMESPACE +#endif // Hadrons_MDistil_DistilVectors_hpp_ diff --git a/Hadrons/Modules/MDistil/LapEvec.cc b/Hadrons/Modules/MDistil/LapEvec.cc new file mode 100644 index 00000000..be6838fa --- /dev/null +++ b/Hadrons/Modules/MDistil/LapEvec.cc @@ -0,0 +1,36 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: Hadrons/Modules/MDistil/LapEvec.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::TLapEvec; diff --git a/Hadrons/Modules/MDistil/LapEvec.hpp b/Hadrons/Modules/MDistil/LapEvec.hpp new file mode 100644 index 00000000..9b5197c2 --- /dev/null +++ b/Hadrons/Modules/MDistil/LapEvec.hpp @@ -0,0 +1,335 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: Hadrons/Modules/MDistil/LapEvec.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_LapEvec_hpp_ +#define Hadrons_MDistil_LapEvec_hpp_ + +#include + +BEGIN_HADRONS_NAMESPACE +BEGIN_MODULE_NAMESPACE(MDistil) + +/****************************************************************************** + + Laplacian eigenvectors - parameters + + Computes the eigenvectors of the 3D-Laplacian, built from stout-smeared + gauge links with the specified number of steps and smearing parameter rho. + The smearing is only applied to the spatial components of the gauge field, + i.e. rho_{4i} = rho_{i4} = rho_{44} = 0. + + Chebyshev-preconditioning is needed for convergence of the nvec lowest + eigenvectors. + + ******************************************************************************/ + +struct StoutParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(StoutParameters, + int, steps, + double, rho) + StoutParameters() = default; + template StoutParameters(Reader& Reader){read(Reader,"StoutSmearing",*this);} +}; + +struct ChebyshevParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(ChebyshevParameters, + int, PolyOrder, + double, alpha, + double, beta) + ChebyshevParameters() = default; + template ChebyshevParameters(Reader& Reader){read(Reader,"Chebyshev",*this);} +}; + +struct LanczosParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(LanczosParameters, + int, Nvec, + int, Nk, + int, Np, + int, MaxIt, + double, resid, + int, IRLLog) + LanczosParameters() = default; + template LanczosParameters(Reader& Reader){read(Reader,"Lanczos",*this);} +}; + +// These are the actual parameters passed to the module during construction + +struct LapEvecPar: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(LapEvecPar + ,std::string, gauge + ,StoutParameters, Stout + ,ChebyshevParameters, Cheby + ,LanczosParameters, Lanczos + ,std::string, FileName) +}; + +/****************************************************************************** + + Laplacian eigenvectors - Module (class) definition + + ******************************************************************************/ + +template +class TLapEvec: public Module +{ +public: + GAUGE_TYPE_ALIASES(GImpl,); + // constructor + TLapEvec(const std::string name); + // destructor + virtual ~TLapEvec(void) {}; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +protected: + std::unique_ptr gridLD; // Owned by me, so I must delete it +}; + +MODULE_REGISTER_TMP(LapEvec, TLapEvec, MDistil); + +/****************************************************************************** + TLapEvec implementation + ******************************************************************************/ + +// constructor ///////////////////////////////////////////////////////////////// +template +TLapEvec::TLapEvec(const std::string name) : Module(name) {} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TLapEvec::getInput(void) +{ + return std::vector{par().gauge}; +} + +template +std::vector TLapEvec::getOutput(void) +{ + return {getName()}; // This is the higher dimensional eigenpack +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TLapEvec::setup(void) +{ + GridCartesian * gridHD = env().getGrid(); + MakeLowerDimGrid(gridLD,gridHD); + const int Ntlocal{gridHD->LocalDimensions()[Tdir]}; + // Temporaries + envTmpLat(GaugeField, "Umu_stout"); + envTmpLat(GaugeField, "Umu_smear"); + envTmp(LatticeGaugeField, "UmuNoTime",1,LatticeGaugeField(gridLD.get())); + envTmp(LatticeColourVector, "src",1,LatticeColourVector(gridLD.get())); + envTmp(std::vector, "eig",1,std::vector(Ntlocal)); + // Output objects + envCreate(LapEvecs, getName(), 1, par().Lanczos.Nvec, gridHD); +} + +/************************************************************************************* + + -Grad^2 (Peardon, 2009, pg 2, equation 3, https://arxiv.org/abs/0905.2160) + Field Type of field the operator will be applied to + GaugeField Gauge field the operator will smear using + + *************************************************************************************/ + +template +class Laplacian3D : public LinearOperatorBase, public LinearFunction { + typedef typename GaugeField::vector_type vCoeff_t; +public: + int nd; // number of spatial dimensions + std::vector > > U; + // Construct this operator given a gauge field and the number of dimensions it should act on + Laplacian3D( GaugeField& gf, int dimSpatial = Tdir ) : nd{dimSpatial} + { + if (dimSpatial<1) + { + HADRONS_ERROR(Range,"Must be at least one spatial dimension"); + } + for (int mu = 0 ; mu < nd ; mu++) + U.push_back(PeekIndex(gf,mu)); + } + + // Apply this operator to "in", return result in "out" + void operator()(const Field& in, Field& out) { + if (nd > in.Grid()->Nd()) + { + HADRONS_ERROR(Range,"nd too large"); + } + conformable( in, out ); + out = ( ( Real ) ( 2 * nd ) ) * in; + Field tmp_(in.Grid()); + typedef typename GaugeField::vector_type vCoeff_t; + for (int mu = 0 ; mu < nd ; mu++) + { + out -= U[mu] * Cshift( in, mu, 1); + tmp_ = adj( U[mu] ) * in; + out -= Cshift(tmp_,mu,-1); + } + } + + void OpDiag (const Field &in, Field &out) { HADRONS_ERROR(Definition, "OpDiag() undefined"); }; + void OpDir (const Field &in, Field &out,int dir,int disp) { HADRONS_ERROR(Definition, "OpDir() undefined"); }; + void Op (const Field &in, Field &out) { HADRONS_ERROR(Definition, "Op() undefined"); }; + void AdjOp (const Field &in, Field &out) { HADRONS_ERROR(Definition, "AdjOp() undefined"); }; + void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2) { HADRONS_ERROR(Definition, "HermOpAndNorm() undefined"); }; + void HermOp(const Field &in, Field &out) { operator()(in,out); }; +}; + +template +class Laplacian3DHerm : public LinearFunction { +public: + OperatorFunction & poly_; + LinearOperatorBase &Linop_; + Laplacian3DHerm(OperatorFunction & poly,LinearOperatorBase& linop) + : poly_{poly}, Linop_{linop} {} + void operator()(const Field& in, Field& out) + { + poly_(Linop_,in,out); + } +}; + +/****************************************************************************** + Calculate low-mode eigenvalues of the Laplacian + ******************************************************************************/ + +// execution /////////////////////////////////////////////////////////////////// +template +void TLapEvec::execute(void) +{ + const ChebyshevParameters &ChebPar{par().Cheby}; + const LanczosParameters &LPar{par().Lanczos}; + + // Disable IRL logging if requested + LOG(Message) << "IRLLog=" << LPar.IRLLog << std::endl; + const int PreviousIRLLogState{GridLogIRL.isActive()}; + GridLogIRL.Active( LPar.IRLLog == 0 ? 0 : 1 ); + + // Stout smearing + envGetTmp(GaugeField, Umu_smear); + Umu_smear = envGet(GaugeField, par().gauge); // The smeared field starts off as the Gauge field + LOG(Message) << "Initial plaquette: " << WilsonLoops::avgPlaquette(Umu_smear) << std::endl; + const StoutParameters &Stout{par().Stout}; + if( Stout.steps ) + { + envGetTmp(GaugeField, Umu_stout); + 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; + } + LOG(Message) << "Smeared plaquette: " << WilsonLoops::avgPlaquette(Umu_smear) << std::endl; + } + + //////////////////////////////////////////////////////////////////////// + // Invert nabla operator separately on each time-slice + //////////////////////////////////////////////////////////////////////// + + auto & eig4d = envGet(LapEvecs, getName() ); + envGetTmp(std::vector, eig); // Eigenpack for each timeslice + envGetTmp(LatticeGaugeField, UmuNoTime); // Gauge field without time dimension + envGetTmp(LatticeColourVector, src); + GridCartesian * gridHD = env().getGrid(); + const int Ntlocal{gridHD->LocalDimensions()[Tdir]}; + const int Ntfirst{gridHD->LocalStarts()[Tdir]}; + uint32_t ConvergenceErrors{0}; + for (int t = 0; t < Ntlocal; t++ ) + { + LOG(Message) << "------------------------------------------------------------" << std::endl; + LOG(Message) << " Compute eigenpack, local timeslice = " << t << " / " << Ntlocal << std::endl; + LOG(Message) << "------------------------------------------------------------" << std::endl; + eig[t].resize(LPar.Nk+LPar.Np,gridLD.get()); + + // Construct smearing operator + ExtractSliceLocal(UmuNoTime,Umu_smear,0,t,Tdir); // switch to 3d/4d objects + Laplacian3D Nabla(UmuNoTime); + LOG(Message) << "Chebyshev preconditioning to order " << ChebPar.PolyOrder + << " with parameters (alpha,beta) = (" << ChebPar.alpha << "," << ChebPar.beta << ")" << std::endl; + Chebyshev Cheb(ChebPar.alpha,ChebPar.beta,ChebPar.PolyOrder); + + // Construct source vector according to Test_dwf_compressed_lanczos.cc + src = 11.0; // NB: This is a dummy parameter and just needs to be non-zero + RealD nn = norm2(src); + nn = Grid::sqrt(nn); + src = src * (1.0/nn); + + Laplacian3DHerm NablaCheby(Cheb,Nabla); + ImplicitlyRestartedLanczos + IRL(NablaCheby,Nabla,LPar.Nvec,LPar.Nk,LPar.Nk+LPar.Np,LPar.resid,LPar.MaxIt); + int Nconv = 0; + IRL.calc(eig[t].eval,eig[t].evec,src,Nconv); + if (Nconv < LPar.Nvec) + { + // NB: Can't assert here since we are processing local slices - i.e. not all nodes would assert + ConvergenceErrors = 1; + LOG(Error) << "MDistil::LapEvec : Not enough eigenvectors converged. If this occurs in practice, we should modify the eigensolver to iterate once more to ensure the second convergence test does not take us below the requested number of eigenvectors" << std::endl; + } + if( Nconv != LPar.Nvec ) + eig[t].resize(LPar.Nvec, gridLD.get()); + RotateEigen( eig[t].evec ); // Rotate the eigenvectors into our phase convention + + for (int i=0;iGlobalSum(ConvergenceErrors); + if(ConvergenceErrors!=0) + { + HADRONS_ERROR(Program,"The eingensolver failed to find enough eigenvectors on at least one node"); + } + // Now write out the 4d eigenvectors + std::string sEigenPackName(par().FileName); + if( !sEigenPackName.empty() ) + { + eig4d.record.solverXml = parString(); + ModuleBase * b{vm().getModule(par().gauge)}; + std::string sOperatorXml{ "" }; + sOperatorXml.append( b->getRegisteredName() ); + sOperatorXml.append( "" ); + sOperatorXml.append( b->parString() ); + sOperatorXml.append( "" ); + eig4d.record.operatorXml = sOperatorXml; + sEigenPackName.append("."); + sEigenPackName.append(std::to_string(vm().getTrajectory())); + eig4d.write(sEigenPackName,false); + } +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MDistil_LapEvec_hpp_ diff --git a/Hadrons/Modules/MDistil/Noises.cc b/Hadrons/Modules/MDistil/Noises.cc new file mode 100644 index 00000000..16576b0d --- /dev/null +++ b/Hadrons/Modules/MDistil/Noises.cc @@ -0,0 +1,36 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: Hadrons/Modules/MDistil/Noises.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::TNoises; diff --git a/Hadrons/Modules/MDistil/Noises.hpp b/Hadrons/Modules/MDistil/Noises.hpp new file mode 100644 index 00000000..bd1a43ff --- /dev/null +++ b/Hadrons/Modules/MDistil/Noises.hpp @@ -0,0 +1,146 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: Hadrons/Modules/MDistil/Noises.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_Noises_hpp_ +#define Hadrons_MDistil_Noises_hpp_ + +#include + +BEGIN_HADRONS_NAMESPACE +BEGIN_MODULE_NAMESPACE(MDistil) + +/****************************************************************************** + * Noises * + ******************************************************************************/ + +class NoisesPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(NoisesPar, + std::string, DistilParams, + std::string, NoiseFileName) +}; + +template +class TNoises: public Module +{ +public: + // constructor + TNoises(const std::string name); + // destructor + virtual ~TNoises(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(Noises, TNoises, MDistil); + +/****************************************************************************** + * TNoises implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TNoises::TNoises(const std::string name) : Module(name) {} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TNoises::getInput(void) +{ + return {par().DistilParams}; +} + +template +std::vector TNoises::getOutput(void) +{ + return {getName()}; +} + +// setup /////////////////////////////////////////////////////////////////////// + +template +void TNoises::setup(void) +{ + const DistilParameters &dp{envGet(DistilParameters, par().DistilParams)}; + const int Nt{env().getDim(Tdir)}; + std::cout << dp.nnoise << dp.nvec << Nt << Ns << std::endl; + envCreate(NoiseTensor, getName(), 1, dp.nnoise, Nt, dp.nvec, Ns); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TNoises::execute(void) +{ + const DistilParameters &dp{envGet(DistilParameters, par().DistilParams)}; + const int Nt{env().getDim(Tdir)}; + const bool full_tdil{ dp.TI == Nt }; + const bool exact_distillation{ full_tdil && dp.LI == dp.nvec }; + + // We use our own seeds so we can specify different noises per quark + Real rn; + auto &noise = envGet(NoiseTensor, getName()); + for (int inoise = 0; inoise < dp.nnoise; inoise++) + { + for (int t = 0; t < Nt; t++) + { + for (int ivec = 0; ivec < dp.nvec; ivec++) + { + for (int is = 0; is < Ns; is++) + { + if (exact_distillation) + { + noise.tensor(inoise, t, ivec, is) = 1.; + } + else + { + random(rngSerial(),rn); + // We could use a greater number of complex roots of unity + // ... but this seems to work well + noise.tensor(inoise, t, ivec, is) = (rn > 0.5) ? -1 : 1; + } + } + } + } + } + if (env().getGrid()->IsBoss()) + { + std::string sName {par().NoiseFileName}; + sName.append("."); + sName.append(std::to_string(vm().getTrajectory())); + noise.write(sName.c_str()); + } +} + +END_MODULE_NAMESPACE +END_HADRONS_NAMESPACE +#endif diff --git a/Hadrons/Modules/MDistil/PerambFromSolve.cc b/Hadrons/Modules/MDistil/PerambFromSolve.cc new file mode 100644 index 00000000..69cd0d35 --- /dev/null +++ b/Hadrons/Modules/MDistil/PerambFromSolve.cc @@ -0,0 +1,36 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: Hadrons/Modules/MDistil/PerambFromSolve.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::TPerambFromSolve; diff --git a/Hadrons/Modules/MDistil/PerambFromSolve.hpp b/Hadrons/Modules/MDistil/PerambFromSolve.hpp new file mode 100644 index 00000000..a0a9b61c --- /dev/null +++ b/Hadrons/Modules/MDistil/PerambFromSolve.hpp @@ -0,0 +1,183 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: Hadrons/Modules/MDistil/PerambFromSolve.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_PerambFromSolve_hpp_ +#define Hadrons_MDistil_PerambFromSolve_hpp_ + +#include + +BEGIN_HADRONS_NAMESPACE +BEGIN_MODULE_NAMESPACE(MDistil) + +/****************************************************************************** + * PerambFromSolve + + This module computes a perambulator from an already completed solve. + Optionally, the number of eigenvectors used in the perambulator and the + parameter LI can be chosen to be lower than the ones in the solve, allowing + for a study of the signal with different values of nvec. + + LI_reduced : value of LI actually used in the computation + nvec_reduced: value of nvec actually used in the computation + LI : value of LI used to compute the 'solve' + nvec : value of nvec used to compute the 'solve' + + ******************************************************************************/ + +class PerambFromSolvePar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(PerambFromSolvePar, + std::string, lapevec, + std::string, PerambFileName, + std::string, solve, + int, nvec_reduced, + int, LI_reduced, + std::string, DistilParams); +}; + +template +class TPerambFromSolve: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl,); + // constructor + TPerambFromSolve(const std::string name); + // destructor + virtual ~TPerambFromSolve(void) {}; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +protected: + std::unique_ptr grid3d; // Owned by me, so I must delete it +}; + +MODULE_REGISTER_TMP(PerambFromSolve, TPerambFromSolve, MDistil); + +/****************************************************************************** + * TPerambFromSolve implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TPerambFromSolve::TPerambFromSolve(const std::string name) : Module(name){} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TPerambFromSolve::getInput(void) +{ + return {par().solve, par().lapevec, par().DistilParams}; +} + +template +std::vector TPerambFromSolve::getOutput(void) +{ + return std::vector{getName()}; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TPerambFromSolve::setup(void) +{ + const DistilParameters & dp{envGet(MDistil::DistilParameters, par().DistilParams)}; + const int Nt{env().getDim(Tdir)}; + const bool full_tdil{ dp.TI == Nt }; + const int Nt_inv{ full_tdil ? 1 : dp.TI }; + MakeLowerDimGrid( grid3d, env().getGrid() ); + const int nvec_reduced{par().nvec_reduced}; + const int LI_reduced{ par().LI_reduced}; + envCreate(PerambTensor, getName(), 1, Nt,nvec_reduced,LI_reduced,dp.nnoise,Nt_inv,dp.SI); + envCreate(NoiseTensor, getName() + "_noise", 1, dp.nnoise, Nt, dp.nvec, Ns ); + envTmp(LatticeColourVector, "result3d_nospin",1,LatticeColourVector(grid3d.get())); + envTmp(LatticeColourVector, "evec3d",1,LatticeColourVector(grid3d.get())); + envTmpLat(LatticeColourVector, "result4d_nospin"); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TPerambFromSolve::execute(void) +{ + GridCartesian * grid4d = env().getGrid(); + const int Ntlocal{grid4d->LocalDimensions()[3]}; + const int Ntfirst{grid4d->LocalStarts()[3]}; + const DistilParameters &dp{envGet(DistilParameters, par().DistilParams)}; + const int Nt{env().getDim(Tdir)}; + const bool full_tdil{ dp.TI == Nt }; + const int Nt_inv{ full_tdil ? 1 : dp.TI }; + const int nvec_reduced{par().nvec_reduced}; + const int LI_reduced{ par().LI_reduced}; + auto &perambulator = envGet(PerambTensor, getName()); + auto &solve = envGet(std::vector, par().solve); + auto &epack = envGet(Grid::Hadrons::EigenPack, par().lapevec); + + envGetTmp(LatticeColourVector, result4d_nospin); + envGetTmp(LatticeColourVector, result3d_nospin); + envGetTmp(LatticeColourVector, evec3d); + + for (int inoise = 0; inoise < dp.nnoise; inoise++) + { + for (int dk = 0; dk < LI_reduced; dk++) + { + for (int dt = 0; dt < Nt_inv; dt++) + { + for (int ds = 0; ds < dp.SI; ds++) + { + for (int is = 0; is < Ns; is++) + { + result4d_nospin = peekSpin(solve[inoise+dp.nnoise*(dk+dp.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, 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; + } + } + } + } + } + } + } + if(grid4d->IsBoss()) + { + std::string sPerambName{par().PerambFileName}; + sPerambName.append( "." ); + sPerambName.append( std::to_string(vm().getTrajectory())); + perambulator.write(sPerambName.c_str()); + } +} + +END_MODULE_NAMESPACE +END_HADRONS_NAMESPACE +#endif // Hadrons_MDistil_PerambFromSolve_hpp_ diff --git a/Hadrons/Modules/MDistil/Perambulator.cc b/Hadrons/Modules/MDistil/Perambulator.cc new file mode 100644 index 00000000..e710698e --- /dev/null +++ b/Hadrons/Modules/MDistil/Perambulator.cc @@ -0,0 +1,57 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: Hadrons/Modules/MDistil/Perambulator.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::TPerambulator; + +BEGIN_HADRONS_NAMESPACE + +// Global constants for distillation + +#ifdef HAVE_HDF5 +extern const std::string NamedTensorFileExtension{".h5"}; +#else +extern const std::string NamedTensorFileExtension{".dat"}; +#endif + +BEGIN_MODULE_NAMESPACE(MDistil) + +const std::string NoiseTensor::Name__{"Noises"}; +const std::array NoiseTensor::DefaultIndexNames__{"nNoise", "nT", "nVec", "nS"}; + +const std::string PerambTensor::Name__{"Perambulator"}; +const std::array PerambTensor::DefaultIndexNames__{"nT", "nVec", "LI", "nNoise", "nT_inv", "SI"}; + +END_MODULE_NAMESPACE +END_HADRONS_NAMESPACE diff --git a/Hadrons/Modules/MDistil/Perambulator.hpp b/Hadrons/Modules/MDistil/Perambulator.hpp new file mode 100644 index 00000000..fb9d16dc --- /dev/null +++ b/Hadrons/Modules/MDistil/Perambulator.hpp @@ -0,0 +1,263 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: Hadrons/Modules/MDistil/Perambulator.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_Perambulator_hpp_ +#define Hadrons_MDistil_Perambulator_hpp_ + +#include + +BEGIN_HADRONS_NAMESPACE +BEGIN_MODULE_NAMESPACE(MDistil) + +/****************************************************************************** + * Perambulator * + ******************************************************************************/ + +class PerambulatorPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(PerambulatorPar, + std::string, lapevec, + std::string, solver, + std::string, noise, + std::string, PerambFileName, + std::string, UnsmearedSinkFileName, + std::string, DistilParams); +}; + +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); +protected: + std::unique_ptr grid3d; // Owned by me, so I must delete it + unsigned int Ls_; +}; + +MODULE_REGISTER_TMP(Perambulator, TPerambulator, MDistil); + +/****************************************************************************** + * TPerambulator implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TPerambulator::TPerambulator(const std::string name) : Module(name) {} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TPerambulator::getInput(void) +{ + return {par().lapevec, par().solver, par().noise, par().DistilParams}; +} + +template +std::vector TPerambulator::getOutput(void) +{ + return {getName(), getName() + "_unsmeared_sink"}; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TPerambulator::setup(void) +{ + MakeLowerDimGrid(grid3d, env().getGrid()); + const DistilParameters &dp = envGet(DistilParameters, par().DistilParams); + const int Nt{env().getDim(Tdir)}; + const bool full_tdil{ dp.TI == Nt }; + const int Nt_inv{ full_tdil ? 1 : dp.TI }; + + envCreate(PerambTensor, getName(), 1, Nt, dp.nvec, dp.LI, dp.nnoise, Nt_inv, dp.SI); + envCreate(std::vector, getName() + "_unsmeared_sink", 1, + dp.nnoise*dp.LI*Ns*Nt_inv, envGetGrid(FermionField)); + + envTmpLat(LatticeSpinColourVector, "dist_source"); + envTmpLat(LatticeSpinColourVector, "source4d"); + envTmp(LatticeSpinColourVector, "source3d",1,LatticeSpinColourVector(grid3d.get())); + envTmp(LatticeColourVector, "source3d_nospin",1,LatticeColourVector(grid3d.get())); + envTmpLat(LatticeSpinColourVector, "result4d"); + envTmpLat(LatticeColourVector, "result4d_nospin"); + envTmp(LatticeColourVector, "result3d_nospin",1,LatticeColourVector(grid3d.get())); + envTmp(LatticeColourVector, "evec3d",1,LatticeColourVector(grid3d.get())); + + Ls_ = env().getObjectLs(par().solver); + envTmpLat(FermionField, "v4dtmp"); + envTmpLat(FermionField, "v5dtmp", Ls_); + envTmpLat(FermionField, "v5dtmp_sol", Ls_); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TPerambulator::execute(void) +{ + const DistilParameters &dp{ envGet(DistilParameters, par().DistilParams) }; + const int Nt{env().getDim(Tdir)}; + const bool full_tdil{ dp.TI == Nt }; + const int Nt_inv{ full_tdil ? 1 : dp.TI }; + + auto &solver=envGet(Solver, par().solver); + auto &mat = solver.getFMat(); + envGetTmp(FermionField, v4dtmp); + envGetTmp(FermionField, v5dtmp); + envGetTmp(FermionField, v5dtmp_sol); + auto &noise = envGet(NoiseTensor, par().noise); + auto &perambulator = envGet(PerambTensor, getName()); + auto &epack = envGet(LapEvecs, par().lapevec); + auto &unsmeared_sink = envGet(std::vector, getName() + "_unsmeared_sink"); + envGetTmp(LatticeSpinColourVector, dist_source); + envGetTmp(LatticeSpinColourVector, source4d); + envGetTmp(LatticeSpinColourVector, source3d); + envGetTmp(LatticeColourVector, source3d_nospin); + envGetTmp(LatticeSpinColourVector, result4d); + envGetTmp(LatticeColourVector, result4d_nospin); + envGetTmp(LatticeColourVector, result3d_nospin); + envGetTmp(LatticeColourVector, evec3d); + GridCartesian * const grid4d{ env().getGrid() }; // Owned by environment (so I won't delete it) + const int Ntlocal{grid4d->LocalDimensions()[3]}; + const int Ntfirst{grid4d->LocalStarts()[3]}; + const std::string UnsmearedSinkFileName{ par().UnsmearedSinkFileName }; + + for (int inoise = 0; inoise < dp.nnoise; inoise++) + { + for (int dk = 0; dk < dp.LI; dk++) + { + for (int dt = 0; dt < Nt_inv; dt++) + { + for (int ds = 0; ds < dp.SI; ds++) + { + 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; + evec3d = 0; + for (int it = dt; it < Nt; it += dp.TI) + { + const int t_inv{full_tdil ? dp.tsrc : it}; + if( t_inv >= Ntfirst && t_inv < Ntfirst + Ntlocal ) + { + for (int ik = dk; ik < dp.nvec; ik += dp.LI) + { + for (int is = ds; is < Ns; is += dp.SI) + { + ExtractSliceLocal(evec3d,epack.evec[ik],0,t_inv-Ntfirst,Tdir); + 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; + } + } + } + } + result4d=0; + v4dtmp = dist_source; + if (Ls_ == 1) + solver(result4d, v4dtmp); + else + { + mat.ImportPhysicalFermionSource(v4dtmp, v5dtmp); + solver(v5dtmp_sol, v5dtmp); + mat.ExportPhysicalFermionSolution(v5dtmp_sol, v4dtmp); + result4d = v4dtmp; + } + if (!UnsmearedSinkFileName.empty()) + unsmeared_sink[inoise+dp.nnoise*(dk+dp.LI*(dt+Nt_inv*ds))] = result4d; + for (int is = 0; is < Ns; is++) + { + result4d_nospin = peekSpin(result4d,is); + for (int t = Ntfirst; t < Ntfirst + Ntlocal; t++) + { + ExtractSliceLocal(result3d_nospin,result4d_nospin,0,t-Ntfirst,Tdir); + for (int ivec = 0; ivec < dp.nvec; ivec++) + { + ExtractSliceLocal(evec3d,epack.evec[ivec],0,t-Ntfirst,Tdir); + pokeSpin(perambulator.tensor(t, ivec, dk, inoise,dt,ds),static_cast(innerProduct(evec3d, result3d_nospin)),is); + } + } + } + } + } + } + } + // Now share my timeslice data with other members of the grid + const int NumSlices{grid4d->_processors[Tdir] / grid3d->_processors[Tdir]}; + if (NumSlices > 1) + { + LOG(Debug) << "Sharing perambulator data with other nodes" << std::endl; + const int MySlice {grid4d->_processor_coor[Tdir]}; + const int SliceCount {static_cast(perambulator.tensor.size()/NumSlices)}; + PerambTensor::Scalar * const MyData {perambulator.tensor.data()+MySlice*SliceCount}; + Coordinate coor(Nd); + for (int i = 0 ; i < Tdir ; i++) coor[i] = grid4d->_processor_coor[i]; + std::vector reqs(0); + for (int i = 1; i < NumSlices ; i++) + { + coor[Tdir] = (MySlice+i)%NumSlices; + const int SendRank { grid4d->RankFromProcessorCoor(coor) }; + const int RecvSlice { ( MySlice - i + NumSlices ) % NumSlices }; + coor[Tdir] = RecvSlice; + const auto RecvRank = grid4d->RankFromProcessorCoor(coor); + grid4d->SendToRecvFromBegin(reqs,MyData,SendRank, perambulator.tensor.data() + + RecvSlice*SliceCount,RecvRank,SliceCount*sizeof(PerambTensor::Scalar)); + } + grid4d->SendToRecvFromComplete(reqs); + } + + // Save the perambulator to disk from the boss node + if (grid4d->IsBoss()) + { + std::string sPerambName {par().PerambFileName}; + sPerambName.append("."); + sPerambName.append(std::to_string(vm().getTrajectory())); + perambulator.write(sPerambName.c_str()); + } + + //Save the unsmeared sinks if filename specified + if (!UnsmearedSinkFileName.empty()) + { + LOG(Message) << "Writing unsmeared sink to " << UnsmearedSinkFileName << std::endl; + A2AVectorsIo::write(UnsmearedSinkFileName, unsmeared_sink, false, vm().getTrajectory()); + } +} + +END_MODULE_NAMESPACE +END_HADRONS_NAMESPACE + +#endif // Hadrons_MDistil_Perambulator_hpp_ diff --git a/Hadrons/Modules/MIO/LoadDistilNoise.cc b/Hadrons/Modules/MIO/LoadDistilNoise.cc new file mode 100644 index 00000000..da532830 --- /dev/null +++ b/Hadrons/Modules/MIO/LoadDistilNoise.cc @@ -0,0 +1,7 @@ +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MIO; + +template class Grid::Hadrons::MIO::TLoadDistilNoise; diff --git a/Hadrons/Modules/MIO/LoadDistilNoise.hpp b/Hadrons/Modules/MIO/LoadDistilNoise.hpp new file mode 100644 index 00000000..d23a0d02 --- /dev/null +++ b/Hadrons/Modules/MIO/LoadDistilNoise.hpp @@ -0,0 +1,111 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: Hadrons/Modules/MDistil/LoadDistilNoise.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_MIO_LoadDistilNoise_hpp_ +#define Hadrons_MIO_LoadDistilNoise_hpp_ + +#include + +BEGIN_HADRONS_NAMESPACE +BEGIN_MODULE_NAMESPACE(MIO) + +/****************************************************************************** + * LoadDistilNoise * + ******************************************************************************/ + +class LoadDistilNoisePar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(LoadDistilNoisePar, + std::string, NoiseFileName, + std::string, DistilParams); +}; + +template +class TLoadDistilNoise: public Module +{ +public: + // constructor + TLoadDistilNoise(const std::string name); + // destructor + virtual ~TLoadDistilNoise(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(LoadDistilNoise, TLoadDistilNoise, MIO); + +/****************************************************************************** + * TLoadDistilNoise implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TLoadDistilNoise::TLoadDistilNoise(const std::string name) : Module(name) {} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TLoadDistilNoise::getInput(void) +{ + return {par().DistilParams}; +} + +template +std::vector TLoadDistilNoise::getOutput(void) +{ + return {getName()}; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TLoadDistilNoise::setup(void) +{ + const MDistil::DistilParameters &dp{envGet(MDistil::DistilParameters, par().DistilParams)}; + const int Nt{env().getDim(Tdir)}; + envCreate(MDistil::NoiseTensor, getName(), 1, dp.nnoise, Nt, dp.nvec, Ns); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TLoadDistilNoise::execute(void) +{ + auto &noises = envGet(MDistil::NoiseTensor, getName()); + std::string sNoiseName{ par().NoiseFileName }; + sNoiseName.append( 1, '.' ); + sNoiseName.append( std::to_string( vm().getTrajectory() ) ); + noises.read(sNoiseName.c_str()); +} + +END_MODULE_NAMESPACE +END_HADRONS_NAMESPACE +#endif diff --git a/Hadrons/Modules/MIO/LoadPerambulator.cc b/Hadrons/Modules/MIO/LoadPerambulator.cc new file mode 100644 index 00000000..870d8da4 --- /dev/null +++ b/Hadrons/Modules/MIO/LoadPerambulator.cc @@ -0,0 +1,36 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: Hadrons/Modules/MDistil/LoadPerambulator.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 MIO; + +template class Grid::Hadrons::MIO::TLoadPerambulator; diff --git a/Hadrons/Modules/MIO/LoadPerambulator.hpp b/Hadrons/Modules/MIO/LoadPerambulator.hpp new file mode 100644 index 00000000..cd5a4cc8 --- /dev/null +++ b/Hadrons/Modules/MIO/LoadPerambulator.hpp @@ -0,0 +1,113 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: Hadrons/Modules/MDistil/LoadPerambulator.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_MIO_LoadPerambulator_hpp_ +#define Hadrons_MIO_LoadPerambulator_hpp_ + +#include + +BEGIN_HADRONS_NAMESPACE +BEGIN_MODULE_NAMESPACE(MIO) + +/****************************************************************************** + * LoadPerambulator * + ******************************************************************************/ + +class LoadPerambulatorPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(LoadPerambulatorPar, + std::string, PerambFileName, + std::string, DistilParams); +}; + +template +class TLoadPerambulator: public Module +{ +public: + // constructor + TLoadPerambulator(const std::string name); + // destructor + virtual ~TLoadPerambulator(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(LoadPerambulator, TLoadPerambulator, MIO); + +/****************************************************************************** + * TLoadPerambulator implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TLoadPerambulator::TLoadPerambulator(const std::string name) : Module(name) {} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TLoadPerambulator::getInput(void) +{ + return {par().DistilParams}; +} + +template +std::vector TLoadPerambulator::getOutput(void) +{ + return {getName()}; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TLoadPerambulator::setup(void) +{ + const MDistil::DistilParameters &dp{envGet(MDistil::DistilParameters, par().DistilParams)}; + const int Nt{env().getDim(Tdir)}; + const bool full_tdil{ dp.TI == Nt }; + const int Nt_inv{ full_tdil ? 1 : dp.TI }; + envCreate(MDistil::PerambTensor, getName(), 1, Nt,dp.nvec,dp.LI,dp.nnoise,Nt_inv,dp.SI); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TLoadPerambulator::execute(void) +{ + auto &perambulator = envGet(MDistil::PerambTensor, getName()); + std::string sPerambName{ par().PerambFileName }; + sPerambName.append( 1, '.' ); + sPerambName.append( std::to_string( vm().getTrajectory() ) ); + perambulator.read(sPerambName.c_str()); +} + +END_MODULE_NAMESPACE +END_HADRONS_NAMESPACE +#endif diff --git a/Hadrons/Modules/MNPR/Amputate.hpp b/Hadrons/Modules/MNPR/Amputate.hpp index 6258e515..7cbe45ef 100644 --- a/Hadrons/Modules/MNPR/Amputate.hpp +++ b/Hadrons/Modules/MNPR/Amputate.hpp @@ -187,7 +187,7 @@ void TAmputate::execute(void) result.Vamp.resize(Gamma::nGamma/2); for( int mu=0; mu < Gamma::nGamma/2; mu++){ Gamma::Algebra gam = mu; - result.Vamp[mu] = 1/12.0*trace(adj(Gamma(mu*2+1))*g5*Sout_inv*g5*vertex[mu]*Sin_inv); + result.Vamp[mu] = static_cast( 1 / 12.0 ) * trace(adj(Gamma(mu*2+1))*g5*Sout_inv*g5*vertex[mu]*Sin_inv)()()(); LOG(Message) << "Vamp[" << mu << "] - " << result.Vamp[mu] << std::endl; } } diff --git a/Hadrons/NamedTensor.hpp b/Hadrons/NamedTensor.hpp new file mode 100644 index 00000000..954de2d8 --- /dev/null +++ b/Hadrons/NamedTensor.hpp @@ -0,0 +1,190 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: Hadrons/NamedTensor.hpp + + Copyright (C) 2015-2019 + + Author: Felix Erben + Author: Michael Marshall + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_NamedTensor_hpp_ +#define Hadrons_NamedTensor_hpp_ + +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + NamedTensor contains: + 1) Name of the tensor. By default, this is the tag name used for save / load + 2) Eigen::Tensor of type Scalar_ and rank NumIndices_ (row-major order) + 3) Name for each index + They can be persisted to / restored from disk. During restore, these validations are performed: + 1) Tensor dimensionality must match + 2) IndexNames are validated against current values + 3) If the tensor has non-zero size, the tensor being loaded must have same extent in each dimension + ******************************************************************************/ + +extern const std::string NamedTensorFileExtension; + +template +class NamedTensor : Serializable +{ +public: + using Scalar = Scalar_; + static constexpr int NumIndices = NumIndices_; + using ET = Eigen::Tensor; + using Index = typename ET::Index; + GRID_SERIALIZABLE_CLASS_MEMBERS(NamedTensor, + ET, tensor, + std::vector, IndexNames ); + + // Name of the object and Index names as set in the constructor + const std::string &Name_; + const std::array &DefaultIndexNames_; + + // Default constructor (assumes tensor will be loaded from file) + NamedTensor(const std::string &Name, + const std::array &indexNames) + : IndexNames{indexNames.begin(), indexNames.end()}, Name_{Name}, DefaultIndexNames_{indexNames} {} + + // Construct a named tensor explicitly specifying size of each dimension + template + NamedTensor(const std::string &Name, + const std::array &indexNames, + Eigen::Index firstDimension, IndexTypes... otherDimensions) + : tensor(firstDimension, otherDimensions...), + IndexNames{indexNames.begin(), indexNames.end()}, Name_{Name}, DefaultIndexNames_{indexNames} + { + if(sizeof...(otherDimensions) + 1 != NumIndices_) + { + HADRONS_ERROR(Argument, "NamedTensor: dimensions != tensor rank"); + } + } + + // Do my index names match the default for my type? + template + bool ValidateIndexNames( const array_or_vector_of_string &CheckNames ) const + { + return IndexNames.size() == CheckNames.size() && std::equal( IndexNames.begin(), IndexNames.end(), CheckNames.begin(), + [](const std::string &s1, const std::string &s2) + { + return s1.size() == s2.size() && std::equal( s1.begin(), s1.end(), s2.begin(), + [](const char & c1, const char & c2) + { return c1 == c2 || std::toupper(c1) == std::toupper(c2); }); // case insensitive + }); + } + bool ValidateIndexNames() const { return ValidateIndexNames(DefaultIndexNames_); } + +#ifdef HAVE_HDF5 + using Default_Reader = Grid::Hdf5Reader; + using Default_Writer = Grid::Hdf5Writer; +#else + using Default_Reader = Grid::BinaryReader; + using Default_Writer = Grid::BinaryWriter; +#endif + + void write(const std::string &FileName, const std::string &Tag) const + { + std::string FileName_{FileName}; + FileName_.append( NamedTensorFileExtension ); + LOG(Message) << "Writing " << Name_ << " to file " << FileName_ << " tag " << Tag << std::endl; + Default_Writer w( FileName_ ); + write( w, Tag, *this ); + } + void write(const std::string &FileName) const { return write(FileName, Name_); } + + // Read tensor. + // Validate: + // 1) index names (if requested) + // 2) index dimensions (if they are non-zero when called) + template void read(Reader &r, bool bValidate, const std::string &Tag) + { + // Grab index names and dimensions + std::vector OldIndexNames{std::move(IndexNames)}; + const typename ET::Dimensions OldDimensions{tensor.dimensions()}; + read(r, Tag, *this); + const typename ET::Dimensions & NewDimensions{tensor.dimensions()}; + for (int i = 0; i < NumIndices_; i++) + if(OldDimensions[i] && OldDimensions[i] != NewDimensions[i]) + { + HADRONS_ERROR(Size,"NamedTensor::read dimension size"); + } + if (bValidate && !ValidateIndexNames(OldIndexNames)) + { + HADRONS_ERROR(Definition,"NamedTensor::read dimension name"); + } + } + template void read(Reader &r, bool bValidate = true) { read(r, bValidate, Name_); } + + inline void read (const std::string &FileName, bool bValidate, const std::string &Tag) + { + Default_Reader r(FileName + NamedTensorFileExtension); + read(r, bValidate, Tag); + } + inline void read (const std::string &FileName, bool bValidate= true) { return read(FileName, bValidate, Name_); } +}; + +/****************************************************************************** + Common elements for distillation + ******************************************************************************/ + +BEGIN_MODULE_NAMESPACE(MDistil) + +//Eigenvectors of the Laplacian +using LapEvecs = Grid::Hadrons::EigenPack; + +// Noise vector (index order: nnoise, nt, nvec, ns) + +class NoiseTensor : public NamedTensor +{ + static const std::string Name__; + static const std::array DefaultIndexNames__; + public: + // Default constructor (assumes tensor will be loaded from file) + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NoiseTensor() : NamedTensor{Name__, DefaultIndexNames__} {} + + // Construct a named tensor explicitly specifying size of each dimension + template + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NoiseTensor(Eigen::Index nNoise, Eigen::Index nT, Eigen::Index nVec, Eigen::Index nS) + : NamedTensor{Name__, DefaultIndexNames__, nNoise, nT, nVec, nS} {} +}; + +class PerambTensor : public NamedTensor +{ + static const std::string Name__; + static const std::array DefaultIndexNames__; + public: + // Default constructor (assumes tensor will be loaded from file) + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE PerambTensor() : NamedTensor{Name__, DefaultIndexNames__} {} + + // Construct a named tensor explicitly specifying size of each dimension + template + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE PerambTensor(Eigen::Index nT, Eigen::Index nVec, Eigen::Index LI, Eigen::Index nNoise, Eigen::Index nT_inv, Eigen::Index SI) + : NamedTensor{Name__, DefaultIndexNames__, nT, nVec, LI, nNoise, nT_inv, SI} {} +}; + +END_MODULE_NAMESPACE +END_HADRONS_NAMESPACE +#endif // Hadrons_NamedTensor_hpp_ diff --git a/Hadrons/modules.inc b/Hadrons/modules.inc index bcf58d00..79bfba0a 100644 --- a/Hadrons/modules.inc +++ b/Hadrons/modules.inc @@ -48,7 +48,9 @@ modules_cc =\ Modules/MNoise/FullVolumeSpinColorDiagonal.cc \ Modules/MIO/LoadA2AVectors.cc \ Modules/MIO/LoadBinary.cc \ + Modules/MIO/LoadPerambulator.cc \ Modules/MIO/LoadA2AMatrixDiskVector.cc \ + Modules/MIO/LoadDistilNoise.cc \ Modules/MIO/LoadNersc.cc \ Modules/MIO/LoadCoarseEigenPack.cc \ Modules/MIO/LoadCosmHol.cc \ @@ -63,6 +65,12 @@ modules_cc =\ Modules/MScalarSUN/TrKinetic.cc \ Modules/MScalarSUN/StochFreeField.cc \ Modules/MScalarSUN/TrMag.cc \ + Modules/MDistil/PerambFromSolve.cc \ + Modules/MDistil/DistilVectors.cc \ + Modules/MDistil/DistilPar.cc \ + Modules/MDistil/LapEvec.cc \ + Modules/MDistil/Perambulator.cc \ + Modules/MDistil/Noises.cc \ Modules/MAction/MobiusDWF.cc \ Modules/MAction/DWF.cc \ Modules/MAction/ScaledDWF.cc \ @@ -126,8 +134,10 @@ modules_hpp =\ Modules/MIO/LoadBinary.hpp \ Modules/MIO/LoadNersc.hpp \ Modules/MIO/LoadEigenPack.hpp \ + Modules/MIO/LoadDistilNoise.hpp \ Modules/MIO/LoadA2AMatrixDiskVector.hpp \ Modules/MIO/LoadCoarseEigenPack.hpp \ + Modules/MIO/LoadPerambulator.hpp \ Modules/MIO/LoadA2AVectors.hpp \ Modules/MIO/LoadCosmHol.hpp \ Modules/MScalarSUN/TwoPointNPR.hpp \ @@ -141,6 +151,13 @@ modules_hpp =\ Modules/MScalarSUN/StochFreeField.hpp \ Modules/MScalarSUN/Div.hpp \ Modules/MScalarSUN/TrKinetic.hpp \ + Modules/MDistil/PerambFromSolve.hpp \ + Modules/MDistil/Distil.hpp \ + Modules/MDistil/DistilVectors.hpp \ + Modules/MDistil/Noises.hpp \ + Modules/MDistil/LapEvec.hpp \ + Modules/MDistil/Perambulator.hpp \ + Modules/MDistil/DistilPar.hpp \ Modules/MAction/ZMobiusDWF.hpp \ Modules/MAction/ScaledDWF.hpp \ Modules/MAction/WilsonClover.hpp \ diff --git a/configure.ac b/configure.ac index 640b2dbf..d451c163 100644 --- a/configure.ac +++ b/configure.ac @@ -136,15 +136,15 @@ case ${ac_SFW_FP16} in esac ############### SUMMIT JSRUN -AC_ARG_ENABLE([jsrun], - [AC_HELP_STRING([--enable-jsrun=yes|no], [enable IBMs jsrun resource manager for SUMMIT])], - [ac_JSRUN=${enable_jsrun}], [ac_JSRUN=no]) -case ${ac_JSRUN} in +AC_ARG_ENABLE([summit], + [AC_HELP_STRING([--enable-summit=yes|no], [enable IBMs jsrun resource manager for SUMMIT])], + [ac_JSRUN=${enable_summit}], [ac_SUMMIT=no]) +case ${ac_SUMMIT} in + no);; yes) AC_DEFINE([GRID_IBM_SUMMIT],[1],[Let JSRUN manage the GPU device allocation]);; - no);; *) - AC_MSG_ERROR(["JSRUN option not supported ${ac_JSRUN}"]);; + AC_DEFINE([GRID_IBM_SUMMIT],[1],[Let JSRUN manage the GPU device allocation]);; esac ############### Intel libraries @@ -255,7 +255,7 @@ AC_ARG_ENABLE([simd],[AC_HELP_STRING([--enable-simd=code], AC_ARG_ENABLE([gen-simd-width], [AS_HELP_STRING([--enable-gen-simd-width=size], - [size (in bytes) of the generic SIMD vectors (default: 32)])], + [size (in bytes) of the generic SIMD vectors (default: 64)])], [ac_gen_simd_width=$enable_gen_simd_width], [ac_gen_simd_width=64]) diff --git a/documentation/GridXcode/readme.md b/documentation/GridXcode/readme.md index 8d9d7ad8..8cbf9112 100644 --- a/documentation/GridXcode/readme.md +++ b/documentation/GridXcode/readme.md @@ -119,13 +119,13 @@ Install [MacPorts][MacPorts] if you haven't done so already, and then install pa These are the `portname`s for mandatory Grid libraries: -* git +* git-flow-avh * gmp * mpfr and these are the `portname`s for optional Grid libraries: -* fftw-3 +* fftw-3-single * hdf5 * lapack * doxygen @@ -369,7 +369,7 @@ Instead: 3. From a terminal session, locate and run your executable using `mpirun` (*the mangled name of the project build products will not be exactly the same as this example*): - `$GridPre/openmpi-3.1.3/bin/mpirun -np 2 ~/Library/Developer/Xcode/DerivedData/HelloGrid-fiyyuveptaqelbbvllomcgjyvghr/Build/Products/Debug/HelloGrid --grid 4.4.4.8 --mpi 1.1.1.2` + `$GridPre/openmpi/bin/mpirun -np 2 ~/Library/Developer/Xcode/DerivedData/HelloGrid-fiyyuveptaqelbbvllomcgjyvghr/Build/Products/Debug/HelloGrid --grid 4.4.4.8 --mpi 1.1.1.2` The Xcode debugger will attach to the first process. diff --git a/tests/hadrons/Test_distil.cc b/tests/hadrons/Test_distil.cc new file mode 100644 index 00000000..0156f340 --- /dev/null +++ b/tests/hadrons/Test_distil.cc @@ -0,0 +1,382 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: Tests/Hadrons/Test_hadrons_distil.cc + + Copyright (C) 2015-2019 + + Author: Felix Erben + Author: Michael Marshall + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ +/* END LEGAL */ + +#include +#include +#include + +using namespace Grid; +using namespace Hadrons; + +// Very simple iterators for Eigen tensors +// The only way I could get these iterators to work is to put the begin() and end() functions in the Eigen namespace +// So if Eigen ever defines these, we'll have a conflict and have to change this +namespace Eigen { + template + inline typename std::enable_if::value, typename EigenIO::Traits::scalar_type *>::type + begin( ET & et ) { return reinterpret_cast::scalar_type *>(et.data()); } + template + inline typename std::enable_if::value, typename EigenIO::Traits::scalar_type *>::type + end( ET & et ) { return begin(et) + et.size() * EigenIO::Traits::count; } +} + +///////////////////////////////////////////////////////////// +// Test creation of laplacian eigenvectors +///////////////////////////////////////////////////////////// + +void test_Global(Application &application) +{ + // global parameters + Application::GlobalPar globalPar; + globalPar.trajCounter.start = 1100; + globalPar.trajCounter.end = 1120; + globalPar.trajCounter.step = 20; + globalPar.runId = "test"; + globalPar.graphFile = ""; + globalPar.scheduleFile = ""; + globalPar.saveSchedule = "false"; + globalPar.parallelWriteMaxRetry = -1; + application.setPar(globalPar); +} + +///////////////////////////////////////////////////////////// +// Create a random gauge with the correct name +///////////////////////////////////////////////////////////// + +std::string test_Gauge(Application &application ) +{ + std::string sGaugeName{ "gauge" }; + application.createModule( sGaugeName ); + return sGaugeName; +} + +///////////////////////////////////////////////////////////// +// Test creation of laplacian eigenvectors +///////////////////////////////////////////////////////////// + +void test_LapEvec(Application &application) +{ + const char szModuleName[] = "LapEvec"; + test_Gauge( application ); + MDistil::LapEvecPar p; + p.gauge = "gauge"; + p.Stout.steps = 3; + p.Stout.rho = 0.2; + p.Cheby.PolyOrder = 11; + p.Cheby.alpha = 0.55; + p.Cheby.beta = 35.5; + p.Lanczos.Nvec = 5; + p.Lanczos.Nk = 6; + p.Lanczos.Np = 2; + p.Lanczos.MaxIt = 1000; + p.Lanczos.resid = 1e-2; + p.Lanczos.IRLLog = 0; + application.createModule(szModuleName,p); +} + +///////////////////////////////////////////////////////////// +// Test creation Solver +///////////////////////////////////////////////////////////// + +std::string SolverName( const char * pSuffix = nullptr ) { + std::string sSolverName{ "CG" }; + if( pSuffix && pSuffix[0] ) { + sSolverName.append( "_" ); + sSolverName.append( pSuffix ); + } + return sSolverName; +} + +std::string test_Solver(Application &application, const char * pSuffix = nullptr ) +{ + std::string sActionName{ "DWF" }; + if( pSuffix && pSuffix[0] ) { + sActionName.append( "_" ); + sActionName.append( pSuffix ); + } + MAction::DWF::Par actionPar; + actionPar.gauge = "gauge"; + actionPar.Ls = 16; + actionPar.M5 = 1.8; + actionPar.mass = 0.005; + actionPar.boundary = "1 1 1 -1"; + actionPar.twist = "0. 0. 0. 0."; + application.createModule( sActionName, actionPar ); + MSolver::RBPrecCG::Par solverPar; + solverPar.action = sActionName; + solverPar.residual = 1.0e-2; + solverPar.maxIteration = 10000; + std::string sSolverName{ SolverName( pSuffix ) }; + application.createModule( sSolverName, solverPar ); + return sSolverName; +} + +///////////////////////////////////////////////////////////// +// DistilParameters +///////////////////////////////////////////////////////////// + +std::string test_DPar(Application &application) { + MDistil::DistilParameters DPar; + DPar.nvec = 5; + DPar.nnoise = 1; + DPar.tsrc = 0; + DPar.LI = 5; + DPar.TI = 8; + DPar.SI = 4; + std::string sDParName{"DPar_l"}; + application.createModule(sDParName,DPar); + return sDParName; +} +///////////////////////////////////////////////////////////// +// Noises +///////////////////////////////////////////////////////////// + +std::string test_Noises(Application &application, const std::string &sNoiseBaseName ) { + MDistil::NoisesPar NoisePar; + NoisePar.DistilParams = "DPar_l"; + NoisePar.NoiseFileName = "noise"; + std::string sNoiseName{"noise"}; + application.createModule(sNoiseName,NoisePar); + return sNoiseName; +} + +///////////////////////////////////////////////////////////// +// Perambulators +///////////////////////////////////////////////////////////// + +std::string PerambulatorName( const char * pszSuffix = nullptr ) +{ + std::string sPerambulatorName{ "Peramb" }; + if( pszSuffix && pszSuffix[0] ) + sPerambulatorName.append( pszSuffix ); + return sPerambulatorName; +} + +void test_LoadPerambulators( Application &application, const char * pszSuffix = nullptr ) +{ + std::string sModuleName{ "Peramb_load" }; + MIO::LoadPerambulator::Par PerambPar; + PerambPar.PerambFileName = "Peramb"; + PerambPar.DistilParams = "DPar_l"; + test_Noises(application, sModuleName); // I want these written after solver stuff + application.createModule( sModuleName, PerambPar ); +} + +void test_Perambulators( Application &application, const char * pszSuffix = nullptr ) +{ + std::string sModuleName{ PerambulatorName( pszSuffix ) }; + // Perambulator parameters + MDistil::Perambulator::Par PerambPar; + PerambPar.lapevec = "LapEvec"; + PerambPar.PerambFileName = sModuleName; + PerambPar.solver = test_Solver( application, pszSuffix ); + PerambPar.DistilParams = "DPar_l"; + PerambPar.noise = "noise"; + test_Noises(application, sModuleName); // I want these written after solver stuff + application.createModule( sModuleName, PerambPar ); +} + +///////////////////////////////////////////////////////////// +// DistilVectors +///////////////////////////////////////////////////////////// + +void test_DistilVectors(Application &application, const char * pszSuffix = nullptr, const char * pszNvec = nullptr ) +{ + std::string sModuleName{"DistilVecs"}; + if( pszSuffix ) + sModuleName.append( pszSuffix ); + std::string sPerambName{"Peramb"}; + if( pszSuffix ) + sPerambName.append( pszSuffix ); + MDistil::DistilVectors::Par DistilVecPar; + DistilVecPar.noise = "noise"; + DistilVecPar.rho = "rho"; + DistilVecPar.phi = "phi"; + DistilVecPar.perambulator = sPerambName; + DistilVecPar.lapevec = "LapEvec"; + DistilVecPar.DistilParams = "DPar_l"; + application.createModule(sModuleName,DistilVecPar); +} + +///////////////////////////////////////////////////////////// +// MesonSink +///////////////////////////////////////////////////////////// + +void test_MesonSink(Application &application) +{ + // DistilVectors parameters + MContraction::A2AMesonField::Par A2AMesonFieldPar; + //A2AMesonFieldPar.left="Peramb_unsmeared_sink"; + A2AMesonFieldPar.left="g5phi"; + A2AMesonFieldPar.right="Peramb_unsmeared_sink"; + A2AMesonFieldPar.output="DistilFields"; + A2AMesonFieldPar.gammas="Identity"; + A2AMesonFieldPar.mom={"0 0 0"}; + A2AMesonFieldPar.cacheBlock=2; + A2AMesonFieldPar.block=4; + application.createModule("DistilMesonSink",A2AMesonFieldPar); +} +///////////////////////////////////////////////////////////// +// MesonFields +///////////////////////////////////////////////////////////// + +void test_MesonField(Application &application, const char * pszFileSuffix, + const char * pszObjectLeft = nullptr, const char * pszObjectRight = nullptr ) +{ + // DistilVectors parameters + if( pszObjectLeft == nullptr ) + pszObjectLeft = pszFileSuffix; + if( pszObjectRight == nullptr ) + pszObjectRight = pszObjectLeft; + MContraction::A2AMesonField::Par A2AMesonFieldPar; + A2AMesonFieldPar.left=""; + A2AMesonFieldPar.right=A2AMesonFieldPar.left; + A2AMesonFieldPar.left.append( pszObjectLeft ); + A2AMesonFieldPar.right.append( pszObjectRight ); + A2AMesonFieldPar.output="MesonSinks"; + A2AMesonFieldPar.output.append( pszFileSuffix ); + A2AMesonFieldPar.gammas="Identity"; + A2AMesonFieldPar.mom={"0 0 0"}; + A2AMesonFieldPar.cacheBlock=2; + A2AMesonFieldPar.block=4; + std::string sObjectName{"DistilMesonField"}; + sObjectName.append( pszFileSuffix ); + application.createModule(sObjectName, A2AMesonFieldPar); +} + +bool bNumber( int &ri, const char * & pstr, bool bGobbleWhiteSpace = true ) +{ + if( bGobbleWhiteSpace ) + while( std::isspace(static_cast(*pstr)) ) + pstr++; + const char * p = pstr; + bool bMinus = false; + char c = * p++; + if( c == '+' ) + c = * p++; + else if( c == '-' ) { + bMinus = true; + c = * p++; + } + int n = c - '0'; + if( n < 0 || n > 9 ) + return false; + while( * p >= '0' && * p <= '9' ) { + n = n * 10 + ( * p ) - '0'; + p++; + } + if( bMinus ) + n *= -1; + ri = n; + pstr = p; + return true; +} + +int main(int argc, char *argv[]) +{ + // Decode command-line parameters. 1st one is which test to run + int iTestNum = -1; + + for(int i = 1 ; i < argc ; i++ ) { + std::cout << "argv[" << i << "]=\"" << argv[i] << "\"" << std::endl; + const char * p = argv[i]; + if( * p == '/' || * p == '-' ) { + p++; + char c = * p++; + switch(toupper(c)) { + case 'T': + if( bNumber( iTestNum, p ) ) { + std::cout << "Test " << iTestNum << " requested"; + if( * p ) + std::cout << " (ignoring trailer \"" << p << "\")"; + std::cout << std::endl; + } + else + std::cout << "Invalid test \"" << &argv[i][2] << "\"" << std::endl; + break; + default: + std::cout << "Ignoring switch \"" << &argv[i][1] << "\"" << std::endl; + break; + } + } + } + + // initialization ////////////////////////////////////////////////////////// + Grid_init(&argc, &argv); + HadronsLogError.Active(GridLogError.isActive()); + HadronsLogWarning.Active(GridLogWarning.isActive()); + HadronsLogMessage.Active(GridLogMessage.isActive()); + HadronsLogIterative.Active(GridLogIterative.isActive()); + HadronsLogDebug.Active(GridLogDebug.isActive()); + LOG(Message) << "Grid initialized" << std::endl; + + // run setup /////////////////////////////////////////////////////////////// + Application application; + + // For now perform free propagator test - replace this with distillation test(s) + LOG(Message) << "====== Creating xml for test " << iTestNum << " ======" << std::endl; + + switch(iTestNum) { + default: // 0 + LOG(Message) << "Computing Meson 2pt-function" << std::endl; + test_Global( application ); + test_LapEvec( application ); + test_DPar( application ); + test_Perambulators( application ); + test_DistilVectors( application ); + test_MesonField( application, "Phi", "phi" ); + test_MesonField( application, "Rho", "rho" ); + break; + case 1: + LOG(Message) << "Computing Meson 2pt-function by loading perambulators" << std::endl; + test_Global( application ); + test_LapEvec( application ); + test_DPar( application ); + test_LoadPerambulators( application ); + test_DistilVectors( application, "_load" ); + test_MesonField( application, "Phi", "phi" ); + test_MesonField( application, "Rho", "rho" ); + break; + } + // execution + static const char XmlFileName[] = "test_distil.xml"; + application.saveParameterFile( XmlFileName ); + + const Grid::Coordinate &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; + Grid_finalize(); + + return EXIT_SUCCESS; +}