From d1f7c6b94e732d7619c1e2862264e2dfe8884475 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Mon, 5 Dec 2016 16:47:29 +0900 Subject: [PATCH] Hadrons: templatisation of the fermion implementation --- extras/Hadrons/Global.hpp | 17 +++- extras/Hadrons/Modules.hpp | 2 +- extras/Hadrons/Modules/MAction/DWF.cc | 82 ---------------- extras/Hadrons/Modules/MAction/DWF.hpp | 66 ++++++++++++- extras/Hadrons/Modules/MAction/Wilson.cc | 76 -------------- extras/Hadrons/Modules/MAction/Wilson.hpp | 64 +++++++++++- extras/Hadrons/Modules/MContraction/Meson.cc | 94 ------------------ extras/Hadrons/Modules/MContraction/Meson.hpp | 81 ++++++++++++++- extras/Hadrons/Modules/MQuark.hpp | 71 -------------- extras/Hadrons/Modules/MSolver/RBPrecCG.cc | 84 ---------------- extras/Hadrons/Modules/MSolver/RBPrecCG.hpp | 68 ++++++++++++- extras/Hadrons/Modules/MSource/Point.cc | 75 -------------- extras/Hadrons/Modules/MSource/Point.hpp | 63 +++++++++++- extras/Hadrons/Modules/MSource/Z2.cc | 88 ----------------- extras/Hadrons/Modules/MSource/Z2.hpp | 75 +++++++++++++- .../Hadrons/Modules/{MQuark.cc => Quark.hpp} | 98 ++++++++++++++----- extras/Hadrons/modules.inc | 13 +-- tests/hadrons/Test_hadrons_mesons.cc | 4 +- 18 files changed, 483 insertions(+), 638 deletions(-) delete mode 100644 extras/Hadrons/Modules/MAction/DWF.cc delete mode 100644 extras/Hadrons/Modules/MAction/Wilson.cc delete mode 100644 extras/Hadrons/Modules/MContraction/Meson.cc delete mode 100644 extras/Hadrons/Modules/MQuark.hpp delete mode 100644 extras/Hadrons/Modules/MSolver/RBPrecCG.cc delete mode 100644 extras/Hadrons/Modules/MSource/Point.cc delete mode 100644 extras/Hadrons/Modules/MSource/Z2.cc rename extras/Hadrons/Modules/{MQuark.cc => Quark.hpp} (58%) diff --git a/extras/Hadrons/Global.hpp b/extras/Hadrons/Global.hpp index 1dc9164a..3107ef12 100644 --- a/extras/Hadrons/Global.hpp +++ b/extras/Hadrons/Global.hpp @@ -56,10 +56,19 @@ using Grid::operator<<; BEGIN_HADRONS_NAMESPACE // type aliases -typedef FermionOperator FMat; -typedef FIMPL::FermionField FermionField; -typedef FIMPL::PropagatorField PropagatorField; -typedef std::function SolverFn; +//typedef FermionOperator FMat; +//typedef FIMPL::FermionField FermionField; +//typedef FIMPL::PropagatorField PropagatorField; +//typedef std::function SolverFn; + +#define TYPE_ALIASES(FImpl, suffix)\ +typedef FermionOperator FMat##suffix; \ +typedef typename FImpl::FermionField FermionField##suffix; \ +typedef typename FImpl::PropagatorField PropagatorField##suffix; \ +typedef typename FImpl::SitePropagator SitePropagator##suffix; \ +typedef typename FImpl::DoubledGaugeField DoubledGaugeField##suffix;\ +typedef std::function SolverFn##suffix; // logger class HadronsLogger: public Logger diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 998c609f..e90d0313 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -4,7 +4,7 @@ #include #include #include -#include #include #include #include +#include diff --git a/extras/Hadrons/Modules/MAction/DWF.cc b/extras/Hadrons/Modules/MAction/DWF.cc deleted file mode 100644 index 36527b58..00000000 --- a/extras/Hadrons/Modules/MAction/DWF.cc +++ /dev/null @@ -1,82 +0,0 @@ -/******************************************************************************* -Grid physics library, www.github.com/paboyle/Grid - -Source file: programs/Hadrons/DWF.cc - -Copyright (C) 2016 - -Author: Antonin Portelli - -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. -*******************************************************************************/ - -#include - -using namespace Grid; -using namespace Hadrons; -using namespace MAction; - -/****************************************************************************** -* DWF implementation * -******************************************************************************/ -// constructor ///////////////////////////////////////////////////////////////// -DWF::DWF(const std::string name) -: Module(name) -{} - -// dependencies/products /////////////////////////////////////////////////////// -std::vector DWF::getInput(void) -{ - std::vector in = {par().gauge}; - - return in; -} - -std::vector DWF::getOutput(void) -{ - std::vector out = {getName()}; - - return out; -} - -// setup /////////////////////////////////////////////////////////////////////// -void DWF::setup(void) -{ - unsigned int size; - - size = 3*env().lattice4dSize(); - env().registerObject(getName(), size, par().Ls); -} - -// execution /////////////////////////////////////////////////////////////////// -void DWF::execute(void) -{ - LOG(Message) << "Setting up domain wall fermion matrix with m= " - << par().mass << ", M5= " << par().M5 << " and Ls= " - << par().Ls << " using gauge field '" << par().gauge << "'" - << std::endl; - env().createGrid(par().Ls); - auto &U = *env().getObject(par().gauge); - auto &g4 = *env().getGrid(); - auto &grb4 = *env().getRbGrid(); - auto &g5 = *env().getGrid(par().Ls); - auto &grb5 = *env().getRbGrid(par().Ls); - FMat *fMatPt = new DomainWallFermion(U, g5, grb5, g4, grb4, - par().mass, par().M5); - env().setObject(getName(), fMatPt); -} diff --git a/extras/Hadrons/Modules/MAction/DWF.hpp b/extras/Hadrons/Modules/MAction/DWF.hpp index 246fafd4..ceec975c 100644 --- a/extras/Hadrons/Modules/MAction/DWF.hpp +++ b/extras/Hadrons/Modules/MAction/DWF.hpp @@ -49,13 +49,16 @@ public: double , M5); }; -class DWF: public Module +template +class TDWF: public Module { +public: + TYPE_ALIASES(FImpl,); public: // constructor - DWF(const std::string name); + TDWF(const std::string name); // destructor - virtual ~DWF(void) = default; + virtual ~TDWF(void) = default; // dependency relation virtual std::vector getInput(void); virtual std::vector getOutput(void); @@ -65,6 +68,63 @@ public: virtual void execute(void); }; +/****************************************************************************** + * DWF template implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TDWF::TDWF(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TDWF::getInput(void) +{ + std::vector in = {par().gauge}; + + return in; +} + +template +std::vector TDWF::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TDWF::setup(void) +{ + unsigned int size; + + size = 3*env().template lattice4dSize(); + env().registerObject(getName(), size, par().Ls); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TDWF::execute(void) +{ + LOG(Message) << "Setting up domain wall fermion matrix with m= " + << par().mass << ", M5= " << par().M5 << " and Ls= " + << par().Ls << " using gauge field '" << par().gauge << "'" + << std::endl; + env().createGrid(par().Ls); + auto &U = *env().template getObject(par().gauge); + auto &g4 = *env().getGrid(); + auto &grb4 = *env().getRbGrid(); + auto &g5 = *env().getGrid(par().Ls); + auto &grb5 = *env().getRbGrid(par().Ls); + FMat *fMatPt = new DomainWallFermion(U, g5, grb5, g4, grb4, + par().mass, par().M5); + env().setObject(getName(), fMatPt); +} + +typedef TDWF DWF; + END_MODULE_NAMESPACE MODULE_REGISTER_NS(DWF, MAction); diff --git a/extras/Hadrons/Modules/MAction/Wilson.cc b/extras/Hadrons/Modules/MAction/Wilson.cc deleted file mode 100644 index 3188a636..00000000 --- a/extras/Hadrons/Modules/MAction/Wilson.cc +++ /dev/null @@ -1,76 +0,0 @@ -/******************************************************************************* -Grid physics library, www.github.com/paboyle/Grid - -Source file: programs/Hadrons/Wilson.cc - -Copyright (C) 2016 - -Author: Antonin Portelli - -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. -*******************************************************************************/ - -#include - -using namespace Grid; -using namespace Hadrons; -using namespace MAction; - -/****************************************************************************** -* Wilson implementation * -******************************************************************************/ -// constructor ///////////////////////////////////////////////////////////////// -Wilson::Wilson(const std::string name) -: Module(name) -{} - -// dependencies/products /////////////////////////////////////////////////////// -std::vector Wilson::getInput(void) -{ - std::vector in = {par().gauge}; - - return in; -} - -std::vector Wilson::getOutput(void) -{ - std::vector out = {getName()}; - - return out; -} - -// setup /////////////////////////////////////////////////////////////////////// -void Wilson::setup(void) -{ - unsigned int size; - - size = 3*env().lattice4dSize(); - env().registerObject(getName(), size); -} - -// execution /////////////////////////////////////////////////////////////////// -void Wilson::execute() -{ - LOG(Message) << "Setting up Wilson fermion matrix with m= " << par().mass - << " using gauge field '" << par().gauge << "'" << std::endl; - auto &U = *env().getObject(par().gauge); - auto &grid = *env().getGrid(); - auto &gridRb = *env().getRbGrid(); - FMat *fMatPt = new WilsonFermion(U, grid, gridRb, par().mass); - env().setObject(getName(), fMatPt); -} diff --git a/extras/Hadrons/Modules/MAction/Wilson.hpp b/extras/Hadrons/Modules/MAction/Wilson.hpp index d3965f09..6a4528ca 100644 --- a/extras/Hadrons/Modules/MAction/Wilson.hpp +++ b/extras/Hadrons/Modules/MAction/Wilson.hpp @@ -1,7 +1,7 @@ /******************************************************************************* Grid physics library, www.github.com/paboyle/Grid -Source file: programs/Hadrons/Wilson.hpp +Source file: programs/Hadrons/TWilson.hpp Copyright (C) 2016 @@ -35,7 +35,7 @@ directory. BEGIN_HADRONS_NAMESPACE /****************************************************************************** - * Wilson quark action * + * TWilson quark action * ******************************************************************************/ BEGIN_MODULE_NAMESPACE(MAction) @@ -47,13 +47,16 @@ public: double , mass); }; -class Wilson: public Module +template +class TWilson: public Module { +public: + TYPE_ALIASES(FImpl,); public: // constructor - Wilson(const std::string name); + TWilson(const std::string name); // destructor - virtual ~Wilson(void) = default; + virtual ~TWilson(void) = default; // dependencies/products virtual std::vector getInput(void); virtual std::vector getOutput(void); @@ -63,6 +66,57 @@ public: virtual void execute(void); }; +/****************************************************************************** + * TWilson template implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TWilson::TWilson(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TWilson::getInput(void) +{ + std::vector in = {par().gauge}; + + return in; +} + +template +std::vector TWilson::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TWilson::setup(void) +{ + unsigned int size; + + size = 3*env().template lattice4dSize(); + env().registerObject(getName(), size); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TWilson::execute() +{ + LOG(Message) << "Setting up TWilson fermion matrix with m= " << par().mass + << " using gauge field '" << par().gauge << "'" << std::endl; + auto &U = *env().template getObject(par().gauge); + auto &grid = *env().getGrid(); + auto &gridRb = *env().getRbGrid(); + FMat *fMatPt = new WilsonFermion(U, grid, gridRb, par().mass); + env().setObject(getName(), fMatPt); +} + +typedef TWilson Wilson; + END_MODULE_NAMESPACE MODULE_REGISTER_NS(Wilson, MAction); diff --git a/extras/Hadrons/Modules/MContraction/Meson.cc b/extras/Hadrons/Modules/MContraction/Meson.cc deleted file mode 100644 index fbdfceb1..00000000 --- a/extras/Hadrons/Modules/MContraction/Meson.cc +++ /dev/null @@ -1,94 +0,0 @@ -/******************************************************************************* -Grid physics library, www.github.com/paboyle/Grid - -Source file: programs/Hadrons/Meson.cc - -Copyright (C) 2015 - -Author: Antonin Portelli - -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. -*******************************************************************************/ - -#include - -using namespace Grid; -using namespace QCD; -using namespace Hadrons; -using namespace MContraction; - -/****************************************************************************** - * Meson implementation * - ******************************************************************************/ -// constructor ///////////////////////////////////////////////////////////////// -Meson::Meson(const std::string name) -: Module(name) -{} - -// dependencies/products /////////////////////////////////////////////////////// -std::vector Meson::getInput(void) -{ - std::vector input = {par().q1, par().q2}; - - return input; -} - -std::vector Meson::getOutput(void) -{ - std::vector output = {getName()}; - - return output; -} - -// execution /////////////////////////////////////////////////////////////////// -void Meson::execute(void) -{ - LOG(Message) << "Computing meson contraction '" << getName() << "' using" - << " quarks '" << par().q1 << "' and '" << par().q2 << "'" - << std::endl; - - XmlWriter writer(par().output); - PropagatorField &q1 = *env().getObject(par().q1); - PropagatorField &q2 = *env().getObject(par().q2); - LatticeComplex c(env().getGrid()); - SpinMatrix g[Ns*Ns], g5; - std::vector buf; - Result result; - - g5 = makeGammaProd(Ns*Ns - 1); - result.corr.resize(Ns*Ns); - for (unsigned int i = 0; i < Ns*Ns; ++i) - { - g[i] = makeGammaProd(i); - } - for (unsigned int iSink = 0; iSink < Ns*Ns; ++iSink) - { - result.corr[iSink].resize(Ns*Ns); - for (unsigned int iSrc = 0; iSrc < Ns*Ns; ++iSrc) - { - c = trace(g[iSink]*q1*g[iSrc]*g5*adj(q2)*g5); - sliceSum(c, buf, Tp); - result.corr[iSink][iSrc].resize(buf.size()); - for (unsigned int t = 0; t < buf.size(); ++t) - { - result.corr[iSink][iSrc][t] = TensorRemove(buf[t]); - } - } - } - write(writer, "meson", result); -} diff --git a/extras/Hadrons/Modules/MContraction/Meson.hpp b/extras/Hadrons/Modules/MContraction/Meson.hpp index 713c8631..fde05fca 100644 --- a/extras/Hadrons/Modules/MContraction/Meson.hpp +++ b/extras/Hadrons/Modules/MContraction/Meson.hpp @@ -1,7 +1,7 @@ /******************************************************************************* Grid physics library, www.github.com/paboyle/Grid -Source file: programs/Hadrons/Meson.hpp +Source file: programs/Hadrons/TMeson.hpp Copyright (C) 2015 @@ -35,7 +35,7 @@ directory. BEGIN_HADRONS_NAMESPACE /****************************************************************************** - * Meson * + * TMeson * ******************************************************************************/ BEGIN_MODULE_NAMESPACE(MContraction) @@ -48,8 +48,12 @@ public: std::string, output); }; -class Meson: public Module +template +class TMeson: public Module { +public: + TYPE_ALIASES(FImpl1, 1); + TYPE_ALIASES(FImpl2, 2); public: class Result: Serializable { @@ -59,9 +63,9 @@ public: }; public: // constructor - Meson(const std::string name); + TMeson(const std::string name); // destructor - virtual ~Meson(void) = default; + virtual ~TMeson(void) = default; // dependencies/products virtual std::vector getInput(void); virtual std::vector getOutput(void); @@ -69,6 +73,73 @@ public: virtual void execute(void); }; +/****************************************************************************** + * TMeson implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TMeson::TMeson(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TMeson::getInput(void) +{ + std::vector input = {par().q1, par().q2}; + + return input; +} + +template +std::vector TMeson::getOutput(void) +{ + std::vector output = {getName()}; + + return output; +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TMeson::execute(void) +{ + LOG(Message) << "Computing meson contraction '" << getName() << "' using" + << " quarks '" << par().q1 << "' and '" << par().q2 << "'" + << std::endl; + + XmlWriter writer(par().output); + PropagatorField1 &q1 = *env().template getObject(par().q1); + PropagatorField2 &q2 = *env().template getObject(par().q2); + LatticeComplex c(env().getGrid()); + SpinMatrix g[Ns*Ns], g5; + std::vector buf; + Result result; + + g5 = makeGammaProd(Ns*Ns - 1); + result.corr.resize(Ns*Ns); + for (unsigned int i = 0; i < Ns*Ns; ++i) + { + g[i] = makeGammaProd(i); + } + for (unsigned int iSink = 0; iSink < Ns*Ns; ++iSink) + { + result.corr[iSink].resize(Ns*Ns); + for (unsigned int iSrc = 0; iSrc < Ns*Ns; ++iSrc) + { + c = trace(g[iSink]*q1*g[iSrc]*g5*adj(q2)*g5); + sliceSum(c, buf, Tp); + result.corr[iSink][iSrc].resize(buf.size()); + for (unsigned int t = 0; t < buf.size(); ++t) + { + result.corr[iSink][iSrc][t] = TensorRemove(buf[t]); + } + } + } + write(writer, "meson", result); +} + +typedef TMeson Meson; + END_MODULE_NAMESPACE MODULE_REGISTER_NS(Meson, MContraction); diff --git a/extras/Hadrons/Modules/MQuark.hpp b/extras/Hadrons/Modules/MQuark.hpp deleted file mode 100644 index ed6dfb0c..00000000 --- a/extras/Hadrons/Modules/MQuark.hpp +++ /dev/null @@ -1,71 +0,0 @@ -/******************************************************************************* -Grid physics library, www.github.com/paboyle/Grid - -Source file: programs/Hadrons/MQuark.hpp - -Copyright (C) 2015 - -Author: Antonin Portelli - -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. -*******************************************************************************/ - -#ifndef Hadrons_MQuark_hpp_ -#define Hadrons_MQuark_hpp_ - -#include -#include -#include - -BEGIN_HADRONS_NAMESPACE - -/****************************************************************************** - * MQuark * - ******************************************************************************/ -class MQuarkPar: Serializable -{ -public: - GRID_SERIALIZABLE_CLASS_MEMBERS(MQuarkPar, - std::string, source, - std::string, solver); -}; - -class MQuark: public Module -{ -public: - // constructor - MQuark(const std::string name); - // destructor - virtual ~MQuark(void) = default; - // dependencies/products - virtual std::vector getInput(void); - virtual std::vector getOutput(void); - // setup - virtual void setup(void); - // execution - virtual void execute(void); -private: - unsigned int Ls_; - SolverFn *solver_{nullptr}; -}; - -MODULE_REGISTER(MQuark); - -END_HADRONS_NAMESPACE - -#endif // Hadrons_MQuark_hpp_ diff --git a/extras/Hadrons/Modules/MSolver/RBPrecCG.cc b/extras/Hadrons/Modules/MSolver/RBPrecCG.cc deleted file mode 100644 index 9cdd05aa..00000000 --- a/extras/Hadrons/Modules/MSolver/RBPrecCG.cc +++ /dev/null @@ -1,84 +0,0 @@ -/******************************************************************************* -Grid physics library, www.github.com/paboyle/Grid - -Source file: programs/Hadrons/RBPrecCG.cc - -Copyright (C) 2016 - -Author: Antonin Portelli - -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. -*******************************************************************************/ - -#include - -using namespace Grid; -using namespace QCD; -using namespace Hadrons; -using namespace MSolver; - -/****************************************************************************** -* RBPrecCG implementation * -******************************************************************************/ -// constructor ///////////////////////////////////////////////////////////////// -RBPrecCG::RBPrecCG(const std::string name) -: Module(name) -{} - -// dependencies/products /////////////////////////////////////////////////////// -std::vector RBPrecCG::getInput(void) -{ - std::vector in = {par().action}; - - return in; -} - -std::vector RBPrecCG::getOutput(void) -{ - std::vector out = {getName()}; - - return out; -} - -// setup /////////////////////////////////////////////////////////////////////// -void RBPrecCG::setup(void) -{ - auto Ls = env().getObjectLs(par().action); - - env().registerObject(getName(), 0, Ls); - env().addOwnership(getName(), par().action); -} - -// execution /////////////////////////////////////////////////////////////////// -void RBPrecCG::execute(void) -{ - auto &mat = *(env().getObject(par().action)); - auto solver = [&mat, this](LatticeFermion &sol, - const LatticeFermion &source) - { - ConjugateGradient cg(par().residual, 10000); - SchurRedBlackDiagMooeeSolve schurSolver(cg); - - schurSolver(mat, source, sol); - }; - - LOG(Message) << "setting up Schur red-black preconditioned CG for" - << " action '" << par().action << "' with residual " - << par().residual << std::endl; - env().setObject(getName(), new SolverFn(solver)); -} diff --git a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp index 14d8c169..57ebb816 100644 --- a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp +++ b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp @@ -1,7 +1,7 @@ /******************************************************************************* Grid physics library, www.github.com/paboyle/Grid -Source file: programs/Hadrons/RBPrecCG.hpp +Source file: programs/Hadrons/TRBPrecCG.hpp Copyright (C) 2016 @@ -47,13 +47,16 @@ public: double , residual); }; -class RBPrecCG: public Module +template +class TRBPrecCG: public Module { +public: + TYPE_ALIASES(FImpl,); public: // constructor - RBPrecCG(const std::string name); + TRBPrecCG(const std::string name); // destructor - virtual ~RBPrecCG(void) = default; + virtual ~TRBPrecCG(void) = default; // dependencies/products virtual std::vector getInput(void); virtual std::vector getOutput(void); @@ -63,6 +66,63 @@ public: virtual void execute(void); }; +/****************************************************************************** + * TRBPrecCG template implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TRBPrecCG::TRBPrecCG(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TRBPrecCG::getInput(void) +{ + std::vector in = {par().action}; + + return in; +} + +template +std::vector TRBPrecCG::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TRBPrecCG::setup(void) +{ + auto Ls = env().getObjectLs(par().action); + + env().registerObject(getName(), 0, Ls); + env().addOwnership(getName(), par().action); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TRBPrecCG::execute(void) +{ + auto &mat = *(env().template getObject(par().action)); + auto solver = [&mat, this](FermionField &sol, const FermionField &source) + { + ConjugateGradient cg(par().residual, 10000); + SchurRedBlackDiagMooeeSolve schurSolver(cg); + + schurSolver(mat, source, sol); + }; + + LOG(Message) << "setting up Schur red-black preconditioned CG for" + << " action '" << par().action << "' with residual " + << par().residual << std::endl; + env().setObject(getName(), new SolverFn(solver)); +} + +typedef TRBPrecCG RBPrecCG; + END_MODULE_NAMESPACE MODULE_REGISTER_NS(RBPrecCG, MSolver); diff --git a/extras/Hadrons/Modules/MSource/Point.cc b/extras/Hadrons/Modules/MSource/Point.cc deleted file mode 100644 index f22f41d9..00000000 --- a/extras/Hadrons/Modules/MSource/Point.cc +++ /dev/null @@ -1,75 +0,0 @@ -/******************************************************************************* -Grid physics library, www.github.com/paboyle/Grid - -Source file: programs/Hadrons/Point.cc - -Copyright (C) 2016 - -Author: Antonin Portelli - -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. -*******************************************************************************/ - -#include - -using namespace Grid; -using namespace Hadrons; -using namespace MSource; - -/****************************************************************************** -* Point implementation * -******************************************************************************/ -// constructor ///////////////////////////////////////////////////////////////// -Point::Point(const std::string name) -: Module(name) -{} - -// dependencies/products /////////////////////////////////////////////////////// -std::vector Point::getInput(void) -{ - std::vector in; - - return in; -} - -std::vector Point::getOutput(void) -{ - std::vector out = {getName()}; - - return out; -} - -// setup /////////////////////////////////////////////////////////////////////// -void Point::setup(void) -{ - env().registerLattice(getName()); -} - -// execution /////////////////////////////////////////////////////////////////// -void Point::execute(void) -{ - std::vector position = strToVec(par().position); - SpinColourMatrix id; - - LOG(Message) << "Creating point source at position [" << par().position - << "]" << std::endl; - PropagatorField &src = *env().createLattice(getName()); - id = 1.; - src = zero; - pokeSite(id, src, position); -} diff --git a/extras/Hadrons/Modules/MSource/Point.hpp b/extras/Hadrons/Modules/MSource/Point.hpp index 1245ed9a..ad0a37b5 100644 --- a/extras/Hadrons/Modules/MSource/Point.hpp +++ b/extras/Hadrons/Modules/MSource/Point.hpp @@ -1,7 +1,7 @@ /******************************************************************************* Grid physics library, www.github.com/paboyle/Grid -Source file: programs/Hadrons/Point.hpp +Source file: programs/Hadrons/TPoint.hpp Copyright (C) 2016 @@ -46,7 +46,7 @@ BEGIN_HADRONS_NAMESPACE */ /****************************************************************************** - * Point * + * TPoint * ******************************************************************************/ BEGIN_MODULE_NAMESPACE(MSource) @@ -57,13 +57,16 @@ public: std::string, position); }; -class Point: public Module +template +class TPoint: public Module { +public: + TYPE_ALIASES(FImpl,); public: // constructor - Point(const std::string name); + TPoint(const std::string name); // destructor - virtual ~Point(void) = default; + virtual ~TPoint(void) = default; // dependency relation virtual std::vector getInput(void); virtual std::vector getOutput(void); @@ -73,6 +76,56 @@ public: virtual void execute(void); }; +/****************************************************************************** + * TPoint template implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TPoint::TPoint(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TPoint::getInput(void) +{ + std::vector in; + + return in; +} + +template +std::vector TPoint::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TPoint::setup(void) +{ + env().template registerLattice(getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TPoint::execute(void) +{ + std::vector position = strToVec(par().position); + typename SitePropagator::scalar_object id; + + LOG(Message) << "Creating point source at position [" << par().position + << "]" << std::endl; + PropagatorField &src = *env().template createLattice(getName()); + id = 1.; + src = zero; + pokeSite(id, src, position); +} + +typedef TPoint Point; + END_MODULE_NAMESPACE MODULE_REGISTER_NS(Point, MSource); diff --git a/extras/Hadrons/Modules/MSource/Z2.cc b/extras/Hadrons/Modules/MSource/Z2.cc deleted file mode 100644 index 90f9f274..00000000 --- a/extras/Hadrons/Modules/MSource/Z2.cc +++ /dev/null @@ -1,88 +0,0 @@ -/******************************************************************************* -Grid physics library, www.github.com/paboyle/Grid - -Source file: programs/Hadrons/Z2.cc - -Copyright (C) 2016 - -Author: Antonin Portelli - -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. -*******************************************************************************/ - -#include - -using namespace Grid; -using namespace Hadrons; -using namespace MSource; - -/****************************************************************************** -* Z2 implementation * -******************************************************************************/ -// constructor ///////////////////////////////////////////////////////////////// -Z2::Z2(const std::string name) -: Module(name) -{} - -// dependencies/products /////////////////////////////////////////////////////// -std::vector Z2::getInput(void) -{ - std::vector in; - - return in; -} - -std::vector Z2::getOutput(void) -{ - std::vector out = {getName()}; - - return out; -} - -// setup /////////////////////////////////////////////////////////////////////// -void Z2::setup(void) -{ - env().registerLattice(getName()); -} - -// execution /////////////////////////////////////////////////////////////////// -void Z2::execute(void) -{ - Lattice> t(env().getGrid()); - LatticeComplex eta(env().getGrid()); - LatticeFermion phi(env().getGrid()); - Complex shift(1., 1.); - - if (par().tA == par().tB) - { - LOG(Message) << "Generating Z_2 wall source at t= " << par().tA - << std::endl; - } - else - { - LOG(Message) << "Generating Z_2 band for " << par().tA << " <= t <= " - << par().tB << std::endl; - } - PropagatorField &src = *env().createLattice(getName()); - LatticeCoordinate(t, Tp); - bernoulli(*env().get4dRng(), eta); - eta = (2.*eta - shift)*(1./::sqrt(2.)); - eta = where((t >= par().tA) and (t <= par().tB), eta, 0.*eta); - src = 1.; - src = src*eta; -} diff --git a/extras/Hadrons/Modules/MSource/Z2.hpp b/extras/Hadrons/Modules/MSource/Z2.hpp index ac569ced..21d0b2bd 100644 --- a/extras/Hadrons/Modules/MSource/Z2.hpp +++ b/extras/Hadrons/Modules/MSource/Z2.hpp @@ -1,7 +1,7 @@ /******************************************************************************* Grid physics library, www.github.com/paboyle/Grid -Source file: programs/Hadrons/Z2.hpp +Source file: programs/Hadrons/TZ2.hpp Copyright (C) 2016 @@ -47,7 +47,7 @@ BEGIN_HADRONS_NAMESPACE */ /****************************************************************************** - * Z2 * + * Z2 stochastic source * ******************************************************************************/ BEGIN_MODULE_NAMESPACE(MSource) @@ -59,13 +59,16 @@ public: unsigned int, tB); }; -class Z2: public Module +template +class TZ2: public Module { +public: + TYPE_ALIASES(FImpl,); public: // constructor - Z2(const std::string name); + TZ2(const std::string name); // destructor - virtual ~Z2(void) = default; + virtual ~TZ2(void) = default; // dependency relation virtual std::vector getInput(void); virtual std::vector getOutput(void); @@ -75,6 +78,68 @@ public: virtual void execute(void); }; +/****************************************************************************** + * TZ2 template implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TZ2::TZ2(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TZ2::getInput(void) +{ + std::vector in; + + return in; +} + +template +std::vector TZ2::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TZ2::setup(void) +{ + env().template registerLattice(getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TZ2::execute(void) +{ + Lattice> t(env().getGrid()); + LatticeComplex eta(env().getGrid()); + Complex shift(1., 1.); + + if (par().tA == par().tB) + { + LOG(Message) << "Generating Z_2 wall source at t= " << par().tA + << std::endl; + } + else + { + LOG(Message) << "Generating Z_2 band for " << par().tA << " <= t <= " + << par().tB << std::endl; + } + PropagatorField &src = *env().template createLattice(getName()); + LatticeCoordinate(t, Tp); + bernoulli(*env().get4dRng(), eta); + eta = (2.*eta - shift)*(1./::sqrt(2.)); + eta = where((t >= par().tA) and (t <= par().tB), eta, 0.*eta); + src = 1.; + src = src*eta; +} + +typedef TZ2 Z2; + END_MODULE_NAMESPACE MODULE_REGISTER_NS(Z2, MSource); diff --git a/extras/Hadrons/Modules/MQuark.cc b/extras/Hadrons/Modules/Quark.hpp similarity index 58% rename from extras/Hadrons/Modules/MQuark.cc rename to extras/Hadrons/Modules/Quark.hpp index 07f25e1c..59274546 100644 --- a/extras/Hadrons/Modules/MQuark.cc +++ b/extras/Hadrons/Modules/Quark.hpp @@ -1,7 +1,7 @@ /******************************************************************************* Grid physics library, www.github.com/paboyle/Grid -Source file: programs/Hadrons/MQuark.cc +Source file: programs/Hadrons/TQuark.hpp Copyright (C) 2015 @@ -25,29 +25,68 @@ See the full license in the file "LICENSE" in the top level distribution directory. *******************************************************************************/ -#include +#ifndef Hadrons_Quark_hpp_ +#define Hadrons_Quark_hpp_ -using namespace Grid; -using namespace QCD; -using namespace Hadrons; +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE /****************************************************************************** - * MQuark implementation * + * TQuark * + ******************************************************************************/ +class QuarkPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(QuarkPar, + std::string, source, + std::string, solver); +}; + +template +class TQuark: public Module +{ +public: + TYPE_ALIASES(FImpl,); +public: + // constructor + TQuark(const std::string name); + // destructor + virtual ~TQuark(void) = default; + // dependencies/products + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +private: + unsigned int Ls_; + SolverFn *solver_{nullptr}; +}; + +/****************************************************************************** + * TQuark implementation * ******************************************************************************/ // constructor ///////////////////////////////////////////////////////////////// -MQuark::MQuark(const std::string name) +template +TQuark::TQuark(const std::string name) : Module(name) {} // dependencies/products /////////////////////////////////////////////////////// -std::vector MQuark::getInput(void) +template +std::vector TQuark::getInput(void) { std::vector in = {par().source, par().solver}; return in; } -std::vector MQuark::getOutput(void) +template +std::vector TQuark::getOutput(void) { std::vector out = {getName(), getName() + "_5d"}; @@ -55,40 +94,42 @@ std::vector MQuark::getOutput(void) } // setup /////////////////////////////////////////////////////////////////////// -void MQuark::setup(void) +template +void TQuark::setup(void) { Ls_ = env().getObjectLs(par().solver); - env().registerLattice(getName()); + env().template registerLattice(getName()); if (Ls_ > 1) { - env().registerLattice(getName() + "_5d", Ls_); + env().template registerLattice(getName() + "_5d", Ls_); } } // execution /////////////////////////////////////////////////////////////////// -void MQuark::execute(void) +template +void TQuark::execute(void) { - LatticeFermion source(env().getGrid(Ls_)), sol(env().getGrid(Ls_)), - tmp(env().getGrid()); - std::string propName = (Ls_ == 1) ? getName() : (getName() + "_5d"); - LOG(Message) << "Computing quark propagator '" << getName() << "'" << std::endl; - LatticePropagator &prop = *env().createLattice(propName); - LatticePropagator &fullSrc = *env().getObject(par().source); - SolverFn &solver = *env().getObject(par().solver); + + FermionField source(env().getGrid(Ls_)), sol(env().getGrid(Ls_)), + tmp(env().getGrid()); + std::string propName = (Ls_ == 1) ? getName() : (getName() + "_5d"); + PropagatorField &prop = *env().template createLattice(propName); + PropagatorField &fullSrc = *env().template getObject(par().source); + SolverFn &solver = *env().template getObject(par().solver); if (Ls_ > 1) { - env().createLattice(getName()); + env().template createLattice(getName()); } - + LOG(Message) << "Inverting using solver '" << par().solver << "' on source '" << par().source << "'" << std::endl; for (unsigned int s = 0; s < Ns; ++s) for (unsigned int c = 0; c < Nc; ++c) { LOG(Message) << "Inversion for spin= " << s << ", color= " << c - << std::endl; + << std::endl; // source conversion for 4D sources if (!env().isObject5d(par().source)) { @@ -124,7 +165,8 @@ void MQuark::execute(void) // create 4D propagators from 5D one if necessary if (Ls_ > 1) { - LatticePropagator &p4d = *env().getObject(getName()); + PropagatorField &p4d = + *env().template getObject(getName()); axpby_ssp_pminus(sol, 0., sol, 1., sol, 0, 0); axpby_ssp_pplus(sol, 0., sol, 1., sol, 0, Ls_-1); @@ -133,3 +175,11 @@ void MQuark::execute(void) } } } + +typedef TQuark Quark; + +MODULE_REGISTER(Quark); + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Quark_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index a460ae86..aacbbf14 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -1,14 +1,7 @@ modules_cc =\ - Modules/MAction/DWF.cc \ - Modules/MAction/Wilson.cc \ - Modules/MContraction/Meson.cc \ Modules/MGauge/Load.cc \ Modules/MGauge/Random.cc \ - Modules/MGauge/Unit.cc \ - Modules/MQuark.cc \ - Modules/MSolver/RBPrecCG.cc \ - Modules/MSource/Point.cc \ - Modules/MSource/Z2.cc + Modules/MGauge/Unit.cc modules_hpp =\ Modules/MAction/DWF.hpp \ @@ -17,8 +10,8 @@ modules_hpp =\ Modules/MGauge/Load.hpp \ Modules/MGauge/Random.hpp \ Modules/MGauge/Unit.hpp \ - Modules/MQuark.hpp \ Modules/MSolver/RBPrecCG.hpp \ Modules/MSource/Point.hpp \ - Modules/MSource/Z2.hpp + Modules/MSource/Z2.hpp \ + Modules/Quark.hpp diff --git a/tests/hadrons/Test_hadrons_mesons.cc b/tests/hadrons/Test_hadrons_mesons.cc index b9561694..75d752b2 100644 --- a/tests/hadrons/Test_hadrons_mesons.cc +++ b/tests/hadrons/Test_hadrons_mesons.cc @@ -82,10 +82,10 @@ int main(int argc, char *argv[]) solverPar); // propagators - MQuark::Par quarkPar; + Quark::Par quarkPar; quarkPar.solver = "CG_" + flavour[i]; quarkPar.source = (flavour[i][0] == 'c') ? "z2" : "pt"; - application.createModule("Q_" + flavour[i], quarkPar); + application.createModule("Q_" + flavour[i], quarkPar); } for (unsigned int i = 0; i < flavour.size(); ++i) for (unsigned int j = i + 1; j < flavour.size(); ++j)