portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-10-25 02:04:48 +01:00
Code Releases Activity

Merge branch 'develop' of https://github.com/paboyle/Grid into develop

Browse Source
This commit is contained in:
fionnoh
2018-09-17 12:19:47 +01:00
parent 7e3647246c 4af6c7e7aa
commit d9de8fd5c9
21 changed files with 639 additions and 376 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
Grid/Makefile.am
View File

@@ -24,7 +24,7 @@ endif
all: version-cache
version-cache:
if [ `git status --porcelain | grep -v '??' | wc -l` -gt 0 ]; then\
@if [ `git status --porcelain | grep -v '??' | wc -l` -gt 0 ]; then\
a="uncommited changes";\
else\
a="clean";\

16
Grid/algorithms/iterative/Deflation.h
View File

@@ -31,21 +31,13 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
namespace Grid {
template<class Field>
class Guesser {
public:
Guesser(void) = default;
virtual ~Guesser(void) = default;
virtual void operator()(const Field &src, Field &guess) = 0;
};
template<class Field>
class ZeroGuesser: public Guesser<Field> {
class ZeroGuesser: public LinearFunction<Field> {
public:
virtual void operator()(const Field &src, Field &guess) { guess = zero; };
};
template<class Field>
class SourceGuesser: public Guesser<Field> {
class SourceGuesser: public LinearFunction<Field> {
public:
virtual void operator()(const Field &src, Field &guess) { guess = src; };
};
@@ -54,7 +46,7 @@ public:
// Fine grid deflation
////////////////////////////////
template<class Field>
class DeflatedGuesser: public Guesser<Field> {
class DeflatedGuesser: public LinearFunction<Field> {
private:
const std::vector<Field> &evec;
const std::vector<RealD> &eval;
@@ -76,7 +68,7 @@ public:
};
template<class FineField, class CoarseField>
class LocalCoherenceDeflatedGuesser: public Guesser<FineField> {
class LocalCoherenceDeflatedGuesser: public LinearFunction<FineField> {
private:
const std::vector<FineField> &subspace;
const std::vector<CoarseField> &evec_coarse;

13
Grid/serialisation/VectorUtils.h
View File

@@ -423,4 +423,17 @@ namespace Grid {
}
}
// helper function to read space-separated values
template <typename T>
std::string vecToStr(const std::vector<T> &v)
{
using Grid::operator<<;
std::ostringstream sstr;
sstr << v;
return sstr.str();
}
#endif

8
Hadrons/Application.cc
View File

@@ -62,10 +62,10 @@ Application::Application(void)
LOG(Message) << std::endl;
LOG(Message) << "** Default parameters (and associated C macros)" << std::endl;
LOG(Message) << "ASCII output precision : " << MACOUT(DEFAULT_ASCII_PREC) << std::endl;
LOG(Message) << "Fermion implementation : " << MACOUTS(FIMPL) << std::endl;
LOG(Message) << "z-Fermion implementation: " << MACOUTS(ZFIMPL) << std::endl;
LOG(Message) << "Scalar implementation : " << MACOUTS(SIMPL) << std::endl;
LOG(Message) << "Gauge implementation : " << MACOUTS(GIMPL) << std::endl;
LOG(Message) << "Fermion implementation : " << MACOUTS(FIMPLBASE) << std::endl;
LOG(Message) << "z-Fermion implementation: " << MACOUTS(ZFIMPLBASE) << std::endl;
LOG(Message) << "Scalar implementation : " << MACOUTS(SIMPLBASE) << std::endl;
LOG(Message) << "Gauge implementation : " << MACOUTS(GIMPLBASE) << std::endl;
LOG(Message) << "Eigenvector base size : "
<< MACOUT(HADRONS_DEFAULT_LANCZOS_NBASIS) << std::endl;
LOG(Message) << "Schur decomposition : " << MACOUTS(HADRONS_DEFAULT_SCHUR) << std::endl;

133
Hadrons/Environment.cc
View File

@@ -43,137 +43,24 @@ HADRONS_ERROR_REF(ObjectDefinition, "no object with address " + std::to_string(a
// constructor /////////////////////////////////////////////////////////////////
Environment::Environment(void)
{
dim_ = GridDefaultLatt();
nd_ = dim_.size();
grid4d_.reset(SpaceTimeGrid::makeFourDimGrid(
dim_, GridDefaultSimd(nd_, vComplex::Nsimd()),
GridDefaultMpi()));
gridRb4d_.reset(SpaceTimeGrid::makeFourDimRedBlackGrid(grid4d_.get()));
dim_ = GridDefaultLatt();
nd_ = dim_.size();
defaultGrid_ = {typeHash<vComplex>(), 1};
grid4d_[defaultGrid_].reset(
SpaceTimeGrid::makeFourDimGrid(dim_,
GridDefaultSimd(nd_, vComplex::Nsimd()),
GridDefaultMpi()));
gridRb4d_[defaultGrid_].reset(
SpaceTimeGrid::makeFourDimRedBlackGrid(grid4d_[defaultGrid_].get()));
vol_ = 1.;
for (auto d: dim_)
{
vol_ *= d;
}
rng4d_.reset(new GridParallelRNG(grid4d_.get()));
rng4d_.reset(new GridParallelRNG(grid4d_[defaultGrid_].get()));
}
// grids ///////////////////////////////////////////////////////////////////////
void Environment::createGrid(const unsigned int Ls)
{
if ((Ls > 1) and (grid5d_.find(Ls) == grid5d_.end()))
{
auto g = getGrid();
grid5d_[Ls].reset(SpaceTimeGrid::makeFiveDimGrid(Ls, g));
gridRb5d_[Ls].reset(SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, g));
}
}
void Environment::createCoarseGrid(const std::vector<int> &blockSize,
const unsigned int Ls)
{
int nd = getNd();
std::vector<int> 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])
{
HADRONS_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)
{
HADRONS_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()));
if (Ls > 1)
{
gridCoarse5d_[key5d].reset(
SpaceTimeGrid::makeFiveDimGrid(cLs, gridCoarse4d_[key4d].get()));
}
}
GridCartesian * Environment::getGrid(const unsigned int Ls) const
{
try
{
if (Ls == 1)
{
return grid4d_.get();
}
else
{
return grid5d_.at(Ls).get();
}
}
catch(std::out_of_range &)
{
HADRONS_ERROR(Definition, "no grid with Ls= " + std::to_string(Ls));
}
}
GridRedBlackCartesian * Environment::getRbGrid(const unsigned int Ls) const
{
try
{
if (Ls == 1)
{
return gridRb4d_.get();
}
else
{
return gridRb5d_.at(Ls).get();
}
}
catch(std::out_of_range &)
{
HADRONS_ERROR(Definition, "no red-black grid with Ls= " + std::to_string(Ls));
}
}
GridCartesian * Environment::getCoarseGrid(
const std::vector<int> &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 &)
{
HADRONS_ERROR(Definition, "no coarse grid with Ls= " + std::to_string(Ls));
}
}
unsigned int Environment::getNd(void) const
{
return nd_;

248
Hadrons/Environment.hpp
View File

@@ -63,6 +63,9 @@ inline Environment & env(void) const\
return Environment::getInstance();\
}
#define DEFINE_ENV_LAMBDA \
auto env = [](void)->Environment &{return Environment::getInstance();}
class Environment
{
SINGLETON(Environment);
@@ -83,15 +86,28 @@ private:
int module{-1};
std::unique_ptr<Object> data{nullptr};
};
typedef std::pair<size_t, unsigned int> FineGridKey;
typedef std::pair<size_t, std::vector<int>> CoarseGridKey;
public:
// grids
template <typename VType = vComplex>
void createGrid(const unsigned int Ls);
template <typename VType = vComplex>
void createCoarseGrid(const std::vector<int> &blockSize,
const unsigned int Ls = 1);
GridCartesian * getGrid(const unsigned int Ls = 1) const;
GridRedBlackCartesian * getRbGrid(const unsigned int Ls = 1) const;
template <typename VType = vComplex>
GridCartesian * getGrid(void) const;
template <typename VType = vComplex>
GridRedBlackCartesian * getRbGrid(void) const;
template <typename VType = vComplex>
GridCartesian * getCoarseGrid(const std::vector<int> &blockSize) const;
template <typename VType = vComplex>
GridCartesian * getGrid(const unsigned int Ls) const;
template <typename VType = vComplex>
GridRedBlackCartesian * getRbGrid(const unsigned int Ls) const;
template <typename VType = vComplex>
GridCartesian * getCoarseGrid(const std::vector<int> &blockSize,
const unsigned int Ls = 1) const;
const unsigned int Ls) const;
std::vector<int> getDim(void) const;
int getDim(const unsigned int mu) const;
unsigned int getNd(void) const;
@@ -154,22 +170,23 @@ public:
void printContent(void) const;
private:
// general
double vol_;
bool protect_{true};
double vol_;
bool protect_{true};
// grids
std::vector<int> dim_;
GridPt grid4d_;
std::map<unsigned int, GridPt> grid5d_;
GridRbPt gridRb4d_;
std::map<unsigned int, GridRbPt> gridRb5d_;
std::map<std::vector<int>, GridPt> gridCoarse4d_;
std::map<std::vector<int>, GridPt> gridCoarse5d_;
unsigned int nd_;
std::vector<int> dim_;
FineGridKey defaultGrid_;
std::map<FineGridKey, GridPt> grid4d_;
std::map<FineGridKey, GridPt> grid5d_;
std::map<FineGridKey, GridRbPt> gridRb4d_;
std::map<FineGridKey, GridRbPt> gridRb5d_;
std::map<CoarseGridKey, GridPt> gridCoarse4d_;
std::map<CoarseGridKey, GridPt> gridCoarse5d_;
unsigned int nd_;
// random number generator
RngPt rng4d_;
RngPt rng4d_;
// object store
std::vector<ObjInfo> object_;
std::map<std::string, unsigned int> objectAddress_;
std::vector<ObjInfo> object_;
std::map<std::string, unsigned int> objectAddress_;
};
/******************************************************************************
@@ -203,6 +220,191 @@ void Holder<T>::reset(T *pt)
/******************************************************************************
* Environment template implementation *
******************************************************************************/
// grids ///////////////////////////////////////////////////////////////////////
template <typename VType>
void Environment::createGrid(const unsigned int Ls)
{
size_t hash = typeHash<VType>();
if (grid4d_.find({hash, 1}) == grid4d_.end())
{
grid4d_[{hash, 1}].reset(
SpaceTimeGrid::makeFourDimGrid(getDim(),
GridDefaultSimd(getNd(), VType::Nsimd()),
GridDefaultMpi()));
gridRb4d_[{hash, 1}].reset(
SpaceTimeGrid::makeFourDimRedBlackGrid(grid4d_[{hash, 1}].get()));
}
if (grid5d_.find({hash, Ls}) == grid5d_.end())
{
auto g = grid4d_[{hash, 1}].get();
grid5d_[{hash, Ls}].reset(SpaceTimeGrid::makeFiveDimGrid(Ls, g));
gridRb5d_[{hash, Ls}].reset(SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, g));
}
}
template <typename VType>
void Environment::createCoarseGrid(const std::vector<int> &blockSize,
const unsigned int Ls)
{
int nd = getNd();
std::vector<int> fineDim = getDim(), coarseDim(nd);
unsigned int cLs;
auto key4d = blockSize, key5d = blockSize;
size_t hash = typeHash<VType>();
createGrid(Ls);
for (int d = 0; d < coarseDim.size(); d++)
{
coarseDim[d] = fineDim[d]/blockSize[d];
if (coarseDim[d]*blockSize[d] != fineDim[d])
{
HADRONS_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)
{
HADRONS_ERROR(Size, "Fine Ls (" + std::to_string(Ls)
+ ") not divisible by coarse Ls ("
+ std::to_string(cLs) + ")");
}
}
else
{
cLs = Ls;
}
key4d.resize(nd);
key5d.push_back(Ls);
CoarseGridKey hkey4d = {hash, key4d}, hkey5d = {hash, key5d};
if (gridCoarse4d_.find(hkey4d) == gridCoarse4d_.end())
{
gridCoarse4d_[hkey4d].reset(
SpaceTimeGrid::makeFourDimGrid(coarseDim,
GridDefaultSimd(nd, VType::Nsimd()), GridDefaultMpi()));
}
if (gridCoarse5d_.find(hkey5d) == gridCoarse5d_.end())
{
gridCoarse5d_[hkey5d].reset(
SpaceTimeGrid::makeFiveDimGrid(cLs, gridCoarse4d_[hkey4d].get()));
}
}
template <typename VType>
GridCartesian * Environment::getGrid(void) const
{
auto it = grid4d_.find({typeHash<VType>(), 1});
if (it != grid4d_.end())
{
return it->second.get();
}
else
{
HADRONS_ERROR(Definition, "no 4D grid for SIMD type '"
+ typeName<VType>() + "'");
}
}
template <typename VType>
GridRedBlackCartesian * Environment::getRbGrid(void) const
{
auto it = gridRb4d_.find({typeHash<VType>(), 1});
if (it != gridRb4d_.end())
{
return it->second.get();
}
else
{
HADRONS_ERROR(Definition, "no 4D red-black grid for SIMD type '"
+ typeName<VType>() + "'");
}
}
template <typename VType>
GridCartesian * Environment::getCoarseGrid(const std::vector<int> &blockSize) const
{
std::vector<int> s = blockSize;
s.resize(getNd());
auto it = gridCoarse4d_.find({typeHash<VType>(), s});
if (it != gridCoarse4d_.end())
{
return it->second.get();
}
else
{
HADRONS_ERROR(Definition, "no 4D coarse grid for SIMD type '"
+ typeName<VType>() + "' and block size "
+ vecToStr(blockSize));
}
}
template <typename VType>
GridCartesian * Environment::getGrid(const unsigned int Ls) const
{
auto it = grid5d_.find({typeHash<VType>(), Ls});
if (it != grid5d_.end())
{
return it->second.get();
}
else
{
HADRONS_ERROR(Definition, "no 5D grid for SIMD type '"
+ typeName<VType>() + "' and Ls = "
+ std::to_string(Ls));
}
}
template <typename VType>
GridRedBlackCartesian * Environment::getRbGrid(const unsigned int Ls) const
{
auto it = gridRb5d_.find({typeHash<VType>(), Ls});
if (it != gridRb5d_.end())
{
return it->second.get();
}
else
{
HADRONS_ERROR(Definition, "no 5D red-black grid for SIMD type '"
+ typeName<VType>() + "' and Ls = "
+ std::to_string(Ls));
}
}
template <typename VType>
GridCartesian * Environment::getCoarseGrid(const std::vector<int> &blockSize,
const unsigned int Ls) const
{
std::vector<int> s = blockSize;
s.push_back(Ls);
auto it = gridCoarse5d_.find({typeHash<VType>(), s});
if (it != gridCoarse5d_.end())
{
return it->second.get();
}
else
{
HADRONS_ERROR(Definition, "no 5D coarse grid for SIMD type '"
+ typeName<VType>() + "', block size "
+ vecToStr(blockSize)
+ " and Ls = " + std::to_string(Ls));
}
}
// general memory management ///////////////////////////////////////////////////
template <typename B, typename T, typename ... Ts>
void Environment::createDerivedObject(const std::string name,
@@ -230,19 +432,19 @@ void Environment::createDerivedObject(const std::string name,
object_[address].Ls = Ls;
object_[address].data.reset(new Holder<B>(new T(std::forward<Ts>(args)...)));
object_[address].size = MemoryProfiler::stats->maxAllocated - initMem;
object_[address].type = &typeid(B);
object_[address].derivedType = &typeid(T);
object_[address].type = typeIdPt<B>();
object_[address].derivedType = typeIdPt<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)) or
(object_[address].derivedType != &typeid(T)))
else if ((object_[address].storage != Storage::cache) or
(object_[address].storage != storage) or
(object_[address].name != name) or
(typeHash(object_[address].type) != typeHash<B>()) or
(typeHash(object_[address].derivedType) != typeHash<T>()))
{
HADRONS_ERROR_REF(ObjectDefinition, "object '" + name + "' already allocated", address);
}

