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/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] ) ;; diff --git a/extras/Hadrons/EigenPack.hpp b/extras/Hadrons/EigenPack.hpp new file mode 100644 index 00000000..d27e35b9 --- /dev/null +++ b/extras/Hadrons/EigenPack.hpp @@ -0,0 +1,218 @@ +/************************************************************************************* + +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 +#include + +BEGIN_HADRONS_NAMESPACE + +// Lanczos type +#ifndef HADRONS_DEFAULT_LANCZOS_NBASIS +#define HADRONS_DEFAULT_LANCZOS_NBASIS 60 +#endif + +template +class EigenPack +{ +public: + typedef F Field; +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: + typedef CoarseF CoarseField; +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..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}; @@ -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; @@ -110,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 @@ -155,6 +165,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_; @@ -176,7 +190,7 @@ Holder::Holder(T *pt) template T & Holder::get(void) const { - return &objPt_.get(); + return *objPt_.get(); } template @@ -216,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(T); + 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(T))) + 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"); } @@ -246,21 +262,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) + "')"); } } @@ -276,6 +308,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/Global.hpp b/extras/Hadrons/Global.hpp index ed8f4f32..12b9a029 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 @@ -187,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) 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/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) diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 04cc7a62..400c4c8c 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -64,6 +64,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include #include #include #include @@ -72,5 +73,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/MScalarSUN/ShiftProbe.hpp b/extras/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp new file mode 100644 index 00000000..8d52327e --- /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"); + envCreateLat(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/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; diff --git a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp index 6e2103ce..80717f32 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 @@ -78,16 +77,6 @@ public: virtual void setup(void); // 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_; }; MODULE_REGISTER_NS(LocalCoherenceLanczos, @@ -104,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 @@ -121,61 +107,12 @@ std::vector TLocalCoherenceLanczos::getInput(void) template std::vector TLocalCoherenceLanczos::getOutput(void) { - std::vector out = {fineName_, coarseName_}; + std::vector out = {getName()}; return out; } // 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 +120,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 = envGetDerived(BasePack, 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 +147,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 = envGetDerived(BasePack, 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()); } } diff --git a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp index 2b914625..8d1b0fc6 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,16 +45,24 @@ 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: 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; public: // constructor TRBPrecCG(const std::string name); @@ -70,37 +79,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,30 +119,60 @@ std::vector TRBPrecCG::getOutput(void) } // setup /////////////////////////////////////////////////////////////////////// -template -void TRBPrecCG::setup(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); - 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"; + GuesserPt guesser{nullptr}; + + if (par().eigenPack.empty()) + { + guesser.reset(new ZeroGuesser()); + } + else + { + try + { + auto &epack = envGetDerived(EPack, CoarseEPack, par().eigenPack); + + guesser.reset(new CoarseGuesser(epack.evec, epack.evecCoarse, + epack.evalCoarse)); + } + catch (Exceptions::Definition &e) + { + auto &epack = envGet(EPack, par().eigenPack); + + guesser.reset(new FineGuesser(epack.evec, epack.eval)); + } + } + 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 eb1f68ba..2b7f9620 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) @@ -632,6 +630,17 @@ 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) + { + std::string msg = ""; + + for (auto &a: freeProg[i]) + { + msg += env().getObjectName(a) + " "; + } + msg += "]"; + LOG(Debug) << std::setw(4) << i + 1 << ": [" << msg << std::endl; + } // program execution LOG(Debug) << "Executing program..." << std::endl; diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index f7c6c414..55f99f66 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -50,6 +50,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 \ @@ -58,6 +59,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/communicator/SharedMemoryMPI.cc b/lib/communicator/SharedMemoryMPI.cc index 9e5d8f15..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 "<(),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/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 b86e250f..b81d1e43 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); + uint64_t offset2 = ftello(File); // std::cout << " now at offset "<=0); //////////////////////////////////////// @@ -568,7 +581,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); } }; 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); 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/lib/serialisation/BaseIO.h b/lib/serialisation/BaseIO.h index 24e1cec7..bc178e0d 100644 --- a/lib/serialisation/BaseIO.h +++ b/lib/serialisation/BaseIO.h @@ -31,161 +31,10 @@ Author: Guido Cossu #define GRID_SERIALISATION_ABSTRACT_READER_H #include +#include +#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 - { - 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 - 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; @@ -205,6 +54,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; }; @@ -224,6 +79,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); @@ -237,203 +98,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) { @@ -467,6 +134,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 @@ -502,7 +190,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) @@ -514,13 +232,83 @@ 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; 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 94ad9736..12625ab8 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 diff --git a/lib/serialisation/VectorUtils.h b/lib/serialisation/VectorUtils.h new file mode 100644 index 00000000..6df9416d --- /dev/null +++ b/lib/serialisation/VectorUtils.h @@ -0,0 +1,384 @@ +#ifndef GRID_SERIALISATION_VECTORUTILS_H +#define GRID_SERIALISATION_VECTORUTILS_H + +#include +#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 + { + 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/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"<, 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; } }; @@ -93,8 +100,30 @@ 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) { + GridSerialRNG rng; + + rng.SeedFixedIntegers(std::vector({42,10,81,9})); + std::cout << "==== basic IO" << std::endl; XmlWriter WR("bother.xml"); @@ -120,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); @@ -131,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; @@ -141,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; @@ -197,68 +218,11 @@ int main(int argc,char **argv) std::cout << flatdv.getVector() << std::endl; std::cout << 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); - } -*/ - - + std::cout << "==== Grid tensor to vector test" << std::endl; + tensorConvTest(rng, SpinColourMatrix); + tensorConvTest(rng, SpinColourVector); + tensorConvTest(rng, ColourMatrix); + tensorConvTest(rng, ColourVector); + tensorConvTest(rng, SpinMatrix); + tensorConvTest(rng, SpinVector); } 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); 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); 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));