From fe6a372f752a4f0040a93b6dfad918ca5b078ca2 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Mon, 8 Oct 2018 15:14:08 +0100 Subject: [PATCH 01/37] Hadrons: fixes and cleaning in the scalar SU(N) part --- Hadrons/Modules.hpp | 1 - Hadrons/Modules/MScalarSUN/TimeMomProbe.cc | 38 --- Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp | 268 -------------------- Hadrons/Modules/MScalarSUN/TrMag.hpp | 3 +- Hadrons/Modules/MScalarSUN/Utils.hpp | 2 +- Hadrons/modules.inc | 2 - 6 files changed, 3 insertions(+), 311 deletions(-) delete mode 100644 Hadrons/Modules/MScalarSUN/TimeMomProbe.cc delete mode 100644 Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp diff --git a/Hadrons/Modules.hpp b/Hadrons/Modules.hpp index 2859f0fc..6d56ea2d 100644 --- a/Hadrons/Modules.hpp +++ b/Hadrons/Modules.hpp @@ -50,7 +50,6 @@ #include #include #include -#include #include #include #include diff --git a/Hadrons/Modules/MScalarSUN/TimeMomProbe.cc b/Hadrons/Modules/MScalarSUN/TimeMomProbe.cc deleted file mode 100644 index 5bfa203a..00000000 --- a/Hadrons/Modules/MScalarSUN/TimeMomProbe.cc +++ /dev/null @@ -1,38 +0,0 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: Hadrons/Modules/MScalarSUN/TimeMomProbe.cc - -Copyright (C) 2015-2018 - -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 MScalarSUN; - -template class Grid::Hadrons::MScalarSUN::TTimeMomProbe>; -template class Grid::Hadrons::MScalarSUN::TTimeMomProbe>; -template class Grid::Hadrons::MScalarSUN::TTimeMomProbe>; -template class Grid::Hadrons::MScalarSUN::TTimeMomProbe>; -template class Grid::Hadrons::MScalarSUN::TTimeMomProbe>; diff --git a/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp b/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp deleted file mode 100644 index a73da8de..00000000 --- a/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp +++ /dev/null @@ -1,268 +0,0 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp - -Copyright (C) 2015-2018 - -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_MScalarSUN_TimeMomProbe_hpp_ -#define Hadrons_MScalarSUN_TimeMomProbe_hpp_ - -#include -#include -#include -#include - -BEGIN_HADRONS_NAMESPACE - -/****************************************************************************** - * n-point functions O(t,p)*tr(phi(t_1,p_1)*...*phi(t_n,p_n)) * - ******************************************************************************/ -BEGIN_MODULE_NAMESPACE(MScalarSUN) - -class TimeMomProbePar: Serializable -{ -public: - GRID_SERIALIZABLE_CLASS_MEMBERS(TimeMomProbePar, - std::string, field, - std::vector, op, - std::vector>, timeMom, - std::string, output); -}; - -class TimeMomProbeResult: Serializable -{ -public: - GRID_SERIALIZABLE_CLASS_MEMBERS(TimeMomProbeResult, - std::string, op, - std::vector>, timeMom, - std::vector, data); -}; - -template -class TTimeMomProbe: public Module -{ -public: - typedef typename SImpl::Field Field; - typedef typename SImpl::SiteField::scalar_object Site; - typedef typename SImpl::ComplexField ComplexField; - typedef std::vector SlicedOp; -public: - // constructor - TTimeMomProbe(const std::string name); - // destructor - virtual ~TTimeMomProbe(void) {}; - // dependency relation - virtual std::vector getInput(void); - virtual std::vector getOutput(void); - // setup - virtual void setup(void); - // execution - virtual void execute(void); -private: - void vectorModulo(std::vector &v); -}; - -MODULE_REGISTER_TMP(TimeMomProbeSU2, TTimeMomProbe>, MScalarSUN); -MODULE_REGISTER_TMP(TimeMomProbeSU3, TTimeMomProbe>, MScalarSUN); -MODULE_REGISTER_TMP(TimeMomProbeSU4, TTimeMomProbe>, MScalarSUN); -MODULE_REGISTER_TMP(TimeMomProbeSU5, TTimeMomProbe>, MScalarSUN); -MODULE_REGISTER_TMP(TimeMomProbeSU6, TTimeMomProbe>, MScalarSUN); - -/****************************************************************************** - * TTimeMomProbe implementation * - ******************************************************************************/ -// constructor ///////////////////////////////////////////////////////////////// -template -TTimeMomProbe::TTimeMomProbe(const std::string name) -: Module(name) -{} - -// dependencies/products /////////////////////////////////////////////////////// -template -std::vector TTimeMomProbe::getInput(void) -{ - std::vector in = par().op; - - in.push_back(par().field); - - return in; -} - -template -std::vector TTimeMomProbe::getOutput(void) -{ - std::vector out; - - return out; -} - -// setup /////////////////////////////////////////////////////////////////////// -template -void TTimeMomProbe::setup(void) -{ - envTmpLat(ComplexField, "ftBuf"); - envTmpLat(Field, "ftMatBuf"); -} - -// execution /////////////////////////////////////////////////////////////////// -// NB: time is direction 0 -template -void TTimeMomProbe::vectorModulo(std::vector &v) -{ - for (unsigned int mu = 0; mu < env().getNd(); ++mu) - { - auto d = env().getDim(mu); - v[mu] = ((v[mu] % d) + d) % d; - } -} - -template -void TTimeMomProbe::execute(void) -{ - const unsigned int nd = env().getNd(); - const unsigned int nt = env().getDim(0); - double partVol = 1.; - std::set> timeMomSet; - std::vector>> timeMom; - std::vector> transferMom; - FFT fft(envGetGrid(Field)); - std::vector dMask(nd, 1); - std::vector result; - std::map> slicedOp; - std::vector slicedProbe; - auto &phi = envGet(Field, par().field); - - envGetTmp(ComplexField, ftBuf); - envGetTmp(Field, ftMatBuf); - dMask[0] = 0; - for (unsigned int mu = 1; mu < nd; ++mu) - { - partVol *= env().getDim(mu); - } - timeMom.resize(par().timeMom.size()); - for (unsigned int p = 0; p < timeMom.size(); ++p) - { - for (auto &tms: par().timeMom[p]) - { - std::vector tm = strToVec(tms); - - timeMom[p].push_back(tm); - timeMomSet.insert(tm); - } - transferMom.push_back(std::vector(nd - 1, 0)); - for (auto &tm: timeMom[p]) - { - for (unsigned int j = 1; j < nd; ++j) - { - transferMom[p][j - 1] -= tm[j]; - } - } - LOG(Message) << "Probe " << p << " (" << timeMom[p].size() << " points) : " << std::endl; - LOG(Message) << " phi(t_i, p_i) for (t_i, p_i) in " << timeMom[p] << std::endl; - LOG(Message) << " operator with momentum " << transferMom[p] << std::endl; - } - LOG(Message) << "FFT: field '" << par().field << "'" << std::endl; - fft.FFT_dim_mask(ftMatBuf, phi, dMask, FFT::forward); - slicedProbe.resize(timeMom.size()); - for (unsigned int p = 0; p < timeMom.size(); ++p) - { - std::vector qt; - - LOG(Message) << "Making probe " << p << std::endl; - slicedProbe[p].resize(nt); - for (unsigned int t = 0; t < nt; ++t) - { - Site acc; - - for (unsigned int i = 0; i < timeMom[p].size(); ++i) - { - Site buf; - - qt = timeMom[p][i]; - qt[0] += t; - vectorModulo(qt); - peekSite(buf, ftMatBuf, qt); - if (i == 0) - { - acc = buf; - } - else - { - acc *= buf; - } - } - slicedProbe[p][t] = TensorRemove(trace(acc)); - } - //std::cout << slicedProbe[p]<< std::endl; - } - for (auto &o: par().op) - { - auto &op = envGet(ComplexField, o); - - slicedOp[o].resize(transferMom.size()); - LOG(Message) << "FFT: operator '" << o << "'" << std::endl; - fft.FFT_dim_mask(ftBuf, op, dMask, FFT::forward); - //std::cout << ftBuf << std::endl; - for (unsigned int p = 0; p < transferMom.size(); ++p) - { - std::vector qt(nd, 0); - - for (unsigned int j = 1; j < nd; ++j) - { - qt[j] = transferMom[p][j - 1]; - } - slicedOp[o][p].resize(nt); - for (unsigned int t = 0; t < nt; ++t) - { - TComplex buf; - - qt[0] = t; - vectorModulo(qt); - peekSite(buf, ftBuf, qt); - slicedOp[o][p][t] = TensorRemove(buf); - } - //std::cout << ftBuf << std::endl; - //std::cout << slicedOp[o][p] << std::endl; - } - } - LOG(Message) << "Making correlators" << std::endl; - for (auto &o: par().op) - for (unsigned int p = 0; p < timeMom.size(); ++p) - { - TimeMomProbeResult r; - - LOG(Message) << " <" << o << " probe_" << p << ">" << std::endl; - r.op = o; - r.timeMom = timeMom[p]; - r.data = makeTwoPoint(slicedOp[o][p], slicedProbe[p], 1./partVol); - result.push_back(r); - } - saveResult(par().output, "timemomprobe", result); -} - -END_MODULE_NAMESPACE - -END_HADRONS_NAMESPACE - -#endif // Hadrons_MScalarSUN_TimeMomProbe_hpp_ diff --git a/Hadrons/Modules/MScalarSUN/TrMag.hpp b/Hadrons/Modules/MScalarSUN/TrMag.hpp index d7d40e88..b9602be3 100644 --- a/Hadrons/Modules/MScalarSUN/TrMag.hpp +++ b/Hadrons/Modules/MScalarSUN/TrMag.hpp @@ -124,7 +124,8 @@ void TTrMag::execute(void) std::vector result; auto &phi = envGet(Field, par().field); - auto m2 = sum(phi), mn = m2; + auto m2 = sum(phi); + auto mn = m2; m2 = -m2*m2; mn = 1.; diff --git a/Hadrons/Modules/MScalarSUN/Utils.hpp b/Hadrons/Modules/MScalarSUN/Utils.hpp index fdc42c9a..7eba5900 100644 --- a/Hadrons/Modules/MScalarSUN/Utils.hpp +++ b/Hadrons/Modules/MScalarSUN/Utils.hpp @@ -103,7 +103,7 @@ std::vector makeTwoPoint(const std::vector &sink, { for (unsigned int t = 0; t < nt; ++t) { - res[dt] += trace(sink[(t+dt)%nt]*source[t]); + res[dt] += trace(sink[(t+dt)%nt]*adj(source[t])); } res[dt] *= factor/static_cast(nt); } diff --git a/Hadrons/modules.inc b/Hadrons/modules.inc index 866dfd91..74138238 100644 --- a/Hadrons/modules.inc +++ b/Hadrons/modules.inc @@ -46,7 +46,6 @@ modules_cc =\ Modules/MAction/ScaledDWF.cc \ Modules/MScalarSUN/TrPhi.cc \ Modules/MScalarSUN/Grad.cc \ - Modules/MScalarSUN/TimeMomProbe.cc \ Modules/MScalarSUN/TrMag.cc \ Modules/MScalarSUN/TrKinetic.cc \ Modules/MScalarSUN/EMT.cc \ @@ -115,7 +114,6 @@ modules_hpp =\ Modules/MScalarSUN/TwoPointNPR.hpp \ Modules/MScalarSUN/ShiftProbe.hpp \ Modules/MScalarSUN/Div.hpp \ - Modules/MScalarSUN/TimeMomProbe.hpp \ Modules/MScalarSUN/TrMag.hpp \ Modules/MScalarSUN/EMT.hpp \ Modules/MScalarSUN/TwoPoint.hpp \ From 936eaac8e1c1d802bdac320933e11ad12437303b Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Mon, 8 Oct 2018 19:00:50 +0100 Subject: [PATCH 02/37] function to get the sha256 string --- Grid/lattice/Lattice_rng.h | 6 +----- Grid/util/Sha.h | 16 +++++++++++++++- 2 files changed, 16 insertions(+), 6 deletions(-) diff --git a/Grid/lattice/Lattice_rng.h b/Grid/lattice/Lattice_rng.h index c6b9e14e..5348538c 100644 --- a/Grid/lattice/Lattice_rng.h +++ b/Grid/lattice/Lattice_rng.h @@ -392,14 +392,10 @@ namespace Grid { void SeedUniqueString(const std::string &s){ std::vector seeds; - std::stringstream sha; seeds = GridChecksum::sha256_seeds(s); - for(int i=0;i &seeds){ diff --git a/Grid/util/Sha.h b/Grid/util/Sha.h index 6cfbe0bd..b0a8cc10 100644 --- a/Grid/util/Sha.h +++ b/Grid/util/Sha.h @@ -38,7 +38,21 @@ public: { return ::crc32(0L,(unsigned char *)data,bytes); } - static inline std::vector sha256(void *data,size_t bytes) + template + static inline std::string sha256_string(const std::vector &hash) + { + std::stringstream sha; + std::string s; + + for(unsigned int i = 0; i < hash.size(); i++) + { + sha << std::hex << static_cast(hash[i]); + } + s = sha.str(); + + return s; + } + static inline std::vector sha256(const void *data,size_t bytes) { std::vector hash(SHA256_DIGEST_LENGTH); SHA256_CTX sha256; From efc0c65056135be0cadd13a24086c23ab44317dc Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Mon, 8 Oct 2018 19:02:00 +0100 Subject: [PATCH 03/37] Hadrons: DiskVector Eigen specialisation with binary I/O and sha256 correctness check --- Hadrons/DiskVector.hpp | 115 ++++++++++++++++++++++++++----- tests/hadrons/Test_diskvector.cc | 16 +++++ 2 files changed, 113 insertions(+), 18 deletions(-) diff --git a/Hadrons/DiskVector.hpp b/Hadrons/DiskVector.hpp index 231854f3..94c3e597 100644 --- a/Hadrons/DiskVector.hpp +++ b/Hadrons/DiskVector.hpp @@ -34,6 +34,12 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include +#ifdef DV_DEBUG +#define DV_DEBUG_MSG(dv, stream) LOG(Debug) << "diskvector " << (dv) << ": " << stream << std::endl +#else +#define DV_DEBUG_MSG(dv, stream) +#endif + BEGIN_HADRONS_NAMESPACE /****************************************************************************** @@ -56,9 +62,7 @@ public: // write to disk and cache T &operator=(const T &obj) const { -#ifdef DV_DEBUG - LOG(Debug) << "diskvector " << &master_ << ": writing to " << i_ << std::endl; -#endif + DV_DEBUG_MSG(&master_, "writing to " << i_); master_.cacheInsert(i_, obj); master_.save(master_.filename(i_), obj); @@ -82,6 +86,8 @@ public: virtual ~DiskVectorBase(void); const T & operator[](const unsigned int i) const; RwAccessHelper operator[](const unsigned int i); + double hitRatio(void) const; + void resetStat(void); private: virtual void load(T &obj, const std::string filename) const = 0; virtual void save(const std::string filename, const T &obj) const = 0; @@ -93,6 +99,7 @@ private: private: std::string dirname_; unsigned int size_, cacheSize_; + double access_{0.}, hit_{0.}; bool clean_; // using pointers to allow modifications when class is const // semantic: const means data unmodified, but cache modification allowed @@ -115,6 +122,7 @@ private: read(reader, basename(filename), obj); } + virtual void save(const std::string filename, const T &obj) const { Writer writer(filename); @@ -123,13 +131,70 @@ private: } }; +/****************************************************************************** + * Specialisation for Eigen matrices * + ******************************************************************************/ +template +using EigenDiskVectorMat = Eigen::Matrix; + +template +class EigenDiskVector: public DiskVectorBase> +{ +public: + using DiskVectorBase>::DiskVectorBase; + typedef EigenDiskVectorMat Matrix; +private: + virtual void load(EigenDiskVectorMat &obj, const std::string filename) const + { + std::ifstream f(filename, std::ios::binary); + std::vector hash(SHA256_DIGEST_LENGTH); + Eigen::Index nRow, nCol; + size_t matSize; + double t; + + f.read(reinterpret_cast(hash.data()), hash.size()*sizeof(unsigned char)); + f.read(reinterpret_cast(&nRow), sizeof(Eigen::Index)); + f.read(reinterpret_cast(&nCol), sizeof(Eigen::Index)); + obj.resize(nRow, nCol); + matSize = nRow*nCol*sizeof(T); + t = -usecond(); + f.read(reinterpret_cast(obj.data()), matSize); + t += usecond(); + DV_DEBUG_MSG(this, "Eigen read " << matSize/t*1.0e6/1024/1024 << " MB/s"); + auto check = GridChecksum::sha256(obj.data(), matSize); + DV_DEBUG_MSG(this, "Eigen sha256 " << GridChecksum::sha256_string(check)); + if (hash != check) + { + HADRONS_ERROR(Io, "checksum failed") + } + } + + virtual void save(const std::string filename, const EigenDiskVectorMat &obj) const + { + std::ofstream f(filename, std::ios::binary); + std::vector hash(SHA256_DIGEST_LENGTH); + Eigen::Index nRow, nCol; + size_t matSize; + double t; + + nRow = obj.rows(); + nCol = obj.cols(); + matSize = nRow*nCol*sizeof(T); + hash = GridChecksum::sha256(obj.data(), matSize); + DV_DEBUG_MSG(this, "Eigen sha256 " << GridChecksum::sha256_string(hash)); + f.write(reinterpret_cast(hash.data()), hash.size()*sizeof(unsigned char)); + f.write(reinterpret_cast(&nRow), sizeof(Eigen::Index)); + f.write(reinterpret_cast(&nCol), sizeof(Eigen::Index)); + t = -usecond(); + f.write(reinterpret_cast(obj.data()), matSize); + t += usecond(); + DV_DEBUG_MSG(this, "Eigen write " << matSize/t*1.0e6/1024/1024 << " MB/s"); + } +}; + /****************************************************************************** * DiskVectorBase implementation * ******************************************************************************/ -#ifdef DV_DEBUG -#define DV_DEBUG_MSG(stream) LOG(Debug) << "diskvector " << this << ": " << stream << std::endl -#endif - template DiskVectorBase::DiskVectorBase(const std::string dirname, const unsigned int size, @@ -160,28 +225,29 @@ DiskVectorBase::~DiskVectorBase(void) template const T & DiskVectorBase::operator[](const unsigned int i) const { - auto &cache = *cachePtr_; - auto &loads = *loadsPtr_; + auto &cache = *cachePtr_; + auto &loads = *loadsPtr_; - DV_DEBUG_MSG("accessing " << i << " (RO)"); + DV_DEBUG_MSG(this, "accessing " << i << " (RO)"); if (i >= size_) { HADRONS_ERROR(Size, "index out of range"); } - + const_cast(access_)++; if (cache.find(i) == cache.end()) { // cache miss - DV_DEBUG_MSG("cache miss"); + DV_DEBUG_MSG(this, "cache miss"); fetch(i); } else { - DV_DEBUG_MSG("cache hit"); + DV_DEBUG_MSG(this, "cache hit"); auto pos = std::find(loads.begin(), loads.end(), i); + const_cast(hit_)++; loads.erase(pos); loads.push_back(i); } @@ -193,7 +259,7 @@ const T & DiskVectorBase::operator[](const unsigned int i) const { msg += std::to_string(p) + " "; } - DV_DEBUG_MSG("in cache: " << msg); + DV_DEBUG_MSG(this, "in cache: " << msg); #endif return cache.at(i); @@ -202,7 +268,7 @@ const T & DiskVectorBase::operator[](const unsigned int i) const template typename DiskVectorBase::RwAccessHelper DiskVectorBase::operator[](const unsigned int i) { - DV_DEBUG_MSG("accessing " << i << " (RW)"); + DV_DEBUG_MSG(this, "accessing " << i << " (RW)"); if (i >= size_) { @@ -212,6 +278,19 @@ typename DiskVectorBase::RwAccessHelper DiskVectorBase::operator[](const u return RwAccessHelper(*this, i); } +template +double DiskVectorBase::hitRatio(void) const +{ + return hit_/access_; +} + +template +void DiskVectorBase::resetStat(void) +{ + access_ = 0.; + hit_ = 0.; +} + template std::string DiskVectorBase::filename(const unsigned int i) const { @@ -226,7 +305,7 @@ void DiskVectorBase::evict(void) const if (cache.size() >= cacheSize_) { - DV_DEBUG_MSG("evicting " << loads.front()); + DV_DEBUG_MSG(this, "evicting " << loads.front()); cache.erase(loads.front()); loads.pop_front(); } @@ -239,7 +318,7 @@ void DiskVectorBase::fetch(const unsigned int i) const auto &loads = *loadsPtr_; struct stat s; - DV_DEBUG_MSG("loading " << i << " from disk"); + DV_DEBUG_MSG(this, "loading " << i << " from disk"); evict(); if(stat(filename(i).c_str(), &s) != 0) @@ -267,7 +346,7 @@ void DiskVectorBase::cacheInsert(const unsigned int i, const T &obj) const { msg += std::to_string(p) + " "; } - DV_DEBUG_MSG("in cache: " << msg); + DV_DEBUG_MSG(this, "in cache: " << msg); #endif } diff --git a/tests/hadrons/Test_diskvector.cc b/tests/hadrons/Test_diskvector.cc index 363ae2ce..10bc4db1 100644 --- a/tests/hadrons/Test_diskvector.cc +++ b/tests/hadrons/Test_diskvector.cc @@ -91,6 +91,22 @@ int main(int argc, char *argv[]) v13r = v[13]; LOG(Message) << "v[13] correct? " << ((v13r == v13w) ? "yes" : "no" ) << std::endl; + LOG(Message) << "hit ratio " << v.hitRatio() << std::endl; + + EigenDiskVector w("eigendiskvector_test", 1000, 4); + EigenDiskVector::Matrix m,n; + + w[2] = EigenDiskVectorMat::Random(2000, 2000); + m = w[2]; + w[3] = EigenDiskVectorMat::Random(2000, 2000); + w[4] = EigenDiskVectorMat::Random(2000, 2000); + w[5] = EigenDiskVectorMat::Random(2000, 2000); + w[6] = EigenDiskVectorMat::Random(2000, 2000); + w[7] = EigenDiskVectorMat::Random(2000, 2000); + n = w[2]; + LOG(Message) << "w[2] correct? " + << ((m == n) ? "yes" : "no" ) << std::endl; + LOG(Message) << "hit ratio " << w.hitRatio() << std::endl; Grid_finalize(); From 49f25e08e8e7c5003c8507a1242ff644ab9a0630 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 11 Oct 2018 12:35:32 +0100 Subject: [PATCH 04/37] PauliVillars based 4D -> 5D reconstruction with Fourier Accelerated PV inverse by Christoph. Differs from the one by Rudy in BFM since it vectorises the twisted 4D solves in pairs. --- Grid/qcd/action/fermion/CayleyFermion5D.cc | 20 ++ Grid/qcd/action/fermion/CayleyFermion5D.h | 8 + Grid/qcd/action/fermion/Fermion.h | 13 +- Grid/qcd/action/fermion/FermionOperatorImpl.h | 1 + .../qcd/action/fermion/FourierAcceleratedPV.h | 237 ++++++++++++++++++ .../action/fermion/PauliVillarsInverters.h | 95 +++++++ Grid/qcd/action/fermion/Reconstruct5Dprop.h | 135 ++++++++++ Grid/qcd/action/fermion/WilsonTMFermion5D.h | 155 ++++++++++++ tests/debug/Test_cayley_cg.cc | 169 ++++++++++++- 9 files changed, 823 insertions(+), 10 deletions(-) create mode 100644 Grid/qcd/action/fermion/FourierAcceleratedPV.h create mode 100644 Grid/qcd/action/fermion/PauliVillarsInverters.h create mode 100644 Grid/qcd/action/fermion/Reconstruct5Dprop.h create mode 100644 Grid/qcd/action/fermion/WilsonTMFermion5D.h diff --git a/Grid/qcd/action/fermion/CayleyFermion5D.cc b/Grid/qcd/action/fermion/CayleyFermion5D.cc index d4ed5a7c..0e0e7bee 100644 --- a/Grid/qcd/action/fermion/CayleyFermion5D.cc +++ b/Grid/qcd/action/fermion/CayleyFermion5D.cc @@ -68,6 +68,26 @@ void CayleyFermion5D::ExportPhysicalFermionSolution(const FermionField &so ExtractSlice(exported4d, tmp, 0, 0); } template +void CayleyFermion5D::P(const FermionField &psi, FermionField &chi) +{ + int Ls= this->Ls; + chi=zero; + for(int s=0;s +void CayleyFermion5D::Pdag(const FermionField &psi, FermionField &chi) +{ + int Ls= this->Ls; + chi=zero; + for(int s=0;s void CayleyFermion5D::ExportPhysicalFermionSource(const FermionField &solution5d,FermionField &exported4d) { int Ls = this->Ls; diff --git a/Grid/qcd/action/fermion/CayleyFermion5D.h b/Grid/qcd/action/fermion/CayleyFermion5D.h index 3bf2a8b6..a84c0e8b 100644 --- a/Grid/qcd/action/fermion/CayleyFermion5D.h +++ b/Grid/qcd/action/fermion/CayleyFermion5D.h @@ -93,6 +93,14 @@ namespace Grid { virtual void ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d); virtual void ImportUnphysicalFermion(const FermionField &solution5d, FermionField &exported4d); + /////////////////////////////////////////////////////////////// + // Support for MADWF tricks + /////////////////////////////////////////////////////////////// + RealD Mass(void) { return mass; }; + void SetMass(RealD _mass) { mass=_mass; } ; + void P(const FermionField &psi, FermionField &chi); + void Pdag(const FermionField &psi, FermionField &chi); + ///////////////////////////////////////////////////// // Instantiate different versions depending on Impl ///////////////////////////////////////////////////// diff --git a/Grid/qcd/action/fermion/Fermion.h b/Grid/qcd/action/fermion/Fermion.h index 2a008cb7..d0200ae7 100644 --- a/Grid/qcd/action/fermion/Fermion.h +++ b/Grid/qcd/action/fermion/Fermion.h @@ -80,12 +80,23 @@ Author: Peter Boyle /////////////////////////////////////////////////////////////////////////////// #include +/////////////////////////////////////////////////////////////////////////////// +// Fourier accelerated Pauli Villars inverse support +/////////////////////////////////////////////////////////////////////////////// +#include + +//////////////////////////////////////////////////////////////////////////////// +// Move this group to a DWF specific tools/algorithms subdir? +//////////////////////////////////////////////////////////////////////////////// +#include +#include +#include + //////////////////////////////////////////////////////////////////////////////////////////////////// // More maintainable to maintain the following typedef list centrally, as more "impl" targets // are added, (e.g. extension for gparity, half precision project in comms etc..) //////////////////////////////////////////////////////////////////////////////////////////////////// - // Cayley 5d namespace Grid { namespace QCD { diff --git a/Grid/qcd/action/fermion/FermionOperatorImpl.h b/Grid/qcd/action/fermion/FermionOperatorImpl.h index e6f6e136..b0ffa90e 100644 --- a/Grid/qcd/action/fermion/FermionOperatorImpl.h +++ b/Grid/qcd/action/fermion/FermionOperatorImpl.h @@ -141,6 +141,7 @@ namespace QCD { //////////////////////////////////////////////////////////////////////// #define INHERIT_FIMPL_TYPES(Impl)\ + typedef Impl Impl_t; \ typedef typename Impl::FermionField FermionField; \ typedef typename Impl::PropagatorField PropagatorField; \ typedef typename Impl::DoubledGaugeField DoubledGaugeField; \ diff --git a/Grid/qcd/action/fermion/FourierAcceleratedPV.h b/Grid/qcd/action/fermion/FourierAcceleratedPV.h new file mode 100644 index 00000000..d8441cf2 --- /dev/null +++ b/Grid/qcd/action/fermion/FourierAcceleratedPV.h @@ -0,0 +1,237 @@ + + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/qcd/action/fermion/FourierAcceleratedPV.h + + Copyright (C) 2015 + +Author: Christoph Lehner +Author: Peter Boyle + + 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 */ +#pragma once +namespace Grid { +namespace QCD { + + template + void get_real_const_bc(M& m, RealD& _b, RealD& _c) { + ComplexD b,c; + b=m.bs[0]; + c=m.cs[0]; + std::cout << GridLogMessage << "b=" << b << ", c=" << c << std::endl; + for (size_t i=1;i +class FourierAcceleratedPV { + public: + + ConjugateGradient &cg; + M& dwfPV; + G& Umu; + GridCartesian* grid5D; + GridRedBlackCartesian* gridRB5D; + int group_in_s; + + FourierAcceleratedPV(M& _dwfPV, G& _Umu, ConjugateGradient &_cg, int _group_in_s = 2) + : dwfPV(_dwfPV), Umu(_Umu), cg(_cg), group_in_s(_group_in_s) + { + assert( dwfPV.FermionGrid()->_fdimensions[0] % (2*group_in_s) == 0); + grid5D = QCD::SpaceTimeGrid::makeFiveDimGrid(2*group_in_s, (GridCartesian*)Umu._grid); + gridRB5D = QCD::SpaceTimeGrid::makeFiveDimRedBlackGrid(2*group_in_s, (GridCartesian*)Umu._grid); + } + + void rotatePV(const Vi& _src, Vi& dst, bool forward) const { + + GridStopWatch gsw1, gsw2; + + typedef typename Vi::scalar_type Coeff_t; + int Ls = dst._grid->_fdimensions[0]; + + Vi _tmp(dst._grid); + double phase = M_PI / (double)Ls; + Coeff_t bzero(0.0,0.0); + + FFT theFFT((GridCartesian*)dst._grid); + + if (!forward) { + gsw1.Start(); + for (int s=0;s_fdimensions[0]; + + GridStopWatch gswT; + gswT.Start(); + + RealD b,c; + get_real_const_bc(dwfPV,b,c); + RealD M5 = dwfPV.M5; + + // U(true) Rightinv TMinv U(false) = Minv + + Vi _src_diag(_dst._grid); + Vi _src_diag_slice(dwfPV.GaugeGrid()); + Vi _dst_diag_slice(dwfPV.GaugeGrid()); + Vi _src_diag_slices(grid5D); + Vi _dst_diag_slices(grid5D); + Vi _dst_diag(_dst._grid); + + rotatePV(_src,_src_diag,false); + + // now do TM solves + Gamma G5(Gamma::Algebra::Gamma5); + + GridStopWatch gswA, gswB; + + gswA.Start(); + + typedef typename M::Impl_t Impl; + //WilsonTMFermion tm(x.Umu,*x.UGridF,*x.UrbGridF,0.0,0.0,solver_outer.parent.par.wparams_f); + std::vector vmass(grid5D->_fdimensions[0],0.0); + std::vector vmu(grid5D->_fdimensions[0],0.0); + + WilsonTMFermion5D tm(Umu,*grid5D,*gridRB5D, + *(GridCartesian*)dwfPV.GaugeGrid(), + *(GridRedBlackCartesian*)dwfPV.GaugeRedBlackGrid(), + vmass,vmu); + + //SchurRedBlackDiagTwoSolve sol(cg); + SchurRedBlackDiagMooeeSolve sol(cg); // same performance as DiagTwo + gswA.Stop(); + + gswB.Start(); + + for (int sgroup=0;sgroup + + 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 */ +#pragma once + +namespace Grid { +namespace QCD { + +template +class PauliVillarsSolverUnprec +{ + public: + ConjugateGradient & CG; + PauliVillarsSolverUnprec( ConjugateGradient &_CG) : CG(_CG){}; + + template + void operator() (Matrix &_Matrix,const Field &src,Field &sol) + { + RealD m = _Matrix.Mass(); + Field A (_Matrix.FermionGrid()); + + MdagMLinearOperator HermOp(_Matrix); + + _Matrix.SetMass(1.0); + _Matrix.Mdag(src,A); + CG(HermOp,A,sol); + _Matrix.SetMass(m); + }; +}; + +template +class PauliVillarsSolverRBprec +{ + public: + SchurSolverType & SchurSolver; + PauliVillarsSolverRBprec( SchurSolverType &_SchurSolver) : SchurSolver(_SchurSolver){}; + + template + void operator() (Matrix &_Matrix,const Field &src,Field &sol) + { + RealD m = _Matrix.Mass(); + Field A (_Matrix.FermionGrid()); + + _Matrix.SetMass(1.0); + SchurSolver(_Matrix,src,sol); + _Matrix.SetMass(m); + }; +}; + +template +class PauliVillarsSolverFourierAccel +{ + public: + GaugeField & Umu; + ConjugateGradient & CG; + + PauliVillarsSolverFourierAccel(GaugeField &_Umu,ConjugateGradient &_CG) : Umu(_Umu), CG(_CG) + { + }; + + template + void operator() (Matrix &_Matrix,const Field &src,Field &sol) + { + FourierAcceleratedPV faPV(_Matrix,Umu,CG) ; + faPV.pvInv(src,sol); + }; +}; + + +} +} diff --git a/Grid/qcd/action/fermion/Reconstruct5Dprop.h b/Grid/qcd/action/fermion/Reconstruct5Dprop.h new file mode 100644 index 00000000..6862c5ee --- /dev/null +++ b/Grid/qcd/action/fermion/Reconstruct5Dprop.h @@ -0,0 +1,135 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/algorithms/iterative/SchurRedBlack.h + + Copyright (C) 2015 + +Author: Peter Boyle + + 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 */ +#pragma once + +namespace Grid { +namespace QCD { + +template class Reconstruct5DfromPhysical { + private: + PVinverter & PauliVillarsSolver; + public: + + ///////////////////////////////////////////////////// + // First cut works, 10 Oct 2018. + // + // Must form a plan to get this into production for Zmobius acceleration + // of the Mobius exact AMA corrections. + // + // TODO : understand absence of contact term in eqns in Hantao's thesis + // sol4 is contact term subtracted, but thesis & Brower's paper suggests not. + // + // Step 1: Localise PV inverse in a routine. [DONE] + // Step 2: Schur based PV inverse [DONE] + // Step 3: Fourier accelerated PV inverse [DONE] + // + ///////////////////////////////////////////////////// + + Reconstruct5DfromPhysical(PVinverter &_PauliVillarsSolver) + : PauliVillarsSolver(_PauliVillarsSolver) + { + }; + + + template + void PV(Matrix &_Matrix,const Field &src,Field &sol) + { + RealD m = _Matrix.Mass(); + _Matrix.SetMass(1.0); + _Matrix.M(src,sol); + _Matrix.SetMass(m); + } + template + void PVdag(Matrix &_Matrix,const Field &src,Field &sol) + { + RealD m = _Matrix.Mass(); + _Matrix.SetMass(1.0); + _Matrix.Mdag(src,sol); + _Matrix.SetMass(m); + } + template + void operator() (Matrix & _Matrix,const Field &sol4,const Field &src4, Field &sol5){ + + int Ls = _Matrix.Ls; + + Field psi4(_Matrix.GaugeGrid()); + Field psi(_Matrix.FermionGrid()); + Field A (_Matrix.FermionGrid()); + Field B (_Matrix.FermionGrid()); + Field c (_Matrix.FermionGrid()); + + typedef typename Matrix::Coeff_t Coeff_t; + + std::cout << GridLogMessage<< " ************************************************" << std::endl; + std::cout << GridLogMessage<< " Reconstruct5Dprop: c.f. MADWF algorithm " << std::endl; + std::cout << GridLogMessage<< " ************************************************" << std::endl; + + /////////////////////////////////////// + //Import source, include Dminus factors + /////////////////////////////////////// + _Matrix.ImportPhysicalFermionSource(src4,B); + + /////////////////////////////////////// + // Set up c from src4 + /////////////////////////////////////// + PauliVillarsSolver(_Matrix,B,A); + _Matrix.Pdag(A,c); + + ////////////////////////////////////// + // Build Pdag PV^-1 Dm P [-sol4,c2,c3... cL] + ////////////////////////////////////// + psi4 = - sol4; + InsertSlice(psi4, psi, 0 , 0); + for (int s=1;s ; NB Christoph did similar in GPT + + 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 */ +#pragma once + +#include +#include + + +namespace Grid { + + namespace QCD { + + template + class WilsonTMFermion5D : public WilsonFermion5D + { + public: + INHERIT_IMPL_TYPES(Impl); + public: + + virtual void Instantiatable(void) {}; + + // Constructors + WilsonTMFermion5D(GaugeField &_Umu, + GridCartesian &Fgrid, + GridRedBlackCartesian &Frbgrid, + GridCartesian &Ugrid, + GridRedBlackCartesian &Urbgrid, + const std::vector _mass, + const std::vector _mu, + const ImplParams &p= ImplParams() + ) : + WilsonFermion5D(_Umu, + Fgrid, + Frbgrid, + Ugrid, + Urbgrid, + 4.0,p) + + { + update(_mass,_mu); + } + + virtual void Meooe(const FermionField &in, FermionField &out) { + if (in.checkerboard == Odd) { + this->DhopEO(in, out, DaggerNo); + } else { + this->DhopOE(in, out, DaggerNo); + } + } + + virtual void MeooeDag(const FermionField &in, FermionField &out) { + if (in.checkerboard == Odd) { + this->DhopEO(in, out, DaggerYes); + } else { + this->DhopOE(in, out, DaggerYes); + } + } + + // allow override for twisted mass and clover + virtual void Mooee(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + //axpibg5x(out,in,a,b); // out = a*in + b*i*G5*in + for (int s=0;s<(int)this->mass.size();s++) { + ComplexD a = 4.0+this->mass[s]; + ComplexD b(0.0,this->mu[s]); + axpbg5y_ssp(out,a,in,b,in,s,s); + } + } + + virtual void MooeeDag(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + for (int s=0;s<(int)this->mass.size();s++) { + ComplexD a = 4.0+this->mass[s]; + ComplexD b(0.0,-this->mu[s]); + axpbg5y_ssp(out,a,in,b,in,s,s); + } + } + virtual void MooeeInv(const FermionField &in, FermionField &out) { + for (int s=0;s<(int)this->mass.size();s++) { + RealD m = this->mass[s]; + RealD tm = this->mu[s]; + RealD mtil = 4.0+this->mass[s]; + RealD sq = mtil*mtil+tm*tm; + ComplexD a = mtil/sq; + ComplexD b(0.0, -tm /sq); + axpbg5y_ssp(out,a,in,b,in,s,s); + } + } + virtual void MooeeInvDag(const FermionField &in, FermionField &out) { + for (int s=0;s<(int)this->mass.size();s++) { + RealD m = this->mass[s]; + RealD tm = this->mu[s]; + RealD mtil = 4.0+this->mass[s]; + RealD sq = mtil*mtil+tm*tm; + ComplexD a = mtil/sq; + ComplexD b(0.0,tm /sq); + axpbg5y_ssp(out,a,in,b,in,s,s); + } + } + + virtual RealD M(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + this->Dhop(in, out, DaggerNo); + FermionField tmp(out._grid); + for (int s=0;s<(int)this->mass.size();s++) { + ComplexD a = 4.0+this->mass[s]; + ComplexD b(0.0,this->mu[s]); + axpbg5y_ssp(tmp,a,in,b,in,s,s); + } + return axpy_norm(out, 1.0, tmp, out); + } + + // needed for fast PV + void update(const std::vector& _mass, const std::vector& _mu) { + assert(_mass.size() == _mu.size()); + assert(_mass.size() == this->FermionGrid()->_fdimensions[0]); + this->mass = _mass; + this->mu = _mu; + } + + private: + std::vector mu; + std::vector mass; + + }; + + typedef WilsonTMFermion5D WilsonTMFermion5DF; + typedef WilsonTMFermion5D WilsonTMFermion5DD; + +}} diff --git a/tests/debug/Test_cayley_cg.cc b/tests/debug/Test_cayley_cg.cc index 57fe7717..fbc978a6 100644 --- a/tests/debug/Test_cayley_cg.cc +++ b/tests/debug/Test_cayley_cg.cc @@ -1,5 +1,4 @@ /************************************************************************************* - Grid physics library, www.github.com/paboyle/Grid Source file: ./tests/Test_cayley_cg.cc @@ -27,6 +26,7 @@ Author: paboyle *************************************************************************************/ /* END LEGAL */ #include +#include using namespace std; using namespace Grid; @@ -46,6 +46,7 @@ struct scal { template void TestCGinversions(What & Ddwf, + LatticeGaugeField &Umu, GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, RealD mass, RealD M5, @@ -75,6 +76,24 @@ void TestCGprec(What & Ddwf, GridParallelRNG *RNG4, GridParallelRNG *RNG5); +template +void TestReconstruct5D(What & Ddwf, + LatticeGaugeField &Umu, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5); + +template +void TestReconstruct5DFA(What & Ddwf, + LatticeGaugeField &Umu, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5); + int main (int argc, char ** argv) { Grid_init(&argc,&argv); @@ -100,46 +119,71 @@ int main (int argc, char ** argv) RealD mass=0.1; RealD M5 =1.8; + std::cout<(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestCGinversions(Ddwf,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5DFA(Ddwf,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); RealD b=1.5;// Scale factor b+c=2, b-c=1 RealD c=0.5; std::vector gamma(Ls,ComplexD(1.0,0.0)); + std::cout<(Dmob,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestCGinversions(Dmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5DFA(Dmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + std::cout<(ZDmob,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestCGinversions(ZDmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5D(ZDmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + std::cout<(Dzolo,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestCGinversions(Dzolo,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5D(Dzolo,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + std::cout<(Dsham,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestCGinversions(Dsham,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5DFA(Dsham,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + std::cout<(Dshamz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestCGinversions(Dshamz,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5D(Dshamz,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + std::cout<(Dov,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestCGinversions(Dov,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5DFA(Dov,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + std::cout<(Dovz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestCGinversions(Dovz,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5D(Dovz,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); Grid_finalize(); } template void TestCGinversions(What & Ddwf, + LatticeGaugeField &Umu, GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, RealD mass, RealD M5, @@ -154,6 +198,7 @@ void TestCGinversions(What & Ddwf, TestCGschur(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5); } + template void TestCGunprec(What & Ddwf, GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, @@ -189,6 +234,112 @@ void TestCGprec(What & Ddwf, CG(HermOpEO,src_o,result_o); } +template +void TestReconstruct5D(What & Ddwf, + LatticeGaugeField & Umu, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5) +{ + LatticeFermion src4 (UGrid); random(*RNG4,src4); + LatticeFermion res4 (UGrid); res4 = zero; + + LatticeFermion src (FGrid); + LatticeFermion src_NE(FGrid); + LatticeFermion result(FGrid); + LatticeFermion result_rec(FGrid); + + MdagMLinearOperator HermOp(Ddwf); + double Resid = 1.0e-12; + ConjugateGradient CG(Resid,10000); + + Ddwf.ImportPhysicalFermionSource(src4,src); + Ddwf.Mdag(src,src_NE); + CG(HermOp,src_NE,result); + + Ddwf.ExportPhysicalFermionSolution(result, res4); + + Ddwf.M(result,src_NE); + src_NE = src_NE - src; + std::cout < SchurSolverType; + typedef PauliVillarsSolverRBprec PVinverter; + SchurSolverType SchurSolver(CG); + PVinverter PVinverse(SchurSolver); + + Reconstruct5DfromPhysical reconstructor(PVinverse); + + reconstructor(Ddwf,res4,src4,result_rec); + + std::cout < +void TestReconstruct5DFA(What & Ddwf, + LatticeGaugeField & Umu, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5) +{ + LatticeFermion src4 (UGrid); random(*RNG4,src4); + LatticeFermion res4 (UGrid); res4 = zero; + + LatticeFermion src (FGrid); + LatticeFermion src_NE(FGrid); + LatticeFermion result(FGrid); + LatticeFermion result_rec(FGrid); + + MdagMLinearOperator HermOp(Ddwf); + double Resid = 1.0e-12; + ConjugateGradient CG(Resid,10000); + + Ddwf.ImportPhysicalFermionSource(src4,src); + Ddwf.Mdag(src,src_NE); + CG(HermOp,src_NE,result); + + Ddwf.ExportPhysicalFermionSolution(result, res4); + + Ddwf.M(result,src_NE); + src_NE = src_NE - src; + std::cout < PVinverter; + PVinverter PVinverse(Umu,CG); + + Reconstruct5DfromPhysical reconstructor(PVinverse); + + reconstructor(Ddwf,res4,src4,result_rec); + + std::cout < void TestCGschur(What & Ddwf, From 03c3d495a26facb4715d33fdcc6d2724309d8506 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Fri, 12 Oct 2018 10:59:33 +0100 Subject: [PATCH 05/37] First cut (non functional NPR code) developed by Julia Kettle --- Hadrons/Modules/MGauge/GaugeFix.cc | 93 +++++++++ Hadrons/Modules/MGauge/GaugeFix.hpp | 81 ++++++++ Hadrons/Modules/MNPR/Amputate.cc | 8 + Hadrons/Modules/MNPR/Amputate.hpp | 199 +++++++++++++++++++ Hadrons/Modules/MNPR/Bilinear.cc | 8 + Hadrons/Modules/MNPR/Bilinear.hpp | 224 ++++++++++++++++++++++ Hadrons/Modules/MNPR/FourQuark.cc | 8 + Hadrons/Modules/MNPR/FourQuark.hpp | 273 +++++++++++++++++++++++++++ Hadrons/Modules/MSource/Momentum.cc | 8 + Hadrons/Modules/MSource/Momentum.hpp | 121 ++++++++++++ 10 files changed, 1023 insertions(+) create mode 100644 Hadrons/Modules/MGauge/GaugeFix.cc create mode 100644 Hadrons/Modules/MGauge/GaugeFix.hpp create mode 100644 Hadrons/Modules/MNPR/Amputate.cc create mode 100644 Hadrons/Modules/MNPR/Amputate.hpp create mode 100644 Hadrons/Modules/MNPR/Bilinear.cc create mode 100644 Hadrons/Modules/MNPR/Bilinear.hpp create mode 100644 Hadrons/Modules/MNPR/FourQuark.cc create mode 100644 Hadrons/Modules/MNPR/FourQuark.hpp create mode 100644 Hadrons/Modules/MSource/Momentum.cc create mode 100644 Hadrons/Modules/MSource/Momentum.hpp diff --git a/Hadrons/Modules/MGauge/GaugeFix.cc b/Hadrons/Modules/MGauge/GaugeFix.cc new file mode 100644 index 00000000..bd57ed07 --- /dev/null +++ b/Hadrons/Modules/MGauge/GaugeFix.cc @@ -0,0 +1,93 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MGauge/Fix.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; + + +/****************************************************************************** +* TGaugeFix implementation * +******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TGaugeFix::TGaugeFix(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TGaugeFix::getInput(void) +{ + std::vector in = {par().gauge}; + return in; +} + +template +std::vector TGaugeFix::getOutput(void) +{ + std::vector out = {getName()}; + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TGaugeFix::setup(void) +{ + envCreateLat(GaugeField, getName()); +} + + +// execution /////////////////////////////////////////////////////////////////// +template +void TGaugeFix::execute(void) +//Loads the gauge and fixes it +{ + + std::cout << "executing" << std::endl; + LOG(Message) << "Fixing the Gauge" << std::endl; + LOG(Message) << par().gauge << std::endl; + auto &U = envGet(GaugeField, par().gauge); + auto &Umu = envGet(GaugeField, getName()); + LOG(Message) << "Gauge Field fetched" << std::endl; + //do we allow maxiter etc to be user set? + Real alpha = par().alpha; + int maxiter = par().maxiter; + Real Omega_tol = par().Omega_tol; + Real Phi_tol = par().Phi_tol; + bool Fourier = par().Fourier; + FourierAcceleratedGaugeFixer::SteepestDescentGaugeFix(U,alpha,maxiter,Omega_tol,Phi_tol,Fourier); + Umu = U; + LOG(Message) << "Gauge Fixed" << std::endl; + +} + +template class Grid::Hadrons::MGauge::TGaugeFix; diff --git a/Hadrons/Modules/MGauge/GaugeFix.hpp b/Hadrons/Modules/MGauge/GaugeFix.hpp new file mode 100644 index 00000000..154237a7 --- /dev/null +++ b/Hadrons/Modules/MGauge/GaugeFix.hpp @@ -0,0 +1,81 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MGauge/Fix.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_MGaugeFix_hpp_ +#define Hadrons_MGaugeFix_hpp_ + +#include +#include +#include +#include +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Fix gauge * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MGauge) + +class GaugeFixPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(GaugeFixPar, + std::string, gauge, + Real, alpha, + int, maxiter, + Real, Omega_tol, + Real, Phi_tol, + bool, Fourier); +}; + +template +class TGaugeFix: public Module +{ +public: + GAUGE_TYPE_ALIASES(GImpl,); +public: + // constructor + TGaugeFix(const std::string name); + // destructor + virtual ~TGaugeFix(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_TMP(GaugeFix, TGaugeFix, MGauge); + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MGaugeFix_hpp_ diff --git a/Hadrons/Modules/MNPR/Amputate.cc b/Hadrons/Modules/MNPR/Amputate.cc new file mode 100644 index 00000000..c4ada9d9 --- /dev/null +++ b/Hadrons/Modules/MNPR/Amputate.cc @@ -0,0 +1,8 @@ +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MNPR; + +template class Grid::Hadrons::MNPR::TAmputate; + diff --git a/Hadrons/Modules/MNPR/Amputate.hpp b/Hadrons/Modules/MNPR/Amputate.hpp new file mode 100644 index 00000000..797247b5 --- /dev/null +++ b/Hadrons/Modules/MNPR/Amputate.hpp @@ -0,0 +1,199 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MNpr/Amputate.hpp + +Copyright (C) 2015 +Copyright (C) 2016 + +Author: Julia Kettle J.R.Kettle-2@sms.ed.ac.uk + +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_Amputate_hpp_ +#define Hadrons_Amputate_hpp_ + +#include +#include +#include +#include +//#include +//#include +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * TAmputate * + Performs bilinear contractions of the type tr[g5*adj(Sout)*g5*G*Sin] + Suitable for non exceptional momenta +******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MNPR) + +class AmputatePar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(AmputatePar, + std::string, Sin, //need to make this a propogator type? + std::string, Sout, //same + std::string, vertex, + std::string, pin, + std::string, pout, + std::string, output, + std::string, input); +}; + +template +class TAmputate: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl1, 1); + FERM_TYPE_ALIASES(FImpl2, 2); + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::vector, Vamp, + ); + }; +public: + // constructor + TAmputate(const std::string name); + // destructor + virtual ~TAmputate(void) = default; + // dependencies/products + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + virtual SpinColourMatrix invertspincolmat(SpinColourMatrix &scmat); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_TMP(Amputate, ARG(TAmputate), MNPR); + +/****************************************************************************** + * TAmputate implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TAmputate::TAmputate(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TAmputate::getInput(void) +{ + std::vector input = {par().Sin, par().Sout, par().vertex}; + + return input; +} + +template +std::vector TAmputate::getOutput(void) +{ + std::vector output = {getName()}; + + + return output; +} + +// Invert spin colour matrix using Eigen +template +SpinColourMatrix TAmputate::invertspincolmat(SpinColourMatrix &scmat) +{ + Eigen::MatrixXcf scmat_2d(Ns*Nc,Ns*Nc); + for(int ic=0; ic +void TAmputate::execute(void) +{ + LOG(Message) << "Computing bilinear amputations '" << getName() << "' using" + << " momentum '" << par().Sin << "' and '" << par().Sout << "'" + << std::endl; + BinaryWriter writer(par().output); + PropagatorField1 &Sin = *env().template getObject(par().Sin); //Do these have the phases taken into account?? Don't think so. FIX + PropagatorField2 &Sout = *env().template getObject(par().Sout); + std::vector pin = strToVec(par().pin), pout = strToVec(par().pout); + std::vector latt_size(pin.begin(), pin.end()); + LatticeComplex pdotxin(env().getGrid()), pdotxout(env().getGrid()), coor(env().getGrid()); + LOG(Message) << "Propagators set up " << std::endl; + std::vector vertex; // Let's read from file here + Gamma g5(Gamma::Algebra::Gamma5); + Result result; + LOG(Message) << "reading file - " << par().input << std::endl; + BinaryReader reader(par().input); + Complex Ci(0.0,1.0); + + std::string svertex; + read(reader,"vertex", vertex); + LOG(Message) << "vertex read" << std::endl; + + pdotxin=zero; + pdotxout=zero; + for (unsigned int mu = 0; mu < 4; ++mu) + { + Real TwoPiL = M_PI * 2.0/ latt_size[mu]; + LatticeCoordinate(coor,mu); + pdotxin = pdotxin +(TwoPiL * pin[mu]) * coor; + pdotxout= pdotxout +(TwoPiL * pout[mu]) * coor; + } + Sin = Sin*exp(-Ci*pdotxin); //phase corrections + Sout = Sout*exp(-Ci*pdotxout); + + SpinColourMatrix Sin_mom = sum(Sin); + SpinColourMatrix Sout_mom = sum(Sout); + LOG(Message) << "summed over lattice" << std::endl; + + LOG(Message) << "Lattice -> spincolourmatrix conversion" << std::endl; + + SpinColourMatrix Sin_inv = invertspincolmat(Sin_mom); + SpinColourMatrix Sout_inv = invertspincolmat(Sout_mom); + LOG(Message) << "Inversions done" << std::endl; + + result.Vamp.resize(Gamma::nGamma/2); + for( int mu=0; mu < Gamma::nGamma/2; mu++){ + Gamma::Algebra gam = mu; + result.Vamp[mu] = 1/12.0*trace(adj(Gamma(mu*2+1))*g5*Sout_inv*g5*vertex[mu]*Sin_inv); + LOG(Message) << "Vamp[" << mu << "] - " << result.Vamp[mu] << std::endl; + } +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Amputate_hpp_ diff --git a/Hadrons/Modules/MNPR/Bilinear.cc b/Hadrons/Modules/MNPR/Bilinear.cc new file mode 100644 index 00000000..01fb08c6 --- /dev/null +++ b/Hadrons/Modules/MNPR/Bilinear.cc @@ -0,0 +1,8 @@ +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MNPR; + +template class Grid::Hadrons::MNPR::TBilinear; + diff --git a/Hadrons/Modules/MNPR/Bilinear.hpp b/Hadrons/Modules/MNPR/Bilinear.hpp new file mode 100644 index 00000000..0dc621b9 --- /dev/null +++ b/Hadrons/Modules/MNPR/Bilinear.hpp @@ -0,0 +1,224 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MContraction/Bilinear.hpp + +Copyright (C) 2015 +Copyright (C) 2016 + +Author: Julia Kettle J.R.Kettle-2@sms.ed.ac.uk + +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_Bilinear_hpp_ +#define Hadrons_Bilinear_hpp_ + +#include +#include +#include +#include +//#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * TBilinear * + Performs bilinear contractions of the type tr[g5*adj(Sout)*g5*G*Sin] + Suitable for non exceptional momenta in Rome-Southampton NPR +******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MNPR) + +class BilinearPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(BilinearPar, + std::string, Sin, + std::string, Sout, + std::string, pin, + std::string, pout, + std::string, output); +}; + +template +class TBilinear: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl1, 1); + FERM_TYPE_ALIASES(FImpl2, 2); + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::vector, bilinear); + }; +public: + // constructor + TBilinear(const std::string name); + // destructor + virtual ~TBilinear(void) = default; + // dependencies/products + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + //LatticeSpinColourMatrix PhaseProps(LatticeSpinColourMatrix S, std::vector p); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_TMP(Bilinear, ARG(TBilinear), MNPR); + +/****************************************************************************** + * TBilinear implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TBilinear::TBilinear(const std::string name) +: Module(name) +{} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TBilinear::setup(void) +{ + //env().template registerLattice(getName()); + //env().template registerObject(getName()); +} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TBilinear::getInput(void) +{ + std::vector input = {par().Sin, par().Sout}; + + return input; +} + +template +std::vector TBilinear::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +/* +/////Phase propagators////////////////////////// +template +LatticeSpinColourMatrix TBilinear::PhaseProps(LatticeSpinColourMatrix S, std::vector p) +{ + GridBase *grid = S._grid; + LatticeComplex pdotx(grid), coor(grid); + std::vector latt_size = grid->_fdimensions; + Complex Ci(0.0,1.0); + pdotx=zero; + for (unsigned int mu = 0; mu < 4; ++mu) + { + Real TwoPiL = M_PI * 2.0/ latt_size[mu]; + LatticeCoordinate(coor,mu); + pdotx = pdotx +(TwoPiL * p[mu]) * coor; + } + S = S*exp(-Ci*pdotx); + return S; +} +*/ +// execution /////////////////////////////////////////////////////////////////// +template +void TBilinear::execute(void) +{ +/************************************************************************** + +Compute the bilinear vertex needed for the NPR. +V(G) = sum_x [ g5 * adj(S'(x,p2)) * g5 * G * S'(x,p1) ]_{si,sj,ci,cj} +G is one of the 16 gamma vertices [I,gmu,g5,g5gmu,sig(mu,nu)] + + * G + / \ + p1/ \p2 + / \ + / \ + +Returns a spin-colour matrix, with indices si,sj, ci,cj + +Conventions: +p1 - incoming momenta +p2 - outgoing momenta +q = (p1-p2) +**************************************************************************/ + + LOG(Message) << "Computing bilinear contractions '" << getName() << "' using" + << " momentum '" << par().Sin << "' and '" << par().Sout << "'" + << std::endl; + + BinaryWriter writer(par().output); + + + // Propogators + LatticeSpinColourMatrix &Sin = *env().template getObject(par().Sin); + LatticeSpinColourMatrix &Sout = *env().template getObject(par().Sout); + LatticeComplex pdotxin(env().getGrid()), pdotxout(env().getGrid()), coor(env().getGrid()); + // momentum on legs + std::vector pin = strToVec(par().pin), pout = strToVec(par().pout); + std::vector latt_size(pin.begin(), pin.end()); + //bilinears + LatticeSpinColourMatrix bilinear_x(env().getGrid()); + SpinColourMatrix bilinear; + Gamma g5(Gamma::Algebra::Gamma5); + Result result; + Complex Ci(0.0,1.0); + + // + + pdotxin=zero; + pdotxout=zero; + for (unsigned int mu = 0; mu < 4; ++mu) + { + Real TwoPiL = M_PI * 2.0/ latt_size[mu]; + LatticeCoordinate(coor,mu); + pdotxin = pdotxin +(TwoPiL * pin[mu]) * coor; + pdotxout= pdotxout +(TwoPiL * pout[mu]) * coor; + } + Sin = Sin*exp(-Ci*pdotxin); //phase corrections + Sout = Sout*exp(-Ci*pdotxout); + + ////Set up gamma vector////////////////////////// + std::vector gammavector; + for( int i=0; i + +using namespace Grid; +using namespace Hadrons; +using namespace MNPR; + +template class Grid::Hadrons::MNPR::TFourQuark; + diff --git a/Hadrons/Modules/MNPR/FourQuark.hpp b/Hadrons/Modules/MNPR/FourQuark.hpp new file mode 100644 index 00000000..d71ddf43 --- /dev/null +++ b/Hadrons/Modules/MNPR/FourQuark.hpp @@ -0,0 +1,273 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MNPR/FourQuark.hpp + +Copyright (C) 2015 +Copyright (C) 2016 + +Author: Julia Kettle J.R.Kettle-2@sms.ed.ac.uk + +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_FourQuark_hpp_ +#define Hadrons_FourQuark_hpp_ + +#include +#include +#include +#include +#include +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * TFourQuark * + Performs fourquark contractions of the type tr[g5*adj(Sout)*g5*G*Sin] + Suitable for non exceptional momenta +******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MNPR) + +class FourQuarkPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(FourQuarkPar, + std::string, Sin, //need to make this a propogator type? + std::string, Sout, //same + std::string, pin, + std::string, pout, + bool, fullbasis, + std::string, output); +}; + +template +class TFourQuark: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl1, 1); + FERM_TYPE_ALIASES(FImpl2, 2); + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::vector, fourquark); + }; +public: + // constructor + TFourQuark(const std::string name); + // destructor + virtual ~TFourQuark(void) = default; + // dependencies/products + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void tensorprod(LatticeSpinColourSpinColourMatrix &lret, LatticeSpinColourMatrix a, LatticeSpinColourMatrix b); + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_TMP(FourQuark, ARG(TFourQuark), MNPR); + +/****************************************************************************** + * TFourQuark implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TFourQuark::TFourQuark(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TFourQuark::getInput(void) +{ + std::vector input = {par().Sin, par().Sout}; + + return input; +} + +template +std::vector TFourQuark::getOutput(void) +{ + std::vector output = {getName()}; + + return output; +} + + +template +void TFourQuark::tensorprod(LatticeSpinColourSpinColourMatrix &lret, LatticeSpinColourMatrix a, LatticeSpinColourMatrix b) +{ +#if 0 + parallel_for(auto site=lret.begin();site +void TFourQuark::setup(void) +{ + envCreateLat(LatticeSpinColourMatrix, getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TFourQuark::execute(void) +{ + +/********************************************************************************* + +TFourQuark : Creates the four quark vertex required for the NPR of four-quark ops + +V_{Gamma_1,Gamma_2} = sum_x [ ( g5 * adj(S'(x,p2)) * g5 * G1 * S'(x,p1) )_ci,cj;si,sj x ( g5 * adj(S'(x,p2)) * g5 * G2 S'(x,p1) )_ck,cl;sk,cl ] + +Create a bilinear vertex for G1 and G2 the spin and colour indices are kept free. Where there are 16 potential Gs. +We then find the outer product of V1 and V2, keeping the spin and colour indices uncontracted +Then this is summed over the lattice coordinate +Result is a SpinColourSpinColourMatrix - with 4 colour and 4 spin indices. +We have up to 256 of these including the offdiag (G1 != G2). + + \ / + \p1 p1/ + \ / + \ / + G1 * * G2 + / \ + / \ + /p2 p2\ + / \ + +*********************************************************************************/ + + + + + LOG(Message) << "Computing fourquark contractions '" << getName() << "' using" + << " momentum '" << par().Sin << "' and '" << par().Sout << "'" + << std::endl; + + BinaryWriter writer(par().output); + + PropagatorField1 &Sin = *env().template getObject(par().Sin); + PropagatorField2 &Sout = *env().template getObject(par().Sout); + std::vector pin = strToVec(par().pin), pout = strToVec(par().pout); + bool fullbasis = par().fullbasis; + Gamma g5(Gamma::Algebra::Gamma5); + Result result; + std::vector latt_size(pin.begin(), pin.end()); + LatticeComplex pdotxin(env().getGrid()), pdotxout(env().getGrid()), coor(env().getGrid()); + LatticeSpinColourMatrix bilinear_mu(env().getGrid()), bilinear_nu(env().getGrid()); + LatticeSpinColourSpinColourMatrix lret(env().getGrid()); + Complex Ci(0.0,1.0); + + //Phase propagators + //Sin = Grid::QCD::PropUtils::PhaseProps(Sin,pin); + //Sout = Grid::QCD::PropUtils::PhaseProps(Sout,pout); + + //find p.x for in and out so phase can be accounted for in propagators + pdotxin=zero; + pdotxout=zero; + for (unsigned int mu = 0; mu < 4; ++mu) + { + Real TwoPiL = M_PI * 2.0/ latt_size[mu]; + LatticeCoordinate(coor,mu); + pdotxin = pdotxin +(TwoPiL * pin[mu]) * coor; + pdotxout= pdotxout +(TwoPiL * pout[mu]) * coor; + } + Sin = Sin*exp(-Ci*pdotxin); //phase corrections + Sout = Sout*exp(-Ci*pdotxout); + + + //Set up Gammas + std::vector gammavector; + for( int i=1; i + +using namespace Grid; +using namespace Hadrons; +using namespace MSource; + +template class Grid::Hadrons::MSource::TMomSource; + diff --git a/Hadrons/Modules/MSource/Momentum.hpp b/Hadrons/Modules/MSource/Momentum.hpp new file mode 100644 index 00000000..7e4b3c09 --- /dev/null +++ b/Hadrons/Modules/MSource/Momentum.hpp @@ -0,0 +1,121 @@ +#ifndef Hadrons_MomSource_hpp_ +#define Hadrons_MomSource_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/* +Plane Wave source +----------------- +src_x = e^i2pi/L * p *position +*/ + +/****************************************************************************** + * TPlane * + ******************************************************************************/ + BEGIN_MODULE_NAMESPACE(MSource) + + class MomSourcePar: Serializable + { + public: +//What is meant by serializable in this context +GRID_SERIALIZABLE_CLASS_MEMBERS(MomSourcePar, +std::string, mom); +}; + + +template +class TMomSource: public Module +{ +public: +FERM_TYPE_ALIASES(FImpl,); +public: +// constructor +TMomSource(const std::string name); +// destructor +virtual ~TMomSource(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_TMP(MomSource, TMomSource, MSource); +//MODULE_REGISTER_NS(MomSource, TMomSource, MSource); + +/****************************************************************************** +* TMomSource template implementation * +******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TMomSource::TMomSource(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TMomSource::getInput(void) +{ + std::vector in; + return in; +} + +template +std::vector TMomSource::getOutput(void) +{ + std::vector out = {getName()}; + return out; +} + + +// setup /////////////////////////////////////////////////////////////////////// +template +void TMomSource::setup(void) +{ + envCreateLat(PropagatorField, getName()); +} + + +//execution////////////////////////////////////////////////////////////////// +template +void TMomSource::execute(void) +{ + LOG(Message) << "Generating planewave momentum source with momentum " << par().mom << std::endl; + //what does this env do? + PropagatorField &src = envGet(PropagatorField, getName()); + Lattice> t(env().getGrid()); + LatticeComplex C(env().getGrid()), coor(env().getGrid()); + std::vector p; + std::vector latt_size(GridDefaultLatt().begin(), GridDefaultLatt().end()); + Complex i(0.0,1.0); + + LOG(Message) << " " << std::endl; + //get the momentum from parameters + p = strToVec(par().mom); + C = zero; + LOG(Message) << "momentum converted from string - " << std::to_string(p[0]) < Date: Fri, 12 Oct 2018 11:00:58 +0100 Subject: [PATCH 06/37] NPR first cut by Julia Kettle --- Grid/qcd/QCD.h | 48 +++++++++++++++++++++++++++++++++------------ Hadrons/Modules.hpp | 5 +++++ Hadrons/modules.inc | 10 ++++++++++ 3 files changed, 51 insertions(+), 12 deletions(-) diff --git a/Grid/qcd/QCD.h b/Grid/qcd/QCD.h index e1e38e82..7e9b2da0 100644 --- a/Grid/qcd/QCD.h +++ b/Grid/qcd/QCD.h @@ -90,17 +90,20 @@ namespace QCD { // That probably makes for GridRedBlack4dCartesian grid. // s,sp,c,spc,lc - template using iSinglet = iScalar > >; - template using iSpinMatrix = iScalar, Ns> >; - template using iColourMatrix = iScalar > > ; - template using iSpinColourMatrix = iScalar, Ns> >; - template using iLorentzColourMatrix = iVector >, Nd > ; - template using iDoubleStoredColourMatrix = iVector >, Nds > ; - template using iSpinVector = iScalar, Ns> >; - template using iColourVector = iScalar > >; - template using iSpinColourVector = iScalar, Ns> >; - template using iHalfSpinVector = iScalar, Nhs> >; - template using iHalfSpinColourVector = iScalar, Nhs> >; + + template using iSinglet = iScalar > >; + template using iSpinMatrix = iScalar, Ns> >; + template using iColourMatrix = iScalar > > ; + template using iSpinColourMatrix = iScalar, Ns> >; + template using iLorentzColourMatrix = iVector >, Nd > ; + template using iDoubleStoredColourMatrix = iVector >, Nds > ; + template using iSpinVector = iScalar, Ns> >; + template using iColourVector = iScalar > >; + template using iSpinColourVector = iScalar, Ns> >; + template using iHalfSpinVector = iScalar, Nhs> >; + template using iHalfSpinColourVector = iScalar, Nhs> >; + template using iSpinColourSpinColourMatrix = iScalar, Ns>, Nc>, Ns> >; + template using iGparitySpinColourVector = iVector, Ns>, Ngp >; template using iGparityHalfSpinColourVector = iVector, Nhs>, Ngp >; @@ -127,10 +130,28 @@ namespace QCD { typedef iSpinColourMatrix SpinColourMatrix; typedef iSpinColourMatrix SpinColourMatrixF; typedef iSpinColourMatrix SpinColourMatrixD; - + typedef iSpinColourMatrix vSpinColourMatrix; typedef iSpinColourMatrix vSpinColourMatrixF; typedef iSpinColourMatrix vSpinColourMatrixD; + + // SpinColourSpinColour matrix + typedef iSpinColourSpinColourMatrix SpinColourSpinColourMatrix; + typedef iSpinColourSpinColourMatrix SpinColourSpinColourMatrixF; + typedef iSpinColourSpinColourMatrix SpinColourSpinColourMatrixD; + + typedef iSpinColourSpinColourMatrix vSpinColourSpinColourMatrix; + typedef iSpinColourSpinColourMatrix vSpinColourSpinColourMatrixF; + typedef iSpinColourSpinColourMatrix vSpinColourSpinColourMatrixD; + + // SpinColourSpinColour matrix + typedef iSpinColourSpinColourMatrix SpinColourSpinColourMatrix; + typedef iSpinColourSpinColourMatrix SpinColourSpinColourMatrixF; + typedef iSpinColourSpinColourMatrix SpinColourSpinColourMatrixD; + + typedef iSpinColourSpinColourMatrix vSpinColourSpinColourMatrix; + typedef iSpinColourSpinColourMatrix vSpinColourSpinColourMatrixF; + typedef iSpinColourSpinColourMatrix vSpinColourSpinColourMatrixD; // LorentzColour typedef iLorentzColourMatrix LorentzColourMatrix; @@ -229,6 +250,9 @@ namespace QCD { typedef Lattice LatticeSpinColourMatrixF; typedef Lattice LatticeSpinColourMatrixD; + typedef Lattice LatticeSpinColourSpinColourMatrix; + typedef Lattice LatticeSpinColourSpinColourMatrixF; + typedef Lattice LatticeSpinColourSpinColourMatrixD; typedef Lattice LatticeLorentzColourMatrix; typedef Lattice LatticeLorentzColourMatrixF; diff --git a/Hadrons/Modules.hpp b/Hadrons/Modules.hpp index 6d56ea2d..d86aed3c 100644 --- a/Hadrons/Modules.hpp +++ b/Hadrons/Modules.hpp @@ -16,6 +16,7 @@ #include #include #include +#include #include #include #include @@ -27,6 +28,7 @@ #include #include #include +#include #include #include #include @@ -40,6 +42,9 @@ #include #include #include +#include +#include +#include #include #include #include diff --git a/Hadrons/modules.inc b/Hadrons/modules.inc index 74138238..6eeb891d 100644 --- a/Hadrons/modules.inc +++ b/Hadrons/modules.inc @@ -11,6 +11,7 @@ modules_cc =\ Modules/MContraction/Gamma3pt.cc \ Modules/MFermion/FreeProp.cc \ Modules/MFermion/GaugeProp.cc \ + Modules/MSource/Momentum.cc \ Modules/MSource/Point.cc \ Modules/MSource/Wall.cc \ Modules/MSource/SeqConserved.cc \ @@ -28,6 +29,7 @@ modules_cc =\ Modules/MGauge/StochEm.cc \ Modules/MGauge/Random.cc \ Modules/MGauge/FundtoHirep.cc \ + Modules/MGauge/GaugeFix.cc \ Modules/MNoise/TimeDilutedSpinColorDiagonal.cc \ Modules/MUtilities/RandomVectors.cc \ Modules/MUtilities/TestSeqGamma.cc \ @@ -38,6 +40,9 @@ modules_cc =\ Modules/MScalar/VPCounterTerms.cc \ Modules/MScalar/ChargedProp.cc \ Modules/MScalar/ScalarVP.cc \ + Modules/MNPR/Amputate.cc \ + Modules/MNPR/Bilinear.cc \ + Modules/MNPR/FourQuark.cc \ Modules/MAction/Wilson.cc \ Modules/MAction/MobiusDWF.cc \ Modules/MAction/ZMobiusDWF.cc \ @@ -80,6 +85,7 @@ modules_hpp =\ Modules/MSource/Wall.hpp \ Modules/MSource/Z2.hpp \ Modules/MSource/SeqConserved.hpp \ + Modules/MSource/Momentum.hpp \ Modules/MSink/Smear.hpp \ Modules/MSink/Point.hpp \ Modules/MSolver/MixedPrecisionRBPrecCG.hpp \ @@ -91,6 +97,7 @@ modules_hpp =\ Modules/MGauge/StoutSmearing.hpp \ Modules/MGauge/Unit.hpp \ Modules/MGauge/Random.hpp \ + Modules/MGauge/GaugeFix.hpp \ Modules/MGauge/FundtoHirep.hpp \ Modules/MGauge/StochEm.hpp \ Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp \ @@ -104,6 +111,9 @@ modules_hpp =\ Modules/MScalar/ScalarVP.hpp \ Modules/MScalar/Scalar.hpp \ Modules/MScalar/ChargedProp.hpp \ + Modules/MNPR/Bilinear.hpp \ + Modules/MNPR/Amputate.hpp \ + Modules/MNPR/FourQuark.hpp \ Modules/MAction/DWF.hpp \ Modules/MAction/MobiusDWF.hpp \ Modules/MAction/Wilson.hpp \ From f0229025e2e513f86258e5a25724446b4e7d3e3c Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sat, 13 Oct 2018 19:55:03 +0100 Subject: [PATCH 07/37] MADWF working across a range of actions --- Grid/algorithms/iterative/ConjugateGradient.h | 14 +- Grid/qcd/action/fermion/CayleyFermion5D.cc | 6 +- Grid/qcd/action/fermion/CayleyFermion5D.h | 11 +- Grid/qcd/action/fermion/Fermion.h | 1 + .../qcd/action/fermion/FourierAcceleratedPV.h | 2 +- Grid/qcd/action/fermion/MADWF.h | 179 ++++++++++++++++++ tests/debug/Test_cayley_cg.cc | 35 +++- 7 files changed, 237 insertions(+), 11 deletions(-) create mode 100644 Grid/qcd/action/fermion/MADWF.h diff --git a/Grid/algorithms/iterative/ConjugateGradient.h b/Grid/algorithms/iterative/ConjugateGradient.h index ff4ba8ac..8a6a4bae 100644 --- a/Grid/algorithms/iterative/ConjugateGradient.h +++ b/Grid/algorithms/iterative/ConjugateGradient.h @@ -150,13 +150,13 @@ class ConjugateGradient : public OperatorFunction { std::cout << GridLogMessage << "\tTrue residual " << true_residual<::SetCoefficientsInternal(RealD zolo_hi,std::vector _gamma; + RealD _zolo_hi; + RealD _b; + RealD _c; + // Cayley form Moebius (tanh and zolotarev) std::vector omega; std::vector bs; // S dependent coeffs diff --git a/Grid/qcd/action/fermion/Fermion.h b/Grid/qcd/action/fermion/Fermion.h index d0200ae7..77a4681f 100644 --- a/Grid/qcd/action/fermion/Fermion.h +++ b/Grid/qcd/action/fermion/Fermion.h @@ -91,6 +91,7 @@ Author: Peter Boyle #include #include #include +#include //////////////////////////////////////////////////////////////////////////////////////////////////// // More maintainable to maintain the following typedef list centrally, as more "impl" targets diff --git a/Grid/qcd/action/fermion/FourierAcceleratedPV.h b/Grid/qcd/action/fermion/FourierAcceleratedPV.h index d8441cf2..d6196eee 100644 --- a/Grid/qcd/action/fermion/FourierAcceleratedPV.h +++ b/Grid/qcd/action/fermion/FourierAcceleratedPV.h @@ -7,7 +7,7 @@ Copyright (C) 2015 -Author: Christoph Lehner +Author: Christoph Lehner (lifted with permission by Peter Boyle, brought back to Grid) Author: Peter Boyle This program is free software; you can redistribute it and/or modify diff --git a/Grid/qcd/action/fermion/MADWF.h b/Grid/qcd/action/fermion/MADWF.h new file mode 100644 index 00000000..f7ffea2f --- /dev/null +++ b/Grid/qcd/action/fermion/MADWF.h @@ -0,0 +1,179 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/algorithms/iterative/MADWF.h + + Copyright (C) 2015 + +Author: Peter Boyle + + 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 */ +#pragma once + +namespace Grid { +namespace QCD { + +template +class MADWF +{ + private: + typedef typename Matrixo::FermionField FermionFieldo; + typedef typename Matrixi::FermionField FermionFieldi; + + PVinverter & PauliVillarsSolvero;// For the outer field + SchurSolver & SchurSolveri; // For the inner approx field + Guesser & Guesseri; // To deflate the inner approx solves + + Matrixo & Mato; // Action object for outer + Matrixi & Mati; // Action object for inner + + RealD target_resid; + int maxiter; + public: + + MADWF(Matrixo &_Mato, + Matrixi &_Mati, + PVinverter &_PauliVillarsSolvero, + SchurSolver &_SchurSolveri, + Guesser & _Guesseri, + RealD resid, + int _maxiter) : + + Mato(_Mato),Mati(_Mati), + SchurSolveri(_SchurSolveri), + PauliVillarsSolvero(_PauliVillarsSolvero),Guesseri(_Guesseri) + { + target_resid=resid; + maxiter =_maxiter; + }; + + void operator() (const FermionFieldo &src4,FermionFieldo &sol5) + { + std::cout << GridLogMessage<< " ************************************************" << std::endl; + std::cout << GridLogMessage<< " MADWF-like algorithm " << std::endl; + std::cout << GridLogMessage<< " ************************************************" << std::endl; + + FermionFieldo c0(Mato.GaugeGrid()); // 4d + FermionFieldo y0(Mato.GaugeGrid()); // 4d + + FermionFieldo A(Mato.FermionGrid()); // Temporary outer + FermionFieldo B(Mato.FermionGrid()); // Temporary outer + FermionFieldo b(Mato.FermionGrid()); // 5d source + + FermionFieldo c(Mato.FermionGrid()); // PVinv source; reused so store + FermionFieldo defect(Mato.FermionGrid()); // 5d source + + FermionFieldi ci(Mati.FermionGrid()); + FermionFieldi yi(Mati.FermionGrid()); + FermionFieldi xi(Mati.FermionGrid()); + FermionFieldi srci(Mati.FermionGrid()); + FermionFieldi Ai(Mati.FermionGrid()); + + RealD m=Mati.Mass(); + + /////////////////////////////////////// + //Import source, include Dminus factors + /////////////////////////////////////// + Mato.ImportPhysicalFermionSource(src4,b); + std::cout << GridLogMessage << " src4 " < HermOp(Ddwf); double Resid = 1.0e-12; + double Residi = 1.0e-6; ConjugateGradient CG(Resid,10000); + ConjugateGradient CGi(Residi,10000); Ddwf.ImportPhysicalFermionSource(src4,src); Ddwf.Mdag(src,src_NE); @@ -271,7 +274,8 @@ void TestReconstruct5D(What & Ddwf, // RBprec PV inverse //////////////////////////// typedef LatticeFermion Field; - typedef SchurRedBlackDiagMooeeSolve SchurSolverType; + typedef SchurRedBlackDiagTwoSolve SchurSolverType; + typedef SchurRedBlackDiagTwoSolve SchurSolverTypei; typedef PauliVillarsSolverRBprec PVinverter; SchurSolverType SchurSolver(CG); PVinverter PVinverse(SchurSolver); @@ -286,6 +290,19 @@ void TestReconstruct5D(What & Ddwf, result_rec = result_rec - result; std::cout < Guess; + MADWF > + madwf(Ddwf,Ddwf,PVinverse,SchurSolveri,Guess,Resid,10); + + madwf(src4,result_madwf); + result_madwf = result_madwf - result; + std::cout < void TestReconstruct5DFA(What & Ddwf, @@ -303,10 +320,13 @@ void TestReconstruct5DFA(What & Ddwf, LatticeFermion src_NE(FGrid); LatticeFermion result(FGrid); LatticeFermion result_rec(FGrid); + LatticeFermion result_madwf(FGrid); MdagMLinearOperator HermOp(Ddwf); double Resid = 1.0e-12; + double Residi = 1.0e-5; ConjugateGradient CG(Resid,10000); + ConjugateGradient CGi(Residi,10000); Ddwf.ImportPhysicalFermionSource(src4,src); Ddwf.Mdag(src,src_NE); @@ -324,6 +344,7 @@ void TestReconstruct5DFA(What & Ddwf, // Fourier accel PV inverse //////////////////////////// typedef LatticeFermion Field; + typedef SchurRedBlackDiagTwoSolve SchurSolverTypei; typedef PauliVillarsSolverFourierAccel PVinverter; PVinverter PVinverse(Umu,CG); @@ -337,6 +358,18 @@ void TestReconstruct5DFA(What & Ddwf, result_rec = result_rec - result; std::cout < Guess; + MADWF > + madwf(Ddwf,Ddwf,PVinverse,SchurSolver,Guess,Resid,10); + + madwf(src4,result_madwf); + result_madwf = result_madwf - result; + std::cout < Date: Sun, 14 Oct 2018 00:22:52 +0100 Subject: [PATCH 08/37] Mixed precision now supported in MADWF --- Grid/qcd/action/fermion/MADWF.h | 24 +++++++++++--- tests/debug/Test_cayley_cg.cc | 57 ++++++++++++++++++++++----------- 2 files changed, 57 insertions(+), 24 deletions(-) diff --git a/Grid/qcd/action/fermion/MADWF.h b/Grid/qcd/action/fermion/MADWF.h index f7ffea2f..064b13a8 100644 --- a/Grid/qcd/action/fermion/MADWF.h +++ b/Grid/qcd/action/fermion/MADWF.h @@ -30,6 +30,17 @@ Author: Peter Boyle namespace Grid { namespace QCD { +template X=0> +inline void convert(const Fieldi &from,Fieldo &to) +{ + precisionChange(to,from); +} +template X=0> +inline void convert(const Fieldi &from,Fieldo &to) +{ + to=from; +} + template class MADWF { @@ -70,8 +81,10 @@ class MADWF std::cout << GridLogMessage<< " MADWF-like algorithm " << std::endl; std::cout << GridLogMessage<< " ************************************************" << std::endl; - FermionFieldo c0(Mato.GaugeGrid()); // 4d - FermionFieldo y0(Mato.GaugeGrid()); // 4d + FermionFieldi c0i(Mati.GaugeGrid()); // 4d + FermionFieldi y0i(Mati.GaugeGrid()); // 4d + FermionFieldo c0 (Mato.GaugeGrid()); // 4d + FermionFieldo y0 (Mato.GaugeGrid()); // 4d FermionFieldo A(Mato.FermionGrid()); // Temporary outer FermionFieldo B(Mato.FermionGrid()); // Temporary outer @@ -110,15 +123,16 @@ class MADWF // Solve the inner system with surface term c0 //////////////////////////////////////////////// ci = zero; - InsertSlice(c0,ci,0, 0); + convert(c0,c0i); // Possible precison change + InsertSlice(c0i,ci,0, 0); // Dwm P y = Dwm x = D(1) P (c0,0,0,0)^T - Mati.P(ci,Ai); Mati.SetMass(1.0); Mati.M(Ai,srci); Mati.SetMass(m); SchurSolveri(Mati,srci,xi,Guesseri); Mati.Pdag(xi,yi); - ExtractSlice(y0, yi, 0 , 0); + ExtractSlice(y0i, yi, 0 , 0); + convert(y0i,y0); // Possible precision change ////////////////////////////////////// // Propagate solution back to outer system diff --git a/tests/debug/Test_cayley_cg.cc b/tests/debug/Test_cayley_cg.cc index 361b07e9..5150268f 100644 --- a/tests/debug/Test_cayley_cg.cc +++ b/tests/debug/Test_cayley_cg.cc @@ -85,8 +85,9 @@ void TestReconstruct5D(What & Ddwf, GridParallelRNG *RNG4, GridParallelRNG *RNG5); -template +template void TestReconstruct5DFA(What & Ddwf, + WhatF & DdwfF, LatticeGaugeField &Umu, GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, @@ -102,19 +103,31 @@ int main (int argc, char ** argv) std::cout< seeds4({1,2,3,4}); std::vector seeds5({5,6,7,8}); GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); LatticeGaugeField Umu(UGrid); + LatticeGaugeFieldF UmuF(UGridF); SU3::HotConfiguration(RNG4,Umu); + precisionChange(UmuF,Umu); std::vector U(4,UGrid); RealD mass=0.1; @@ -123,8 +136,9 @@ int main (int argc, char ** argv) std::cout<(Ddwf,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); - TestReconstruct5DFA(Ddwf,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5DFA(Ddwf,DdwfF,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); RealD b=1.5;// Scale factor b+c=2, b-c=1 RealD c=0.5; @@ -134,8 +148,9 @@ int main (int argc, char ** argv) std::cout<(Dmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); - TestReconstruct5DFA(Dmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5DFA(Dmob,DmobF,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dsham,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); - TestReconstruct5DFA(Dsham,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5DFA(Dsham,DshamF,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dov,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); - TestReconstruct5DFA(Dov,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5DFA(Dov,DovF,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout< +template void TestReconstruct5DFA(What & Ddwf, - LatticeGaugeField & Umu, - GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, - GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, - RealD mass, RealD M5, - GridParallelRNG *RNG4, - GridParallelRNG *RNG5) + WhatF & DdwfF, + LatticeGaugeField & Umu, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5) { LatticeFermion src4 (UGrid); random(*RNG4,src4); LatticeFermion res4 (UGrid); res4 = zero; @@ -326,7 +344,7 @@ void TestReconstruct5DFA(What & Ddwf, double Resid = 1.0e-12; double Residi = 1.0e-5; ConjugateGradient CG(Resid,10000); - ConjugateGradient CGi(Residi,10000); + ConjugateGradient CGi(Residi,10000); Ddwf.ImportPhysicalFermionSource(src4,src); Ddwf.Mdag(src,src_NE); @@ -344,7 +362,8 @@ void TestReconstruct5DFA(What & Ddwf, // Fourier accel PV inverse //////////////////////////// typedef LatticeFermion Field; - typedef SchurRedBlackDiagTwoSolve SchurSolverTypei; + typedef LatticeFermionF FieldF; + typedef SchurRedBlackDiagTwoSolve SchurSolverTypei; typedef PauliVillarsSolverFourierAccel PVinverter; PVinverter PVinverse(Umu,CG); @@ -362,9 +381,9 @@ void TestReconstruct5DFA(What & Ddwf, // Now try MADWF ////////////////////////////// SchurSolverTypei SchurSolver(CGi); - ZeroGuesser Guess; - MADWF > - madwf(Ddwf,Ddwf,PVinverse,SchurSolver,Guess,Resid,10); + ZeroGuesser Guess; + MADWF > + madwf(Ddwf,DdwfF,PVinverse,SchurSolver,Guess,Resid,10); madwf(src4,result_madwf); result_madwf = result_madwf - result; From 959211534198a4ebecba1aafb7f3cad8eae05e7e Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Mon, 15 Oct 2018 15:49:42 +0100 Subject: [PATCH 09/37] Hadrons: NPR and gauge fixing linking fix --- Hadrons/Modules/MGauge/GaugeFix.cc | 57 ---------------------------- Hadrons/Modules/MGauge/GaugeFix.hpp | 56 ++++++++++++++++++++++++++- Hadrons/Modules/MNPR/Amputate.hpp | 2 +- Hadrons/Modules/MNPR/Bilinear.hpp | 2 +- Hadrons/Modules/MNPR/FourQuark.hpp | 2 +- Hadrons/Modules/MSource/Momentum.cc | 2 +- Hadrons/Modules/MSource/Momentum.hpp | 44 ++++++++++----------- 7 files changed, 81 insertions(+), 84 deletions(-) diff --git a/Hadrons/Modules/MGauge/GaugeFix.cc b/Hadrons/Modules/MGauge/GaugeFix.cc index bd57ed07..1cae39ad 100644 --- a/Hadrons/Modules/MGauge/GaugeFix.cc +++ b/Hadrons/Modules/MGauge/GaugeFix.cc @@ -33,61 +33,4 @@ using namespace Grid; using namespace Hadrons; using namespace MGauge; - -/****************************************************************************** -* TGaugeFix implementation * -******************************************************************************/ -// constructor ///////////////////////////////////////////////////////////////// -template -TGaugeFix::TGaugeFix(const std::string name) -: Module(name) -{} - -// dependencies/products /////////////////////////////////////////////////////// -template -std::vector TGaugeFix::getInput(void) -{ - std::vector in = {par().gauge}; - return in; -} - -template -std::vector TGaugeFix::getOutput(void) -{ - std::vector out = {getName()}; - return out; -} - -// setup /////////////////////////////////////////////////////////////////////// -template -void TGaugeFix::setup(void) -{ - envCreateLat(GaugeField, getName()); -} - - -// execution /////////////////////////////////////////////////////////////////// -template -void TGaugeFix::execute(void) -//Loads the gauge and fixes it -{ - - std::cout << "executing" << std::endl; - LOG(Message) << "Fixing the Gauge" << std::endl; - LOG(Message) << par().gauge << std::endl; - auto &U = envGet(GaugeField, par().gauge); - auto &Umu = envGet(GaugeField, getName()); - LOG(Message) << "Gauge Field fetched" << std::endl; - //do we allow maxiter etc to be user set? - Real alpha = par().alpha; - int maxiter = par().maxiter; - Real Omega_tol = par().Omega_tol; - Real Phi_tol = par().Phi_tol; - bool Fourier = par().Fourier; - FourierAcceleratedGaugeFixer::SteepestDescentGaugeFix(U,alpha,maxiter,Omega_tol,Phi_tol,Fourier); - Umu = U; - LOG(Message) << "Gauge Fixed" << std::endl; - -} - template class Grid::Hadrons::MGauge::TGaugeFix; diff --git a/Hadrons/Modules/MGauge/GaugeFix.hpp b/Hadrons/Modules/MGauge/GaugeFix.hpp index 154237a7..e8625b01 100644 --- a/Hadrons/Modules/MGauge/GaugeFix.hpp +++ b/Hadrons/Modules/MGauge/GaugeFix.hpp @@ -62,7 +62,7 @@ public: // constructor TGaugeFix(const std::string name); // destructor - virtual ~TGaugeFix(void) = default; + virtual ~TGaugeFix(void) {}; // dependencies/products virtual std::vector getInput(void); virtual std::vector getOutput(void); @@ -74,6 +74,60 @@ public: MODULE_REGISTER_TMP(GaugeFix, TGaugeFix, MGauge); +/****************************************************************************** +* TGaugeFix implementation * +******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TGaugeFix::TGaugeFix(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TGaugeFix::getInput(void) +{ + std::vector in = {par().gauge}; + return in; +} + +template +std::vector TGaugeFix::getOutput(void) +{ + std::vector out = {getName()}; + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TGaugeFix::setup(void) +{ + envCreateLat(GaugeField, getName()); +} + + +// execution /////////////////////////////////////////////////////////////////// +template +void TGaugeFix::execute(void) +//Loads the gauge and fixes it +{ + std::cout << "executing" << std::endl; + LOG(Message) << "Fixing the Gauge" << std::endl; + LOG(Message) << par().gauge << std::endl; + auto &U = envGet(GaugeField, par().gauge); + auto &Umu = envGet(GaugeField, getName()); + LOG(Message) << "Gauge Field fetched" << std::endl; + //do we allow maxiter etc to be user set? + Real alpha = par().alpha; + int maxiter = par().maxiter; + Real Omega_tol = par().Omega_tol; + Real Phi_tol = par().Phi_tol; + bool Fourier = par().Fourier; + FourierAcceleratedGaugeFixer::SteepestDescentGaugeFix(U,alpha,maxiter,Omega_tol,Phi_tol,Fourier); + Umu = U; + LOG(Message) << "Gauge Fixed" << std::endl; +} + END_MODULE_NAMESPACE END_HADRONS_NAMESPACE diff --git a/Hadrons/Modules/MNPR/Amputate.hpp b/Hadrons/Modules/MNPR/Amputate.hpp index 797247b5..417b81e4 100644 --- a/Hadrons/Modules/MNPR/Amputate.hpp +++ b/Hadrons/Modules/MNPR/Amputate.hpp @@ -75,7 +75,7 @@ public: // constructor TAmputate(const std::string name); // destructor - virtual ~TAmputate(void) = default; + virtual ~TAmputate(void) {}; // dependencies/products virtual std::vector getInput(void); virtual std::vector getOutput(void); diff --git a/Hadrons/Modules/MNPR/Bilinear.hpp b/Hadrons/Modules/MNPR/Bilinear.hpp index 0dc621b9..85d4fda7 100644 --- a/Hadrons/Modules/MNPR/Bilinear.hpp +++ b/Hadrons/Modules/MNPR/Bilinear.hpp @@ -72,7 +72,7 @@ public: // constructor TBilinear(const std::string name); // destructor - virtual ~TBilinear(void) = default; + virtual ~TBilinear(void) {}; // dependencies/products virtual std::vector getInput(void); virtual std::vector getOutput(void); diff --git a/Hadrons/Modules/MNPR/FourQuark.hpp b/Hadrons/Modules/MNPR/FourQuark.hpp index d71ddf43..35ef7cad 100644 --- a/Hadrons/Modules/MNPR/FourQuark.hpp +++ b/Hadrons/Modules/MNPR/FourQuark.hpp @@ -72,7 +72,7 @@ public: // constructor TFourQuark(const std::string name); // destructor - virtual ~TFourQuark(void) = default; + virtual ~TFourQuark(void) {}; // dependencies/products virtual std::vector getInput(void); virtual std::vector getOutput(void); diff --git a/Hadrons/Modules/MSource/Momentum.cc b/Hadrons/Modules/MSource/Momentum.cc index 66280600..2a55f5dc 100644 --- a/Hadrons/Modules/MSource/Momentum.cc +++ b/Hadrons/Modules/MSource/Momentum.cc @@ -4,5 +4,5 @@ using namespace Grid; using namespace Hadrons; using namespace MSource; -template class Grid::Hadrons::MSource::TMomSource; +template class Grid::Hadrons::MSource::TMomentum; diff --git a/Hadrons/Modules/MSource/Momentum.hpp b/Hadrons/Modules/MSource/Momentum.hpp index 7e4b3c09..b6e93bbc 100644 --- a/Hadrons/Modules/MSource/Momentum.hpp +++ b/Hadrons/Modules/MSource/Momentum.hpp @@ -1,5 +1,5 @@ -#ifndef Hadrons_MomSource_hpp_ -#define Hadrons_MomSource_hpp_ +#ifndef Hadrons_Momentum_hpp_ +#define Hadrons_Momentum_hpp_ #include #include @@ -14,29 +14,29 @@ src_x = e^i2pi/L * p *position */ /****************************************************************************** - * TPlane * - ******************************************************************************/ - BEGIN_MODULE_NAMESPACE(MSource) + * Plane Wave source * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MSource) - class MomSourcePar: Serializable - { - public: +class MomentumPar: Serializable +{ +public: //What is meant by serializable in this context -GRID_SERIALIZABLE_CLASS_MEMBERS(MomSourcePar, +GRID_SERIALIZABLE_CLASS_MEMBERS(MomentumPar, std::string, mom); }; template -class TMomSource: public Module +class TMomentum: public Module { public: FERM_TYPE_ALIASES(FImpl,); public: // constructor -TMomSource(const std::string name); +TMomentum(const std::string name); // destructor -virtual ~TMomSource(void) = default; +virtual ~TMomentum(void) {}; // dependency relation virtual std::vector getInput(void); virtual std::vector getOutput(void); @@ -46,28 +46,28 @@ virtual void setup(void); virtual void execute(void); }; -MODULE_REGISTER_TMP(MomSource, TMomSource, MSource); -//MODULE_REGISTER_NS(MomSource, TMomSource, MSource); +MODULE_REGISTER_TMP(Momentum, TMomentum, MSource); +//MODULE_REGISTER_NS(Momentum, TMomentum, MSource); /****************************************************************************** -* TMomSource template implementation * +* TMomentum template implementation * ******************************************************************************/ // constructor ///////////////////////////////////////////////////////////////// template -TMomSource::TMomSource(const std::string name) -: Module(name) +TMomentum::TMomentum(const std::string name) +: Module(name) {} // dependencies/products /////////////////////////////////////////////////////// template -std::vector TMomSource::getInput(void) +std::vector TMomentum::getInput(void) { std::vector in; return in; } template -std::vector TMomSource::getOutput(void) +std::vector TMomentum::getOutput(void) { std::vector out = {getName()}; return out; @@ -76,7 +76,7 @@ std::vector TMomSource::getOutput(void) // setup /////////////////////////////////////////////////////////////////////// template -void TMomSource::setup(void) +void TMomentum::setup(void) { envCreateLat(PropagatorField, getName()); } @@ -84,7 +84,7 @@ void TMomSource::setup(void) //execution////////////////////////////////////////////////////////////////// template -void TMomSource::execute(void) +void TMomentum::execute(void) { LOG(Message) << "Generating planewave momentum source with momentum " << par().mom << std::endl; //what does this env do? @@ -118,4 +118,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_MomSource_hpp_ +#endif // Hadrons_Momentum_hpp_ From 2f368c33fc7db569d202f6e53edc9595b3ab03b6 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Mon, 15 Oct 2018 15:51:45 +0100 Subject: [PATCH 10/37] Hadrons: copyright update --- .../Modules/MContraction/A2AAslashField.cc | 27 ++++++++++++++++++ .../Modules/MContraction/A2AAslashField.hpp | 27 ++++++++++++++++++ Hadrons/Modules/MGauge/GaugeFix.cc | 6 ++-- Hadrons/Modules/MGauge/GaugeFix.hpp | 6 ++-- Hadrons/Modules/MNPR/Amputate.cc | 27 ++++++++++++++++++ Hadrons/Modules/MNPR/Amputate.hpp | 7 +++-- Hadrons/Modules/MNPR/Bilinear.cc | 27 ++++++++++++++++++ Hadrons/Modules/MNPR/Bilinear.hpp | 7 +++-- Hadrons/Modules/MNPR/FourQuark.cc | 27 ++++++++++++++++++ Hadrons/Modules/MNPR/FourQuark.hpp | 7 +++-- Hadrons/Modules/MSource/Momentum.cc | 28 +++++++++++++++++++ Hadrons/Modules/MSource/Momentum.hpp | 28 +++++++++++++++++++ Hadrons/TimerArray.cc | 27 ++++++++++++++++++ Hadrons/TimerArray.hpp | 27 ++++++++++++++++++ Hadrons/Utilities/EigenPackCast.cc | 27 ++++++++++++++++++ Hadrons/Utilities/HadronsXmlRun.cc | 2 +- 16 files changed, 291 insertions(+), 16 deletions(-) diff --git a/Hadrons/Modules/MContraction/A2AAslashField.cc b/Hadrons/Modules/MContraction/A2AAslashField.cc index eb76932e..65c49198 100644 --- a/Hadrons/Modules/MContraction/A2AAslashField.cc +++ b/Hadrons/Modules/MContraction/A2AAslashField.cc @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MContraction/A2AAslashField.cc + +Copyright (C) 2015-2018 + +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; diff --git a/Hadrons/Modules/MContraction/A2AAslashField.hpp b/Hadrons/Modules/MContraction/A2AAslashField.hpp index 792ff6ee..cc8b1559 100644 --- a/Hadrons/Modules/MContraction/A2AAslashField.hpp +++ b/Hadrons/Modules/MContraction/A2AAslashField.hpp @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MContraction/A2AAslashField.hpp + +Copyright (C) 2015-2018 + +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_MContraction_A2AAslashField_hpp_ #define Hadrons_MContraction_A2AAslashField_hpp_ diff --git a/Hadrons/Modules/MGauge/GaugeFix.cc b/Hadrons/Modules/MGauge/GaugeFix.cc index 1cae39ad..53aa16da 100644 --- a/Hadrons/Modules/MGauge/GaugeFix.cc +++ b/Hadrons/Modules/MGauge/GaugeFix.cc @@ -2,12 +2,12 @@ Grid physics library, www.github.com/paboyle/Grid -Source file: extras/Hadrons/Modules/MGauge/Fix.cc +Source file: Hadrons/Modules/MGauge/GaugeFix.cc -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli +Author: Peter Boyle 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 diff --git a/Hadrons/Modules/MGauge/GaugeFix.hpp b/Hadrons/Modules/MGauge/GaugeFix.hpp index e8625b01..ece8c19d 100644 --- a/Hadrons/Modules/MGauge/GaugeFix.hpp +++ b/Hadrons/Modules/MGauge/GaugeFix.hpp @@ -2,12 +2,12 @@ Grid physics library, www.github.com/paboyle/Grid -Source file: extras/Hadrons/Modules/MGauge/Fix.hpp +Source file: Hadrons/Modules/MGauge/GaugeFix.hpp -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli +Author: Peter Boyle 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 diff --git a/Hadrons/Modules/MNPR/Amputate.cc b/Hadrons/Modules/MNPR/Amputate.cc index c4ada9d9..8d55accb 100644 --- a/Hadrons/Modules/MNPR/Amputate.cc +++ b/Hadrons/Modules/MNPR/Amputate.cc @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MNPR/Amputate.cc + +Copyright (C) 2015-2018 + +Author: Peter Boyle + +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; diff --git a/Hadrons/Modules/MNPR/Amputate.hpp b/Hadrons/Modules/MNPR/Amputate.hpp index 417b81e4..953bb354 100644 --- a/Hadrons/Modules/MNPR/Amputate.hpp +++ b/Hadrons/Modules/MNPR/Amputate.hpp @@ -2,12 +2,13 @@ Grid physics library, www.github.com/paboyle/Grid -Source file: extras/Hadrons/Modules/MNpr/Amputate.hpp +Source file: Hadrons/Modules/MNPR/Amputate.hpp -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 +Author: Antonin Portelli Author: Julia Kettle J.R.Kettle-2@sms.ed.ac.uk +Author: Peter Boyle 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 diff --git a/Hadrons/Modules/MNPR/Bilinear.cc b/Hadrons/Modules/MNPR/Bilinear.cc index 01fb08c6..865739dd 100644 --- a/Hadrons/Modules/MNPR/Bilinear.cc +++ b/Hadrons/Modules/MNPR/Bilinear.cc @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MNPR/Bilinear.cc + +Copyright (C) 2015-2018 + +Author: Peter Boyle + +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; diff --git a/Hadrons/Modules/MNPR/Bilinear.hpp b/Hadrons/Modules/MNPR/Bilinear.hpp index 85d4fda7..66cd22a6 100644 --- a/Hadrons/Modules/MNPR/Bilinear.hpp +++ b/Hadrons/Modules/MNPR/Bilinear.hpp @@ -2,12 +2,13 @@ Grid physics library, www.github.com/paboyle/Grid -Source file: extras/Hadrons/Modules/MContraction/Bilinear.hpp +Source file: Hadrons/Modules/MNPR/Bilinear.hpp -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 +Author: Antonin Portelli Author: Julia Kettle J.R.Kettle-2@sms.ed.ac.uk +Author: Peter Boyle 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 diff --git a/Hadrons/Modules/MNPR/FourQuark.cc b/Hadrons/Modules/MNPR/FourQuark.cc index 6e951434..9f14ff0f 100644 --- a/Hadrons/Modules/MNPR/FourQuark.cc +++ b/Hadrons/Modules/MNPR/FourQuark.cc @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MNPR/FourQuark.cc + +Copyright (C) 2015-2018 + +Author: Peter Boyle + +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; diff --git a/Hadrons/Modules/MNPR/FourQuark.hpp b/Hadrons/Modules/MNPR/FourQuark.hpp index 35ef7cad..f961a366 100644 --- a/Hadrons/Modules/MNPR/FourQuark.hpp +++ b/Hadrons/Modules/MNPR/FourQuark.hpp @@ -2,12 +2,13 @@ Grid physics library, www.github.com/paboyle/Grid -Source file: extras/Hadrons/Modules/MNPR/FourQuark.hpp +Source file: Hadrons/Modules/MNPR/FourQuark.hpp -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 +Author: Antonin Portelli Author: Julia Kettle J.R.Kettle-2@sms.ed.ac.uk +Author: Peter Boyle 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 diff --git a/Hadrons/Modules/MSource/Momentum.cc b/Hadrons/Modules/MSource/Momentum.cc index 2a55f5dc..9bcf65ae 100644 --- a/Hadrons/Modules/MSource/Momentum.cc +++ b/Hadrons/Modules/MSource/Momentum.cc @@ -1,3 +1,31 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MSource/Momentum.cc + +Copyright (C) 2015-2018 + +Author: Antonin Portelli +Author: Peter Boyle + +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; diff --git a/Hadrons/Modules/MSource/Momentum.hpp b/Hadrons/Modules/MSource/Momentum.hpp index b6e93bbc..8336ea5c 100644 --- a/Hadrons/Modules/MSource/Momentum.hpp +++ b/Hadrons/Modules/MSource/Momentum.hpp @@ -1,3 +1,31 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MSource/Momentum.hpp + +Copyright (C) 2015-2018 + +Author: Antonin Portelli +Author: Peter Boyle + +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_Momentum_hpp_ #define Hadrons_Momentum_hpp_ diff --git a/Hadrons/TimerArray.cc b/Hadrons/TimerArray.cc index c99ed41f..2b85bc7e 100644 --- a/Hadrons/TimerArray.cc +++ b/Hadrons/TimerArray.cc @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/TimerArray.cc + +Copyright (C) 2015-2018 + +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; diff --git a/Hadrons/TimerArray.hpp b/Hadrons/TimerArray.hpp index 477f11cd..77cc2b8c 100644 --- a/Hadrons/TimerArray.hpp +++ b/Hadrons/TimerArray.hpp @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/TimerArray.hpp + +Copyright (C) 2015-2018 + +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_TimerArray_hpp_ #define Hadrons_TimerArray_hpp_ diff --git a/Hadrons/Utilities/EigenPackCast.cc b/Hadrons/Utilities/EigenPackCast.cc index 7247c92d..8f2755e3 100644 --- a/Hadrons/Utilities/EigenPackCast.cc +++ b/Hadrons/Utilities/EigenPackCast.cc @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Utilities/EigenPackCast.cc + +Copyright (C) 2015-2018 + +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 diff --git a/Hadrons/Utilities/HadronsXmlRun.cc b/Hadrons/Utilities/HadronsXmlRun.cc index c6ea7dc5..9b8fc901 100644 --- a/Hadrons/Utilities/HadronsXmlRun.cc +++ b/Hadrons/Utilities/HadronsXmlRun.cc @@ -2,7 +2,7 @@ Grid physics library, www.github.com/paboyle/Grid -Source file: Hadrons/HadronsXmlRun.cc +Source file: Hadrons/Utilities/HadronsXmlRun.cc Copyright (C) 2015-2018 From 291bc2a1f086925b3778aac73dc0f47f0b0543fb Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Mon, 15 Oct 2018 17:25:08 +0100 Subject: [PATCH 11/37] IO benchmark on a list of directories --- benchmarks/Benchmark_IO.cc | 102 ++++++---------------------- benchmarks/Benchmark_IO.hpp | 107 ++++++++++++++++++++++++++++++ benchmarks/Benchmark_IO_vs_dir.cc | 63 ++++++++++++++++++ 3 files changed, 191 insertions(+), 81 deletions(-) create mode 100644 benchmarks/Benchmark_IO.hpp create mode 100644 benchmarks/Benchmark_IO_vs_dir.cc diff --git a/benchmarks/Benchmark_IO.cc b/benchmarks/Benchmark_IO.cc index 479ae037..58e0340d 100644 --- a/benchmarks/Benchmark_IO.cc +++ b/benchmarks/Benchmark_IO.cc @@ -1,108 +1,48 @@ -#include -#ifdef HAVE_LIME -using namespace std; -using namespace Grid; -using namespace Grid::QCD; +#include "Benchmark_IO.hpp" -#define MSG cout << GridLogMessage -#define SEP \ -"=============================================================================" #ifndef BENCH_IO_LMAX #define BENCH_IO_LMAX 40 #endif -typedef function WriterFn; -typedef function ReaderFn; +using namespace Grid; +using namespace QCD; -string filestem(const int l) +std::string filestem(const int l) { - return "iobench_l" + to_string(l); -} - -void limeWrite(const string filestem, LatticeFermion &vec) -{ - emptyUserRecord record; - ScidacWriter binWriter(vec._grid->IsBoss()); - - binWriter.open(filestem + ".bin"); - binWriter.writeScidacFieldRecord(vec, record); - binWriter.close(); -} - -void limeRead(LatticeFermion &vec, const string filestem) -{ - emptyUserRecord record; - ScidacReader binReader; - - binReader.open(filestem + ".bin"); - binReader.readScidacFieldRecord(vec, record); - binReader.close(); -} - -void writeBenchmark(const int l, const WriterFn &write) -{ - auto mpi = GridDefaultMpi(); - auto simd = GridDefaultSimd(Nd, vComplex::Nsimd()); - vector latt = {l*mpi[0], l*mpi[1], l*mpi[2], l*mpi[3]}; - unique_ptr gPt(SpaceTimeGrid::makeFourDimGrid(latt, simd, mpi)); - GridCartesian *g = gPt.get(); - GridParallelRNG rng(g); - LatticeFermion vec(g); - emptyUserRecord record; - ScidacWriter binWriter(g->IsBoss()); - - cout << "-- Local volume " << l << "^4" << endl; - random(rng, vec); - write(filestem(l), vec); -} - -void readBenchmark(const int l, const ReaderFn &read) -{ - auto mpi = GridDefaultMpi(); - auto simd = GridDefaultSimd(Nd, vComplex::Nsimd()); - vector latt = {l*mpi[0], l*mpi[1], l*mpi[2], l*mpi[3]}; - unique_ptr gPt(SpaceTimeGrid::makeFourDimGrid(latt, simd, mpi)); - GridCartesian *g = gPt.get(); - LatticeFermion vec(g); - emptyUserRecord record; - ScidacReader binReader; - - cout << "-- Local volume " << l << "^4" << endl; - read(vec, filestem(l)); + return "iobench_l" + std::to_string(l); } int main (int argc, char ** argv) { Grid_init(&argc,&argv); - auto simd = GridDefaultSimd(Nd,vComplex::Nsimd()); - auto mpi = GridDefaultMpi(); - int64_t threads = GridThread::GetThreads(); - MSG << "Grid is setup to use " << threads << " threads" << endl; - MSG << SEP << endl; - MSG << "Benchmark Lime write" << endl; - MSG << SEP << endl; + MSG << "Grid is setup to use " << threads << " threads" << std::endl; + MSG << SEP << std::endl; + MSG << "Benchmark Lime write" << std::endl; + MSG << SEP << std::endl; for (int l = 4; l <= BENCH_IO_LMAX; l += 2) { - writeBenchmark(l, limeWrite); + auto mpi = GridDefaultMpi(); + std::vector latt = {l*mpi[0], l*mpi[1], l*mpi[2], l*mpi[3]}; + + std::cout << "-- Local volume " << l << "^4" << std::endl; + writeBenchmark(latt, filestem(l), limeWrite); } - MSG << "Benchmark Lime read" << endl; - MSG << SEP << endl; + MSG << "Benchmark Lime read" << std::endl; + MSG << SEP << std::endl; for (int l = 4; l <= BENCH_IO_LMAX; l += 2) { - readBenchmark(l, limeRead); + auto mpi = GridDefaultMpi(); + std::vector latt = {l*mpi[0], l*mpi[1], l*mpi[2], l*mpi[3]}; + + std::cout << "-- Local volume " << l << "^4" << std::endl; + readBenchmark(latt, filestem(l), limeRead); } Grid_finalize(); return EXIT_SUCCESS; } -#else -int main (int argc, char ** argv) -{ - return EXIT_SUCCESS; -} -#endif diff --git a/benchmarks/Benchmark_IO.hpp b/benchmarks/Benchmark_IO.hpp new file mode 100644 index 00000000..9565c928 --- /dev/null +++ b/benchmarks/Benchmark_IO.hpp @@ -0,0 +1,107 @@ +#ifndef Benchmark_IO_hpp_ +#define Benchmark_IO_hpp_ + +#include + +#define MSG std::cout << GridLogMessage +#define SEP \ +"=============================================================================" + +namespace Grid { + +template +using WriterFn = std::function ; +template +using ReaderFn = std::function; + +template +void limeWrite(const std::string filestem, Field &vec) +{ + emptyUserRecord record; + QCD::ScidacWriter binWriter(vec._grid->IsBoss()); + + binWriter.open(filestem + ".bin"); + binWriter.writeScidacFieldRecord(vec, record); + binWriter.close(); +} + +template +void limeRead(Field &vec, const std::string filestem) +{ + emptyUserRecord record; + QCD::ScidacReader binReader; + + binReader.open(filestem + ".bin"); + binReader.readScidacFieldRecord(vec, record); + binReader.close(); +} + +inline void makeGrid(std::shared_ptr &gPt, + const std::shared_ptr &gBasePt, + const unsigned int Ls = 1, const bool rb = false) +{ + if (rb) + { + if (Ls > 1) + { + gPt.reset(QCD::SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, gBasePt.get())); + } + else + { + gPt.reset(QCD::SpaceTimeGrid::makeFourDimRedBlackGrid(gBasePt.get())); + } + } + else + { + if (Ls > 1) + { + gPt.reset(QCD::SpaceTimeGrid::makeFiveDimGrid(Ls, gBasePt.get())); + } + else + { + gPt = gBasePt; + } + } +} + +template +void writeBenchmark(const std::vector &latt, const std::string filename, + const WriterFn &write, + const unsigned int Ls = 1, const bool rb = false) +{ + auto mpi = GridDefaultMpi(); + auto simd = GridDefaultSimd(latt.size(), Field::vector_type::Nsimd()); + std::shared_ptr gBasePt(QCD::SpaceTimeGrid::makeFourDimGrid(latt, simd, mpi)); + std::shared_ptr gPt; + + makeGrid(gPt, gBasePt, Ls, rb); + + GridBase *g = gPt.get(); + GridParallelRNG rng(g); + Field vec(g); + + random(rng, vec); + write(filename, vec); +} + +template +void readBenchmark(const std::vector &latt, const std::string filename, + const ReaderFn &read, + const unsigned int Ls = 1, const bool rb = false) +{ + auto mpi = GridDefaultMpi(); + auto simd = GridDefaultSimd(latt.size(), Field::vector_type::Nsimd()); + std::shared_ptr gBasePt(QCD::SpaceTimeGrid::makeFourDimGrid(latt, simd, mpi)); + std::shared_ptr gPt; + + makeGrid(gPt, gBasePt, Ls, rb); + + GridBase *g = gPt.get(); + Field vec(g); + + read(vec, filename); +} + +} + +#endif // Benchmark_IO_hpp_ diff --git a/benchmarks/Benchmark_IO_vs_dir.cc b/benchmarks/Benchmark_IO_vs_dir.cc new file mode 100644 index 00000000..c22550c3 --- /dev/null +++ b/benchmarks/Benchmark_IO_vs_dir.cc @@ -0,0 +1,63 @@ +#include "Benchmark_IO.hpp" + +#define MSG std::cout << GridLogMessage +#define SEP \ +"=============================================================================" +#ifndef BENCH_IO_LMAX +#define BENCH_IO_LMAX 40 +#endif + +using namespace Grid; +using namespace QCD; + +int main (int argc, char ** argv) +{ + std::vector dir; + unsigned int Ls; + bool rb; + if (argc < 4) + { + std::cerr << "usage: " << argv[0] << " [ ... ] [Grid options]"; + std::cerr << std::endl; + } + Ls = std::stoi(argv[1]); + rb = (std::string(argv[2]) == "1"); + for (unsigned int i = 3; i < argc; ++i) + { + std::string a = argv[i]; + + if (a[0] != '-') + { + dir.push_back(std::string(argv[i])); + } + else + { + break; + } + } + Grid_init(&argc,&argv); + + + int64_t threads = GridThread::GetThreads(); + MSG << "Grid is setup to use " << threads << " threads" << std::endl; + MSG << SEP << std::endl; + MSG << "Benchmark Lime write" << std::endl; + MSG << SEP << std::endl; + for (auto &d: dir) + { + MSG << "-- Directory " << d << std::endl; + writeBenchmark(GridDefaultLatt(), d + "/ioBench", limeWrite, Ls, rb); + } + + MSG << "Benchmark Lime read" << std::endl; + MSG << SEP << std::endl; + for (auto &d: dir) + { + MSG << "-- Directory " << d << std::endl; + readBenchmark(GridDefaultLatt(), d + "/ioBench", limeRead, Ls, rb); + } + + Grid_finalize(); + + return EXIT_SUCCESS; +} From 3023287fd98cebc086d5f9bcd82e316394200722 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 16 Oct 2018 14:44:14 +0100 Subject: [PATCH 12/37] Hadrons: 3-index RO access to Eigen disk vector --- Hadrons/DiskVector.hpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/Hadrons/DiskVector.hpp b/Hadrons/DiskVector.hpp index 94c3e597..0b4bffe1 100644 --- a/Hadrons/DiskVector.hpp +++ b/Hadrons/DiskVector.hpp @@ -143,6 +143,12 @@ class EigenDiskVector: public DiskVectorBase> public: using DiskVectorBase>::DiskVectorBase; typedef EigenDiskVectorMat Matrix; +public: + T operator()(const unsigned int i, const Eigen::Index j, + const Eigen::Index k) const + { + return (*this)[i](j, k); + } private: virtual void load(EigenDiskVectorMat &obj, const std::string filename) const { From 109c74bed83be98962babe9ded39c1870249d01d Mon Sep 17 00:00:00 2001 From: Vera Guelpers Date: Tue, 16 Oct 2018 14:56:12 +0100 Subject: [PATCH 13/37] Hadrons: full volume noise source for A2A --- Hadrons/DilutedNoise.hpp | 75 +++++++++++ Hadrons/Modules.hpp | 1 + .../MNoise/FullVolumeSpinColorDiagonal.cc | 36 ++++++ .../MNoise/FullVolumeSpinColorDiagonal.hpp | 121 ++++++++++++++++++ Hadrons/modules.inc | 2 + 5 files changed, 235 insertions(+) create mode 100644 Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.cc create mode 100644 Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp diff --git a/Hadrons/DilutedNoise.hpp b/Hadrons/DilutedNoise.hpp index 1f8186ae..9440eb74 100644 --- a/Hadrons/DilutedNoise.hpp +++ b/Hadrons/DilutedNoise.hpp @@ -7,6 +7,7 @@ Source file: Hadrons/DilutedNoise.hpp Copyright (C) 2015-2018 Author: Antonin Portelli +Author: Vera Guelpers 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 @@ -55,6 +56,7 @@ public: GridCartesian *getGrid(void) const; // generate noise (pure virtual) virtual void generateNoise(GridParallelRNG &rng) = 0; + virtual void generateNoise(GridParallelRNG &rng, unsigned int n_src) = 0; private: std::vector noise_; GridCartesian *grid_; @@ -72,10 +74,28 @@ public: virtual ~TimeDilutedSpinColorDiagonalNoise(void) = default; // generate noise virtual void generateNoise(GridParallelRNG &rng); + virtual void generateNoise(GridParallelRNG &rng, unsigned int n_src); private: unsigned int nt_; }; +template +class FullVolumeSpinColorDiagonalNoise: public DilutedNoise +{ +public: + typedef typename FImpl::FermionField FermionField; +public: + // constructor/destructor + FullVolumeSpinColorDiagonalNoise(GridCartesian *g, unsigned int n_src); + virtual ~FullVolumeSpinColorDiagonalNoise(void) = default; + // generate noise + virtual void generateNoise(GridParallelRNG &rng); + virtual void generateNoise(GridParallelRNG &rng, unsigned int n_src); +public: + unsigned int n_src; +}; + + /****************************************************************************** * DilutedNoise template implementation * ******************************************************************************/ @@ -186,6 +206,61 @@ void TimeDilutedSpinColorDiagonalNoise::generateNoise(GridParallelRNG &rn } } +template +void TimeDilutedSpinColorDiagonalNoise::generateNoise(GridParallelRNG &rng, unsigned int n_src) +{ + assert(0); +} + +/****************************************************************************** + * FullVolumeSpinColorDiagonalNoise template implementation * + ******************************************************************************/ +template +FullVolumeSpinColorDiagonalNoise:: +FullVolumeSpinColorDiagonalNoise(GridCartesian *g, unsigned int n_src) +: DilutedNoise(g) +{ + this->resize(n_src*Ns*FImpl::Dimension); +} + +template +void FullVolumeSpinColorDiagonalNoise::generateNoise(GridParallelRNG &rng) +{ + assert(0); +} + +template +void FullVolumeSpinColorDiagonalNoise::generateNoise(GridParallelRNG &rng,unsigned int n_src) +{ + typedef decltype(peekColour((*this)[0], 0)) SpinField; + + auto &noise = *this; + auto g = this->getGrid(); + auto nd = g->GlobalDimensions().size(); + auto nc = FImpl::Dimension; + Complex shift(1., 1.); + LatticeComplex eta(g); + SpinField etas(g); + unsigned int i = 0; + + bernoulli(rng, eta); + eta = (2.*eta - shift)*(1./::sqrt(2.)); + for (unsigned int n = 0; n < n_src; ++n) + { + for (unsigned int s = 0; s < Ns; ++s) + { + etas = zero; + pokeSpin(etas, eta, s); + for (unsigned int c = 0; c < nc; ++c) + { + noise[i] = zero; + pokeColour(noise[i], etas, c); + i++; + } + } + } +} + END_HADRONS_NAMESPACE #endif // Hadrons_DilutedNoise_hpp_ diff --git a/Hadrons/Modules.hpp b/Hadrons/Modules.hpp index d86aed3c..682cb87a 100644 --- a/Hadrons/Modules.hpp +++ b/Hadrons/Modules.hpp @@ -32,6 +32,7 @@ #include #include #include +#include #include #include #include diff --git a/Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.cc b/Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.cc new file mode 100644 index 00000000..b909bf92 --- /dev/null +++ b/Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.cc @@ -0,0 +1,36 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.cc + +Copyright (C) 2015-2018 + +Author: Antonin Portelli +Author: Vera Guelpers + +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 MNoise; + +template class Grid::Hadrons::MNoise::TFullVolumeSpinColorDiagonal; +template class Grid::Hadrons::MNoise::TFullVolumeSpinColorDiagonal; diff --git a/Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp b/Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp new file mode 100644 index 00000000..3a4cc639 --- /dev/null +++ b/Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp @@ -0,0 +1,121 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp + +Copyright (C) 2015-2018 + +Author: Antonin Portelli +Author: Vera Guelpers + +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_MNoise_FullVolumeSpinColorDiagonal_hpp_ +#define Hadrons_MNoise_FullVolumeSpinColorDiagonal_hpp_ + +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Generate full volume spin-color diagonal noise * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MNoise) + +class FullVolumeSpinColorDiagonalPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(FullVolumeSpinColorDiagonalPar, + unsigned int, nsrc); +}; + +template +class TFullVolumeSpinColorDiagonal: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl,); +public: + // constructor + TFullVolumeSpinColorDiagonal(const std::string name); + // destructor + virtual ~TFullVolumeSpinColorDiagonal(void) {}; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_TMP(FullVolumeSpinColorDiagonal, TFullVolumeSpinColorDiagonal, MNoise); +MODULE_REGISTER_TMP(ZFullVolumeSpinColorDiagonal, TFullVolumeSpinColorDiagonal, MNoise); + +/****************************************************************************** + * TFullVolumeSpinColorDiagonal implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TFullVolumeSpinColorDiagonal::TFullVolumeSpinColorDiagonal(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TFullVolumeSpinColorDiagonal::getInput(void) +{ + std::vector in; + + return in; +} + +template +std::vector TFullVolumeSpinColorDiagonal::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TFullVolumeSpinColorDiagonal::setup(void) +{ + envCreateDerived(DilutedNoise, + FullVolumeSpinColorDiagonalNoise, + getName(), 1, envGetGrid(FermionField), par().nsrc); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TFullVolumeSpinColorDiagonal::execute(void) +{ + auto &noise = envGet(DilutedNoise, getName()); + LOG(Message) << "Generating full volume, spin-color diagonal noise" << std::endl; + noise.generateNoise(rng4d(),par().nsrc); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MNoise_FullVolumeSpinColorDiagonal_hpp_ diff --git a/Hadrons/modules.inc b/Hadrons/modules.inc index 6eeb891d..c7c8cc67 100644 --- a/Hadrons/modules.inc +++ b/Hadrons/modules.inc @@ -31,6 +31,7 @@ modules_cc =\ Modules/MGauge/FundtoHirep.cc \ Modules/MGauge/GaugeFix.cc \ Modules/MNoise/TimeDilutedSpinColorDiagonal.cc \ + Modules/MNoise/FullVolumeSpinColorDiagonal.cc \ Modules/MUtilities/RandomVectors.cc \ Modules/MUtilities/TestSeqGamma.cc \ Modules/MUtilities/PrecisionCast.cc \ @@ -101,6 +102,7 @@ modules_hpp =\ Modules/MGauge/FundtoHirep.hpp \ Modules/MGauge/StochEm.hpp \ Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp \ + Modules/MNoise/FullVolumeSpinColorDiagonal.hpp \ Modules/MUtilities/PrecisionCast.hpp \ Modules/MUtilities/RandomVectors.hpp \ Modules/MUtilities/TestSeqGamma.hpp \ From 0ba3d469c784c2a7bab5fc33f38a0ad9f1bb5020 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Wed, 17 Oct 2018 20:27:34 +0100 Subject: [PATCH 14/37] Benchmark IO in single and double precision --- benchmarks/Benchmark_IO_vs_dir.cc | 26 +++++++++++++++++++++----- 1 file changed, 21 insertions(+), 5 deletions(-) diff --git a/benchmarks/Benchmark_IO_vs_dir.cc b/benchmarks/Benchmark_IO_vs_dir.cc index c22550c3..15f6194f 100644 --- a/benchmarks/Benchmark_IO_vs_dir.cc +++ b/benchmarks/Benchmark_IO_vs_dir.cc @@ -3,9 +3,6 @@ #define MSG std::cout << GridLogMessage #define SEP \ "=============================================================================" -#ifndef BENCH_IO_LMAX -#define BENCH_IO_LMAX 40 -#endif using namespace Grid; using namespace QCD; @@ -41,7 +38,7 @@ int main (int argc, char ** argv) int64_t threads = GridThread::GetThreads(); MSG << "Grid is setup to use " << threads << " threads" << std::endl; MSG << SEP << std::endl; - MSG << "Benchmark Lime write" << std::endl; + MSG << "Benchmark double precision Lime write" << std::endl; MSG << SEP << std::endl; for (auto &d: dir) { @@ -49,7 +46,8 @@ int main (int argc, char ** argv) writeBenchmark(GridDefaultLatt(), d + "/ioBench", limeWrite, Ls, rb); } - MSG << "Benchmark Lime read" << std::endl; + MSG << SEP << std::endl; + MSG << "Benchmark double precision Lime read" << std::endl; MSG << SEP << std::endl; for (auto &d: dir) { @@ -57,6 +55,24 @@ int main (int argc, char ** argv) readBenchmark(GridDefaultLatt(), d + "/ioBench", limeRead, Ls, rb); } + MSG << SEP << std::endl; + MSG << "Benchmark single precision Lime write" << std::endl; + MSG << SEP << std::endl; + for (auto &d: dir) + { + MSG << "-- Directory " << d << std::endl; + writeBenchmark(GridDefaultLatt(), d + "/ioBench", limeWrite, Ls, rb); + } + + MSG << SEP << std::endl; + MSG << "Benchmark single precision Lime read" << std::endl; + MSG << SEP << std::endl; + for (auto &d: dir) + { + MSG << "-- Directory " << d << std::endl; + readBenchmark(GridDefaultLatt(), d + "/ioBench", limeRead, Ls, rb); + } + Grid_finalize(); return EXIT_SUCCESS; From b1c3cbe35e77fb10b36cb993a9010d3ef78b8872 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 18 Oct 2018 19:44:58 +0100 Subject: [PATCH 15/37] Hadrons: A2A vectors I/O --- Hadrons/A2AVectors.hpp | 122 ++++++++++++++++++++++++- Hadrons/Modules.hpp | 1 + Hadrons/Modules/MIO/LoadA2AVectors.cc | 7 ++ Hadrons/Modules/MIO/LoadA2AVectors.hpp | 93 +++++++++++++++++++ Hadrons/Modules/MSolver/A2AVectors.hpp | 11 ++- Hadrons/modules.inc | 4 +- 6 files changed, 235 insertions(+), 3 deletions(-) create mode 100644 Hadrons/Modules/MIO/LoadA2AVectors.cc create mode 100644 Hadrons/Modules/MIO/LoadA2AVectors.hpp diff --git a/Hadrons/A2AVectors.hpp b/Hadrons/A2AVectors.hpp index 423a5f08..ff32ff5e 100644 --- a/Hadrons/A2AVectors.hpp +++ b/Hadrons/A2AVectors.hpp @@ -36,7 +36,7 @@ See the full license in the file "LICENSE" in the top level distribution directo BEGIN_HADRONS_NAMESPACE /****************************************************************************** - * Classes to generate V & W all-to-all vectors * + * Class to generate V & W all-to-all vectors * ******************************************************************************/ template class A2AVectorsSchurDiagTwo @@ -70,6 +70,42 @@ private: SchurDiagTwoOperator op_; }; +/****************************************************************************** + * Methods for V & W all-to-all vectors I/O * + ******************************************************************************/ +class A2AVectorsIo +{ +public: + struct Record: Serializable + { + GRID_SERIALIZABLE_CLASS_MEMBERS(Record, + unsigned int, index); + Record(void): index(0) {} + }; +public: + template + static void write(const std::string fileStem, std::vector &vec, + const bool multiFile, const int trajectory = -1); + template + static void read(std::vector &vec, const std::string fileStem, + const bool multiFile, const int trajectory = -1); +private: + static inline std::string vecFilename(const std::string stem, const int traj, + const bool multiFile) + { + std::string t = (traj < 0) ? "" : ("." + std::to_string(traj)); + + if (multiFile) + { + return stem + t; + } + else + { + return stem + t + ".bin"; + } + } +}; + /****************************************************************************** * A2AVectorsSchurDiagTwo template implementation * ******************************************************************************/ @@ -217,6 +253,90 @@ void A2AVectorsSchurDiagTwo::makeHighModeW5D(FermionField &wout_4d, } } +/****************************************************************************** + * all-to-all vectors I/O template implementation * + ******************************************************************************/ +template +void A2AVectorsIo::write(const std::string fileStem, std::vector &vec, + const bool multiFile, const int trajectory) +{ + Record record; + GridBase *grid = vec[0]._grid; + ScidacWriter binWriter(grid->IsBoss()); + std::string filename = vecFilename(fileStem, multiFile, trajectory); + + if (multiFile) + { + std::string fullFilename; + + for (unsigned int i = 0; i < vec.size(); ++i) + { + fullFilename = filename + "/elem" + std::to_string(i) + ".bin"; + + LOG(Message) << "Writing vector " << i << std::endl; + makeFileDir(fullFilename, grid); + binWriter.open(fullFilename); + record.index = i; + binWriter.writeScidacFieldRecord(vec[i], record); + binWriter.close(); + } + } + else + { + makeFileDir(filename, grid); + binWriter.open(filename); + for (unsigned int i = 0; i < vec.size(); ++i) + { + LOG(Message) << "Writing vector " << i << std::endl; + record.index = i; + binWriter.writeScidacFieldRecord(vec[i], record); + } + binWriter.close(); + } +} + +template +void A2AVectorsIo::read(std::vector &vec, const std::string fileStem, + const bool multiFile, const int trajectory) +{ + Record record; + ScidacReader binReader; + std::string filename = vecFilename(fileStem, multiFile, trajectory); + + if (multiFile) + { + std::string fullFilename; + + for (unsigned int i = 0; i < vec.size(); ++i) + { + fullFilename = filename + "/elem" + std::to_string(i) + ".bin"; + + LOG(Message) << "Reading vector " << i << std::endl; + binReader.open(fullFilename); + binReader.readScidacFieldRecord(vec[i], record); + binReader.close(); + if (record.index != i) + { + HADRONS_ERROR(Io, "vector index mismatch"); + } + } + } + else + { + binReader.open(filename); + for (unsigned int i = 0; i < vec.size(); ++i) + { + LOG(Message) << "Reading vector " << i << std::endl; + binReader.readScidacFieldRecord(vec[i], record); + if (record.index != i) + { + HADRONS_ERROR(Io, "vector index mismatch"); + } + } + binReader.close(); + } +} + END_HADRONS_NAMESPACE #endif // A2A_Vectors_hpp_ diff --git a/Hadrons/Modules.hpp b/Hadrons/Modules.hpp index d86aed3c..a866c5cd 100644 --- a/Hadrons/Modules.hpp +++ b/Hadrons/Modules.hpp @@ -65,6 +65,7 @@ #include #include #include +#include #include #include #include diff --git a/Hadrons/Modules/MIO/LoadA2AVectors.cc b/Hadrons/Modules/MIO/LoadA2AVectors.cc new file mode 100644 index 00000000..4951b441 --- /dev/null +++ b/Hadrons/Modules/MIO/LoadA2AVectors.cc @@ -0,0 +1,7 @@ +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MIO; + +template class Grid::Hadrons::MIO::TLoadA2AVectors; diff --git a/Hadrons/Modules/MIO/LoadA2AVectors.hpp b/Hadrons/Modules/MIO/LoadA2AVectors.hpp new file mode 100644 index 00000000..cfb0cac1 --- /dev/null +++ b/Hadrons/Modules/MIO/LoadA2AVectors.hpp @@ -0,0 +1,93 @@ +#ifndef Hadrons_MIO_LoadA2AVectors_hpp_ +#define Hadrons_MIO_LoadA2AVectors_hpp_ + +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Module to load all-to-all vectors * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MIO) + +class LoadA2AVectorsPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(LoadA2AVectorsPar, + std::string, filestem, + bool, multiFile, + unsigned int, size); +}; + +template +class TLoadA2AVectors: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl,); +public: + // constructor + TLoadA2AVectors(const std::string name); + // destructor + virtual ~TLoadA2AVectors(void) {}; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_TMP(LoadA2AVectors, TLoadA2AVectors, MIO); + +/****************************************************************************** + * TLoadA2AVectors implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TLoadA2AVectors::TLoadA2AVectors(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TLoadA2AVectors::getInput(void) +{ + std::vector in; + + return in; +} + +template +std::vector TLoadA2AVectors::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TLoadA2AVectors::setup(void) +{ + envCreate(std::vector, getName(), 1, par().size, + envGetGrid(FermionField)); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TLoadA2AVectors::execute(void) +{ + auto &vec = envGet(std::vector, getName()); + + A2AVectorsIo::read(vec, par().filestem, par().multiFile, vm().getTrajectory()); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MIO_LoadA2AVectors_hpp_ diff --git a/Hadrons/Modules/MSolver/A2AVectors.hpp b/Hadrons/Modules/MSolver/A2AVectors.hpp index ee7cfa30..47ae6aee 100644 --- a/Hadrons/Modules/MSolver/A2AVectors.hpp +++ b/Hadrons/Modules/MSolver/A2AVectors.hpp @@ -51,7 +51,9 @@ public: std::string, noise, std::string, action, std::string, eigenPack, - std::string, solver); + std::string, solver, + std::string, output, + bool, multiFile); }; template @@ -236,6 +238,13 @@ void TA2AVectors::execute(void) } stopTimer("W high mode"); } + + // I/O if necessary + if (!par().output.empty()) + { + A2AVectorsIo::write(par().output + "_w", w, par().multiFile, vm().getTrajectory()); + A2AVectorsIo::write(par().output + "_v", v, par().multiFile, vm().getTrajectory()); + } } END_MODULE_NAMESPACE diff --git a/Hadrons/modules.inc b/Hadrons/modules.inc index 6eeb891d..c716b470 100644 --- a/Hadrons/modules.inc +++ b/Hadrons/modules.inc @@ -64,7 +64,8 @@ modules_cc =\ Modules/MIO/LoadBinary.cc \ Modules/MIO/LoadNersc.cc \ Modules/MIO/LoadCoarseEigenPack.cc \ - Modules/MIO/LoadCosmHol.cc + Modules/MIO/LoadCosmHol.cc \ + Modules/MIO/LoadA2AVectors.cc modules_hpp =\ Modules/MContraction/Baryon.hpp \ @@ -134,6 +135,7 @@ modules_hpp =\ Modules/MScalarSUN/TrKinetic.hpp \ Modules/MIO/LoadEigenPack.hpp \ Modules/MIO/LoadNersc.hpp \ + Modules/MIO/LoadA2AVectors.hpp \ Modules/MIO/LoadCosmHol.hpp \ Modules/MIO/LoadCoarseEigenPack.hpp \ Modules/MIO/LoadBinary.hpp From a3ace57e014fb7bb0f9da2304fe0f12983ef3f47 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 18 Oct 2018 19:46:11 +0100 Subject: [PATCH 16/37] Hadrons copyright update --- Hadrons/Modules/MIO/LoadA2AVectors.cc | 27 ++++++++++++++++++++++++++ Hadrons/Modules/MIO/LoadA2AVectors.hpp | 27 ++++++++++++++++++++++++++ Hadrons/Modules/MNPR/Amputate.cc | 1 + Hadrons/Modules/MNPR/Bilinear.cc | 1 + Hadrons/Modules/MNPR/FourQuark.cc | 1 + 5 files changed, 57 insertions(+) diff --git a/Hadrons/Modules/MIO/LoadA2AVectors.cc b/Hadrons/Modules/MIO/LoadA2AVectors.cc index 4951b441..7a40a6f5 100644 --- a/Hadrons/Modules/MIO/LoadA2AVectors.cc +++ b/Hadrons/Modules/MIO/LoadA2AVectors.cc @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MIO/LoadA2AVectors.cc + +Copyright (C) 2015-2018 + +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; diff --git a/Hadrons/Modules/MIO/LoadA2AVectors.hpp b/Hadrons/Modules/MIO/LoadA2AVectors.hpp index cfb0cac1..5b194c16 100644 --- a/Hadrons/Modules/MIO/LoadA2AVectors.hpp +++ b/Hadrons/Modules/MIO/LoadA2AVectors.hpp @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MIO/LoadA2AVectors.hpp + +Copyright (C) 2015-2018 + +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_MIO_LoadA2AVectors_hpp_ #define Hadrons_MIO_LoadA2AVectors_hpp_ diff --git a/Hadrons/Modules/MNPR/Amputate.cc b/Hadrons/Modules/MNPR/Amputate.cc index 8d55accb..ec7c5940 100644 --- a/Hadrons/Modules/MNPR/Amputate.cc +++ b/Hadrons/Modules/MNPR/Amputate.cc @@ -6,6 +6,7 @@ Source file: Hadrons/Modules/MNPR/Amputate.cc Copyright (C) 2015-2018 +Author: Antonin Portelli Author: Peter Boyle This program is free software; you can redistribute it and/or modify diff --git a/Hadrons/Modules/MNPR/Bilinear.cc b/Hadrons/Modules/MNPR/Bilinear.cc index 865739dd..c5b38d37 100644 --- a/Hadrons/Modules/MNPR/Bilinear.cc +++ b/Hadrons/Modules/MNPR/Bilinear.cc @@ -6,6 +6,7 @@ Source file: Hadrons/Modules/MNPR/Bilinear.cc Copyright (C) 2015-2018 +Author: Antonin Portelli Author: Peter Boyle This program is free software; you can redistribute it and/or modify diff --git a/Hadrons/Modules/MNPR/FourQuark.cc b/Hadrons/Modules/MNPR/FourQuark.cc index 9f14ff0f..1943c25b 100644 --- a/Hadrons/Modules/MNPR/FourQuark.cc +++ b/Hadrons/Modules/MNPR/FourQuark.cc @@ -6,6 +6,7 @@ Source file: Hadrons/Modules/MNPR/FourQuark.cc Copyright (C) 2015-2018 +Author: Antonin Portelli Author: Peter Boyle This program is free software; you can redistribute it and/or modify From 21304e21398e74f3d2c556bfe583e53b99be9279 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 18 Oct 2018 19:58:50 +0100 Subject: [PATCH 17/37] Hadrons: fix to allow single-prec build again --- Hadrons/Modules/MAction/DWF.cc | 2 ++ Hadrons/Modules/MAction/DWF.hpp | 2 ++ Hadrons/Modules/MAction/MobiusDWF.cc | 2 ++ Hadrons/Modules/MAction/MobiusDWF.hpp | 2 ++ Hadrons/Modules/MAction/ScaledDWF.cc | 2 ++ Hadrons/Modules/MAction/ScaledDWF.hpp | 2 ++ Hadrons/Modules/MAction/Wilson.cc | 2 ++ Hadrons/Modules/MAction/Wilson.hpp | 2 ++ Hadrons/Modules/MAction/WilsonClover.cc | 2 ++ Hadrons/Modules/MAction/WilsonClover.hpp | 2 ++ Hadrons/Modules/MAction/ZMobiusDWF.cc | 2 ++ Hadrons/Modules/MAction/ZMobiusDWF.hpp | 2 ++ 12 files changed, 24 insertions(+) diff --git a/Hadrons/Modules/MAction/DWF.cc b/Hadrons/Modules/MAction/DWF.cc index 6cda84ac..38d25cb9 100644 --- a/Hadrons/Modules/MAction/DWF.cc +++ b/Hadrons/Modules/MAction/DWF.cc @@ -32,4 +32,6 @@ using namespace Hadrons; using namespace MAction; template class Grid::Hadrons::MAction::TDWF; +#ifdef GRID_DEFAULT_PRECISION_DOUBLE template class Grid::Hadrons::MAction::TDWF; +#endif diff --git a/Hadrons/Modules/MAction/DWF.hpp b/Hadrons/Modules/MAction/DWF.hpp index 67cfeb1b..257782a1 100644 --- a/Hadrons/Modules/MAction/DWF.hpp +++ b/Hadrons/Modules/MAction/DWF.hpp @@ -73,7 +73,9 @@ protected: }; MODULE_REGISTER_TMP(DWF, TDWF, MAction); +#ifdef GRID_DEFAULT_PRECISION_DOUBLE MODULE_REGISTER_TMP(DWFF, TDWF, MAction); +#endif /****************************************************************************** * DWF template implementation * diff --git a/Hadrons/Modules/MAction/MobiusDWF.cc b/Hadrons/Modules/MAction/MobiusDWF.cc index 41683920..879452d8 100644 --- a/Hadrons/Modules/MAction/MobiusDWF.cc +++ b/Hadrons/Modules/MAction/MobiusDWF.cc @@ -32,4 +32,6 @@ using namespace Hadrons; using namespace MAction; template class Grid::Hadrons::MAction::TMobiusDWF; +#ifdef GRID_DEFAULT_PRECISION_DOUBLE template class Grid::Hadrons::MAction::TMobiusDWF; +#endif diff --git a/Hadrons/Modules/MAction/MobiusDWF.hpp b/Hadrons/Modules/MAction/MobiusDWF.hpp index 22964c9a..01223267 100644 --- a/Hadrons/Modules/MAction/MobiusDWF.hpp +++ b/Hadrons/Modules/MAction/MobiusDWF.hpp @@ -72,7 +72,9 @@ public: }; MODULE_REGISTER_TMP(MobiusDWF, TMobiusDWF, MAction); +#ifdef GRID_DEFAULT_PRECISION_DOUBLE MODULE_REGISTER_TMP(MobiusDWFF, TMobiusDWF, MAction); +#endif /****************************************************************************** * TMobiusDWF implementation * diff --git a/Hadrons/Modules/MAction/ScaledDWF.cc b/Hadrons/Modules/MAction/ScaledDWF.cc index 70eafac5..7008bf5d 100644 --- a/Hadrons/Modules/MAction/ScaledDWF.cc +++ b/Hadrons/Modules/MAction/ScaledDWF.cc @@ -32,4 +32,6 @@ using namespace Hadrons; using namespace MAction; template class Grid::Hadrons::MAction::TScaledDWF; +#ifdef GRID_DEFAULT_PRECISION_DOUBLE template class Grid::Hadrons::MAction::TScaledDWF; +#endif diff --git a/Hadrons/Modules/MAction/ScaledDWF.hpp b/Hadrons/Modules/MAction/ScaledDWF.hpp index c890f0e9..f7a09eef 100644 --- a/Hadrons/Modules/MAction/ScaledDWF.hpp +++ b/Hadrons/Modules/MAction/ScaledDWF.hpp @@ -71,7 +71,9 @@ public: }; MODULE_REGISTER_TMP(ScaledDWF, TScaledDWF, MAction); +#ifdef GRID_DEFAULT_PRECISION_DOUBLE MODULE_REGISTER_TMP(ScaledDWFF, TScaledDWF, MAction); +#endif /****************************************************************************** * TScaledDWF implementation * diff --git a/Hadrons/Modules/MAction/Wilson.cc b/Hadrons/Modules/MAction/Wilson.cc index 4ce3e7ef..1e801ed6 100644 --- a/Hadrons/Modules/MAction/Wilson.cc +++ b/Hadrons/Modules/MAction/Wilson.cc @@ -32,4 +32,6 @@ using namespace Hadrons; using namespace MAction; template class Grid::Hadrons::MAction::TWilson; +#ifdef GRID_DEFAULT_PRECISION_DOUBLE template class Grid::Hadrons::MAction::TWilson; +#endif diff --git a/Hadrons/Modules/MAction/Wilson.hpp b/Hadrons/Modules/MAction/Wilson.hpp index a9327c2f..093449e6 100644 --- a/Hadrons/Modules/MAction/Wilson.hpp +++ b/Hadrons/Modules/MAction/Wilson.hpp @@ -71,7 +71,9 @@ protected: }; MODULE_REGISTER_TMP(Wilson, TWilson, MAction); +#ifdef GRID_DEFAULT_PRECISION_DOUBLE MODULE_REGISTER_TMP(WilsonF, TWilson, MAction); +#endif /****************************************************************************** * TWilson template implementation * diff --git a/Hadrons/Modules/MAction/WilsonClover.cc b/Hadrons/Modules/MAction/WilsonClover.cc index 2c5c0e66..eed1582c 100644 --- a/Hadrons/Modules/MAction/WilsonClover.cc +++ b/Hadrons/Modules/MAction/WilsonClover.cc @@ -32,4 +32,6 @@ using namespace Hadrons; using namespace MAction; template class Grid::Hadrons::MAction::TWilsonClover; +#ifdef GRID_DEFAULT_PRECISION_DOUBLE template class Grid::Hadrons::MAction::TWilsonClover; +#endif diff --git a/Hadrons/Modules/MAction/WilsonClover.hpp b/Hadrons/Modules/MAction/WilsonClover.hpp index 349abe84..0b78bb55 100644 --- a/Hadrons/Modules/MAction/WilsonClover.hpp +++ b/Hadrons/Modules/MAction/WilsonClover.hpp @@ -75,7 +75,9 @@ public: }; MODULE_REGISTER_TMP(WilsonClover, TWilsonClover, MAction); +#ifdef GRID_DEFAULT_PRECISION_DOUBLE MODULE_REGISTER_TMP(WilsonCloverF, TWilsonClover, MAction); +#endif /****************************************************************************** * TWilsonClover template implementation * diff --git a/Hadrons/Modules/MAction/ZMobiusDWF.cc b/Hadrons/Modules/MAction/ZMobiusDWF.cc index ef8e4799..609b76cc 100644 --- a/Hadrons/Modules/MAction/ZMobiusDWF.cc +++ b/Hadrons/Modules/MAction/ZMobiusDWF.cc @@ -32,4 +32,6 @@ using namespace Hadrons; using namespace MAction; template class Grid::Hadrons::MAction::TZMobiusDWF; +#ifdef GRID_DEFAULT_PRECISION_DOUBLE template class Grid::Hadrons::MAction::TZMobiusDWF; +#endif diff --git a/Hadrons/Modules/MAction/ZMobiusDWF.hpp b/Hadrons/Modules/MAction/ZMobiusDWF.hpp index f7959127..40b04d7f 100644 --- a/Hadrons/Modules/MAction/ZMobiusDWF.hpp +++ b/Hadrons/Modules/MAction/ZMobiusDWF.hpp @@ -73,7 +73,9 @@ public: }; MODULE_REGISTER_TMP(ZMobiusDWF, TZMobiusDWF, MAction); +#ifdef GRID_DEFAULT_PRECISION_DOUBLE MODULE_REGISTER_TMP(ZMobiusDWFF, TZMobiusDWF, MAction); +#endif /****************************************************************************** * TZMobiusDWF implementation * From 2d3916418efb66ba7e43a15646eeabbb7a2497af Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 18 Oct 2018 23:45:13 +0100 Subject: [PATCH 18/37] Hadrons: more precision fix --- Hadrons/Modules/MIO/LoadEigenPack.cc | 2 ++ Hadrons/Modules/MIO/LoadEigenPack.hpp | 2 ++ 2 files changed, 4 insertions(+) diff --git a/Hadrons/Modules/MIO/LoadEigenPack.cc b/Hadrons/Modules/MIO/LoadEigenPack.cc index f0ae63c3..28fdeb01 100644 --- a/Hadrons/Modules/MIO/LoadEigenPack.cc +++ b/Hadrons/Modules/MIO/LoadEigenPack.cc @@ -32,4 +32,6 @@ using namespace Hadrons; using namespace MIO; template class Grid::Hadrons::MIO::TLoadEigenPack>; +#ifdef GRID_DEFAULT_PRECISION_DOUBLE template class Grid::Hadrons::MIO::TLoadEigenPack>; +#endif diff --git a/Hadrons/Modules/MIO/LoadEigenPack.hpp b/Hadrons/Modules/MIO/LoadEigenPack.hpp index 2302b15d..016675c9 100644 --- a/Hadrons/Modules/MIO/LoadEigenPack.hpp +++ b/Hadrons/Modules/MIO/LoadEigenPack.hpp @@ -72,7 +72,9 @@ public: }; MODULE_REGISTER_TMP(LoadFermionEigenPack, TLoadEigenPack>, MIO); +#ifdef GRID_DEFAULT_PRECISION_DOUBLE MODULE_REGISTER_TMP(LoadFermionEigenPackIo32, ARG(TLoadEigenPack>), MIO); +#endif /****************************************************************************** * TLoadEigenPack implementation * From 2e2e5ce596c72b22632e019a3255ed4d63d8d5de Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 19 Oct 2018 20:36:32 +0100 Subject: [PATCH 19/37] SciDAC I/O print data checksums --- Grid/parallelIO/IldgIO.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Grid/parallelIO/IldgIO.h b/Grid/parallelIO/IldgIO.h index 16951a26..518a84a9 100644 --- a/Grid/parallelIO/IldgIO.h +++ b/Grid/parallelIO/IldgIO.h @@ -233,7 +233,8 @@ class GridLimeReader : public BinaryIO { // std::cout << " ReadLatticeObject from offset "< munge; BinaryIO::readLatticeObject< vobj, sobj >(field, filename, munge, offset, format,nersc_csum,scidac_csuma,scidac_csumb); - + std::cout << GridLogMessage << "SciDAC checksum A " << std::hex << scidac_csuma << std::dec << std::endl; + std::cout << GridLogMessage << "SciDAC checksum B " << std::hex << scidac_csumb << std::dec << std::endl; ///////////////////////////////////////////// // Insist checksum is next record ///////////////////////////////////////////// From 1982cc58dd5737e1649bfbcdc7d1c249544fa70f Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Sun, 21 Oct 2018 01:20:05 +0100 Subject: [PATCH 20/37] Hadrons: A2A vectors I/O filename fix --- Hadrons/A2AVectors.hpp | 4 ++-- Hadrons/Modules/MSolver/A2AVectors.hpp | 6 +++++- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/Hadrons/A2AVectors.hpp b/Hadrons/A2AVectors.hpp index ff32ff5e..f55eb6d7 100644 --- a/Hadrons/A2AVectors.hpp +++ b/Hadrons/A2AVectors.hpp @@ -263,7 +263,7 @@ void A2AVectorsIo::write(const std::string fileStem, std::vector &vec, Record record; GridBase *grid = vec[0]._grid; ScidacWriter binWriter(grid->IsBoss()); - std::string filename = vecFilename(fileStem, multiFile, trajectory); + std::string filename = vecFilename(fileStem, trajectory, multiFile); if (multiFile) { @@ -301,7 +301,7 @@ void A2AVectorsIo::read(std::vector &vec, const std::string fileStem, { Record record; ScidacReader binReader; - std::string filename = vecFilename(fileStem, multiFile, trajectory); + std::string filename = vecFilename(fileStem, trajectory, multiFile); if (multiFile) { diff --git a/Hadrons/Modules/MSolver/A2AVectors.hpp b/Hadrons/Modules/MSolver/A2AVectors.hpp index 47ae6aee..f9980ee3 100644 --- a/Hadrons/Modules/MSolver/A2AVectors.hpp +++ b/Hadrons/Modules/MSolver/A2AVectors.hpp @@ -242,8 +242,12 @@ void TA2AVectors::execute(void) // I/O if necessary if (!par().output.empty()) { - A2AVectorsIo::write(par().output + "_w", w, par().multiFile, vm().getTrajectory()); + startTimer("V I/O"); A2AVectorsIo::write(par().output + "_v", v, par().multiFile, vm().getTrajectory()); + stopTimer("V I/O"); + startTimer("W I/O"); + A2AVectorsIo::write(par().output + "_w", w, par().multiFile, vm().getTrajectory()); + stopTimer("W I/O"); } } From 6b559d68aab01e8e969f600bb9b233261a31f9c2 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Mon, 22 Oct 2018 11:10:18 +0100 Subject: [PATCH 21/37] Hadrons: eigenpack converter can do test reads --- Hadrons/Utilities/EigenPackCast.cc | 36 +++++++++++++++++++++++++----- 1 file changed, 31 insertions(+), 5 deletions(-) diff --git a/Hadrons/Utilities/EigenPackCast.cc b/Hadrons/Utilities/EigenPackCast.cc index 8f2755e3..8c3a2b3f 100644 --- a/Hadrons/Utilities/EigenPackCast.cc +++ b/Hadrons/Utilities/EigenPackCast.cc @@ -35,7 +35,7 @@ using namespace Hadrons; template void convert(const std::string outFilename, const std::string inFilename, const unsigned int Ls, const bool rb, const unsigned int size, - const bool multiFile) + const bool multiFile, const bool testRead) { assert(outFilename != inFilename); @@ -102,6 +102,7 @@ void convert(const std::string outFilename, const std::string inFilename, LOG(Message) << "Out type : " << typeName() << std::endl; LOG(Message) << "#vectors : " << size << std::endl; LOG(Message) << "Multifile : " << (multiFile ? "yes" : "no") << std::endl; + LOG(Message) << "Test read : " << (testRead ? "yes" : "no") << std::endl; if (multiFile) { for(unsigned int k = 0; k < size; ++k) @@ -112,6 +113,8 @@ void convert(const std::string outFilename, const std::string inFilename, LOG(Message) << "==== Converting vector " << k << std::endl; LOG(Message) << "In : " << inV << std::endl; LOG(Message) << "Out: " << outV << std::endl; + // conversion + LOG(Message) << "-- Doing conversion" << std::endl; makeFileDir(outV, gOut); binWriter.open(outV); binReader.open(inV); @@ -121,10 +124,20 @@ void convert(const std::string outFilename, const std::string inFilename, EigenPackIo::writeElement(binWriter, bufIn, eval, k, &bufOut, &testIn); binWriter.close(); binReader.close(); + // read test + if (testRead) + { + LOG(Message) << "-- Test read" << std::endl; + binReader.open(outV); + EigenPackIo::readElement(bufOut, eval, k, binReader); + binReader.close(); + } } } else { + // conversion + LOG(Message) << "-- Doing conversion" << std::endl; makeFileDir(outFilename, gOut); binWriter.open(outFilename); binReader.open(inFilename); @@ -137,6 +150,18 @@ void convert(const std::string outFilename, const std::string inFilename, } binWriter.close(); binReader.close(); + // read test + if (testRead) + { + LOG(Message) << "-- Test read" << std::endl; + binReader.open(outFilename); + EigenPackIo::readHeader(record, binReader); + for(unsigned int k = 0; k < size; ++k) + { + EigenPackIo::readElement(bufOut, eval, k, binReader); + } + binReader.close(); + } } } @@ -154,11 +179,11 @@ int main(int argc, char *argv[]) // parse command line std::string outFilename, inFilename; unsigned int size, Ls; - bool rb, multiFile; + bool rb, multiFile, testRead; - if (argc < 7) + if (argc < 8) { - std::cerr << "usage: " << argv[0] << " <#vector> [Grid options]"; + std::cerr << "usage: " << argv[0] << " <#vector> [Grid options]"; std::cerr << std::endl; std::exit(EXIT_FAILURE); } @@ -168,6 +193,7 @@ int main(int argc, char *argv[]) rb = (std::string(argv[4]) != "0"); size = std::stoi(std::string(argv[5])); multiFile = (std::string(argv[6]) != "0"); + testRead = (std::string(argv[7]) != "0"); // initialization Grid_init(&argc, &argv); @@ -176,7 +202,7 @@ int main(int argc, char *argv[]) // execution try { - convert(outFilename, inFilename, Ls, rb, size, multiFile); + convert(outFilename, inFilename, Ls, rb, size, multiFile, testRead); } catch (const std::exception& e) { From febe41cc1d5832869fbc71fd97beb6e3a46a5065 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 23 Oct 2018 12:48:15 +0100 Subject: [PATCH 22/37] Hadrons: improvement on PR #176 --- Hadrons/DilutedNoise.hpp | 29 ++++--------------- .../MNoise/FullVolumeSpinColorDiagonal.hpp | 2 +- 2 files changed, 7 insertions(+), 24 deletions(-) diff --git a/Hadrons/DilutedNoise.hpp b/Hadrons/DilutedNoise.hpp index 9440eb74..1eb7ff40 100644 --- a/Hadrons/DilutedNoise.hpp +++ b/Hadrons/DilutedNoise.hpp @@ -56,7 +56,6 @@ public: GridCartesian *getGrid(void) const; // generate noise (pure virtual) virtual void generateNoise(GridParallelRNG &rng) = 0; - virtual void generateNoise(GridParallelRNG &rng, unsigned int n_src) = 0; private: std::vector noise_; GridCartesian *grid_; @@ -74,7 +73,6 @@ public: virtual ~TimeDilutedSpinColorDiagonalNoise(void) = default; // generate noise virtual void generateNoise(GridParallelRNG &rng); - virtual void generateNoise(GridParallelRNG &rng, unsigned int n_src); private: unsigned int nt_; }; @@ -90,9 +88,8 @@ public: virtual ~FullVolumeSpinColorDiagonalNoise(void) = default; // generate noise virtual void generateNoise(GridParallelRNG &rng); - virtual void generateNoise(GridParallelRNG &rng, unsigned int n_src); -public: - unsigned int n_src; +private: + unsigned int nSrc_; }; @@ -206,31 +203,17 @@ void TimeDilutedSpinColorDiagonalNoise::generateNoise(GridParallelRNG &rn } } -template -void TimeDilutedSpinColorDiagonalNoise::generateNoise(GridParallelRNG &rng, unsigned int n_src) -{ - assert(0); -} - /****************************************************************************** * FullVolumeSpinColorDiagonalNoise template implementation * ******************************************************************************/ template FullVolumeSpinColorDiagonalNoise:: -FullVolumeSpinColorDiagonalNoise(GridCartesian *g, unsigned int n_src) -: DilutedNoise(g) -{ - this->resize(n_src*Ns*FImpl::Dimension); -} +FullVolumeSpinColorDiagonalNoise(GridCartesian *g, unsigned int nSrc) +: DilutedNoise(g, nSrc*Ns*FImpl::Dimension), nSrc_(nSrc) +{} template void FullVolumeSpinColorDiagonalNoise::generateNoise(GridParallelRNG &rng) -{ - assert(0); -} - -template -void FullVolumeSpinColorDiagonalNoise::generateNoise(GridParallelRNG &rng,unsigned int n_src) { typedef decltype(peekColour((*this)[0], 0)) SpinField; @@ -245,7 +228,7 @@ void FullVolumeSpinColorDiagonalNoise::generateNoise(GridParallelRNG &rng bernoulli(rng, eta); eta = (2.*eta - shift)*(1./::sqrt(2.)); - for (unsigned int n = 0; n < n_src; ++n) + for (unsigned int n = 0; n < nSrc_; ++n) { for (unsigned int s = 0; s < Ns; ++s) { diff --git a/Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp b/Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp index 3a4cc639..93990882 100644 --- a/Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp +++ b/Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp @@ -111,7 +111,7 @@ void TFullVolumeSpinColorDiagonal::execute(void) { auto &noise = envGet(DilutedNoise, getName()); LOG(Message) << "Generating full volume, spin-color diagonal noise" << std::endl; - noise.generateNoise(rng4d(),par().nsrc); + noise.generateNoise(rng4d()); } END_MODULE_NAMESPACE From 8f514ae55002cc267a823f300768ada6d8488102 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Mon, 5 Nov 2018 11:41:10 +0000 Subject: [PATCH 23/37] Hadrons: Lanczos 32bit IO --- .../Modules/MSolver/LocalCoherenceLanczos.cc | 5 ++- .../Modules/MSolver/LocalCoherenceLanczos.hpp | 34 +++++++++++-------- 2 files changed, 23 insertions(+), 16 deletions(-) diff --git a/Hadrons/Modules/MSolver/LocalCoherenceLanczos.cc b/Hadrons/Modules/MSolver/LocalCoherenceLanczos.cc index fa73e5b6..dacc871f 100644 --- a/Hadrons/Modules/MSolver/LocalCoherenceLanczos.cc +++ b/Hadrons/Modules/MSolver/LocalCoherenceLanczos.cc @@ -33,4 +33,7 @@ using namespace MSolver; template class Grid::Hadrons::MSolver::TLocalCoherenceLanczos; template class Grid::Hadrons::MSolver::TLocalCoherenceLanczos; - +#ifdef GRID_DEFAULT_PRECISION_DOUBLE +template class Grid::Hadrons::MSolver::TLocalCoherenceLanczos; +template class Grid::Hadrons::MSolver::TLocalCoherenceLanczos; +#endif diff --git a/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp b/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp index 7242eaf9..492ff39e 100644 --- a/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp +++ b/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp @@ -55,17 +55,17 @@ public: bool, multiFile); }; -template +template class TLocalCoherenceLanczos: public Module { public: FERM_TYPE_ALIASES(FImpl,); typedef LocalCoherenceLanczos LCL; - typedef BaseFermionEigenPack BasePack; - typedef CoarseFermionEigenPack CoarsePack; - typedef HADRONS_DEFAULT_SCHUR_OP SchurFMat; + nBasis> LCL; + typedef BaseFermionEigenPack BasePack; + typedef CoarseFermionEigenPack CoarsePack; + typedef HADRONS_DEFAULT_SCHUR_OP SchurFMat; public: // constructor TLocalCoherenceLanczos(const std::string name); @@ -82,27 +82,31 @@ public: MODULE_REGISTER_TMP(LocalCoherenceLanczos, ARG(TLocalCoherenceLanczos), MSolver); MODULE_REGISTER_TMP(ZLocalCoherenceLanczos, ARG(TLocalCoherenceLanczos), MSolver); +#ifdef GRID_DEFAULT_PRECISION_DOUBLE +MODULE_REGISTER_TMP(LocalCoherenceLanczosIo32, ARG(TLocalCoherenceLanczos), MSolver); +MODULE_REGISTER_TMP(ZLocalCoherenceLanczosIo32, ARG(TLocalCoherenceLanczos), MSolver); +#endif /****************************************************************************** * TLocalCoherenceLanczos implementation * ******************************************************************************/ // constructor ///////////////////////////////////////////////////////////////// -template -TLocalCoherenceLanczos::TLocalCoherenceLanczos(const std::string name) +template +TLocalCoherenceLanczos::TLocalCoherenceLanczos(const std::string name) : Module(name) {} // dependencies/products /////////////////////////////////////////////////////// -template -std::vector TLocalCoherenceLanczos::getInput(void) +template +std::vector TLocalCoherenceLanczos::getInput(void) { std::vector in = {par().action}; return in; } -template -std::vector TLocalCoherenceLanczos::getOutput(void) +template +std::vector TLocalCoherenceLanczos::getOutput(void) { std::vector out = {getName()}; @@ -110,8 +114,8 @@ std::vector TLocalCoherenceLanczos::getOutput(void) } // setup /////////////////////////////////////////////////////////////////////// -template -void TLocalCoherenceLanczos::setup(void) +template +void TLocalCoherenceLanczos::setup(void) { LOG(Message) << "Setting up local coherence Lanczos eigensolver for" << " action '" << par().action << "' (" << nBasis @@ -138,8 +142,8 @@ void TLocalCoherenceLanczos::setup(void) } // execution /////////////////////////////////////////////////////////////////// -template -void TLocalCoherenceLanczos::execute(void) +template +void TLocalCoherenceLanczos::execute(void) { auto &finePar = par().fineParams; auto &coarsePar = par().coarseParams; From 85699daef29603eba224a03a2e618cbf56cd8a83 Mon Sep 17 00:00:00 2001 From: Vera Guelpers Date: Tue, 6 Nov 2018 10:27:18 +0000 Subject: [PATCH 24/37] Hadrons: Module to electrify a gauge field --- Hadrons/Modules.hpp | 1 + Hadrons/Modules/MGauge/Electrify.cc | 34 ++++++ Hadrons/Modules/MGauge/Electrify.hpp | 151 +++++++++++++++++++++++++++ Hadrons/modules.inc | 2 + 4 files changed, 188 insertions(+) create mode 100644 Hadrons/Modules/MGauge/Electrify.cc create mode 100644 Hadrons/Modules/MGauge/Electrify.hpp diff --git a/Hadrons/Modules.hpp b/Hadrons/Modules.hpp index da57ab2b..b4d55f99 100644 --- a/Hadrons/Modules.hpp +++ b/Hadrons/Modules.hpp @@ -27,6 +27,7 @@ #include #include #include +#include #include #include #include diff --git a/Hadrons/Modules/MGauge/Electrify.cc b/Hadrons/Modules/MGauge/Electrify.cc new file mode 100644 index 00000000..1feea9ec --- /dev/null +++ b/Hadrons/Modules/MGauge/Electrify.cc @@ -0,0 +1,34 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MGauge/Electrify.cc + +Copyright (C) 2015-2018 + +Author: Vera Guelpers + +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; + +template class Grid::Hadrons::MGauge::TElectrify; diff --git a/Hadrons/Modules/MGauge/Electrify.hpp b/Hadrons/Modules/MGauge/Electrify.hpp new file mode 100644 index 00000000..533ac704 --- /dev/null +++ b/Hadrons/Modules/MGauge/Electrify.hpp @@ -0,0 +1,151 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MGauge/Electrify.hpp + +Copyright (C) 2015-2018 + +Author: Vera Guelpers + +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_MGauge_Electrify_hpp_ +#define Hadrons_MGauge_Electrify_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Electrify gauge * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MGauge) + +/**************************************************************************** +* Electrify a gauge field: +* +* Ue_mu(x) = U_mu(x)*exp(ieqA_mu(x)) +* +* with +* +* - gauge: U_mu(x): gauge field +* - photon: A_mu(x): electromagnetic photon field +* - e: value for the elementary charge +* - q: charge in units of e +* +*****************************************************************************/ + + +class ElectrifyPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(ElectrifyPar, + std::string, gauge, + std::string, photon, + double, e, + double, charge); +}; + +template +class TElectrify: public Module +{ +public: + GAUGE_TYPE_ALIASES(GImpl,); +public: + typedef PhotonR::GaugeField EmField; +public: + // constructor + TElectrify(const std::string name); + // destructor + virtual ~TElectrify(void) {}; + // dependencies/products + virtual std::vector getInput(void); + virtual std::vector getOutput(void); +protected: + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_TMP(Electrify, TElectrify, MGauge); + +/****************************************************************************** +* TElectrify implementation * +******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TElectrify::TElectrify(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TElectrify::getInput(void) +{ + std::vector in = {par().gauge, par().photon}; + + return in; +} + +template +std::vector TElectrify::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TElectrify::setup(void) +{ + envCreateLat(GaugeField, getName()); + envTmpLat(LatticeComplex, "eiAmu"); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TElectrify::execute(void) +{ + LOG(Message) << "Electrify the gauge field " << par().gauge << " using the photon field " + << par().photon << " with charge e*q= " << par().e << "*" << par().charge << std::endl; + + auto &Ue = envGet(GaugeField, getName()); + auto &U = envGet(GaugeField, par().gauge); + auto &A = envGet(EmField, par().photon); + envGetTmp(LatticeComplex, eiAmu); + + Complex i(0.0,1.0); + + for(unsigned int mu = 0; mu < env().getNd(); mu++) + { + eiAmu = exp(i * par().e * par().charge * PeekIndex(A, mu)); + PokeIndex(Ue, PeekIndex(U, mu) * eiAmu, mu); + } +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MGauge_Electrify_hpp_ diff --git a/Hadrons/modules.inc b/Hadrons/modules.inc index 157fe1a2..1573685b 100644 --- a/Hadrons/modules.inc +++ b/Hadrons/modules.inc @@ -25,6 +25,7 @@ modules_cc =\ Modules/MSolver/LocalCoherenceLanczos.cc \ Modules/MGauge/StoutSmearing.cc \ Modules/MGauge/Unit.cc \ + Modules/MGauge/Electrify.cc \ Modules/MGauge/UnitEm.cc \ Modules/MGauge/StochEm.cc \ Modules/MGauge/Random.cc \ @@ -98,6 +99,7 @@ modules_hpp =\ Modules/MGauge/UnitEm.hpp \ Modules/MGauge/StoutSmearing.hpp \ Modules/MGauge/Unit.hpp \ + Modules/MGauge/Electrify.hpp \ Modules/MGauge/Random.hpp \ Modules/MGauge/GaugeFix.hpp \ Modules/MGauge/FundtoHirep.hpp \ From 4a47b1187637890b89939186fa7d70a80aa585fe Mon Sep 17 00:00:00 2001 From: Azusa Yamaguchi Date: Tue, 6 Nov 2018 12:49:05 +0000 Subject: [PATCH 25/37] Block CG improvements to develop --- .../iterative/BlockConjugateGradient.h | 446 +++++++++++------- Grid/algorithms/iterative/ConjugateGradient.h | 16 +- tests/solver/Test_dwf_mrhs_cg_mpi.cc | 9 +- tests/solver/Test_mobius_bcg.cc | 220 +++++++++ tests/solver/Test_mobius_bcg_nosplit.cc | 144 ++++++ tests/solver/Test_mobius_bcg_phys_nosplit.cc | 148 ++++++ tests/solver/Test_mobius_bcg_prec_nosplit.cc | 144 ++++++ tests/solver/Test_staggered_block_cg_prec.cc | 65 ++- 8 files changed, 1003 insertions(+), 189 deletions(-) create mode 100644 tests/solver/Test_mobius_bcg.cc create mode 100644 tests/solver/Test_mobius_bcg_nosplit.cc create mode 100644 tests/solver/Test_mobius_bcg_phys_nosplit.cc create mode 100644 tests/solver/Test_mobius_bcg_prec_nosplit.cc diff --git a/Grid/algorithms/iterative/BlockConjugateGradient.h b/Grid/algorithms/iterative/BlockConjugateGradient.h index e0eeddcb..9a1e8b34 100644 --- a/Grid/algorithms/iterative/BlockConjugateGradient.h +++ b/Grid/algorithms/iterative/BlockConjugateGradient.h @@ -33,7 +33,7 @@ directory namespace Grid { -enum BlockCGtype { BlockCG, BlockCGrQ, CGmultiRHS }; +enum BlockCGtype { BlockCG, BlockCGrQ, CGmultiRHS, BlockCGVec, BlockCGrQVec }; ////////////////////////////////////////////////////////////////////////// // Block conjugate gradient. Dimension zero should be the block direction @@ -42,7 +42,6 @@ template class BlockConjugateGradient : public OperatorFunction { public: - typedef typename Field::scalar_type scomplex; int blockDim ; @@ -54,21 +53,15 @@ class BlockConjugateGradient : public OperatorFunction { RealD Tolerance; Integer MaxIterations; Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion + Integer PrintInterval; //GridLogMessages or Iterative BlockConjugateGradient(BlockCGtype cgtype,int _Orthog,RealD tol, Integer maxit, bool err_on_no_conv = true) - : Tolerance(tol), CGtype(cgtype), blockDim(_Orthog), MaxIterations(maxit), ErrorOnNoConverge(err_on_no_conv) + : Tolerance(tol), CGtype(cgtype), blockDim(_Orthog), MaxIterations(maxit), ErrorOnNoConverge(err_on_no_conv),PrintInterval(100) {}; //////////////////////////////////////////////////////////////////////////////////////////////////// // Thin QR factorisation (google it) //////////////////////////////////////////////////////////////////////////////////////////////////// -void ThinQRfact (Eigen::MatrixXcd &m_rr, - Eigen::MatrixXcd &C, - Eigen::MatrixXcd &Cinv, - Field & Q, - const Field & R) -{ - int Orthog = blockDim; // First dimension is block dim; this is an assumption //////////////////////////////////////////////////////////////////////////////////////////////////// //Dimensions // R_{ferm x Nblock} = Q_{ferm x Nblock} x C_{Nblock x Nblock} -> ferm x Nblock @@ -85,22 +78,20 @@ void ThinQRfact (Eigen::MatrixXcd &m_rr, // Cdag C = Rdag R ; passes. // QdagQ = 1 ; passes //////////////////////////////////////////////////////////////////////////////////////////////////// +void ThinQRfact (Eigen::MatrixXcd &m_rr, + Eigen::MatrixXcd &C, + Eigen::MatrixXcd &Cinv, + Field & Q, + const Field & R) +{ + int Orthog = blockDim; // First dimension is block dim; this is an assumption sliceInnerProductMatrix(m_rr,R,R,Orthog); // Force manifest hermitian to avoid rounding related m_rr = 0.5*(m_rr+m_rr.adjoint()); -#if 0 - std::cout << " Calling Cholesky ldlt on m_rr " << m_rr < & Q, + const std::vector & R) +{ + InnerProductMatrix(m_rr,R,R); + + m_rr = 0.5*(m_rr+m_rr.adjoint()); + + Eigen::MatrixXcd L = m_rr.llt().matrixL(); + + C = L.adjoint(); + Cinv = C.inverse(); + + MulMatrix(Q,Cinv,R); +} + //////////////////////////////////////////////////////////////////////////////////////////////////// // Call one of several implementations //////////////////////////////////////////////////////////////////////////////////////////////////// @@ -119,14 +129,20 @@ void operator()(LinearOperatorBase &Linop, const Field &Src, Field &Psi) { if ( CGtype == BlockCGrQ ) { BlockCGrQsolve(Linop,Src,Psi); - } else if (CGtype == BlockCG ) { - BlockCGsolve(Linop,Src,Psi); } else if (CGtype == CGmultiRHS ) { CGmultiRHSsolve(Linop,Src,Psi); } else { assert(0); } } +void operator()(LinearOperatorBase &Linop, const std::vector &Src, std::vector &Psi) +{ + if ( CGtype == BlockCGrQVec ) { + BlockCGrQsolveVec(Linop,Src,Psi); + } else { + assert(0); + } +} //////////////////////////////////////////////////////////////////////////// // BlockCGrQ implementation: @@ -139,7 +155,8 @@ void BlockCGrQsolve(LinearOperatorBase &Linop, const Field &B, Field &X) { int Orthog = blockDim; // First dimension is block dim; this is an assumption Nblock = B._grid->_fdimensions[Orthog]; - +/* FAKE */ + Nblock=8; std::cout< &Linop, const Field &B, Field &X) std::cout << GridLogMessage<<"BlockCGrQ algorithm initialisation " < Thin QR factorisation (google it) - Linop.HermOp(X, AD); tmp = B - AD; - //std::cout << GridLogMessage << " initial tmp " << norm2(tmp)<< std::endl; + ThinQRfact (m_rr, m_C, m_Cinv, Q, tmp); - //std::cout << GridLogMessage << " initial Q " << norm2(Q)<< std::endl; - //std::cout << GridLogMessage << " m_rr " << m_rr< &Linop, const Field &B, Field &X) MatrixTimer.Start(); Linop.HermOp(D, Z); MatrixTimer.Stop(); - //std::cout << GridLogMessage << " norm2 Z " < &Linop, const Field &B, Field &X) //7. D = Q + D S^dag m_tmp = m_S.adjoint(); + sliceMaddTimer.Start(); sliceMaddMatrix(D,m_tmp,D,Q,Orthog); sliceMaddTimer.Stop(); @@ -317,152 +328,6 @@ void BlockCGrQsolve(LinearOperatorBase &Linop, const Field &B, Field &X) IterationsToComplete = k; } ////////////////////////////////////////////////////////////////////////// -// Block conjugate gradient; Original O'Leary Dimension zero should be the block direction -////////////////////////////////////////////////////////////////////////// -void BlockCGsolve(LinearOperatorBase &Linop, const Field &Src, Field &Psi) -{ - int Orthog = blockDim; // First dimension is block dim; this is an assumption - Nblock = Src._grid->_fdimensions[Orthog]; - - std::cout< residuals(Nblock); - std::vector ssq(Nblock); - - sliceNorm(ssq,Src,Orthog); - RealD sssum=0; - for(int b=0;b max_resid ) max_resid = rr; - } - - if ( max_resid < Tolerance*Tolerance ) { - - SolverTimer.Stop(); - - std::cout << GridLogMessage<<"BlockCG converged in "< &Linop, const Field &Src, Field & IterationsToComplete = k; } +void InnerProductMatrix(Eigen::MatrixXcd &m , const std::vector &X, const std::vector &Y){ + for(int b=0;b &AP, Eigen::MatrixXcd &m , const std::vector &X,const std::vector &Y,RealD scale=1.0){ + // Should make this cache friendly with site outermost, parallel_for + // Deal with case AP aliases with either Y or X + std::vector tmp(Nblock,X[0]); + for(int b=0;b &AP, Eigen::MatrixXcd &m , const std::vector &X){ + // Should make this cache friendly with site outermost, parallel_for + for(int b=0;b &P){ + double nn = 0.0; + for(int b=0;b &Linop, const std::vector &B, std::vector &X) +{ + Nblock = B.size(); + assert(Nblock == X.size()); + + std::cout< tmp(Nblock,Fake); + std::vector Q(Nblock,Fake); + std::vector D(Nblock,Fake); + std::vector Z(Nblock,Fake); + std::vector AD(Nblock,Fake); + + Eigen::MatrixXcd m_DZ = Eigen::MatrixXcd::Identity(Nblock,Nblock); + Eigen::MatrixXcd m_M = Eigen::MatrixXcd::Identity(Nblock,Nblock); + Eigen::MatrixXcd m_rr = Eigen::MatrixXcd::Zero(Nblock,Nblock); + + Eigen::MatrixXcd m_C = Eigen::MatrixXcd::Zero(Nblock,Nblock); + Eigen::MatrixXcd m_Cinv = Eigen::MatrixXcd::Zero(Nblock,Nblock); + Eigen::MatrixXcd m_S = Eigen::MatrixXcd::Zero(Nblock,Nblock); + Eigen::MatrixXcd m_Sinv = Eigen::MatrixXcd::Zero(Nblock,Nblock); + + Eigen::MatrixXcd m_tmp = Eigen::MatrixXcd::Identity(Nblock,Nblock); + Eigen::MatrixXcd m_tmp1 = Eigen::MatrixXcd::Identity(Nblock,Nblock); + + // Initial residual computation & set up + std::vector residuals(Nblock); + std::vector ssq(Nblock); + + RealD sssum=0; + for(int b=0;b Thin QR factorisation (google it) + * for k: + * Z = AD + * M = [D^dag Z]^{-1} + * X = X + D MC + * QS = Q - ZM + * D = Q + D S^dag + * C = S C + */ + /////////////////////////////////////// + // Initial block: initial search dir is guess + /////////////////////////////////////// + std::cout << GridLogMessage<<"BlockCGrQvec algorithm initialisation " < Thin QR factorisation (google it) + for(int b=0;b max_resid ) max_resid = rr; + } + + std::cout << GridLogIterative << "\t Block Iteration "< { LinalgTimer.Stop(); std::cout << GridLogIterative << "ConjugateGradient: Iteration " << k - << " residual " << cp << " target " << rsq << std::endl; + << " residual^2 " << sqrt(cp/ssq) << " target " << Tolerance << std::endl; // Stopping condition if (cp <= rsq) { @@ -150,13 +150,13 @@ class ConjugateGradient : public OperatorFunction { std::cout << GridLogMessage << "\tTrue residual " << true_residual< HermOp(Ddwf); MdagMLinearOperator HermOpCk(Dchk); - ConjugateGradient CG((1.0e-2),10000); + ConjugateGradient CG((stp),10000); s_res = zero; CG(HermOp,s_src,s_res); @@ -227,5 +228,11 @@ int main (int argc, char ** argv) std::cout << GridLogMessage<<" resid["< BCGV (BlockCGVec,blockDim,stp,10000); + BCGV.PrintInterval=10; + BCGV(HermOpCk,src,result); + Grid_finalize(); } diff --git a/tests/solver/Test_mobius_bcg.cc b/tests/solver/Test_mobius_bcg.cc new file mode 100644 index 00000000..e59cb7e0 --- /dev/null +++ b/tests/solver/Test_mobius_bcg.cc @@ -0,0 +1,220 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_dwf_mrhs_cg.cc + + Copyright (C) 2015 + +Author: Peter Boyle + + 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 std; +using namespace Grid; +using namespace Grid::QCD; + +int main (int argc, char ** argv) +{ + typedef typename MobiusFermionR::FermionField FermionField; + typedef typename MobiusFermionR::ComplexField ComplexField; + typename MobiusFermionR::ImplParams params; + + const int Ls=12; + + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + std::vector mpi_split (mpi_layout.size(),1); + std::vector split_coor (mpi_layout.size(),1); + std::vector split_dim (mpi_layout.size(),1); + + std::vector boundary_phases(Nd,1.); + boundary_phases[Nd-1]=-1.; + params.boundary_phases = boundary_phases; + + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), + GridDefaultSimd(Nd,vComplex::Nsimd()), + GridDefaultMpi()); + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + GridRedBlackCartesian * rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + + ///////////////////////////////////////////// + // Split into 1^4 mpi communicators + ///////////////////////////////////////////// + + for(int i=0;i> mpi_split[k]; + } + break; + } + } + + + double stp = 1.e-8; + int nrhs = 1; + int me; + for(int i=0;i seeds({1,2,3,4}); + + std::vector src(nrhs,FGrid); + std::vector src_chk(nrhs,FGrid); + std::vector result(nrhs,FGrid); + FermionField tmp(FGrid); + std::cout << GridLogMessage << "Made the Fermion Fields"< HermOp(Ddwf); + MdagMLinearOperator HermOpCk(Dchk); + ConjugateGradient CG((stp),100000); + s_res = zero; + + CG(HermOp,s_src,s_res); + + std::cout << GridLogMessage << " split residual norm "< iterations(nrhs,0); + iterations[me] = CG.IterationsToComplete; + + for(int n=0;nGlobalSum(iterations[n]); + std::cout << GridLogMessage<<" Rank "< BCGV (BlockCGrQVec,blockDim,stp,100000); + { + BCGV(HermOpCk,src,result); + } + + + + Grid_finalize(); +} diff --git a/tests/solver/Test_mobius_bcg_nosplit.cc b/tests/solver/Test_mobius_bcg_nosplit.cc new file mode 100644 index 00000000..f3ed621f --- /dev/null +++ b/tests/solver/Test_mobius_bcg_nosplit.cc @@ -0,0 +1,144 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_dwf_mrhs_cg.cc + + Copyright (C) 2015 + +Author: Peter Boyle + + 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 std; +using namespace Grid; +using namespace Grid::QCD; + +int main (int argc, char ** argv) +{ + typedef typename DomainWallFermionR::FermionField FermionField; + typedef typename DomainWallFermionR::ComplexField ComplexField; + typename DomainWallFermionR::ImplParams params; + + const int Ls=16; + + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + + std::vector boundary_phases(Nd,1.); + boundary_phases[Nd-1]=-1.; + params.boundary_phases = boundary_phases; + + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), + GridDefaultSimd(Nd,vComplex::Nsimd()), + GridDefaultMpi()); + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + GridRedBlackCartesian * rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + + double stp = 1.e-8; + int nrhs = 2; + + /////////////////////////////////////////////// + // Set up the problem as a 4d spreadout job + /////////////////////////////////////////////// + std::vector seeds({1,2,3,4}); + + std::vector src(nrhs,FGrid); + std::vector src_chk(nrhs,FGrid); + std::vector result(nrhs,FGrid); + FermionField tmp(FGrid); + std::cout << GridLogMessage << "Made the Fermion Fields"< HermOp(Ddwf); + ConjugateGradient CG((stp),100000); + + for(int rhs=0;rhs<1;rhs++){ + result[rhs] = zero; + CG(HermOp,src[rhs],result[rhs]); + } + + for(int rhs=0;rhs<1;rhs++){ + std::cout << " Result["< + +#include +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +int main (int argc, char ** argv) +{ + typedef typename DomainWallFermionR::FermionField FermionField; + typedef typename DomainWallFermionR::ComplexField ComplexField; + typename DomainWallFermionR::ImplParams params; + + const int Ls=16; + + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + + std::vector boundary_phases(Nd,1.); + boundary_phases[Nd-1]=-1.; + params.boundary_phases = boundary_phases; + + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), + GridDefaultSimd(Nd,vComplex::Nsimd()), + GridDefaultMpi()); + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + GridRedBlackCartesian * rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + + double stp = 1.e-8; + int nrhs = 2; + + /////////////////////////////////////////////// + // Set up the problem as a 4d spreadout job + /////////////////////////////////////////////// + std::vector seeds({1,2,3,4}); + + std::vector src4(nrhs,UGrid); + std::vector src(nrhs,FGrid); + std::vector src_chk(nrhs,FGrid); + std::vector result(nrhs,FGrid); + FermionField tmp(FGrid); + std::cout << GridLogMessage << "Made the Fermion Fields"< HermOp(Ddwf); + ConjugateGradient CG((stp),100000); + + for(int rhs=0;rhs<1;rhs++){ + result[rhs] = zero; + // CG(HermOp,src[rhs],result[rhs]); + } + + for(int rhs=0;rhs<1;rhs++){ + std::cout << " Result["< + +#include +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +int main (int argc, char ** argv) +{ + typedef typename DomainWallFermionR::FermionField FermionField; + typedef typename DomainWallFermionR::ComplexField ComplexField; + typename DomainWallFermionR::ImplParams params; + + const int Ls=16; + + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + + std::vector boundary_phases(Nd,1.); + boundary_phases[Nd-1]=-1.; + params.boundary_phases = boundary_phases; + + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), + GridDefaultSimd(Nd,vComplex::Nsimd()), + GridDefaultMpi()); + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + GridRedBlackCartesian * rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + + double stp = 1.e-8; + int nrhs = 4; + + /////////////////////////////////////////////// + // Set up the problem as a 4d spreadout job + /////////////////////////////////////////////// + std::vector seeds({1,2,3,4}); + + std::vector src(nrhs,FGrid); + std::vector src_chk(nrhs,FGrid); + std::vector result(nrhs,FGrid); + FermionField tmp(FGrid); + std::cout << GridLogMessage << "Made the Fermion Fields"< HermOp(Ddwf); + ConjugateGradient CG((stp),100000); + + for(int rhs=0;rhs<1;rhs++){ + result[rhs] = zero; + CG(HermOp,src[rhs],result[rhs]); + } + + for(int rhs=0;rhs<1;rhs++){ + std::cout << " Result["< src_v (Ls,UrbGrid); + std::vector result_v(Ls,UrbGrid); + for(int s=0;s void RedBlackSolution(Matrix & _Matrix,const Field &sol_o, const Field &src_e,Field &sol) + { + GridBase *grid = _Matrix.RedBlackGrid(); + GridBase *fgrid= _Matrix.Grid(); + + Field tmp(grid); + Field sol_e(grid); /////////////////////////////////////////////////// // sol_e = M_ee^-1 * ( src_e - Meo sol_o )... /////////////////////////////////////////////////// _Matrix.Meooe(sol_o,tmp); assert( tmp.checkerboard ==Even); - src_e = src_e-tmp; assert( src_e.checkerboard ==Even); - _Matrix.MooeeInv(src_e,sol_e); assert( sol_e.checkerboard ==Even); + tmp = src_e-tmp; assert( src_e.checkerboard ==Even); + _Matrix.MooeeInv(tmp,sol_e); assert( sol_e.checkerboard ==Even); - setCheckerboard(out,sol_e); assert( sol_e.checkerboard ==Even); - setCheckerboard(out,sol_o); assert( sol_o.checkerboard ==Odd ); + setCheckerboard(sol,sol_e); assert( sol_e.checkerboard ==Even); + setCheckerboard(sol,sol_o); assert( sol_o.checkerboard ==Odd ); + }; + + template + void operator()(Matrix &_Matrix, const std::vector &in, std::vector &out) + { + ZeroGuesser guess; + (*this)(_Matrix,in,out,guess); + } + + template + void operator()(Matrix &_Matrix, const std::vector &in, std::vector &out,Guesser &guess) + { + GridBase *grid = _Matrix.RedBlackGrid(); + GridBase *fgrid= _Matrix.Grid(); + int nblock = in.size(); + + std::vector src_o(nblock,grid); + std::vector sol_o(nblock,grid); + + std::vector guess_save; + + Field resid(fgrid); + Field tmp(grid); + + //////////////////////////////////////////////// + // Prepare RedBlack source + //////////////////////////////////////////////// + for(int b=0;b _HermOpEO(_Matrix); + _HermitianRBSolver(_HermOpEO,src_o,sol_o); + + //////////////////////////////////////////////// + // A2A boolean behavioural control & reconstruct other checkerboard + //////////////////////////////////////////////// + for(int b=0;b + void operator() (Matrix & _Matrix,const Field &in, Field &out,Guesser &guess){ + + // FIXME CGdiagonalMee not implemented virtual function + // FIXME use CBfactorise to control schur decomp + GridBase *grid = _Matrix.RedBlackGrid(); + GridBase *fgrid= _Matrix.Grid(); + + Field resid(fgrid); + Field src_o(grid); + Field src_e(grid); + Field sol_o(grid); + + + //////////////////////////////////////////////// + // RedBlack source + //////////////////////////////////////////////// + RedBlackSource(_Matrix,in,src_e,src_o); + + //////////////////////////////// + // Construct the guess + //////////////////////////////// + Field tmp(grid); + guess(src_o,sol_o); + + Field guess_save(grid); + guess_save = sol_o; + + std::cout< _HermOpEO(_Matrix); + _HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(tmp.checkerboard==Odd); + + //////////////////////////////////////////////// + // Fionn A2A boolean behavioural control + //////////////////////////////////////////////// + if (subGuess) tmp = sol_o-guess_save; + else tmp = sol_o; + + //////////////////////////////////////////////// + // MooeeInv + //////////////////////////////////////////////// + _Matrix.MooeeInv(tmp,sol_o); assert( sol_o.checkerboard ==Odd); + + /////////////////////////////////////////////////// + // RedBlack solution needs the even source + /////////////////////////////////////////////////// + RedBlackSolution(_Matrix,sol_o,src_e,out); // Verify the unprec residual if ( ! subGuess ) { diff --git a/tests/solver/Test_mobius_bcg_prec_nosplit.cc b/tests/solver/Test_mobius_bcg_prec_nosplit.cc index 347322c8..63078613 100644 --- a/tests/solver/Test_mobius_bcg_prec_nosplit.cc +++ b/tests/solver/Test_mobius_bcg_prec_nosplit.cc @@ -83,7 +83,7 @@ int main (int argc, char ** argv) LatticeGaugeField Umu(UGrid); - int conf = 0; + int conf = 2; if(conf==0) { FieldMetaData header; std::string file("./lat.in"); @@ -131,9 +131,12 @@ int main (int argc, char ** argv) for(int s=0;s BCGV (BlockCGrQVec,blockDim,stp,100000); + + { - BCGV(HermOp,src,result); + BlockConjugateGradient BCGV (BlockCGrQVec,blockDim,stp,100000); + SchurRedBlackDiagTwoSolve SchurSolver(BCGV); + SchurSolver(Ddwf,src,result); } for(int rhs=0;rhs Date: Tue, 6 Nov 2018 22:18:17 +0000 Subject: [PATCH 27/37] Patch to broken assertion --- Grid/algorithms/iterative/SchurRedBlack.h | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Grid/algorithms/iterative/SchurRedBlack.h b/Grid/algorithms/iterative/SchurRedBlack.h index cf6b5a13..b779cfcb 100644 --- a/Grid/algorithms/iterative/SchurRedBlack.h +++ b/Grid/algorithms/iterative/SchurRedBlack.h @@ -474,16 +474,17 @@ namespace Grid { Field tmp(grid); guess(src_o,sol_o); + // std::cout< _HermOpEO(_Matrix); - _HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(tmp.checkerboard==Odd); + _HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Odd); //////////////////////////////////////////////// // Fionn A2A boolean behavioural control From 6f421c7a6f42a8e41bb6c9aa4360618ca9391a18 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Wed, 7 Nov 2018 12:26:56 +0000 Subject: [PATCH 28/37] Block solver in the SchurRedBlack plus timing report cleaner --- Grid/algorithms/SparseMatrix.h | 8 + Grid/algorithms/iterative/SchurRedBlack.h | 570 ++++++++-------------- Grid/log/Log.h | 4 +- Grid/perfmon/Timer.h | 30 +- Grid/qcd/action/fermion/FermionOperator.h | 5 - 5 files changed, 243 insertions(+), 374 deletions(-) diff --git a/Grid/algorithms/SparseMatrix.h b/Grid/algorithms/SparseMatrix.h index 1611a6f4..9b0f7659 100644 --- a/Grid/algorithms/SparseMatrix.h +++ b/Grid/algorithms/SparseMatrix.h @@ -55,6 +55,14 @@ namespace Grid { template class CheckerBoardedSparseMatrixBase : public SparseMatrixBase { public: virtual GridBase *RedBlackGrid(void)=0; + + ////////////////////////////////////////////////////////////////////// + // Query the even even properties to make algorithmic decisions + ////////////////////////////////////////////////////////////////////// + virtual RealD Mass(void) { return 0.0; }; + virtual int ConstEE(void) { return 0; }; // Disable assumptions unless overridden + virtual int isTrivialEE(void) { return 0; }; // by a derived class that knows better + // half checkerboard operaions virtual void Meooe (const Field &in, Field &out)=0; virtual void Mooee (const Field &in, Field &out)=0; diff --git a/Grid/algorithms/iterative/SchurRedBlack.h b/Grid/algorithms/iterative/SchurRedBlack.h index b779cfcb..fdb17a98 100644 --- a/Grid/algorithms/iterative/SchurRedBlack.h +++ b/Grid/algorithms/iterative/SchurRedBlack.h @@ -86,229 +86,23 @@ Author: Peter Boyle */ namespace Grid { + /////////////////////////////////////////////////////////////////////////////////////////////////////// + // Use base class to share code + /////////////////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////////////////// // Take a matrix and form a Red Black solver calling a Herm solver // Use of RB info prevents making SchurRedBlackSolve conform to standard interface /////////////////////////////////////////////////////////////////////////////////////////////////////// - // Now make the norm reflect extra factor of Mee - template class SchurRedBlackStaggeredSolve { - private: + template class SchurRedBlackBase { + protected: + typedef CheckerBoardedSparseMatrixBase Matrix; OperatorFunction & _HermitianRBSolver; int CBfactorise; bool subGuess; public: - ///////////////////////////////////////////////////// - // Wrap the usual normal equations Schur trick - ///////////////////////////////////////////////////// - SchurRedBlackStaggeredSolve(OperatorFunction &HermitianRBSolver, const bool initSubGuess = false) : - _HermitianRBSolver(HermitianRBSolver) - { - CBfactorise=0; - subtractGuess(initSubGuess); - }; - void subtractGuess(const bool initSubGuess) - { - subGuess = initSubGuess; - } - bool isSubtractGuess(void) - { - return subGuess; - } - - template - void operator() (Matrix & _Matrix,const Field &in, Field &out){ - ZeroGuesser guess; - (*this)(_Matrix,in,out,guess); - } - template - void operator() (Matrix & _Matrix,const Field &in, Field &out, Guesser &guess){ - - // FIXME CGdiagonalMee not implemented virtual function - // FIXME use CBfactorise to control schur decomp - GridBase *grid = _Matrix.RedBlackGrid(); - GridBase *fgrid= _Matrix.Grid(); - - SchurStaggeredOperator _HermOpEO(_Matrix); - - Field src_e(grid); - Field src_o(grid); - Field sol_e(grid); - Field sol_o(grid); - Field tmp(grid); - Field Mtmp(grid); - Field resid(fgrid); - - std::cout << GridLogMessage << " SchurRedBlackStaggeredSolve " < using SchurRedBlackStagSolve = SchurRedBlackStaggeredSolve; - - /////////////////////////////////////////////////////////////////////////////////////////////////////// - // Take a matrix and form a Red Black solver calling a Herm solver - // Use of RB info prevents making SchurRedBlackSolve conform to standard interface - /////////////////////////////////////////////////////////////////////////////////////////////////////// - template class SchurRedBlackDiagMooeeSolve { - private: - OperatorFunction & _HermitianRBSolver; - int CBfactorise; - bool subGuess; - public: - - ///////////////////////////////////////////////////// - // Wrap the usual normal equations Schur trick - ///////////////////////////////////////////////////// - SchurRedBlackDiagMooeeSolve(OperatorFunction &HermitianRBSolver,int cb=0, const bool initSubGuess = false) : _HermitianRBSolver(HermitianRBSolver) - { - CBfactorise=cb; - subtractGuess(initSubGuess); - }; - void subtractGuess(const bool initSubGuess) - { - subGuess = initSubGuess; - } - bool isSubtractGuess(void) - { - return subGuess; - } - template - void operator() (Matrix & _Matrix,const Field &in, Field &out){ - ZeroGuesser guess; - (*this)(_Matrix,in,out,guess); - } - template - void operator() (Matrix & _Matrix,const Field &in, Field &out,Guesser &guess){ - - // FIXME CGdiagonalMee not implemented virtual function - // FIXME use CBfactorise to control schur decomp - GridBase *grid = _Matrix.RedBlackGrid(); - GridBase *fgrid= _Matrix.Grid(); - - SchurDiagMooeeOperator _HermOpEO(_Matrix); - - Field src_e(grid); - Field src_o(grid); - Field sol_e(grid); - Field sol_o(grid); - Field tmp(grid); - Field Mtmp(grid); - Field resid(fgrid); - - pickCheckerboard(Even,src_e,in); - pickCheckerboard(Odd ,src_o,in); - pickCheckerboard(Even,sol_e,out); - pickCheckerboard(Odd ,sol_o,out); - - ///////////////////////////////////////////////////// - // src_o = Mdag * (source_o - Moe MeeInv source_e) - ///////////////////////////////////////////////////// - _Matrix.MooeeInv(src_e,tmp); assert( tmp.checkerboard ==Even); - _Matrix.Meooe (tmp,Mtmp); assert( Mtmp.checkerboard ==Odd); - tmp=src_o-Mtmp; assert( tmp.checkerboard ==Odd); - - // get the right MpcDag - _HermOpEO.MpcDag(tmp,src_o); assert(src_o.checkerboard ==Odd); - - ////////////////////////////////////////////////////////////// - // Call the red-black solver - ////////////////////////////////////////////////////////////// - std::cout< class SchurRedBlackDiagTwoSolve { - private: - OperatorFunction & _HermitianRBSolver; - int CBfactorise; - bool subGuess; - public: - - ///////////////////////////////////////////////////// - // Wrap the usual normal equations Schur trick - ///////////////////////////////////////////////////// - SchurRedBlackDiagTwoSolve(OperatorFunction &HermitianRBSolver, const bool initSubGuess = false) : - _HermitianRBSolver(HermitianRBSolver) + SchurRedBlackBase(OperatorFunction &HermitianRBSolver, const bool initSubGuess = false) : + _HermitianRBSolver(HermitianRBSolver) { CBfactorise = 0; subtractGuess(initSubGuess); @@ -322,63 +116,20 @@ namespace Grid { return subGuess; } - template + ///////////////////////////////////////////////////////////// + // Shared code + ///////////////////////////////////////////////////////////// void operator() (Matrix & _Matrix,const Field &in, Field &out){ ZeroGuesser guess; (*this)(_Matrix,in,out,guess); } - - template void RedBlackSource(Matrix & _Matrix,const Field &src, Field &src_e,Field &src_o) - { - GridBase *grid = _Matrix.RedBlackGrid(); - GridBase *fgrid= _Matrix.Grid(); - - SchurDiagTwoOperator _HermOpEO(_Matrix); - - Field tmp(grid); - Field Mtmp(grid); - - pickCheckerboard(Even,src_e,src); - pickCheckerboard(Odd ,src_o,src); - - ///////////////////////////////////////////////////// - // src_o = Mdag * (source_o - Moe MeeInv source_e) - ///////////////////////////////////////////////////// - _Matrix.MooeeInv(src_e,tmp); assert( tmp.checkerboard ==Even); - _Matrix.Meooe (tmp,Mtmp); assert( Mtmp.checkerboard ==Odd); - tmp=src_o-Mtmp; assert( tmp.checkerboard ==Odd); - - // get the right MpcDag - _HermOpEO.MpcDag(tmp,src_o); assert(src_o.checkerboard ==Odd); - } - - template void RedBlackSolution(Matrix & _Matrix,const Field &sol_o, const Field &src_e,Field &sol) - { - GridBase *grid = _Matrix.RedBlackGrid(); - GridBase *fgrid= _Matrix.Grid(); - - Field tmp(grid); - Field sol_e(grid); - - /////////////////////////////////////////////////// - // sol_e = M_ee^-1 * ( src_e - Meo sol_o )... - /////////////////////////////////////////////////// - _Matrix.Meooe(sol_o,tmp); assert( tmp.checkerboard ==Even); - tmp = src_e-tmp; assert( src_e.checkerboard ==Even); - _Matrix.MooeeInv(tmp,sol_e); assert( sol_e.checkerboard ==Even); - - setCheckerboard(sol,sol_e); assert( sol_e.checkerboard ==Even); - setCheckerboard(sol,sol_o); assert( sol_o.checkerboard ==Odd ); - }; - - template void operator()(Matrix &_Matrix, const std::vector &in, std::vector &out) { ZeroGuesser guess; (*this)(_Matrix,in,out,guess); } - template + template void operator()(Matrix &_Matrix, const std::vector &in, std::vector &out,Guesser &guess) { GridBase *grid = _Matrix.RedBlackGrid(); @@ -414,19 +165,15 @@ namespace Grid { ////////////////////////////////////////////////////////////// // Call the block solver ////////////////////////////////////////////////////////////// - std::cout< _HermOpEO(_Matrix); - _HermitianRBSolver(_HermOpEO,src_o,sol_o); + std::cout< + template void operator() (Matrix & _Matrix,const Field &in, Field &out,Guesser &guess){ // FIXME CGdiagonalMee not implemented virtual function @@ -462,7 +208,6 @@ namespace Grid { Field src_e(grid); Field sol_o(grid); - //////////////////////////////////////////////// // RedBlack source //////////////////////////////////////////////// @@ -474,28 +219,18 @@ namespace Grid { Field tmp(grid); guess(src_o,sol_o); - // std::cout< _HermOpEO(_Matrix); - _HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Odd); + RedBlackSolve(_Matrix,src_o,sol_o); //////////////////////////////////////////////// // Fionn A2A boolean behavioural control //////////////////////////////////////////////// - if (subGuess) tmp = sol_o-guess_save; - else tmp = sol_o; - - //////////////////////////////////////////////// - // MooeeInv - //////////////////////////////////////////////// - _Matrix.MooeeInv(tmp,sol_o); assert( sol_o.checkerboard ==Odd); + if (subGuess) sol_o= sol_o-guess_save; /////////////////////////////////////////////////// // RedBlack solution needs the even source @@ -509,68 +244,182 @@ namespace Grid { RealD ns = norm2(in); RealD nr = norm2(resid); - std::cout< &src_o, std::vector &sol_o)=0; + }; - /////////////////////////////////////////////////////////////////////////////////////////////////////// - // Take a matrix and form a Red Black solver calling a Herm solver - // Use of RB info prevents making SchurRedBlackSolve conform to standard interface - /////////////////////////////////////////////////////////////////////////////////////////////////////// - template class SchurRedBlackDiagTwoMixed { - private: - LinearFunction & _HermitianRBSolver; - int CBfactorise; - bool subGuess; + + template class SchurRedBlackStaggeredSolve : public SchurRedBlackBase { public: + typedef CheckerBoardedSparseMatrixBase Matrix; + + SchurRedBlackStaggeredSolve(OperatorFunction &HermitianRBSolver, const bool initSubGuess = false) + : SchurRedBlackBase (HermitianRBSolver,initSubGuess) + { + } + + ////////////////////////////////////////////////////// + // Override RedBlack specialisation + ////////////////////////////////////////////////////// + virtual void RedBlackSource(Matrix & _Matrix,const Field &src, Field &src_e,Field &src_o) + { + GridBase *grid = _Matrix.RedBlackGrid(); + GridBase *fgrid= _Matrix.Grid(); + + Field tmp(grid); + Field Mtmp(grid); + + pickCheckerboard(Even,src_e,src); + pickCheckerboard(Odd ,src_o,src); + + ///////////////////////////////////////////////////// + // src_o = (source_o - Moe MeeInv source_e) + ///////////////////////////////////////////////////// + _Matrix.MooeeInv(src_e,tmp); assert( tmp.checkerboard ==Even); + _Matrix.Meooe (tmp,Mtmp); assert( Mtmp.checkerboard ==Odd); + tmp=src_o-Mtmp; assert( tmp.checkerboard ==Odd); + + _Matrix.Mooee(tmp,src_o); // Extra factor of "m" in source from dumb choice of matrix norm. + } + virtual void RedBlackSolution(Matrix & _Matrix,const Field &sol_o, const Field &src_e_c,Field &sol) + { + GridBase *grid = _Matrix.RedBlackGrid(); + GridBase *fgrid= _Matrix.Grid(); + + Field tmp(grid); + Field sol_e(grid); + Field src_e(grid); + + src_e = src_e_c; // Const correctness + + /////////////////////////////////////////////////// + // sol_e = M_ee^-1 * ( src_e - Meo sol_o )... + /////////////////////////////////////////////////// + _Matrix.Meooe(sol_o,tmp); assert( tmp.checkerboard ==Even); + src_e = src_e-tmp; assert( src_e.checkerboard ==Even); + _Matrix.MooeeInv(src_e,sol_e); assert( sol_e.checkerboard ==Even); + + setCheckerboard(sol,sol_e); assert( sol_e.checkerboard ==Even); + setCheckerboard(sol,sol_o); assert( sol_o.checkerboard ==Odd ); + } + virtual void RedBlackSolve (Matrix & _Matrix,const Field &src_o, Field &sol_o) + { + SchurStaggeredOperator _HermOpEO(_Matrix); + this->_HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Odd); + }; + virtual void RedBlackSolve (Matrix & _Matrix,const std::vector &src_o, std::vector &sol_o) + { + SchurStaggeredOperator _HermOpEO(_Matrix); + this->_HermitianRBSolver(_HermOpEO,src_o,sol_o); + } + }; + template using SchurRedBlackStagSolve = SchurRedBlackStaggeredSolve; + + /////////////////////////////////////////////////////////////////////////////////////////////////////// + // Site diagonal has Mooee on it. + /////////////////////////////////////////////////////////////////////////////////////////////////////// + template class SchurRedBlackDiagMooeeSolve : public SchurRedBlackBase { + public: + typedef CheckerBoardedSparseMatrixBase Matrix; + + SchurRedBlackDiagMooeeSolve(OperatorFunction &HermitianRBSolver, const bool initSubGuess = false) + : SchurRedBlackBase (HermitianRBSolver,initSubGuess) {}; + + + ////////////////////////////////////////////////////// + // Override RedBlack specialisation + ////////////////////////////////////////////////////// + virtual void RedBlackSource(Matrix & _Matrix,const Field &src, Field &src_e,Field &src_o) + { + GridBase *grid = _Matrix.RedBlackGrid(); + GridBase *fgrid= _Matrix.Grid(); + + Field tmp(grid); + Field Mtmp(grid); + + pickCheckerboard(Even,src_e,src); + pickCheckerboard(Odd ,src_o,src); + + ///////////////////////////////////////////////////// + // src_o = Mdag * (source_o - Moe MeeInv source_e) + ///////////////////////////////////////////////////// + _Matrix.MooeeInv(src_e,tmp); assert( tmp.checkerboard ==Even); + _Matrix.Meooe (tmp,Mtmp); assert( Mtmp.checkerboard ==Odd); + tmp=src_o-Mtmp; assert( tmp.checkerboard ==Odd); + + // get the right MpcDag + SchurDiagMooeeOperator _HermOpEO(_Matrix); + _HermOpEO.MpcDag(tmp,src_o); assert(src_o.checkerboard ==Odd); + + } + virtual void RedBlackSolution(Matrix & _Matrix,const Field &sol_o, const Field &src_e,Field &sol) + { + GridBase *grid = _Matrix.RedBlackGrid(); + GridBase *fgrid= _Matrix.Grid(); + + Field tmp(grid); + Field sol_e(grid); + Field src_e_i(grid); + /////////////////////////////////////////////////// + // sol_e = M_ee^-1 * ( src_e - Meo sol_o )... + /////////////////////////////////////////////////// + _Matrix.Meooe(sol_o,tmp); assert( tmp.checkerboard ==Even); + src_e_i = src_e-tmp; assert( src_e_i.checkerboard ==Even); + _Matrix.MooeeInv(src_e_i,sol_e); assert( sol_e.checkerboard ==Even); + + setCheckerboard(sol,sol_e); assert( sol_e.checkerboard ==Even); + setCheckerboard(sol,sol_o); assert( sol_o.checkerboard ==Odd ); + } + virtual void RedBlackSolve (Matrix & _Matrix,const Field &src_o, Field &sol_o) + { + SchurDiagMooeeOperator _HermOpEO(_Matrix); + this->_HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Odd); + }; + virtual void RedBlackSolve (Matrix & _Matrix,const std::vector &src_o, std::vector &sol_o) + { + SchurDiagMooeeOperator _HermOpEO(_Matrix); + this->_HermitianRBSolver(_HermOpEO,src_o,sol_o); + } + }; + + /////////////////////////////////////////////////////////////////////////////////////////////////////// + // Site diagonal is identity, right preconditioned by Mee^inv + // ( 1 - Meo Moo^inv Moe Mee^inv ) phi =( 1 - Meo Moo^inv Moe Mee^inv ) Mee psi = = eta = eta + //=> psi = MeeInv phi + /////////////////////////////////////////////////////////////////////////////////////////////////////// + template class SchurRedBlackDiagTwoSolve : public SchurRedBlackBase { + public: + typedef CheckerBoardedSparseMatrixBase Matrix; ///////////////////////////////////////////////////// // Wrap the usual normal equations Schur trick ///////////////////////////////////////////////////// - SchurRedBlackDiagTwoMixed(LinearFunction &HermitianRBSolver, const bool initSubGuess = false) : - _HermitianRBSolver(HermitianRBSolver) - { - CBfactorise=0; - subtractGuess(initSubGuess); - }; - void subtractGuess(const bool initSubGuess) - { - subGuess = initSubGuess; - } - bool isSubtractGuess(void) - { - return subGuess; - } + SchurRedBlackDiagTwoSolve(OperatorFunction &HermitianRBSolver, const bool initSubGuess = false) + : SchurRedBlackBase(HermitianRBSolver,initSubGuess) {}; - template - void operator() (Matrix & _Matrix,const Field &in, Field &out){ - ZeroGuesser guess; - (*this)(_Matrix,in,out,guess); - } - template - void operator() (Matrix & _Matrix,const Field &in, Field &out,Guesser &guess){ - - // FIXME CGdiagonalMee not implemented virtual function - // FIXME use CBfactorise to control schur decomp + virtual void RedBlackSource(Matrix & _Matrix,const Field &src, Field &src_e,Field &src_o) + { GridBase *grid = _Matrix.RedBlackGrid(); GridBase *fgrid= _Matrix.Grid(); SchurDiagTwoOperator _HermOpEO(_Matrix); - - Field src_e(grid); - Field src_o(grid); - Field sol_e(grid); - Field sol_o(grid); + Field tmp(grid); Field Mtmp(grid); - Field resid(fgrid); - pickCheckerboard(Even,src_e,in); - pickCheckerboard(Odd ,src_o,in); - pickCheckerboard(Even,sol_e,out); - pickCheckerboard(Odd ,sol_o,out); + pickCheckerboard(Even,src_e,src); + pickCheckerboard(Odd ,src_o,src); ///////////////////////////////////////////////////// // src_o = Mdag * (source_o - Moe MeeInv source_e) @@ -581,43 +430,44 @@ namespace Grid { // get the right MpcDag _HermOpEO.MpcDag(tmp,src_o); assert(src_o.checkerboard ==Odd); + } - ////////////////////////////////////////////////////////////// - // Call the red-black solver - ////////////////////////////////////////////////////////////// - std::cout< _HermOpEO(_Matrix); + this->_HermitianRBSolver(_HermOpEO,src_o,sol_o); + }; + virtual void RedBlackSolve (Matrix & _Matrix,const std::vector &src_o, std::vector &sol_o) + { + SchurDiagTwoOperator _HermOpEO(_Matrix); + this->_HermitianRBSolver(_HermOpEO,src_o,sol_o); + } }; - } #endif diff --git a/Grid/log/Log.h b/Grid/log/Log.h index b58c5d16..5d97ee5a 100644 --- a/Grid/log/Log.h +++ b/Grid/log/Log.h @@ -146,9 +146,11 @@ public: if ( log.timestamp ) { log.StopWatch->Stop(); GridTime now = log.StopWatch->Elapsed(); + if ( log.timing_mode==1 ) log.StopWatch->Reset(); log.StopWatch->Start(); - stream << log.evidence()<< std::setw(6)< GridTimePoint; -typedef std::chrono::milliseconds GridMillisecs; -typedef std::chrono::microseconds GridTime; -typedef std::chrono::microseconds GridUsecs; -inline std::ostream& operator<< (std::ostream & stream, const std::chrono::milliseconds & time) +typedef std::chrono::seconds GridSecs; +typedef std::chrono::milliseconds GridMillisecs; +typedef std::chrono::microseconds GridUsecs; +typedef std::chrono::microseconds GridTime; + +inline std::ostream& operator<< (std::ostream & stream, const GridSecs & time) { - stream << time.count()<<" ms"; + stream << time.count()<<" s"; return stream; } -inline std::ostream& operator<< (std::ostream & stream, const std::chrono::microseconds & time) +inline std::ostream& operator<< (std::ostream & stream, const GridMillisecs & now) { - stream << time.count()<<" usec"; + GridSecs second(1); + auto secs = now/second ; + auto subseconds = now%second ; + stream << secs<<"."< Date: Wed, 7 Nov 2018 13:32:46 +0000 Subject: [PATCH 29/37] Add CXX to grid-config --- configure.ac | 2 ++ grid-config.in | 4 ++++ 2 files changed, 6 insertions(+) diff --git a/configure.ac b/configure.ac index d362ed5a..ad2b8bba 100644 --- a/configure.ac +++ b/configure.ac @@ -485,6 +485,7 @@ DX_INIT_DOXYGEN([$PACKAGE_NAME], [doxygen.cfg]) ############### Ouput cwd=`pwd -P`; cd ${srcdir}; abs_srcdir=`pwd -P`; cd ${cwd} +GRID_CXX="$CXX" GRID_CXXFLAGS="$AM_CXXFLAGS $CXXFLAGS" GRID_LDFLAGS="$AM_LDFLAGS $LDFLAGS" GRID_LIBS=$LIBS @@ -497,6 +498,7 @@ AM_LDFLAGS="-L${cwd}/Grid $AM_LDFLAGS" AC_SUBST([AM_CFLAGS]) AC_SUBST([AM_CXXFLAGS]) AC_SUBST([AM_LDFLAGS]) +AC_SUBST([GRID_CXX]) AC_SUBST([GRID_CXXFLAGS]) AC_SUBST([GRID_LDFLAGS]) AC_SUBST([GRID_LIBS]) diff --git a/grid-config.in b/grid-config.in index bd340846..f39b01bd 100755 --- a/grid-config.in +++ b/grid-config.in @@ -61,6 +61,10 @@ while test $# -gt 0; do echo @GRID_CXXFLAGS@ ;; + --cxx) + echo @GRID_CXX@ + ;; + --ldflags) echo @GRID_LDFLAGS@ ;; From e0a79a5bbfe11ee5ea327a4cf1f30d4a283384dc Mon Sep 17 00:00:00 2001 From: Vera Guelpers Date: Wed, 7 Nov 2018 15:01:22 +0000 Subject: [PATCH 30/37] Hadrons: PR#177: Electrify gauge: Single Precision fix --- Hadrons/Modules/MGauge/Electrify.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Hadrons/Modules/MGauge/Electrify.hpp b/Hadrons/Modules/MGauge/Electrify.hpp index 533ac704..c5e6735d 100644 --- a/Hadrons/Modules/MGauge/Electrify.hpp +++ b/Hadrons/Modules/MGauge/Electrify.hpp @@ -139,7 +139,7 @@ void TElectrify::execute(void) for(unsigned int mu = 0; mu < env().getNd(); mu++) { - eiAmu = exp(i * par().e * par().charge * PeekIndex(A, mu)); + eiAmu = exp(i * (Real)(par().e * par().charge) * PeekIndex(A, mu)); PokeIndex(Ue, PeekIndex(U, mu) * eiAmu, mu); } } From 839605c45c4da742f074f807d481d3d05761a27c Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Wed, 7 Nov 2018 23:38:46 +0000 Subject: [PATCH 31/37] Verbose reduce --- Grid/communicator/Communicator_mpi3.cc | 6 ------ Grid/communicator/SharedMemoryMPI.cc | 6 +++--- 2 files changed, 3 insertions(+), 9 deletions(-) diff --git a/Grid/communicator/Communicator_mpi3.cc b/Grid/communicator/Communicator_mpi3.cc index 424b7973..e6800afe 100644 --- a/Grid/communicator/Communicator_mpi3.cc +++ b/Grid/communicator/Communicator_mpi3.cc @@ -50,8 +50,6 @@ void CartesianCommunicator::Init(int *argc, char ***argv) assert(0); } - Grid_quiesce_nodes(); - // Never clean up as done once. MPI_Comm_dup (MPI_COMM_WORLD,&communicator_world); @@ -124,10 +122,8 @@ CartesianCommunicator::CartesianCommunicator(const std::vector &processors, // split the communicator ////////////////////////////////////////////////////////////////////////////////////////////////////// // int Nparent = parent._processors ; - // std::cout << " splitting from communicator "< &processors, int Nchild = Nparent/childsize; assert (childsize * Nchild == Nparent); - // std::cout << " child size "< ccoor(_ndimension); // coor within subcommunicator std::vector scoor(_ndimension); // coor of split within parent std::vector ssize(_ndimension); // coor of split within parent diff --git a/Grid/communicator/SharedMemoryMPI.cc b/Grid/communicator/SharedMemoryMPI.cc index b425d97b..d3b094ad 100644 --- a/Grid/communicator/SharedMemoryMPI.cc +++ b/Grid/communicator/SharedMemoryMPI.cc @@ -413,7 +413,7 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags) assert(((uint64_t)ptr&0x3F)==0); close(fd); WorldShmCommBufs[r] =ptr; - std::cout << "Set WorldShmCommBufs["< Date: Wed, 7 Nov 2018 23:39:18 +0000 Subject: [PATCH 32/37] Allow shrinking machine in orthog direction for extract slice local --- Grid/lattice/Lattice_transfer.h | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/Grid/lattice/Lattice_transfer.h b/Grid/lattice/Lattice_transfer.h index f988f310..4eca50d6 100644 --- a/Grid/lattice/Lattice_transfer.h +++ b/Grid/lattice/Lattice_transfer.h @@ -464,8 +464,10 @@ void InsertSliceLocal(const Lattice &lowDim, Lattice & higherDim,int assert(orthog>=0); for(int d=0;d_processors[d] == hg->_processors[d]); - assert(lg->_ldimensions[d] == hg->_ldimensions[d]); + if ( d!=orthog ) { + assert(lg->_processors[d] == hg->_processors[d]); + assert(lg->_ldimensions[d] == hg->_ldimensions[d]); + } } // the above should guarantee that the operations are local @@ -499,8 +501,10 @@ void ExtractSliceLocal(Lattice &lowDim, Lattice & higherDim,int slic assert(orthog>=0); for(int d=0;d_processors[d] == hg->_processors[d]); - assert(lg->_ldimensions[d] == hg->_ldimensions[d]); + if ( d!=orthog ) { + assert(lg->_processors[d] == hg->_processors[d]); + assert(lg->_ldimensions[d] == hg->_ldimensions[d]); + } } // the above should guarantee that the operations are local From 68c13045d6ce9982f9bc760c085e29f2c3eb8f2e Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Wed, 7 Nov 2018 23:40:15 +0000 Subject: [PATCH 33/37] Added a test for Felix and Michael to look at --- tests/debug/Test_split_laplacian.cc | 104 ++++++++++++++++++++++++++++ 1 file changed, 104 insertions(+) create mode 100644 tests/debug/Test_split_laplacian.cc diff --git a/tests/debug/Test_split_laplacian.cc b/tests/debug/Test_split_laplacian.cc new file mode 100644 index 00000000..174dc1b1 --- /dev/null +++ b/tests/debug/Test_split_laplacian.cc @@ -0,0 +1,104 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_dwf_mrhs_cg.cc + + Copyright (C) 2015 + +Author: Peter Boyle + + 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 std; +using namespace Grid; +using namespace Grid::QCD; + + +int main (int argc, char ** argv) +{ + typedef LatticeComplex ComplexField; + + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + int nd = latt_size.size(); + int ndm1 = nd-1; + + std::vector simd_layout = GridDefaultSimd(nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + std::vector mpi_split (mpi_layout.size(),1); + + std::cout << " Full " << GridCmdVectorIntToString(latt_size) << " subgrid" < latt_m = latt_size; latt_m[nd-1] = 1; + std::vector mpi_m = mpi_layout; mpi_m [nd-1] = 1; + std::vector simd_m = GridDefaultSimd(ndm1,vComplex::Nsimd()); simd_m.push_back(1); + + + std::cout << " Requesting " << GridCmdVectorIntToString(latt_m)<< " subgrid" <GlobalSum(tmp); + std::cout << GridLogMessage<< " Full nodes "<< tmp <GlobalSum(tmp); + std::cout << GridLogMessage<< " Split nodes "<< tmp <Barrier(); + + auto local_latt = GridN->LocalDimensions(); + + Full_cpy = zero; + std::vector seeds({1,2,3,4}); + GridParallelRNG RNG(GridN); RNG.SeedFixedIntegers(seeds); + + random(RNG,Full); + for(int t=0;t Date: Thu, 8 Nov 2018 08:58:09 +0000 Subject: [PATCH 34/37] Hadrons: sequential Aslash insertion and propagator on A2A vector --- Hadrons/Modules.hpp | 1 + Hadrons/Modules/MSolver/A2AAslashVector.cc | 35 ++++ Hadrons/Modules/MSolver/A2AAslashVector.hpp | 189 ++++++++++++++++++++ Hadrons/modules.inc | 2 + 4 files changed, 227 insertions(+) create mode 100644 Hadrons/Modules/MSolver/A2AAslashVector.cc create mode 100644 Hadrons/Modules/MSolver/A2AAslashVector.hpp diff --git a/Hadrons/Modules.hpp b/Hadrons/Modules.hpp index da57ab2b..27a889bf 100644 --- a/Hadrons/Modules.hpp +++ b/Hadrons/Modules.hpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include #include diff --git a/Hadrons/Modules/MSolver/A2AAslashVector.cc b/Hadrons/Modules/MSolver/A2AAslashVector.cc new file mode 100644 index 00000000..27344397 --- /dev/null +++ b/Hadrons/Modules/MSolver/A2AAslashVector.cc @@ -0,0 +1,35 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MSolver/A2AAslashVector.cc + +Copyright (C) 2015-2018 + +Author: Vera Guelpers + +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 MSolver; + +template class Grid::Hadrons::MSolver::TA2AAslashVector; +template class Grid::Hadrons::MSolver::TA2AAslashVector; diff --git a/Hadrons/Modules/MSolver/A2AAslashVector.hpp b/Hadrons/Modules/MSolver/A2AAslashVector.hpp new file mode 100644 index 00000000..c99a6f6b --- /dev/null +++ b/Hadrons/Modules/MSolver/A2AAslashVector.hpp @@ -0,0 +1,189 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MSolver/A2AAslashVector.hpp + +Copyright (C) 2015-2018 + +Author: Vera Guelpers + +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_MSolver_A2AAslashVector_hpp_ +#define Hadrons_MSolver_A2AAslashVector_hpp_ + +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Create all-to-all V & W vectors * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MSolver) + +/**************************************************************************** +* Calculate a sequential propagator on an insertion of i*g_mu*A_mu +* on an A2A vector +* +* vv_i(y) = S(y,x) * i * g_mu*A_mu(x) * v_i(x) +* +* with +* +* - v: A2A vector v_i(x) +* - photon: A_mu(x): electromagnetic photon field +* - solver: the solver for calculating the sequential propagator +* +*****************************************************************************/ + + +class A2AAslashVectorPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(A2AAslashVectorPar, + std::string, v, + std::string, photon, + std::string, solver); +}; + +template +class TA2AAslashVector : public Module +{ +public: + FERM_TYPE_ALIASES(FImpl,); + SOLVER_TYPE_ALIASES(FImpl,); +public: + typedef PhotonR::GaugeField EmField; +public: + // constructor + TA2AAslashVector(const std::string name); + // destructor + virtual ~TA2AAslashVector(void) {}; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +private: + unsigned int Ls_; +}; + +MODULE_REGISTER_TMP(A2AAslashVector,TA2AAslashVector, MSolver); +MODULE_REGISTER_TMP(ZA2AAslashVector,TA2AAslashVector, MSolver); + +/****************************************************************************** + * TA2AAslashVector implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TA2AAslashVector::TA2AAslashVector(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TA2AAslashVector::getInput(void) +{ + std::vector in = {par().v, par().photon, par().solver}; + + return in; +} + +template +std::vector TA2AAslashVector::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TA2AAslashVector::setup(void) +{ + Ls_ = env().getObjectLs(par().solver); + auto &vvector = envGet(std::vector, par().v); + unsigned int Nmodes = vvector.size(); + envCreate(std::vector, getName(), 1, + Nmodes, envGetGrid(FermionField)); + + envTmpLat(FermionField, "v4dtmp"); + envTmpLat(FermionField, "v5dtmp", Ls_); + envTmpLat(FermionField, "v5dtmp_sol", Ls_); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TA2AAslashVector::execute(void) +{ + auto &solver = envGet(Solver, par().solver); + auto &stoch_photon = envGet(EmField, par().photon); + auto &vvector = envGet(std::vector, par().v); + auto &Aslashv = envGet(std::vector, getName()); + unsigned int Nmodes = vvector.size(); + auto &mat = solver.getFMat(); + envGetTmp(FermionField, v4dtmp); + envGetTmp(FermionField, v5dtmp); + envGetTmp(FermionField, v5dtmp_sol); + + Complex ci(0.0,1.0); + + + startTimer("Seq Aslash"); + + LOG(Message) << "Calculate Sequential propagator on Aslash * v with the A2A vector " << par().v + << " and the photon field " << par().photon << std::endl; + + + for(unsigned int i=0; i(stoch_photon, mu) * (gmu * vvector[i]); + } + stopTimer("Multiply Aslash"); + + if (Ls_ == 1) + { + solver(Aslashv[i], v4dtmp); + } + else + { + mat.ImportPhysicalFermionSource(v4dtmp, v5dtmp); + solver(v5dtmp_sol, v5dtmp); + mat.ExportPhysicalFermionSolution(v5dtmp_sol, v4dtmp); + Aslashv[i] = v4dtmp; + } + } + + stopTimer("Seq Aslash"); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MSolver_A2AAslashVector_hpp_ diff --git a/Hadrons/modules.inc b/Hadrons/modules.inc index 157fe1a2..4ee09bb6 100644 --- a/Hadrons/modules.inc +++ b/Hadrons/modules.inc @@ -20,6 +20,7 @@ modules_cc =\ Modules/MSink/Point.cc \ Modules/MSink/Smear.cc \ Modules/MSolver/A2AVectors.cc \ + Modules/MSolver/A2AAslashVector.cc \ Modules/MSolver/RBPrecCG.cc \ Modules/MSolver/MixedPrecisionRBPrecCG.cc \ Modules/MSolver/LocalCoherenceLanczos.cc \ @@ -95,6 +96,7 @@ modules_hpp =\ Modules/MSolver/Guesser.hpp \ Modules/MSolver/RBPrecCG.hpp \ Modules/MSolver/A2AVectors.hpp \ + Modules/MSolver/A2AAslashVector.hpp \ Modules/MGauge/UnitEm.hpp \ Modules/MGauge/StoutSmearing.hpp \ Modules/MGauge/Unit.hpp \ From f3f24b3017eea98336c1c88783d51c7fd9c767c3 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 8 Nov 2018 12:55:25 +0000 Subject: [PATCH 35/37] Optional Twisted BC's added, in "DoubleStore" for WilsonImpl. Untested but doesn't affect answers when twists are all zero. The zero is the default behaviour for ImplParams. --- Grid/qcd/action/ActionParams.h | 7 +++++-- Grid/qcd/action/fermion/FermionOperatorImpl.h | 20 ++++++++++++++++--- 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/Grid/qcd/action/ActionParams.h b/Grid/qcd/action/ActionParams.h index d25b60a9..76868abc 100644 --- a/Grid/qcd/action/ActionParams.h +++ b/Grid/qcd/action/ActionParams.h @@ -44,12 +44,15 @@ namespace QCD { struct WilsonImplParams { bool overlapCommsCompute; + std::vector twist_n_2pi_L; std::vector boundary_phases; WilsonImplParams() : overlapCommsCompute(false) { boundary_phases.resize(Nd, 1.0); + twist_n_2pi_L.resize(Nd, 0.0); }; - WilsonImplParams(const std::vector phi) - : boundary_phases(phi), overlapCommsCompute(false) {} + WilsonImplParams(const std::vector phi) : boundary_phases(phi), overlapCommsCompute(false) { + twist_n_2pi_L.resize(Nd, 0.0); + } }; struct StaggeredImplParams { diff --git a/Grid/qcd/action/fermion/FermionOperatorImpl.h b/Grid/qcd/action/fermion/FermionOperatorImpl.h index b0ffa90e..721004e1 100644 --- a/Grid/qcd/action/fermion/FermionOperatorImpl.h +++ b/Grid/qcd/action/fermion/FermionOperatorImpl.h @@ -240,16 +240,30 @@ namespace QCD { GaugeLinkField tmp(GaugeGrid); Lattice > coor(GaugeGrid); + //////////////////////////////////////////////////// + // apply any boundary phase or twists + //////////////////////////////////////////////////// for (int mu = 0; mu < Nd; mu++) { - auto pha = Params.boundary_phases[mu]; - scalar_type phase( real(pha),imag(pha) ); + ////////// boundary phase ///////////// + auto pha = Params.boundary_phases[mu]; + scalar_type phase( real(pha),imag(pha) ); - int Lmu = GaugeGrid->GlobalDimensions()[mu] - 1; + int L = GaugeGrid->GlobalDimensions()[mu]; + int Lmu = L - 1; LatticeCoordinate(coor, mu); U = PeekIndex(Umu, mu); + + // apply any twists + RealD theta = Params.twist_n_2pi_L[mu] * 2*M_PI / L; + if ( theta != 0.0) { + scalar_type twphase(::cos(theta),::sin(theta)); + U = twphase*U; + std::cout << GridLogMessage << " Twist ["<(Uds, tmp, mu); From 8e0d2f3402e2cc259fd84382627b590759c4f9a3 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Mon, 12 Nov 2018 17:16:18 +0000 Subject: [PATCH 36/37] Hadrons: support for twisted boundary conditions --- Hadrons/Modules/MAction/DWF.hpp | 8 +++++--- Hadrons/Modules/MAction/MobiusDWF.hpp | 8 +++++--- Hadrons/Modules/MAction/ScaledDWF.hpp | 8 +++++--- Hadrons/Modules/MAction/Wilson.hpp | 9 ++++++--- Hadrons/Modules/MAction/WilsonClover.hpp | 8 +++++--- Hadrons/Modules/MAction/ZMobiusDWF.hpp | 8 +++++--- tests/hadrons/Test_QED.cc | 2 ++ tests/hadrons/Test_hadrons.hpp | 2 ++ tests/hadrons/Test_hadrons_meson_3pt.cc | 2 ++ tests/hadrons/Test_hadrons_spectrum.cc | 2 ++ 10 files changed, 39 insertions(+), 18 deletions(-) diff --git a/Hadrons/Modules/MAction/DWF.hpp b/Hadrons/Modules/MAction/DWF.hpp index 257782a1..a7104b42 100644 --- a/Hadrons/Modules/MAction/DWF.hpp +++ b/Hadrons/Modules/MAction/DWF.hpp @@ -49,7 +49,8 @@ public: unsigned int, Ls, double , mass, double , M5, - std::string , boundary); + std::string , boundary, + std::string , twist); }; template @@ -119,8 +120,9 @@ void TDWF::setup(void) auto &grb4 = *envGetRbGrid(FermionField); auto &g5 = *envGetGrid(FermionField, par().Ls); auto &grb5 = *envGetRbGrid(FermionField, par().Ls); - std::vector boundary = strToVec(par().boundary); - typename DomainWallFermion::ImplParams implParams(boundary); + typename DomainWallFermion::ImplParams implParams; + implParams.boundary_phases = strToVec(par().boundary); + implParams.twist_n_2pi_L = strToVec(par().twist); envCreateDerived(FMat, DomainWallFermion, getName(), par().Ls, U, g5, grb5, g4, grb4, par().mass, par().M5, implParams); } diff --git a/Hadrons/Modules/MAction/MobiusDWF.hpp b/Hadrons/Modules/MAction/MobiusDWF.hpp index 01223267..0ba9c4c3 100644 --- a/Hadrons/Modules/MAction/MobiusDWF.hpp +++ b/Hadrons/Modules/MAction/MobiusDWF.hpp @@ -49,7 +49,8 @@ public: double , M5, double , b, double , c, - std::string , boundary); + std::string , boundary, + std::string , twist); }; template @@ -119,8 +120,9 @@ void TMobiusDWF::setup(void) auto &grb4 = *envGetRbGrid(FermionField); auto &g5 = *envGetGrid(FermionField, par().Ls); auto &grb5 = *envGetRbGrid(FermionField, par().Ls); - std::vector boundary = strToVec(par().boundary); - typename MobiusFermion::ImplParams implParams(boundary); + typename MobiusFermion::ImplParams implParams; + implParams.boundary_phases = strToVec(par().boundary); + implParams.twist_n_2pi_L = strToVec(par().twist); envCreateDerived(FMat, MobiusFermion, getName(), par().Ls, U, g5, grb5, g4, grb4, par().mass, par().M5, par().b, par().c, implParams); diff --git a/Hadrons/Modules/MAction/ScaledDWF.hpp b/Hadrons/Modules/MAction/ScaledDWF.hpp index f7a09eef..3b8066e6 100644 --- a/Hadrons/Modules/MAction/ScaledDWF.hpp +++ b/Hadrons/Modules/MAction/ScaledDWF.hpp @@ -48,7 +48,8 @@ public: double , mass, double , M5, double , scale, - std::string , boundary); + std::string , boundary, + std::string , twist); }; template @@ -118,8 +119,9 @@ void TScaledDWF::setup(void) auto &grb4 = *envGetRbGrid(FermionField); auto &g5 = *envGetGrid(FermionField, par().Ls); auto &grb5 = *envGetRbGrid(FermionField, par().Ls); - std::vector boundary = strToVec(par().boundary); - typename MobiusFermion::ImplParams implParams(boundary); + typename ScaledShamirFermion::ImplParams implParams; + implParams.boundary_phases = strToVec(par().boundary); + implParams.twist_n_2pi_L = strToVec(par().twist); envCreateDerived(FMat, ScaledShamirFermion, getName(), par().Ls, U, g5, grb5, g4, grb4, par().mass, par().M5, par().scale, implParams); diff --git a/Hadrons/Modules/MAction/Wilson.hpp b/Hadrons/Modules/MAction/Wilson.hpp index 093449e6..b5e53837 100644 --- a/Hadrons/Modules/MAction/Wilson.hpp +++ b/Hadrons/Modules/MAction/Wilson.hpp @@ -47,7 +47,9 @@ public: GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonPar, std::string, gauge, double , mass, - std::string, boundary); + std::string, boundary, + std::string, string, + std::string, twist); }; template @@ -113,8 +115,9 @@ void TWilson::setup(void) auto &U = envGet(GaugeField, par().gauge); auto &grid = *envGetGrid(FermionField); auto &gridRb = *envGetRbGrid(FermionField); - std::vector boundary = strToVec(par().boundary); - typename WilsonFermion::ImplParams implParams(boundary); + typename WilsonFermion::ImplParams implParams; + implParams.boundary_phases = strToVec(par().boundary); + implParams.twist_n_2pi_L = strToVec(par().twist); envCreateDerived(FMat, WilsonFermion, getName(), 1, U, grid, gridRb, par().mass, implParams); } diff --git a/Hadrons/Modules/MAction/WilsonClover.hpp b/Hadrons/Modules/MAction/WilsonClover.hpp index 0b78bb55..ad301380 100644 --- a/Hadrons/Modules/MAction/WilsonClover.hpp +++ b/Hadrons/Modules/MAction/WilsonClover.hpp @@ -51,7 +51,8 @@ public: double , csw_r, double , csw_t, WilsonAnisotropyCoefficients ,clover_anisotropy, - std::string, boundary + std::string, boundary, + std::string, twist ); }; @@ -119,8 +120,9 @@ void TWilsonClover::setup(void) auto &U = envGet(GaugeField, par().gauge); auto &grid = *envGetGrid(FermionField); auto &gridRb = *envGetRbGrid(FermionField); - std::vector boundary = strToVec(par().boundary); - typename WilsonCloverFermion::ImplParams implParams(boundary); + typename WilsonCloverFermion::ImplParams implParams; + implParams.boundary_phases = strToVec(par().boundary); + implParams.twist_n_2pi_L = strToVec(par().twist); envCreateDerived(FMat, WilsonCloverFermion, getName(), 1, U, grid, gridRb, par().mass, par().csw_r, par().csw_t, par().clover_anisotropy, implParams); diff --git a/Hadrons/Modules/MAction/ZMobiusDWF.hpp b/Hadrons/Modules/MAction/ZMobiusDWF.hpp index 40b04d7f..12ad82ea 100644 --- a/Hadrons/Modules/MAction/ZMobiusDWF.hpp +++ b/Hadrons/Modules/MAction/ZMobiusDWF.hpp @@ -50,7 +50,8 @@ public: double , b, double , c, std::vector>, omega, - std::string , boundary); + std::string , boundary, + std::string , twist); }; template @@ -127,8 +128,9 @@ void TZMobiusDWF::setup(void) auto &g5 = *envGetGrid(FermionField, par().Ls); auto &grb5 = *envGetRbGrid(FermionField, par().Ls); auto omega = par().omega; - std::vector boundary = strToVec(par().boundary); - typename ZMobiusFermion::ImplParams implParams(boundary); + typename ZMobiusFermion::ImplParams implParams; + implParams.boundary_phases = strToVec(par().boundary); + implParams.twist_n_2pi_L = strToVec(par().twist); envCreateDerived(FMat, ZMobiusFermion, getName(), par().Ls, U, g5, grb5, g4, grb4, par().mass, par().M5, omega, par().b, par().c, implParams); diff --git a/tests/hadrons/Test_QED.cc b/tests/hadrons/Test_QED.cc index fb6f8198..9053c6ed 100644 --- a/tests/hadrons/Test_QED.cc +++ b/tests/hadrons/Test_QED.cc @@ -72,6 +72,7 @@ int main(int argc, char *argv[]) // set fermion boundary conditions to be periodic space, antiperiodic time. std::string boundary = "1 1 1 -1"; + std::string twist = "0. 0. 0. 0."; //stochastic photon field MGauge::StochEm::Par photonPar; @@ -90,6 +91,7 @@ int main(int argc, char *argv[]) actionPar.M5 = 1.8; actionPar.mass = mass[i]; actionPar.boundary = boundary; + actionPar.twist = "0. 0. 0. 0."; application.createModule("DWF_" + flavour[i], actionPar); diff --git a/tests/hadrons/Test_hadrons.hpp b/tests/hadrons/Test_hadrons.hpp index 1ee72356..7f07bea5 100644 --- a/tests/hadrons/Test_hadrons.hpp +++ b/tests/hadrons/Test_hadrons.hpp @@ -126,6 +126,7 @@ inline void makeWilsonAction(Application &application, std::string actionName, actionPar.gauge = gaugeField; actionPar.mass = mass; actionPar.boundary = boundary; + actionPar.twist = "0. 0. 0. 0."; application.createModule(actionName, actionPar); } } @@ -154,6 +155,7 @@ inline void makeDWFAction(Application &application, std::string actionName, actionPar.M5 = M5; actionPar.mass = mass; actionPar.boundary = boundary; + actionPar.twist = "0. 0. 0. 0."; application.createModule(actionName, actionPar); } } diff --git a/tests/hadrons/Test_hadrons_meson_3pt.cc b/tests/hadrons/Test_hadrons_meson_3pt.cc index 5b7cda22..741e2a7c 100644 --- a/tests/hadrons/Test_hadrons_meson_3pt.cc +++ b/tests/hadrons/Test_hadrons_meson_3pt.cc @@ -66,6 +66,7 @@ int main(int argc, char *argv[]) // set fermion boundary conditions to be periodic space, antiperiodic time. std::string boundary = "1 1 1 -1"; + std::string twist = "0. 0. 0. 0."; // sink MSink::Point::Par sinkPar; @@ -80,6 +81,7 @@ int main(int argc, char *argv[]) actionPar.M5 = 1.8; actionPar.mass = mass[i]; actionPar.boundary = boundary; + actionPar.twist = twist; application.createModule("DWF_" + flavour[i], actionPar); // solvers diff --git a/tests/hadrons/Test_hadrons_spectrum.cc b/tests/hadrons/Test_hadrons_spectrum.cc index b57deb38..af1dccd7 100644 --- a/tests/hadrons/Test_hadrons_spectrum.cc +++ b/tests/hadrons/Test_hadrons_spectrum.cc @@ -72,6 +72,7 @@ int main(int argc, char *argv[]) // set fermion boundary conditions to be periodic space, antiperiodic time. std::string boundary = "1 1 1 -1"; + std::string twist = "0. 0. 0. 0."; for (unsigned int i = 0; i < flavour.size(); ++i) { @@ -82,6 +83,7 @@ int main(int argc, char *argv[]) actionPar.M5 = 1.8; actionPar.mass = mass[i]; actionPar.boundary = boundary; + actionPar.twist = twist; application.createModule("DWF_" + flavour[i], actionPar); // solvers From 995f20e45d40be9a762856f1bfec89a3f71287e4 Mon Sep 17 00:00:00 2001 From: Vera Guelpers Date: Tue, 13 Nov 2018 14:54:48 +0000 Subject: [PATCH 37/37] Hadrons: some renamings --- Hadrons/Modules/MGauge/Electrify.hpp | 10 +++++----- Hadrons/Modules/MSolver/A2AAslashVector.hpp | 20 ++++++++++---------- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/Hadrons/Modules/MGauge/Electrify.hpp b/Hadrons/Modules/MGauge/Electrify.hpp index c5e6735d..58d65eba 100644 --- a/Hadrons/Modules/MGauge/Electrify.hpp +++ b/Hadrons/Modules/MGauge/Electrify.hpp @@ -48,7 +48,7 @@ BEGIN_MODULE_NAMESPACE(MGauge) * with * * - gauge: U_mu(x): gauge field -* - photon: A_mu(x): electromagnetic photon field +* - emField: A_mu(x): electromagnetic photon field * - e: value for the elementary charge * - q: charge in units of e * @@ -60,7 +60,7 @@ class ElectrifyPar: Serializable public: GRID_SERIALIZABLE_CLASS_MEMBERS(ElectrifyPar, std::string, gauge, - std::string, photon, + std::string, emField, double, e, double, charge); }; @@ -102,7 +102,7 @@ TElectrify::TElectrify(const std::string name) template std::vector TElectrify::getInput(void) { - std::vector in = {par().gauge, par().photon}; + std::vector in = {par().gauge, par().emField}; return in; } @@ -128,11 +128,11 @@ template void TElectrify::execute(void) { LOG(Message) << "Electrify the gauge field " << par().gauge << " using the photon field " - << par().photon << " with charge e*q= " << par().e << "*" << par().charge << std::endl; + << par().emField << " with charge e*q= " << par().e << "*" << par().charge << std::endl; auto &Ue = envGet(GaugeField, getName()); auto &U = envGet(GaugeField, par().gauge); - auto &A = envGet(EmField, par().photon); + auto &A = envGet(EmField, par().emField); envGetTmp(LatticeComplex, eiAmu); Complex i(0.0,1.0); diff --git a/Hadrons/Modules/MSolver/A2AAslashVector.hpp b/Hadrons/Modules/MSolver/A2AAslashVector.hpp index c99a6f6b..3142d1e7 100644 --- a/Hadrons/Modules/MSolver/A2AAslashVector.hpp +++ b/Hadrons/Modules/MSolver/A2AAslashVector.hpp @@ -48,8 +48,8 @@ BEGIN_MODULE_NAMESPACE(MSolver) * * with * -* - v: A2A vector v_i(x) -* - photon: A_mu(x): electromagnetic photon field +* - vector: A2A vector v_i(x) +* - emField: A_mu(x): electromagnetic photon field * - solver: the solver for calculating the sequential propagator * *****************************************************************************/ @@ -59,8 +59,8 @@ class A2AAslashVectorPar: Serializable { public: GRID_SERIALIZABLE_CLASS_MEMBERS(A2AAslashVectorPar, - std::string, v, - std::string, photon, + std::string, vector, + std::string, emField, std::string, solver); }; @@ -104,7 +104,7 @@ TA2AAslashVector::TA2AAslashVector(const std::string name) template std::vector TA2AAslashVector::getInput(void) { - std::vector in = {par().v, par().photon, par().solver}; + std::vector in = {par().vector, par().emField, par().solver}; return in; } @@ -122,7 +122,7 @@ template void TA2AAslashVector::setup(void) { Ls_ = env().getObjectLs(par().solver); - auto &vvector = envGet(std::vector, par().v); + auto &vvector = envGet(std::vector, par().vector); unsigned int Nmodes = vvector.size(); envCreate(std::vector, getName(), 1, Nmodes, envGetGrid(FermionField)); @@ -137,8 +137,8 @@ template void TA2AAslashVector::execute(void) { auto &solver = envGet(Solver, par().solver); - auto &stoch_photon = envGet(EmField, par().photon); - auto &vvector = envGet(std::vector, par().v); + auto &stoch_photon = envGet(EmField, par().emField); + auto &vvector = envGet(std::vector, par().vector); auto &Aslashv = envGet(std::vector, getName()); unsigned int Nmodes = vvector.size(); auto &mat = solver.getFMat(); @@ -151,8 +151,8 @@ void TA2AAslashVector::execute(void) startTimer("Seq Aslash"); - LOG(Message) << "Calculate Sequential propagator on Aslash * v with the A2A vector " << par().v - << " and the photon field " << par().photon << std::endl; + LOG(Message) << "Calculate Sequential propagator on Aslash * v with the A2A vector " << par().vector + << " and the photon field " << par().emField << std::endl; for(unsigned int i=0; i