From 1dc86efd26244d59c4db442ca6ee2853ac2d488d Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 5 Mar 2018 12:22:18 +0000 Subject: [PATCH 01/29] Finalize protection --- lib/communicator/SharedMemoryMPI.cc | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/lib/communicator/SharedMemoryMPI.cc b/lib/communicator/SharedMemoryMPI.cc index 9e5d8f15..45edbb07 100644 --- a/lib/communicator/SharedMemoryMPI.cc +++ b/lib/communicator/SharedMemoryMPI.cc @@ -401,7 +401,10 @@ void *SharedMemory::ShmBufferTranslate(int rank,void * local_p) } SharedMemory::~SharedMemory() { - MPI_Comm_free(&ShmComm); + int MPI_is_finalised; MPI_Finalized(&MPI_is_finalised); + if ( !MPI_is_finalised ) { + MPI_Comm_free(&ShmComm); + } }; } From c399c2b44dea7e6cc4ca6ee34adcd1a86b07c338 Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 5 Mar 2018 12:55:41 +0000 Subject: [PATCH 02/29] Guido broke the charge conjugate plaquette action with premature optimisation. This sector of the code does not matter for anything other than Guido's quenched HMC studies, and any plaq specific optimisations should be retained in a private branch instead of destroying the code simplicity. --- lib/qcd/action/gauge/WilsonGaugeAction.h | 12 ++++-------- lib/qcd/utils/WilsonLoops.h | 5 +++-- tests/forces/Test_gp_rect_force.cc | 4 ++-- tests/forces/Test_gpwilson_force.cc | 2 +- 4 files changed, 10 insertions(+), 13 deletions(-) diff --git a/lib/qcd/action/gauge/WilsonGaugeAction.h b/lib/qcd/action/gauge/WilsonGaugeAction.h index 1ea780b7..77c2424c 100644 --- a/lib/qcd/action/gauge/WilsonGaugeAction.h +++ b/lib/qcd/action/gauge/WilsonGaugeAction.h @@ -71,18 +71,14 @@ class WilsonGaugeAction : public Action { RealD factor = 0.5 * beta / RealD(Nc); - //GaugeLinkField Umu(U._grid); + GaugeLinkField Umu(U._grid); GaugeLinkField dSdU_mu(U._grid); for (int mu = 0; mu < Nd; mu++) { - //Umu = PeekIndex(U, mu); + Umu = PeekIndex(U, mu); // Staple in direction mu - //WilsonLoops::Staple(dSdU_mu, U, mu); - //dSdU_mu = Ta(Umu * dSdU_mu) * factor; - - - WilsonLoops::StapleMult(dSdU_mu, U, mu); - dSdU_mu = Ta(dSdU_mu) * factor; + WilsonLoops::Staple(dSdU_mu, U, mu); + dSdU_mu = Ta(Umu * dSdU_mu) * factor; PokeIndex(dSdU, dSdU_mu, mu); } diff --git a/lib/qcd/utils/WilsonLoops.h b/lib/qcd/utils/WilsonLoops.h index cdd76ecc..6cf34e0c 100644 --- a/lib/qcd/utils/WilsonLoops.h +++ b/lib/qcd/utils/WilsonLoops.h @@ -212,6 +212,7 @@ public: // For the force term +/* static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) { GridBase *grid = Umu._grid; std::vector U(Nd, grid); @@ -225,7 +226,7 @@ static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) { for (int nu = 0; nu < Nd; nu++) { if (nu != mu) { - // this is ~10% faster than the Staple + // this is ~10% faster than the Staple -- PAB: so what it gives the WRONG answers for other BC's! tmp1 = Cshift(U[nu], mu, 1); tmp2 = Cshift(U[mu], nu, 1); staple += tmp1* adj(U[nu]*tmp2); @@ -235,7 +236,7 @@ static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) { } staple = U[mu]*staple; } - +*/ ////////////////////////////////////////////////// // the sum over all staples on each site ////////////////////////////////////////////////// diff --git a/tests/forces/Test_gp_rect_force.cc b/tests/forces/Test_gp_rect_force.cc index bb35c77a..6b3349e0 100644 --- a/tests/forces/Test_gp_rect_force.cc +++ b/tests/forces/Test_gp_rect_force.cc @@ -59,8 +59,8 @@ int main (int argc, char ** argv) double beta = 1.0; double c1 = 0.331; - //GparityPlaqPlusRectangleActionR Action(beta,c1); - ConjugateWilsonGaugeActionR Action(beta); + ConjugatePlaqPlusRectangleActionR Action(beta,c1); + // ConjugateWilsonGaugeActionR Action(beta); //WilsonGaugeActionR Action(beta); ComplexD S = Action.S(U); diff --git a/tests/forces/Test_gpwilson_force.cc b/tests/forces/Test_gpwilson_force.cc index ebde61a5..e52ed7ee 100644 --- a/tests/forces/Test_gpwilson_force.cc +++ b/tests/forces/Test_gpwilson_force.cc @@ -91,7 +91,7 @@ int main (int argc, char ** argv) //////////////////////////////////// // Modify the gauge field a little //////////////////////////////////// - RealD dt = 0.0001; + RealD dt = 0.01; LatticeColourMatrix mommu(UGrid); LatticeColourMatrix forcemu(UGrid); From cd51b9af99e0226a22112592682eb3bfe6c63fc2 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Mon, 5 Mar 2018 19:58:13 +0000 Subject: [PATCH 03/29] Torture yourself with namespace lookup 101 --- extras/Hadrons/Global.hpp | 7 +++- lib/serialisation/BaseIO.h | 82 +++++++++++++++++++------------------- 2 files changed, 46 insertions(+), 43 deletions(-) diff --git a/extras/Hadrons/Global.hpp b/extras/Hadrons/Global.hpp index ed8f4f32..e9f5933b 100644 --- a/extras/Hadrons/Global.hpp +++ b/extras/Hadrons/Global.hpp @@ -43,12 +43,15 @@ See the full license in the file "LICENSE" in the top level distribution directo namespace Grid {\ using namespace QCD;\ namespace Hadrons {\ -using Grid::operator<<; +using Grid::operator<<;\ +using Grid::operator>>; #define END_HADRONS_NAMESPACE }} #define BEGIN_MODULE_NAMESPACE(name)\ namespace name {\ -using Grid::operator<<; +using Grid::operator<<;\ +using Grid::operator>>; + #define END_MODULE_NAMESPACE } /* the 'using Grid::operator<<;' statement prevents a very nasty compilation diff --git a/lib/serialisation/BaseIO.h b/lib/serialisation/BaseIO.h index 24e1cec7..d2db8bfc 100644 --- a/lib/serialisation/BaseIO.h +++ b/lib/serialisation/BaseIO.h @@ -33,42 +33,6 @@ Author: Guido Cossu #include namespace Grid { - // Vector IO utilities /////////////////////////////////////////////////////// - // helper function to read space-separated values - template - std::vector strToVec(const std::string s) - { - std::istringstream sstr(s); - T buf; - std::vector v; - - while(!sstr.eof()) - { - sstr >> buf; - v.push_back(buf); - } - - return v; - } - - // output to streams for vectors - template < class T > - inline std::ostream & operator<<(std::ostream &os, const std::vector &v) - { - os << "["; - for (auto &x: v) - { - os << x << " "; - } - if (v.size() > 0) - { - os << "\b"; - } - os << "]"; - - return os; - } - // Vector element trait ////////////////////////////////////////////////////// template struct element @@ -151,15 +115,15 @@ namespace Grid { do { is.get(c); - } while (c != '<' && !is.eof()); - if (c == '<') + } while (c != '{' && !is.eof()); + if (c == '{') { int start = is.tellg(); do { is.get(c); - } while (c != '>' && !is.eof()); - if (c == '>') + } while (c != '}' && !is.eof()); + if (c == '}') { int end = is.tellg(); int psize = end - start - 1; @@ -182,7 +146,43 @@ namespace Grid { template inline std::ostream & operator<<(std::ostream &os, const std::pair &p) { - os << "<" << p.first << " " << p.second << ">"; + os << "{" << p.first << " " << p.second << "}"; + return os; + } + + // Vector IO utilities /////////////////////////////////////////////////////// + // helper function to read space-separated values + template + std::vector strToVec(const std::string s) + { + std::istringstream sstr(s); + T buf; + std::vector v; + + while(!sstr.eof()) + { + sstr >> buf; + v.push_back(buf); + } + + return v; + } + + // output to streams for vectors + template < class T > + inline std::ostream & operator<<(std::ostream &os, const std::vector &v) + { + os << "["; + for (auto &x: v) + { + os << x << " "; + } + if (v.size() > 0) + { + os << "\b"; + } + os << "]"; + return os; } From dda6c69d5b32aa7d108b225928b76b8e6748b265 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Mon, 5 Mar 2018 19:58:40 +0000 Subject: [PATCH 04/29] Hadrons: scalar SU(N) shift probes --- extras/Hadrons/Modules.hpp | 1 + .../Hadrons/Modules/MScalarSUN/ShiftProbe.hpp | 169 ++++++++++++++++++ extras/Hadrons/modules.inc | 1 + 3 files changed, 171 insertions(+) create mode 100644 extras/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 695dfd02..ee53faa8 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -61,6 +61,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include #include #include #include diff --git a/extras/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp b/extras/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp new file mode 100644 index 00000000..89059180 --- /dev/null +++ b/extras/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp @@ -0,0 +1,169 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MScalarSUN/ShiftProbe.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_ShiftProbe_hpp_ +#define Hadrons_MScalarSUN_ShiftProbe_hpp_ + +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Ward identity phi^n probe with fields at different positions * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MScalarSUN) + +typedef std::pair ShiftPair; + +class ShiftProbePar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(ShiftProbePar, + std::string, field, + std::string, shifts, + std::string, output); +}; + +template +class TShiftProbe: public Module +{ +public: + + typedef typename SImpl::Field Field; + typedef typename SImpl::ComplexField ComplexField; + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::string, op, + Complex , value); + }; +public: + // constructor + TShiftProbe(const std::string name); + // destructor + virtual ~TShiftProbe(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(ShiftProbeSU2, TShiftProbe>, MScalarSUN); +MODULE_REGISTER_NS(ShiftProbeSU3, TShiftProbe>, MScalarSUN); +MODULE_REGISTER_NS(ShiftProbeSU4, TShiftProbe>, MScalarSUN); +MODULE_REGISTER_NS(ShiftProbeSU5, TShiftProbe>, MScalarSUN); +MODULE_REGISTER_NS(ShiftProbeSU6, TShiftProbe>, MScalarSUN); + +/****************************************************************************** + * TShiftProbe implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TShiftProbe::TShiftProbe(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TShiftProbe::getInput(void) +{ + std::vector in = {par().field}; + + return in; +} + +template +std::vector TShiftProbe::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TShiftProbe::setup(void) +{ + envTmpLat(Field, "acc"); + envCacheLat(ComplexField, getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TShiftProbe::execute(void) +{ + LOG(Message) << "Creating shift probe for shifts " << par().shifts + << std::endl; + + std::vector shift; + unsigned int sign; + auto &phi = envGet(Field, par().field); + auto &probe = envGet(ComplexField, getName()); + + shift = strToVec(par().shifts); + if (shift.size() % 2 != 0) + { + HADRON_ERROR(Size, "the number of shifts is odd"); + } + sign = (shift.size() % 4 == 0) ? 1 : -1; + for (auto &s: shift) + { + if (s.first >= env().getNd()) + { + HADRON_ERROR(Size, "dimension to large for shift <" + + std::to_string(s.first) + " " + + std::to_string(s.second) + ">" ); + } + } + envGetTmp(Field, acc); + acc = 1.; + for (unsigned int i = 0; i < shift.size(); ++i) + { + if (shift[i].second == 0) + { + acc *= phi; + } + else + { + acc *= Cshift(phi, shift[i].first, shift[i].second); + } + } + probe = sign*trace(acc); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MScalarSUN_ShiftProbe_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 44df3091..c0be8aef 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -44,6 +44,7 @@ modules_hpp =\ Modules/MAction/Wilson.hpp \ Modules/MAction/WilsonClover.hpp \ Modules/MAction/ZMobiusDWF.hpp \ + Modules/MScalarSUN/ShiftProbe.hpp \ Modules/MScalarSUN/Div.hpp \ Modules/MScalarSUN/TrMag.hpp \ Modules/MScalarSUN/EMT.hpp \ From 485c5db0fe28b04c867caf33c879c58f9b924d96 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 6 Mar 2018 19:22:03 +0000 Subject: [PATCH 05/29] conversion of Grid tensors to nested std::vector in preparation for tensor serialisation --- lib/serialisation/BaseIO.h | 35 +++++++++++++++ tests/IO/Test_serialisation.cc | 78 +++++++--------------------------- 2 files changed, 51 insertions(+), 62 deletions(-) diff --git a/lib/serialisation/BaseIO.h b/lib/serialisation/BaseIO.h index 24e1cec7..1af5acc6 100644 --- a/lib/serialisation/BaseIO.h +++ b/lib/serialisation/BaseIO.h @@ -31,6 +31,7 @@ Author: Guido Cossu #define GRID_SERIALISATION_ABSTRACT_READER_H #include +#include namespace Grid { // Vector IO utilities /////////////////////////////////////////////////////// @@ -69,6 +70,40 @@ namespace Grid { return os; } + // convert Grid scalar tensors to std::vectors + template + void tensorToVec(V &v, const T& t) + { + v = t; + } + + template + void tensorToVec(V &v, const iScalar& t) + { + tensorToVec(v, t._internal); + } + + template + void tensorToVec(std::vector &v, const iVector& t) + { + v.resize(N); + for (unsigned int i = 0; i < N; i++) + { + tensorToVec(v[i], t._internal[i]); + } + } + + template + void tensorToVec(std::vector> &v, const iMatrix& t) + { + v.resize(N, std::vector(N)); + for (unsigned int i = 0; i < N; i++) + for (unsigned int j = 0; j < N; j++) + { + tensorToVec(v[i][j], t._internal[i][j]); + } + } + // Vector element trait ////////////////////////////////////////////////////// template struct element diff --git a/tests/IO/Test_serialisation.cc b/tests/IO/Test_serialisation.cc index 82638ad9..cdafd5c0 100644 --- a/tests/IO/Test_serialisation.cc +++ b/tests/IO/Test_serialisation.cc @@ -197,68 +197,22 @@ int main(int argc,char **argv) std::cout << flatdv.getVector() << std::endl; std::cout << std::endl; + std::cout << "==== Grid tensor to vector test" << std::endl; - std::cout << ".:::::: Testing JSON classes "<< std::endl; - - - { - JSONWriter JW("bother.json"); - - // test basic type writing - myenum a = myenum::red; - push(JW,"BasicTypes"); - write(JW,std::string("i16"),i16); - write(JW,"myenum",a); - write(JW,"u16",u16); - write(JW,"i32",i32); - write(JW,"u32",u32); - write(JW,"i64",i64); - write(JW,"u64",u64); - write(JW,"f",f); - write(JW,"d",d); - write(JW,"b",b); - pop(JW); - - - // test serializable class writing - myclass obj(1234); // non-trivial constructor - std::cout << obj << std::endl; - std::cout << "-- serialisable class writing to 'bother.json'..." << std::endl; - write(JW,"obj",obj); - JW.write("obj2", obj); - - - std::vector vec; - vec.push_back(myclass(1234)); - vec.push_back(myclass(5678)); - vec.push_back(myclass(3838)); - write(JW, "objvec", vec); - - } - - - { - JSONReader RD("bother.json"); - myclass jcopy1; - std::vector jveccopy1; - read(RD,"obj",jcopy1); - read(RD,"objvec", jveccopy1); - std::cout << "Loaded (JSON) -----------------" << std::endl; - std::cout << jcopy1 << std::endl << jveccopy1 << std::endl; - } - - -/* - // This is still work in progress - { - // Testing the next element function - JSONReader RD("test.json"); - RD.push("grid"); - RD.push("Observable"); - std::string name; - read(RD,"name", name); - } -*/ - + GridSerialRNG rng; + SpinColourMatrix scm; + SpinColourVector scv; + std::vector>>> scmv; + std::vector> scvv; + rng.SeedFixedIntegers(std::vector({42,10,81,9})); + random(rng, scm); + random(rng, scv); + std::cout << "Test spin-color matrix: " << scm << std::endl; + std::cout << "Test spin-color vector: " << scv << std::endl; + std::cout << "Converting to std::vector" << std::endl; + tensorToVec(scmv, scm); + tensorToVec(scvv, scv); + std::cout << "Spin-color matrix: " << scmv << std::endl; + std::cout << "Spin-color vector: " << scvv << std::endl; } From 8b14096990ff0fe1969ace0bad933ff3dbbac8fc Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Wed, 7 Mar 2018 15:12:18 +0000 Subject: [PATCH 06/29] Conversion of Grid tensors to std::vector made more elegant, also pair syntax changed to (x y) to avoid issues with JSON/XML --- lib/serialisation/BaseIO.h | 128 ++++++++++++++++++++++++------------- 1 file changed, 83 insertions(+), 45 deletions(-) diff --git a/lib/serialisation/BaseIO.h b/lib/serialisation/BaseIO.h index 1af5acc6..5b37e1fb 100644 --- a/lib/serialisation/BaseIO.h +++ b/lib/serialisation/BaseIO.h @@ -34,74 +34,76 @@ Author: Guido Cossu #include namespace Grid { - // Vector IO utilities /////////////////////////////////////////////////////// - // helper function to read space-separated values + // Grid scalar tensors to nested std::vectors ////////////////////////////////// template - std::vector strToVec(const std::string s) + struct TensorToVec { - std::istringstream sstr(s); - T buf; - std::vector v; - - while(!sstr.eof()) - { - sstr >> buf; - v.push_back(buf); - } - - return v; - } - - // output to streams for vectors - template < class T > - inline std::ostream & operator<<(std::ostream &os, const std::vector &v) + typedef T type; + }; + + template + struct TensorToVec> { - os << "["; - for (auto &x: v) - { - os << x << " "; - } - if (v.size() > 0) - { - os << "\b"; - } - os << "]"; - - return os; - } - - // convert Grid scalar tensors to std::vectors - template - void tensorToVec(V &v, const T& t) + typedef TensorToVec::type type; + }; + + template + struct TensorToVec> + { + typedef TensorToVec::type type; + }; + + template + struct TensorToVec> + { + typedef std::vector::type> type; + }; + + template + struct TensorToVec> + { + typedef std::vector::type>> type; + }; + + template + TensorToVec::type tensorToVec(const T &t) { v = t; } template - void tensorToVec(V &v, const iScalar& t) + TensorToVec>::type tensorToVec(V &v, const iScalar& t) { - tensorToVec(v, t._internal); + return tensorToVec(t._internal); } - template - void tensorToVec(std::vector &v, const iVector& t) + template + TensorToVec>::type tensorToVec(const iVector& t) { + TensorToVec>::type v; + v.resize(N); for (unsigned int i = 0; i < N; i++) { - tensorToVec(v[i], t._internal[i]); + v[i] = tensorToVec(t._internal[i]); } + + return v; } - template - void tensorToVec(std::vector> &v, const iMatrix& t) + template + TensorToVec>::type tensorToVec(const iMatrix& t) { + TensorToVec>::type v; + v.resize(N, std::vector(N)); for (unsigned int i = 0; i < N; i++) for (unsigned int j = 0; j < N; j++) { - tensorToVec(v[i][j], t._internal[i][j]); + v[i][j] = tensorToVec(t._internal[i][j]); } + + return v; } // Vector element trait ////////////////////////////////////////////////////// @@ -217,7 +219,43 @@ namespace Grid { template inline std::ostream & operator<<(std::ostream &os, const std::pair &p) { - os << "<" << p.first << " " << p.second << ">"; + os << "{" << p.first << " " << p.second << "}"; + return os; + } + + // Vector IO utilities /////////////////////////////////////////////////////// + // helper function to read space-separated values + template + std::vector strToVec(const std::string s) + { + std::istringstream sstr(s); + T buf; + std::vector v; + + while(!sstr.eof()) + { + sstr >> buf; + v.push_back(buf); + } + + return v; + } + + // output to streams for vectors + template < class T > + inline std::ostream & operator<<(std::ostream &os, const std::vector &v) + { + os << "["; + for (auto &x: v) + { + os << x << " "; + } + if (v.size() > 0) + { + os << "\b"; + } + os << "]"; + return os; } From 90dbe03e1764c28f3afac7f53576dd4249f07ea8 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Wed, 7 Mar 2018 15:12:18 +0000 Subject: [PATCH 07/29] Conversion of Grid tensors to std::vector made more elegant, also pair syntax changed to (x y) to avoid issues with JSON/XML --- lib/serialisation/BaseIO.h | 141 ++++++++++++++++++++------------- tests/IO/Test_serialisation.cc | 6 +- 2 files changed, 90 insertions(+), 57 deletions(-) diff --git a/lib/serialisation/BaseIO.h b/lib/serialisation/BaseIO.h index 1af5acc6..0a919aab 100644 --- a/lib/serialisation/BaseIO.h +++ b/lib/serialisation/BaseIO.h @@ -34,74 +34,73 @@ Author: Guido Cossu #include namespace Grid { - // Vector IO utilities /////////////////////////////////////////////////////// - // helper function to read space-separated values + // Grid scalar tensors to nested std::vectors ////////////////////////////////// template - std::vector strToVec(const std::string s) + struct TensorToVec { - std::istringstream sstr(s); - T buf; - std::vector v; - - while(!sstr.eof()) - { - sstr >> buf; - v.push_back(buf); - } - - return v; - } - - // output to streams for vectors - template < class T > - inline std::ostream & operator<<(std::ostream &os, const std::vector &v) + typedef T type; + }; + + template + struct TensorToVec> { - os << "["; - for (auto &x: v) - { - os << x << " "; - } - if (v.size() > 0) - { - os << "\b"; - } - os << "]"; - - return os; - } - - // convert Grid scalar tensors to std::vectors - template - void tensorToVec(V &v, const T& t) + typedef typename TensorToVec::type type; + }; + + template + struct TensorToVec> { - v = t; + typedef typename std::vector::type> type; + }; + + template + struct TensorToVec> + { + typedef typename std::vector::type>> type; + }; + + template + typename TensorToVec::type tensorToVec(const T &t) + { + return t; } - template - void tensorToVec(V &v, const iScalar& t) + template + typename TensorToVec>::type tensorToVec(const iScalar& t) { - tensorToVec(v, t._internal); + return tensorToVec(t._internal); } - template - void tensorToVec(std::vector &v, const iVector& t) + template + typename TensorToVec>::type tensorToVec(const iVector& t) { + typename TensorToVec>::type v; + v.resize(N); for (unsigned int i = 0; i < N; i++) { - tensorToVec(v[i], t._internal[i]); + v[i] = tensorToVec(t._internal[i]); } + + return v; } - template - void tensorToVec(std::vector> &v, const iMatrix& t) + template + typename TensorToVec>::type tensorToVec(const iMatrix& t) { - v.resize(N, std::vector(N)); + typename TensorToVec>::type v; + + v.resize(N); for (unsigned int i = 0; i < N; i++) - for (unsigned int j = 0; j < N; j++) { - tensorToVec(v[i][j], t._internal[i][j]); + v[i].resize(N); + for (unsigned int j = 0; j < N; j++) + { + v[i][j] = tensorToVec(t._internal[i][j]); + } } + + return v; } // Vector element trait ////////////////////////////////////////////////////// @@ -186,15 +185,15 @@ namespace Grid { do { is.get(c); - } while (c != '<' && !is.eof()); - if (c == '<') + } while (c != '(' && !is.eof()); + if (c == '(') { int start = is.tellg(); do { is.get(c); - } while (c != '>' && !is.eof()); - if (c == '>') + } while (c != ')' && !is.eof()); + if (c == ')') { int end = is.tellg(); int psize = end - start - 1; @@ -217,7 +216,43 @@ namespace Grid { template inline std::ostream & operator<<(std::ostream &os, const std::pair &p) { - os << "<" << p.first << " " << p.second << ">"; + os << "(" << p.first << " " << p.second << ")"; + return os; + } + + // Vector IO utilities /////////////////////////////////////////////////////// + // helper function to read space-separated values + template + std::vector strToVec(const std::string s) + { + std::istringstream sstr(s); + T buf; + std::vector v; + + while(!sstr.eof()) + { + sstr >> buf; + v.push_back(buf); + } + + return v; + } + + // output to streams for vectors + template < class T > + inline std::ostream & operator<<(std::ostream &os, const std::vector &v) + { + os << "["; + for (auto &x: v) + { + os << x << " "; + } + if (v.size() > 0) + { + os << "\b"; + } + os << "]"; + return os; } diff --git a/tests/IO/Test_serialisation.cc b/tests/IO/Test_serialisation.cc index cdafd5c0..d4b89652 100644 --- a/tests/IO/Test_serialisation.cc +++ b/tests/IO/Test_serialisation.cc @@ -202,8 +202,6 @@ int main(int argc,char **argv) GridSerialRNG rng; SpinColourMatrix scm; SpinColourVector scv; - std::vector>>> scmv; - std::vector> scvv; rng.SeedFixedIntegers(std::vector({42,10,81,9})); random(rng, scm); @@ -211,8 +209,8 @@ int main(int argc,char **argv) std::cout << "Test spin-color matrix: " << scm << std::endl; std::cout << "Test spin-color vector: " << scv << std::endl; std::cout << "Converting to std::vector" << std::endl; - tensorToVec(scmv, scm); - tensorToVec(scvv, scv); + auto scmv = tensorToVec(scm); + auto scvv = tensorToVec(scv); std::cout << "Spin-color matrix: " << scmv << std::endl; std::cout << "Spin-color vector: " << scvv << std::endl; } From 5e8af396fd2855d4cbc84931a1e0b99b86dbbb03 Mon Sep 17 00:00:00 2001 From: Dan H Date: Wed, 7 Mar 2018 13:11:51 -0500 Subject: [PATCH 08/29] Add print of the current git hash on Grid init. --- .gitignore | 1 + Makefile.am | 4 ++++ lib/util/Init.cc | 7 +++++++ 3 files changed, 12 insertions(+) diff --git a/.gitignore b/.gitignore index dc59879f..49295fc6 100644 --- a/.gitignore +++ b/.gitignore @@ -123,6 +123,7 @@ make-bin-BUCK.sh ##################### lib/qcd/spin/gamma-gen/*.h lib/qcd/spin/gamma-gen/*.cc +lib/version.h # vs code editor files # ######################## diff --git a/Makefile.am b/Makefile.am index 3a65cf1b..d507bf08 100644 --- a/Makefile.am +++ b/Makefile.am @@ -5,6 +5,10 @@ include $(top_srcdir)/doxygen.inc bin_SCRIPTS=grid-config +BUILT_SOURCES = version.h + +version.h: + echo "`git log -n 1 --format=format:"#define GITHASH \\"%H:%d\\"%n" HEAD`" > $(srcdir)/lib/version.h .PHONY: bench check tests doxygen-run doxygen-doc $(DX_PS_GOAL) $(DX_PDF_GOAL) diff --git a/lib/util/Init.cc b/lib/util/Init.cc index fb3d7a1e..b4ac14b7 100644 --- a/lib/util/Init.cc +++ b/lib/util/Init.cc @@ -49,6 +49,7 @@ Author: paboyle #include #include +#include #include @@ -288,6 +289,12 @@ void Grid_init(int *argc,char ***argv) std::cout << "but WITHOUT ANY WARRANTY; without even the implied warranty of"< Date: Thu, 8 Mar 2018 09:50:39 +0000 Subject: [PATCH 09/29] std::vector to tensor conversion + test units --- lib/serialisation/BaseIO.h | 32 +++++++++++++++++++++++++++++++ tests/IO/Test_serialisation.cc | 35 +++++++++++++++++++++++----------- 2 files changed, 56 insertions(+), 11 deletions(-) diff --git a/lib/serialisation/BaseIO.h b/lib/serialisation/BaseIO.h index 0a919aab..d129b9e5 100644 --- a/lib/serialisation/BaseIO.h +++ b/lib/serialisation/BaseIO.h @@ -103,6 +103,38 @@ namespace Grid { return v; } + template + void vecToTensor(T &t, const typename TensorToVec::type &v) + { + t = v; + } + + + template + void vecToTensor(iScalar &t, const typename TensorToVec>::type &v) + { + vecToTensor(t._internal, v); + } + + template + void vecToTensor(iVector &t, const typename TensorToVec>::type &v) + { + for (unsigned int i = 0; i < N; i++) + { + vecToTensor(t._internal[i], v[i]); + } + } + + template + void vecToTensor(iMatrix &t, const typename TensorToVec>::type &v) + { + for (unsigned int i = 0; i < N; i++) + for (unsigned int j = 0; j < N; j++) + { + vecToTensor(t._internal[i][j], v[i][j]); + } + } + // Vector element trait ////////////////////////////////////////////////////// template struct element diff --git a/tests/IO/Test_serialisation.cc b/tests/IO/Test_serialisation.cc index d4b89652..93007e44 100644 --- a/tests/IO/Test_serialisation.cc +++ b/tests/IO/Test_serialisation.cc @@ -93,6 +93,24 @@ void ioTest(const std::string &filename, const O &object, const std::string &nam if (!good) exit(EXIT_FAILURE); } +template +void tensorConvTestFn(GridSerialRNG &rng, const std::string label) +{ + T t, ft; + Real n; + bool good; + + random(rng, t); + auto tv = tensorToVec(t); + vecToTensor(ft, tv); + n = norm2(t - ft); + good = (n == 0); + std::cout << label << " norm 2 diff: " << n << " -- " + << (good ? "success" : "failure") << std::endl; +} + +#define tensorConvTest(rng, type) tensorConvTestFn(rng, #type) + int main(int argc,char **argv) { std::cout << "==== basic IO" << std::endl; @@ -200,17 +218,12 @@ int main(int argc,char **argv) std::cout << "==== Grid tensor to vector test" << std::endl; GridSerialRNG rng; - SpinColourMatrix scm; - SpinColourVector scv; rng.SeedFixedIntegers(std::vector({42,10,81,9})); - random(rng, scm); - random(rng, scv); - std::cout << "Test spin-color matrix: " << scm << std::endl; - std::cout << "Test spin-color vector: " << scv << std::endl; - std::cout << "Converting to std::vector" << std::endl; - auto scmv = tensorToVec(scm); - auto scvv = tensorToVec(scv); - std::cout << "Spin-color matrix: " << scmv << std::endl; - std::cout << "Spin-color vector: " << scvv << std::endl; + tensorConvTest(rng, SpinColourMatrix); + tensorConvTest(rng, SpinColourVector); + tensorConvTest(rng, ColourMatrix); + tensorConvTest(rng, ColourVector); + tensorConvTest(rng, SpinMatrix); + tensorConvTest(rng, SpinVector); } From c49be8988be95f37c05741f6e807e41707c847b1 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 8 Mar 2018 09:51:22 +0000 Subject: [PATCH 10/29] Grid tensor serialisation --- lib/serialisation/Hdf5IO.h | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/lib/serialisation/Hdf5IO.h b/lib/serialisation/Hdf5IO.h index 94ad9736..9140435d 100644 --- a/lib/serialisation/Hdf5IO.h +++ b/lib/serialisation/Hdf5IO.h @@ -5,6 +5,7 @@ #include #include #include +#include #include "Hdf5Type.h" #ifndef H5_NO_NAMESPACE @@ -37,6 +38,12 @@ namespace Grid template typename std::enable_if>::is_number, void>::type writeDefault(const std::string &s, const std::vector &x); + template + void writeDefault(const std::string &s, const iScalar &t); + template + void writeDefault(const std::string &s, const iVector &t); + template + void writeDefault(const std::string &s, const iMatrix &t); private: template void writeSingleAttribute(const U &x, const std::string &name, @@ -147,6 +154,24 @@ namespace Grid } pop(); } + + template + void Hdf5Writer::writeDefault(const std::string &s, const iScalar &t) + { + writeDefault(s, tensorToVec(t)); + } + + template + void Hdf5Writer::writeDefault(const std::string &s, const iVector &t) + { + writeDefault(s, tensorToVec(t)); + } + + template + void Hdf5Writer::writeDefault(const std::string &s, const iMatrix &t) + { + writeDefault(s, tensorToVec(t)); + } // Reader template implementation //////////////////////////////////////////// template From 360cface3349f923ef00abaabc064e4e4af94c2b Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 8 Mar 2018 19:12:03 +0000 Subject: [PATCH 11/29] Grid tensor serialisation fully implemented and tested --- lib/lattice/Lattice_comparison_utils.h | 2 +- lib/serialisation/BaseIO.h | 538 +++++++------------------ lib/serialisation/Hdf5IO.h | 24 -- lib/serialisation/VectorUtils.h | 336 +++++++++++++++ lib/tensors/Tensor_logical.h | 33 ++ tests/IO/Test_serialisation.cc | 27 +- 6 files changed, 519 insertions(+), 441 deletions(-) create mode 100644 lib/serialisation/VectorUtils.h diff --git a/lib/lattice/Lattice_comparison_utils.h b/lib/lattice/Lattice_comparison_utils.h index 14a19383..9580d4d2 100644 --- a/lib/lattice/Lattice_comparison_utils.h +++ b/lib/lattice/Lattice_comparison_utils.h @@ -198,7 +198,7 @@ namespace Grid { typedef typename vsimd::scalar_type scalar;\ return Comparison(functor(),lhs,rhs);\ }\ - template\ + template = 0>\ inline vInteger operator op(const iScalar &lhs,const iScalar &rhs)\ { \ return lhs._internal op rhs._internal; \ diff --git a/lib/serialisation/BaseIO.h b/lib/serialisation/BaseIO.h index d129b9e5..c9b3fb9e 100644 --- a/lib/serialisation/BaseIO.h +++ b/lib/serialisation/BaseIO.h @@ -32,178 +32,9 @@ Author: Guido Cossu #include #include +#include namespace Grid { - // Grid scalar tensors to nested std::vectors ////////////////////////////////// - template - struct TensorToVec - { - typedef T type; - }; - - template - struct TensorToVec> - { - typedef typename TensorToVec::type type; - }; - - template - struct TensorToVec> - { - typedef typename std::vector::type> type; - }; - - template - struct TensorToVec> - { - typedef typename std::vector::type>> type; - }; - - template - typename TensorToVec::type tensorToVec(const T &t) - { - return t; - } - - template - typename TensorToVec>::type tensorToVec(const iScalar& t) - { - return tensorToVec(t._internal); - } - - template - typename TensorToVec>::type tensorToVec(const iVector& t) - { - typename TensorToVec>::type v; - - v.resize(N); - for (unsigned int i = 0; i < N; i++) - { - v[i] = tensorToVec(t._internal[i]); - } - - return v; - } - - template - typename TensorToVec>::type tensorToVec(const iMatrix& t) - { - typename TensorToVec>::type v; - - v.resize(N); - for (unsigned int i = 0; i < N; i++) - { - v[i].resize(N); - for (unsigned int j = 0; j < N; j++) - { - v[i][j] = tensorToVec(t._internal[i][j]); - } - } - - return v; - } - - template - void vecToTensor(T &t, const typename TensorToVec::type &v) - { - t = v; - } - - - template - void vecToTensor(iScalar &t, const typename TensorToVec>::type &v) - { - vecToTensor(t._internal, v); - } - - template - void vecToTensor(iVector &t, const typename TensorToVec>::type &v) - { - for (unsigned int i = 0; i < N; i++) - { - vecToTensor(t._internal[i], v[i]); - } - } - - template - void vecToTensor(iMatrix &t, const typename TensorToVec>::type &v) - { - for (unsigned int i = 0; i < N; i++) - for (unsigned int j = 0; j < N; j++) - { - vecToTensor(t._internal[i][j], v[i][j]); - } - } - - // Vector element trait ////////////////////////////////////////////////////// - template - struct element - { - typedef T type; - static constexpr bool is_number = false; - }; - - template - struct element> - { - typedef typename element::type type; - static constexpr bool is_number = std::is_arithmetic::value - or is_complex::value - or element::is_number; - }; - - // Vector flattening utility class //////////////////////////////////////////// - // Class to flatten a multidimensional std::vector - template - class Flatten - { - public: - typedef typename element::type Element; - public: - explicit Flatten(const V &vector); - const V & getVector(void); - const std::vector & getFlatVector(void); - const std::vector & getDim(void); - private: - void accumulate(const Element &e); - template - void accumulate(const W &v); - void accumulateDim(const Element &e); - template - void accumulateDim(const W &v); - private: - const V &vector_; - std::vector flatVector_; - std::vector dim_; - }; - - // Class to reconstruct a multidimensional std::vector - template - class Reconstruct - { - public: - typedef typename element::type Element; - public: - Reconstruct(const std::vector &flatVector, - const std::vector &dim); - const V & getVector(void); - const std::vector & getFlatVector(void); - const std::vector & getDim(void); - private: - void fill(std::vector &v); - template - void fill(W &v); - void resize(std::vector &v, const unsigned int dim); - template - void resize(W &v, const unsigned int dim); - private: - V vector_; - const std::vector &flatVector_; - std::vector dim_; - size_t ind_{0}; - unsigned int dimInd_{0}; - }; - // Pair IO utilities ///////////////////////////////////////////////////////// // helper function to parse input in the format "" template @@ -252,42 +83,6 @@ namespace Grid { return os; } - // Vector IO utilities /////////////////////////////////////////////////////// - // helper function to read space-separated values - template - std::vector strToVec(const std::string s) - { - std::istringstream sstr(s); - T buf; - std::vector v; - - while(!sstr.eof()) - { - sstr >> buf; - v.push_back(buf); - } - - return v; - } - - // output to streams for vectors - template < class T > - inline std::ostream & operator<<(std::ostream &os, const std::vector &v) - { - os << "["; - for (auto &x: v) - { - os << x << " "; - } - if (v.size() > 0) - { - os << "\b"; - } - os << "]"; - - return os; - } - // Abstract writer/reader classes //////////////////////////////////////////// // static polymorphism implemented using CRTP idiom class Serializable; @@ -307,6 +102,12 @@ namespace Grid { template typename std::enable_if::value, void>::type write(const std::string& s, const U &output); + template + void write(const std::string &s, const iScalar &output); + template + void write(const std::string &s, const iVector &output); + template + void write(const std::string &s, const iMatrix &output); private: T *upcast; }; @@ -326,6 +127,12 @@ namespace Grid { template typename std::enable_if::value, void>::type read(const std::string& s, U &output); + template + void read(const std::string &s, iScalar &output); + template + void read(const std::string &s, iVector &output); + template + void read(const std::string &s, iMatrix &output); protected: template void fromString(U &output, const std::string &s); @@ -339,203 +146,9 @@ namespace Grid { }; template struct isWriter { static const bool value = false; - }; - - - - // Generic writer interface - // serializable base class - class Serializable - { - public: - template - static inline void write(Writer &WR,const std::string &s, - const Serializable &obj) - {} - - template - static inline void read(Reader &RD,const std::string &s, - Serializable &obj) - {} - - friend inline std::ostream & operator<<(std::ostream &os, - const Serializable &obj) - { - return os; - } }; - - // Flatten class template implementation ///////////////////////////////////// - template - void Flatten::accumulate(const Element &e) - { - flatVector_.push_back(e); - } - - template - template - void Flatten::accumulate(const W &v) - { - for (auto &e: v) - { - accumulate(e); - } - } - - template - void Flatten::accumulateDim(const Element &e) {}; - - template - template - void Flatten::accumulateDim(const W &v) - { - dim_.push_back(v.size()); - accumulateDim(v[0]); - } - - template - Flatten::Flatten(const V &vector) - : vector_(vector) - { - accumulate(vector_); - accumulateDim(vector_); - } - - template - const V & Flatten::getVector(void) - { - return vector_; - } - - template - const std::vector::Element> & - Flatten::getFlatVector(void) - { - return flatVector_; - } - - template - const std::vector & Flatten::getDim(void) - { - return dim_; - } - - // Reconstruct class template implementation ///////////////////////////////// - template - void Reconstruct::fill(std::vector &v) - { - for (auto &e: v) - { - e = flatVector_[ind_++]; - } - } - - template - template - void Reconstruct::fill(W &v) - { - for (auto &e: v) - { - fill(e); - } - } - - template - void Reconstruct::resize(std::vector &v, const unsigned int dim) - { - v.resize(dim_[dim]); - } - - template - template - void Reconstruct::resize(W &v, const unsigned int dim) - { - v.resize(dim_[dim]); - for (auto &e: v) - { - resize(e, dim + 1); - } - } - - template - Reconstruct::Reconstruct(const std::vector &flatVector, - const std::vector &dim) - : flatVector_(flatVector) - , dim_(dim) - { - resize(vector_, 0); - fill(vector_); - } - - template - const V & Reconstruct::getVector(void) - { - return vector_; - } - - template - const std::vector::Element> & - Reconstruct::getFlatVector(void) - { - return flatVector_; - } - - template - const std::vector & Reconstruct::getDim(void) - { - return dim_; - } - - // Generic writer interface ////////////////////////////////////////////////// - template - inline void push(Writer &w, const std::string &s) { - w.push(s); - } - - template - inline void push(Writer &w, const char *s) - { - w.push(std::string(s)); - } - - template - inline void pop(Writer &w) - { - w.pop(); - } - - template - inline void write(Writer &w, const std::string& s, const U &output) - { - w.write(s, output); - } - - // Generic reader interface - template - inline bool push(Reader &r, const std::string &s) - { - return r.push(s); - } - - template - inline bool push(Reader &r, const char *s) - { - return r.push(std::string(s)); - } - - template - inline void pop(Reader &r) - { - r.pop(); - } - - template - inline void read(Reader &r, const std::string &s, U &output) - { - r.read(s, output); - } - - // Writer template implementation //////////////////////////////////////////// + + // Writer template implementation template Writer::Writer(void) { @@ -569,6 +182,27 @@ namespace Grid { { upcast->writeDefault(s, output); } + + template + template + void Writer::write(const std::string &s, const iScalar &output) + { + upcast->writeDefault(s, tensorToVec(output)); + } + + template + template + void Writer::write(const std::string &s, const iVector &output) + { + upcast->writeDefault(s, tensorToVec(output)); + } + + template + template + void Writer::write(const std::string &s, const iMatrix &output) + { + upcast->writeDefault(s, tensorToVec(output)); + } // Reader template implementation template @@ -604,7 +238,37 @@ namespace Grid { { upcast->readDefault(s, output); } + + template + template + void Reader::read(const std::string &s, iScalar &output) + { + typename TensorToVec>::type v; + + upcast->readDefault(s, v); + vecToTensor(output, v); + } + + template + template + void Reader::read(const std::string &s, iVector &output) + { + typename TensorToVec>::type v; + + upcast->readDefault(s, v); + vecToTensor(output, v); + } + template + template + void Reader::read(const std::string &s, iMatrix &output) + { + typename TensorToVec>::type v; + + upcast->readDefault(s, v); + vecToTensor(output, v); + } + template template void Reader::fromString(U &output, const std::string &s) @@ -623,6 +287,76 @@ namespace Grid { abort(); } } + + // serializable base class /////////////////////////////////////////////////// + class Serializable + { + public: + template + static inline void write(Writer &WR,const std::string &s, + const Serializable &obj) + {} + + template + static inline void read(Reader &RD,const std::string &s, + Serializable &obj) + {} + + friend inline std::ostream & operator<<(std::ostream &os, + const Serializable &obj) + { + return os; + } + }; + + // Generic writer interface ////////////////////////////////////////////////// + template + inline void push(Writer &w, const std::string &s) { + w.push(s); + } + + template + inline void push(Writer &w, const char *s) + { + w.push(std::string(s)); + } + + template + inline void pop(Writer &w) + { + w.pop(); + } + + template + inline void write(Writer &w, const std::string& s, const U &output) + { + w.write(s, output); + } + + // Generic reader interface ////////////////////////////////////////////////// + template + inline bool push(Reader &r, const std::string &s) + { + return r.push(s); + } + + template + inline bool push(Reader &r, const char *s) + { + return r.push(std::string(s)); + } + + template + inline void pop(Reader &r) + { + r.pop(); + } + + template + inline void read(Reader &r, const std::string &s, U &output) + { + r.read(s, output); + } } #endif diff --git a/lib/serialisation/Hdf5IO.h b/lib/serialisation/Hdf5IO.h index 9140435d..12625ab8 100644 --- a/lib/serialisation/Hdf5IO.h +++ b/lib/serialisation/Hdf5IO.h @@ -38,12 +38,6 @@ namespace Grid template typename std::enable_if>::is_number, void>::type writeDefault(const std::string &s, const std::vector &x); - template - void writeDefault(const std::string &s, const iScalar &t); - template - void writeDefault(const std::string &s, const iVector &t); - template - void writeDefault(const std::string &s, const iMatrix &t); private: template void writeSingleAttribute(const U &x, const std::string &name, @@ -154,24 +148,6 @@ namespace Grid } pop(); } - - template - void Hdf5Writer::writeDefault(const std::string &s, const iScalar &t) - { - writeDefault(s, tensorToVec(t)); - } - - template - void Hdf5Writer::writeDefault(const std::string &s, const iVector &t) - { - writeDefault(s, tensorToVec(t)); - } - - template - void Hdf5Writer::writeDefault(const std::string &s, const iMatrix &t) - { - writeDefault(s, tensorToVec(t)); - } // Reader template implementation //////////////////////////////////////////// template diff --git a/lib/serialisation/VectorUtils.h b/lib/serialisation/VectorUtils.h new file mode 100644 index 00000000..f5c76b84 --- /dev/null +++ b/lib/serialisation/VectorUtils.h @@ -0,0 +1,336 @@ +#ifndef GRID_SERIALISATION_VECTORUTILS_H +#define GRID_SERIALISATION_VECTORUTILS_H + +#include +#include + +namespace Grid { + // Grid scalar tensors to nested std::vectors ////////////////////////////////// + template + struct TensorToVec + { + typedef T type; + }; + + template + struct TensorToVec> + { + typedef typename TensorToVec::type type; + }; + + template + struct TensorToVec> + { + typedef typename std::vector::type> type; + }; + + template + struct TensorToVec> + { + typedef typename std::vector::type>> type; + }; + + template + typename TensorToVec::type tensorToVec(const T &t) + { + return t; + } + + template + typename TensorToVec>::type tensorToVec(const iScalar& t) + { + return tensorToVec(t._internal); + } + + template + typename TensorToVec>::type tensorToVec(const iVector& t) + { + typename TensorToVec>::type v; + + v.resize(N); + for (unsigned int i = 0; i < N; i++) + { + v[i] = tensorToVec(t._internal[i]); + } + + return v; + } + + template + typename TensorToVec>::type tensorToVec(const iMatrix& t) + { + typename TensorToVec>::type v; + + v.resize(N); + for (unsigned int i = 0; i < N; i++) + { + v[i].resize(N); + for (unsigned int j = 0; j < N; j++) + { + v[i][j] = tensorToVec(t._internal[i][j]); + } + } + + return v; + } + + template + void vecToTensor(T &t, const typename TensorToVec::type &v) + { + t = v; + } + + + template + void vecToTensor(iScalar &t, const typename TensorToVec>::type &v) + { + vecToTensor(t._internal, v); + } + + template + void vecToTensor(iVector &t, const typename TensorToVec>::type &v) + { + for (unsigned int i = 0; i < N; i++) + { + vecToTensor(t._internal[i], v[i]); + } + } + + template + void vecToTensor(iMatrix &t, const typename TensorToVec>::type &v) + { + for (unsigned int i = 0; i < N; i++) + for (unsigned int j = 0; j < N; j++) + { + vecToTensor(t._internal[i][j], v[i][j]); + } + } + + // Vector element trait ////////////////////////////////////////////////////// + template + struct element + { + typedef T type; + static constexpr bool is_number = false; + }; + + template + struct element> + { + typedef typename element::type type; + static constexpr bool is_number = std::is_arithmetic::value + or is_complex::value + or element::is_number; + }; + + // Vector flattening utility class //////////////////////////////////////////// + // Class to flatten a multidimensional std::vector + template + class Flatten + { + public: + typedef typename element::type Element; + public: + explicit Flatten(const V &vector); + const V & getVector(void); + const std::vector & getFlatVector(void); + const std::vector & getDim(void); + private: + void accumulate(const Element &e); + template + void accumulate(const W &v); + void accumulateDim(const Element &e); + template + void accumulateDim(const W &v); + private: + const V &vector_; + std::vector flatVector_; + std::vector dim_; + }; + + // Class to reconstruct a multidimensional std::vector + template + class Reconstruct + { + public: + typedef typename element::type Element; + public: + Reconstruct(const std::vector &flatVector, + const std::vector &dim); + const V & getVector(void); + const std::vector & getFlatVector(void); + const std::vector & getDim(void); + private: + void fill(std::vector &v); + template + void fill(W &v); + void resize(std::vector &v, const unsigned int dim); + template + void resize(W &v, const unsigned int dim); + private: + V vector_; + const std::vector &flatVector_; + std::vector dim_; + size_t ind_{0}; + unsigned int dimInd_{0}; + }; + + // Flatten class template implementation + template + void Flatten::accumulate(const Element &e) + { + flatVector_.push_back(e); + } + + template + template + void Flatten::accumulate(const W &v) + { + for (auto &e: v) + { + accumulate(e); + } + } + + template + void Flatten::accumulateDim(const Element &e) {}; + + template + template + void Flatten::accumulateDim(const W &v) + { + dim_.push_back(v.size()); + accumulateDim(v[0]); + } + + template + Flatten::Flatten(const V &vector) + : vector_(vector) + { + accumulate(vector_); + accumulateDim(vector_); + } + + template + const V & Flatten::getVector(void) + { + return vector_; + } + + template + const std::vector::Element> & + Flatten::getFlatVector(void) + { + return flatVector_; + } + + template + const std::vector & Flatten::getDim(void) + { + return dim_; + } + + // Reconstruct class template implementation + template + void Reconstruct::fill(std::vector &v) + { + for (auto &e: v) + { + e = flatVector_[ind_++]; + } + } + + template + template + void Reconstruct::fill(W &v) + { + for (auto &e: v) + { + fill(e); + } + } + + template + void Reconstruct::resize(std::vector &v, const unsigned int dim) + { + v.resize(dim_[dim]); + } + + template + template + void Reconstruct::resize(W &v, const unsigned int dim) + { + v.resize(dim_[dim]); + for (auto &e: v) + { + resize(e, dim + 1); + } + } + + template + Reconstruct::Reconstruct(const std::vector &flatVector, + const std::vector &dim) + : flatVector_(flatVector) + , dim_(dim) + { + resize(vector_, 0); + fill(vector_); + } + + template + const V & Reconstruct::getVector(void) + { + return vector_; + } + + template + const std::vector::Element> & + Reconstruct::getFlatVector(void) + { + return flatVector_; + } + + template + const std::vector & Reconstruct::getDim(void) + { + return dim_; + } + + // Vector IO utilities /////////////////////////////////////////////////////// + // helper function to read space-separated values + template + std::vector strToVec(const std::string s) + { + std::istringstream sstr(s); + T buf; + std::vector v; + + while(!sstr.eof()) + { + sstr >> buf; + v.push_back(buf); + } + + return v; + } + + // output to streams for vectors + template < class T > + inline std::ostream & operator<<(std::ostream &os, const std::vector &v) + { + os << "["; + for (auto &x: v) + { + os << x << " "; + } + if (v.size() > 0) + { + os << "\b"; + } + os << "]"; + + return os; + } +} + +#endif \ No newline at end of file diff --git a/lib/tensors/Tensor_logical.h b/lib/tensors/Tensor_logical.h index 7ab3668b..58b2b03b 100644 --- a/lib/tensors/Tensor_logical.h +++ b/lib/tensors/Tensor_logical.h @@ -55,5 +55,38 @@ LOGICAL_BINOP(&); LOGICAL_BINOP(||); LOGICAL_BINOP(&&); +template +strong_inline bool operator==(const iScalar &t1, const iScalar &t2) +{ + return (t1._internal == t2._internal); +} + +template +strong_inline bool operator==(const iVector &t1, const iVector &t2) +{ + bool res = true; + + for (unsigned int i = 0; i < N; ++i) + { + res = (res && (t1._internal[i] == t2._internal[i])); + } + + return res; +} + +template +strong_inline bool operator==(const iMatrix &t1, const iMatrix &t2) +{ + bool res = true; + + for (unsigned int i = 0; i < N; ++i) + for (unsigned int j = 0; j < N; ++j) + { + res = (res && (t1._internal[i][j] == t2._internal[i][j])); + } + + return res; +} + } #endif diff --git a/tests/IO/Test_serialisation.cc b/tests/IO/Test_serialisation.cc index 93007e44..bca4d01c 100644 --- a/tests/IO/Test_serialisation.cc +++ b/tests/IO/Test_serialisation.cc @@ -45,7 +45,8 @@ public: bool , b, std::vector, array, std::vector >, twodimarray, - std::vector > >, cmplx3darray + std::vector > >, cmplx3darray, + SpinColourMatrix, scm ); myclass() {} myclass(int i) @@ -59,6 +60,12 @@ public: y=2*i; b=true; name="bother said pooh"; + scm()(0, 1)(2, 1) = 2.356; + scm()(3, 0)(1, 1) = 1.323; + scm()(2, 1)(0, 1) = 5.3336; + scm()(0, 2)(1, 1) = 6.336; + scm()(2, 1)(2, 2) = 7.344; + scm()(1, 1)(2, 0) = 8.3534; } }; @@ -113,6 +120,10 @@ void tensorConvTestFn(GridSerialRNG &rng, const std::string label) int main(int argc,char **argv) { + GridSerialRNG rng; + + rng.SeedFixedIntegers(std::vector({42,10,81,9})); + std::cout << "==== basic IO" << std::endl; XmlWriter WR("bother.xml"); @@ -138,7 +149,7 @@ int main(int argc,char **argv) std::cout << "-- serialisable class writing to 'bother.xml'..." << std::endl; write(WR,"obj",obj); WR.write("obj2", obj); - vec.push_back(myclass(1234)); + vec.push_back(obj); vec.push_back(myclass(5678)); vec.push_back(myclass(3838)); pair = std::make_pair(myenum::red, myenum::blue); @@ -149,8 +160,6 @@ int main(int argc,char **argv) std::cout << "-- serialisable class comparison:" << std::endl; std::cout << "vec[0] == obj: " << ((vec[0] == obj) ? "true" : "false") << std::endl; std::cout << "vec[1] == obj: " << ((vec[1] == obj) ? "true" : "false") << std::endl; - - write(WR, "objpair", pair); std::cout << "-- pair writing to std::cout:" << std::endl; std::cout << pair << std::endl; @@ -159,26 +168,20 @@ int main(int argc,char **argv) //// XML ioTest("iotest.xml", obj, "XML (object) "); ioTest("iotest.xml", vec, "XML (vector of objects)"); - ioTest("iotest.xml", pair, "XML (pair of objects)"); //// binary ioTest("iotest.bin", obj, "binary (object) "); ioTest("iotest.bin", vec, "binary (vector of objects)"); - ioTest("iotest.bin", pair, "binary (pair of objects)"); //// text ioTest("iotest.dat", obj, "text (object) "); ioTest("iotest.dat", vec, "text (vector of objects)"); - ioTest("iotest.dat", pair, "text (pair of objects)"); //// text ioTest("iotest.json", obj, "JSON (object) "); ioTest("iotest.json", vec, "JSON (vector of objects)"); - ioTest("iotest.json", pair, "JSON (pair of objects)"); //// HDF5 -#undef HAVE_HDF5 #ifdef HAVE_HDF5 ioTest("iotest.h5", obj, "HDF5 (object) "); ioTest("iotest.h5", vec, "HDF5 (vector of objects)"); - ioTest("iotest.h5", pair, "HDF5 (pair of objects)"); #endif std::cout << "\n==== vector flattening/reconstruction" << std::endl; @@ -216,10 +219,6 @@ int main(int argc,char **argv) std::cout << std::endl; std::cout << "==== Grid tensor to vector test" << std::endl; - - GridSerialRNG rng; - - rng.SeedFixedIntegers(std::vector({42,10,81,9})); tensorConvTest(rng, SpinColourMatrix); tensorConvTest(rng, SpinColourVector); tensorConvTest(rng, ColourMatrix); From 2f849ee25227dde692d93c2d571051d18e8068f4 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 8 Mar 2018 23:34:00 +0000 Subject: [PATCH 12/29] declaration fix --- lib/serialisation/BaseIO.h | 48 --------------------------------- lib/serialisation/VectorUtils.h | 48 +++++++++++++++++++++++++++++++++ 2 files changed, 48 insertions(+), 48 deletions(-) diff --git a/lib/serialisation/BaseIO.h b/lib/serialisation/BaseIO.h index c9b3fb9e..298bff87 100644 --- a/lib/serialisation/BaseIO.h +++ b/lib/serialisation/BaseIO.h @@ -35,54 +35,6 @@ Author: Guido Cossu #include namespace Grid { - // Pair IO utilities ///////////////////////////////////////////////////////// - // helper function to parse input in the format "" - template - inline std::istream & operator>>(std::istream &is, std::pair &buf) - { - T1 buf1; - T2 buf2; - char c; - - // Search for "pair" delimiters. - do - { - is.get(c); - } while (c != '(' && !is.eof()); - if (c == '(') - { - int start = is.tellg(); - do - { - is.get(c); - } while (c != ')' && !is.eof()); - if (c == ')') - { - int end = is.tellg(); - int psize = end - start - 1; - - // Only read data between pair limiters. - is.seekg(start); - std::string tmpstr(psize, ' '); - is.read(&tmpstr[0], psize); - std::istringstream temp(tmpstr); - temp >> buf1 >> buf2; - buf = std::make_pair(buf1, buf2); - is.seekg(end); - } - } - is.peek(); - return is; - } - - // output to streams for pairs - template - inline std::ostream & operator<<(std::ostream &os, const std::pair &p) - { - os << "(" << p.first << " " << p.second << ")"; - return os; - } - // Abstract writer/reader classes //////////////////////////////////////////// // static polymorphism implemented using CRTP idiom class Serializable; diff --git a/lib/serialisation/VectorUtils.h b/lib/serialisation/VectorUtils.h index f5c76b84..6df9416d 100644 --- a/lib/serialisation/VectorUtils.h +++ b/lib/serialisation/VectorUtils.h @@ -5,6 +5,54 @@ #include namespace Grid { + // Pair IO utilities ///////////////////////////////////////////////////////// + // helper function to parse input in the format "" + template + inline std::istream & operator>>(std::istream &is, std::pair &buf) + { + T1 buf1; + T2 buf2; + char c; + + // Search for "pair" delimiters. + do + { + is.get(c); + } while (c != '(' && !is.eof()); + if (c == '(') + { + int start = is.tellg(); + do + { + is.get(c); + } while (c != ')' && !is.eof()); + if (c == ')') + { + int end = is.tellg(); + int psize = end - start - 1; + + // Only read data between pair limiters. + is.seekg(start); + std::string tmpstr(psize, ' '); + is.read(&tmpstr[0], psize); + std::istringstream temp(tmpstr); + temp >> buf1 >> buf2; + buf = std::make_pair(buf1, buf2); + is.seekg(end); + } + } + is.peek(); + return is; + } + + // output to streams for pairs + template + inline std::ostream & operator<<(std::ostream &os, const std::pair &p) + { + os << "(" << p.first << " " << p.second << ")"; + return os; + } + // Grid scalar tensors to nested std::vectors ////////////////////////////////// template struct TensorToVec From 70ec2faa98a86f904e1f63234e4ce1a602d939c5 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 9 Mar 2018 19:53:55 +0000 Subject: [PATCH 13/29] Hadrons: maximum iteration specified for tests and error if 0 --- extras/Hadrons/Modules/MSolver/RBPrecCG.hpp | 8 +++++++- tests/hadrons/Test_QED.cc | 5 +++-- tests/hadrons/Test_hadrons.hpp | 5 +++-- tests/hadrons/Test_hadrons_2AS_spectrum.cc | 5 +++-- tests/hadrons/Test_hadrons_meson_3pt.cc | 5 +++-- tests/hadrons/Test_hadrons_spectrum.cc | 5 +++-- tests/hadrons/Test_hadrons_wilsonFund.cc | 5 +++-- 7 files changed, 25 insertions(+), 13 deletions(-) diff --git a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp index 2b914625..4bff910e 100644 --- a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp +++ b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp @@ -111,9 +111,15 @@ std::vector TRBPrecCG::getOutput(void) template void TRBPrecCG::setup(void) { + if (par().maxIteration == 0) + { + HADRON_ERROR(Argument, "zero maximum iteration"); + } + LOG(Message) << "setting up Schur red-black preconditioned CG for" << " action '" << par().action << "' with residual " - << par().residual << std::endl; + << par().residual << ", maximum iteration " + << par().maxIteration << std::endl; auto Ls = env().getObjectLs(par().action); auto &mat = envGet(FMat, par().action); diff --git a/tests/hadrons/Test_QED.cc b/tests/hadrons/Test_QED.cc index 44dd0578..3377bf3c 100644 --- a/tests/hadrons/Test_QED.cc +++ b/tests/hadrons/Test_QED.cc @@ -93,8 +93,9 @@ int main(int argc, char *argv[]) // solvers MSolver::RBPrecCG::Par solverPar; - solverPar.action = "DWF_" + flavour[i]; - solverPar.residual = 1.0e-8; + solverPar.action = "DWF_" + flavour[i]; + solverPar.residual = 1.0e-8; + solverPar.maxIteration = 10000; application.createModule("CG_" + flavour[i], solverPar); diff --git a/tests/hadrons/Test_hadrons.hpp b/tests/hadrons/Test_hadrons.hpp index 6188a580..67124d6c 100644 --- a/tests/hadrons/Test_hadrons.hpp +++ b/tests/hadrons/Test_hadrons.hpp @@ -176,8 +176,9 @@ inline void makeRBPrecCGSolver(Application &application, std::string &solverName if (!(VirtualMachine::getInstance().hasModule(solverName))) { MSolver::RBPrecCG::Par solverPar; - solverPar.action = actionName; - solverPar.residual = residual; + solverPar.action = actionName; + solverPar.residual = residual; + solverPar.maxIteration = 10000; application.createModule(solverName, solverPar); } diff --git a/tests/hadrons/Test_hadrons_2AS_spectrum.cc b/tests/hadrons/Test_hadrons_2AS_spectrum.cc index 2f519834..b3906730 100644 --- a/tests/hadrons/Test_hadrons_2AS_spectrum.cc +++ b/tests/hadrons/Test_hadrons_2AS_spectrum.cc @@ -106,8 +106,9 @@ int main(int argc, char *argv[]) // solvers MSolver::RBPrecCG2AS::Par solverPar; - solverPar.action = "WilsonClover2AS_" + flavour[i]; - solverPar.residual = 1.0e-8; + solverPar.action = "WilsonClover2AS_" + flavour[i]; + solverPar.residual = 1.0e-8; + solverPar.maxIteration = 10000; application.createModule("CG_" + flavour[i], solverPar); diff --git a/tests/hadrons/Test_hadrons_meson_3pt.cc b/tests/hadrons/Test_hadrons_meson_3pt.cc index 382c39d4..b12ef472 100644 --- a/tests/hadrons/Test_hadrons_meson_3pt.cc +++ b/tests/hadrons/Test_hadrons_meson_3pt.cc @@ -82,8 +82,9 @@ int main(int argc, char *argv[]) // solvers MSolver::RBPrecCG::Par solverPar; - solverPar.action = "DWF_" + flavour[i]; - solverPar.residual = 1.0e-8; + solverPar.action = "DWF_" + flavour[i]; + solverPar.residual = 1.0e-8; + solverPar.maxIteration = 10000; application.createModule("CG_" + flavour[i], solverPar); } diff --git a/tests/hadrons/Test_hadrons_spectrum.cc b/tests/hadrons/Test_hadrons_spectrum.cc index 801674f7..682dcc9b 100644 --- a/tests/hadrons/Test_hadrons_spectrum.cc +++ b/tests/hadrons/Test_hadrons_spectrum.cc @@ -84,8 +84,9 @@ int main(int argc, char *argv[]) // solvers MSolver::RBPrecCG::Par solverPar; - solverPar.action = "DWF_" + flavour[i]; - solverPar.residual = 1.0e-8; + solverPar.action = "DWF_" + flavour[i]; + solverPar.residual = 1.0e-8; + solverPar.maxIteration = 10000; application.createModule("CG_" + flavour[i], solverPar); diff --git a/tests/hadrons/Test_hadrons_wilsonFund.cc b/tests/hadrons/Test_hadrons_wilsonFund.cc index 083d3b8c..36e751b6 100644 --- a/tests/hadrons/Test_hadrons_wilsonFund.cc +++ b/tests/hadrons/Test_hadrons_wilsonFund.cc @@ -96,8 +96,9 @@ int main(int argc, char *argv[]) // solvers MSolver::RBPrecCG::Par solverPar; - solverPar.action = "WilsonClover_" + flavour[i]; - solverPar.residual = 1.0e-8; + solverPar.action = "WilsonClover_" + flavour[i]; + solverPar.residual = 1.0e-8; + solverPar.maxIteration = 10000; application.createModule("CG_" + flavour[i], solverPar); From b801e1fcd637b402e381de98261386f11d0d8da4 Mon Sep 17 00:00:00 2001 From: paboyle Date: Fri, 9 Mar 2018 20:44:10 +0000 Subject: [PATCH 14/29] fclose should be called through a call to close() --- lib/parallelIO/IldgIO.h | 1 - 1 file changed, 1 deletion(-) diff --git a/lib/parallelIO/IldgIO.h b/lib/parallelIO/IldgIO.h index b86e250f..b0bd7e2c 100644 --- a/lib/parallelIO/IldgIO.h +++ b/lib/parallelIO/IldgIO.h @@ -568,7 +568,6 @@ class IldgWriter : public ScidacWriter { writeLimeIldgLFN(header.ildg_lfn); // rec writeLimeLatticeBinaryObject(Umu,std::string(ILDG_BINARY_DATA)); // Closes message with checksum // limeDestroyWriter(LimeW); - fclose(File); } }; From e485a07133f0e8694f3e51e8167c84848f673060 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 9 Mar 2018 21:56:01 +0000 Subject: [PATCH 15/29] Hadrons: garbage collector debug output --- extras/Hadrons/VirtualMachine.cc | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/extras/Hadrons/VirtualMachine.cc b/extras/Hadrons/VirtualMachine.cc index eb1f68ba..8d253805 100644 --- a/extras/Hadrons/VirtualMachine.cc +++ b/extras/Hadrons/VirtualMachine.cc @@ -632,6 +632,15 @@ void VirtualMachine::executeProgram(const Program &p) const // build garbage collection schedule LOG(Debug) << "Building garbage collection schedule..." << std::endl; freeProg = makeGarbageSchedule(p); + for (unsigned int i = 0; i < freeProg.size(); ++i) + { + LOG(Debug) << std::setw(4) << i + 1 << ": ["; + for (auto &a: freeProg[i]) + { + std::cout << env().getObjectName(a) << " "; + } + std::cout << "]" << std::endl; + } // program execution LOG(Debug) << "Executing program..." << std::endl; From 229977c9553a4b6f9757fddec33dfe8b43b2cd7d Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 9 Mar 2018 21:56:27 +0000 Subject: [PATCH 16/29] Hadrons: minor memory fix for ShiftProbe module --- extras/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/extras/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp b/extras/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp index 89059180..8d52327e 100644 --- a/extras/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp +++ b/extras/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp @@ -116,7 +116,7 @@ template void TShiftProbe::setup(void) { envTmpLat(Field, "acc"); - envCacheLat(ComplexField, getName()); + envCreateLat(ComplexField, getName()); } // execution /////////////////////////////////////////////////////////////////// From f57afe2079aad747b7ccd2d8f1835f2b58cb523b Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 13 Mar 2018 13:51:09 +0000 Subject: [PATCH 17/29] Hadrons: much cleaner eigenpack implementation, to be tested --- extras/Hadrons/EigenPack.hpp | 213 ++++++++++++++++++ extras/Hadrons/Environment.cc | 99 +++++++- extras/Hadrons/Environment.hpp | 10 + extras/Hadrons/LanczosUtils.hpp | 115 ---------- extras/Hadrons/Makefile.am | 2 +- .../Modules/MSolver/LocalCoherenceLanczos.hpp | 142 ++++-------- 6 files changed, 362 insertions(+), 219 deletions(-) create mode 100644 extras/Hadrons/EigenPack.hpp delete mode 100644 extras/Hadrons/LanczosUtils.hpp diff --git a/extras/Hadrons/EigenPack.hpp b/extras/Hadrons/EigenPack.hpp new file mode 100644 index 00000000..eeb04d17 --- /dev/null +++ b/extras/Hadrons/EigenPack.hpp @@ -0,0 +1,213 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/EigenPack.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_EigenPack_hpp_ +#define Hadrons_EigenPack_hpp_ + +#include +#include + +BEGIN_HADRONS_NAMESPACE + +// Lanczos type +#ifndef HADRONS_DEFAULT_LANCZOS_NBASIS +#define HADRONS_DEFAULT_LANCZOS_NBASIS 60 +#endif + +template +class EigenPack +{ +public: + std::vector eval; + std::vector evec; +public: + EigenPack(void) = default; + virtual ~EigenPack(void) = default; + + EigenPack(const size_t size, GridBase *grid) + { + resize(size, grid); + } + + void resize(const size_t size, GridBase *grid) + { + eval.resize(size); + evec.resize(size, grid); + } + + virtual void read(const std::string fileStem, const int traj = -1) + { + std::string evecFilename, evalFilename; + + makeFilenames(evecFilename, evalFilename, fileStem, traj); + XmlReader xmlReader(evalFilename); + basicRead(evec, evecFilename, evec.size()); + LOG(Message) << "Reading " << eval.size() << " eigenvalues from '" + << evalFilename << "'" << std::endl; + Grid::read(xmlReader, "evals", eval); + } + + virtual void write(const std::string fileStem, const int traj = -1) + { + std::string evecFilename, evalFilename; + + makeFilenames(evecFilename, evalFilename, fileStem, traj); + XmlWriter xmlWriter(evalFilename); + basicWrite(evecFilename, evec, evec.size()); + LOG(Message) << "Writing " << eval.size() << " eigenvalues to '" + << evalFilename << "'" << std::endl; + Grid::write(xmlWriter, "evals", eval); + } +protected: + void makeFilenames(std::string &evecFilename, std::string &evalFilename, + const std::string stem, const int traj = -1) + { + std::string t = (traj < 0) ? "" : ("." + std::to_string(traj)); + + evecFilename = stem + "_evec" + t + ".bin"; + evalFilename = stem + "_eval" + t + ".xml"; + } + + template + static void basicRead(std::vector &evec, const std::string filename, + const unsigned int size) + { + emptyUserRecord record; + ScidacReader binReader; + + binReader.open(filename); + for(int k = 0; k < size; ++k) + { + binReader.readScidacFieldRecord(evec[k], record); + } + binReader.close(); + } + + template + static void basicWrite(const std::string filename, std::vector &evec, + const unsigned int size) + { + emptyUserRecord record; + ScidacWriter binWriter; + + binWriter.open(filename); + for(int k = 0; k < size; ++k) + { + binWriter.writeScidacFieldRecord(evec[k], record); + } + binWriter.close(); + } +}; + +template +class CoarseEigenPack: public EigenPack +{ +public: + std::vector evalCoarse; + std::vector evecCoarse; +public: + CoarseEigenPack(void) = default; + virtual ~CoarseEigenPack(void) = default; + + CoarseEigenPack(const size_t sizeFine, const size_t sizeCoarse, + GridBase *gridFine, GridBase *gridCoarse) + { + resize(sizeFine, sizeCoarse, gridFine, gridCoarse); + } + + void resize(const size_t sizeFine, const size_t sizeCoarse, + GridBase *gridFine, GridBase *gridCoarse) + { + EigenPack::resize(sizeFine, gridFine); + evalCoarse.resize(sizeCoarse); + evecCoarse.resize(sizeCoarse, gridCoarse); + } + + virtual void read(const std::string fileStem, const int traj = -1) + { + std::string evecFineFilename, evalFineFilename; + std::string evecCoarseFilename, evalCoarseFilename; + + this->makeFilenames(evecFineFilename, evalFineFilename, + fileStem + "_fine", traj); + this->makeFilenames(evecCoarseFilename, evalCoarseFilename, + fileStem + "_coarse", traj); + XmlReader xmlFineReader(evalFineFilename); + XmlReader xmlCoarseReader(evalCoarseFilename); + LOG(Message) << "Reading " << this->evec.size() << " fine eigenvectors from '" + << evecFineFilename << "'" << std::endl; + this->basicRead(this->evec, evecFineFilename, this->evec.size()); + LOG(Message) << "Reading " << evecCoarse.size() << " coarse eigenvectors from '" + << evecCoarseFilename << "'" << std::endl; + this->basicRead(evecCoarse, evecCoarseFilename, evecCoarse.size()); + LOG(Message) << "Reading " << this->eval.size() << " fine eigenvalues from '" + << evalFineFilename << "'" << std::endl; + Grid::read(xmlFineReader, "evals", this->eval); + LOG(Message) << "Reading " << evalCoarse.size() << " coarse eigenvalues from '" + << evalCoarseFilename << "'" << std::endl; + Grid::read(xmlCoarseReader, "evals", evalCoarse); + } + + virtual void write(const std::string fileStem, const int traj = -1) + { + std::string evecFineFilename, evalFineFilename; + std::string evecCoarseFilename, evalCoarseFilename; + + this->makeFilenames(evecFineFilename, evalFineFilename, + fileStem + "_fine", traj); + this->makeFilenames(evecCoarseFilename, evalCoarseFilename, + fileStem + "_coarse", traj); + XmlWriter xmlFineWriter(evalFineFilename); + XmlWriter xmlCoarseWriter(evalCoarseFilename); + LOG(Message) << "Writing " << this->evec.size() << " fine eigenvectors to '" + << evecFineFilename << "'" << std::endl; + this->basicWrite(evecFineFilename, this->evec, this->evec.size()); + LOG(Message) << "Writing " << evecCoarse.size() << " coarse eigenvectors to '" + << evecCoarseFilename << "'" << std::endl; + this->basicWrite(evecCoarseFilename, evecCoarse, evecCoarse.size()); + LOG(Message) << "Writing " << this->eval.size() << " fine eigenvalues to '" + << evalFineFilename << "'" << std::endl; + Grid::write(xmlFineWriter, "evals", this->eval); + LOG(Message) << "Writing " << evalCoarse.size() << " coarse eigenvalues to '" + << evalCoarseFilename << "'" << std::endl; + Grid::write(xmlCoarseWriter, "evals", evalCoarse); + } +}; + +template +using FermionEigenPack = EigenPack; + +template +using CoarseFermionEigenPack = CoarseEigenPack< + typename FImpl::FermionField, + typename LocalCoherenceLanczos::CoarseField>; + +END_HADRONS_NAMESPACE + +#endif // Hadrons_EigenPack_hpp_ diff --git a/extras/Hadrons/Environment.cc b/extras/Hadrons/Environment.cc index 6554122e..9c9618f7 100644 --- a/extras/Hadrons/Environment.cc +++ b/extras/Hadrons/Environment.cc @@ -61,7 +61,7 @@ Environment::Environment(void) // grids /////////////////////////////////////////////////////////////////////// void Environment::createGrid(const unsigned int Ls) { - if (grid5d_.find(Ls) == grid5d_.end()) + if ((Ls > 1) and (grid5d_.find(Ls) == grid5d_.end())) { auto g = getGrid(); @@ -70,6 +70,53 @@ void Environment::createGrid(const unsigned int Ls) } } +void Environment::createCoarseGrid(const std::vector &blockSize, + const unsigned int Ls) +{ + int nd = getNd(); + std::vector fineDim = getDim(), coarseDim; + unsigned int cLs; + auto key4d = blockSize, key5d = blockSize; + + createGrid(Ls); + coarseDim.resize(nd); + for (int d = 0; d < coarseDim.size(); d++) + { + coarseDim[d] = fineDim[d]/blockSize[d]; + if (coarseDim[d]*blockSize[d] != fineDim[d]) + { + HADRON_ERROR(Size, "Fine dimension " + std::to_string(d) + + " (" + std::to_string(fineDim[d]) + + ") not divisible by coarse dimension (" + + std::to_string(coarseDim[d]) + ")"); + } + } + if (blockSize.size() > nd) + { + cLs = Ls/blockSize[nd]; + if (cLs*blockSize[nd] != Ls) + { + HADRON_ERROR(Size, "Fine Ls (" + std::to_string(Ls) + + ") not divisible by coarse Ls (" + + std::to_string(cLs) + ")"); + } + key4d.resize(nd); + key5d.push_back(Ls); + } + gridCoarse4d_[key4d].reset( + SpaceTimeGrid::makeFourDimGrid(coarseDim, + GridDefaultSimd(nd, vComplex::Nsimd()), GridDefaultMpi())); + gridCoarseRb4d_[key4d].reset( + SpaceTimeGrid::makeFourDimRedBlackGrid(gridCoarse4d_[key4d].get())); + if (Ls > 1) + { + gridCoarse5d_[key5d].reset( + SpaceTimeGrid::makeFiveDimGrid(cLs, gridCoarse4d_[key4d].get())); + gridCoarseRb5d_[key5d].reset( + SpaceTimeGrid::makeFiveDimRedBlackGrid(cLs, gridCoarse4d_[key4d].get())); + } +} + GridCartesian * Environment::getGrid(const unsigned int Ls) const { try @@ -104,7 +151,55 @@ GridRedBlackCartesian * Environment::getRbGrid(const unsigned int Ls) const } catch(std::out_of_range &) { - HADRON_ERROR(Definition, "no red-black 5D grid with Ls= " + std::to_string(Ls)); + HADRON_ERROR(Definition, "no red-black grid with Ls= " + std::to_string(Ls)); + } +} + +GridCartesian * Environment::getCoarseGrid( + const std::vector &blockSize, const unsigned int Ls) const +{ + auto key = blockSize; + + try + { + if (Ls == 1) + { + key.resize(getNd()); + return gridCoarse4d_.at(key).get(); + } + else + { + key.push_back(Ls); + return gridCoarse5d_.at(key).get(); + } + } + catch(std::out_of_range &) + { + HADRON_ERROR(Definition, "no coarse grid with Ls= " + std::to_string(Ls)); + } +} + +GridRedBlackCartesian * Environment::getRbCoarseGrid( + const std::vector &blockSize, const unsigned int Ls) const +{ + auto key = blockSize; + + try + { + if (Ls == 1) + { + key.resize(getNd()); + return gridCoarseRb4d_.at(key).get(); + } + else + { + key.push_back(Ls); + return gridCoarseRb5d_.at(key).get(); + } + } + catch(std::out_of_range &) + { + HADRON_ERROR(Definition, "no coarse red-black grid with Ls= " + std::to_string(Ls)); } } diff --git a/extras/Hadrons/Environment.hpp b/extras/Hadrons/Environment.hpp index e9bfffe1..3b1d45f8 100644 --- a/extras/Hadrons/Environment.hpp +++ b/extras/Hadrons/Environment.hpp @@ -86,8 +86,14 @@ private: public: // grids void createGrid(const unsigned int Ls); + void createCoarseGrid(const std::vector &blockSize, + const unsigned int Ls = 1); GridCartesian * getGrid(const unsigned int Ls = 1) const; GridRedBlackCartesian * getRbGrid(const unsigned int Ls = 1) const; + GridCartesian * getCoarseGrid(const std::vector &blockSize, + const unsigned int Ls = 1) const; + GridRedBlackCartesian * getRbCoarseGrid(const std::vector &blockSize, + const unsigned int Ls = 1) const; std::vector getDim(void) const; int getDim(const unsigned int mu) const; unsigned long int getLocalVolume(void) const; @@ -155,6 +161,10 @@ private: std::map grid5d_; GridRbPt gridRb4d_; std::map gridRb5d_; + std::map, GridPt> gridCoarse4d_; + std::map, GridRbPt> gridCoarseRb4d_; + std::map, GridPt> gridCoarse5d_; + std::map, GridRbPt> gridCoarseRb5d_; unsigned int nd_; // random number generator RngPt rng4d_; diff --git a/extras/Hadrons/LanczosUtils.hpp b/extras/Hadrons/LanczosUtils.hpp deleted file mode 100644 index a080da4b..00000000 --- a/extras/Hadrons/LanczosUtils.hpp +++ /dev/null @@ -1,115 +0,0 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: extras/Hadrons/LanczosUtils.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_LanczosUtils_hpp_ -#define Hadrons_LanczosUtils_hpp_ - -#include -#include - -BEGIN_HADRONS_NAMESPACE - -// Lanczos type -#ifndef HADRONS_DEFAULT_LANCZOS_NBASIS -#define HADRONS_DEFAULT_LANCZOS_NBASIS 60 -#endif - -template -struct EigenPack -{ - typedef T VectorType; - std::vector eval; - std::vector evec; - - EigenPack(void) = default; - - EigenPack(const size_t size, GridBase *grid) - { - resize(size, grid); - } - - void resize(const size_t size, GridBase *grid) - { - eval.resize(size); - evec.resize(size, grid); - } - - void read(const std::string fileStem) - { - std::string evecFilename = fileStem + "_evec.bin"; - std::string evalFilename = fileStem + "_eval.xml"; - emptyUserRecord record; - ScidacReader binReader; - XmlReader xmlReader(evalFilename); - - LOG(Message) << "Reading " << evec.size() << " eigenvectors from '" - << evecFilename << "'" << std::endl; - binReader.open(evecFilename); - for(int k = 0; k < evec.size(); ++k) - { - binReader.readScidacFieldRecord(evec[k], record); - } - binReader.close(); - LOG(Message) << "Reading " << eval.size() << " eigenvalues from '" - << evalFilename << "'" << std::endl; - Grid::read(xmlReader, "evals", eval); - } - - void write(const std::string fileStem) - { - std::string evecFilename = fileStem + "_evec.bin"; - std::string evalFilename = fileStem + "_eval.xml"; - emptyUserRecord record; - ScidacWriter binWriter; - XmlWriter xmlWriter(evalFilename); - - LOG(Message) << "Writing " << evec.size() << " eigenvectors to '" - << evecFilename << "'" << std::endl; - binWriter.open(fileStem + "_evec.bin"); - for(int k = 0; k < evec.size(); ++k) - { - binWriter.writeScidacFieldRecord(evec[k], record); - } - binWriter.close(); - LOG(Message) << "Writing " << eval.size() << " eigenvalues to '" - << evalFilename << "'" << std::endl; - Grid::write(xmlWriter, "evals", eval); - } -}; - -template -using FineEigenPack = EigenPack; - -template -using CoarseEigenPack = EigenPack< - typename LocalCoherenceLanczos::CoarseField>; - -END_HADRONS_NAMESPACE - -#endif // Hadrons_LanczosUtils_hpp_ diff --git a/extras/Hadrons/Makefile.am b/extras/Hadrons/Makefile.am index 3b393cd1..3b945124 100644 --- a/extras/Hadrons/Makefile.am +++ b/extras/Hadrons/Makefile.am @@ -21,7 +21,7 @@ nobase_libHadrons_a_HEADERS = \ GeneticScheduler.hpp \ Global.hpp \ Graph.hpp \ - LanczosUtils.hpp \ + EigenPack.hpp \ Module.hpp \ Modules.hpp \ ModuleFactory.hpp \ diff --git a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp index 6e2103ce..3cc17bdc 100644 --- a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp +++ b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp @@ -31,7 +31,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include -#include +#include BEGIN_HADRONS_NAMESPACE @@ -45,8 +45,7 @@ class LocalCoherenceLanczosPar: Serializable public: GRID_SERIALIZABLE_CLASS_MEMBERS(LocalCoherenceLanczosPar, std::string, action, - int, doFine, - int, doCoarse, + bool, doCoarse, LanczosParams, fineParams, LanczosParams, coarseParams, ChebyParams, smoother, @@ -63,8 +62,8 @@ public: typedef LocalCoherenceLanczos LCL; - typedef FineEigenPack FinePack; - typedef CoarseEigenPack CoarsePack; + typedef FermionEigenPack BasePack; + typedef CoarseFermionEigenPack CoarsePack; typedef HADRONS_DEFAULT_SCHUR_OP SchurFMat; public: // constructor @@ -79,15 +78,7 @@ public: // execution virtual void execute(void); private: - void makeCoarseGrid(void); -private: - std::vector coarseDim_; - int Ls_, cLs_{1}; - std::unique_ptr coarseGrid4_{nullptr}; - std::unique_ptr coarseGrid_{nullptr}; - std::unique_ptr coarseGrid4Rb_{nullptr}; - std::unique_ptr coarseGridRb_{nullptr}; - std::string fineName_, coarseName_; + std::string fineName_, coarseName_; }; MODULE_REGISTER_NS(LocalCoherenceLanczos, @@ -127,55 +118,6 @@ std::vector TLocalCoherenceLanczos::getOutput(void) } // setup /////////////////////////////////////////////////////////////////////// -template -void TLocalCoherenceLanczos::makeCoarseGrid(void) -{ - int nd = env().getNd(); - std::vector blockSize = strToVec(par().blockSize); - auto fineDim = env().getDim(); - - Ls_ = env().getObjectLs(par().action); - env().createGrid(Ls_); - coarseDim_.resize(nd); - for (int d = 0; d < coarseDim_.size(); d++) - { - coarseDim_[d] = fineDim[d]/blockSize[d]; - if (coarseDim_[d]*blockSize[d] != fineDim[d]) - { - HADRON_ERROR(Size, "Fine dimension " + std::to_string(d) - + " (" + std::to_string(fineDim[d]) - + ") not divisible by coarse dimension (" - + std::to_string(coarseDim_[d]) + ")"); - } - } - if (blockSize.size() > nd) - { - cLs_ = Ls_/blockSize[nd]; - if (cLs_*blockSize[nd] != Ls_) - { - HADRON_ERROR(Size, "Fine Ls (" + std::to_string(Ls_) - + ") not divisible by coarse Ls (" - + std::to_string(cLs_) + ")"); - } - } - if (Ls_ > 1) - { - coarseGrid4_.reset(SpaceTimeGrid::makeFourDimGrid( - coarseDim_, GridDefaultSimd(nd, vComplex::Nsimd()), - GridDefaultMpi())); - coarseGrid4Rb_.reset(SpaceTimeGrid::makeFourDimRedBlackGrid(coarseGrid4_.get())); - coarseGrid_.reset(SpaceTimeGrid::makeFiveDimGrid(cLs_, coarseGrid4_.get())); - coarseGridRb_.reset(SpaceTimeGrid::makeFiveDimRedBlackGrid(cLs_, coarseGrid4_.get())); - } - else - { - coarseGrid_.reset(SpaceTimeGrid::makeFourDimGrid( - coarseDim_, GridDefaultSimd(nd, vComplex::Nsimd()), - GridDefaultMpi())); - coarseGridRb_.reset(SpaceTimeGrid::makeFourDimRedBlackGrid(coarseGrid_.get())); - } -} - template void TLocalCoherenceLanczos::setup(void) { @@ -183,19 +125,25 @@ void TLocalCoherenceLanczos::setup(void) << " action '" << par().action << "' (" << nBasis << " eigenvectors)..." << std::endl; - if (!coarseGrid_) - { - makeCoarseGrid(); - } - LOG(Message) << "Coarse grid: " << coarseGrid_->GlobalDimensions() << std::endl; - envCreate(FinePack, fineName_, Ls_, par().fineParams.Nm, env().getRbGrid(Ls_)); - envCreate(CoarsePack, coarseName_, Ls_, par().coarseParams.Nm, coarseGridRb_.get()); - auto &fine = envGet(FinePack, fineName_); - auto &coarse = envGet(CoarsePack, coarseName_); - envTmp(SchurFMat, "mat", Ls_, envGet(FMat, par().action)); + unsigned int Ls = env().getObjectLs(par().action); + auto blockSize = strToVec(par().blockSize); + + env().createCoarseGrid(blockSize, Ls); + + auto cg = env().getCoarseGrid(blockSize, Ls); + auto cgrb = env().getRbCoarseGrid(blockSize, Ls); + int cNm = (par().doCoarse) ? par().coarseParams.Nm : 0; + + LOG(Message) << "Coarse grid: " << cg->GlobalDimensions() << std::endl; + envCreateDerived(BasePack, CoarsePack, getName(), Ls, + par().fineParams.Nm, cNm, env().getRbGrid(Ls), cgrb); + + auto &epack = envGet(CoarsePack, getName()); + + envTmp(SchurFMat, "mat", Ls, envGet(FMat, par().action)); envGetTmp(SchurFMat, mat); - envTmp(LCL, "solver", Ls_, env().getRbGrid(Ls_), coarseGridRb_.get(), mat, - Odd, fine.evec, coarse.evec, fine.eval, coarse.eval); + envTmp(LCL, "solver", Ls, env().getRbGrid(Ls), cgrb, mat, + Odd, epack.evec, epack.evecCoarse, epack.eval, epack.evalCoarse); } // execution /////////////////////////////////////////////////////////////////// @@ -204,41 +152,33 @@ void TLocalCoherenceLanczos::execute(void) { auto &finePar = par().fineParams; auto &coarsePar = par().coarseParams; - auto &fine = envGet(FinePack, fineName_); - auto &coarse = envGet(CoarsePack, coarseName_); + auto &epack = envGet(CoarsePack, getName()); envGetTmp(LCL, solver); - if (par().doFine) - { - LOG(Message) << "Performing fine grid IRL -- Nstop= " - << finePar.Nstop << ", Nk= " << finePar.Nk << ", Nm= " - << finePar.Nm << std::endl; - solver.calcFine(finePar.Cheby, finePar.Nstop, finePar.Nk, finePar.Nm, - finePar.resid,finePar.MaxIt, finePar.betastp, - finePar.MinRes); - solver.testFine(finePar.resid*100.0); - LOG(Message) << "Orthogonalising" << std::endl; - solver.Orthogonalise(); - if (!par().output.empty()) - { - fine.write(par().output + "_fine"); - } - } + LOG(Message) << "Performing fine grid IRL -- Nstop= " + << finePar.Nstop << ", Nk= " << finePar.Nk << ", Nm= " + << finePar.Nm << std::endl; + solver.calcFine(finePar.Cheby, finePar.Nstop, finePar.Nk, finePar.Nm, + finePar.resid,finePar.MaxIt, finePar.betastp, + finePar.MinRes); + solver.testFine(finePar.resid*100.0); if (par().doCoarse) { + LOG(Message) << "Orthogonalising" << std::endl; + solver.Orthogonalise(); LOG(Message) << "Performing coarse grid IRL -- Nstop= " - << coarsePar.Nstop << ", Nk= " << coarsePar.Nk << ", Nm= " - << coarsePar.Nm << std::endl; + << coarsePar.Nstop << ", Nk= " << coarsePar.Nk << ", Nm= " + << coarsePar.Nm << std::endl; solver.calcCoarse(coarsePar.Cheby, par().smoother, par().coarseRelaxTol, - coarsePar.Nstop, coarsePar.Nk, coarsePar.Nm, + coarsePar.Nstop, coarsePar.Nk, coarsePar.Nm, coarsePar.resid, coarsePar.MaxIt, coarsePar.betastp, coarsePar.MinRes); solver.testCoarse(coarsePar.resid*100.0, par().smoother, - par().coarseRelaxTol); - if (!par().output.empty()) - { - coarse.write(par().output + "_coarse"); - } + par().coarseRelaxTol); + } + if (!par().output.empty()) + { + epack.write(par().output, vm().getTrajectory()); } } From b85f987b0b29756d320ee17056d15d1e39d5cf44 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 13 Mar 2018 16:09:22 +0000 Subject: [PATCH 18/29] Hadrons: error message channel verbose during profiling --- extras/Hadrons/VirtualMachine.cc | 2 -- 1 file changed, 2 deletions(-) diff --git a/extras/Hadrons/VirtualMachine.cc b/extras/Hadrons/VirtualMachine.cc index 8d253805..ba963276 100644 --- a/extras/Hadrons/VirtualMachine.cc +++ b/extras/Hadrons/VirtualMachine.cc @@ -381,7 +381,6 @@ void VirtualMachine::makeMemoryProfile(void) env().protectObjects(false); GridLogMessage.Active(false); HadronsLogMessage.Active(false); - HadronsLogError.Active(false); for (auto it = program.rbegin(); it != program.rend(); ++it) { auto a = *it; @@ -397,7 +396,6 @@ void VirtualMachine::makeMemoryProfile(void) env().protectObjects(protect); GridLogMessage.Active(gmsg); HadronsLogMessage.Active(hmsg); - HadronsLogError.Active(err); LOG(Debug) << "Memory profile:" << std::endl; LOG(Debug) << "----------------" << std::endl; for (unsigned int a = 0; a < profile_.module.size(); ++a) From 78f8d475283269f43b7d26b5711edf8459a19bb5 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 13 Mar 2018 16:10:16 +0000 Subject: [PATCH 19/29] Hadrons: environment access to derived objects --- extras/Hadrons/Environment.hpp | 48 ++++++++++++++++++++++++++++------ extras/Hadrons/Module.hpp | 3 +++ 2 files changed, 43 insertions(+), 8 deletions(-) diff --git a/extras/Hadrons/Environment.hpp b/extras/Hadrons/Environment.hpp index 3b1d45f8..637962b1 100644 --- a/extras/Hadrons/Environment.hpp +++ b/extras/Hadrons/Environment.hpp @@ -116,6 +116,10 @@ public: Ts && ... args); void setObjectModule(const unsigned int objAddress, const int modAddress); + template + T * getDerivedObject(const unsigned int address) const; + template + T * getDerivedObject(const std::string name) const; template T * getObject(const unsigned int address) const; template @@ -186,7 +190,7 @@ Holder::Holder(T *pt) template T & Holder::get(void) const { - return &objPt_.get(); + return *objPt_.get(); } template @@ -231,7 +235,7 @@ void Environment::createDerivedObject(const std::string name, object_[address].Ls = Ls; object_[address].data.reset(new Holder(new T(std::forward(args)...))); object_[address].size = MemoryProfiler::stats->maxAllocated - initMem; - object_[address].type = &typeid(T); + object_[address].type = &typeid(B); if (MemoryProfiler::stats == &memStats) { MemoryProfiler::stats = nullptr; @@ -241,7 +245,7 @@ void Environment::createDerivedObject(const std::string name, else if ((object_[address].storage != Storage::cache) or (object_[address].storage != storage) or (object_[address].name != name) or - (object_[address].type != &typeid(T))) + (object_[address].type != &typeid(B))) { HADRON_ERROR(Definition, "object '" + name + "' already allocated"); } @@ -256,21 +260,37 @@ void Environment::createObject(const std::string name, createDerivedObject(name, storage, Ls, std::forward(args)...); } -template -T * Environment::getObject(const unsigned int address) const +template +T * Environment::getDerivedObject(const unsigned int address) const { if (hasObject(address)) { if (hasCreatedObject(address)) { - if (auto h = dynamic_cast *>(object_[address].data.get())) + if (auto h = dynamic_cast *>(object_[address].data.get())) { - return h->getPt(); + if (&typeid(T) == &typeid(B)) + { + return dynamic_cast(h->getPt()); + } + else + { + if (auto hder = dynamic_cast(h->getPt())) + { + return hder; + } + else + { + HADRON_ERROR(Definition, "object with address " + std::to_string(address) + + " cannot be casted to '" + typeName(&typeid(T)) + + "' (has type '" + typeName(&typeid(h->get())) + "')"); + } + } } else { HADRON_ERROR(Definition, "object with address " + std::to_string(address) + - " does not have type '" + typeName(&typeid(T)) + + " does not have type '" + typeName(&typeid(B)) + "' (has type '" + getObjectType(address) + "')"); } } @@ -286,6 +306,18 @@ T * Environment::getObject(const unsigned int address) const } } +template +T * Environment::getDerivedObject(const std::string name) const +{ + return getDerivedObject(getObjectAddress(name)); +} + +template +T * Environment::getObject(const unsigned int address) const +{ + return getDerivedObject(address); +} + template T * Environment::getObject(const std::string name) const { diff --git a/extras/Hadrons/Module.hpp b/extras/Hadrons/Module.hpp index 018a26f7..85c27472 100644 --- a/extras/Hadrons/Module.hpp +++ b/extras/Hadrons/Module.hpp @@ -91,6 +91,9 @@ static ns##mod##ModuleRegistrar ns##mod##ModuleRegistrarInstance; #define envGet(type, name)\ *env().template getObject(name) +#define envGetDerived(base, type, name)\ +*env().template getDerivedObject(name) + #define envGetTmp(type, var)\ type &var = *env().template getObject(getName() + "_tmp_" + #var) From 2d4d70d3ecda9cf817f36be7649796489ed35c7b Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 13 Mar 2018 16:10:36 +0000 Subject: [PATCH 20/29] Hadrons: LCL fixes --- .../Modules/MSolver/LocalCoherenceLanczos.hpp | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp index 3cc17bdc..80717f32 100644 --- a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp +++ b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp @@ -77,8 +77,6 @@ public: virtual void setup(void); // execution virtual void execute(void); -private: - std::string fineName_, coarseName_; }; MODULE_REGISTER_NS(LocalCoherenceLanczos, @@ -95,10 +93,7 @@ MODULE_REGISTER_NS(ZLocalCoherenceLanczos, template TLocalCoherenceLanczos::TLocalCoherenceLanczos(const std::string name) : Module(name) -{ - fineName_ = getName() + "_fine"; - coarseName_ = getName() + "_coarse"; -} +{} // dependencies/products /////////////////////////////////////////////////////// template @@ -112,7 +107,7 @@ std::vector TLocalCoherenceLanczos::getInput(void) template std::vector TLocalCoherenceLanczos::getOutput(void) { - std::vector out = {fineName_, coarseName_}; + std::vector out = {getName()}; return out; } @@ -138,7 +133,7 @@ void TLocalCoherenceLanczos::setup(void) envCreateDerived(BasePack, CoarsePack, getName(), Ls, par().fineParams.Nm, cNm, env().getRbGrid(Ls), cgrb); - auto &epack = envGet(CoarsePack, getName()); + auto &epack = envGetDerived(BasePack, CoarsePack, getName()); envTmp(SchurFMat, "mat", Ls, envGet(FMat, par().action)); envGetTmp(SchurFMat, mat); @@ -152,7 +147,7 @@ void TLocalCoherenceLanczos::execute(void) { auto &finePar = par().fineParams; auto &coarsePar = par().coarseParams; - auto &epack = envGet(CoarsePack, getName()); + auto &epack = envGetDerived(BasePack, CoarsePack, getName()); envGetTmp(LCL, solver); LOG(Message) << "Performing fine grid IRL -- Nstop= " From 72344d1418ee165c0acf9c6cf4b3a373efadc085 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 13 Mar 2018 17:10:54 +0000 Subject: [PATCH 21/29] Hadrons: change default Schur convention to DiagTwo --- extras/Hadrons/Global.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/extras/Hadrons/Global.hpp b/extras/Hadrons/Global.hpp index e9f5933b..12b9a029 100644 --- a/extras/Hadrons/Global.hpp +++ b/extras/Hadrons/Global.hpp @@ -190,7 +190,7 @@ name + "." + std::to_string(vm().getTrajectory()) + "." + resultFileExt // default Schur convention #ifndef HADRONS_DEFAULT_SCHUR -#define HADRONS_DEFAULT_SCHUR DiagMooee +#define HADRONS_DEFAULT_SCHUR DiagTwo #endif #define _HADRONS_SCHUR_OP_(conv) Schur##conv##Operator #define HADRONS_SCHUR_OP(conv) _HADRONS_SCHUR_OP_(conv) From d516938707bf4fea9da60f7333126de8be77f992 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Wed, 14 Mar 2018 14:54:25 +0000 Subject: [PATCH 22/29] Hadrons: eigen packs I/O and deflation interface --- extras/Hadrons/EigenPack.hpp | 19 ++- extras/Hadrons/Environment.hpp | 22 +-- extras/Hadrons/Modules.hpp | 2 + extras/Hadrons/Modules/MGauge/FundtoHirep.cc | 12 +- .../Modules/MIO/LoadCoarseEigenPack.hpp | 126 ++++++++++++++++++ extras/Hadrons/Modules/MIO/LoadEigenPack.hpp | 121 +++++++++++++++++ extras/Hadrons/Modules/MSolver/RBPrecCG.hpp | 80 ++++++++--- extras/Hadrons/VirtualMachine.cc | 8 +- extras/Hadrons/modules.inc | 2 + lib/algorithms/iterative/Deflation.h | 27 ++-- lib/algorithms/iterative/SchurRedBlack.h | 8 +- lib/serialisation/BaseIO.h | 2 +- 12 files changed, 368 insertions(+), 61 deletions(-) create mode 100644 extras/Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp create mode 100644 extras/Hadrons/Modules/MIO/LoadEigenPack.hpp diff --git a/extras/Hadrons/EigenPack.hpp b/extras/Hadrons/EigenPack.hpp index eeb04d17..d27e35b9 100644 --- a/extras/Hadrons/EigenPack.hpp +++ b/extras/Hadrons/EigenPack.hpp @@ -29,6 +29,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #define Hadrons_EigenPack_hpp_ #include +#include #include BEGIN_HADRONS_NAMESPACE @@ -38,12 +39,14 @@ BEGIN_HADRONS_NAMESPACE #define HADRONS_DEFAULT_LANCZOS_NBASIS 60 #endif -template +template class EigenPack { +public: + typedef F Field; public: std::vector eval; - std::vector evec; + std::vector evec; public: EigenPack(void) = default; virtual ~EigenPack(void) = default; @@ -123,12 +126,14 @@ protected: } }; -template -class CoarseEigenPack: public EigenPack +template +class CoarseEigenPack: public EigenPack { public: - std::vector evalCoarse; - std::vector evecCoarse; + typedef CoarseF CoarseField; +public: + std::vector evalCoarse; + std::vector evecCoarse; public: CoarseEigenPack(void) = default; virtual ~CoarseEigenPack(void) = default; @@ -142,7 +147,7 @@ public: void resize(const size_t sizeFine, const size_t sizeCoarse, GridBase *gridFine, GridBase *gridCoarse) { - EigenPack::resize(sizeFine, gridFine); + EigenPack::resize(sizeFine, gridFine); evalCoarse.resize(sizeCoarse); evecCoarse.resize(sizeCoarse, gridCoarse); } diff --git a/extras/Hadrons/Environment.hpp b/extras/Hadrons/Environment.hpp index 637962b1..0fb81250 100644 --- a/extras/Hadrons/Environment.hpp +++ b/extras/Hadrons/Environment.hpp @@ -78,7 +78,7 @@ private: Size size{0}; Storage storage{Storage::object}; unsigned int Ls{0}; - const std::type_info *type{nullptr}; + const std::type_info *type{nullptr}, *derivedType{nullptr}; std::string name; int module{-1}; std::unique_ptr data{nullptr}; @@ -230,22 +230,24 @@ void Environment::createDerivedObject(const std::string name, { MemoryProfiler::stats = &memStats; } - size_t initMem = MemoryProfiler::stats->currentlyAllocated; - object_[address].storage = storage; - object_[address].Ls = Ls; + size_t initMem = MemoryProfiler::stats->currentlyAllocated; + object_[address].storage = storage; + object_[address].Ls = Ls; object_[address].data.reset(new Holder(new T(std::forward(args)...))); - object_[address].size = MemoryProfiler::stats->maxAllocated - initMem; - object_[address].type = &typeid(B); + object_[address].size = MemoryProfiler::stats->maxAllocated - initMem; + object_[address].type = &typeid(B); + object_[address].derivedType = &typeid(T); if (MemoryProfiler::stats == &memStats) { MemoryProfiler::stats = nullptr; } } // object already exists, no error if it is a cache, error otherwise - else if ((object_[address].storage != Storage::cache) or - (object_[address].storage != storage) or - (object_[address].name != name) or - (object_[address].type != &typeid(B))) + else if ((object_[address].storage != Storage::cache) or + (object_[address].storage != storage) or + (object_[address].name != name) or + (object_[address].type != &typeid(B)) or + (object_[address].derivedType != &typeid(T))) { HADRON_ERROR(Definition, "object '" + name + "' already allocated"); } diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index ee53faa8..528faecd 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -70,5 +70,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include #include +#include #include diff --git a/extras/Hadrons/Modules/MGauge/FundtoHirep.cc b/extras/Hadrons/Modules/MGauge/FundtoHirep.cc index 31c5a34d..ba8e9330 100644 --- a/extras/Hadrons/Modules/MGauge/FundtoHirep.cc +++ b/extras/Hadrons/Modules/MGauge/FundtoHirep.cc @@ -42,7 +42,8 @@ TFundtoHirep::TFundtoHirep(const std::string name) template std::vector TFundtoHirep::getInput(void) { - std::vector in; + std::vector in = {par().gaugeconf}; + return in; } @@ -50,6 +51,7 @@ template std::vector TFundtoHirep::getOutput(void) { std::vector out = {getName()}; + return out; } @@ -57,19 +59,19 @@ std::vector TFundtoHirep::getOutput(void) template void TFundtoHirep::setup(void) { - envCreateLat(typename Rep::LatticeField, getName()); + envCreateLat(Rep::LatticeField, getName()); } // execution /////////////////////////////////////////////////////////////////// template void TFundtoHirep::execute(void) { - auto &U = *env().template getObject(par().gaugeconf); LOG(Message) << "Transforming Representation" << std::endl; + auto &U = envGet(LatticeGaugeField, par().gaugeconf); + auto &URep = envGet(Rep::LatticeField, getName()); + Rep TargetRepresentation(U._grid); TargetRepresentation.update_representation(U); - - auto &URep = envGet(typename Rep::LatticeField, getName()); URep = TargetRepresentation.U; } diff --git a/extras/Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp b/extras/Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp new file mode 100644 index 00000000..77c3a7ee --- /dev/null +++ b/extras/Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp @@ -0,0 +1,126 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MIO/LoadCoarseEigenPack.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_LoadCoarseEigenPack_hpp_ +#define Hadrons_MIO_LoadCoarseEigenPack_hpp_ + +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Load local coherence eigen vectors/values package * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MIO) + +class LoadCoarseEigenPackPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(LoadCoarseEigenPackPar, + std::string, filestem, + unsigned int, sizeFine, + unsigned int, sizeCoarse, + unsigned int, Ls, + std::vector, blockSize); +}; + +template +class TLoadCoarseEigenPack: public Module +{ +public: + typedef CoarseEigenPack BasePack; +public: + // constructor + TLoadCoarseEigenPack(const std::string name); + // destructor + virtual ~TLoadCoarseEigenPack(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(LoadCoarseFermionEigenPack, + ARG(TLoadCoarseEigenPack>), MIO); + +/****************************************************************************** + * TLoadCoarseEigenPack implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TLoadCoarseEigenPack::TLoadCoarseEigenPack(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TLoadCoarseEigenPack::getInput(void) +{ + std::vector in; + + return in; +} + +template +std::vector TLoadCoarseEigenPack::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TLoadCoarseEigenPack::setup(void) +{ + env().createGrid(par().Ls); + env().createCoarseGrid(par().blockSize, par().Ls); + envCreateDerived(BasePack, Pack, getName(), par().Ls, par().sizeFine, + par().sizeCoarse, env().getRbGrid(par().Ls), + env().getCoarseGrid(par().blockSize, par().Ls)); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TLoadCoarseEigenPack::execute(void) +{ + auto &epack = envGetDerived(BasePack, Pack, getName()); + + epack.read(par().filestem, vm().getTrajectory()); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MIO_LoadCoarseEigenPack_hpp_ diff --git a/extras/Hadrons/Modules/MIO/LoadEigenPack.hpp b/extras/Hadrons/Modules/MIO/LoadEigenPack.hpp new file mode 100644 index 00000000..bcc3f22b --- /dev/null +++ b/extras/Hadrons/Modules/MIO/LoadEigenPack.hpp @@ -0,0 +1,121 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MIO/LoadEigenPack.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_LoadEigenPack_hpp_ +#define Hadrons_MIO_LoadEigenPack_hpp_ + +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Load eigen vectors/values package * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MIO) + +class LoadEigenPackPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(LoadEigenPackPar, + std::string, filestem, + unsigned int, size, + unsigned int, Ls); +}; + +template +class TLoadEigenPack: public Module +{ +public: + typedef EigenPack BasePack; +public: + // constructor + TLoadEigenPack(const std::string name); + // destructor + virtual ~TLoadEigenPack(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(LoadFermionEigenPack, TLoadEigenPack>, MIO); + +/****************************************************************************** + * TLoadEigenPack implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TLoadEigenPack::TLoadEigenPack(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TLoadEigenPack::getInput(void) +{ + std::vector in; + + return in; +} + +template +std::vector TLoadEigenPack::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TLoadEigenPack::setup(void) +{ + env().createGrid(par().Ls); + envCreateDerived(BasePack, Pack, getName(), par().Ls, par().size, + env().getRbGrid(par().Ls)); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TLoadEigenPack::execute(void) +{ + auto &epack = envGetDerived(BasePack, Pack, getName()); + + epack.read(par().filestem, vm().getTrajectory()); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MIO_LoadEigenPack_hpp_ diff --git a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp index 4bff910e..4af81537 100644 --- a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp +++ b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp @@ -32,6 +32,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include BEGIN_HADRONS_NAMESPACE @@ -44,15 +45,22 @@ class RBPrecCGPar: Serializable { public: GRID_SERIALIZABLE_CLASS_MEMBERS(RBPrecCGPar , - std::string , action, - unsigned int , maxIteration, - double , residual); + std::string , action, + unsigned int, maxIteration, + double , residual, + std::string , eigenPack); }; -template +template class TRBPrecCG: public Module { public: + typedef FermionEigenPack EPack; + typedef CoarseFermionEigenPack CoarseEPack; + typedef DeflatedGuesser FineGuesser; + typedef LocalCoherenceDeflatedGuesser< + typename FImpl::FermionField, + typename CoarseEPack::CoarseField> CoarseGuesser; FGS_TYPE_ALIASES(FImpl,); public: // constructor @@ -70,37 +78,39 @@ protected: virtual void execute(void); }; -MODULE_REGISTER_NS(RBPrecCG, TRBPrecCG, MSolver); -MODULE_REGISTER_NS(ZRBPrecCG, TRBPrecCG, MSolver); +MODULE_REGISTER_NS(RBPrecCG, + ARG(TRBPrecCG), MSolver); +MODULE_REGISTER_NS(ZRBPrecCG, + ARG(TRBPrecCG), MSolver); /****************************************************************************** * TRBPrecCG template implementation * ******************************************************************************/ // constructor ///////////////////////////////////////////////////////////////// -template -TRBPrecCG::TRBPrecCG(const std::string name) +template +TRBPrecCG::TRBPrecCG(const std::string name) : Module(name) {} // dependencies/products /////////////////////////////////////////////////////// -template -std::vector TRBPrecCG::getInput(void) +template +std::vector TRBPrecCG::getInput(void) { std::vector in = {}; return in; } -template -std::vector TRBPrecCG::getReference(void) +template +std::vector TRBPrecCG::getReference(void) { std::vector ref = {par().action}; return ref; } -template -std::vector TRBPrecCG::getOutput(void) +template +std::vector TRBPrecCG::getOutput(void) { std::vector out = {getName()}; @@ -108,8 +118,8 @@ std::vector TRBPrecCG::getOutput(void) } // setup /////////////////////////////////////////////////////////////////////// -template -void TRBPrecCG::setup(void) +template +void TRBPrecCG::setup(void) { if (par().maxIteration == 0) { @@ -121,23 +131,49 @@ void TRBPrecCG::setup(void) << par().residual << ", maximum iteration " << par().maxIteration << std::endl; - auto Ls = env().getObjectLs(par().action); - auto &mat = envGet(FMat, par().action); - auto solver = [&mat, this](FermionField &sol, const FermionField &source) + auto Ls = env().getObjectLs(par().action); + auto &mat = envGet(FMat, par().action); + std::string guesserName = getName() + "_guesser"; + + if (par().eigenPack.empty()) + { + env().template createDerivedObject, ZeroGuesser> + (guesserName, Environment::Storage::object, Ls); + } + else + { + try + { + auto &epack = envGetDerived(EPack, CoarseEPack, par().eigenPack); + + envCreateDerived(Guesser, CoarseGuesser, + guesserName, Ls, epack.evec, epack.evecCoarse, + epack.evalCoarse); + } + catch (Exceptions::Definition &e) + { + auto &epack = envGet(EPack, par().eigenPack); + + envCreateDerived(Guesser, FineGuesser, + guesserName, Ls, epack.evec, epack.eval); + } + } + auto &guesser = envGet(Guesser, guesserName); + auto solver = [&mat, &guesser, this](FermionField &sol, const FermionField &source) { ConjugateGradient cg(par().residual, par().maxIteration); HADRONS_DEFAULT_SCHUR_SOLVE schurSolver(cg); - schurSolver(mat, source, sol); + schurSolver(mat, source, sol, guesser); }; envCreate(SolverFn, getName(), Ls, solver); } // execution /////////////////////////////////////////////////////////////////// -template -void TRBPrecCG::execute(void) +template +void TRBPrecCG::execute(void) {} END_MODULE_NAMESPACE diff --git a/extras/Hadrons/VirtualMachine.cc b/extras/Hadrons/VirtualMachine.cc index ba963276..2b7f9620 100644 --- a/extras/Hadrons/VirtualMachine.cc +++ b/extras/Hadrons/VirtualMachine.cc @@ -632,12 +632,14 @@ void VirtualMachine::executeProgram(const Program &p) const freeProg = makeGarbageSchedule(p); for (unsigned int i = 0; i < freeProg.size(); ++i) { - LOG(Debug) << std::setw(4) << i + 1 << ": ["; + std::string msg = ""; + for (auto &a: freeProg[i]) { - std::cout << env().getObjectName(a) << " "; + msg += env().getObjectName(a) + " "; } - std::cout << "]" << std::endl; + msg += "]"; + LOG(Debug) << std::setw(4) << i + 1 << ": [" << msg << std::endl; } // program execution diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index c0be8aef..5932cc89 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -53,6 +53,8 @@ modules_hpp =\ Modules/MScalarSUN/Utils.hpp \ Modules/MScalarSUN/TransProj.hpp \ Modules/MScalarSUN/TrKinetic.hpp \ + Modules/MIO/LoadEigenPack.hpp \ Modules/MIO/LoadNersc.hpp \ + Modules/MIO/LoadCoarseEigenPack.hpp \ Modules/MIO/LoadBinary.hpp diff --git a/lib/algorithms/iterative/Deflation.h b/lib/algorithms/iterative/Deflation.h index b2239c55..316afe90 100644 --- a/lib/algorithms/iterative/Deflation.h +++ b/lib/algorithms/iterative/Deflation.h @@ -30,22 +30,31 @@ Author: Peter Boyle namespace Grid { -struct ZeroGuesser { +template +class Guesser { public: - template - void operator()(const Field &src,Field &guess) { guess = Zero(); }; + Guesser(void) = default; + virtual ~Guesser(void) = default; + virtual void operator()(const Field &src, Field &guess) = 0; }; -struct SourceGuesser { + +template +class ZeroGuesser: public Guesser { public: - template - void operator()(const Field &src,Field &guess) { guess = src; }; + virtual void operator()(const Field &src, Field &guess) { guess = zero; }; +}; + +template +class SourceGuesser: public Guesser { +public: + virtual void operator()(const Field &src, Field &guess) { guess = src; }; }; //////////////////////////////// // Fine grid deflation //////////////////////////////// template -struct DeflatedGuesser { +class DeflatedGuesser: public Guesser { private: const std::vector &evec; const std::vector &eval; @@ -54,7 +63,7 @@ public: DeflatedGuesser(const std::vector & _evec,const std::vector & _eval) : evec(_evec), eval(_eval) {}; - void operator()(const Field &src,Field &guess) { + virtual void operator()(const Field &src,Field &guess) { guess = zero; assert(evec.size()==eval.size()); auto N = evec.size(); @@ -66,7 +75,7 @@ public: }; template -class LocalCoherenceDeflatedGuesser { +class LocalCoherenceDeflatedGuesser: public Guesser { private: const std::vector &subspace; const std::vector &evec_coarse; diff --git a/lib/algorithms/iterative/SchurRedBlack.h b/lib/algorithms/iterative/SchurRedBlack.h index fac2030f..091330b2 100644 --- a/lib/algorithms/iterative/SchurRedBlack.h +++ b/lib/algorithms/iterative/SchurRedBlack.h @@ -108,7 +108,7 @@ namespace Grid { template void operator() (Matrix & _Matrix,const Field &in, Field &out){ - ZeroGuesser guess; + ZeroGuesser guess; (*this)(_Matrix,in,out,guess); } template @@ -195,7 +195,7 @@ namespace Grid { }; template void operator() (Matrix & _Matrix,const Field &in, Field &out){ - ZeroGuesser guess; + ZeroGuesser guess; (*this)(_Matrix,in,out,guess); } template @@ -280,7 +280,7 @@ namespace Grid { template void operator() (Matrix & _Matrix,const Field &in, Field &out){ - ZeroGuesser guess; + ZeroGuesser guess; (*this)(_Matrix,in,out,guess); } template @@ -365,7 +365,7 @@ namespace Grid { template void operator() (Matrix & _Matrix,const Field &in, Field &out){ - ZeroGuesser guess; + ZeroGuesser guess; (*this)(_Matrix,in,out,guess); } template diff --git a/lib/serialisation/BaseIO.h b/lib/serialisation/BaseIO.h index 298bff87..bc178e0d 100644 --- a/lib/serialisation/BaseIO.h +++ b/lib/serialisation/BaseIO.h @@ -232,7 +232,7 @@ namespace Grid { { is >> std::boolalpha >> output; } - catch(std::istringstream::failure &e) + catch(std::ios_base::failure &e) { std::cerr << "numerical conversion failure on '" << s << "' "; std::cerr << "(typeid: " << typeid(U).name() << ")" << std::endl; From d5ce66f6ab2c44a12def7b6d26df80d6e646b1fb Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Fri, 16 Mar 2018 21:37:03 +0000 Subject: [PATCH 23/29] Extra SHM option --- configure.ac | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/configure.ac b/configure.ac index 3a6a2960..aced6a9c 100644 --- a/configure.ac +++ b/configure.ac @@ -340,7 +340,7 @@ case ${ac_PRECISION} in esac ###################### Shared memory allocation technique under MPI3 -AC_ARG_ENABLE([shm],[AC_HELP_STRING([--enable-shm=shmopen|hugetlbfs], +AC_ARG_ENABLE([shm],[AC_HELP_STRING([--enable-shm=shmopen|hugetlbfs|shmnone], [Select SHM allocation technique])],[ac_SHM=${enable_shm}],[ac_SHM=shmopen]) case ${ac_SHM} in @@ -349,6 +349,10 @@ case ${ac_SHM} in AC_DEFINE([GRID_MPI3_SHMOPEN],[1],[GRID_MPI3_SHMOPEN] ) ;; + shmnone) + AC_DEFINE([GRID_MPI3_SHM_NONE],[1],[GRID_MPI3_SHM_NONE] ) + ;; + hugetlbfs) AC_DEFINE([GRID_MPI3_SHMMMAP],[1],[GRID_MPI3_SHMMMAP] ) ;; From 01568b0e62d94ac5d79da8c3edd1db9eea5928e3 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Fri, 16 Mar 2018 21:54:28 +0000 Subject: [PATCH 24/29] Add a new SHM option --- lib/communicator/SharedMemoryMPI.cc | 44 ++++++++++++++++++++++++++++- 1 file changed, 43 insertions(+), 1 deletion(-) diff --git a/lib/communicator/SharedMemoryMPI.cc b/lib/communicator/SharedMemoryMPI.cc index 45edbb07..1fa84dfb 100644 --- a/lib/communicator/SharedMemoryMPI.cc +++ b/lib/communicator/SharedMemoryMPI.cc @@ -226,6 +226,48 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags) }; #endif // MMAP +#ifdef GRID_MPI3_SHM_NONE +void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags) +{ + std::cout << "SharedMemoryAllocate "<< bytes<< " MMAP anonymous implementation "< Date: Fri, 16 Mar 2018 21:54:56 +0000 Subject: [PATCH 25/29] 4GB clean the offsets in parallel IO for multifile records --- lib/parallelIO/BinaryIO.h | 34 +++++++++++++++------------- lib/parallelIO/IldgIO.h | 47 +++++++++++++++++++++++++-------------- lib/parallelIO/NerscIO.h | 10 ++++----- 3 files changed, 53 insertions(+), 38 deletions(-) diff --git a/lib/parallelIO/BinaryIO.h b/lib/parallelIO/BinaryIO.h index b40a75af..39acf0e0 100644 --- a/lib/parallelIO/BinaryIO.h +++ b/lib/parallelIO/BinaryIO.h @@ -91,7 +91,7 @@ class BinaryIO { typedef typename vobj::scalar_object sobj; GridBase *grid = lat._grid; - int lsites = grid->lSites(); + uint64_t lsites = grid->lSites(); std::vector scalardata(lsites); unvectorizeToLexOrdArray(scalardata,lat); @@ -160,7 +160,9 @@ class BinaryIO { /* * Scidac csum is rather more heavyweight + * FIXME -- 128^3 x 256 x 16 will overflow. */ + int global_site; Lexicographic::CoorFromIndex(coor,local_site,local_vol); @@ -261,7 +263,7 @@ class BinaryIO { GridBase *grid, std::vector &iodata, std::string file, - Integer offset, + uint64_t offset, const std::string &format, int control, uint32_t &nersc_csum, uint32_t &scidac_csuma, @@ -523,7 +525,7 @@ class BinaryIO { static inline void readLatticeObject(Lattice &Umu, std::string file, munger munge, - Integer offset, + uint64_t offset, const std::string &format, uint32_t &nersc_csum, uint32_t &scidac_csuma, @@ -533,7 +535,7 @@ class BinaryIO { typedef typename vobj::Realified::scalar_type word; word w=0; GridBase *grid = Umu._grid; - int lsites = grid->lSites(); + uint64_t lsites = grid->lSites(); std::vector scalardata(lsites); std::vector iodata(lsites); // Munge, checksum, byte order in here @@ -544,7 +546,7 @@ class BinaryIO { GridStopWatch timer; timer.Start(); - parallel_for(int x=0;xBarrier(); @@ -560,7 +562,7 @@ class BinaryIO { static inline void writeLatticeObject(Lattice &Umu, std::string file, munger munge, - Integer offset, + uint64_t offset, const std::string &format, uint32_t &nersc_csum, uint32_t &scidac_csuma, @@ -569,7 +571,7 @@ class BinaryIO { typedef typename vobj::scalar_object sobj; typedef typename vobj::Realified::scalar_type word; word w=0; GridBase *grid = Umu._grid; - int lsites = grid->lSites(); + uint64_t lsites = grid->lSites(); std::vector scalardata(lsites); std::vector iodata(lsites); // Munge, checksum, byte order in here @@ -580,7 +582,7 @@ class BinaryIO { GridStopWatch timer; timer.Start(); unvectorizeToLexOrdArray(scalardata,Umu); - parallel_for(int x=0;xBarrier(); timer.Stop(); @@ -597,7 +599,7 @@ class BinaryIO { static inline void readRNG(GridSerialRNG &serial, GridParallelRNG ¶llel, std::string file, - Integer offset, + uint64_t offset, uint32_t &nersc_csum, uint32_t &scidac_csuma, uint32_t &scidac_csumb) @@ -610,8 +612,8 @@ class BinaryIO { std::string format = "IEEE32BIG"; GridBase *grid = parallel._grid; - int gsites = grid->gSites(); - int lsites = grid->lSites(); + uint64_t gsites = grid->gSites(); + uint64_t lsites = grid->lSites(); uint32_t nersc_csum_tmp = 0; uint32_t scidac_csuma_tmp = 0; @@ -626,7 +628,7 @@ class BinaryIO { nersc_csum,scidac_csuma,scidac_csumb); timer.Start(); - parallel_for(int lidx=0;lidx tmp(RngStateCount); std::copy(iodata[lidx].begin(),iodata[lidx].end(),tmp.begin()); parallel.SetState(tmp,lidx); @@ -659,7 +661,7 @@ class BinaryIO { static inline void writeRNG(GridSerialRNG &serial, GridParallelRNG ¶llel, std::string file, - Integer offset, + uint64_t offset, uint32_t &nersc_csum, uint32_t &scidac_csuma, uint32_t &scidac_csumb) @@ -670,8 +672,8 @@ class BinaryIO { typedef std::array RNGstate; GridBase *grid = parallel._grid; - int gsites = grid->gSites(); - int lsites = grid->lSites(); + uint64_t gsites = grid->gSites(); + uint64_t lsites = grid->lSites(); uint32_t nersc_csum_tmp; uint32_t scidac_csuma_tmp; @@ -684,7 +686,7 @@ class BinaryIO { timer.Start(); std::vector iodata(lsites); - parallel_for(int lidx=0;lidx tmp(RngStateCount); parallel.GetState(tmp,lidx); std::copy(tmp.begin(),tmp.end(),iodata[lidx].begin()); diff --git a/lib/parallelIO/IldgIO.h b/lib/parallelIO/IldgIO.h index b0bd7e2c..8655b24c 100644 --- a/lib/parallelIO/IldgIO.h +++ b/lib/parallelIO/IldgIO.h @@ -337,6 +337,20 @@ class GridLimeWriter : public BinaryIO { template void writeLimeLatticeBinaryObject(Lattice &field,std::string record_name) { + //////////////////////////////////////////////////////////////////// + // NB: FILE and iostream are jointly writing disjoint sequences in the + // the same file through different file handles (integer units). + // + // These are both buffered, so why I think this code is right is as follows. + // + // i) write record header to FILE *File, telegraphing the size; flush + // ii) ftello reads the offset from FILE *File . + // iii) iostream / MPI Open independently seek this offset. Write sequence direct to disk. + // Closes iostream and flushes. + // iv) fseek on FILE * to end of this disjoint section. + // v) Continue writing scidac record. + //////////////////////////////////////////////////////////////////// + //////////////////////////////////////////// // Create record header //////////////////////////////////////////// @@ -350,25 +364,24 @@ class GridLimeWriter : public BinaryIO { // std::cout << "W Gsites " <_gsites<(); BinarySimpleMunger munge; - BinaryIO::writeLatticeObject(field, filename, munge, offset, format,nersc_csum,scidac_csuma,scidac_csumb); - // fseek(File,0,SEEK_END); offset = ftello(File);std::cout << " offset now "<(field, filename, munge, offset1, format,nersc_csum,scidac_csuma,scidac_csumb); + + /////////////////////////////////////////// + // Wind forward and close the record + /////////////////////////////////////////// + fseek(File,0,SEEK_END); + unt64_t offset2 = ftello(File); // std::cout << " now at offset "<=0); //////////////////////////////////////// diff --git a/lib/parallelIO/NerscIO.h b/lib/parallelIO/NerscIO.h index 786839f2..e2c2bc39 100644 --- a/lib/parallelIO/NerscIO.h +++ b/lib/parallelIO/NerscIO.h @@ -57,7 +57,7 @@ namespace Grid { // for the header-reader static inline int readHeader(std::string file,GridBase *grid, FieldMetaData &field) { - int offset=0; + uint64_t offset=0; std::map header; std::string line; @@ -139,7 +139,7 @@ namespace Grid { typedef Lattice > GaugeField; GridBase *grid = Umu._grid; - int offset = readHeader(file,Umu._grid,header); + uint64_t offset = readHeader(file,Umu._grid,header); FieldMetaData clone(header); @@ -236,7 +236,7 @@ namespace Grid { GaugeStatistics(Umu,header); MachineCharacteristics(header); - int offset; + uint64_t offset; truncate(file); @@ -278,7 +278,7 @@ namespace Grid { header.plaquette=0.0; MachineCharacteristics(header); - int offset; + uint64_t offset; #ifdef RNG_RANLUX header.floating_point = std::string("UINT64"); @@ -313,7 +313,7 @@ namespace Grid { GridBase *grid = parallel._grid; - int offset = readHeader(file,grid,header); + uint64_t offset = readHeader(file,grid,header); FieldMetaData clone(header); From e1dcfd35538cf0e7ebe39d0c5f55c7427b921d38 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Fri, 16 Mar 2018 23:10:47 +0000 Subject: [PATCH 26/29] typo fix --- lib/parallelIO/IldgIO.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/parallelIO/IldgIO.h b/lib/parallelIO/IldgIO.h index 8655b24c..b81d1e43 100644 --- a/lib/parallelIO/IldgIO.h +++ b/lib/parallelIO/IldgIO.h @@ -378,7 +378,7 @@ class GridLimeWriter : public BinaryIO { // Wind forward and close the record /////////////////////////////////////////// fseek(File,0,SEEK_END); - unt64_t offset2 = ftello(File); // std::cout << " now at offset "< Date: Sat, 17 Mar 2018 09:35:01 +0000 Subject: [PATCH 27/29] Drop RB on coarse space ; that was a mistake --- tests/lanczos/Test_dwf_compressed_lanczos_reorg.cc | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/lanczos/Test_dwf_compressed_lanczos_reorg.cc b/tests/lanczos/Test_dwf_compressed_lanczos_reorg.cc index 3dff4b90..b55b66d9 100644 --- a/tests/lanczos/Test_dwf_compressed_lanczos_reorg.cc +++ b/tests/lanczos/Test_dwf_compressed_lanczos_reorg.cc @@ -180,7 +180,6 @@ int main (int argc, char ** argv) { GridCartesian * CoarseGrid4 = SpaceTimeGrid::makeFourDimGrid(coarseLatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); GridRedBlackCartesian * CoarseGrid4rb = SpaceTimeGrid::makeFourDimRedBlackGrid(CoarseGrid4); GridCartesian * CoarseGrid5 = SpaceTimeGrid::makeFiveDimGrid(cLs,CoarseGrid4); - GridRedBlackCartesian * CoarseGrid5rb = SpaceTimeGrid::makeFourDimRedBlackGrid(CoarseGrid5); // Gauge field LatticeGaugeField Umu(UGrid); @@ -206,7 +205,7 @@ int main (int argc, char ** argv) { const int nbasis= 60; assert(nbasis==Ns1); - LocalCoherenceLanczosScidac _LocalCoherenceLanczos(FrbGrid,CoarseGrid5rb,HermOp,Odd); + LocalCoherenceLanczosScidac _LocalCoherenceLanczos(FrbGrid,CoarseGrid5,HermOp,Odd); std::cout << GridLogMessage << "Constructed LocalCoherenceLanczos" << std::endl; assert( (Params.doFine)||(Params.doFineRead)); From cbc73a3fd176dd3b63e7bfafaf7a8dc35e7f27af Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Mon, 19 Mar 2018 13:11:38 +0000 Subject: [PATCH 28/29] Hadrons: CG guesser fix --- extras/Hadrons/Modules/MSolver/RBPrecCG.hpp | 27 ++++++++++----------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp index 4af81537..8d1b0fc6 100644 --- a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp +++ b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp @@ -55,13 +55,14 @@ template class TRBPrecCG: public Module { public: + FGS_TYPE_ALIASES(FImpl,); typedef FermionEigenPack EPack; typedef CoarseFermionEigenPack CoarseEPack; + typedef std::shared_ptr> GuesserPt; typedef DeflatedGuesser FineGuesser; typedef LocalCoherenceDeflatedGuesser< typename FImpl::FermionField, typename CoarseEPack::CoarseField> CoarseGuesser; - FGS_TYPE_ALIASES(FImpl,); public: // constructor TRBPrecCG(const std::string name); @@ -131,41 +132,39 @@ void TRBPrecCG::setup(void) << par().residual << ", maximum iteration " << par().maxIteration << std::endl; - auto Ls = env().getObjectLs(par().action); - auto &mat = envGet(FMat, par().action); + auto Ls = env().getObjectLs(par().action); + auto &mat = envGet(FMat, par().action); std::string guesserName = getName() + "_guesser"; + GuesserPt guesser{nullptr}; if (par().eigenPack.empty()) { - env().template createDerivedObject, ZeroGuesser> - (guesserName, Environment::Storage::object, Ls); + guesser.reset(new ZeroGuesser()); } else { try { auto &epack = envGetDerived(EPack, CoarseEPack, par().eigenPack); - - envCreateDerived(Guesser, CoarseGuesser, - guesserName, Ls, epack.evec, epack.evecCoarse, - epack.evalCoarse); + + guesser.reset(new CoarseGuesser(epack.evec, epack.evecCoarse, + epack.evalCoarse)); } catch (Exceptions::Definition &e) { auto &epack = envGet(EPack, par().eigenPack); - envCreateDerived(Guesser, FineGuesser, - guesserName, Ls, epack.evec, epack.eval); + guesser.reset(new FineGuesser(epack.evec, epack.eval)); } } - auto &guesser = envGet(Guesser, guesserName); - auto solver = [&mat, &guesser, this](FermionField &sol, const FermionField &source) + auto solver = [&mat, guesser, this](FermionField &sol, + const FermionField &source) { ConjugateGradient cg(par().residual, par().maxIteration); HADRONS_DEFAULT_SCHUR_SOLVE schurSolver(cg); - schurSolver(mat, source, sol, guesser); + schurSolver(mat, source, sol, *guesser); }; envCreate(SolverFn, getName(), Ls, solver); } From 62702dbcb805543ed1beb66fed073f5cc2a79e85 Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Mon, 19 Mar 2018 17:56:53 +0000 Subject: [PATCH 29/29] Fixing bug in the Point sink causing NaNs --- extras/Hadrons/Modules/MSink/Point.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/extras/Hadrons/Modules/MSink/Point.hpp b/extras/Hadrons/Modules/MSink/Point.hpp index c5f6eff0..ee824c03 100644 --- a/extras/Hadrons/Modules/MSink/Point.hpp +++ b/extras/Hadrons/Modules/MSink/Point.hpp @@ -128,7 +128,7 @@ void TPoint::execute(void) envGetTmp(LatticeComplex, coor); p = strToVec(par().mom); ph = zero; - for(unsigned int mu = 0; mu < env().getNd(); mu++) + for(unsigned int mu = 0; mu < p.size(); mu++) { LatticeCoordinate(coor, mu); ph = ph + (p[mu]/env().getGrid()->_fdimensions[mu])*coor;