5
Hadrons/Global.cc
View File

@@ -72,6 +72,11 @@ void Hadrons::initLogger(void)
}
// type utilities //////////////////////////////////////////////////////////////
size_t Hadrons::typeHash(const std::type_info *info)
{
return info->hash_code();
}
constexpr unsigned int maxNameSize = 1024u;
std::string Hadrons::typeName(const std::type_info *info)

52
Hadrons/Global.hpp
View File

@@ -62,18 +62,36 @@ using Grid::operator>>;
#define END_MODULE_NAMESPACE }
#ifndef FIMPL
#define FIMPL WilsonImplR
#define _HADRONS_IMPL(impl, sub) impl##sub
#define HADRONS_IMPL(impl, sub) _HADRONS_IMPL(impl, sub)
#ifndef FIMPLBASE
#define FIMPLBASE WilsonImpl
#endif
#ifndef ZFIMPL
#define ZFIMPL ZWilsonImplR
#define FIMPL HADRONS_IMPL(FIMPLBASE, R)
#define FIMPLF HADRONS_IMPL(FIMPLBASE, F)
#define FIMPLD HADRONS_IMPL(FIMPLBASE, D)
#ifndef ZFIMPLBASE
#define ZFIMPLBASE ZWilsonImpl
#endif
#ifndef SIMPL
#define SIMPL ScalarImplCR
#define ZFIMPL HADRONS_IMPL(ZFIMPLBASE, R)
#define ZFIMPLF HADRONS_IMPL(ZFIMPLBASE, F)
#define ZFIMPLD HADRONS_IMPL(ZFIMPLBASE, D)
#ifndef SIMPLBASE
#define SIMPLBASE ScalarImplC
#endif
#ifndef GIMPL
#define GIMPL PeriodicGimplR
#define SIMPL HADRONS_IMPL(SIMPLBASE, R)
#define SIMPLF HADRONS_IMPL(SIMPLBASE, F)
#define SIMPLD HADRONS_IMPL(SIMPLBASE, D)
#ifndef GIMPLBASE
#define GIMPLBASE PeriodicGimpl
#endif
#define GIMPL HADRONS_IMPL(GIMPLBASE, R)
#define GIMPLF HADRONS_IMPL(GIMPLBASE, F)
#define GIMPLD HADRONS_IMPL(GIMPLBASE, D)
BEGIN_HADRONS_NAMESPACE
@@ -155,14 +173,28 @@ const std::type_info * typeIdPt(const T &x)
return &typeid(x);
}
std::string typeName(const std::type_info *info);
template <typename T>
const std::type_info * typeIdPt(void)
{
return &typeid(T);
}
size_t typeHash(const std::type_info *info);
template <typename T>
size_t typeHash(const T &x)
{
return typeHash(typeIdPt(x));
}
template <typename T>
size_t typeHash(void)
{
return typeHash(typeIdPt<T>());
}
std::string typeName(const std::type_info *info);
template <typename T>
std::string typeName(const T &x)
{

2
Hadrons/Modules.hpp
View File

@@ -45,7 +45,9 @@ See the full license in the file "LICENSE" in the top level distribution directo
#include <Hadrons/Modules/MSource/SeqConserved.hpp>
#include <Hadrons/Modules/MSink/Smear.hpp>
#include <Hadrons/Modules/MSink/Point.hpp>
#include <Hadrons/Modules/MSolver/MixedPrecisionRBPrecCG.hpp>
#include <Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp>
#include <Hadrons/Modules/MSolver/Guesser.hpp>
#include <Hadrons/Modules/MSolver/RBPrecCG.hpp>
#include <Hadrons/Modules/MSolver/A2AVectors.hpp>
#include <Hadrons/Modules/MGauge/UnitEm.hpp>

85
Hadrons/Modules/MSolver/Guesser.hpp Normal file
View File

@@ -0,0 +1,85 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MSolver/Guesser.hpp
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_MSolver_Guesser_hpp_
#define Hadrons_MSolver_Guesser_hpp_
#include <Hadrons/Global.hpp>
#include <Hadrons/EigenPack.hpp>
BEGIN_HADRONS_NAMESPACE
BEGIN_MODULE_NAMESPACE(MSolver)
template <typename FImpl, int nBasis>
std::shared_ptr<LinearFunction<typename FImpl::FermionField>>
makeGuesser(const std::string epackName)
{
typedef typename FImpl::FermionField FermionField;
typedef FermionEigenPack<FImpl> EPack;
typedef CoarseFermionEigenPack<FImpl, nBasis> CoarseEPack;
typedef DeflatedGuesser<FermionField> FineGuesser;
typedef LocalCoherenceDeflatedGuesser<
FermionField, typename CoarseEPack::CoarseField> CoarseGuesser;
std::shared_ptr<LinearFunction<typename FImpl::FermionField>> guesserPt;
DEFINE_ENV_LAMBDA;
if (epackName.empty())
{
guesserPt.reset(new ZeroGuesser<FermionField>());
}
else
{
try
{
auto &epack = envGetDerived(EPack, CoarseEPack, epackName);
LOG(Message) << "using low-mode deflation with coarse eigenpack '"
<< epackName << "' ("
<< epack.evecCoarse.size() << " modes)" << std::endl;
guesserPt.reset(new CoarseGuesser(epack.evec, epack.evecCoarse,
epack.evalCoarse));
}
catch (Exceptions::ObjectType &e)
{
auto &epack = envGet(EPack, epackName);
LOG(Message) << "using low-mode deflation with eigenpack '"
<< epackName << "' ("
<< epack.evec.size() << " modes)" << std::endl;
guesserPt.reset(new FineGuesser(epack.evec, epack.eval));
}
}
return guesserPt;
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif

35
Hadrons/Modules/MSolver/MixedPrecisionRBPrecCG.cc Normal file
View File

@@ -0,0 +1,35 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MSolver/MixedPrecisionRBPrecCG.cc
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Hadrons/Modules/MSolver/MixedPrecisionRBPrecCG.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MSolver;
template class Grid::Hadrons::MSolver::TMixedPrecisionRBPrecCG<FIMPLF, FIMPLD, HADRONS_DEFAULT_LANCZOS_NBASIS>;
template class Grid::Hadrons::MSolver::TMixedPrecisionRBPrecCG<ZFIMPLF, ZFIMPLD, HADRONS_DEFAULT_LANCZOS_NBASIS>;

197
Hadrons/Modules/MSolver/MixedPrecisionRBPrecCG.hpp Normal file
View File

@@ -0,0 +1,197 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MSolver/MixedPrecisionRBPrecCG.hpp
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_MSolver_MixedPrecisionRBPrecCG_hpp_
#define Hadrons_MSolver_MixedPrecisionRBPrecCG_hpp_
#include <Hadrons/Global.hpp>
#include <Hadrons/Module.hpp>
#include <Hadrons/ModuleFactory.hpp>
#include <Hadrons/Solver.hpp>
#include <Hadrons/EigenPack.hpp>
#include <Hadrons/Modules/MSolver/Guesser.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* Mixed precision schur red-black preconditioned CG *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MSolver)
class MixedPrecisionRBPrecCGPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(MixedPrecisionRBPrecCGPar,
std::string , innerAction,
std::string , outerAction,
unsigned int, maxInnerIteration,
unsigned int, maxOuterIteration,
double , residual,
std::string , eigenPack);
};
template <typename FImplInner, typename FImplOuter, int nBasis>
class TMixedPrecisionRBPrecCG: public Module<MixedPrecisionRBPrecCGPar>
{
public:
FG_TYPE_ALIASES(FImplInner, Inner);
FG_TYPE_ALIASES(FImplOuter, Outer);
SOLVER_TYPE_ALIASES(FImplOuter,);
typedef HADRONS_DEFAULT_SCHUR_OP<FMatInner, FermionFieldInner> SchurFMatInner;
typedef HADRONS_DEFAULT_SCHUR_OP<FMatOuter, FermionFieldOuter> SchurFMatOuter;
private:
template <typename Field>
class OperatorFunctionWrapper: public OperatorFunction<Field>
{
public:
OperatorFunctionWrapper(LinearFunction<Field> &fn): fn_(fn) {};
virtual ~OperatorFunctionWrapper(void) = default;
virtual void operator()(LinearOperatorBase<Field> &op,
const Field &in, Field &out)
{
fn_(in, out);
}
private:
LinearFunction<Field> &fn_;
};
public:
// constructor
TMixedPrecisionRBPrecCG(const std::string name);
// destructor
virtual ~TMixedPrecisionRBPrecCG(void) {};
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getReference(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
};
MODULE_REGISTER_TMP(MixedPrecisionRBPrecCG,
ARG(TMixedPrecisionRBPrecCG<FIMPLF, FIMPLD, HADRONS_DEFAULT_LANCZOS_NBASIS>), MSolver);
MODULE_REGISTER_TMP(ZMixedPrecisionRBPrecCG,
ARG(TMixedPrecisionRBPrecCG<ZFIMPLF, ZFIMPLD, HADRONS_DEFAULT_LANCZOS_NBASIS>), MSolver);
/******************************************************************************
* TMixedPrecisionRBPrecCG implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImplInner, typename FImplOuter, int nBasis>
TMixedPrecisionRBPrecCG<FImplInner, FImplOuter, nBasis>
::TMixedPrecisionRBPrecCG(const std::string name)
: Module<MixedPrecisionRBPrecCGPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImplInner, typename FImplOuter, int nBasis>
std::vector<std::string> TMixedPrecisionRBPrecCG<FImplInner, FImplOuter, nBasis>
::getInput(void)
{
std::vector<std::string> in;
return in;
}
template <typename FImplInner, typename FImplOuter, int nBasis>
std::vector<std::string> TMixedPrecisionRBPrecCG<FImplInner, FImplOuter, nBasis>
::getReference(void)
{
std::vector<std::string> ref = {par().innerAction, par().outerAction};
if (!par().eigenPack.empty())
{
ref.push_back(par().eigenPack);
}
return ref;
}
template <typename FImplInner, typename FImplOuter, int nBasis>
std::vector<std::string> TMixedPrecisionRBPrecCG<FImplInner, FImplOuter, nBasis>
::getOutput(void)
{
std::vector<std::string> out = {getName(), getName() + "_subtract"};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImplInner, typename FImplOuter, int nBasis>
void TMixedPrecisionRBPrecCG<FImplInner, FImplOuter, nBasis>
::setup(void)
{
LOG(Message) << "Setting up Schur red-black preconditioned mixed-precision "
<< "CG for inner/outer action '" << par().innerAction
<< "'/'" << par().outerAction << "', residual "
<< par().residual << ", and maximum inner/outer iteration "
<< par().maxInnerIteration << "/" << par().maxOuterIteration
<< std::endl;
auto Ls = env().getObjectLs(par().innerAction);
auto &imat = envGet(FMatInner, par().innerAction);
auto &omat = envGet(FMatOuter, par().outerAction);
auto guesserPt = makeGuesser<FImplOuter, nBasis>(par().eigenPack);
auto makeSolver = [&imat, &omat, guesserPt, Ls, this](bool subGuess)
{
return [&imat, &omat, guesserPt, subGuess, Ls, this]
(FermionFieldOuter &sol, const FermionFieldOuter &source)
{
typedef typename FermionFieldInner::vector_type VTypeInner;
SchurFMatInner simat(imat);
SchurFMatOuter somat(omat);
MixedPrecisionConjugateGradient<FermionFieldOuter, FermionFieldInner>
mpcg(par().residual, par().maxInnerIteration,
par().maxOuterIteration,
env().template getGrid<VTypeInner>(Ls),
simat, somat);
OperatorFunctionWrapper<FermionFieldOuter> wmpcg(mpcg);
HADRONS_DEFAULT_SCHUR_SOLVE<FermionFieldOuter> schurSolver(wmpcg);
schurSolver.subtractGuess(subGuess);
schurSolver(omat, source, sol, *guesserPt);
};
};
auto solver = makeSolver(false);
envCreate(Solver, getName(), Ls, solver, omat);
auto solver_subtract = makeSolver(true);
envCreate(Solver, getName() + "_subtract", Ls, solver_subtract, omat);
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImplInner, typename FImplOuter, int nBasis>
void TMixedPrecisionRBPrecCG<FImplInner, FImplOuter, nBasis>
::execute(void)
{}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MSolver_MixedPrecisionRBPrecCG_hpp_

47
Hadrons/Modules/MSolver/RBPrecCG.hpp
View File

@@ -35,6 +35,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
#include <Hadrons/ModuleFactory.hpp>
#include <Hadrons/Solver.hpp>
#include <Hadrons/EigenPack.hpp>
#include <Hadrons/Modules/MSolver/Guesser.hpp>
BEGIN_HADRONS_NAMESPACE
@@ -59,13 +60,6 @@ class TRBPrecCG: public Module<RBPrecCGPar>
public:
FG_TYPE_ALIASES(FImpl,);
SOLVER_TYPE_ALIASES(FImpl,);
typedef FermionEigenPack<FImpl> EPack;
typedef CoarseFermionEigenPack<FImpl, nBasis> CoarseEPack;
typedef std::shared_ptr<Guesser<FermionField>> GuesserPt;
typedef DeflatedGuesser<typename FImpl::FermionField> FineGuesser;
typedef LocalCoherenceDeflatedGuesser<
typename FImpl::FermionField,
typename CoarseEPack::CoarseField> CoarseGuesser;
public:
// constructor
TRBPrecCG(const std::string name);
@@ -138,45 +132,18 @@ void TRBPrecCG<FImpl, nBasis>::setup(void)
<< par().residual << ", maximum iteration "
<< par().maxIteration << std::endl;
auto Ls = env().getObjectLs(par().action);
auto &mat = envGet(FMat, par().action);
std::string guesserName = getName() + "_guesser";
GuesserPt guesser{nullptr};
auto Ls = env().getObjectLs(par().action);
auto &mat = envGet(FMat, par().action);
auto guesserPt = makeGuesser<FImpl, nBasis>(par().eigenPack);
if (par().eigenPack.empty())
{
guesser.reset(new ZeroGuesser<FermionField>());
}
else
{
try
{
auto &epack = envGetDerived(EPack, CoarseEPack, par().eigenPack);
LOG(Message) << "using low-mode deflation with coarse eigenpack '"
<< par().eigenPack << "' ("
<< epack.evecCoarse.size() << " modes)" << std::endl;
guesser.reset(new CoarseGuesser(epack.evec, epack.evecCoarse,
epack.evalCoarse));
}
catch (Exceptions::ObjectType &e)
{
auto &epack = envGet(EPack, par().eigenPack);
LOG(Message) << "using low-mode deflation with eigenpack '"
<< par().eigenPack << "' ("
<< epack.evec.size() << " modes)" << std::endl;
guesser.reset(new FineGuesser(epack.evec, epack.eval));
}
}
auto makeSolver = [&mat, guesser, this](bool subGuess) {
return [&mat, guesser, subGuess, this](FermionField &sol,
auto makeSolver = [&mat, guesserPt, this](bool subGuess) {
return [&mat, guesserPt, subGuess, this](FermionField &sol,
const FermionField &source) {
ConjugateGradient<FermionField> cg(par().residual,
par().maxIteration);
HADRONS_DEFAULT_SCHUR_SOLVE<FermionField> schurSolver(cg);
schurSolver.subtractGuess(subGuess);
schurSolver(mat, source, sol, *guesser);
schurSolver(mat, source, sol, *guesserPt);
};
};
auto solver = makeSolver(false);

3
Hadrons/modules.inc
View File

@@ -19,6 +19,7 @@ modules_cc =\
Modules/MSink/Smear.cc \
Modules/MSolver/A2AVectors.cc \
Modules/MSolver/RBPrecCG.cc \
Modules/MSolver/MixedPrecisionRBPrecCG.cc \
Modules/MSolver/LocalCoherenceLanczos.cc \
Modules/MGauge/StoutSmearing.cc \
Modules/MGauge/Unit.cc \
@@ -80,7 +81,9 @@ modules_hpp =\
Modules/MSource/SeqConserved.hpp \
Modules/MSink/Smear.hpp \
Modules/MSink/Point.hpp \
Modules/MSolver/MixedPrecisionRBPrecCG.hpp \
Modules/MSolver/LocalCoherenceLanczos.hpp \
Modules/MSolver/Guesser.hpp \
Modules/MSolver/RBPrecCG.hpp \
Modules/MSolver/A2AVectors.hpp \
Modules/MGauge/UnitEm.hpp \

2
tests/hadrons/Test_QED.cc
View File

@@ -57,7 +57,7 @@ int main(int argc, char *argv[])
globalPar.trajCounter.start = 1500;
globalPar.trajCounter.end = 1520;
globalPar.trajCounter.step = 20;
globalPar.seed = "1 2 3 4";
globalPar.runId = "test";
application.setPar(globalPar);
// gauge field
application.createModule<MGauge::Unit>("gauge");

2
tests/hadrons/Test_free_prop.cc
View File

@@ -59,7 +59,7 @@ int main(int argc, char *argv[])
globalPar.trajCounter.start = 1500;
globalPar.trajCounter.end = 1520;
globalPar.trajCounter.step = 20;
globalPar.seed = "1 2 3 4";
globalPar.runId = "test";
application.setPar(globalPar);
// gauge field
application.createModule<MGauge::Unit>("gauge");

2
tests/hadrons/Test_hadrons.hpp
View File

@@ -51,7 +51,7 @@ using namespace Hadrons;
globalPar.trajCounter.start = 1500; \
globalPar.trajCounter.end = 1520; \
globalPar.trajCounter.step = 20; \
globalPar.seed = "1 2 3 4"; \
globalPar.runId = "test"; \
globalPar.genetic.maxGen = 1000; \
globalPar.genetic.maxCstGen = 200; \
globalPar.genetic.popSize = 20; \

2
tests/hadrons/Test_hadrons_meson_3pt.cc
View File

@@ -54,7 +54,7 @@ int main(int argc, char *argv[])
globalPar.trajCounter.start = 1500;
globalPar.trajCounter.end = 1520;
globalPar.trajCounter.step = 20;
globalPar.seed = "1 2 3 4";
globalPar.runId = "test";
globalPar.genetic.maxGen = 1000;
globalPar.genetic.maxCstGen = 200;
globalPar.genetic.popSize = 20;

156
tests/hadrons/Test_hadrons_quark.cc
View File

@@ -1,156 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Tests/Hadrons/Test_hadrons_quark.cc
Copyright (C) 2015-2018
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include "Test_hadrons.hpp"
using namespace Grid;
using namespace QCD;
using namespace Hadrons;
/*******************************************************************************
* Unit test functions within Quark module.
******************************************************************************/
// Alternative 4D & 5D projections
template<class vobj>
inline void make_4D_with_gammas(Lattice<vobj> &in_5d, Lattice<vobj> &out_4d, int Ls)
{
GridBase *_grid(out_4d._grid);
Lattice<vobj> tmp(_grid);
Gamma G5(Gamma::Algebra::Gamma5);
ExtractSlice(tmp, in_5d, 0, 0);
out_4d = 0.5 * (tmp - G5*tmp);
ExtractSlice(tmp, in_5d, Ls - 1, 0);
out_4d += 0.5 * (tmp + G5*tmp);
}
template<class vobj>
inline void make_5D_with_gammas(Lattice<vobj> &in_4d, Lattice<vobj> &out_5d, int Ls)
{
out_5d = zero;
Gamma G5(Gamma::Algebra::Gamma5);
GridBase *_grid(in_4d._grid);
Lattice<vobj> tmp(_grid);
tmp = 0.5 * (in_4d + G5*in_4d);
InsertSlice(tmp, out_5d, 0, 0);
tmp = 0.5 * (in_4d - G5*in_4d);
InsertSlice(tmp, out_5d, Ls - 1, 0);
}
int main(int argc, char **argv)
{
/***************************************************************************
* Initialisation.
**************************************************************************/
Grid_init(&argc, &argv);
std::vector<int> latt_size = GridDefaultLatt();
std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
std::vector<int> mpi_layout = GridDefaultMpi();
const int Ls = 8;
GridCartesian UGrid(latt_size,simd_layout,mpi_layout);
GridCartesian *FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls, &UGrid);
GridSerialRNG sRNG;
GridParallelRNG pRNG(&UGrid);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
GridParallelRNG rng4(&UGrid);
GridParallelRNG rng5(FGrid);
rng4.SeedFixedIntegers(seeds4);
rng5.SeedFixedIntegers(seeds5);
/***************************************************************************
* Build a 4D random source, and convert it to 5D.
**************************************************************************/
LatticeFermion test4(&UGrid);
LatticeFermion test5(FGrid);
LatticeFermion check5(FGrid);
gaussian(rng4, test4);
make_5D(test4, test5, Ls);
make_5D_with_gammas(test4, check5, Ls);
test5 -= check5;
std::cout << "4D -> 5D comparison, diff = " << Grid::sqrt(norm2(test5)) << std::endl;
/***************************************************************************
* Build a 5D random source, and project down to 4D.
**************************************************************************/
LatticeFermion check4(&UGrid);
gaussian(rng5, test5);
check5 = test5;
make_4D(test5, test4, Ls);
make_4D_with_gammas(check5, check4, Ls);
test4 -= check4;
std::cout << "5D -> 4D comparison, diff = " << Grid::sqrt(norm2(test4)) << std::endl;
/***************************************************************************
* Convert a propagator to a fermion & back.
**************************************************************************/
LatticeFermion ferm(&UGrid);
LatticePropagator prop(&UGrid), ref(&UGrid);
gaussian(rng4, prop);
// Define variables for sanity checking a single site.
typename SpinColourVector::scalar_object fermSite;
typename SpinColourMatrix::scalar_object propSite;
std::vector<int> site(Nd, 0);
for (int s = 0; s < Ns; ++s)
for (int c = 0; c < Nc; ++c)
{
ref = prop;
PropToFerm<WilsonImplR>(ferm, prop, s, c);
FermToProp<WilsonImplR>(prop, ferm, s, c);
std::cout << "Spin = " << s << ", Colour = " << c << std::endl;
ref -= prop;
std::cout << "Prop->Ferm->Prop test, diff = " << Grid::sqrt(norm2(ref)) << std::endl;
peekSite(fermSite, ferm, site);
peekSite(propSite, prop, site);
for (int s2 = 0; s2 < Ns; ++s2)
for (int c2 = 0; c2 < Nc; ++c2)
{
if (propSite()(s2, s)(c2, c) != fermSite()(s2)(c2))
{
std::cout << propSite()(s2, s)(c2, c) << " != "
<< fermSite()(s2)(c2) << " for spin = " << s2
<< ", col = " << c2 << std::endl;
}
}
}
Grid_finalize();
return EXIT_SUCCESS;
}

2
tests/hadrons/Test_hadrons_spectrum.cc
View File

@@ -53,7 +53,7 @@ int main(int argc, char *argv[])
globalPar.trajCounter.start = 1500;
globalPar.trajCounter.end = 1520;
globalPar.trajCounter.step = 20;
globalPar.seed = "1 2 3 4";
globalPar.runId = "test";
application.setPar(globalPar);
// gauge field
application.createModule<MGauge::Unit>("gauge");

3
tests/hadrons/Test_hadrons_wilsonFund.cc
View File

@@ -55,8 +55,7 @@ int main(int argc, char *argv[])
globalPar.trajCounter.start = 309;
globalPar.trajCounter.end = 310;
globalPar.trajCounter.step = 1;
globalPar.seed = "1 2 3 4";
globalPar.runId = "test";
application.setPar(globalPar);
// gauge field