diff --git a/.gitignore b/.gitignore index db26c111..d3e1c2a8 100644 --- a/.gitignore +++ b/.gitignore @@ -9,6 +9,7 @@ ################ *~ *# +*.sublime-* # Precompiled Headers # ####################### diff --git a/Makefile.am b/Makefile.am index 818f0983..d3401c48 100644 --- a/Makefile.am +++ b/Makefile.am @@ -1,5 +1,5 @@ # additional include paths necessary to compile the C++ library -SUBDIRS = lib benchmarks tests +SUBDIRS = lib benchmarks tests extras include $(top_srcdir)/doxygen.inc diff --git a/benchmarks/Benchmark_mooee.cc b/benchmarks/Benchmark_mooee.cc index dfaea627..1e51c9d2 100644 --- a/benchmarks/Benchmark_mooee.cc +++ b/benchmarks/Benchmark_mooee.cc @@ -113,6 +113,36 @@ int main (int argc, char ** argv) std::cout<Barrier(); \ + t0=usecond(); \ + for(int i=0;iBarrier(); \ + zDw.CayleyReport(); \ + std::cout<Barrier(); \ + t0=usecond(); \ + for(int i=0;iBarrier(); \ + Dw.CayleyReport(); \ + std::cout< gamma(Ls,std::complex(1.0,0.0)); + ZMobiusFermionVec5dR zDw(Umu,*sFGrid,*sFrbGrid,*sUGrid,*sUrbGrid,mass,M5,gamma,b,c); + std::cout<Barrier(); @@ -173,10 +209,13 @@ int main (int argc, char ** argv) BENCH_DW_MEO(Dhop ,src,result); BENCH_DW_MEO(DhopEO ,src_o,r_e); - BENCH_DW(Meooe ,src_o,r_e); + BENCH_DW_SSC(Meooe ,src_o,r_e); BENCH_DW(Mooee ,src_o,r_o); BENCH_DW(MooeeInv,src_o,r_o); + BENCH_ZDW(Mooee ,src_o,r_o); + BENCH_ZDW(MooeeInv,src_o,r_o); + } Grid_finalize(); diff --git a/configure.ac b/configure.ac index 5c0805aa..27ff79ac 100644 --- a/configure.ac +++ b/configure.ac @@ -102,6 +102,13 @@ case ${ac_MKL} in AC_DEFINE([USE_MKL], [1], [Define to 1 if you use the Intel MKL]);; esac +############### HDF5 +AC_ARG_WITH([hdf5], + [AS_HELP_STRING([--with-hdf5=prefix], + [try this for a non-standard install prefix of the HDF5 library])], + [AM_CXXFLAGS="-I$with_hdf5/include $AM_CXXFLAGS"] + [AM_LDFLAGS="-L$with_hdf5/lib $AM_LDFLAGS"]) + ############### first-touch AC_ARG_ENABLE([numa], [AC_HELP_STRING([--enable-numa=yes|no|prefix], [enable first touch numa opt])], @@ -148,6 +155,12 @@ AC_SEARCH_LIBS([fftw_execute], [fftw3], [AC_DEFINE([HAVE_FFTW], [1], [Define to 1 if you have the `FFTW' library])] [have_fftw=true]) +AC_SEARCH_LIBS([H5Fopen], [hdf5_cpp], + [AC_DEFINE([HAVE_HDF5], [1], [Define to 1 if you have the `HDF5' library])] + [have_hdf5=true] + [LIBS="${LIBS} -lhdf5"], [], [-lhdf5]) +AM_CONDITIONAL(BUILD_HDF5, [ test "${have_hdf5}X" == "trueX" ]) + CXXFLAGS=$CXXFLAGS_CPY LDFLAGS=$LDFLAGS_CPY @@ -393,10 +406,13 @@ AC_CONFIG_FILES(tests/IO/Makefile) AC_CONFIG_FILES(tests/core/Makefile) AC_CONFIG_FILES(tests/debug/Makefile) AC_CONFIG_FILES(tests/forces/Makefile) +AC_CONFIG_FILES(tests/hadrons/Makefile) AC_CONFIG_FILES(tests/hmc/Makefile) AC_CONFIG_FILES(tests/solver/Makefile) AC_CONFIG_FILES(tests/qdpxx/Makefile) AC_CONFIG_FILES(benchmarks/Makefile) +AC_CONFIG_FILES(extras/Makefile) +AC_CONFIG_FILES(extras/Hadrons/Makefile) AC_OUTPUT git_commit=`cd $srcdir && ./scripts/configure.commit` @@ -425,6 +441,7 @@ GMP : `if test "x$have_gmp" = xtrue; then echo yes; else LAPACK : ${ac_LAPACK} FFTW : `if test "x$have_fftw" = xtrue; then echo yes; else echo no; fi` LIME (ILDG support) : `if test "x$have_lime" = xtrue; then echo yes; else echo no; fi` +HDF5 : `if test "x$have_hdf5" = xtrue; then echo yes; else echo no; fi` build DOXYGEN documentation : `if test "$DX_FLAG_doc" = '1'; then echo yes; else echo no; fi` ----- BUILD FLAGS ------------------------------------- CXXFLAGS: diff --git a/extras/Hadrons/Application.cc b/extras/Hadrons/Application.cc new file mode 100644 index 00000000..62674f30 --- /dev/null +++ b/extras/Hadrons/Application.cc @@ -0,0 +1,317 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Application.cc + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#include +#include + +using namespace Grid; +using namespace QCD; +using namespace Hadrons; + +#define BIG_SEP "===============" +#define SEP "---------------" + +/****************************************************************************** + * Application implementation * + ******************************************************************************/ +// constructors //////////////////////////////////////////////////////////////// +Application::Application(void) +{ + LOG(Message) << "Modules available:" << std::endl; + auto list = ModuleFactory::getInstance().getBuilderList(); + for (auto &m: list) + { + LOG(Message) << " " << m << std::endl; + } + auto dim = GridDefaultLatt(), mpi = GridDefaultMpi(), loc(dim); + locVol_ = 1; + for (unsigned int d = 0; d < dim.size(); ++d) + { + loc[d] /= mpi[d]; + locVol_ *= loc[d]; + } + LOG(Message) << "Global lattice: " << dim << std::endl; + LOG(Message) << "MPI partition : " << mpi << std::endl; + LOG(Message) << "Local lattice : " << loc << std::endl; +} + +Application::Application(const Application::GlobalPar &par) +: Application() +{ + setPar(par); +} + +Application::Application(const std::string parameterFileName) +: Application() +{ + parameterFileName_ = parameterFileName; +} + +// environment shortcut //////////////////////////////////////////////////////// +Environment & Application::env(void) const +{ + return Environment::getInstance(); +} + +// access ////////////////////////////////////////////////////////////////////// +void Application::setPar(const Application::GlobalPar &par) +{ + par_ = par; + env().setSeed(strToVec(par_.seed)); +} + +const Application::GlobalPar & Application::getPar(void) +{ + return par_; +} + +// execute ///////////////////////////////////////////////////////////////////// +void Application::run(void) +{ + if (!parameterFileName_.empty() and (env().getNModule() == 0)) + { + parseParameterFile(parameterFileName_); + } + if (!scheduled_) + { + schedule(); + } + printSchedule(); + configLoop(); +} + +// parse parameter file //////////////////////////////////////////////////////// +class ObjectId: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(ObjectId, + std::string, name, + std::string, type); +}; + +void Application::parseParameterFile(const std::string parameterFileName) +{ + XmlReader reader(parameterFileName); + GlobalPar par; + ObjectId id; + + LOG(Message) << "Building application from '" << parameterFileName << "'..." << std::endl; + read(reader, "parameters", par); + setPar(par); + push(reader, "modules"); + push(reader, "module"); + do + { + read(reader, "id", id); + env().createModule(id.name, id.type, reader); + } while (reader.nextElement("module")); + pop(reader); + pop(reader); +} + +void Application::saveParameterFile(const std::string parameterFileName) +{ + XmlWriter writer(parameterFileName); + ObjectId id; + const unsigned int nMod = env().getNModule(); + + LOG(Message) << "Saving application to '" << parameterFileName << "'..." << std::endl; + write(writer, "parameters", getPar()); + push(writer, "modules"); + for (unsigned int i = 0; i < nMod; ++i) + { + push(writer, "module"); + id.name = env().getModuleName(i); + id.type = env().getModule(i)->getRegisteredName(); + write(writer, "id", id); + env().getModule(i)->saveParameters(writer, "options"); + pop(writer); + } + pop(writer); + pop(writer); +} + +// schedule computation //////////////////////////////////////////////////////// +#define MEM_MSG(size)\ +sizeString((size)*locVol_) << " (" << sizeString(size) << "/site)" + +#define DEFINE_MEMPEAK \ +auto memPeak = [this](const std::vector &program)\ +{\ + unsigned int memPeak;\ + bool msg;\ + \ + msg = HadronsLogMessage.isActive();\ + HadronsLogMessage.Active(false);\ + env().dryRun(true);\ + memPeak = env().executeProgram(program);\ + env().dryRun(false);\ + env().freeAll();\ + HadronsLogMessage.Active(true);\ + \ + return memPeak;\ +} + +void Application::schedule(void) +{ + DEFINE_MEMPEAK; + + // build module dependency graph + LOG(Message) << "Building module graph..." << std::endl; + auto graph = env().makeModuleGraph(); + auto con = graph.getConnectedComponents(); + + // constrained topological sort using a genetic algorithm + LOG(Message) << "Scheduling computation..." << std::endl; + LOG(Message) << " #module= " << graph.size() << std::endl; + LOG(Message) << " population size= " << par_.genetic.popSize << std::endl; + LOG(Message) << " max. generation= " << par_.genetic.maxGen << std::endl; + LOG(Message) << " max. cst. generation= " << par_.genetic.maxCstGen << std::endl; + LOG(Message) << " mutation rate= " << par_.genetic.mutationRate << std::endl; + + unsigned int k = 0, gen, prevPeak, nCstPeak = 0; + std::random_device rd; + GeneticScheduler::Parameters par; + + par.popSize = par_.genetic.popSize; + par.mutationRate = par_.genetic.mutationRate; + par.seed = rd(); + memPeak_ = 0; + CartesianCommunicator::BroadcastWorld(0, &(par.seed), sizeof(par.seed)); + for (unsigned int i = 0; i < con.size(); ++i) + { + GeneticScheduler scheduler(con[i], memPeak, par); + + gen = 0; + do + { + LOG(Debug) << "Generation " << gen << ":" << std::endl; + scheduler.nextGeneration(); + if (gen != 0) + { + if (prevPeak == scheduler.getMinValue()) + { + nCstPeak++; + } + else + { + nCstPeak = 0; + } + } + + prevPeak = scheduler.getMinValue(); + if (gen % 10 == 0) + { + LOG(Iterative) << "Generation " << gen << ": " + << MEM_MSG(scheduler.getMinValue()) << std::endl; + } + + gen++; + } while ((gen < par_.genetic.maxGen) + and (nCstPeak < par_.genetic.maxCstGen)); + auto &t = scheduler.getMinSchedule(); + if (scheduler.getMinValue() > memPeak_) + { + memPeak_ = scheduler.getMinValue(); + } + for (unsigned int j = 0; j < t.size(); ++j) + { + program_.push_back(t[j]); + } + } + scheduled_ = true; +} + +void Application::saveSchedule(const std::string filename) +{ + TextWriter writer(filename); + std::vector program; + + if (!scheduled_) + { + HADRON_ERROR("Computation not scheduled"); + } + LOG(Message) << "Saving current schedule to '" << filename << "'..." + << std::endl; + for (auto address: program_) + { + program.push_back(env().getModuleName(address)); + } + write(writer, "schedule", program); +} + +void Application::loadSchedule(const std::string filename) +{ + DEFINE_MEMPEAK; + + TextReader reader(filename); + std::vector program; + + LOG(Message) << "Loading schedule from '" << filename << "'..." + << std::endl; + read(reader, "schedule", program); + program_.clear(); + for (auto &name: program) + { + program_.push_back(env().getModuleAddress(name)); + } + scheduled_ = true; + memPeak_ = memPeak(program_); +} + +void Application::printSchedule(void) +{ + if (!scheduled_) + { + HADRON_ERROR("Computation not scheduled"); + } + LOG(Message) << "Schedule (memory peak: " << MEM_MSG(memPeak_) << "):" + << std::endl; + for (unsigned int i = 0; i < program_.size(); ++i) + { + LOG(Message) << std::setw(4) << i + 1 << ": " + << env().getModuleName(program_[i]) << std::endl; + } +} + +// loop on configurations ////////////////////////////////////////////////////// +void Application::configLoop(void) +{ + auto range = par_.trajCounter; + + for (unsigned int t = range.start; t < range.end; t += range.step) + { + LOG(Message) << BIG_SEP << " Starting measurement for trajectory " << t + << " " << BIG_SEP << std::endl; + env().setTrajectory(t); + env().executeProgram(program_); + } + LOG(Message) << BIG_SEP << " End of measurement " << BIG_SEP << std::endl; + env().freeAll(); +} diff --git a/extras/Hadrons/Application.hpp b/extras/Hadrons/Application.hpp new file mode 100644 index 00000000..fce9b6eb --- /dev/null +++ b/extras/Hadrons/Application.hpp @@ -0,0 +1,132 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Application.hpp + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_Application_hpp_ +#define Hadrons_Application_hpp_ + +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Main program manager * + ******************************************************************************/ +class Application +{ +public: + class TrajRange: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(TrajRange, + unsigned int, start, + unsigned int, end, + unsigned int, step); + }; + class GeneticPar: Serializable + { + public: + GeneticPar(void): + popSize{20}, maxGen{1000}, maxCstGen{100}, mutationRate{.1} {}; + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(GeneticPar, + unsigned int, popSize, + unsigned int, maxGen, + unsigned int, maxCstGen, + double , mutationRate); + }; + class GlobalPar: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(GlobalPar, + TrajRange, trajCounter, + GeneticPar, genetic, + std::string, seed); + }; +public: + // constructors + Application(void); + Application(const GlobalPar &par); + Application(const std::string parameterFileName); + // destructor + virtual ~Application(void) = default; + // access + void setPar(const GlobalPar &par); + const GlobalPar & getPar(void); + // module creation + template + void createModule(const std::string name); + template + void createModule(const std::string name, const typename M::Par &par); + // execute + void run(void); + // XML parameter file I/O + void parseParameterFile(const std::string parameterFileName); + void saveParameterFile(const std::string parameterFileName); + // schedule computation + void schedule(void); + void saveSchedule(const std::string filename); + void loadSchedule(const std::string filename); + void printSchedule(void); + // loop on configurations + void configLoop(void); +private: + // environment shortcut + Environment & env(void) const; +private: + long unsigned int locVol_; + std::string parameterFileName_{""}; + GlobalPar par_; + std::vector program_; + Environment::Size memPeak_; + bool scheduled_{false}; +}; + +/****************************************************************************** + * Application template implementation * + ******************************************************************************/ +// module creation ///////////////////////////////////////////////////////////// +template +void Application::createModule(const std::string name) +{ + env().createModule(name); +} + +template +void Application::createModule(const std::string name, + const typename M::Par &par) +{ + env().createModule(name, par); +} + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Application_hpp_ diff --git a/extras/Hadrons/Environment.cc b/extras/Hadrons/Environment.cc new file mode 100644 index 00000000..37f2a3d7 --- /dev/null +++ b/extras/Hadrons/Environment.cc @@ -0,0 +1,743 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Environment.cc + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#include +#include +#include + +using namespace Grid; +using namespace QCD; +using namespace Hadrons; + +/****************************************************************************** + * Environment implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +Environment::Environment(void) +{ + nd_ = GridDefaultLatt().size(); + grid4d_.reset(SpaceTimeGrid::makeFourDimGrid( + GridDefaultLatt(), GridDefaultSimd(nd_, vComplex::Nsimd()), + GridDefaultMpi())); + gridRb4d_.reset(SpaceTimeGrid::makeFourDimRedBlackGrid(grid4d_.get())); + auto loc = getGrid()->LocalDimensions(); + locVol_ = 1; + for (unsigned int d = 0; d < loc.size(); ++d) + { + locVol_ *= loc[d]; + } + rng4d_.reset(new GridParallelRNG(grid4d_.get())); +} + +// dry run ///////////////////////////////////////////////////////////////////// +void Environment::dryRun(const bool isDry) +{ + dryRun_ = isDry; +} + +bool Environment::isDryRun(void) const +{ + return dryRun_; +} + +// trajectory number /////////////////////////////////////////////////////////// +void Environment::setTrajectory(const unsigned int traj) +{ + traj_ = traj; +} + +unsigned int Environment::getTrajectory(void) const +{ + return traj_; +} + +// grids /////////////////////////////////////////////////////////////////////// +void Environment::createGrid(const unsigned int Ls) +{ + if (grid5d_.find(Ls) == grid5d_.end()) + { + auto g = getGrid(); + + grid5d_[Ls].reset(SpaceTimeGrid::makeFiveDimGrid(Ls, g)); + gridRb5d_[Ls].reset(SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, g)); + } +} + +GridCartesian * Environment::getGrid(const unsigned int Ls) const +{ + try + { + if (Ls == 1) + { + return grid4d_.get(); + } + else + { + return grid5d_.at(Ls).get(); + } + } + catch(std::out_of_range &) + { + HADRON_ERROR("no grid with Ls= " << Ls); + } +} + +GridRedBlackCartesian * Environment::getRbGrid(const unsigned int Ls) const +{ + try + { + if (Ls == 1) + { + return gridRb4d_.get(); + } + else + { + return gridRb5d_.at(Ls).get(); + } + } + catch(std::out_of_range &) + { + HADRON_ERROR("no red-black 5D grid with Ls= " << Ls); + } +} + +unsigned int Environment::getNd(void) const +{ + return nd_; +} + +// random number generator ///////////////////////////////////////////////////// +void Environment::setSeed(const std::vector &seed) +{ + rng4d_->SeedFixedIntegers(seed); +} + +GridParallelRNG * Environment::get4dRng(void) const +{ + return rng4d_.get(); +} + +// module management /////////////////////////////////////////////////////////// +void Environment::pushModule(Environment::ModPt &pt) +{ + std::string name = pt->getName(); + + if (!hasModule(name)) + { + std::vector inputAddress; + unsigned int address; + ModuleInfo m; + + m.data = std::move(pt); + m.type = typeIdPt(*m.data.get()); + m.name = name; + auto input = m.data->getInput(); + for (auto &in: input) + { + if (!hasObject(in)) + { + addObject(in , -1); + } + m.input.push_back(objectAddress_[in]); + } + auto output = m.data->getOutput(); + module_.push_back(std::move(m)); + address = static_cast(module_.size() - 1); + moduleAddress_[name] = address; + for (auto &out: output) + { + if (!hasObject(out)) + { + addObject(out, address); + } + else + { + if (object_[objectAddress_[out]].module < 0) + { + object_[objectAddress_[out]].module = address; + } + else + { + HADRON_ERROR("object '" + out + + "' is already produced by module '" + + module_[object_[getObjectAddress(out)].module].name + + "' (while pushing module '" + name + "')"); + } + } + } + } + else + { + HADRON_ERROR("module '" + name + "' already exists"); + } +} + +unsigned int Environment::getNModule(void) const +{ + return module_.size(); +} + +void Environment::createModule(const std::string name, const std::string type, + XmlReader &reader) +{ + auto &factory = ModuleFactory::getInstance(); + auto pt = factory.create(type, name); + + pt->parseParameters(reader, "options"); + pushModule(pt); +} + +ModuleBase * Environment::getModule(const unsigned int address) const +{ + if (hasModule(address)) + { + return module_[address].data.get(); + } + else + { + HADRON_ERROR("no module with address " + std::to_string(address)); + } +} + +ModuleBase * Environment::getModule(const std::string name) const +{ + return getModule(getModuleAddress(name)); +} + +unsigned int Environment::getModuleAddress(const std::string name) const +{ + if (hasModule(name)) + { + return moduleAddress_.at(name); + } + else + { + HADRON_ERROR("no module with name '" + name + "'"); + } +} + +std::string Environment::getModuleName(const unsigned int address) const +{ + if (hasModule(address)) + { + return module_[address].name; + } + else + { + HADRON_ERROR("no module with address " + std::to_string(address)); + } +} + +std::string Environment::getModuleType(const unsigned int address) const +{ + if (hasModule(address)) + { + return typeName(module_[address].type); + } + else + { + HADRON_ERROR("no module with address " + std::to_string(address)); + } +} + +std::string Environment::getModuleType(const std::string name) const +{ + return getModuleType(getModuleAddress(name)); +} + +bool Environment::hasModule(const unsigned int address) const +{ + return (address < module_.size()); +} + +bool Environment::hasModule(const std::string name) const +{ + return (moduleAddress_.find(name) != moduleAddress_.end()); +} + +Graph Environment::makeModuleGraph(void) const +{ + Graph moduleGraph; + + for (unsigned int i = 0; i < module_.size(); ++i) + { + moduleGraph.addVertex(i); + for (auto &j: module_[i].input) + { + moduleGraph.addEdge(object_[j].module, i); + } + } + + return moduleGraph; +} + +#define BIG_SEP "===============" +#define SEP "---------------" +#define MEM_MSG(size)\ +sizeString((size)*locVol_) << " (" << sizeString(size) << "/site)" + +Environment::Size +Environment::executeProgram(const std::vector &p) +{ + Size memPeak = 0, sizeBefore, sizeAfter; + std::vector> freeProg; + bool continueCollect, nothingFreed; + + // build garbage collection schedule + freeProg.resize(p.size()); + for (unsigned int i = 0; i < object_.size(); ++i) + { + auto pred = [i, this](const unsigned int j) + { + auto &in = module_[j].input; + auto it = std::find(in.begin(), in.end(), i); + + return (it != in.end()) or (j == object_[i].module); + }; + auto it = std::find_if(p.rbegin(), p.rend(), pred); + if (it != p.rend()) + { + freeProg[p.rend() - it - 1].insert(i); + } + } + + // program execution + for (unsigned int i = 0; i < p.size(); ++i) + { + // execute module + if (!isDryRun()) + { + LOG(Message) << SEP << " Measurement step " << i+1 << "/" + << p.size() << " (module '" << module_[p[i]].name + << "') " << SEP << std::endl; + } + (*module_[p[i]].data)(); + sizeBefore = getTotalSize(); + // print used memory after execution + if (!isDryRun()) + { + LOG(Message) << "Allocated objects: " << MEM_MSG(sizeBefore) + << std::endl; + } + if (sizeBefore > memPeak) + { + memPeak = sizeBefore; + } + // garbage collection for step i + if (!isDryRun()) + { + LOG(Message) << "Garbage collection..." << std::endl; + } + nothingFreed = true; + do + { + continueCollect = false; + auto toFree = freeProg[i]; + for (auto &j: toFree) + { + // continue garbage collection while there are still + // objects without owners + continueCollect = continueCollect or !hasOwners(j); + if(freeObject(j)) + { + // if an object has been freed, remove it from + // the garbage collection schedule + freeProg[i].erase(j); + nothingFreed = false; + } + } + } while (continueCollect); + // any remaining objects in step i garbage collection schedule + // is scheduled for step i + 1 + if (i + 1 < p.size()) + { + for (auto &j: freeProg[i]) + { + freeProg[i + 1].insert(j); + } + } + // print used memory after garbage collection if necessary + if (!isDryRun()) + { + sizeAfter = getTotalSize(); + if (sizeBefore != sizeAfter) + { + LOG(Message) << "Allocated objects: " << MEM_MSG(sizeAfter) + << std::endl; + } + else + { + LOG(Message) << "Nothing to free" << std::endl; + } + } + } + + return memPeak; +} + +Environment::Size Environment::executeProgram(const std::vector &p) +{ + std::vector pAddress; + + for (auto &n: p) + { + pAddress.push_back(getModuleAddress(n)); + } + + return executeProgram(pAddress); +} + +// general memory management /////////////////////////////////////////////////// +void Environment::addObject(const std::string name, const int moduleAddress) +{ + if (!hasObject(name)) + { + ObjInfo info; + + info.name = name; + info.module = moduleAddress; + object_.push_back(std::move(info)); + objectAddress_[name] = static_cast(object_.size() - 1); + } + else + { + HADRON_ERROR("object '" + name + "' already exists"); + } +} + +void Environment::registerObject(const unsigned int address, + const unsigned int size, const unsigned int Ls) +{ + if (!hasRegisteredObject(address)) + { + if (hasObject(address)) + { + object_[address].size = size; + object_[address].Ls = Ls; + object_[address].isRegistered = true; + } + else + { + HADRON_ERROR("no object with address " + std::to_string(address)); + } + } + else + { + HADRON_ERROR("object with address " + std::to_string(address) + + " already registered"); + } +} + +void Environment::registerObject(const std::string name, + const unsigned int size, const unsigned int Ls) +{ + if (!hasObject(name)) + { + addObject(name); + } + registerObject(getObjectAddress(name), size, Ls); +} + +unsigned int Environment::getObjectAddress(const std::string name) const +{ + if (hasObject(name)) + { + return objectAddress_.at(name); + } + else + { + HADRON_ERROR("no object with name '" + name + "'"); + } +} + +std::string Environment::getObjectName(const unsigned int address) const +{ + if (hasObject(address)) + { + return object_[address].name; + } + else + { + HADRON_ERROR("no object with address " + std::to_string(address)); + } +} + +std::string Environment::getObjectType(const unsigned int address) const +{ + if (hasRegisteredObject(address)) + { + return typeName(object_[address].type); + } + else if (hasObject(address)) + { + HADRON_ERROR("object with address " + std::to_string(address) + + " exists but is not registered"); + } + else + { + HADRON_ERROR("no object with address " + std::to_string(address)); + } +} + +std::string Environment::getObjectType(const std::string name) const +{ + return getObjectType(getObjectAddress(name)); +} + +Environment::Size Environment::getObjectSize(const unsigned int address) const +{ + if (hasRegisteredObject(address)) + { + return object_[address].size; + } + else if (hasObject(address)) + { + HADRON_ERROR("object with address " + std::to_string(address) + + " exists but is not registered"); + } + else + { + HADRON_ERROR("no object with address " + std::to_string(address)); + } +} + +Environment::Size Environment::getObjectSize(const std::string name) const +{ + return getObjectSize(getObjectAddress(name)); +} + +unsigned int Environment::getObjectLs(const unsigned int address) const +{ + if (hasRegisteredObject(address)) + { + return object_[address].Ls; + } + else if (hasObject(address)) + { + HADRON_ERROR("object with address " + std::to_string(address) + + " exists but is not registered"); + } + else + { + HADRON_ERROR("no object with address " + std::to_string(address)); + } +} + +unsigned int Environment::getObjectLs(const std::string name) const +{ + return getObjectLs(getObjectAddress(name)); +} + +bool Environment::hasObject(const unsigned int address) const +{ + return (address < object_.size()); +} + +bool Environment::hasObject(const std::string name) const +{ + auto it = objectAddress_.find(name); + + return ((it != objectAddress_.end()) and hasObject(it->second)); +} + +bool Environment::hasRegisteredObject(const unsigned int address) const +{ + if (hasObject(address)) + { + return object_[address].isRegistered; + } + else + { + return false; + } +} + +bool Environment::hasRegisteredObject(const std::string name) const +{ + if (hasObject(name)) + { + return hasRegisteredObject(getObjectAddress(name)); + } + else + { + return false; + } +} + +bool Environment::hasCreatedObject(const unsigned int address) const +{ + if (hasObject(address)) + { + return (object_[address].data != nullptr); + } + else + { + return false; + } +} + +bool Environment::hasCreatedObject(const std::string name) const +{ + if (hasObject(name)) + { + return hasCreatedObject(getObjectAddress(name)); + } + else + { + return false; + } +} + +bool Environment::isObject5d(const unsigned int address) const +{ + return (getObjectLs(address) > 1); +} + +bool Environment::isObject5d(const std::string name) const +{ + return (getObjectLs(name) > 1); +} + +Environment::Size Environment::getTotalSize(void) const +{ + Environment::Size size = 0; + + for (auto &o: object_) + { + if (o.isRegistered) + { + size += o.size; + } + } + + return size; +} + +void Environment::addOwnership(const unsigned int owner, + const unsigned int property) +{ + if (hasObject(property)) + { + object_[property].owners.insert(owner); + } + else + { + HADRON_ERROR("no object with address " + std::to_string(property)); + } + if (hasObject(owner)) + { + object_[owner].properties.insert(property); + } + else + { + HADRON_ERROR("no object with address " + std::to_string(owner)); + } +} + +void Environment::addOwnership(const std::string owner, + const std::string property) +{ + addOwnership(getObjectAddress(owner), getObjectAddress(property)); +} + +bool Environment::hasOwners(const unsigned int address) const +{ + + if (hasObject(address)) + { + return (!object_[address].owners.empty()); + } + else + { + HADRON_ERROR("no object with address " + std::to_string(address)); + } +} + +bool Environment::hasOwners(const std::string name) const +{ + return hasOwners(getObjectAddress(name)); +} + +bool Environment::freeObject(const unsigned int address) +{ + if (!hasOwners(address)) + { + if (!isDryRun() and object_[address].isRegistered) + { + LOG(Message) << "Destroying object '" << object_[address].name + << "'" << std::endl; + } + for (auto &p: object_[address].properties) + { + object_[p].owners.erase(address); + } + object_[address].size = 0; + object_[address].Ls = 0; + object_[address].isRegistered = false; + object_[address].type = nullptr; + object_[address].owners.clear(); + object_[address].properties.clear(); + object_[address].data.reset(nullptr); + + return true; + } + else + { + return false; + } +} + +bool Environment::freeObject(const std::string name) +{ + return freeObject(getObjectAddress(name)); +} + +void Environment::freeAll(void) +{ + for (unsigned int i = 0; i < object_.size(); ++i) + { + freeObject(i); + } +} + +void Environment::printContent(void) +{ + LOG(Message) << "Modules: " << std::endl; + for (unsigned int i = 0; i < module_.size(); ++i) + { + LOG(Message) << std::setw(4) << i << ": " + << getModuleName(i) << std::endl; + } + LOG(Message) << "Objects: " << std::endl; + for (unsigned int i = 0; i < object_.size(); ++i) + { + LOG(Message) << std::setw(4) << i << ": " + << getObjectName(i) << std::endl; + } +} diff --git a/extras/Hadrons/Environment.hpp b/extras/Hadrons/Environment.hpp new file mode 100644 index 00000000..2628e5a0 --- /dev/null +++ b/extras/Hadrons/Environment.hpp @@ -0,0 +1,385 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Environment.hpp + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_Environment_hpp_ +#define Hadrons_Environment_hpp_ + +#include +#include + +#ifndef SITE_SIZE_TYPE +#define SITE_SIZE_TYPE unsigned int +#endif + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Global environment * + ******************************************************************************/ +// forward declaration of Module +class ModuleBase; + +class Object +{ +public: + Object(void) = default; + virtual ~Object(void) = default; +}; + +template +class Holder: public Object +{ +public: + Holder(void) = default; + Holder(T *pt); + virtual ~Holder(void) = default; + T & get(void) const; + T * getPt(void) const; + void reset(T *pt); +private: + std::unique_ptr objPt_{nullptr}; +}; + +class Environment +{ + SINGLETON(Environment); +public: + typedef SITE_SIZE_TYPE Size; + typedef std::unique_ptr ModPt; + typedef std::unique_ptr GridPt; + typedef std::unique_ptr GridRbPt; + typedef std::unique_ptr RngPt; + typedef std::unique_ptr LatticePt; +private: + struct ModuleInfo + { + const std::type_info *type{nullptr}; + std::string name; + ModPt data{nullptr}; + std::vector input; + }; + struct ObjInfo + { + Size size{0}; + unsigned int Ls{0}; + bool isRegistered{false}; + const std::type_info *type{nullptr}; + std::string name; + int module{-1}; + std::set owners, properties; + std::unique_ptr data{nullptr}; + }; +public: + // dry run + void dryRun(const bool isDry); + bool isDryRun(void) const; + // trajectory number + void setTrajectory(const unsigned int traj); + unsigned int getTrajectory(void) const; + // grids + void createGrid(const unsigned int Ls); + GridCartesian * getGrid(const unsigned int Ls = 1) const; + GridRedBlackCartesian * getRbGrid(const unsigned int Ls = 1) const; + unsigned int getNd(void) const; + // random number generator + void setSeed(const std::vector &seed); + GridParallelRNG * get4dRng(void) const; + // module management + void pushModule(ModPt &pt); + template + void createModule(const std::string name); + template + void createModule(const std::string name, + const typename M::Par &par); + void createModule(const std::string name, + const std::string type, + XmlReader &reader); + unsigned int getNModule(void) const; + ModuleBase * getModule(const unsigned int address) const; + ModuleBase * getModule(const std::string name) const; + template + M * getModule(const unsigned int address) const; + template + M * getModule(const std::string name) const; + unsigned int getModuleAddress(const std::string name) const; + std::string getModuleName(const unsigned int address) const; + std::string getModuleType(const unsigned int address) const; + std::string getModuleType(const std::string name) const; + bool hasModule(const unsigned int address) const; + bool hasModule(const std::string name) const; + Graph makeModuleGraph(void) const; + Size executeProgram(const std::vector &p); + Size executeProgram(const std::vector &p); + // general memory management + void addObject(const std::string name, + const int moduleAddress = -1); + void registerObject(const unsigned int address, + const unsigned int size, + const unsigned int Ls = 1); + void registerObject(const std::string name, + const unsigned int size, + const unsigned int Ls = 1); + template + unsigned int lattice4dSize(void) const; + template + void registerLattice(const unsigned int address, + const unsigned int Ls = 1); + template + void registerLattice(const std::string name, + const unsigned int Ls = 1); + template + void setObject(const unsigned int address, T *object); + template + void setObject(const std::string name, T *object); + template + T * getObject(const unsigned int address) const; + template + T * getObject(const std::string name) const; + template + T * createLattice(const unsigned int address); + template + T * createLattice(const std::string name); + unsigned int getObjectAddress(const std::string name) const; + std::string getObjectName(const unsigned int address) const; + std::string getObjectType(const unsigned int address) const; + std::string getObjectType(const std::string name) const; + Size getObjectSize(const unsigned int address) const; + Size getObjectSize(const std::string name) const; + unsigned int getObjectLs(const unsigned int address) const; + unsigned int getObjectLs(const std::string name) const; + bool hasObject(const unsigned int address) const; + bool hasObject(const std::string name) const; + bool hasRegisteredObject(const unsigned int address) const; + bool hasRegisteredObject(const std::string name) const; + bool hasCreatedObject(const unsigned int address) const; + bool hasCreatedObject(const std::string name) const; + bool isObject5d(const unsigned int address) const; + bool isObject5d(const std::string name) const; + Environment::Size getTotalSize(void) const; + void addOwnership(const unsigned int owner, + const unsigned int property); + void addOwnership(const std::string owner, + const std::string property); + bool hasOwners(const unsigned int address) const; + bool hasOwners(const std::string name) const; + bool freeObject(const unsigned int address); + bool freeObject(const std::string name); + void freeAll(void); + void printContent(void); +private: + // general + bool dryRun_{false}; + unsigned int traj_, locVol_; + // grids + GridPt grid4d_; + std::map grid5d_; + GridRbPt gridRb4d_; + std::map gridRb5d_; + unsigned int nd_; + // random number generator + RngPt rng4d_; + // module and related maps + std::vector module_; + std::map moduleAddress_; + // lattice store + std::map lattice_; + // object store + std::vector object_; + std::map objectAddress_; +}; + +/****************************************************************************** + * Holder template implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +Holder::Holder(T *pt) +: objPt_(pt) +{} + +// access ////////////////////////////////////////////////////////////////////// +template +T & Holder::get(void) const +{ + return &objPt_.get(); +} + +template +T * Holder::getPt(void) const +{ + return objPt_.get(); +} + +template +void Holder::reset(T *pt) +{ + objPt_.reset(pt); +} + +/****************************************************************************** + * Environment template implementation * + ******************************************************************************/ +// module management /////////////////////////////////////////////////////////// +template +void Environment::createModule(const std::string name) +{ + ModPt pt(new M(name)); + + pushModule(pt); +} + +template +void Environment::createModule(const std::string name, + const typename M::Par &par) +{ + ModPt pt(new M(name)); + + static_cast(pt.get())->setPar(par); + pushModule(pt); +} + +template +M * Environment::getModule(const unsigned int address) const +{ + if (auto *pt = dynamic_cast(getModule(address))) + { + return pt; + } + else + { + HADRON_ERROR("module '" + module_[address].name + + "' does not have type " + typeid(M).name() + + "(object type: " + getModuleType(address) + ")"); + } +} + +template +M * Environment::getModule(const std::string name) const +{ + return getModule(getModuleAddress(name)); +} + +template +unsigned int Environment::lattice4dSize(void) const +{ + return sizeof(typename T::vector_object)/getGrid()->Nsimd(); +} + +template +void Environment::registerLattice(const unsigned int address, + const unsigned int Ls) +{ + createGrid(Ls); + registerObject(address, Ls*lattice4dSize(), Ls); +} + +template +void Environment::registerLattice(const std::string name, const unsigned int Ls) +{ + createGrid(Ls); + registerObject(name, Ls*lattice4dSize(), Ls); +} + +template +void Environment::setObject(const unsigned int address, T *object) +{ + if (hasRegisteredObject(address)) + { + object_[address].data.reset(new Holder(object)); + object_[address].type = &typeid(T); + } + else if (hasObject(address)) + { + HADRON_ERROR("object with address " + std::to_string(address) + + " exists but is not registered"); + } + else + { + HADRON_ERROR("no object with address " + std::to_string(address)); + } +} + +template +void Environment::setObject(const std::string name, T *object) +{ + setObject(getObjectAddress(name), object); +} + +template +T * Environment::getObject(const unsigned int address) const +{ + if (hasRegisteredObject(address)) + { + if (auto h = dynamic_cast *>(object_[address].data.get())) + { + return h->getPt(); + } + else + { + HADRON_ERROR("object with address " + std::to_string(address) + + " does not have type '" + typeid(T).name() + + "' (has type '" + getObjectType(address) + "')"); + } + } + else if (hasObject(address)) + { + HADRON_ERROR("object with address " + std::to_string(address) + + " exists but is not registered"); + } + else + { + HADRON_ERROR("no object with address " + std::to_string(address)); + } +} + +template +T * Environment::getObject(const std::string name) const +{ + return getObject(getObjectAddress(name)); +} + +template +T * Environment::createLattice(const unsigned int address) +{ + GridCartesian *g = getGrid(getObjectLs(address)); + + setObject(address, new T(g)); + + return getObject(address); +} + +template +T * Environment::createLattice(const std::string name) +{ + return createLattice(getObjectAddress(name)); +} + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Environment_hpp_ diff --git a/extras/Hadrons/Factory.hpp b/extras/Hadrons/Factory.hpp new file mode 100644 index 00000000..da86acae --- /dev/null +++ b/extras/Hadrons/Factory.hpp @@ -0,0 +1,106 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Factory.hpp + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_Factory_hpp_ +#define Hadrons_Factory_hpp_ + +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * abstract factory class * + ******************************************************************************/ +template +class Factory +{ +public: + typedef std::function(const std::string)> Func; +public: + // constructor + Factory(void) = default; + // destructor + virtual ~Factory(void) = default; + // registration + void registerBuilder(const std::string type, const Func &f); + // get builder list + std::vector getBuilderList(void) const; + // factory + std::unique_ptr create(const std::string type, + const std::string name) const; +private: + std::map builder_; +}; + +/****************************************************************************** + * template implementation * + ******************************************************************************/ +// registration //////////////////////////////////////////////////////////////// +template +void Factory::registerBuilder(const std::string type, const Func &f) +{ + builder_[type] = f; +} + +// get module list ///////////////////////////////////////////////////////////// +template +std::vector Factory::getBuilderList(void) const +{ + std::vector list; + + for (auto &b: builder_) + { + list.push_back(b.first); + } + + return list; +} + +// factory ///////////////////////////////////////////////////////////////////// +template +std::unique_ptr Factory::create(const std::string type, + const std::string name) const +{ + Func func; + + try + { + func = builder_.at(type); + } + catch (std::out_of_range &) + { + HADRON_ERROR("object of type '" + type + "' unknown"); + } + + return func(name); +} + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Factory_hpp_ diff --git a/extras/Hadrons/GeneticScheduler.hpp b/extras/Hadrons/GeneticScheduler.hpp new file mode 100644 index 00000000..d0c52596 --- /dev/null +++ b/extras/Hadrons/GeneticScheduler.hpp @@ -0,0 +1,329 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/GeneticScheduler.hpp + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_GeneticScheduler_hpp_ +#define Hadrons_GeneticScheduler_hpp_ + +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Scheduler based on a genetic algorithm * + ******************************************************************************/ +template +class GeneticScheduler +{ +public: + typedef std::vector Gene; + typedef std::pair GenePair; + typedef std::function ObjFunc; + struct Parameters + { + double mutationRate; + unsigned int popSize, seed; + }; +public: + // constructor + GeneticScheduler(Graph &graph, const ObjFunc &func, + const Parameters &par); + // destructor + virtual ~GeneticScheduler(void) = default; + // access + const Gene & getMinSchedule(void); + int getMinValue(void); + // breed a new generation + void nextGeneration(void); + // heuristic benchmarks + void benchmarkCrossover(const unsigned int nIt); + // print population + friend std::ostream & operator<<(std::ostream &out, + const GeneticScheduler &s) + { + out << "["; + for (auto &p: s.population_) + { + out << p.first << ", "; + } + out << "\b\b]"; + + return out; + } +private: + // evolution steps + void initPopulation(void); + void doCrossover(void); + void doMutation(void); + // genetic operators + GenePair selectPair(void); + void crossover(Gene &c1, Gene &c2, const Gene &p1, const Gene &p2); + void mutation(Gene &m, const Gene &c); + +private: + Graph &graph_; + const ObjFunc &func_; + const Parameters par_; + std::multimap population_; + std::mt19937 gen_; +}; + +/****************************************************************************** + * template implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +GeneticScheduler::GeneticScheduler(Graph &graph, const ObjFunc &func, + const Parameters &par) +: graph_(graph) +, func_(func) +, par_(par) +{ + gen_.seed(par_.seed); +} + +// access ////////////////////////////////////////////////////////////////////// +template +const typename GeneticScheduler::Gene & +GeneticScheduler::getMinSchedule(void) +{ + return population_.begin()->second; +} + +template +int GeneticScheduler::getMinValue(void) +{ + return population_.begin()->first; +} + +// breed a new generation ////////////////////////////////////////////////////// +template +void GeneticScheduler::nextGeneration(void) +{ + // random initialization of the population if necessary + if (population_.size() != par_.popSize) + { + initPopulation(); + } + LOG(Debug) << "Starting population:\n" << *this << std::endl; + + // random mutations + //PARALLEL_FOR_LOOP + for (unsigned int i = 0; i < par_.popSize; ++i) + { + doMutation(); + } + LOG(Debug) << "After mutations:\n" << *this << std::endl; + + // mating + //PARALLEL_FOR_LOOP + for (unsigned int i = 0; i < par_.popSize/2; ++i) + { + doCrossover(); + } + LOG(Debug) << "After mating:\n" << *this << std::endl; + + // grim reaper + auto it = population_.begin(); + + std::advance(it, par_.popSize); + population_.erase(it, population_.end()); + LOG(Debug) << "After grim reaper:\n" << *this << std::endl; +} + +// evolution steps ///////////////////////////////////////////////////////////// +template +void GeneticScheduler::initPopulation(void) +{ + population_.clear(); + for (unsigned int i = 0; i < par_.popSize; ++i) + { + auto p = graph_.topoSort(gen_); + + population_.insert(std::make_pair(func_(p), p)); + } +} + +template +void GeneticScheduler::doCrossover(void) +{ + auto p = selectPair(); + Gene &p1 = *(p.first), &p2 = *(p.second); + Gene c1, c2; + + crossover(c1, c2, p1, p2); + PARALLEL_CRITICAL + { + population_.insert(std::make_pair(func_(c1), c1)); + population_.insert(std::make_pair(func_(c2), c2)); + } +} + +template +void GeneticScheduler::doMutation(void) +{ + std::uniform_real_distribution mdis(0., 1.); + std::uniform_int_distribution pdis(0, population_.size() - 1); + + if (mdis(gen_) < par_.mutationRate) + { + Gene m; + auto it = population_.begin(); + + std::advance(it, pdis(gen_)); + mutation(m, it->second); + PARALLEL_CRITICAL + { + population_.insert(std::make_pair(func_(m), m)); + } + } +} + +// genetic operators /////////////////////////////////////////////////////////// +template +typename GeneticScheduler::GenePair GeneticScheduler::selectPair(void) +{ + std::vector prob; + unsigned int ind; + Gene *p1, *p2; + + for (auto &c: population_) + { + prob.push_back(1./c.first); + } + do + { + double probCpy; + + std::discrete_distribution dis1(prob.begin(), prob.end()); + auto rIt = population_.begin(); + ind = dis1(gen_); + std::advance(rIt, ind); + p1 = &(rIt->second); + probCpy = prob[ind]; + prob[ind] = 0.; + std::discrete_distribution dis2(prob.begin(), prob.end()); + rIt = population_.begin(); + std::advance(rIt, dis2(gen_)); + p2 = &(rIt->second); + prob[ind] = probCpy; + } while (p1 == p2); + + return std::make_pair(p1, p2); +} + +template +void GeneticScheduler::crossover(Gene &c1, Gene &c2, const Gene &p1, + const Gene &p2) +{ + Gene buf; + std::uniform_int_distribution dis(0, p1.size() - 1); + unsigned int cut = dis(gen_); + + c1.clear(); + buf = p2; + for (unsigned int i = 0; i < cut; ++i) + { + c1.push_back(p1[i]); + buf.erase(std::find(buf.begin(), buf.end(), p1[i])); + } + for (unsigned int i = 0; i < buf.size(); ++i) + { + c1.push_back(buf[i]); + } + c2.clear(); + buf = p2; + for (unsigned int i = cut; i < p1.size(); ++i) + { + buf.erase(std::find(buf.begin(), buf.end(), p1[i])); + } + for (unsigned int i = 0; i < buf.size(); ++i) + { + c2.push_back(buf[i]); + } + for (unsigned int i = cut; i < p1.size(); ++i) + { + c2.push_back(p1[i]); + } +} + +template +void GeneticScheduler::mutation(Gene &m, const Gene &c) +{ + Gene buf; + std::uniform_int_distribution dis(0, c.size() - 1); + unsigned int cut = dis(gen_); + Graph g1 = graph_, g2 = graph_; + + for (unsigned int i = 0; i < cut; ++i) + { + g1.removeVertex(c[i]); + } + for (unsigned int i = cut; i < c.size(); ++i) + { + g2.removeVertex(c[i]); + } + if (g1.size() > 0) + { + buf = g1.topoSort(gen_); + } + if (g2.size() > 0) + { + m = g2.topoSort(gen_); + } + for (unsigned int i = cut; i < c.size(); ++i) + { + m.push_back(buf[i - cut]); + } +} + +template +void GeneticScheduler::benchmarkCrossover(const unsigned int nIt) +{ + Gene p1, p2, c1, c2; + double neg = 0., eq = 0., pos = 0., total; + int improvement; + + LOG(Message) << "Benchmarking crossover..." << std::endl; + for (unsigned int i = 0; i < nIt; ++i) + { + p1 = graph_.topoSort(gen_); + p2 = graph_.topoSort(gen_); + crossover(c1, c2, p1, p2); + improvement = (func_(c1) + func_(c2) - func_(p1) - func_(p2))/2; + if (improvement < 0) neg++; else if (improvement == 0) eq++; else pos++; + } + total = neg + eq + pos; + LOG(Message) << " -: " << neg/total << " =: " << eq/total + << " +: " << pos/total << std::endl; +} + +END_HADRONS_NAMESPACE + +#endif // Hadrons_GeneticScheduler_hpp_ diff --git a/extras/Hadrons/Global.cc b/extras/Hadrons/Global.cc new file mode 100644 index 00000000..7b0b8fb6 --- /dev/null +++ b/extras/Hadrons/Global.cc @@ -0,0 +1,82 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Global.cc + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#include + +using namespace Grid; +using namespace QCD; +using namespace Hadrons; + +HadronsLogger Hadrons::HadronsLogError(1,"Error"); +HadronsLogger Hadrons::HadronsLogWarning(1,"Warning"); +HadronsLogger Hadrons::HadronsLogMessage(1,"Message"); +HadronsLogger Hadrons::HadronsLogIterative(1,"Iterative"); +HadronsLogger Hadrons::HadronsLogDebug(1,"Debug"); + +// pretty size formatting ////////////////////////////////////////////////////// +std::string Hadrons::sizeString(long unsigned int bytes) + +{ + constexpr unsigned int bufSize = 256; + const char *suffixes[7] = {"", "K", "M", "G", "T", "P", "E"}; + char buf[256]; + long unsigned int s = 0; + double count = bytes; + + while (count >= 1024 && s < 7) + { + s++; + count /= 1024; + } + if (count - floor(count) == 0.0) + { + snprintf(buf, bufSize, "%d %sB", (int)count, suffixes[s]); + } + else + { + snprintf(buf, bufSize, "%.1f %sB", count, suffixes[s]); + } + + return std::string(buf); +} + +// type utilities ////////////////////////////////////////////////////////////// +constexpr unsigned int maxNameSize = 1024u; + +std::string Hadrons::typeName(const std::type_info *info) +{ + char *buf; + std::string name; + + buf = abi::__cxa_demangle(info->name(), nullptr, nullptr, nullptr); + name = buf; + free(buf); + + return name; +} diff --git a/extras/Hadrons/Global.hpp b/extras/Hadrons/Global.hpp new file mode 100644 index 00000000..81afab13 --- /dev/null +++ b/extras/Hadrons/Global.hpp @@ -0,0 +1,150 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Global.hpp + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_Global_hpp_ +#define Hadrons_Global_hpp_ + +#include +#include +#include +#include + +#define BEGIN_HADRONS_NAMESPACE \ +namespace Grid {\ +using namespace QCD;\ +namespace Hadrons {\ +using Grid::operator<<; +#define END_HADRONS_NAMESPACE }} + +#define BEGIN_MODULE_NAMESPACE(name)\ +namespace name {\ +using Grid::operator<<; +#define END_MODULE_NAMESPACE } + +/* the 'using Grid::operator<<;' statement prevents a very nasty compilation + * error with GCC 5 (clang & GCC 6 compile fine without it). + */ + +// FIXME: find a way to do that in a more general fashion +#ifndef FIMPL +#define FIMPL WilsonImplR +#endif + +BEGIN_HADRONS_NAMESPACE + +// type aliases +#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 +{ +public: + HadronsLogger(int on, std::string nm): Logger("Hadrons", on, nm, + GridLogColours, "BLACK"){}; +}; + +#define LOG(channel) std::cout << HadronsLog##channel +#define HADRON_ERROR(msg)\ +LOG(Error) << msg << " (" << __FUNCTION__ << " at " << __FILE__ << ":"\ + << __LINE__ << ")" << std::endl;\ +abort(); + +#define DEBUG_VAR(var) LOG(Debug) << #var << "= " << (var) << std::endl; + +extern HadronsLogger HadronsLogError; +extern HadronsLogger HadronsLogWarning; +extern HadronsLogger HadronsLogMessage; +extern HadronsLogger HadronsLogIterative; +extern HadronsLogger HadronsLogDebug; + +// singleton pattern +#define SINGLETON(name)\ +public:\ + name(const name &e) = delete;\ + void operator=(const name &e) = delete;\ + static name & getInstance(void)\ + {\ + static name e;\ + return e;\ + }\ +private:\ + name(void); + +#define SINGLETON_DEFCTOR(name)\ +public:\ + name(const name &e) = delete;\ + void operator=(const name &e) = delete;\ + static name & getInstance(void)\ + {\ + static name e;\ + return e;\ + }\ +private:\ + name(void) = default; + +// pretty size formating +std::string sizeString(long unsigned int bytes); + +// type utilities +template +const std::type_info * typeIdPt(const T &x) +{ + return &typeid(x); +} + +std::string typeName(const std::type_info *info); + +template +const std::type_info * typeIdPt(void) +{ + return &typeid(T); +} + +template +std::string typeName(const T &x) +{ + return typeName(typeIdPt(x)); +} + +template +std::string typeName(void) +{ + return typeName(typeIdPt()); +} + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Global_hpp_ diff --git a/extras/Hadrons/Graph.hpp b/extras/Hadrons/Graph.hpp new file mode 100644 index 00000000..df255517 --- /dev/null +++ b/extras/Hadrons/Graph.hpp @@ -0,0 +1,760 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Graph.hpp + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_Graph_hpp_ +#define Hadrons_Graph_hpp_ + +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Oriented graph class * + ******************************************************************************/ +// I/O for edges +template +std::ostream & operator<<(std::ostream &out, const std::pair &e) +{ + out << "\"" << e.first << "\" -> \"" << e.second << "\""; + + return out; +} + +// main class +template +class Graph +{ +public: + typedef std::pair Edge; +public: + // constructor + Graph(void); + // destructor + virtual ~Graph(void) = default; + // access + void addVertex(const T &value); + void addEdge(const Edge &e); + void addEdge(const T &start, const T &end); + std::vector getVertices(void) const; + void removeVertex(const T &value); + void removeEdge(const Edge &e); + void removeEdge(const T &start, const T &end); + unsigned int size(void) const; + // tests + bool gotValue(const T &value) const; + // graph topological manipulations + std::vector getAdjacentVertices(const T &value) const; + std::vector getChildren(const T &value) const; + std::vector getParents(const T &value) const; + std::vector getRoots(void) const; + std::vector> getConnectedComponents(void) const; + std::vector topoSort(void); + template + std::vector topoSort(Gen &gen); + std::vector> allTopoSort(void); + // I/O + friend std::ostream & operator<<(std::ostream &out, const Graph &g) + { + out << "{"; + for (auto &e: g.edgeSet_) + { + out << e << ", "; + } + if (g.edgeSet_.size() != 0) + { + out << "\b\b"; + } + out << "}"; + + return out; + } +private: + // vertex marking + void mark(const T &value, const bool doMark = true); + void markAll(const bool doMark = true); + void unmark(const T &value); + void unmarkAll(void); + bool isMarked(const T &value) const; + const T * getFirstMarked(const bool isMarked = true) const; + template + const T * getRandomMarked(const bool isMarked, Gen &gen); + const T * getFirstUnmarked(void) const; + template + const T * getRandomUnmarked(Gen &gen); + // prune marked/unmarked vertices + void removeMarked(const bool isMarked = true); + void removeUnmarked(void); + // depth-first search marking + void depthFirstSearch(void); + void depthFirstSearch(const T &root); +private: + std::map isMarked_; + std::set edgeSet_; +}; + +// build depedency matrix from topological sorts +template +std::map> +makeDependencyMatrix(const std::vector> &topSort); + +/****************************************************************************** + * template implementation * + ****************************************************************************** + * in all the following V is the number of vertex and E is the number of edge + * in the worst case E = V^2 + */ + +// constructor ///////////////////////////////////////////////////////////////// +template +Graph::Graph(void) +{} + +// access ////////////////////////////////////////////////////////////////////// +// complexity: log(V) +template +void Graph::addVertex(const T &value) +{ + isMarked_[value] = false; +} + +// complexity: O(log(V)) +template +void Graph::addEdge(const Edge &e) +{ + addVertex(e.first); + addVertex(e.second); + edgeSet_.insert(e); +} + +// complexity: O(log(V)) +template +void Graph::addEdge(const T &start, const T &end) +{ + addEdge(Edge(start, end)); +} + +template +std::vector Graph::getVertices(void) const +{ + std::vector vertex; + + for (auto &v: isMarked_) + { + vertex.push_back(v.first); + } + + return vertex; +} + +// complexity: O(V*log(V)) +template +void Graph::removeVertex(const T &value) +{ + // remove vertex from the mark table + auto vIt = isMarked_.find(value); + + if (vIt != isMarked_.end()) + { + isMarked_.erase(vIt); + } + else + { + HADRON_ERROR("vertex " << value << " does not exists"); + } + + // remove all edges containing the vertex + auto pred = [&value](const Edge &e) + { + return ((e.first == value) or (e.second == value)); + }; + auto eIt = find_if(edgeSet_.begin(), edgeSet_.end(), pred); + + while (eIt != edgeSet_.end()) + { + edgeSet_.erase(eIt); + eIt = find_if(edgeSet_.begin(), edgeSet_.end(), pred); + } +} + +// complexity: O(log(V)) +template +void Graph::removeEdge(const Edge &e) +{ + auto eIt = edgeSet_.find(e); + + if (eIt != edgeSet_.end()) + { + edgeSet_.erase(eIt); + } + else + { + HADRON_ERROR("edge " << e << " does not exists"); + } +} + +// complexity: O(log(V)) +template +void Graph::removeEdge(const T &start, const T &end) +{ + removeEdge(Edge(start, end)); +} + +// complexity: O(1) +template +unsigned int Graph::size(void) const +{ + return isMarked_.size(); +} + +// tests /////////////////////////////////////////////////////////////////////// +// complexity: O(log(V)) +template +bool Graph::gotValue(const T &value) const +{ + auto it = isMarked_.find(value); + + if (it == isMarked_.end()) + { + return false; + } + else + { + return true; + } +} + +// vertex marking ////////////////////////////////////////////////////////////// +// complexity: O(log(V)) +template +void Graph::mark(const T &value, const bool doMark) +{ + if (gotValue(value)) + { + isMarked_[value] = doMark; + } + else + { + HADRON_ERROR("vertex " << value << " does not exists"); + } +} + +// complexity: O(V*log(V)) +template +void Graph::markAll(const bool doMark) +{ + for (auto &v: isMarked_) + { + mark(v.first, doMark); + } +} + +// complexity: O(log(V)) +template +void Graph::unmark(const T &value) +{ + mark(value, false); +} + +// complexity: O(V*log(V)) +template +void Graph::unmarkAll(void) +{ + markAll(false); +} + +// complexity: O(log(V)) +template +bool Graph::isMarked(const T &value) const +{ + if (gotValue(value)) + { + return isMarked_.at(value); + } + else + { + HADRON_ERROR("vertex " << value << " does not exists"); + + return false; + } +} + +// complexity: O(log(V)) +template +const T * Graph::getFirstMarked(const bool isMarked) const +{ + auto pred = [&isMarked](const std::pair &v) + { + return (v.second == isMarked); + }; + auto vIt = std::find_if(isMarked_.begin(), isMarked_.end(), pred); + + if (vIt != isMarked_.end()) + { + return &(vIt->first); + } + else + { + return nullptr; + } +} + +// complexity: O(log(V)) +template +template +const T * Graph::getRandomMarked(const bool isMarked, Gen &gen) +{ + auto pred = [&isMarked](const std::pair &v) + { + return (v.second == isMarked); + }; + std::uniform_int_distribution dis(0, size() - 1); + auto rIt = isMarked_.begin(); + + std::advance(rIt, dis(gen)); + auto vIt = std::find_if(rIt, isMarked_.end(), pred); + if (vIt != isMarked_.end()) + { + return &(vIt->first); + } + else + { + vIt = std::find_if(isMarked_.begin(), rIt, pred); + if (vIt != rIt) + { + return &(vIt->first); + } + else + { + return nullptr; + } + } +} + +// complexity: O(log(V)) +template +const T * Graph::getFirstUnmarked(void) const +{ + return getFirstMarked(false); +} + +// complexity: O(log(V)) +template +template +const T * Graph::getRandomUnmarked(Gen &gen) +{ + return getRandomMarked(false, gen); +} + +// prune marked/unmarked vertices ////////////////////////////////////////////// +// complexity: O(V^2*log(V)) +template +void Graph::removeMarked(const bool isMarked) +{ + auto isMarkedCopy = isMarked_; + + for (auto &v: isMarkedCopy) + { + if (v.second == isMarked) + { + removeVertex(v.first); + } + } +} + +// complexity: O(V^2*log(V)) +template +void Graph::removeUnmarked(void) +{ + removeMarked(false); +} + +// depth-first search marking ////////////////////////////////////////////////// +// complexity: O(V*log(V)) +template +void Graph::depthFirstSearch(void) +{ + depthFirstSearch(isMarked_.begin()->first); +} + +// complexity: O(V*log(V)) +template +void Graph::depthFirstSearch(const T &root) +{ + std::vector adjacentVertex; + + mark(root); + adjacentVertex = getAdjacentVertices(root); + for (auto &v: adjacentVertex) + { + if (!isMarked(v)) + { + depthFirstSearch(v); + } + } +} + +// graph topological manipulations ///////////////////////////////////////////// +// complexity: O(V*log(V)) +template +std::vector Graph::getAdjacentVertices(const T &value) const +{ + std::vector adjacentVertex; + + auto pred = [&value](const Edge &e) + { + return ((e.first == value) or (e.second == value)); + }; + auto eIt = find_if(edgeSet_.begin(), edgeSet_.end(), pred); + + while (eIt != edgeSet_.end()) + { + if (eIt->first == value) + { + adjacentVertex.push_back((*eIt).second); + } + else if (eIt->second == value) + { + adjacentVertex.push_back((*eIt).first); + } + eIt = find_if(++eIt, edgeSet_.end(), pred); + } + + return adjacentVertex; +} + +// complexity: O(V*log(V)) +template +std::vector Graph::getChildren(const T &value) const +{ + std::vector child; + + auto pred = [&value](const Edge &e) + { + return (e.first == value); + }; + auto eIt = find_if(edgeSet_.begin(), edgeSet_.end(), pred); + + while (eIt != edgeSet_.end()) + { + child.push_back((*eIt).second); + eIt = find_if(++eIt, edgeSet_.end(), pred); + } + + return child; +} + +// complexity: O(V*log(V)) +template +std::vector Graph::getParents(const T &value) const +{ + std::vector parent; + + auto pred = [&value](const Edge &e) + { + return (e.second == value); + }; + auto eIt = find_if(edgeSet_.begin(), edgeSet_.end(), pred); + + while (eIt != edgeSet_.end()) + { + parent.push_back((*eIt).first); + eIt = find_if(++eIt, edgeSet_.end(), pred); + } + + return parent; +} + +// complexity: O(V^2*log(V)) +template +std::vector Graph::getRoots(void) const +{ + std::vector root; + + for (auto &v: isMarked_) + { + auto parent = getParents(v.first); + + if (parent.size() == 0) + { + root.push_back(v.first); + } + } + + return root; +} + +// complexity: O(V^2*log(V)) +template +std::vector> Graph::getConnectedComponents(void) const +{ + std::vector> res; + Graph copy(*this); + + while (copy.size() > 0) + { + copy.depthFirstSearch(); + res.push_back(copy); + res.back().removeUnmarked(); + res.back().unmarkAll(); + copy.removeMarked(); + copy.unmarkAll(); + } + + return res; +} + +// topological sort using a directed DFS algorithm +// complexity: O(V*log(V)) +template +std::vector Graph::topoSort(void) +{ + std::stack buf; + std::vector res; + const T *vPt; + std::map tmpMarked(isMarked_); + + // visit function + std::function visit = [&](const T &v) + { + if (tmpMarked.at(v)) + { + HADRON_ERROR("cannot topologically sort a cyclic graph"); + } + if (!isMarked(v)) + { + std::vector child = getChildren(v); + + tmpMarked[v] = true; + for (auto &c: child) + { + visit(c); + } + mark(v); + tmpMarked[v] = false; + buf.push(v); + } + }; + + // reset temporary marks + for (auto &v: tmpMarked) + { + tmpMarked.at(v.first) = false; + } + + // loop on unmarked vertices + unmarkAll(); + vPt = getFirstUnmarked(); + while (vPt) + { + visit(*vPt); + vPt = getFirstUnmarked(); + } + unmarkAll(); + + // create result vector + while (!buf.empty()) + { + res.push_back(buf.top()); + buf.pop(); + } + + return res; +} + +// random version of the topological sort +// complexity: O(V*log(V)) +template +template +std::vector Graph::topoSort(Gen &gen) +{ + std::stack buf; + std::vector res; + const T *vPt; + std::map tmpMarked(isMarked_); + + // visit function + std::function visit = [&](const T &v) + { + if (tmpMarked.at(v)) + { + HADRON_ERROR("cannot topologically sort a cyclic graph"); + } + if (!isMarked(v)) + { + std::vector child = getChildren(v); + + tmpMarked[v] = true; + std::shuffle(child.begin(), child.end(), gen); + for (auto &c: child) + { + visit(c); + } + mark(v); + tmpMarked[v] = false; + buf.push(v); + } + }; + + // reset temporary marks + for (auto &v: tmpMarked) + { + tmpMarked.at(v.first) = false; + } + + // loop on unmarked vertices + unmarkAll(); + vPt = getRandomUnmarked(gen); + while (vPt) + { + visit(*vPt); + vPt = getRandomUnmarked(gen); + } + unmarkAll(); + + // create result vector + while (!buf.empty()) + { + res.push_back(buf.top()); + buf.pop(); + } + + return res; +} + +// generate all possible topological sorts +// Y. L. Varol & D. Rotem, Comput. J. 24(1), pp. 83–84, 1981 +// http://comjnl.oupjournals.org/cgi/doi/10.1093/comjnl/24.1.83 +// complexity: O(V*log(V)) (from the paper, but really ?) +template +std::vector> Graph::allTopoSort(void) +{ + std::vector> res; + std::map> iMat; + + // create incidence matrix + for (auto &v1: isMarked_) + for (auto &v2: isMarked_) + { + iMat[v1.first][v2.first] = false; + } + for (auto &v: isMarked_) + { + auto cVec = getChildren(v.first); + + for (auto &c: cVec) + { + iMat[v.first][c] = true; + } + } + + // generate initial topological sort + res.push_back(topoSort()); + + // generate all other topological sorts by permutation + std::vector p = res[0]; + const unsigned int n = size(); + std::vector loc(n); + unsigned int i, k, k1; + T obj_k, obj_k1; + bool isFinal; + + for (unsigned int j = 0; j < n; ++j) + { + loc[j] = j; + } + i = 0; + while (i < n-1) + { + k = loc[i]; + k1 = k + 1; + obj_k = p[k]; + if (k1 >= n) + { + isFinal = true; + obj_k1 = obj_k; + } + else + { + isFinal = false; + obj_k1 = p[k1]; + } + if (iMat[res[0][i]][obj_k1] or isFinal) + { + for (unsigned int l = k; l >= i + 1; --l) + { + p[l] = p[l-1]; + } + p[i] = obj_k; + loc[i] = i; + i++; + } + else + { + p[k] = obj_k1; + p[k1] = obj_k; + loc[i] = k1; + i = 0; + res.push_back(p); + } + } + + return res; +} + +// build depedency matrix from topological sorts /////////////////////////////// +// complexity: something like O(V^2*log(V!)) +template +std::map> +makeDependencyMatrix(const std::vector> &topSort) +{ + std::map> m; + const std::vector &vList = topSort[0]; + + for (auto &v1: vList) + for (auto &v2: vList) + { + bool dep = true; + + for (auto &t: topSort) + { + auto i1 = std::find(t.begin(), t.end(), v1); + auto i2 = std::find(t.begin(), t.end(), v2); + + dep = dep and (i1 - i2 > 0); + if (!dep) break; + } + m[v1][v2] = dep; + } + + return m; +} + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Graph_hpp_ diff --git a/extras/Hadrons/HadronsXmlRun.cc b/extras/Hadrons/HadronsXmlRun.cc new file mode 100644 index 00000000..0dff8f9a --- /dev/null +++ b/extras/Hadrons/HadronsXmlRun.cc @@ -0,0 +1,80 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/HadronsXmlRun.cc + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#include + +using namespace Grid; +using namespace QCD; +using namespace Hadrons; + +int main(int argc, char *argv[]) +{ + // parse command line + std::string parameterFileName, scheduleFileName = ""; + + if (argc < 2) + { + std::cerr << "usage: " << argv[0] << " [] [Grid options]"; + std::cerr << std::endl; + std::exit(EXIT_FAILURE); + } + parameterFileName = argv[1]; + if (argc > 2) + { + if (argv[2][0] != '-') + { + scheduleFileName = argv[2]; + } + } + + // 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; + + // execution + Application application(parameterFileName); + + application.parseParameterFile(parameterFileName); + if (!scheduleFileName.empty()) + { + application.loadSchedule(scheduleFileName); + } + application.run(); + + // epilogue + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +} diff --git a/extras/Hadrons/HadronsXmlSchedule.cc b/extras/Hadrons/HadronsXmlSchedule.cc new file mode 100644 index 00000000..a8ca9a63 --- /dev/null +++ b/extras/Hadrons/HadronsXmlSchedule.cc @@ -0,0 +1,72 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/HadronsXmlSchedule.cc + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#include + +using namespace Grid; +using namespace QCD; +using namespace Hadrons; + +int main(int argc, char *argv[]) +{ + // parse command line + std::string parameterFileName, scheduleFileName; + + if (argc < 3) + { + std::cerr << "usage: " << argv[0] << " [Grid options]"; + std::cerr << std::endl; + std::exit(EXIT_FAILURE); + } + parameterFileName = argv[1]; + scheduleFileName = argv[2]; + + // 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; + + // execution + Application application; + + application.parseParameterFile(parameterFileName); + application.schedule(); + application.printSchedule(); + application.saveSchedule(scheduleFileName); + + // epilogue + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +} diff --git a/extras/Hadrons/Makefile.am b/extras/Hadrons/Makefile.am new file mode 100644 index 00000000..9cb23600 --- /dev/null +++ b/extras/Hadrons/Makefile.am @@ -0,0 +1,29 @@ +lib_LIBRARIES = libHadrons.a +bin_PROGRAMS = HadronsXmlRun HadronsXmlSchedule + +include modules.inc + +libHadrons_a_SOURCES = \ + $(modules_cc) \ + Application.cc \ + Environment.cc \ + Global.cc \ + Module.cc +libHadrons_adir = $(pkgincludedir)/Hadrons +nobase_libHadrons_a_HEADERS = \ + $(modules_hpp) \ + Application.hpp \ + Environment.hpp \ + Factory.hpp \ + GeneticScheduler.hpp \ + Global.hpp \ + Graph.hpp \ + Module.hpp \ + Modules.hpp \ + ModuleFactory.hpp + +HadronsXmlRun_SOURCES = HadronsXmlRun.cc +HadronsXmlRun_LDADD = libHadrons.a -lGrid + +HadronsXmlSchedule_SOURCES = HadronsXmlSchedule.cc +HadronsXmlSchedule_LDADD = libHadrons.a -lGrid diff --git a/extras/Hadrons/Module.cc b/extras/Hadrons/Module.cc new file mode 100644 index 00000000..2549a931 --- /dev/null +++ b/extras/Hadrons/Module.cc @@ -0,0 +1,71 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Module.cc + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#include + +using namespace Grid; +using namespace QCD; +using namespace Hadrons; + +/****************************************************************************** + * ModuleBase implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +ModuleBase::ModuleBase(const std::string name) +: name_(name) +, env_(Environment::getInstance()) +{} + +// access ////////////////////////////////////////////////////////////////////// +std::string ModuleBase::getName(void) const +{ + return name_; +} + +Environment & ModuleBase::env(void) const +{ + return env_; +} + +// get factory registration name if available +std::string ModuleBase::getRegisteredName(void) +{ + HADRON_ERROR("module '" + getName() + "' has a type not registered" + + " in the factory"); +} + +// execution /////////////////////////////////////////////////////////////////// +void ModuleBase::operator()(void) +{ + setup(); + if (!env().isDryRun()) + { + execute(); + } +} diff --git a/extras/Hadrons/Module.hpp b/extras/Hadrons/Module.hpp new file mode 100644 index 00000000..071e254a --- /dev/null +++ b/extras/Hadrons/Module.hpp @@ -0,0 +1,198 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Module.hpp + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_Module_hpp_ +#define Hadrons_Module_hpp_ + +#include +#include + +BEGIN_HADRONS_NAMESPACE + +// module registration macros +#define MODULE_REGISTER(mod, base)\ +class mod: public base\ +{\ +public:\ + typedef base Base;\ + using Base::Base;\ + virtual std::string getRegisteredName(void)\ + {\ + return std::string(#mod);\ + }\ +};\ +class mod##ModuleRegistrar\ +{\ +public:\ + mod##ModuleRegistrar(void)\ + {\ + ModuleFactory &modFac = ModuleFactory::getInstance();\ + modFac.registerBuilder(#mod, [&](const std::string name)\ + {\ + return std::unique_ptr(new mod(name));\ + });\ + }\ +};\ +static mod##ModuleRegistrar mod##ModuleRegistrarInstance; + +#define MODULE_REGISTER_NS(mod, base, ns)\ +class mod: public base\ +{\ +public:\ + typedef base Base;\ + using Base::Base;\ + virtual std::string getRegisteredName(void)\ + {\ + return std::string(#ns "::" #mod);\ + }\ +};\ +class ns##mod##ModuleRegistrar\ +{\ +public:\ + ns##mod##ModuleRegistrar(void)\ + {\ + ModuleFactory &modFac = ModuleFactory::getInstance();\ + modFac.registerBuilder(#ns "::" #mod, [&](const std::string name)\ + {\ + return std::unique_ptr(new ns::mod(name));\ + });\ + }\ +};\ +static ns##mod##ModuleRegistrar ns##mod##ModuleRegistrarInstance; + +#define ARG(...) __VA_ARGS__ + +/****************************************************************************** + * Module class * + ******************************************************************************/ +// base class +class ModuleBase +{ +public: + // constructor + ModuleBase(const std::string name); + // destructor + virtual ~ModuleBase(void) = default; + // access + std::string getName(void) const; + Environment &env(void) const; + // get factory registration name if available + virtual std::string getRegisteredName(void); + // dependencies/products + virtual std::vector getInput(void) = 0; + virtual std::vector getOutput(void) = 0; + // parse parameters + virtual void parseParameters(XmlReader &reader, const std::string name) = 0; + virtual void saveParameters(XmlWriter &writer, const std::string name) = 0; + // setup + virtual void setup(void) {}; + // execution + void operator()(void); + virtual void execute(void) = 0; +private: + std::string name_; + Environment &env_; +}; + +// derived class, templating the parameter class +template +class Module: public ModuleBase +{ +public: + typedef P Par; +public: + // constructor + Module(const std::string name); + // destructor + virtual ~Module(void) = default; + // parse parameters + virtual void parseParameters(XmlReader &reader, const std::string name); + virtual void saveParameters(XmlWriter &writer, const std::string name); + // parameter access + const P & par(void) const; + void setPar(const P &par); +private: + P par_; +}; + +// no parameter type +class NoPar {}; + +template <> +class Module: public ModuleBase +{ +public: + // constructor + Module(const std::string name): ModuleBase(name) {}; + // destructor + virtual ~Module(void) = default; + // parse parameters (do nothing) + virtual void parseParameters(XmlReader &reader, const std::string name) {}; + virtual void saveParameters(XmlWriter &writer, const std::string name) + { + push(writer, "options"); + pop(writer); + }; +}; + +/****************************************************************************** + * Template implementation * + ******************************************************************************/ +template +Module

