From f709329d969cc832fce71e8861d3f345e3f2a2c3 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Wed, 17 Oct 2018 20:26:48 +0100 Subject: [PATCH] Hadrons: first version of a contractor utility --- Hadrons/Global.hpp | 1 + Hadrons/Utilities/Contractor.cc | 215 +++++++++++++++++++++++++++++++ Hadrons/Utilities/Contractor.hpp | 12 ++ Hadrons/Utilities/Makefile.am | 5 +- 4 files changed, 232 insertions(+), 1 deletion(-) create mode 100644 Hadrons/Utilities/Contractor.cc create mode 100644 Hadrons/Utilities/Contractor.hpp diff --git a/Hadrons/Global.hpp b/Hadrons/Global.hpp index e4003e73..d8047f09 100644 --- a/Hadrons/Global.hpp +++ b/Hadrons/Global.hpp @@ -32,6 +32,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include +#include #include #include diff --git a/Hadrons/Utilities/Contractor.cc b/Hadrons/Utilities/Contractor.cc new file mode 100644 index 00000000..ee7b6f2b --- /dev/null +++ b/Hadrons/Utilities/Contractor.cc @@ -0,0 +1,215 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Utilities/Contractor.cc + +Copyright (C) 2015-2018 + + +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 */ +#define DV_DEBUG + +#include +#include +#include + +using namespace Grid; +using namespace QCD; +using namespace Hadrons; + +namespace Contractor +{ + class GlobalPar: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(GlobalPar, + unsigned int, nt, + std::string, diskVectorDir, + std::string, output); + }; + + class A2AMatrixPar: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(A2AMatrixPar, + std::string, file, + std::string, dataset, + unsigned int, cacheSize, + std::string, name); + }; + + class ProductPar: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(ProductPar, + std::string, terms, + std::vector, timeRange); + }; +} + +struct ContractorPar +{ + Contractor::GlobalPar global; + std::vector a2aMatrix; + std::vector product; +}; + +std::set parseTimeRange(const std::string str) +{ + std::regex rex("([0-9]+)|(([0-9]+)\\.\\.([0-9]+))"); + std::smatch sm; + std::vector rstr = strToVec(str); + std::set t; + + for (auto &s: rstr) + { + std::regex_match(s, sm, rex); + if (sm[1].matched) + { + t.insert(std::stoi(sm[1].str())); + } + else if (sm[2].matched) + { + unsigned int ta, tb; + + ta = std::stoi(sm[3].str()); + tb = std::stoi(sm[4].str()); + for (unsigned int ti = ta; ti <= tb; ++ti) + { + t.insert(ti); + } + } + } + + return t; +} + +int main(int argc, char* argv[]) +{ + // parse command line + std::string parFilename; + + if (argc != 2) + { + std::cerr << "usage: " << argv[0] << " "; + std::cerr << std::endl; + + return EXIT_FAILURE; + } + parFilename = argv[1]; + + // parse parameter file + ContractorPar par; + + unsigned int nMat, nCont; + { + XmlReader reader(parFilename); + read(reader, "global", par.global); + read(reader, "a2aMatrix", par.a2aMatrix); + read(reader, "product", par.product); + } + nMat = par.a2aMatrix.size(); + nCont = par.product.size(); + + // create diskvectors + std::map> a2aMat; + unsigned int cacheSize; + + for (auto &p: par.a2aMatrix) + { + std::string dirName = par.global.diskVectorDir + "/" + p.name; + + a2aMat.emplace(p.name, EigenDiskVector(dirName, par.global.nt, p.cacheSize)); + } + + // load data + for (unsigned int i = 0; i < a2aMat.size(); ++i) + { + auto &p = par.a2aMatrix[i]; + double t, size; + + std::cout << "-- Loading '" << p.file << "'..." << std::endl; + + A2AMatrixIo a2aIo(p.file, p.dataset, par.global.nt); + + a2aIo.load(a2aMat.at(p.name), &t); + std::cout << "Read " << a2aIo.getSize() << " bytes in " << t << " usec, " << a2aIo.getSize()/t*1.0e6/1024/1024 << " MB/s" << std::endl; + } + + // contract + EigenDiskVector::Matrix buf; + + for (auto &p: par.product) + { + std::vector term = strToVec(p.terms); + std::vector> times; + + if (term.size() != p.timeRange.size()) + { + HADRONS_ERROR(Size, "number of terms (" + std::to_string(term.size()) + + ") different from number of time ranges (" + + std::to_string(p.timeRange.size()) + ")"); + } + for (auto &s: p.timeRange) + { + times.push_back(parseTimeRange(s)); + } + + // test: just 2-pt function for now + if (term.size() != 2) + { + HADRONS_ERROR(Implementation, "only 2-pt function implemented"); + } + + std::vector corr(par.global.nt); + double tTrace = 0., flops, bytes; + +#define TIME_MOD(t) (((t) + par.global.nt) % par.global.nt) + for (unsigned int t = 0; t < par.global.nt; ++t) + { + corr[t] = ComplexD(0., 0.); + } + for (auto &trans: times[1]) + for (auto &t0: times[0]) + { + unsigned int tt0 = TIME_MOD(t0 + trans); + EigenDiskVector::Matrix m0 = a2aMat.at(term[0])[tt0]; + EigenDiskVector::Matrix m1 = a2aMat.at(term[1])[trans]; + + std::cout << "tr(" << term[0] << "[" << tt0 << "]*" + << term[1] << "[" << trans << "])"; + tTrace = -usecond(); + corr[t0] += (m0*m1).trace(); + tTrace += usecond(); + flops = 6.*m0.rows()*m0.cols() + 2.*(m0.rows()*m0.cols() - 1.); + bytes = 2.*m0.rows()*m0.cols()*sizeof(ComplexD); + std::cout << " -- Perf " << flops/tTrace/1.0e3 << " GFlop/s " + << bytes/tTrace*1.0e6/1024/1024/1024 << " GB/s "<< std::endl; + } + for (unsigned int t = 0; t < par.global.nt; ++t) + { + std::cout << t << " " << corr[t] << std::endl; + } +#undef TIME_MOD + } + + + return EXIT_SUCCESS; +} diff --git a/Hadrons/Utilities/Contractor.hpp b/Hadrons/Utilities/Contractor.hpp new file mode 100644 index 00000000..9640c7c8 --- /dev/null +++ b/Hadrons/Utilities/Contractor.hpp @@ -0,0 +1,12 @@ +#ifndef Hadrons_Contractor_hpp_ +#define Hadrons_Contractor_hpp_ + +#include + +BEGIN_HADRONS_NAMESPACE + + + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Contractor_hpp_ diff --git a/Hadrons/Utilities/Makefile.am b/Hadrons/Utilities/Makefile.am index 529def24..1b257c32 100644 --- a/Hadrons/Utilities/Makefile.am +++ b/Hadrons/Utilities/Makefile.am @@ -1,4 +1,4 @@ -bin_PROGRAMS = HadronsXmlRun HadronsFermionEP64To32 +bin_PROGRAMS = HadronsXmlRun HadronsFermionEP64To32 HadronsContractor HadronsXmlRun_SOURCES = HadronsXmlRun.cc HadronsXmlRun_LDADD = ../libHadrons.a ../../Grid/libGrid.a @@ -6,3 +6,6 @@ HadronsXmlRun_LDADD = ../libHadrons.a ../../Grid/libGrid.a HadronsFermionEP64To32_SOURCES = EigenPackCast.cc HadronsFermionEP64To32_CXXFLAGS = $(AM_CXXFLAGS) -DFIN=WilsonImplD::FermionField -DFOUT=WilsonImplF::FermionField HadronsFermionEP64To32_LDADD = ../libHadrons.a ../../Grid/libGrid.a + +HadronsContractor_SOURCES = Contractor.cc +HadronsContractor_LDADD = ../libHadrons.a ../../Grid/libGrid.a