::Module(const std::string name) +: ModuleBase(name) +{} + +template +void Module

::parseParameters(XmlReader &reader, const std::string name) +{ + read(reader, name, par_); +} + +template +void Module

::saveParameters(XmlWriter &writer, const std::string name) +{ + write(writer, name, par_); +} + +template +const P & Module

::par(void) const +{ + return par_; +} + +template +void Module

::setPar(const P &par) +{ + par_ = par; +} + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Module_hpp_ diff --git a/extras/Hadrons/ModuleFactory.hpp b/extras/Hadrons/ModuleFactory.hpp new file mode 100644 index 00000000..48ab305c --- /dev/null +++ b/extras/Hadrons/ModuleFactory.hpp @@ -0,0 +1,49 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/ModuleFactory.hpp + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_ModuleFactory_hpp_ +#define Hadrons_ModuleFactory_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * ModuleFactory * + ******************************************************************************/ +class ModuleFactory: public Factory +{ + SINGLETON_DEFCTOR(ModuleFactory) +}; + +END_HADRONS_NAMESPACE + +#endif // Hadrons_ModuleFactory_hpp_ diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp new file mode 100644 index 00000000..77ae08b7 --- /dev/null +++ b/extras/Hadrons/Modules.hpp @@ -0,0 +1,40 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules.hpp + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include diff --git a/extras/Hadrons/Modules/MAction/DWF.hpp b/extras/Hadrons/Modules/MAction/DWF.hpp new file mode 100644 index 00000000..49861e3e --- /dev/null +++ b/extras/Hadrons/Modules/MAction/DWF.hpp @@ -0,0 +1,134 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MAction/DWF.hpp + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_DWF_hpp_ +#define Hadrons_DWF_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Domain wall quark action * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MAction) + +class DWFPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(DWFPar, + std::string, gauge, + unsigned int, Ls, + double , mass, + double , M5); +}; + +template +class TDWF: public Module +{ +public: + TYPE_ALIASES(FImpl,); +public: + // constructor + TDWF(const std::string name); + // destructor + virtual ~TDWF(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(DWF, TDWF, MAction); + +/****************************************************************************** + * 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 = 2*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); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_DWF_hpp_ diff --git a/extras/Hadrons/Modules/MAction/Wilson.hpp b/extras/Hadrons/Modules/MAction/Wilson.hpp new file mode 100644 index 00000000..6ffa997d --- /dev/null +++ b/extras/Hadrons/Modules/MAction/Wilson.hpp @@ -0,0 +1,126 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MAction/Wilson.hpp + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_Wilson_hpp_ +#define Hadrons_Wilson_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * TWilson quark action * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MAction) + +class WilsonPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonPar, + std::string, gauge, + double , mass); +}; + +template +class TWilson: public Module +{ +public: + TYPE_ALIASES(FImpl,); +public: + // constructor + TWilson(const std::string name); + // destructor + virtual ~TWilson(void) = default; + // dependencies/products + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(Wilson, TWilson, MAction); + +/****************************************************************************** + * 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 = 2*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); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Wilson_hpp_ diff --git a/extras/Hadrons/Modules/MContraction/Baryon.hpp b/extras/Hadrons/Modules/MContraction/Baryon.hpp new file mode 100644 index 00000000..6bccda2e --- /dev/null +++ b/extras/Hadrons/Modules/MContraction/Baryon.hpp @@ -0,0 +1,131 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MContraction/Baryon.hpp + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_Baryon_hpp_ +#define Hadrons_Baryon_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Baryon * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MContraction) + +class BaryonPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(BaryonPar, + std::string, q1, + std::string, q2, + std::string, q3, + std::string, output); +}; + +template +class TBaryon: public Module +{ +public: + TYPE_ALIASES(FImpl1, 1); + TYPE_ALIASES(FImpl2, 2); + TYPE_ALIASES(FImpl3, 3); + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::vector>>, corr); + }; +public: + // constructor + TBaryon(const std::string name); + // destructor + virtual ~TBaryon(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(Baryon, ARG(TBaryon), MContraction); + +/****************************************************************************** + * TBaryon implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TBaryon::TBaryon(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TBaryon::getInput(void) +{ + std::vector input = {par().q1, par().q2, par().q3}; + + return input; +} + +template +std::vector TBaryon::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TBaryon::execute(void) +{ + LOG(Message) << "Computing baryon contractions '" << getName() << "' using" + << " quarks '" << par().q1 << "', '" << par().q2 << "', and '" + << par().q3 << "'" << std::endl; + + XmlWriter writer(par().output); + PropagatorField1 &q1 = *env().template getObject(par().q1); + PropagatorField2 &q2 = *env().template getObject(par().q2); + PropagatorField3 &q3 = *env().template getObject(par().q2); + LatticeComplex c(env().getGrid()); + Result result; + + // FIXME: do contractions + + write(writer, "meson", result); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Baryon_hpp_ diff --git a/extras/Hadrons/Modules/MContraction/Meson.hpp b/extras/Hadrons/Modules/MContraction/Meson.hpp new file mode 100644 index 00000000..dc4f60ee --- /dev/null +++ b/extras/Hadrons/Modules/MContraction/Meson.hpp @@ -0,0 +1,148 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MContraction/Meson.hpp + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_Meson_hpp_ +#define Hadrons_Meson_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * TMeson * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MContraction) + +class MesonPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(MesonPar, + std::string, q1, + std::string, q2, + std::string, output); +}; + +template +class TMeson: public Module +{ +public: + TYPE_ALIASES(FImpl1, 1); + TYPE_ALIASES(FImpl2, 2); + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::vector>>, corr); + }; +public: + // constructor + TMeson(const std::string name); + // destructor + virtual ~TMeson(void) = default; + // dependencies/products + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(Meson, ARG(TMeson), MContraction); + +/****************************************************************************** + * 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 contractions '" << 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); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Meson_hpp_ diff --git a/extras/Hadrons/Modules/MGauge/Load.cc b/extras/Hadrons/Modules/MGauge/Load.cc new file mode 100644 index 00000000..e5ee8abb --- /dev/null +++ b/extras/Hadrons/Modules/MGauge/Load.cc @@ -0,0 +1,78 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MGauge/Load.cc + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MGauge; + +/****************************************************************************** +* TLoad implementation * +******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +TLoad::TLoad(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +std::vector TLoad::getInput(void) +{ + std::vector in; + + return in; +} + +std::vector TLoad::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +void TLoad::setup(void) +{ + env().registerLattice(getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +void TLoad::execute(void) +{ + NerscField header; + std::string fileName = par().file + "." + + std::to_string(env().getTrajectory()); + + LOG(Message) << "Loading NERSC configuration from file '" << fileName + << "'" << std::endl; + LatticeGaugeField &U = *env().createLattice(getName()); + NerscIO::readConfiguration(U, header, fileName); + LOG(Message) << "NERSC header:" << std::endl; + dump_nersc_header(header, LOG(Message)); +} diff --git a/extras/Hadrons/Modules/MGauge/Load.hpp b/extras/Hadrons/Modules/MGauge/Load.hpp new file mode 100644 index 00000000..c41f9b8c --- /dev/null +++ b/extras/Hadrons/Modules/MGauge/Load.hpp @@ -0,0 +1,73 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MGauge/Load.hpp + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_Load_hpp_ +#define Hadrons_Load_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Load a NERSC configuration * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MGauge) + +class LoadPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(LoadPar, + std::string, file); +}; + +class TLoad: public Module +{ +public: + // constructor + TLoad(const std::string name); + // destructor + virtual ~TLoad(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(Load, TLoad, MGauge); + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Load_hpp_ diff --git a/extras/Hadrons/Modules/MGauge/Random.cc b/extras/Hadrons/Modules/MGauge/Random.cc new file mode 100644 index 00000000..c10fdfc3 --- /dev/null +++ b/extras/Hadrons/Modules/MGauge/Random.cc @@ -0,0 +1,69 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MGauge/Random.cc + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MGauge; + +/****************************************************************************** +* TRandom implementation * +******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +TRandom::TRandom(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +std::vector TRandom::getInput(void) +{ + return std::vector(); +} + +std::vector TRandom::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +void TRandom::setup(void) +{ + env().registerLattice(getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +void TRandom::execute(void) +{ + LOG(Message) << "Generating random gauge configuration" << std::endl; + LatticeGaugeField &U = *env().createLattice(getName()); + SU3::HotConfiguration(*env().get4dRng(), U); +} diff --git a/extras/Hadrons/Modules/MGauge/Random.hpp b/extras/Hadrons/Modules/MGauge/Random.hpp new file mode 100644 index 00000000..e3fbcf1a --- /dev/null +++ b/extras/Hadrons/Modules/MGauge/Random.hpp @@ -0,0 +1,66 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MGauge/Random.hpp + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_Random_hpp_ +#define Hadrons_Random_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Random gauge * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MGauge) + +class TRandom: public Module +{ +public: + // constructor + TRandom(const std::string name); + // destructor + virtual ~TRandom(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(Random, TRandom, MGauge); + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Random_hpp_ diff --git a/extras/Hadrons/Modules/MGauge/Unit.cc b/extras/Hadrons/Modules/MGauge/Unit.cc new file mode 100644 index 00000000..18d75c59 --- /dev/null +++ b/extras/Hadrons/Modules/MGauge/Unit.cc @@ -0,0 +1,69 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MGauge/Unit.cc + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MGauge; + +/****************************************************************************** +* TUnit implementation * +******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +TUnit::TUnit(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +std::vector TUnit::getInput(void) +{ + return std::vector(); +} + +std::vector TUnit::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +void TUnit::setup(void) +{ + env().registerLattice(getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +void TUnit::execute(void) +{ + LOG(Message) << "Creating unit gauge configuration" << std::endl; + LatticeGaugeField &U = *env().createLattice(getName()); + SU3::ColdConfiguration(*env().get4dRng(), U); +} diff --git a/extras/Hadrons/Modules/MGauge/Unit.hpp b/extras/Hadrons/Modules/MGauge/Unit.hpp new file mode 100644 index 00000000..2ff10bfd --- /dev/null +++ b/extras/Hadrons/Modules/MGauge/Unit.hpp @@ -0,0 +1,66 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MGauge/Unit.hpp + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_Unit_hpp_ +#define Hadrons_Unit_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Unit gauge * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MGauge) + +class TUnit: public Module +{ +public: + // constructor + TUnit(const std::string name); + // destructor + virtual ~TUnit(void) = default; + // dependencies/products + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(Unit, TUnit, MGauge); + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Unit_hpp_ diff --git a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp new file mode 100644 index 00000000..d7220271 --- /dev/null +++ b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp @@ -0,0 +1,132 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MSolver/RBPrecCG.hpp + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_RBPrecCG_hpp_ +#define Hadrons_RBPrecCG_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Schur red-black preconditioned CG * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MSolver) + +class RBPrecCGPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(RBPrecCGPar, + std::string, action, + double , residual); +}; + +template +class TRBPrecCG: public Module +{ +public: + TYPE_ALIASES(FImpl,); +public: + // constructor + TRBPrecCG(const std::string name); + // destructor + virtual ~TRBPrecCG(void) = default; + // dependencies/products + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(RBPrecCG, TRBPrecCG, MSolver); + +/****************************************************************************** + * 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)); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_RBPrecCG_hpp_ diff --git a/extras/Hadrons/Modules/MSource/Point.hpp b/extras/Hadrons/Modules/MSource/Point.hpp new file mode 100644 index 00000000..a0ecbc2a --- /dev/null +++ b/extras/Hadrons/Modules/MSource/Point.hpp @@ -0,0 +1,135 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MSource/Point.hpp + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_Point_hpp_ +#define Hadrons_Point_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/* + + Point source + ------------ + * src_x = delta_x,position + + * options: + - position: space-separated integer sequence (e.g. "0 1 1 0") + + */ + +/****************************************************************************** + * TPoint * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MSource) + +class PointPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(PointPar, + std::string, position); +}; + +template +class TPoint: public Module +{ +public: + TYPE_ALIASES(FImpl,); +public: + // constructor + TPoint(const std::string name); + // destructor + virtual ~TPoint(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(Point, TPoint, MSource); + +/****************************************************************************** + * 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); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Point_hpp_ diff --git a/extras/Hadrons/Modules/MSource/SeqGamma.hpp b/extras/Hadrons/Modules/MSource/SeqGamma.hpp new file mode 100644 index 00000000..611b0108 --- /dev/null +++ b/extras/Hadrons/Modules/MSource/SeqGamma.hpp @@ -0,0 +1,164 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MSource/SeqGamma.hpp + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_SeqGamma_hpp_ +#define Hadrons_SeqGamma_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/* + + Sequential source + ----------------------------- + * src_x = q_x * theta(x_3 - tA) * theta(tB - x_3) * gamma * exp(i x.mom) + + * options: + - q: input propagator (string) + - tA: begin timeslice (integer) + - tB: end timesilce (integer) + - gamma: gamma product to insert (integer) + - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.") + + */ + +/****************************************************************************** + * SeqGamma * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MSource) + +class SeqGammaPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(SeqGammaPar, + std::string, q, + unsigned int, tA, + unsigned int, tB, + unsigned int, gamma, + std::string, mom); +}; + +template +class TSeqGamma: public Module +{ +public: + TYPE_ALIASES(FImpl,); +public: + // constructor + TSeqGamma(const std::string name); + // destructor + virtual ~TSeqGamma(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(SeqGamma, TSeqGamma, MSource); + +/****************************************************************************** + * TSeqGamma implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TSeqGamma::TSeqGamma(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TSeqGamma::getInput(void) +{ + std::vector in = {par().q}; + + return in; +} + +template +std::vector TSeqGamma::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TSeqGamma::setup(void) +{ + env().template registerLattice(getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TSeqGamma::execute(void) +{ + if (par().tA == par().tB) + { + LOG(Message) << "Generating gamma_" << par().gamma + << " sequential source at t= " << par().tA << std::endl; + } + else + { + LOG(Message) << "Generating gamma_" << par().gamma + << " sequential source for " + << par().tA << " <= t <= " << par().tB << std::endl; + } + PropagatorField &src = *env().template createLattice(getName()); + PropagatorField &q = *env().template getObject(par().q); + Lattice> t(env().getGrid()); + LatticeComplex ph(env().getGrid()), coor(env().getGrid()); + SpinMatrix g; + std::vector p; + Complex i(0.0,1.0); + + g = makeGammaProd(par().gamma); + p = strToVec(par().mom); + ph = zero; + for(unsigned int mu = 0; mu < env().getNd(); mu++) + { + LatticeCoordinate(coor, mu); + ph = ph + p[mu]*coor; + } + ph = exp(i*ph); + LatticeCoordinate(t, Tp); + src = where((t >= par().tA) and (t <= par().tB), g*ph*q, 0.*q); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_SeqGamma_hpp_ diff --git a/extras/Hadrons/Modules/MSource/Z2.hpp b/extras/Hadrons/Modules/MSource/Z2.hpp new file mode 100644 index 00000000..cd5727be --- /dev/null +++ b/extras/Hadrons/Modules/MSource/Z2.hpp @@ -0,0 +1,151 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MSource/Z2.hpp + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_Z2_hpp_ +#define Hadrons_Z2_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/* + + Z_2 stochastic source + ----------------------------- + * src_x = eta_x * theta(x_3 - tA) * theta(tB - x_3) + + the eta_x are independent uniform random numbers in {+/- 1 +/- i} + + * options: + - tA: begin timeslice (integer) + - tB: end timesilce (integer) + + */ + +/****************************************************************************** + * Z2 stochastic source * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MSource) + +class Z2Par: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Z2Par, + unsigned int, tA, + unsigned int, tB); +}; + +template +class TZ2: public Module +{ +public: + TYPE_ALIASES(FImpl,); +public: + // constructor + TZ2(const std::string name); + // destructor + virtual ~TZ2(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(Z2, TZ2, MSource); + +/****************************************************************************** + * 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; +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Z2_hpp_ diff --git a/extras/Hadrons/Modules/Quark.hpp b/extras/Hadrons/Modules/Quark.hpp new file mode 100644 index 00000000..e441a096 --- /dev/null +++ b/extras/Hadrons/Modules/Quark.hpp @@ -0,0 +1,185 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/Quark.hpp + +Copyright (C) 2015 +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_Quark_hpp_ +#define Hadrons_Quark_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * 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}; +}; + +MODULE_REGISTER(Quark, TQuark); + +/****************************************************************************** + * TQuark implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TQuark::TQuark(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TQuark::getInput(void) +{ + std::vector in = {par().source, par().solver}; + + return in; +} + +template +std::vector TQuark::getOutput(void) +{ + std::vector out = {getName(), getName() + "_5d"}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TQuark::setup(void) +{ + Ls_ = env().getObjectLs(par().solver); + env().template registerLattice(getName()); + if (Ls_ > 1) + { + env().template registerLattice(getName() + "_5d", Ls_); + } +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TQuark::execute(void) +{ + LOG(Message) << "Computing quark propagator '" << getName() << "'" + << std::endl; + + 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().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; + // source conversion for 4D sources + if (!env().isObject5d(par().source)) + { + if (Ls_ == 1) + { + PropToFerm(source, fullSrc, s, c); + } + else + { + source = zero; + PropToFerm(tmp, fullSrc, s, c); + InsertSlice(tmp, source, 0, 0); + InsertSlice(tmp, source, Ls_-1, 0); + axpby_ssp_pplus(source, 0., source, 1., source, 0, 0); + axpby_ssp_pminus(source, 0., source, 1., source, Ls_-1, Ls_-1); + } + } + // source conversion for 5D sources + else + { + if (Ls_ != env().getObjectLs(par().source)) + { + HADRON_ERROR("Ls mismatch between quark action and source"); + } + else + { + PropToFerm(source, fullSrc, s, c); + } + } + sol = zero; + solver(sol, source); + FermToProp(prop, sol, s, c); + // create 4D propagators from 5D one if necessary + if (Ls_ > 1) + { + 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); + ExtractSlice(tmp, sol, 0, 0); + FermToProp(p4d, tmp, s, c); + } + } +} + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Quark_hpp_ diff --git a/extras/Hadrons/Modules/templates/Module.cc.template b/extras/Hadrons/Modules/templates/Module.cc.template new file mode 100644 index 00000000..0c509d6d --- /dev/null +++ b/extras/Hadrons/Modules/templates/Module.cc.template @@ -0,0 +1,39 @@ +#include + +using namespace Grid; +using namespace Hadrons; + +/****************************************************************************** +* T___FILEBASENAME___ implementation * +******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +T___FILEBASENAME___::T___FILEBASENAME___(const std::string name) +: Module<___FILEBASENAME___Par>(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +std::vector T___FILEBASENAME___::getInput(void) +{ + std::vector in; + + return in; +} + +std::vector T___FILEBASENAME___::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +void T___FILEBASENAME___::setup(void) +{ + +} + +// execution /////////////////////////////////////////////////////////////////// +void T___FILEBASENAME___::execute(void) +{ + +} diff --git a/extras/Hadrons/Modules/templates/Module.hpp.template b/extras/Hadrons/Modules/templates/Module.hpp.template new file mode 100644 index 00000000..fb43260f --- /dev/null +++ b/extras/Hadrons/Modules/templates/Module.hpp.template @@ -0,0 +1,40 @@ +#ifndef Hadrons____FILEBASENAME____hpp_ +#define Hadrons____FILEBASENAME____hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * ___FILEBASENAME___ * + ******************************************************************************/ +class ___FILEBASENAME___Par: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(___FILEBASENAME___Par, + unsigned int, i); +}; + +class T___FILEBASENAME___: public Module<___FILEBASENAME___Par> +{ +public: + // constructor + T___FILEBASENAME___(const std::string name); + // destructor + virtual ~T___FILEBASENAME___(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER(___FILEBASENAME___, T___FILEBASENAME___); + +END_HADRONS_NAMESPACE + +#endif // Hadrons____FILEBASENAME____hpp_ diff --git a/extras/Hadrons/Modules/templates/Module_in_NS.cc.template b/extras/Hadrons/Modules/templates/Module_in_NS.cc.template new file mode 100644 index 00000000..8b2a0ec0 --- /dev/null +++ b/extras/Hadrons/Modules/templates/Module_in_NS.cc.template @@ -0,0 +1,40 @@ +#include + +using namespace Grid; +using namespace Hadrons; +using namespace ___NAMESPACE___; + +/****************************************************************************** +* T___FILEBASENAME___ implementation * +******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +T___FILEBASENAME___::T___FILEBASENAME___(const std::string name) +: Module<___FILEBASENAME___Par>(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +std::vector T___FILEBASENAME___::getInput(void) +{ + std::vector in; + + return in; +} + +std::vector T___FILEBASENAME___::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +void T___FILEBASENAME___::setup(void) +{ + +} + +// execution /////////////////////////////////////////////////////////////////// +void T___FILEBASENAME___::execute(void) +{ + +} diff --git a/extras/Hadrons/Modules/templates/Module_in_NS.hpp.template b/extras/Hadrons/Modules/templates/Module_in_NS.hpp.template new file mode 100644 index 00000000..ece2bb58 --- /dev/null +++ b/extras/Hadrons/Modules/templates/Module_in_NS.hpp.template @@ -0,0 +1,44 @@ +#ifndef Hadrons____FILEBASENAME____hpp_ +#define Hadrons____FILEBASENAME____hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * ___FILEBASENAME___ * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(___NAMESPACE___) + +class ___FILEBASENAME___Par: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(___FILEBASENAME___Par, + unsigned int, i); +}; + +class T___FILEBASENAME___: public Module<___FILEBASENAME___Par> +{ +public: + // constructor + T___FILEBASENAME___(const std::string name); + // destructor + virtual ~T___FILEBASENAME___(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(___FILEBASENAME___, T___FILEBASENAME___, ___NAMESPACE___); + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons____FILEBASENAME____hpp_ diff --git a/extras/Hadrons/Modules/templates/Module_tmp.hpp.template b/extras/Hadrons/Modules/templates/Module_tmp.hpp.template new file mode 100644 index 00000000..2ee053a9 --- /dev/null +++ b/extras/Hadrons/Modules/templates/Module_tmp.hpp.template @@ -0,0 +1,81 @@ +#ifndef Hadrons____FILEBASENAME____hpp_ +#define Hadrons____FILEBASENAME____hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * ___FILEBASENAME___ * + ******************************************************************************/ +class ___FILEBASENAME___Par: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(___FILEBASENAME___Par, + unsigned int, i); +}; + +template +class T___FILEBASENAME___: public Module<___FILEBASENAME___Par> +{ +public: + // constructor + T___FILEBASENAME___(const std::string name); + // destructor + virtual ~T___FILEBASENAME___(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER(___FILEBASENAME___, T___FILEBASENAME___); + +/****************************************************************************** + * T___FILEBASENAME___ implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +T___FILEBASENAME___::T___FILEBASENAME___(const std::string name) +: Module<___FILEBASENAME___Par>(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector T___FILEBASENAME___::getInput(void) +{ + std::vector in; + + return in; +} + +template +std::vector T___FILEBASENAME___::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void T___FILEBASENAME___::setup(void) +{ + +} + +// execution /////////////////////////////////////////////////////////////////// +template +void T___FILEBASENAME___::execute(void) +{ + +} + +END_HADRONS_NAMESPACE + +#endif // Hadrons____FILEBASENAME____hpp_ diff --git a/extras/Hadrons/Modules/templates/Module_tmp_in_NS.hpp.template b/extras/Hadrons/Modules/templates/Module_tmp_in_NS.hpp.template new file mode 100644 index 00000000..a330652d --- /dev/null +++ b/extras/Hadrons/Modules/templates/Module_tmp_in_NS.hpp.template @@ -0,0 +1,85 @@ +#ifndef Hadrons____FILEBASENAME____hpp_ +#define Hadrons____FILEBASENAME____hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * ___FILEBASENAME___ * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(___NAMESPACE___) + +class ___FILEBASENAME___Par: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(___FILEBASENAME___Par, + unsigned int, i); +}; + +template +class T___FILEBASENAME___: public Module<___FILEBASENAME___Par> +{ +public: + // constructor + T___FILEBASENAME___(const std::string name); + // destructor + virtual ~T___FILEBASENAME___(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(___FILEBASENAME___, T___FILEBASENAME___, ___NAMESPACE___); + +/****************************************************************************** + * T___FILEBASENAME___ implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +T___FILEBASENAME___::T___FILEBASENAME___(const std::string name) +: Module<___FILEBASENAME___Par>(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector T___FILEBASENAME___::getInput(void) +{ + std::vector in; + + return in; +} + +template +std::vector T___FILEBASENAME___::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void T___FILEBASENAME___::setup(void) +{ + +} + +// execution /////////////////////////////////////////////////////////////////// +template +void T___FILEBASENAME___::execute(void) +{ + +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons____FILEBASENAME____hpp_ diff --git a/extras/Hadrons/add_module.sh b/extras/Hadrons/add_module.sh new file mode 100755 index 00000000..d5d23ea4 --- /dev/null +++ b/extras/Hadrons/add_module.sh @@ -0,0 +1,31 @@ +#!/usr/bin/env bash + +if (( $# != 1 && $# != 2)); then + echo "usage: `basename $0` []" 1>&2 + exit 1 +fi +NAME=$1 +NS=$2 + +if (( $# == 1 )); then + if [ -e "Modules/${NAME}.cc" ] || [ -e "Modules/${NAME}.hpp" ]; then + echo "error: files Modules/${NAME}.* already exists" 1>&2 + exit 1 + fi + sed "s/___FILEBASENAME___/${NAME}/g" Modules/templates/Module.cc.template > Modules/${NAME}.cc + sed "s/___FILEBASENAME___/${NAME}/g" Modules/templates/Module.hpp.template > Modules/${NAME}.hpp +elif (( $# == 2 )); then + mkdir -p Modules/${NS} + if [ -e "Modules/${NS}/${NAME}.cc" ] || [ -e "Modules/${NS}/${NAME}.hpp" ]; then + echo "error: files Modules/${NS}/${NAME}.* already exists" 1>&2 + exit 1 + fi + TMPCC=".${NS}.${NAME}.tmp.cc" + TMPHPP=".${NS}.${NAME}.tmp.hpp" + sed "s/___FILEBASENAME___/${NAME}/g" Modules/templates/Module_in_NS.cc.template > ${TMPCC} + sed "s/___FILEBASENAME___/${NAME}/g" Modules/templates/Module_in_NS.hpp.template > ${TMPHPP} + sed "s/___NAMESPACE___/${NS}/g" ${TMPCC} > Modules/${NS}/${NAME}.cc + sed "s/___NAMESPACE___/${NS}/g" ${TMPHPP} > Modules/${NS}/${NAME}.hpp + rm -f ${TMPCC} ${TMPHPP} +fi +./make_module_list.sh diff --git a/extras/Hadrons/add_module_template.sh b/extras/Hadrons/add_module_template.sh new file mode 100755 index 00000000..0069fcea --- /dev/null +++ b/extras/Hadrons/add_module_template.sh @@ -0,0 +1,28 @@ +#!/usr/bin/env bash + +if (( $# != 1 && $# != 2)); then + echo "usage: `basename $0` []" 1>&2 + exit 1 +fi +NAME=$1 +NS=$2 + +if (( $# == 1 )); then + if [ -e "Modules/${NAME}.cc" ] || [ -e "Modules/${NAME}.hpp" ]; then + echo "error: files Modules/${NAME}.* already exists" 1>&2 + exit 1 + fi + sed "s/___FILEBASENAME___/${NAME}/g" Modules/templates/Module_tmp.hpp.template > Modules/${NAME}.hpp +elif (( $# == 2 )); then + mkdir -p Modules/${NS} + if [ -e "Modules/${NS}/${NAME}.cc" ] || [ -e "Modules/${NS}/${NAME}.hpp" ]; then + echo "error: files Modules/${NS}/${NAME}.* already exists" 1>&2 + exit 1 + fi + TMPCC=".${NS}.${NAME}.tmp.cc" + TMPHPP=".${NS}.${NAME}.tmp.hpp" + sed "s/___FILEBASENAME___/${NAME}/g" Modules/templates/Module_tmp_in_NS.hpp.template > ${TMPHPP} + sed "s/___NAMESPACE___/${NS}/g" ${TMPHPP} > Modules/${NS}/${NAME}.hpp + rm -f ${TMPCC} ${TMPHPP} +fi +./make_module_list.sh diff --git a/extras/Hadrons/make_module_list.sh b/extras/Hadrons/make_module_list.sh new file mode 100755 index 00000000..ddc56ff6 --- /dev/null +++ b/extras/Hadrons/make_module_list.sh @@ -0,0 +1,12 @@ +#!/usr/bin/env bash + +echo 'modules_cc =\' > modules.inc +find Modules -name '*.cc' -type f -print | sed 's/^/ /;$q;s/$/ \\/' >> modules.inc +echo '' >> modules.inc +echo 'modules_hpp =\' >> modules.inc +find Modules -name '*.hpp' -type f -print | sed 's/^/ /;$q;s/$/ \\/' >> modules.inc +echo '' >> modules.inc +rm -f Modules.hpp +for f in `find Modules -name '*.hpp'`; do + echo "#include " >> Modules.hpp +done diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc new file mode 100644 index 00000000..4251ffa3 --- /dev/null +++ b/extras/Hadrons/modules.inc @@ -0,0 +1,19 @@ +modules_cc =\ + Modules/MGauge/Load.cc \ + Modules/MGauge/Random.cc \ + Modules/MGauge/Unit.cc + +modules_hpp =\ + Modules/MAction/DWF.hpp \ + Modules/MAction/Wilson.hpp \ + Modules/MContraction/Baryon.hpp \ + Modules/MContraction/Meson.hpp \ + Modules/MGauge/Load.hpp \ + Modules/MGauge/Random.hpp \ + Modules/MGauge/Unit.hpp \ + Modules/MSolver/RBPrecCG.hpp \ + Modules/MSource/Point.hpp \ + Modules/MSource/SeqGamma.hpp \ + Modules/MSource/Z2.hpp \ + Modules/Quark.hpp + diff --git a/extras/Makefile.am b/extras/Makefile.am new file mode 100644 index 00000000..d8c2b675 --- /dev/null +++ b/extras/Makefile.am @@ -0,0 +1 @@ +SUBDIRS = Hadrons \ No newline at end of file diff --git a/lib/AlignedAllocator.cc b/lib/AlignedAllocator.cc new file mode 100644 index 00000000..9df4ec1c --- /dev/null +++ b/lib/AlignedAllocator.cc @@ -0,0 +1,65 @@ + + + +#include + +namespace Grid { + +int PointerCache::victim; + + PointerCache::PointerCacheEntry PointerCache::Entries[PointerCache::Ncache]; + +void *PointerCache::Insert(void *ptr,size_t bytes) { + + if (bytes < 4096 ) return NULL; + +#ifdef _OPENMP + assert(omp_in_parallel()==0); +#endif + void * ret = NULL; + int v = -1; + + for(int e=0;e namespace Grid { + class PointerCache { + private: + + static const int Ncache=8; + static int victim; + + typedef struct { + void *address; + size_t bytes; + int valid; + } PointerCacheEntry; + + static PointerCacheEntry Entries[Ncache]; + + public: + + + static void *Insert(void *ptr,size_t bytes) ; + static void *Lookup(size_t bytes) ; + + }; + //////////////////////////////////////////////////////////////////// // A lattice of something, but assume the something is SIMDized. //////////////////////////////////////////////////////////////////// + template class alignedAllocator { public: @@ -66,27 +89,27 @@ public: pointer allocate(size_type __n, const void* _p= 0) { + size_type bytes = __n*sizeof(_Tp); + + _Tp *ptr = (_Tp *) PointerCache::Lookup(bytes); + #ifdef HAVE_MM_MALLOC_H - _Tp * ptr = (_Tp *) _mm_malloc(__n*sizeof(_Tp),128); + if ( ptr == (_Tp *) NULL ) ptr = (_Tp *) _mm_malloc(bytes,128); #else - _Tp * ptr = (_Tp *) memalign(128,__n*sizeof(_Tp)); + if ( ptr == (_Tp *) NULL ) ptr = (_Tp *) memalign(128,bytes); #endif - _Tp tmp; -#ifdef GRID_NUMA -#pragma omp parallel for schedule(static) - for(int i=0;i<__n;i++){ - ptr[i]=tmp; - } -#endif return ptr; } - void deallocate(pointer __p, size_type) { + void deallocate(pointer __p, size_type __n) { + size_type bytes = __n * sizeof(_Tp); + pointer __freeme = (pointer)PointerCache::Insert((void *)__p,bytes); + #ifdef HAVE_MM_MALLOC_H - _mm_free((void *)__p); + if ( __freeme ) _mm_free((void *)__freeme); #else - free((void *)__p); + if ( __freeme ) free((void *)__freeme); #endif } void construct(pointer __p, const _Tp& __val) { }; diff --git a/lib/Grid.h b/lib/Grid.h index 718b2135..c2f6e00f 100644 --- a/lib/Grid.h +++ b/lib/Grid.h @@ -59,13 +59,13 @@ Author: paboyle /////////////////// // Grid headers /////////////////// -#include #include "Config.h" #include #include #include #include #include +#include #include #include #include diff --git a/lib/Hadrons b/lib/Hadrons new file mode 120000 index 00000000..1f422592 --- /dev/null +++ b/lib/Hadrons @@ -0,0 +1 @@ +../extras/Hadrons \ No newline at end of file diff --git a/lib/Log.h b/lib/Log.h index d7422ca9..8dce0acd 100644 --- a/lib/Log.h +++ b/lib/Log.h @@ -110,8 +110,8 @@ public: friend std::ostream& operator<< (std::ostream& stream, Logger& log){ if ( log.active ) { - stream << log.background()<< log.topName << log.background()<< " : "; - stream << log.colour() < >& table,const Lattice { PARALLEL_FOR_LOOP for(int i=0;i #endif #define PARALLEL_NESTED_LOOP2 _Pragma("omp parallel for collapse(2)") #define PARALLEL_REGION _Pragma("omp parallel") +#define PARALLEL_CRITICAL _Pragma("omp critical") #else #define PARALLEL_FOR_LOOP #define PARALLEL_FOR_LOOP_INTERN #define PARALLEL_NESTED_LOOP2 #define PARALLEL_REGION +#define PARALLEL_CRITICAL #endif namespace Grid { diff --git a/lib/lattice/Lattice_transfer.h b/lib/lattice/Lattice_transfer.h index cc4617de..a49b1b5f 100644 --- a/lib/lattice/Lattice_transfer.h +++ b/lib/lattice/Lattice_transfer.h @@ -386,7 +386,7 @@ void InsertSlice(Lattice &lowDim,Lattice & higherDim,int slice, int } // the above should guarantee that the operations are local - //PARALLEL_FOR_LOOP + PARALLEL_FOR_LOOP for(int idx=0;idxlSites();idx++){ std::vector lcoor(nl); std::vector hcoor(nh); @@ -428,7 +428,7 @@ void ExtractSlice(Lattice &lowDim, Lattice & higherDim,int slice, in } } // the above should guarantee that the operations are local - //PARALLEL_FOR_LOOP + PARALLEL_FOR_LOOP for(int idx=0;idxlSites();idx++){ std::vector lcoor(nl); std::vector hcoor(nh); diff --git a/lib/qcd/action/fermion/CayleyFermion5D.cc b/lib/qcd/action/fermion/CayleyFermion5D.cc index b8e98dce..781380e5 100644 --- a/lib/qcd/action/fermion/CayleyFermion5D.cc +++ b/lib/qcd/action/fermion/CayleyFermion5D.cc @@ -29,6 +29,7 @@ Author: paboyle *************************************************************************************/ /* END LEGAL */ +#include #include @@ -48,18 +49,18 @@ namespace QCD { FourDimGrid, FourDimRedBlackGrid,_M5,p), mass(_mass) - { } + { + } template void CayleyFermion5D::Dminus(const FermionField &psi, FermionField &chi) { int Ls=this->Ls; - FermionField tmp(psi._grid); - this->DW(psi,tmp,DaggerNo); + this->DW(psi,this->tmp(),DaggerNo); for(int s=0;stmp(),s,s);// chi = (1-c[s] D_W) psi } } @@ -87,8 +88,8 @@ template void CayleyFermion5D::CayleyReport(void) std::cout << GridLogMessage << "CayleyFermion5D Number of MooeeInv Calls : " << MooeeInvCalls << std::endl; std::cout << GridLogMessage << "CayleyFermion5D ComputeTime/Calls : " << MooeeInvTime / MooeeInvCalls << " us" << std::endl; - // Flops = 9*12*Ls*vol/2 - RealD mflops = 9.0*12*volume*MooeeInvCalls/MooeeInvTime/2; // 2 for red black counting + // Flops = MADD * Ls *Ls *4dvol * spin/colour/complex + RealD mflops = 2.0*24*this->Ls*volume*MooeeInvCalls/MooeeInvTime/2; // 2 for red black counting std::cout << GridLogMessage << "Average mflops/s per call : " << mflops << std::endl; std::cout << GridLogMessage << "Average mflops/s per call per rank : " << mflops/NP << std::endl; } @@ -110,12 +111,11 @@ template void CayleyFermion5D::DminusDag(const FermionField &psi, FermionField &chi) { int Ls=this->Ls; - FermionField tmp(psi._grid); - this->DW(psi,tmp,DaggerYes); + this->DW(psi,this->tmp(),DaggerYes); for(int s=0;stmp(),s,s);// chi = (1-c[s] D_W) psi } } template @@ -138,6 +138,7 @@ void CayleyFermion5D::Meooe5D (const FermionField &psi, FermionField &D lower[0] =-mass*lower[0]; M5D(psi,psi,Din,lower,diag,upper); } +// FIXME Redunant with the above routine; check this and eliminate template void CayleyFermion5D::Meo5D (const FermionField &psi, FermionField &chi) { int Ls=this->Ls; @@ -259,36 +260,33 @@ template void CayleyFermion5D::Meooe (const FermionField &psi, FermionField &chi) { int Ls=this->Ls; - FermionField tmp(psi._grid); - Meooe5D(psi,tmp); + Meooe5D(psi,this->tmp()); if ( psi.checkerboard == Odd ) { - this->DhopEO(tmp,chi,DaggerNo); + this->DhopEO(this->tmp(),chi,DaggerNo); } else { - this->DhopOE(tmp,chi,DaggerNo); + this->DhopOE(this->tmp(),chi,DaggerNo); } } template void CayleyFermion5D::MeooeDag (const FermionField &psi, FermionField &chi) { - FermionField tmp(psi._grid); // Apply 4d dslash if ( psi.checkerboard == Odd ) { - this->DhopEO(psi,tmp,DaggerYes); + this->DhopEO(psi,this->tmp(),DaggerYes); } else { - this->DhopOE(psi,tmp,DaggerYes); + this->DhopOE(psi,this->tmp(),DaggerYes); } - MeooeDag5D(tmp,chi); + MeooeDag5D(this->tmp(),chi); } template void CayleyFermion5D::Mdir (const FermionField &psi, FermionField &chi,int dir,int disp){ - FermionField tmp(psi._grid); - Meo5D(psi,tmp); + Meo5D(psi,this->tmp()); // Apply 4d dslash fragment - this->DhopDir(tmp,chi,dir,disp); + this->DhopDir(this->tmp(),chi,dir,disp); } // force terms; five routines; default to Dhop on diagonal template @@ -459,9 +457,91 @@ void CayleyFermion5D::SetCoefficientsInternal(RealD zolo_hi,std::vectorMooeeInternalCompute(0,inv,MatpInv,MatmInv); + this->MooeeInternalCompute(1,inv,MatpInvDag,MatmInvDag); + } +template +void CayleyFermion5D::MooeeInternalCompute(int dag, int inv, + Vector > & Matp, + Vector > & Matm) +{ + int Ls=this->Ls; + + GridBase *grid = this->FermionRedBlackGrid(); + int LLs = grid->_rdimensions[0]; + + if ( LLs == Ls ) return; // Not vectorised in 5th direction + + Eigen::MatrixXcd Pplus = Eigen::MatrixXcd::Zero(Ls,Ls); + Eigen::MatrixXcd Pminus = Eigen::MatrixXcd::Zero(Ls,Ls); + + for(int s=0;s::iscomplex() ) { + sp[l] = PplusMat (l*istride+s1*ostride,s2); + sm[l] = PminusMat(l*istride+s1*ostride,s2); + } else { + // if real + scalar_type tmp; + tmp = PplusMat (l*istride+s1*ostride,s2); + sp[l] = scalar_type(tmp.real(),tmp.real()); + tmp = PminusMat(l*istride+s1*ostride,s2); + sm[l] = scalar_type(tmp.real(),tmp.real()); + } + } + Matp[LLs*s2+s1] = Vp; + Matm[LLs*s2+s1] = Vm; + }} +} + FermOpTemplateInstantiate(CayleyFermion5D); GparityFermOpTemplateInstantiate(CayleyFermion5D); diff --git a/lib/qcd/action/fermion/CayleyFermion5D.h b/lib/qcd/action/fermion/CayleyFermion5D.h index 6fb58234..86255be6 100644 --- a/lib/qcd/action/fermion/CayleyFermion5D.h +++ b/lib/qcd/action/fermion/CayleyFermion5D.h @@ -33,6 +33,31 @@ namespace Grid { namespace QCD { + template struct switcheroo { + static inline int iscomplex() { return 0; } + + template + static inline vec mult(vec a, vec b) { + return real_mult(a,b); + } + }; + template<> struct switcheroo { + static inline int iscomplex() { return 1; } + + template + static inline vec mult(vec a, vec b) { + return a*b; + } + }; + template<> struct switcheroo { + static inline int iscomplex() { return 1; } + template + static inline vec mult(vec a, vec b) { + return a*b; + } + }; + + template class CayleyFermion5D : public WilsonFermion5D { @@ -75,7 +100,19 @@ namespace Grid { std::vector &lower, std::vector &diag, std::vector &upper); + void MooeeInternal(const FermionField &in, FermionField &out,int dag,int inv); + void MooeeInternalCompute(int dag, int inv, Vector > & Matp, Vector > & Matm); + + void MooeeInternalAsm(const FermionField &in, FermionField &out, + int LLs, int site, + Vector > &Matp, + Vector > &Matm); + void MooeeInternalZAsm(const FermionField &in, FermionField &out, + int LLs, int site, + Vector > &Matp, + Vector > &Matm); + virtual void Instantiatable(void)=0; @@ -112,6 +149,12 @@ namespace Grid { std::vector ueem; std::vector dee; + // Matrices of 5d ee inverse params + Vector > MatpInv; + Vector > MatmInv; + Vector > MatpInvDag; + Vector > MatmInvDag; + // Constructors CayleyFermion5D(GaugeField &_Umu, GridCartesian &FiveDimGrid, diff --git a/lib/qcd/action/fermion/CayleyFermion5Dvec.cc b/lib/qcd/action/fermion/CayleyFermion5Dvec.cc index 35a10de2..ed742ea3 100644 --- a/lib/qcd/action/fermion/CayleyFermion5Dvec.cc +++ b/lib/qcd/action/fermion/CayleyFermion5Dvec.cc @@ -29,13 +29,12 @@ Author: paboyle *************************************************************************************/ /* END LEGAL */ -#include + #include namespace Grid { -namespace QCD { - /* +namespace QCD { /* * Dense matrix versions of routines */ template @@ -126,7 +125,6 @@ PARALLEL_FOR_LOOP for(int v=0;v(hp_00.v); hp_01.v = Optimization::Rotate::tRotate<2>(hp_01.v); @@ -165,42 +160,20 @@ PARALLEL_FOR_LOOP hm_12.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_12.v); } - /* - if ( ss==0) std::cout << " dphi_00 " <::mult(d[v]()()(), phi[ss+v]()(0)(0)) + switcheroo::mult(l[v]()()(),hm_00); + Simd p_01 = switcheroo::mult(d[v]()()(), phi[ss+v]()(0)(1)) + switcheroo::mult(l[v]()()(),hm_01); + Simd p_02 = switcheroo::mult(d[v]()()(), phi[ss+v]()(0)(2)) + switcheroo::mult(l[v]()()(),hm_02); + Simd p_10 = switcheroo::mult(d[v]()()(), phi[ss+v]()(1)(0)) + switcheroo::mult(l[v]()()(),hm_10); + Simd p_11 = switcheroo::mult(d[v]()()(), phi[ss+v]()(1)(1)) + switcheroo::mult(l[v]()()(),hm_11); + Simd p_12 = switcheroo::mult(d[v]()()(), phi[ss+v]()(1)(2)) + switcheroo::mult(l[v]()()(),hm_12); + Simd p_20 = switcheroo::mult(d[v]()()(), phi[ss+v]()(2)(0)) + switcheroo::mult(u[v]()()(),hp_00); + Simd p_21 = switcheroo::mult(d[v]()()(), phi[ss+v]()(2)(1)) + switcheroo::mult(u[v]()()(),hp_01); + Simd p_22 = switcheroo::mult(d[v]()()(), phi[ss+v]()(2)(2)) + switcheroo::mult(u[v]()()(),hp_02); + Simd p_30 = switcheroo::mult(d[v]()()(), phi[ss+v]()(3)(0)) + switcheroo::mult(u[v]()()(),hp_10); + Simd p_31 = switcheroo::mult(d[v]()()(), phi[ss+v]()(3)(1)) + switcheroo::mult(u[v]()()(),hp_11); + Simd p_32 = switcheroo::mult(d[v]()()(), phi[ss+v]()(3)(2)) + switcheroo::mult(u[v]()()(),hp_12); - - // if ( ss==0){ - /* - std::cout << ss<<" "<< v<< " good "<< chi[ss+v]()(0)(0) << " bad "<::M5Ddag(const FermionField &psi, M5Dtime-=usecond(); PARALLEL_FOR_LOOP for(int ss=0;ssoSites();ss+=LLs){ // adds LLs - +#if 0 alignas(64) SiteHalfSpinor hp; alignas(64) SiteHalfSpinor hm; alignas(64) SiteSpinor fp; @@ -287,9 +260,504 @@ PARALLEL_FOR_LOOP chi[ss+v] = chi[ss+v] +l[v]*fm; } +#else + for(int v=0;v(hp_00.v); + hp_01.v = Optimization::Rotate::tRotate<2>(hp_01.v); + hp_02.v = Optimization::Rotate::tRotate<2>(hp_02.v); + hp_10.v = Optimization::Rotate::tRotate<2>(hp_10.v); + hp_11.v = Optimization::Rotate::tRotate<2>(hp_11.v); + hp_12.v = Optimization::Rotate::tRotate<2>(hp_12.v); + } + if ( vm>=v ) { + hm_00.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_00.v); + hm_01.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_01.v); + hm_02.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_02.v); + hm_10.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_10.v); + hm_11.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_11.v); + hm_12.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_12.v); + } + + Simd p_00 = switcheroo::mult(d[v]()()(), phi[ss+v]()(0)(0)) + switcheroo::mult(u[v]()()(),hp_00); + Simd p_01 = switcheroo::mult(d[v]()()(), phi[ss+v]()(0)(1)) + switcheroo::mult(u[v]()()(),hp_01); + Simd p_02 = switcheroo::mult(d[v]()()(), phi[ss+v]()(0)(2)) + switcheroo::mult(u[v]()()(),hp_02); + Simd p_10 = switcheroo::mult(d[v]()()(), phi[ss+v]()(1)(0)) + switcheroo::mult(u[v]()()(),hp_10); + Simd p_11 = switcheroo::mult(d[v]()()(), phi[ss+v]()(1)(1)) + switcheroo::mult(u[v]()()(),hp_11); + Simd p_12 = switcheroo::mult(d[v]()()(), phi[ss+v]()(1)(2)) + switcheroo::mult(u[v]()()(),hp_12); + + Simd p_20 = switcheroo::mult(d[v]()()(), phi[ss+v]()(2)(0)) + switcheroo::mult(l[v]()()(),hm_00); + Simd p_21 = switcheroo::mult(d[v]()()(), phi[ss+v]()(2)(1)) + switcheroo::mult(l[v]()()(),hm_01); + Simd p_22 = switcheroo::mult(d[v]()()(), phi[ss+v]()(2)(2)) + switcheroo::mult(l[v]()()(),hm_02); + Simd p_30 = switcheroo::mult(d[v]()()(), phi[ss+v]()(3)(0)) + switcheroo::mult(l[v]()()(),hm_10); + Simd p_31 = switcheroo::mult(d[v]()()(), phi[ss+v]()(3)(1)) + switcheroo::mult(l[v]()()(),hm_11); + Simd p_32 = switcheroo::mult(d[v]()()(), phi[ss+v]()(3)(2)) + switcheroo::mult(l[v]()()(),hm_12); + + vstream(chi[ss+v]()(0)(0),p_00); + vstream(chi[ss+v]()(0)(1),p_01); + vstream(chi[ss+v]()(0)(2),p_02); + vstream(chi[ss+v]()(1)(0),p_10); + vstream(chi[ss+v]()(1)(1),p_11); + vstream(chi[ss+v]()(1)(2),p_12); + vstream(chi[ss+v]()(2)(0),p_20); + vstream(chi[ss+v]()(2)(1),p_21); + vstream(chi[ss+v]()(2)(2),p_22); + vstream(chi[ss+v]()(3)(0),p_30); + vstream(chi[ss+v]()(3)(1),p_31); + vstream(chi[ss+v]()(3)(2),p_32); + } +#endif } M5Dtime+=usecond(); } + + +#ifdef AVX512 +#include +#include +#include +#endif + +template +void CayleyFermion5D::MooeeInternalAsm(const FermionField &psi, FermionField &chi, + int LLs, int site, + Vector > &Matp, + Vector > &Matm) +{ +#ifndef AVX512 + { + SiteHalfSpinor BcastP; + SiteHalfSpinor BcastM; + SiteHalfSpinor SiteChiP; + SiteHalfSpinor SiteChiM; + + // Ls*Ls * 2 * 12 * vol flops + for(int s1=0;s1); + for(int s1=0;s1 +void CayleyFermion5D::MooeeInternalZAsm(const FermionField &psi, FermionField &chi, + int LLs, int site, Vector > &Matp, Vector > &Matm) +{ +#ifndef AVX512 + { + SiteHalfSpinor BcastP; + SiteHalfSpinor BcastM; + SiteHalfSpinor SiteChiP; + SiteHalfSpinor SiteChiM; + + // Ls*Ls * 2 * 12 * vol flops + for(int s1=0;s1); + for(int s1=0;s1 void CayleyFermion5D::MooeeInternal(const FermionField &psi, FermionField &chi,int dag, int inv) { @@ -299,108 +767,41 @@ void CayleyFermion5D::MooeeInternal(const FermionField &psi, FermionField chi.checkerboard=psi.checkerboard; - Eigen::MatrixXcd Pplus = Eigen::MatrixXcd::Zero(Ls,Ls); - Eigen::MatrixXcd Pminus = Eigen::MatrixXcd::Zero(Ls,Ls); + Vector > Matp; + Vector > Matm; + Vector > *_Matp; + Vector > *_Matm; - for(int s=0;s > Matp(Ls*LLs); - Vector > Matm(Ls*LLs); + assert(_Matp->size()==Ls*LLs); - for(int s2=0;s2 SitePplus(LLs); - Vector SitePminus(LLs); - Vector SiteChiP(LLs); - Vector SiteChiM(LLs); - Vector SiteChi(LLs); - - SiteHalfSpinor BcastP; - SiteHalfSpinor BcastM; - -#pragma omp for - for(auto site=0;site::iscomplex() ) { + PARALLEL_FOR_LOOP + for(auto site=0;site::MooeeInternal(const Fermion template void CayleyFermion5D::MooeeInternal(const FermionField &psi, FermionField &chi,int dag, int inv); template void CayleyFermion5D::MooeeInternal(const FermionField &psi, FermionField &chi,int dag, int inv); + }} diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h index 742c6e08..676a0e83 100644 --- a/lib/qcd/action/fermion/FermionOperator.h +++ b/lib/qcd/action/fermion/FermionOperator.h @@ -48,6 +48,8 @@ namespace Grid { FermionOperator(const ImplParams &p= ImplParams()) : Impl(p) {}; + virtual FermionField &tmp(void) = 0; + GridBase * Grid(void) { return FermionGrid(); }; // this is all the linalg routines need to know GridBase * RedBlackGrid(void) { return FermionRedBlackGrid(); }; diff --git a/lib/qcd/action/fermion/FermionOperatorImpl.h b/lib/qcd/action/fermion/FermionOperatorImpl.h index f658872d..46837d0c 100644 --- a/lib/qcd/action/fermion/FermionOperatorImpl.h +++ b/lib/qcd/action/fermion/FermionOperatorImpl.h @@ -48,10 +48,12 @@ namespace QCD { // typedef typename XXX GaugeField; // typedef typename XXX GaugeActField; // typedef typename XXX FermionField; + // typedef typename XXX PropagatorField; // typedef typename XXX DoubledGaugeField; // typedef typename XXX SiteSpinor; - // typedef typename XXX SiteHalfSpinor; - // typedef typename XXX Compressor; + // typedef typename XXX SitePropagator; + // typedef typename XXX SiteHalfSpinor; + // typedef typename XXX Compressor; // // and Methods: // void ImportGauge(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu) @@ -94,14 +96,16 @@ namespace QCD { //////////////////////////////////////////////////////////////////////// #define INHERIT_FIMPL_TYPES(Impl)\ - typedef typename Impl::FermionField FermionField; \ - typedef typename Impl::DoubledGaugeField DoubledGaugeField; \ - typedef typename Impl::SiteSpinor SiteSpinor; \ - typedef typename Impl::SiteHalfSpinor SiteHalfSpinor; \ - typedef typename Impl::Compressor Compressor; \ - typedef typename Impl::StencilImpl StencilImpl; \ - typedef typename Impl::ImplParams ImplParams; \ - typedef typename Impl::Coeff_t Coeff_t; + typedef typename Impl::FermionField FermionField; \ + typedef typename Impl::PropagatorField PropagatorField; \ + typedef typename Impl::DoubledGaugeField DoubledGaugeField; \ + typedef typename Impl::SiteSpinor SiteSpinor; \ + typedef typename Impl::SitePropagator SitePropagator; \ + typedef typename Impl::SiteHalfSpinor SiteHalfSpinor; \ + typedef typename Impl::Compressor Compressor; \ + typedef typename Impl::StencilImpl StencilImpl; \ + typedef typename Impl::ImplParams ImplParams; \ + typedef typename Impl::Coeff_t Coeff_t; \ #define INHERIT_IMPL_TYPES(Base) \ INHERIT_GIMPL_TYPES(Base) \ @@ -127,14 +131,17 @@ namespace QCD { INHERIT_GIMPL_TYPES(Gimpl); template using iImplSpinor = iScalar, Ns> >; + template using iImplPropagator = iScalar, Ns> >; template using iImplHalfSpinor = iScalar, Nhs> >; template using iImplDoubledGaugeField = iVector >, Nds>; typedef iImplSpinor SiteSpinor; + typedef iImplPropagator SitePropagator; typedef iImplHalfSpinor SiteHalfSpinor; typedef iImplDoubledGaugeField SiteDoubledGaugeField; typedef Lattice FermionField; + typedef Lattice PropagatorField; typedef Lattice DoubledGaugeField; typedef WilsonCompressor Compressor; @@ -216,14 +223,17 @@ class DomainWallVec5dImpl : public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres INHERIT_GIMPL_TYPES(Gimpl); template using iImplSpinor = iScalar, Ns> >; + template using iImplPropagator = iScalar, Ns> >; template using iImplHalfSpinor = iScalar, Nhs> >; template using iImplDoubledGaugeField = iVector >, Nds>; template using iImplGaugeField = iVector >, Nd>; template using iImplGaugeLink = iScalar > >; typedef iImplSpinor SiteSpinor; + typedef iImplPropagator SitePropagator; typedef iImplHalfSpinor SiteHalfSpinor; typedef Lattice FermionField; + typedef Lattice PropagatorField; // Make the doubled gauge field a *scalar* typedef iImplDoubledGaugeField SiteDoubledGaugeField; // This is a scalar @@ -352,14 +362,17 @@ class GparityWilsonImpl : public ConjugateGaugeImpl using iImplSpinor = iVector, Ns>, Ngp>; + template using iImplPropagator = iVector, Ns>, Ngp >; template using iImplHalfSpinor = iVector, Nhs>, Ngp>; template using iImplDoubledGaugeField = iVector >, Nds>, Ngp>; typedef iImplSpinor SiteSpinor; + typedef iImplPropagator SitePropagator; typedef iImplHalfSpinor SiteHalfSpinor; typedef iImplDoubledGaugeField SiteDoubledGaugeField; typedef Lattice FermionField; + typedef Lattice PropagatorField; typedef Lattice DoubledGaugeField; typedef WilsonCompressor Compressor; diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index e75584bc..2811e5c7 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -61,7 +61,9 @@ WilsonFermion::WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid, LebesgueEvenOdd(_cbgrid), Umu(&Fgrid), UmuEven(&Hgrid), - UmuOdd(&Hgrid) { + UmuOdd(&Hgrid), + _tmp(&Hgrid) +{ // Allocate the required comms buffer ImportGauge(_Umu); } diff --git a/lib/qcd/action/fermion/WilsonFermion.h b/lib/qcd/action/fermion/WilsonFermion.h index 40fbd1bf..933be732 100644 --- a/lib/qcd/action/fermion/WilsonFermion.h +++ b/lib/qcd/action/fermion/WilsonFermion.h @@ -58,6 +58,9 @@ class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { GridBase *FermionGrid(void) { return _grid; } GridBase *FermionRedBlackGrid(void) { return _cbgrid; } + FermionField _tmp; + FermionField &tmp(void) { return _tmp; } + ////////////////////////////////////////////////////////////////// // override multiply; cut number routines if pass dagger argument // and also make interface more uniformly consistent diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index d930e87c..16594e0f 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -61,7 +61,8 @@ WilsonFermion5D::WilsonFermion5D(GaugeField &_Umu, UmuEven(_FourDimRedBlackGrid), UmuOdd (_FourDimRedBlackGrid), Lebesgue(_FourDimGrid), - LebesgueEvenOdd(_FourDimRedBlackGrid) + LebesgueEvenOdd(_FourDimRedBlackGrid), + _tmp(&FiveDimRedBlackGrid) { if (Impl::LsVectorised) { diff --git a/lib/qcd/action/fermion/WilsonFermion5D.h b/lib/qcd/action/fermion/WilsonFermion5D.h index ffb5c58e..fb4fa925 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.h +++ b/lib/qcd/action/fermion/WilsonFermion5D.h @@ -74,6 +74,9 @@ namespace QCD { typedef WilsonKernels Kernels; PmuStat stat; + FermionField _tmp; + FermionField &tmp(void) { return _tmp; } + void Report(void); void ZeroCounters(void); double DhopCalls; diff --git a/lib/serialisation/BaseIO.h b/lib/serialisation/BaseIO.h index 01b21784..887f7a05 100644 --- a/lib/serialisation/BaseIO.h +++ b/lib/serialisation/BaseIO.h @@ -32,6 +32,7 @@ Author: Peter Boyle #include namespace Grid { + // Vector IO utilities /////////////////////////////////////////////////////// // helper function to read space-separated values template std::vector strToVec(const std::string s) @@ -67,12 +68,78 @@ namespace Grid { return os; } - class Serializable {}; + // Vector element trait ////////////////////////////////////////////////////// + template + struct element + { + typedef T type; + static constexpr bool is_number = false; + }; - - - + template + struct element> + { + typedef typename element::type type; + static constexpr bool is_number = std::is_arithmetic::value + or is_complex::value + or element::is_number; + }; + + // Vector flatening utility class //////////////////////////////////////////// + // Class to flatten a multidimensional std::vector + template + class Flatten + { + public: + typedef typename element::type Element; + public: + explicit Flatten(const V &vector); + const V & getVector(void); + const std::vector & getFlatVector(void); + const std::vector & getDim(void); + private: + void accumulate(const Element &e); + template + void accumulate(const W &v); + void accumulateDim(const Element &e); + template + void accumulateDim(const W &v); + private: + const V &vector_; + std::vector flatVector_; + std::vector dim_; + }; + + // Class to reconstruct a multidimensional std::vector + template + class Reconstruct + { + public: + typedef typename element::type Element; + public: + Reconstruct(const std::vector &flatVector, + const std::vector &dim); + const V & getVector(void); + const std::vector & getFlatVector(void); + const std::vector & getDim(void); + private: + void fill(std::vector &v); + template + void fill(W &v); + void resize(std::vector &v, const unsigned int dim); + template + void resize(W &v, const unsigned int dim); + private: + V vector_; + const std::vector &flatVector_; + std::vector dim_; + size_t ind_{0}; + unsigned int dimInd_{0}; + }; + + // Abstract writer/reader classes //////////////////////////////////////////// // static polymorphism implemented using CRTP idiom + class Serializable; // Static abstract writer template @@ -87,12 +154,7 @@ namespace Grid { typename std::enable_if::value, void>::type write(const std::string& s, const U &output); template - typename std::enable_if::value, void>::type - write(const std::string& s, const U &output); - template - typename std::enable_if< - !(std::is_base_of::value or std::is_enum::value), - void>::type + typename std::enable_if::value, void>::type write(const std::string& s, const U &output); private: T *upcast; @@ -111,12 +173,7 @@ namespace Grid { typename std::enable_if::value, void>::type read(const std::string& s, U &output); template - typename std::enable_if::value, void>::type - read(const std::string& s, U &output); - template - typename std::enable_if< - !(std::is_base_of::value or std::is_enum::value), - void>::type + typename std::enable_if::value, void>::type read(const std::string& s, U &output); protected: template @@ -125,16 +182,151 @@ namespace Grid { T *upcast; }; - // type traits - // What is the vtype - template struct isReader { - static const bool value = false; - }; - template struct isWriter { - static const bool value = false; - }; - // Generic writer interface + // serializable base class + class Serializable + { + public: + template + static inline void write(Writer &WR,const std::string &s, + const Serializable &obj) + {} + + template + static inline void read(Reader &RD,const std::string &s, + Serializable &obj) + {} + + friend inline std::ostream & operator<<(std::ostream &os, + const Serializable &obj) + { + return os; + } + }; + + // Flatten class template implementation ///////////////////////////////////// + template + void Flatten::accumulate(const Element &e) + { + flatVector_.push_back(e); + } + + template + template + void Flatten::accumulate(const W &v) + { + for (auto &e: v) + { + accumulate(e); + } + } + + template + void Flatten::accumulateDim(const Element &e) {}; + + template + template + void Flatten::accumulateDim(const W &v) + { + dim_.push_back(v.size()); + accumulateDim(v[0]); + } + + template + Flatten::Flatten(const V &vector) + : vector_(vector) + { + accumulate(vector_); + accumulateDim(vector_); + } + + template + const V & Flatten::getVector(void) + { + return vector_; + } + + template + const std::vector::Element> & + Flatten::getFlatVector(void) + { + return flatVector_; + } + + template + const std::vector & Flatten::getDim(void) + { + return dim_; + } + + // Reconstruct class template implementation ///////////////////////////////// + template + void Reconstruct::fill(std::vector &v) + { + for (auto &e: v) + { + e = flatVector_[ind_++]; + } + } + + template + template + void Reconstruct::fill(W &v) + { + for (auto &e: v) + { + fill(e); + } + } + + template + void Reconstruct::resize(std::vector &v, const unsigned int dim) + { + v.resize(dim_[dim]); + } + + template + template + void Reconstruct::resize(W &v, const unsigned int dim) + { + v.resize(dim_[dim]); + for (auto &e: v) + { + resize(e, dim + 1); + } + } + + template + Reconstruct::Reconstruct(const std::vector &flatVector, + const std::vector &dim) + : flatVector_(flatVector) + , dim_(dim) + { + resize(vector_, 0); + fill(vector_); + } + + template + const V & Reconstruct::getVector(void) + { + return vector_; + } + + template + const std::vector::Element> & + Reconstruct::getFlatVector(void) + { + return flatVector_; + } + + template + const std::vector & Reconstruct::getDim(void) + { + return dim_; + } + + // Generic writer interface ////////////////////////////////////////////////// +>>>>>>> develop template inline void push(Writer &w, const std::string &s) { w.push(s); @@ -212,23 +404,13 @@ namespace Grid { template template - typename std::enable_if::value, void>::type - Writer::write(const std::string &s, const U &output) - { - EnumIO::write(*this, s, output); - } - - template - template - typename std::enable_if< - !(std::is_base_of::value or std::is_enum::value), - void>::type + typename std::enable_if::value, void>::type Writer::write(const std::string &s, const U &output) { upcast->writeDefault(s, output); } - // Reader template implementation //////////////////////////////////////////// + // Reader template implementation template Reader::Reader(void) { @@ -257,17 +439,7 @@ namespace Grid { template template - typename std::enable_if::value, void>::type - Reader::read(const std::string &s, U &output) - { - EnumIO::read(*this, s, output); - } - - template - template - typename std::enable_if< - !(std::is_base_of::value or std::is_enum::value), - void>::type + typename std::enable_if::value, void>::type Reader::read(const std::string &s, U &output) { upcast->readDefault(s, output); @@ -291,7 +463,6 @@ namespace Grid { abort(); } } - } #endif diff --git a/lib/serialisation/Hdf5IO.cc b/lib/serialisation/Hdf5IO.cc new file mode 100644 index 00000000..c5313495 --- /dev/null +++ b/lib/serialisation/Hdf5IO.cc @@ -0,0 +1,103 @@ +#include + +using namespace Grid; +#ifndef H5_NO_NAMESPACE +using namespace H5NS; +#endif + +// Writer implementation /////////////////////////////////////////////////////// +Hdf5Writer::Hdf5Writer(const std::string &fileName) +: fileName_(fileName) +, file_(fileName.c_str(), H5F_ACC_TRUNC) +{ + group_ = file_.openGroup("/"); + writeSingleAttribute(dataSetThres_, HDF5_GRID_GUARD "dataset_threshold", + Hdf5Type::type()); +} + +void Hdf5Writer::push(const std::string &s) +{ + group_ = group_.createGroup(s); + path_.push_back(s); +} + +void Hdf5Writer::pop(void) +{ + path_.pop_back(); + if (path_.empty()) + { + group_ = file_.openGroup("/"); + } + else + { + auto binOp = [](const std::string &a, const std::string &b)->std::string + { + return a + "/" + b; + }; + + group_ = group_.openGroup(std::accumulate(path_.begin(), path_.end(), + std::string(""), binOp)); + } +} + +template <> +void Hdf5Writer::writeDefault(const std::string &s, const std::string &x) +{ + StrType strType(PredType::C_S1, x.size()); + + writeSingleAttribute(*(x.data()), s, strType); +} + +void Hdf5Writer::writeDefault(const std::string &s, const char *x) +{ + std::string sx(x); + + writeDefault(s, sx); +} + +// Reader implementation /////////////////////////////////////////////////////// +Hdf5Reader::Hdf5Reader(const std::string &fileName) +: fileName_(fileName) +, file_(fileName.c_str(), H5F_ACC_RDONLY) +{ + group_ = file_.openGroup("/"); + readSingleAttribute(dataSetThres_, HDF5_GRID_GUARD "dataset_threshold", + Hdf5Type::type()); +} + +void Hdf5Reader::push(const std::string &s) +{ + group_ = group_.openGroup(s); + path_.push_back(s); +} + +void Hdf5Reader::pop(void) +{ + path_.pop_back(); + if (path_.empty()) + { + group_ = file_.openGroup("/"); + } + else + { + auto binOp = [](const std::string &a, const std::string &b)->std::string + { + return a + "/" + b; + }; + + group_ = group_.openGroup(std::accumulate(path_.begin(), path_.end(), + std::string(""), binOp)); + } +} + +template <> +void Hdf5Reader::readDefault(const std::string &s, std::string &x) +{ + Attribute attribute; + + attribute = group_.openAttribute(s); + StrType strType = attribute.getStrType(); + + x.resize(strType.getSize()); + attribute.read(strType, &(x[0])); +} diff --git a/lib/serialisation/Hdf5IO.h b/lib/serialisation/Hdf5IO.h new file mode 100644 index 00000000..2f891cd4 --- /dev/null +++ b/lib/serialisation/Hdf5IO.h @@ -0,0 +1,242 @@ +#ifndef GRID_SERIALISATION_HDF5_H +#define GRID_SERIALISATION_HDF5_H + +#include +#include +#include +#include +#include "Hdf5Type.h" + +#ifndef H5_NO_NAMESPACE +#define H5NS H5 +#endif + +// default thresold above which datasets are used instead of attributes +#ifndef HDF5_DEF_DATASET_THRES +#define HDF5_DEF_DATASET_THRES 6u +#endif + +// name guard for Grid metadata +#define HDF5_GRID_GUARD "_Grid_" + +namespace Grid +{ + class Hdf5Writer: public Writer + { + public: + Hdf5Writer(const std::string &fileName); + virtual ~Hdf5Writer(void) = default; + void push(const std::string &s); + void pop(void); + void writeDefault(const std::string &s, const char *x); + template + void writeDefault(const std::string &s, const U &x); + template + typename std::enable_if>::is_number, void>::type + writeDefault(const std::string &s, const std::vector &x); + template + typename std::enable_if>::is_number, void>::type + writeDefault(const std::string &s, const std::vector &x); + private: + template + void writeSingleAttribute(const U &x, const std::string &name, + const H5NS::DataType &type); + private: + std::string fileName_; + std::vector path_; + H5NS::H5File file_; + H5NS::Group group_; + unsigned int dataSetThres_{HDF5_DEF_DATASET_THRES}; + }; + + class Hdf5Reader: public Reader + { + public: + Hdf5Reader(const std::string &fileName); + virtual ~Hdf5Reader(void) = default; + void push(const std::string &s); + void pop(void); + template + void readDefault(const std::string &s, U &output); + template + typename std::enable_if>::is_number, void>::type + readDefault(const std::string &s, std::vector &x); + template + typename std::enable_if>::is_number, void>::type + readDefault(const std::string &s, std::vector &x); + private: + template + void readSingleAttribute(U &x, const std::string &name, + const H5NS::DataType &type); + private: + std::string fileName_; + std::vector path_; + H5NS::H5File file_; + H5NS::Group group_; + unsigned int dataSetThres_; + }; + + // Writer template implementation //////////////////////////////////////////// + template + void Hdf5Writer::writeSingleAttribute(const U &x, const std::string &name, + const H5NS::DataType &type) + { + H5NS::Attribute attribute; + hsize_t attrDim = 1; + H5NS::DataSpace attrSpace(1, &attrDim); + + attribute = group_.createAttribute(name, type, attrSpace); + attribute.write(type, &x); + } + + template + void Hdf5Writer::writeDefault(const std::string &s, const U &x) + { + writeSingleAttribute(x, s, Hdf5Type::type()); + } + + template <> + void Hdf5Writer::writeDefault(const std::string &s, const std::string &x); + + template + typename std::enable_if>::is_number, void>::type + Hdf5Writer::writeDefault(const std::string &s, const std::vector &x) + { + // alias to element type + typedef typename element>::type Element; + + // flatten the vector and getting dimensions + Flatten> flat(x); + std::vector dim; + const auto &flatx = flat.getFlatVector(); + + for (auto &d: flat.getDim()) + { + dim.push_back(d); + } + + // write to file + H5NS::DataSpace dataSpace(dim.size(), dim.data()); + + if (flatx.size() > dataSetThres_) + { + H5NS::DataSet dataSet; + + dataSet = group_.createDataSet(s, Hdf5Type::type(), dataSpace); + dataSet.write(flatx.data(), Hdf5Type::type()); + } + else + { + H5NS::Attribute attribute; + + attribute = group_.createAttribute(s, Hdf5Type::type(), dataSpace); + attribute.write(Hdf5Type::type(), flatx.data()); + } + } + + template + typename std::enable_if>::is_number, void>::type + Hdf5Writer::writeDefault(const std::string &s, const std::vector &x) + { + push(s); + writeSingleAttribute(x.size(), HDF5_GRID_GUARD "vector_size", + Hdf5Type::type()); + for (hsize_t i = 0; i < x.size(); ++i) + { + write(s + "_" + std::to_string(i), x[i]); + } + pop(); + } + + // Reader template implementation //////////////////////////////////////////// + template + void Hdf5Reader::readSingleAttribute(U &x, const std::string &name, + const H5NS::DataType &type) + { + H5NS::Attribute attribute; + + attribute = group_.openAttribute(name); + attribute.read(type, &x); + } + + template + void Hdf5Reader::readDefault(const std::string &s, U &output) + { + readSingleAttribute(output, s, Hdf5Type::type()); + } + + template <> + void Hdf5Reader::readDefault(const std::string &s, std::string &x); + + template + typename std::enable_if>::is_number, void>::type + Hdf5Reader::readDefault(const std::string &s, std::vector &x) + { + // alias to element type + typedef typename element>::type Element; + + // read the dimensions + H5NS::DataSpace dataSpace; + std::vector hdim; + std::vector dim; + hsize_t size = 1; + + if (group_.attrExists(s)) + { + dataSpace = group_.openAttribute(s).getSpace(); + } + else + { + dataSpace = group_.openDataSet(s).getSpace(); + } + hdim.resize(dataSpace.getSimpleExtentNdims()); + dataSpace.getSimpleExtentDims(hdim.data()); + for (auto &d: hdim) + { + dim.push_back(d); + size *= d; + } + + // read the flat vector + std::vector buf(size); + + if (size > dataSetThres_) + { + H5NS::DataSet dataSet; + + dataSet = group_.openDataSet(s); + dataSet.read(buf.data(), Hdf5Type::type()); + } + else + { + H5NS::Attribute attribute; + + attribute = group_.openAttribute(s); + attribute.read(Hdf5Type::type(), buf.data()); + } + + // reconstruct the multidimensional vector + Reconstruct> r(buf, dim); + + x = r.getVector(); + } + + template + typename std::enable_if>::is_number, void>::type + Hdf5Reader::readDefault(const std::string &s, std::vector &x) + { + uint64_t size; + + push(s); + readSingleAttribute(size, HDF5_GRID_GUARD "vector_size", + Hdf5Type::type()); + x.resize(size); + for (hsize_t i = 0; i < x.size(); ++i) + { + read(s + "_" + std::to_string(i), x[i]); + } + pop(); + } +} + +#endif diff --git a/lib/serialisation/Hdf5Type.h b/lib/serialisation/Hdf5Type.h new file mode 100644 index 00000000..8634f35b --- /dev/null +++ b/lib/serialisation/Hdf5Type.h @@ -0,0 +1,77 @@ +#ifndef GRID_SERIALISATION_HDF5_TYPE_H +#define GRID_SERIALISATION_HDF5_TYPE_H + +#include +#include +#include + +#ifndef H5_NO_NAMESPACE +#define H5NS H5 +#endif + +#define HDF5_NATIVE_TYPE(predType, cType)\ +template <>\ +class Hdf5Type\ +{\ +public:\ + static inline const H5NS::DataType & type(void)\ + {\ + return H5NS::PredType::predType;\ + }\ + static constexpr bool isNative = true;\ +}; + +#define DEFINE_HDF5_NATIVE_TYPES \ +HDF5_NATIVE_TYPE(NATIVE_B8, bool);\ +HDF5_NATIVE_TYPE(NATIVE_CHAR, char);\ +HDF5_NATIVE_TYPE(NATIVE_SCHAR, signed char);\ +HDF5_NATIVE_TYPE(NATIVE_UCHAR, unsigned char);\ +HDF5_NATIVE_TYPE(NATIVE_SHORT, short);\ +HDF5_NATIVE_TYPE(NATIVE_USHORT, unsigned short);\ +HDF5_NATIVE_TYPE(NATIVE_INT, int);\ +HDF5_NATIVE_TYPE(NATIVE_UINT, unsigned int);\ +HDF5_NATIVE_TYPE(NATIVE_LONG, long);\ +HDF5_NATIVE_TYPE(NATIVE_ULONG, unsigned long);\ +HDF5_NATIVE_TYPE(NATIVE_LLONG, long long);\ +HDF5_NATIVE_TYPE(NATIVE_ULLONG, unsigned long long);\ +HDF5_NATIVE_TYPE(NATIVE_FLOAT, float);\ +HDF5_NATIVE_TYPE(NATIVE_DOUBLE, double);\ +HDF5_NATIVE_TYPE(NATIVE_LDOUBLE, long double); + +namespace Grid +{ + template class Hdf5Type + { + public: + static constexpr bool isNative = false; + }; + + DEFINE_HDF5_NATIVE_TYPES; + + template + class Hdf5Type> + { + public: + static inline const H5NS::DataType & type(void) + { + if (typePtr_ == nullptr) + { + typePtr_.reset(new H5NS::CompType(sizeof(std::complex))); + typePtr_->insertMember("re", 0, Hdf5Type::type()); + typePtr_->insertMember("im", sizeof(R), Hdf5Type::type()); + } + + return *typePtr_; + } + static constexpr bool isNative = false; + private: + static std::unique_ptr typePtr_; + }; + + template + std::unique_ptr Hdf5Type>::typePtr_ = nullptr; +} + +#undef HDF5_NATIVE_TYPE + +#endif /* GRID_SERIALISATION_HDF5_TYPE_H */ diff --git a/lib/serialisation/MacroMagic.h b/lib/serialisation/MacroMagic.h index c78bba0c..8b027f30 100644 --- a/lib/serialisation/MacroMagic.h +++ b/lib/serialisation/MacroMagic.h @@ -109,40 +109,36 @@ THE SOFTWARE. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// #define GRID_MACRO_MEMBER(A,B) A B; +#define GRID_MACRO_COMP_MEMBER(A,B) result = (result and (lhs. B == rhs. B)); #define GRID_MACRO_OS_WRITE_MEMBER(A,B) os<< #A <<" "#B <<" = "<< obj. B <<" ; " <\ - static inline void write(Writer &WR,const std::string &s, const cname &obj){ \ - push(WR,s);\ - GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_WRITE_MEMBER,__VA_ARGS__)) \ - pop(WR);\ - } \ - \ - \ - template \ - static inline void read(Reader &RD,const std::string &s, cname &obj){ \ - push(RD,s);\ - GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_READ_MEMBER,__VA_ARGS__)) \ - pop(RD);\ - } \ - \ - \ - friend inline std::ostream & operator << (std::ostream &os, const cname &obj ) { \ - os<<"class "<<#cname<<" {"<\ +static inline void write(Writer &WR,const std::string &s, const cname &obj){ \ + push(WR,s);\ + GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_WRITE_MEMBER,__VA_ARGS__)) \ + pop(WR);\ +}\ +template \ +static inline void read(Reader &RD,const std::string &s, cname &obj){ \ + push(RD,s);\ + GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_READ_MEMBER,__VA_ARGS__)) \ + pop(RD);\ +}\ +friend inline std::ostream & operator << (std::ostream &os, const cname &obj ) { \ + os<<"class "<<#cname<<" {"<::type #define GRID_MACRO_ENUMVAL(A,B) A = B, @@ -150,44 +146,52 @@ THE SOFTWARE. #define GRID_MACRO_ENUMTEST(A,B) else if (buf == #A) {obj = GRID_ENUM_TYPE(obj)::A;} #define GRID_MACRO_ENUMCASEIO(A,B) case GRID_ENUM_TYPE(obj)::A: os << #A; break; -namespace Grid { - template - class EnumIO {}; -} - #define GRID_SERIALIZABLE_ENUM(name,undefname,...)\ - enum class name {\ - GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMVAL,__VA_ARGS__))\ - undefname = -1\ +class name: public Grid::Serializable\ +{\ +public:\ + enum EnumType\ + {\ + GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMVAL,__VA_ARGS__))\ + undefname = -1\ };\ +public:\ + name(void): value_(undefname) {};\ + name(EnumType value): value_(value) {};\ + template \ + static inline void write(Grid::Writer &WR,const std::string &s, const name &obj)\ + {\ + switch (obj.value_)\ + {\ + GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMCASE,__VA_ARGS__))\ + default: Grid::write(WR,s,#undefname); break;\ + }\ + }\ \ - template<>\ - class EnumIO {\ - public:\ - template \ - static inline void write(Writer &WR,const std::string &s, const name &obj){ \ - switch (obj) {\ - GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMCASE,__VA_ARGS__))\ - default: Grid::write(WR,s,#undefname); break;\ - }\ - }\ - \ - template \ - static inline void read(Reader &RD,const std::string &s, name &obj){ \ - std::string buf;\ - Grid::read(RD, s, buf);\ - if (buf == #undefname) {obj = name::undefname;}\ - GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMTEST,__VA_ARGS__))\ - else {obj = name::undefname;}\ - }\ - };\ - \ - inline std::ostream & operator << (std::ostream &os, const name &obj ) { \ + template \ + static inline void read(Grid::Reader &RD,const std::string &s, name &obj)\ + {\ + std::string buf;\ + Grid::read(RD, s, buf);\ + if (buf == #undefname) {obj = name::undefname;}\ + GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMTEST,__VA_ARGS__))\ + else {obj = name::undefname;}\ + }\ + inline operator EnumType(void) const\ + {\ + return value_;\ + }\ + inline friend std::ostream & operator<<(std::ostream &os, const name &obj)\ + {\ switch (obj) {\ - GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMCASEIO,__VA_ARGS__))\ - default: os << #undefname; break;\ + GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMCASEIO,__VA_ARGS__))\ + default: os << #undefname; break;\ }\ return os;\ - }; + }\ +private:\ + EnumType value_;\ +}; + #endif diff --git a/lib/serialisation/Serialisation.h b/lib/serialisation/Serialisation.h index cd86e252..0e1c7531 100644 --- a/lib/serialisation/Serialisation.h +++ b/lib/serialisation/Serialisation.h @@ -37,6 +37,11 @@ Author: Peter Boyle #include "TextIO.h" #include "XmlIO.h" #include "JSON_IO.h" + +#ifdef HAVE_HDF5 +#include "Hdf5IO.h" +#endif + ////////////////////////////////////////// // Todo: ////////////////////////////////////////// diff --git a/lib/simd/Grid_avx.h b/lib/simd/Grid_avx.h index ac6ec9f4..e9082e30 100644 --- a/lib/simd/Grid_avx.h +++ b/lib/simd/Grid_avx.h @@ -204,6 +204,29 @@ namespace Optimization { } }; + struct MultRealPart{ + inline __m256 operator()(__m256 a, __m256 b){ + __m256 ymm0; + ymm0 = _mm256_shuffle_ps(a,a,_MM_SELECT_FOUR_FOUR(2,2,0,0)); // ymm0 <- ar ar, + return _mm256_mul_ps(ymm0,b); // ymm0 <- ar bi, ar br + } + inline __m256d operator()(__m256d a, __m256d b){ + __m256d ymm0; + ymm0 = _mm256_shuffle_pd(a,a,0x0); // ymm0 <- ar ar, ar,ar b'00,00 + return _mm256_mul_pd(ymm0,b); // ymm0 <- ar bi, ar br + } + }; + struct MaddRealPart{ + inline __m256 operator()(__m256 a, __m256 b, __m256 c){ + __m256 ymm0 = _mm256_moveldup_ps(a); // ymm0 <- ar ar, + return _mm256_add_ps(_mm256_mul_ps( ymm0, b),c); + } + inline __m256d operator()(__m256d a, __m256d b, __m256d c){ + __m256d ymm0 = _mm256_shuffle_pd( a, a, 0x0 ); + return _mm256_add_pd(_mm256_mul_pd( ymm0, b),c); + } + }; + struct MultComplex{ // Complex float inline __m256 operator()(__m256 a, __m256 b){ @@ -618,7 +641,9 @@ namespace Optimization { typedef Optimization::Sub SubSIMD; typedef Optimization::Div DivSIMD; typedef Optimization::Mult MultSIMD; - typedef Optimization::MultComplex MultComplexSIMD; + typedef Optimization::MultComplex MultComplexSIMD; + typedef Optimization::MultRealPart MultRealPartSIMD; + typedef Optimization::MaddRealPart MaddRealPartSIMD; typedef Optimization::Conj ConjSIMD; typedef Optimization::TimesMinusI TimesMinusISIMD; typedef Optimization::TimesI TimesISIMD; diff --git a/lib/simd/Grid_avx512.h b/lib/simd/Grid_avx512.h index d6531d57..ebf99e16 100644 --- a/lib/simd/Grid_avx512.h +++ b/lib/simd/Grid_avx512.h @@ -189,6 +189,29 @@ namespace Optimization { // 2mul,4 mac +add+sub = 8 flop type insns // 3shuf + 2 (+shuf) = 5/6 simd perm and 1/2 the load. + struct MultRealPart{ + inline __m512 operator()(__m512 a, __m512 b){ + __m512 ymm0; + ymm0 = _mm512_moveldup_ps(a); // ymm0 <- ar ar, + return _mm512_mul_ps(ymm0,b); // ymm0 <- ar bi, ar br + } + inline __m512d operator()(__m512d a, __m512d b){ + __m512d ymm0; + ymm0 = _mm512_shuffle_pd(a,a,0x00); // ymm0 <- ar ar, ar,ar b'00,00 + return _mm512_mul_pd(ymm0,b); // ymm0 <- ar bi, ar br + } + }; + struct MaddRealPart{ + inline __m512 operator()(__m512 a, __m512 b, __m512 c){ + __m512 ymm0 = _mm512_moveldup_ps(a); // ymm0 <- ar ar, + return _mm512_fmadd_ps( ymm0, b, c); + } + inline __m512d operator()(__m512d a, __m512d b, __m512d c){ + __m512d ymm0 = _mm512_shuffle_pd( a, a, 0x00 ); + return _mm512_fmadd_pd( ymm0, b, c); + } + }; + struct MultComplex{ // Complex float inline __m512 operator()(__m512 a, __m512 b){ @@ -501,6 +524,8 @@ namespace Optimization { typedef Optimization::Mult MultSIMD; typedef Optimization::Div DivSIMD; typedef Optimization::MultComplex MultComplexSIMD; + typedef Optimization::MultRealPart MultRealPartSIMD; + typedef Optimization::MaddRealPart MaddRealPartSIMD; typedef Optimization::Conj ConjSIMD; typedef Optimization::TimesMinusI TimesMinusISIMD; typedef Optimization::TimesI TimesISIMD; diff --git a/lib/simd/Grid_generic.h b/lib/simd/Grid_generic.h index 62c78afb..91e9cda2 100644 --- a/lib/simd/Grid_generic.h +++ b/lib/simd/Grid_generic.h @@ -224,6 +224,21 @@ namespace Optimization { #define cmul(a, b, c, i)\ c[i] = a[i]*b[i] - a[i+1]*b[i+1];\ c[i+1] = a[i]*b[i+1] + a[i+1]*b[i]; + + struct MultRealPart{ + template + inline vec operator()(vec a, vec b){ + vec out; + + VECTOR_FOR(i, W::c, 1) + { + out.v[2*i] = a[2*i]*b[2*i]; + out.v[2*i+1] = a[2*i]*b[2*i+1]; + } + return out; + }; + }; + struct MultComplex{ // Complex @@ -456,6 +471,7 @@ namespace Optimization { typedef Optimization::Div DivSIMD; typedef Optimization::Mult MultSIMD; typedef Optimization::MultComplex MultComplexSIMD; + typedef Optimization::MultRealPart MultRealPartSIMD; typedef Optimization::Conj ConjSIMD; typedef Optimization::TimesMinusI TimesMinusISIMD; typedef Optimization::TimesI TimesISIMD; diff --git a/lib/simd/Grid_qpx.h b/lib/simd/Grid_qpx.h index bc86291d..99a9ea68 100644 --- a/lib/simd/Grid_qpx.h +++ b/lib/simd/Grid_qpx.h @@ -220,6 +220,14 @@ namespace Optimization { } }; + struct MultRealPart{ + // Complex double + inline vector4double operator()(vector4double a, vector4double b){ + // return vec_xmul(b, a); + return vec_xmul(a, b); + } + FLOAT_WRAP_2(operator(), inline) + }; struct MultComplex{ // Complex double inline vector4double operator()(vector4double a, vector4double b){ @@ -430,6 +438,7 @@ typedef Optimization::Sub SubSIMD; typedef Optimization::Mult MultSIMD; typedef Optimization::Div DivSIMD; typedef Optimization::MultComplex MultComplexSIMD; +typedef Optimization::MultRealPart MultRealPartSIMD; typedef Optimization::Conj ConjSIMD; typedef Optimization::TimesMinusI TimesMinusISIMD; typedef Optimization::TimesI TimesISIMD; diff --git a/lib/simd/Grid_sse4.h b/lib/simd/Grid_sse4.h index 560eda11..943756b2 100644 --- a/lib/simd/Grid_sse4.h +++ b/lib/simd/Grid_sse4.h @@ -177,6 +177,29 @@ namespace Optimization { } }; + struct MultRealPart{ + inline __m128 operator()(__m128 a, __m128 b){ + __m128 ymm0; + ymm0 = _mm_shuffle_ps(a,a,_MM_SELECT_FOUR_FOUR(2,2,0,0)); // ymm0 <- ar ar, + return _mm_mul_ps(ymm0,b); // ymm0 <- ar bi, ar br + } + inline __m128d operator()(__m128d a, __m128d b){ + __m128d ymm0; + ymm0 = _mm_shuffle_pd(a,a,0x0); // ymm0 <- ar ar, ar,ar b'00,00 + return _mm_mul_pd(ymm0,b); // ymm0 <- ar bi, ar br + } + }; + struct MaddRealPart{ + inline __m128 operator()(__m128 a, __m128 b, __m128 c){ + __m128 ymm0 = _mm_shuffle_ps(a,a,_MM_SELECT_FOUR_FOUR(2,2,0,0)); // ymm0 <- ar ar, + return _mm_add_ps(_mm_mul_ps( ymm0, b),c); + } + inline __m128d operator()(__m128d a, __m128d b, __m128d c){ + __m128d ymm0 = _mm_shuffle_pd( a, a, 0x0 ); + return _mm_add_pd(_mm_mul_pd( ymm0, b),c); + } + }; + struct MultComplex{ // Complex float inline __m128 operator()(__m128 a, __m128 b){ @@ -325,9 +348,11 @@ namespace Optimization { } } +#ifndef _mm_alignr_epi64 #define _mm_alignr_epi32(a,b,n) _mm_alignr_epi8(a,b,(n*4)%16) #define _mm_alignr_epi64(a,b,n) _mm_alignr_epi8(a,b,(n*8)%16) - +#endif + template static inline __m128 tRotate(__m128 in){ return (__m128)_mm_alignr_epi32((__m128i)in,(__m128i)in,n); }; template static inline __m128d tRotate(__m128d in){ return (__m128d)_mm_alignr_epi64((__m128i)in,(__m128i)in,n); }; @@ -415,6 +440,8 @@ namespace Optimization { typedef Optimization::Div DivSIMD; typedef Optimization::Mult MultSIMD; typedef Optimization::MultComplex MultComplexSIMD; + typedef Optimization::MultRealPart MultRealPartSIMD; + typedef Optimization::MaddRealPart MaddRealPartSIMD; typedef Optimization::Conj ConjSIMD; typedef Optimization::TimesMinusI TimesMinusISIMD; typedef Optimization::TimesI TimesISIMD; diff --git a/lib/simd/Grid_vector_types.h b/lib/simd/Grid_vector_types.h index fc217cc5..46e074db 100644 --- a/lib/simd/Grid_vector_types.h +++ b/lib/simd/Grid_vector_types.h @@ -101,6 +101,11 @@ template using IfNotInteger = Invoke +Out trinary(Input1 src_1, Input2 src_2, Input3 src_3, Operation op) { + return op(src_1, src_2, src_3); +} + template Out binary(Input1 src_1, Input2 src_2, Operation op) { return op(src_1, src_2); @@ -178,6 +183,7 @@ class Grid_simd { const Grid_simd *__restrict__ r) { *y = (*l) * (*r); } + friend inline void sub(Grid_simd *__restrict__ y, const Grid_simd *__restrict__ l, const Grid_simd *__restrict__ r) { @@ -188,7 +194,6 @@ class Grid_simd { const Grid_simd *__restrict__ r) { *y = (*l) + (*r); } - friend inline void mac(Grid_simd *__restrict__ y, const Scalar_type *__restrict__ a, const Grid_simd *__restrict__ x) { @@ -260,7 +265,7 @@ class Grid_simd { } //////////////////////////// - // opreator scalar * simd + // operator scalar * simd //////////////////////////// friend inline Grid_simd operator*(const Scalar_type &a, Grid_simd b) { Grid_simd va; @@ -446,6 +451,11 @@ inline void vbroadcast(Grid_simd &ret,const Grid_simd &src,int lane){ S* typepun =(S*) &src; vsplat(ret,typepun[lane]); } +template =0> +inline void rbroadcast(Grid_simd &ret,const Grid_simd &src,int lane){ + S* typepun =(S*) &src; + ret.v = unary(real(typepun[lane]), VsplatSIMD()); +} /////////////////////// // Splat @@ -462,6 +472,10 @@ template inline void vsplat(Grid_simd &ret, EnableIf, S> c) { vsplat(ret, real(c), imag(c)); } +template +inline void rsplat(Grid_simd &ret, EnableIf, S> c) { + vsplat(ret, real(c), real(c)); +} // if real fill with a, if complex fill with a in the real part (first function // above) @@ -563,6 +577,21 @@ inline Grid_simd operator-(Grid_simd a, Grid_simd b) { return ret; }; +// Distinguish between complex types and others +template = 0> +inline Grid_simd real_mult(Grid_simd a, Grid_simd b) { + Grid_simd ret; + ret.v = binary(a.v, b.v, MultRealPartSIMD()); + return ret; +}; +template = 0> +inline Grid_simd real_madd(Grid_simd a, Grid_simd b, Grid_simd c) { + Grid_simd ret; + ret.v = trinary(a.v, b.v, c.v, MaddRealPartSIMD()); + return ret; +}; + + // Distinguish between complex types and others template = 0> inline Grid_simd operator*(Grid_simd a, Grid_simd b) { diff --git a/lib/simd/Intel512avx.h b/lib/simd/Intel512avx.h index 19157db4..7b5964ad 100644 --- a/lib/simd/Intel512avx.h +++ b/lib/simd/Intel512avx.h @@ -95,10 +95,14 @@ Author: paboyle #define VIDUPd(SRC,DEST) "vpshufd $0xee," #SRC"," #DEST ";\n" // 32 bit level: 3,2,3,2 #define VIDUPf(SRC,DEST) "vmovshdup " #SRC ", " #DEST ";\n" -#define VBCASTRDUPd(OFF,A,DEST) "vbroadcastsd (" #OFF "*16+0)(" #A ")," #DEST ";\n" -#define VBCASTIDUPd(OFF,A,DEST) "vbroadcastsd (" #OFF "*16+8)(" #A ")," #DEST ";\n" -#define VBCASTRDUPf(OFF,PTR,DEST) "vbroadcastss (" #OFF "*8 +0)(" #PTR "), " #DEST ";\n" -#define VBCASTIDUPf(OFF,PTR,DEST) "vbroadcastss (" #OFF "*8 +4)(" #PTR "), " #DEST ";\n" +#define VBCASTRDUPd(OFF,A,DEST) "vbroadcastsd (" #OFF "*16+0)(" #A ")," #DEST ";\n" +#define VBCASTIDUPd(OFF,A,DEST) "vbroadcastsd (" #OFF "*16+8)(" #A ")," #DEST ";\n" +#define VBCASTRDUPf(OFF,PTR,DEST) "vbroadcastss (" #OFF "*8 +0)(" #PTR "), " #DEST ";\n" +#define VBCASTIDUPf(OFF,PTR,DEST) "vbroadcastss (" #OFF "*8 +4)(" #PTR "), " #DEST ";\n" +#define VBCASTCDUPf(OFF,A,DEST) "vbroadcastsd (" #OFF "*64 )(" #A ")," #DEST ";\n" +#define VBCASTZDUPf(OFF,A,DEST) "vbroadcastf32x4 (" #OFF "*64 )(" #A ")," #DEST ";\n" +#define VBCASTCDUP(OFF,A,DEST) VBCASTCDUPf(OFF,A,DEST) +#define VBCASTZDUP(OFF,A,DEST) VBCASTZDUPf(OFF,A,DEST) #define VMADDSUBf(A,B,accum) "vfmaddsub231ps " #A "," #B "," #accum ";\n" #define VMADDSUBd(A,B,accum) "vfmaddsub231pd " #A "," #B "," #accum ";\n" @@ -106,11 +110,15 @@ Author: paboyle #define VMADDSUBMEMd(O,P,B,accum) "vfmaddsub231pd " #O"*64("#P "),"#B "," #accum ";\n" +#define VMADDRDUPf(O,P,B,accum) "vfmadd231ps (" #O"*8+0)("#P "){1to16},"#B "," #accum ";\n" +#define VMADDIDUPf(O,P,B,accum) "vfmadd231ps (" #O"*8+4)("#P "){1to16},"#B "," #accum ";\n" #define VMADDSUBRDUPf(O,P,B,accum) "vfmaddsub231ps (" #O"*8+0)("#P "){1to16},"#B "," #accum ";\n" #define VMADDSUBIDUPf(O,P,B,accum) "vfmaddsub231ps (" #O"*8+4)("#P "){1to16},"#B "," #accum ";\n" #define VMULRDUPf(O,P,B,accum) "vmulps (" #O"*8+0)("#P "){1to16},"#B "," #accum ";\n" #define VMULIDUPf(O,P,B,accum) "vmulps (" #O"*8+4)("#P "){1to16},"#B "," #accum ";\n" +#define VMADDRDUPd(O,P,B,accum) "vfmadd231pd (" #O"*16+0)("#P "){1to8},"#B "," #accum ";\n" +#define VMADDIDUPd(O,P,B,accum) "vfmadd231pd (" #O"*16+8)("#P "){1to8},"#B "," #accum ";\n" #define VMADDSUBRDUPd(O,P,B,accum) "vfmaddsub231pd (" #O"*16+0)("#P "){1to8},"#B "," #accum ";\n" #define VMADDSUBIDUPd(O,P,B,accum) "vfmaddsub231pd (" #O"*16+8)("#P "){1to8},"#B "," #accum ";\n" #define VMULRDUPd(O,P,B,accum) "vmulpd (" #O"*16+0)("#P "){1to8},"#B "," #accum ";\n" diff --git a/lib/simd/Intel512common.h b/lib/simd/Intel512common.h index cfa20c26..e69e541c 100644 --- a/lib/simd/Intel512common.h +++ b/lib/simd/Intel512common.h @@ -87,7 +87,8 @@ Author: paboyle VACCTIMESMINUSI1d(A,ACC,tmp) \ VACCTIMESMINUSI2d(A,ACC,tmp) -#define LOAD64i(A,ptr) __asm__ ( "movq %0, %" #A : : "r"(ptr) : #A ); +#define LOAD64a(A,ptr) "movq %0, %" #A : : "r"(ptr) : #A +#define LOAD64i(A,ptr) __asm__ ( LOAD64a(A,ptr)); #define LOAD64(A,ptr) LOAD64i(A,ptr) #define VMOVf(A,DEST) "vmovaps " #A ", " #DEST ";\n" @@ -108,8 +109,8 @@ Author: paboyle //"vprefetche0 "#O"*64("#A");\n" "vprefetche1 ("#O"+12)*64("#A");\n" // "clevict0 "#O"*64("#A");\n" -#define VLOADf(OFF,PTR,DEST) "vmovaps " #OFF "*64(" #PTR "), " #DEST ";\n" -#define VLOADd(OFF,PTR,DEST) "vmovapd " #OFF "*64(" #PTR "), " #DEST ";\n" +#define VLOADf(OFF,PTR,DEST) "vmovups " #OFF "*64(" #PTR "), " #DEST ";\n" +#define VLOADd(OFF,PTR,DEST) "vmovupd " #OFF "*64(" #PTR "), " #DEST ";\n" #define VADDf(A,B,DEST) "vaddps " #A "," #B "," #DEST ";\n" #define VADDd(A,B,DEST) "vaddpd " #A "," #B "," #DEST ";\n" @@ -143,8 +144,8 @@ Author: paboyle #define VSTOREf(OFF,PTR,SRC) "vmovntps " #SRC "," #OFF "*64(" #PTR ")" ";\n" #define VSTOREd(OFF,PTR,SRC) "vmovntpd " #SRC "," #OFF "*64(" #PTR ")" ";\n" #else -#define VSTOREf(OFF,PTR,SRC) "vmovaps " #SRC "," #OFF "*64(" #PTR ")" ";\n" -#define VSTOREd(OFF,PTR,SRC) "vmovapd " #SRC "," #OFF "*64(" #PTR ")" ";\n" +#define VSTOREf(OFF,PTR,SRC) "vmovups " #SRC "," #OFF "*64(" #PTR ")" ";\n" +#define VSTOREd(OFF,PTR,SRC) "vmovupd " #SRC "," #OFF "*64(" #PTR ")" ";\n" #endif // Swaps Re/Im ; could unify this with IMCI diff --git a/lib/simd/Intel512double.h b/lib/simd/Intel512double.h index 224c593d..632b5639 100644 --- a/lib/simd/Intel512double.h +++ b/lib/simd/Intel512double.h @@ -144,10 +144,12 @@ Author: paboyle #define VMADDSUBMEM(O,P,B,accum) VMADDSUBMEMd(O,P,B,accum) #define VMADDMEM(O,P,B,accum) VMADDMEMd(O,P,B,accum) #define VMULMEM(O,P,B,accum) VMULMEMd(O,P,B,accum) +#undef VMADDRDUP #undef VMADDSUBRDUP #undef VMADDSUBIDUP #undef VMULRDUP #undef VMULIDUP +#define VMADDRDUP(O,P,B,accum) VMADDRDUPd(O,P,B,accum) #define VMADDSUBRDUP(O,P,B,accum) VMADDSUBRDUPd(O,P,B,accum) #define VMADDSUBIDUP(O,P,B,accum) VMADDSUBIDUPd(O,P,B,accum) #define VMULRDUP(O,P,B,accum) VMULRDUPd(O,P,B,accum) diff --git a/lib/simd/Intel512single.h b/lib/simd/Intel512single.h index 3fa47668..ed135651 100644 --- a/lib/simd/Intel512single.h +++ b/lib/simd/Intel512single.h @@ -144,10 +144,12 @@ Author: paboyle #define VMADDMEM(O,P,B,accum) VMADDMEMf(O,P,B,accum) #define VMULMEM(O,P,B,accum) VMULMEMf(O,P,B,accum) +#undef VMADDRDUP #undef VMADDSUBRDUP #undef VMADDSUBIDUP #undef VMULRDUP #undef VMULIDUP +#define VMADDRDUP(O,P,B,accum) VMADDRDUPf(O,P,B,accum) #define VMADDSUBRDUP(O,P,B,accum) VMADDSUBRDUPf(O,P,B,accum) #define VMADDSUBIDUP(O,P,B,accum) VMADDSUBIDUPf(O,P,B,accum) #define VMULRDUP(O,P,B,accum) VMULRDUPf(O,P,B,accum) diff --git a/lib/stencil/Lebesgue.cc b/lib/stencil/Lebesgue.cc index c30c14c7..c83975a9 100644 --- a/lib/stencil/Lebesgue.cc +++ b/lib/stencil/Lebesgue.cc @@ -65,7 +65,7 @@ void LebesgueOrder::CartesianBlocking(void) { _LebesgueReorder.resize(0); - std::cout << GridLogDebug << " CartesianBlocking "; + // std::cout << GridLogDebug << " CartesianBlocking "; // for(int d=0;d Make.inc echo >> Make.inc echo CCFILES=$CCFILES >> Make.inc @@ -24,10 +23,11 @@ for subdir in $dirs; do echo "tests: ${TESTLIST} ${SUB}" > Make.inc echo ${PREF}_PROGRAMS = ${TESTLIST} >> Make.inc echo >> Make.inc + HADLINK=`[ $subdir = './hadrons' ] && echo '-lHadrons '` for f in $TESTS; do BNAME=`basename $f .cc` - echo ${BNAME}_SOURCES=$f >> Make.inc - echo ${BNAME}_LDADD=-lGrid>> Make.inc + echo ${BNAME}_SOURCES=$f >> Make.inc + echo ${BNAME}_LDADD=${HADLINK}-lGrid >> Make.inc echo >> Make.inc done if [ $subdir != '.' ]; then diff --git a/tests/IO/Test_serialisation.cc b/tests/IO/Test_serialisation.cc index 3bd83848..212b653f 100644 --- a/tests/IO/Test_serialisation.cc +++ b/tests/IO/Test_serialisation.cc @@ -29,132 +29,155 @@ Author: Peter Boyle /* END LEGAL */ #include -namespace Grid { - - GRID_SERIALIZABLE_ENUM(myenum, undef, red, 1, blue, 2, green, 3); - - class myclass: Serializable { - public: - - GRID_SERIALIZABLE_CLASS_MEMBERS(myclass, - myenum, e, - std::vector, ve, - std::string, name, - int, x, - double, y, - bool , b, - std::vector, array, - std::vector>, twodimarray, - ); - - myclass() {} - myclass(int i) - : array(4,5.1), twodimarray(3,std::vector(2,1.23456)), ve({myenum::blue, myenum::red}) - { - e=myenum::red; - x=i; - y=2*i; - b=true; - name="bother said pooh"; - } - }; - - -} using namespace Grid; -int16_t i16 = 1; +GRID_SERIALIZABLE_ENUM(myenum, undef, red, 1, blue, 2, green, 3); + +class myclass: Serializable { +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(myclass, + myenum, e, + std::vector, ve, + std::string, name, + int, x, + double, y, + bool , b, + std::vector, array, + std::vector>, twodimarray, + std::vector>>, cmplx3darray + ); + myclass() {} + myclass(int i) + : array(4,5.1) + , twodimarray(3,std::vector(5, 1.23456)) + , cmplx3darray(3,std::vector>(5, std::vector(7, Complex(1.2, 3.4)))) + , ve(2, myenum::blue) + { + e=myenum::red; + x=i; + y=2*i; + b=true; + name="bother said pooh"; + } +}; + +int16_t i16 = 1; uint16_t u16 = 2; -int32_t i32 = 3; +int32_t i32 = 3; uint32_t u32 = 4; -int64_t i64 = 5; +int64_t i64 = 5; uint64_t u64 = 6; -float f = M_PI; -double d = 2*M_PI; -bool b = false; +float f = M_PI; +double d = 2*M_PI; +bool b = false; + +template +void ioTest(const std::string &filename, const O &object, const std::string &name) +{ + // writer needs to be destroyed so that writing physically happens + { + W writer(filename); + + write(writer, "testobject", object); + } + + R reader(filename); + O buf; + bool good; + + read(reader, "testobject", buf); + good = (object == buf); + std::cout << name << " IO test: " << (good ? "success" : "failure"); + std::cout << std::endl; + if (!good) exit(EXIT_FAILURE); +} int main(int argc,char **argv) { - { - XmlWriter WR("bother.xml"); - - // test basic type writing - push(WR,"BasicTypes"); - write(WR,std::string("i16"),i16); - write(WR,"u16",u16); - write(WR,"i32",i32); - write(WR,"u32",u32); - write(WR,"i64",i64); - write(WR,"u64",u64); - write(WR,"f",f); - write(WR,"d",d); - write(WR,"b",b); - pop(WR); - - // test serializable class writing - myclass obj(1234); // non-trivial constructor - write(WR,"obj",obj); - WR.write("obj2", obj); - std::cout << obj << std::endl; - - std::vector vec; - vec.push_back(myclass(1234)); - vec.push_back(myclass(5678)); - vec.push_back(myclass(3838)); - write(WR, "objvec", vec); - }; + std::cout << "==== basic IO" << std::endl; + XmlWriter WR("bother.xml"); + + // test basic type writing + std::cout << "-- basic writing to 'bother.xml'..." << std::endl; + push(WR,"BasicTypes"); + write(WR,std::string("i16"),i16); + write(WR,"u16",u16); + write(WR,"i32",i32); + write(WR,"u32",u32); + write(WR,"i64",i64); + write(WR,"u64",u64); + write(WR,"f",f); + write(WR,"d",d); + write(WR,"b",b); + pop(WR); + + // test serializable class writing + myclass obj(1234); // non-trivial constructor + std::vector vec; + + std::cout << "-- serialisable class writing to 'bother.xml'..." << std::endl; + write(WR,"obj",obj); + WR.write("obj2", obj); + vec.push_back(myclass(1234)); + vec.push_back(myclass(5678)); + vec.push_back(myclass(3838)); + write(WR, "objvec", vec); + std::cout << "-- serialisable class writing to std::cout:" << std::endl; + std::cout << obj << std::endl; + std::cout << "-- serialisable class comparison:" << std::endl; + std::cout << "vec[0] == obj: " << ((vec[0] == obj) ? "true" : "false") << std::endl; + std::cout << "vec[1] == obj: " << ((vec[1] == obj) ? "true" : "false") << std::endl; // read tests - myclass copy1, copy2, copy3; - std::vector veccopy1, veccopy2, veccopy3; + std::cout << "\n==== IO self-consistency tests" << std::endl; //// XML - { - XmlReader RD("bother.xml"); - read(RD,"obj",copy1); - read(RD,"objvec", veccopy1); - std::cout << "Loaded (XML) -----------------" << std::endl; - std::cout << copy1 << std::endl << veccopy1 << std::endl; - } + ioTest("iotest.xml", obj, "XML (object) "); + ioTest("iotest.xml", vec, "XML (vector of objects)"); //// binary - { - BinaryWriter BWR("bother.bin"); - write(BWR,"discard",copy1 ); - write(BWR,"discard",veccopy1 ); - } - { - BinaryReader BRD("bother.bin"); - read (BRD,"discard",copy2 ); - read (BRD,"discard",veccopy2 ); - std::cout << "Loaded (bin) -----------------" << std::endl; - std::cout << copy2 << std::endl << veccopy2 << std::endl; - } + ioTest("iotest.bin", obj, "binary (object) "); + ioTest("iotest.bin", vec, "binary (vector of objects)"); //// text - { - TextWriter TWR("bother.txt"); - write(TWR,"discard",copy1 ); - write(TWR,"discard",veccopy1 ); - } - { - TextReader TRD("bother.txt"); - read (TRD,"discard",copy3 ); - read (TRD,"discard",veccopy3 ); - std::cout << "Loaded (txt) -----------------" << std::endl; - std::cout << copy3 << std::endl << veccopy3 << std::endl; - } + ioTest("iotest.dat", obj, "text (object) "); + ioTest("iotest.dat", vec, "text (vector of objects)"); + //// HDF5 +#ifdef HAVE_HDF5 + ioTest("iotest.h5", obj, "HDF5 (object) "); + ioTest("iotest.h5", vec, "HDF5 (vector of objects)"); +#endif - std::vector iv = strToVec("1 2 2 4"); - std::vector sv = strToVec("bli bla blu"); + std::cout << "\n==== vector flattening/reconstruction" << std::endl; + typedef std::vector>> vec3d; - for (auto &e: iv) + vec3d dv, buf; + double d = 0.; + + dv.resize(4); + for (auto &v1: dv) { - std::cout << e << " "; - } - std::cout << std::endl; - for (auto &e: sv) - { - std::cout << e << " "; + v1.resize(3); + for (auto &v2: v1) + { + v2.resize(5); + for (auto &x: v2) + { + x = d++; + } + } } + std::cout << "original 3D vector:" << std::endl; + std::cout << dv << std::endl; + + Flatten flatdv(dv); + + std::cout << "\ndimensions:" << std::endl; + std::cout << flatdv.getDim() << std::endl; + std::cout << "\nflattened vector:" << std::endl; + std::cout << flatdv.getFlatVector() << std::endl; + + Reconstruct rec(flatdv.getFlatVector(), flatdv.getDim()); + std::cout << "\nreconstructed vector:" << std::endl; + std::cout << flatdv.getVector() << std::endl; std::cout << std::endl; @@ -211,4 +234,6 @@ int main(int argc,char **argv) read(RD,"name", name); } + + } diff --git a/tests/debug/Test_cayley_even_odd_vec.cc b/tests/Test_cayley_even_odd_vec.cc similarity index 100% rename from tests/debug/Test_cayley_even_odd_vec.cc rename to tests/Test_cayley_even_odd_vec.cc diff --git a/tests/hadrons/Makefile.am b/tests/hadrons/Makefile.am new file mode 100644 index 00000000..c8ec1612 --- /dev/null +++ b/tests/hadrons/Makefile.am @@ -0,0 +1,3 @@ +AM_LDFLAGS += -L../../extras/Hadrons + +include Make.inc diff --git a/tests/hadrons/Test_hadrons_meson_3pt.cc b/tests/hadrons/Test_hadrons_meson_3pt.cc new file mode 100644 index 00000000..56d0efa7 --- /dev/null +++ b/tests/hadrons/Test_hadrons_meson_3pt.cc @@ -0,0 +1,170 @@ +/******************************************************************************* + Grid physics library, www.github.com/paboyle/Grid + + Source file: tests/hadrons/Test_hadrons_meson_3pt.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 Hadrons; + +int main(int argc, char *argv[]) +{ + // 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; + std::vector flavour = {"l", "s", "c1", "c2", "c3"}; + std::vector mass = {.01, .04, .2 , .25 , .3 }; + unsigned int nt = GridDefaultLatt()[Tp]; + + // global parameters + Application::GlobalPar globalPar; + globalPar.trajCounter.start = 1500; + globalPar.trajCounter.end = 1520; + globalPar.trajCounter.step = 20; + globalPar.seed = "1 2 3 4"; + globalPar.genetic.maxGen = 1000; + globalPar.genetic.maxCstGen = 200; + globalPar.genetic.popSize = 20; + globalPar.genetic.mutationRate = .1; + application.setPar(globalPar); + + // gauge field + application.createModule("gauge"); + for (unsigned int i = 0; i < flavour.size(); ++i) + { + // actions + MAction::DWF::Par actionPar; + actionPar.gauge = "gauge"; + actionPar.Ls = 12; + actionPar.M5 = 1.8; + actionPar.mass = mass[i]; + application.createModule("DWF_" + flavour[i], actionPar); + + // solvers + MSolver::RBPrecCG::Par solverPar; + solverPar.action = "DWF_" + flavour[i]; + solverPar.residual = 1.0e-8; + application.createModule("CG_" + flavour[i], + solverPar); + } + for (unsigned int t = 0; t < nt; t += 1) + { + std::string srcName; + std::vector qName; + std::vector> seqName; + + // Z2 source + MSource::Z2::Par z2Par; + z2Par.tA = t; + z2Par.tB = t; + srcName = "z2_" + std::to_string(t); + application.createModule(srcName, z2Par); + for (unsigned int i = 0; i < flavour.size(); ++i) + { + // sequential sources + MSource::SeqGamma::Par seqPar; + qName.push_back("QZ2_" + flavour[i] + "_" + std::to_string(t)); + seqPar.q = qName[i]; + seqPar.tA = (t + nt/4) % nt; + seqPar.tB = (t + nt/4) % nt; + seqPar.mom = "1. 0. 0. 0."; + seqName.push_back(std::vector(Nd)); + for (unsigned int mu = 0; mu < Nd; ++mu) + { + seqPar.gamma = 0x1 << mu; + seqName[i][mu] = "G" + std::to_string(seqPar.gamma) + + "_" + std::to_string(seqPar.tA) + "-" + + qName[i]; + application.createModule(seqName[i][mu], seqPar); + } + + // propagators + Quark::Par quarkPar; + quarkPar.solver = "CG_" + flavour[i]; + quarkPar.source = srcName; + application.createModule(qName[i], quarkPar); + for (unsigned int mu = 0; mu < Nd; ++mu) + { + quarkPar.source = seqName[i][mu]; + seqName[i][mu] = "Q_" + flavour[i] + "-" + seqName[i][mu]; + application.createModule(seqName[i][mu], quarkPar); + } + } + + // contractions + MContraction::Meson::Par mesPar; + for (unsigned int i = 0; i < flavour.size(); ++i) + for (unsigned int j = i; j < flavour.size(); ++j) + { + mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j]; + mesPar.q1 = qName[i]; + mesPar.q2 = qName[j]; + application.createModule("meson_Z2_" + + std::to_string(t) + + "_" + + flavour[i] + + flavour[j], + mesPar); + } + for (unsigned int i = 0; i < flavour.size(); ++i) + for (unsigned int j = 0; j < flavour.size(); ++j) + for (unsigned int mu = 0; mu < Nd; ++mu) + { + MContraction::Meson::Par mesPar; + + mesPar.output = "3pt/Z2_" + flavour[i] + flavour[j] + "_" + + std::to_string(mu); + mesPar.q1 = qName[i]; + mesPar.q2 = seqName[j][mu]; + application.createModule("3pt_Z2_" + + std::to_string(t) + + "_" + + flavour[i] + + flavour[j] + + "_" + + std::to_string(mu), + mesPar); + } + } + + // execution + application.saveParameterFile("meson3pt.xml"); + application.run(); + + // epilogue + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +} diff --git a/tests/hadrons/Test_hadrons_spectrum.cc b/tests/hadrons/Test_hadrons_spectrum.cc new file mode 100644 index 00000000..ea6e4479 --- /dev/null +++ b/tests/hadrons/Test_hadrons_spectrum.cc @@ -0,0 +1,132 @@ +/******************************************************************************* + Grid physics library, www.github.com/paboyle/Grid + + Source file: tests/hadrons/Test_hadrons_spectrum.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 Hadrons; + +int main(int argc, char *argv[]) +{ + // 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; + std::vector flavour = {"l", "s", "c1", "c2", "c3"}; + std::vector mass = {.01, .04, .2 , .25 , .3 }; + + // global parameters + Application::GlobalPar globalPar; + globalPar.trajCounter.start = 1500; + globalPar.trajCounter.end = 1520; + globalPar.trajCounter.step = 20; + globalPar.seed = "1 2 3 4"; + application.setPar(globalPar); + // gauge field + application.createModule("gauge"); + // sources + MSource::Z2::Par z2Par; + z2Par.tA = 0; + z2Par.tB = 0; + application.createModule("z2", z2Par); + MSource::Point::Par ptPar; + ptPar.position = "0 0 0 0"; + application.createModule("pt", ptPar); + for (unsigned int i = 0; i < flavour.size(); ++i) + { + // actions + MAction::DWF::Par actionPar; + actionPar.gauge = "gauge"; + actionPar.Ls = 12; + actionPar.M5 = 1.8; + actionPar.mass = mass[i]; + application.createModule("DWF_" + flavour[i], actionPar); + + // solvers + MSolver::RBPrecCG::Par solverPar; + solverPar.action = "DWF_" + flavour[i]; + solverPar.residual = 1.0e-8; + application.createModule("CG_" + flavour[i], + solverPar); + + // propagators + Quark::Par quarkPar; + quarkPar.solver = "CG_" + flavour[i]; + quarkPar.source = "pt"; + application.createModule("Qpt_" + flavour[i], quarkPar); + quarkPar.source = "z2"; + application.createModule("QZ2_" + flavour[i], quarkPar); + } + for (unsigned int i = 0; i < flavour.size(); ++i) + for (unsigned int j = i; j < flavour.size(); ++j) + { + MContraction::Meson::Par mesPar; + + mesPar.output = "mesons/pt_" + flavour[i] + flavour[j]; + mesPar.q1 = "Qpt_" + flavour[i]; + mesPar.q2 = "Qpt_" + flavour[j]; + application.createModule("meson_pt_" + + flavour[i] + flavour[j], + mesPar); + mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j]; + mesPar.q1 = "QZ2_" + flavour[i]; + mesPar.q2 = "QZ2_" + flavour[j]; + application.createModule("meson_Z2_" + + flavour[i] + flavour[j], + mesPar); + } + for (unsigned int i = 0; i < flavour.size(); ++i) + for (unsigned int j = i; j < flavour.size(); ++j) + for (unsigned int k = j; k < flavour.size(); ++k) + { + MContraction::Baryon::Par barPar; + + barPar.output = "baryons/pt_" + flavour[i] + flavour[j] + flavour[k]; + barPar.q1 = "Qpt_" + flavour[i]; + barPar.q2 = "Qpt_" + flavour[j]; + barPar.q3 = "Qpt_" + flavour[k]; + application.createModule( + "baryon_pt_" + flavour[i] + flavour[j] + flavour[k], barPar); + } + + // execution + application.saveParameterFile("spectrum.xml"); + application.run(); + + // epilogue + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +}