portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-09 23:45:36 +00:00
Code Releases Activity

merged in develop

Browse Source
This commit is contained in:
ferben 2019-12-04 17:12:46 +00:00
commit cd9fd80a5d
39 changed files with 2641 additions and 68 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
Grid/communicator/SharedMemoryMPI.cc
View File

@ -162,11 +162,8 @@ static inline int divides(int a,int b)
void GlobalSharedMemory::GetShmDims(const Coordinate &WorldDims,Coordinate &ShmDims)
{
////////////////////////////////////////////////////////////////
// Assert power of two shm_size.
// Powers of 2,3,5 only in prime decomposition for now
////////////////////////////////////////////////////////////////
int log2size = Log2Size(WorldShmSize,MAXLOG2RANKSPERNODE);
assert(log2size != -1);
int ndimension = WorldDims.size();
ShmDims=Coordinate(ndimension,1);
@ -177,7 +174,8 @@ void GlobalSharedMemory::GetShmDims(const Coordinate &WorldDims,Coordinate &ShmD
while(AutoShmSize != WorldShmSize) {
for(int p=0;p<primes.size();p++) {
int prime=primes[p];
if ( divides(prime,WorldDims[dim]/ShmDims[dim]) ) {
if ( divides(prime,WorldDims[dim]/ShmDims[dim])
&& divides(prime,WorldShmSize/AutoShmSize) ) {
AutoShmSize*=prime;
ShmDims[dim]*=prime;
break;
@ -308,7 +306,6 @@ void GlobalSharedMemory::OptimalCommunicatorHypercube(const Coordinate &processo
}
void GlobalSharedMemory::OptimalCommunicatorSharedMemory(const Coordinate &processors,Grid_MPI_Comm & optimal_comm)
{
////////////////////////////////////////////////////////////////
// Identify subblock of ranks on node spreading across dims
// in a maximally symmetrical way
@ -435,10 +432,13 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
// e.g. DGX1, supermicro board,
//////////////////////////////////////////////////////////////////////////////////////////////////////////
// cudaDeviceGetP2PAttribute(&perfRank, cudaDevP2PAttrPerformanceRank, device1, device2);
#ifdef GRID_IBM_SUMMIT
std::cout << header << "flag IBM_SUMMIT disabled CUDA set device: ensure jsrun is used correctly" <<std::endl;
// IBM Jsrun makes cuda Device numbering screwy and not match rank
std::cout << "IBM Summit or similar - NOT setting device to WorldShmRank"<<std::endl;
#else
cudaSetDevice(WorldShmRank);
std::cout << "setting device to WorldShmRank"<<std::endl;
cudaSetDevice(WorldShmRank);
#endif
///////////////////////////////////////////////////////////////////////////////////////////////////////////
// Each MPI rank should allocate our own buffer
@ -466,7 +466,7 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
// If it is me, pass around the IPC access key
//////////////////////////////////////////////////
cudaIpcMemHandle_t handle;
if ( r==WorldShmRank ) {
err = cudaIpcGetMemHandle(&handle,ShmCommBuf);
if ( err != cudaSuccess) {
@ -735,6 +735,24 @@ void SharedMemory::SetCommunicator(Grid_MPI_Comm comm)
std::vector<int> ranks(size); for(int r=0;r<size;r++) ranks[r]=r;
MPI_Group_translate_ranks (FullGroup,size,&ranks[0],ShmGroup, &ShmRanks[0]);
#ifdef GRID_IBM_SUMMIT
// Hide the shared memory path between sockets
// if even number of nodes
if ( (ShmSize & 0x1)==0 ) {
int SocketSize = ShmSize/2;
int mySocket = ShmRank/SocketSize;
for(int r=0;r<size;r++){
int hisRank=ShmRanks[r];
if ( hisRank!= MPI_UNDEFINED ) {
int hisSocket=hisRank/SocketSize;
if ( hisSocket != mySocket ) {
ShmRanks[r] = MPI_UNDEFINED;
}
}
}
}
#endif
SharedMemoryTest();
}
//////////////////////////////////////////////////////////////////

10
Grid/perfmon/PerfCount.h
View File

@ -44,8 +44,13 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
#include <sys/syscall.h>
#endif
#ifdef __x86_64__
#ifdef GRID_NVCC
accelerator_inline uint64_t __rdtsc(void) { return 0; }
accelerator_inline uint64_t __rdpmc(int ) { return 0; }
#else
#include <x86intrin.h>
#endif
#endif
NAMESPACE_BEGIN(Grid);
@ -89,13 +94,8 @@ inline uint64_t cyclecount(void){
return tmp;
}
#elif defined __x86_64__
#ifdef GRID_NVCC
accelerator_inline uint64_t __rdtsc(void) { return 0; }
#endif
inline uint64_t cyclecount(void){
return __rdtsc();
// unsigned int dummy;
// return __rdtscp(&dummy);
}
#else

2
Grid/qcd/action/fermion/instantiation/ImprovedStaggeredFermionInstantiation.cc
View File

@ -26,7 +26,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include <Grid.h>
#include <Grid/Grid.h>
NAMESPACE_BEGIN(Grid);

2
Grid/qcd/action/fermion/instantiation/ImprovedStaggeredFermionInstantiation.cc.master
View File

@ -26,7 +26,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include <Grid.h>
#include <Grid/Grid.h>
#include <Grid/qcd/action/fermion/implementation/ImprovedStaggeredFermionImplementation.h>
NAMESPACE_BEGIN(Grid);

2
Grid/qcd/action/fermion/instantiation/StaggeredImplD/ImprovedStaggeredFermionInstantiationStaggeredImplD.cc
View File

@ -26,7 +26,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include <Grid.h>
#include <Grid/Grid.h>
#include <Grid/qcd/action/fermion/implementation/ImprovedStaggeredFermionImplementation.h>
NAMESPACE_BEGIN(Grid);

2
Grid/qcd/action/fermion/instantiation/StaggeredImplF/ImprovedStaggeredFermionInstantiationStaggeredImplF.cc
View File

@ -26,7 +26,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include <Grid.h>
#include <Grid/Grid.h>
#include <Grid/qcd/action/fermion/implementation/ImprovedStaggeredFermionImplementation.h>
NAMESPACE_BEGIN(Grid);

92
Grid/qcd/smearing/StoutSmearing.h
View File

@ -1,5 +1,34 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/smearing/StoutSmearing.h
Copyright (C) 2019
Author: unknown
Author: Felix Erben <ferben@ed.ac.uk>
Author: Michael Marshall <Michael.Marshall@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/*
@file stoutSmear.hpp
@file StoutSmearing.h
@brief Declares Stout smearing class
*/
#pragma once
@ -9,19 +38,43 @@ NAMESPACE_BEGIN(Grid);
/*! @brief Stout smearing of link variable. */
template <class Gimpl>
class Smear_Stout : public Smear<Gimpl> {
private:
const Smear<Gimpl>* SmearBase;
private:
int OrthogDim = -1;
const std::vector<double> SmearRho;
// Smear<Gimpl>* ownership semantics:
// Smear<Gimpl>* passed in to constructor are owned by caller, so we don't delete them here
// Smear<Gimpl>* created within constructor need to be deleted as part of the destructor
const std::unique_ptr<Smear<Gimpl>> OwnedBase; // deleted at destruction
const Smear<Gimpl>* SmearBase; // Not owned by this object, so not deleted at destruction
// only anticipated to be used from default constructor
inline static std::vector<double> rho3D(double rho, int orthogdim){
std::vector<double> rho3d(Nd*Nd);
for (int mu=0; mu<Nd; mu++)
for (int nu=0; nu<Nd; nu++)
rho3d[mu + Nd * nu] = (mu == nu || mu == orthogdim || nu == orthogdim) ? 0.0 : rho;
return rho3d;
};
public:
INHERIT_GIMPL_TYPES(Gimpl)
Smear_Stout(Smear<Gimpl>* base) : SmearBase(base) {
assert(Nc == 3);// "Stout smearing currently implemented only for Nc==3");
/*! Stout smearing with base explicitly specified */
Smear_Stout(Smear<Gimpl>* base) : SmearBase{base} {
assert(Nc == 3 && "Stout smearing currently implemented only for Nc==3");
}
/*! Default constructor */
Smear_Stout(double rho = 1.0) : SmearBase(new Smear_APE<Gimpl>(rho)) {
assert(Nc == 3);// "Stout smearing currently implemented only for Nc==3");
/*! Construct stout smearing object from explicitly specified rho matrix */
Smear_Stout(const std::vector<double>& rho_)
: OwnedBase{new Smear_APE<Gimpl>(rho_)}, SmearBase{OwnedBase.get()} {
std::cout << GridLogDebug << "Stout smearing constructor : Smear_Stout(const std::vector<double>& " << rho_ << " )" << std::endl
assert(Nc == 3 && "Stout smearing currently implemented only for Nc==3");
}
/*! Default constructor. rho is constant in all directions, optionally except for orthogonal dimension */
Smear_Stout(double rho = 1.0, int orthogdim = -1)
: OrthogDim{orthogdim}, SmearRho{ rho3D(rho,orthogdim) }, OwnedBase{ new Smear_APE<Gimpl>(SmearRho) }, SmearBase{OwnedBase.get()} {
assert(Nc == 3 && "Stout smearing currently implemented only for Nc==3");
}
~Smear_Stout() {} // delete SmearBase...
@ -36,12 +89,16 @@ public:
SmearBase->smear(C, U);
for (int mu = 0; mu < Nd; mu++) {
tmp = peekLorentz(C, mu);
Umu = peekLorentz(U, mu);
iq_mu = Ta(
tmp *
adj(Umu)); // iq_mu = Ta(Omega_mu) to match the signs with the paper
exponentiate_iQ(tmp, iq_mu);
if( mu == OrthogDim )
tmp = 1.0; // Don't smear in the orthogonal direction
else {
tmp = peekLorentz(C, mu);
Umu = peekLorentz(U, mu);
iq_mu = Ta(
tmp *
adj(Umu)); // iq_mu = Ta(Omega_mu) to match the signs with the paper
exponentiate_iQ(tmp, iq_mu);
}
pokeLorentz(u_smr, tmp * Umu, mu); // u_smr = exp(iQ_mu)*U_mu
}
std::cout << GridLogDebug << "Stout smearing completed\n";
@ -80,6 +137,7 @@ public:
iQ2 = iQ * iQ;
iQ3 = iQ * iQ2;
//We should check sgn(c0) here already and then apply eq (34) from 0311018
set_uw(u, w, iQ2, iQ3);
set_fj(f0, f1, f2, u, w);
@ -139,9 +197,8 @@ public:
}
LatticeComplex func_xi0(const LatticeComplex& w) const {
// Define a function to do the check
// if( w < 1e-4 ) std::cout << GridLogWarning<< "[Smear_stout] w too small:
// "<< w <<"\n";
// Definition from arxiv 0311018
//if (abs(w) < 0.05) {w2 = w*w; return 1.0 - w2/6.0 * (1.0-w2/20.0 * (1.0-w2/42.0));}
return sin(w) / w;
}
@ -154,4 +211,3 @@ public:
};
NAMESPACE_END(Grid);

2
Grid/qcd/utils/A2Autils.h
View File

@ -260,7 +260,7 @@ void A2Autils<FImpl>::MesonField(TensorType &mat,
int ij_dx = m+Nmom*i + Nmom*Lblock * j + Nmom*Lblock * Rblock * lt;
for(int mu=0;mu<Ngamma;mu++){
// this is a bit slow
mat(m,mu,t,i,j) = trace(lsSum[ij_dx]*Gamma(gammas[mu]));
mat(m,mu,t,i,j) = trace(lsSum[ij_dx]*Gamma(gammas[mu]))()()();
}
}
}

39
Grid/qcd/utils/BaryonUtils.h
View File

@ -45,8 +45,8 @@ public:
typedef Lattice<iSpinMatrix<typename FImpl::Simd>> SpinMatrixField;
typedef typename SpinMatrixField::vector_object sobj;
static constexpr int epsilon[6][3] = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}};
static constexpr int epsilon_sgn[6]= {1,1,1,-1,-1,-1};
static const int epsilon[6][3] ;
static const Complex epsilon_sgn[6];
private:
template <class mobj, class robj>
@ -149,10 +149,15 @@ public:
SpinMatrixField &stn_corr);
};
template <class FImpl>
constexpr int BaryonUtils<FImpl>::epsilon[6][3];
template <class FImpl>
constexpr int BaryonUtils<FImpl>::epsilon_sgn[6];
template <class FImpl>
const int BaryonUtils<FImpl>::epsilon[6][3] = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}};
template <class FImpl>
const Complex BaryonUtils<FImpl>::epsilon_sgn[6] = {Complex(1),
Complex(1),
Complex(1),
Complex(-1),
Complex(-1),
Complex(-1)};
template <class FImpl>
template <class mobj, class robj>
@ -189,7 +194,7 @@ void BaryonUtils<FImpl>::baryon_site(const mobj &D1,
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
for (int beta_left=0; beta_left<Ns; beta_left++){
for (int gamma_left=0; gamma_left<Ns; gamma_left++){
result()()() += static_cast<Complex>(epsilon_sgn[ie_left] * epsilon_sgn[ie_right]) * pD1()(gamma_left,gamma_left)(c_right,c_left)*D2g()(alpha_right,beta_left)(a_right,a_left)*gD3()(alpha_right,beta_left)(b_right,b_left);
result()()() += epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1()(gamma_left,gamma_left)(c_right,c_left)*D2g()(alpha_right,beta_left)(a_right,a_left)*gD3()(alpha_right,beta_left)(b_right,b_left);
}}}
}
//This is the \delta_{456}^{231} part
@ -198,7 +203,7 @@ void BaryonUtils<FImpl>::baryon_site(const mobj &D1,
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
for (int beta_left=0; beta_left<Ns; beta_left++){
for (int gamma_left=0; gamma_left<Ns; gamma_left++){
result()()() += static_cast<Complex>(epsilon_sgn[ie_left] * epsilon_sgn[ie_right]) * pD1g()(gamma_left,beta_left)(c_right,a_left)*D2()(alpha_right,beta_left)(a_right,b_left)*gD3()(alpha_right,gamma_left)(b_right,c_left);
result()()() += epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1g()(gamma_left,beta_left)(c_right,a_left)*D2()(alpha_right,beta_left)(a_right,b_left)*gD3()(alpha_right,gamma_left)(b_right,c_left);
}}}
}
//This is the \delta_{456}^{312} part
@ -207,7 +212,7 @@ void BaryonUtils<FImpl>::baryon_site(const mobj &D1,
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
for (int beta_left=0; beta_left<Ns; beta_left++){
for (int gamma_left=0; gamma_left<Ns; gamma_left++){
result()()() += static_cast<Complex>(epsilon_sgn[ie_left] * epsilon_sgn[ie_right]) * pD1()(gamma_left,beta_left)(c_right,b_left)*D2()(alpha_right,gamma_left)(a_right,c_left)*gD3g()(alpha_right,beta_left)(b_right,a_left);
result()()() += epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1()(gamma_left,beta_left)(c_right,b_left)*D2()(alpha_right,gamma_left)(a_right,c_left)*gD3g()(alpha_right,beta_left)(b_right,a_left);
}}}
}
//This is the \delta_{456}^{132} part
@ -216,7 +221,7 @@ void BaryonUtils<FImpl>::baryon_site(const mobj &D1,
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
for (int beta_left=0; beta_left<Ns; beta_left++){
for (int gamma_left=0; gamma_left<Ns; gamma_left++){
result()()() -= static_cast<Complex>(epsilon_sgn[ie_left] * epsilon_sgn[ie_right]) * pD1()(gamma_left,gamma_left)(c_right,c_left)*D2()(alpha_right,beta_left)(a_right,b_left)*gD3g()(alpha_right,beta_left)(b_right,a_left);
result()()() -= epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1()(gamma_left,gamma_left)(c_right,c_left)*D2()(alpha_right,beta_left)(a_right,b_left)*gD3g()(alpha_right,beta_left)(b_right,a_left);
}}}
}
//This is the \delta_{456}^{321} part
@ -225,7 +230,7 @@ void BaryonUtils<FImpl>::baryon_site(const mobj &D1,
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
for (int beta_left=0; beta_left<Ns; beta_left++){
for (int gamma_left=0; gamma_left<Ns; gamma_left++){
result()()() -= static_cast<Complex>(epsilon_sgn[ie_left] * epsilon_sgn[ie_right]) * pD1()(gamma_left,beta_left)(c_right,b_left)*D2g()(alpha_right,beta_left)(a_right,a_left)*gD3()(alpha_right,gamma_left)(b_right,c_left);
result()()() -= epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1()(gamma_left,beta_left)(c_right,b_left)*D2g()(alpha_right,beta_left)(a_right,a_left)*gD3()(alpha_right,gamma_left)(b_right,c_left);
}}}
}
//This is the \delta_{456}^{213} part
@ -234,7 +239,7 @@ void BaryonUtils<FImpl>::baryon_site(const mobj &D1,
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
for (int beta_left=0; beta_left<Ns; beta_left++){
for (int gamma_left=0; gamma_left<Ns; gamma_left++){
result()()() -= static_cast<Complex>(epsilon_sgn[ie_left] * epsilon_sgn[ie_right]) * pD1g()(gamma_left,beta_left)(c_right,a_left)*D2()(alpha_right,gamma_left)(a_right,c_left)*gD3()(alpha_right,beta_left)(b_right,b_left);
result()()() -= epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1g()(gamma_left,beta_left)(c_right,a_left)*D2()(alpha_right,gamma_left)(a_right,c_left)*gD3()(alpha_right,beta_left)(b_right,b_left);
}}}
}
}
@ -414,10 +419,10 @@ void BaryonUtils<FImpl>::Sigma_to_Nucleon_Q1_NonEye_site(const mobj &Du_xi,
for (int tau2=0; tau2<Ns; tau2++){
for (int j=0; j<Nc; j++){
auto DuGH_at_aj = DuGH()(alpha_s,tau2)(a_s,j);
auto ee_GDGDDG_a = static_cast<Complex>(epsilon_sgn[ie_n] * epsilon_sgn[ie_s]) * GDsGDd_ab_bb * DuGH_at_aj;
auto ee_GDGDDG_a = epsilon_sgn[ie_n] * epsilon_sgn[ie_s] * GDsGDd_ab_bb * DuGH_at_aj;
for (int gamma_s=0; gamma_s<Ns; gamma_s++){
auto DuGH_gt_cj = DuGH()(gamma_s,tau2)(c_s,j);
auto ee_GDGDDG_c = static_cast<Complex>(epsilon_sgn[ie_n] * epsilon_sgn[ie_s]) * GDsGDd_ab_bb * DuGH_gt_cj;
auto ee_GDGDDG_c = epsilon_sgn[ie_n] * epsilon_sgn[ie_s] * GDsGDd_ab_bb * DuGH_gt_cj;
for (int gamma_n=0; gamma_n<Ns; gamma_n++){
result()(gamma_s,gamma_n)() += ee_GDGDDG_a * DuG()(gamma_s, beta_n)(c_s,a_n) * adjDu()(tau2,gamma_n)(j,c_n);
result()(gamma_s,gamma_n)() += ee_GDGDDG_c * adjDuG()(tau2, beta_n)(j,a_n) * Du_spec()(alpha_s,gamma_n)(a_s,c_n);
@ -466,7 +471,7 @@ void BaryonUtils<FImpl>::Sigma_to_Nucleon_Q2_Eye_site(const mobj &Dq_loop,
auto GDsG_at_bi = GDsG()(alpha_s,tau)(b_s,i);
for (int beta_n=0; beta_n<Ns; beta_n++){
auto DqGDd_tb_ib = DqGDd()(tau,beta_n)(i,b_n);
auto ee_GDGDGD = static_cast<Complex>(epsilon_sgn[ie_n] * epsilon_sgn[ie_s]) * GDsG_at_bi * DqGDd_tb_ib;
auto ee_GDGDGD = epsilon_sgn[ie_n] * epsilon_sgn[ie_s] * GDsG_at_bi * DqGDd_tb_ib;
for (int gamma_s=0; gamma_s<Ns; gamma_s++){
for (int gamma_n=0; gamma_n<Ns; gamma_n++){
result()(gamma_s,gamma_n)() -= ee_GDGDGD * DuG()(alpha_s, beta_n)(a_s,a_n) * Du_spec()(gamma_s,gamma_n)(c_s,c_n);
@ -515,10 +520,10 @@ void BaryonUtils<FImpl>::Sigma_to_Nucleon_Q2_NonEye_site(const mobj &Du_xi,
auto GDsG_at_bi = GDsG()(alpha_s,tau)(b_s,i);
for (int beta_n=0; beta_n<Ns; beta_n++){
auto DuGDd_ab_ab = DuGDd()(alpha_s,beta_n)(a_s,b_n);
auto ee_GDGDGD_a = static_cast<Complex>(epsilon_sgn[ie_n] * epsilon_sgn[ie_s]) * GDsG_at_bi * DuGDd_ab_ab;
auto ee_GDGDGD_a = epsilon_sgn[ie_n] * epsilon_sgn[ie_s] * GDsG_at_bi * DuGDd_ab_ab;
for (int gamma_s=0; gamma_s<Ns; gamma_s++){
auto DuGDd_gb_cb = DuGDd()(gamma_s,beta_n)(c_s,b_n);
auto ee_GDGDGD_c = static_cast<Complex>(epsilon_sgn[ie_n] * epsilon_sgn[ie_s]) * GDsG_at_bi * DuGDd_gb_cb;
auto ee_GDGDGD_c = epsilon_sgn[ie_n] * epsilon_sgn[ie_s] * GDsG_at_bi * DuGDd_gb_cb;
for (int gamma_n=0; gamma_n<Ns; gamma_n++){
result()(gamma_s,gamma_n)() -= ee_GDGDGD_a * DuG()(gamma_s, beta_n)(c_s,a_n) * adjDu()(tau,gamma_n)(i,c_n);
result()(gamma_s,gamma_n)() -= ee_GDGDGD_c * adjDuG()(tau, beta_n)(i,a_n) * Du_spec()(alpha_s,gamma_n)(a_s,c_n);

4
Grid/stencil/Stencil.h
View File

@ -1233,7 +1233,7 @@ public:
};
void Report(void) {
#define AVERAGE(A) _grid->GlobalSum(A);A/=NP;
#define AVERAGE(A)
#define PRINTIT(A) AVERAGE(A); std::cout << GridLogMessage << " Stencil " << #A << " "<< A/calls<<std::endl;
RealD NP = _grid->_Nprocessors;
RealD NN = _grid->NodeCount();
@ -1281,11 +1281,13 @@ public:
std::cout << GridLogMessage << " Stencil SHM mem " << (membytes)/gatheralltime/1000. << " GB/s per rank"<<std::endl;
std::cout << GridLogMessage << " Stencil SHM mem " << (membytes)/gatheralltime/1000.*NP/NN << " GB/s per node"<<std::endl;
}
/*
PRINTIT(mpi3synctime);
PRINTIT(mpi3synctime_g);
PRINTIT(shmmergetime);
PRINTIT(splicetime);
PRINTIT(nosplicetime);
*/
}
#undef PRINTIT
#undef AVERAGE

10
Hadrons/Environment.cc
View File

@ -84,6 +84,16 @@ GridParallelRNG * Environment::get4dRng(void)
return rng4d_.get();
}
GridSerialRNG * Environment::getSerialRng(void)
{
if (rngSerial_ == nullptr)
{
rngSerial_.reset(new GridSerialRNG());
}
return rngSerial_.get();
}
// general memory management ///////////////////////////////////////////////////
void Environment::addObject(const std::string name, const int moduleAddress)
{

3
Hadrons/Environment.hpp
View File

@ -74,6 +74,7 @@ public:
typedef std::unique_ptr<GridCartesian> GridPt;
typedef std::unique_ptr<GridRedBlackCartesian> GridRbPt;
typedef std::unique_ptr<GridParallelRNG> RngPt;
typedef std::unique_ptr<GridSerialRNG> SerialRngPt;
enum class Storage {object, cache, temporary};
private:
struct ObjInfo
@ -114,6 +115,7 @@ public:
double getVolume(void) const;
// random number generator
GridParallelRNG * get4dRng(void);
GridSerialRNG * getSerialRng(void);
// general memory management
void addObject(const std::string name,
const int moduleAddress = -1);
@ -183,6 +185,7 @@ private:
unsigned int nd_;
// random number generator
RngPt rng4d_{nullptr};
SerialRngPt rngSerial_{nullptr};
// object store
std::vector<ObjInfo> object_;
std::map<std::string, unsigned int> objectAddress_;

15
Hadrons/Module.cc
View File

@ -93,3 +93,18 @@ GridParallelRNG & ModuleBase::rng4d(void)
return r;
}
GridSerialRNG & ModuleBase::rngSerial(void)
{
auto &r = *env().getSerialRng();
if (makeSeedString() != seed_)
{
seed_ = makeSeedString();
LOG(Message) << "Seeding Serial RNG " << &r << " with string '"
<< seed_ << "'" << std::endl;
r.SeedUniqueString(seed_);
}
return r;
}

2
Hadrons/Module.hpp
View File

@ -1,7 +1,6 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Module.hpp
Copyright (C) 2015-2019
@ -196,6 +195,7 @@ protected:
DEFINE_VM_ALIAS;
// RNG seeded from module string
GridParallelRNG &rng4d(void);
GridSerialRNG &rngSerial(void);
private:
std::string makeSeedString(void);
private:

9
Hadrons/Modules.hpp
View File

@ -50,8 +50,10 @@
#include <Hadrons/Modules/MIO/LoadBinary.hpp>
#include <Hadrons/Modules/MIO/LoadNersc.hpp>
#include <Hadrons/Modules/MIO/LoadEigenPack.hpp>
#include <Hadrons/Modules/MIO/LoadDistilNoise.hpp>
#include <Hadrons/Modules/MIO/LoadA2AMatrixDiskVector.hpp>
#include <Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp>
#include <Hadrons/Modules/MIO/LoadPerambulator.hpp>
#include <Hadrons/Modules/MIO/LoadA2AVectors.hpp>
#include <Hadrons/Modules/MIO/LoadCosmHol.hpp>
#include <Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp>
@ -65,6 +67,13 @@
#include <Hadrons/Modules/MScalarSUN/StochFreeField.hpp>
#include <Hadrons/Modules/MScalarSUN/Div.hpp>
#include <Hadrons/Modules/MScalarSUN/TrKinetic.hpp>
#include <Hadrons/Modules/MDistil/PerambFromSolve.hpp>
#include <Hadrons/Modules/MDistil/Distil.hpp>
#include <Hadrons/Modules/MDistil/DistilVectors.hpp>
#include <Hadrons/Modules/MDistil/Noises.hpp>
#include <Hadrons/Modules/MDistil/LapEvec.hpp>
#include <Hadrons/Modules/MDistil/Perambulator.hpp>
#include <Hadrons/Modules/MDistil/DistilPar.hpp>
#include <Hadrons/Modules/MAction/ZMobiusDWF.hpp>
#include <Hadrons/Modules/MAction/ScaledDWF.hpp>
#include <Hadrons/Modules/MAction/WilsonClover.hpp>

2
Hadrons/Modules/MContraction/Meson.hpp
View File

@ -199,7 +199,7 @@ void TMeson<FImpl1, FImpl2>::execute(void)
Gamma gSnk(gammaList[i].first);
Gamma gSrc(gammaList[i].second);
for (unsigned int t = 0; t < buf.size(); ++t)
for (unsigned int t = 0; t < nt; ++t)
{
result[i].corr[t] = TensorRemove(trace(mesonConnected(q1[t], q2[t], gSnk, gSrc)));
}

124
Hadrons/Modules/MDistil/Distil.hpp Normal file
View File

@ -0,0 +1,124 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MDistil/Distil.hpp
Copyright (C) 2015-2019
Author: Felix Erben <ferben@ed.ac.uk>
Author: Michael Marshall <Michael.Marshall@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_MDistil_Distil_hpp_
#define Hadrons_MDistil_Distil_hpp_
#include <Hadrons/NamedTensor.hpp>
#include <Hadrons/Module.hpp>
#include <Hadrons/ModuleFactory.hpp>
#include <Hadrons/Solver.hpp>
#include <Hadrons/A2AVectors.hpp>
#include <Hadrons/DilutedNoise.hpp>
BEGIN_HADRONS_NAMESPACE
BEGIN_MODULE_NAMESPACE(MDistil)
/******************************************************************************
Distillation code that is common across modules
Documentation on how to use this code available at
* https://aportelli.github.io/Hadrons-doc/#/mdistil *
Notation for (stochastic) DistilParameters taken from 1104.3870:
TI is interlaced dilution in time (corresponding to Nt = time-dimension of the lattice)
LI is interlaced dilution in laplacian-eigenvector space (corresponding to nvec)
SI is interlaced dilution in spin (corresponding to Ns, taken from Grid, usually Ns=4)
This code automatically computes perambulators using exact distillation if
* (TI,LI,SI) = (Nt,nvec,Ns) *
In this case, nnoise=1 and Noises is set to an array of values =1 as well.
tsrc then specifies the only timeslice on which the sources are supported.
(( for stochastic distillation, the vaue of tsrc has no meaning in this code ))
******************************************************************************/
struct DistilParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(DistilParameters,
int, nvec,
int, nnoise,
int, tsrc,
int, TI,
int, LI,
int, SI )
};
/******************************************************************************
Make a lower dimensional grid in preparation for local slice operations
******************************************************************************/
inline void MakeLowerDimGrid( std::unique_ptr<GridCartesian> &up, GridCartesian * gridHD )
{
int nd{static_cast<int>(gridHD->_ndimension)};
Coordinate latt_size = gridHD->_gdimensions;
latt_size[nd-1] = 1;
Coordinate simd_layout = GridDefaultSimd(nd-1, vComplex::Nsimd());
simd_layout.push_back( 1 );
Coordinate mpi_layout = gridHD->_processors;
mpi_layout[nd-1] = 1;
up.reset( new GridCartesian(latt_size,simd_layout,mpi_layout,*gridHD) );
}
/*************************************************************************************
Rotate eigenvectors into our phase convention
First component of first eigenvector is real and positive
*************************************************************************************/
inline void RotateEigen(std::vector<LatticeColourVector> & evec)
{
ColourVector cv0;
auto grid = evec[0].Grid();
Coordinate siteFirst(grid->Nd(),0);
peekSite(cv0, evec[0], siteFirst);
const std::complex<Real> cplx0{cv0()()(0).real(), cv0()()(0).imag()};
if( cplx0.imag() == 0 )
LOG(Message) << "RotateEigen() : Site 0 : " << cplx0 << " => already meets phase convention" << std::endl;
else
{
const Real cplx0_mag{ std::abs(cplx0) };
const std::complex<Real> std_phase{std::conj(cplx0/cplx0_mag)};
LOG(Message) << "RotateEigen() : Site 0 : |" << cplx0 << "|=" << cplx0_mag
<< " => phase=" << (std::arg(std_phase) / M_PI) << " pi" << std::endl;
{
const Grid::Complex phase{std_phase.real(),std_phase.imag()};
for( int k = 0 ; k < evec.size() ; k++ )
evec[k] *= phase;
// Get rid of the rounding error in imaginary phase on the very first site
peekSite(cv0, evec[0], siteFirst);
cv0()()(0).imag(0); // this should be zero after the phase multiply - force it to be so
pokeSite(cv0, evec[0], siteFirst);
}
}
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif

7
Hadrons/Modules/MDistil/DistilPar.cc Normal file
View File

@ -0,0 +1,7 @@
#include <Hadrons/Modules/MDistil/DistilPar.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MDistil;
template class Grid::Hadrons::MDistil::TDistilPar<FIMPL>;

97
Hadrons/Modules/MDistil/DistilPar.hpp Normal file
View File

@ -0,0 +1,97 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MDistil/DistilPar.hpp
Copyright (C) 2019
Author: Felix Erben <ferben@ed.ac.uk>
Author: Michael Marshall <Michael.Marshall@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_MDistil_DistilPar_hpp_
#define Hadrons_MDistil_DistilPar_hpp_
#include <Hadrons/Modules/MDistil/Distil.hpp>
BEGIN_HADRONS_NAMESPACE
BEGIN_MODULE_NAMESPACE(MDistil)
/******************************************************************************
* DistilPar *
******************************************************************************/
template <typename FImpl>
class TDistilPar: public Module<DistilParameters>
{
public:
// constructor
TDistilPar(const std::string name);
// destructor
virtual ~TDistilPar(void) {};
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
};
MODULE_REGISTER_TMP(DistilPar, TDistilPar<FIMPL>, MDistil);
/******************************************************************************
* TDistilPar implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TDistilPar<FImpl>::TDistilPar(const std::string name) : Module<DistilParameters>(name) {}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TDistilPar<FImpl>::getInput(void)
{
return {};
}
template <typename FImpl>
std::vector<std::string> TDistilPar<FImpl>::getOutput(void)
{
return {getName()};
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TDistilPar<FImpl>::setup(void)
{
envCreate(DistilParameters, getName(), 1, par() );
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TDistilPar<FImpl>::execute(void)
{
// Nothing to do. setup() created and initialised the output object
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif

36
Hadrons/Modules/MDistil/DistilVectors.cc Normal file
View File

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

243
Hadrons/Modules/MDistil/DistilVectors.hpp Normal file
View File

@ -0,0 +1,243 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MDistil/DistilVectors.hpp
Copyright (C) 2019
Author: Felix Erben <ferben@ed.ac.uk>
Author: Michael Marshall <Michael.Marshall@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_MDistil_DistilVectors_hpp_
#define Hadrons_MDistil_DistilVectors_hpp_
#include <Hadrons/Modules/MDistil/Distil.hpp>
BEGIN_HADRONS_NAMESPACE
BEGIN_MODULE_NAMESPACE(MDistil)
/******************************************************************************
* DistilVectors *
* (Create rho and/or phi vectors) *
******************************************************************************/
class DistilVectorsPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(DistilVectorsPar,
std::string, noise,
std::string, perambulator,
std::string, lapevec,
std::string, rho,
std::string, phi,
std::string, DistilParams);
};
template <typename FImpl>
class TDistilVectors: public Module<DistilVectorsPar>
{
public:
FERM_TYPE_ALIASES(FImpl,);
// constructor
TDistilVectors(const std::string name);
// destructor
virtual ~TDistilVectors(void) {};
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
protected:
std::unique_ptr<GridCartesian> grid3d; // Owned by me, so I must delete it
public:
// These variables contain parameters
std::string RhoName;
std::string PhiName;
};
MODULE_REGISTER_TMP(DistilVectors, TDistilVectors<FIMPL>, MDistil);
/******************************************************************************
* TDistilVectors implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TDistilVectors<FImpl>::TDistilVectors(const std::string name) : Module<DistilVectorsPar>(name) {}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TDistilVectors<FImpl>::getInput(void)
{
return {par().noise,par().perambulator,par().lapevec,par().DistilParams};
}
template <typename FImpl>
std::vector<std::string> TDistilVectors<FImpl>::getOutput(void)
{
RhoName = par().rho;
PhiName = par().phi;
if (RhoName.empty() && PhiName.empty())
{
HADRONS_ERROR(Argument,"No output specified");
}
std::vector<std::string> out;
if (!RhoName.empty())
out.push_back(RhoName);
if (!PhiName.empty())
out.push_back(PhiName);
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TDistilVectors<FImpl>::setup(void)
{
// We expect the perambulator to have been created with these indices
auto &perambulator = envGet(PerambTensor, par().perambulator);
if (!perambulator.ValidateIndexNames())
{
HADRONS_ERROR(Range,"Perambulator index names bad");
}
const DistilParameters &dp{envGet(DistilParameters, par().DistilParams)};
const int Nt{env().getDim(Tdir)};
const bool full_tdil{ dp.TI == Nt };
const int Nt_inv{ full_tdil ? 1 : dp.TI };
if (!RhoName.empty())
envCreate(std::vector<FermionField>, RhoName, 1, dp.nnoise*dp.LI*dp.SI*Nt_inv, envGetGrid(FermionField));
if (!PhiName.empty())
envCreate(std::vector<FermionField>, PhiName, 1, dp.nnoise*dp.LI*dp.SI*Nt_inv, envGetGrid(FermionField));
Coordinate latt_size = GridDefaultLatt();
Coordinate mpi_layout = GridDefaultMpi();
Coordinate simd_layout_3 = GridDefaultSimd(Nd-1, vComplex::Nsimd());
latt_size[Nd-1] = 1;
simd_layout_3.push_back( 1 );
mpi_layout[Nd-1] = 1;
GridCartesian * const grid4d{env().getGrid()};
MakeLowerDimGrid(grid3d, grid4d);
envTmp(LatticeSpinColourVector, "source4d",1,LatticeSpinColourVector(grid4d));
envTmp(LatticeSpinColourVector, "source3d",1,LatticeSpinColourVector(grid3d.get()));
envTmp(LatticeColourVector, "source3d_nospin",1,LatticeColourVector(grid3d.get()));
envTmp(LatticeSpinColourVector, "sink3d",1,LatticeSpinColourVector(grid3d.get()));
envTmp(LatticeColourVector, "evec3d",1,LatticeColourVector(grid3d.get()));
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TDistilVectors<FImpl>::execute(void)
{
auto &noise = envGet(NoiseTensor, par().noise);
auto &perambulator = envGet(PerambTensor, par().perambulator);
auto &epack = envGet(Grid::Hadrons::EigenPack<LatticeColourVector>, par().lapevec);
const DistilParameters &dp{envGet(DistilParameters, par().DistilParams)};
envGetTmp(LatticeSpinColourVector, source4d);
envGetTmp(LatticeSpinColourVector, source3d);
envGetTmp(LatticeColourVector, source3d_nospin);
envGetTmp(LatticeSpinColourVector, sink3d);
envGetTmp(LatticeColourVector, evec3d);
GridCartesian * const grid4d{env().getGrid()};
const int Ntlocal{ grid4d->LocalDimensions()[3] };
const int Ntfirst{ grid4d->LocalStarts()[3] };
const int Nt{env().getDim(Tdir)};
const bool full_tdil{ dp.TI == Nt };
const int Nt_inv{ full_tdil ? 1 : dp.TI };
int vecindex;
if (!RhoName.empty())
{
auto &rho = envGet(std::vector<FermionField>, RhoName);
for (int inoise = 0; inoise < dp.nnoise; inoise++)
{
for (int dk = 0; dk < dp.LI; dk++)
{
for (int dt = 0; dt < Nt_inv; dt++)
{
for (int ds = 0; ds < dp.SI; ds++)
{
vecindex = inoise + dp.nnoise * (dk + dp.LI * (ds + dp.SI * dt));
rho[vecindex] = 0;
for (int it = dt; it < Nt; it += dp.TI)
{
const int t_inv{full_tdil ? dp.tsrc : it};
if (t_inv >= Ntfirst && t_inv < Ntfirst + Ntlocal)
{
for (int ik = dk; ik < dp.nvec; ik += dp.LI)
{
for (int is = ds; is < Ns; is += dp.SI)
{
ExtractSliceLocal(evec3d,epack.evec[ik],0,t_inv-Ntfirst,Tdir);
source3d_nospin = evec3d * noise.tensor(inoise, t_inv, ik, is);
source3d=0;
pokeSpin(source3d,source3d_nospin,is);
source4d=0;
InsertSliceLocal(source3d,source4d,0,t_inv-Ntfirst,Tdir);
rho[vecindex] += source4d;
}
}
}
}
}
}
}
}
}
if (!PhiName.empty())
{
auto &phi = envGet(std::vector<FermionField>, PhiName);
for (int inoise = 0; inoise < dp.nnoise; inoise++)
{
for (int dk = 0; dk < dp.LI; dk++)
{
for (int dt = 0; dt < Nt_inv; dt++)
{
for (int ds = 0; ds < dp.SI; ds++)
{
vecindex = inoise + dp.nnoise * (dk + dp.LI * (ds + dp.SI * dt));
phi[vecindex] = 0;
for (int t = Ntfirst; t < Ntfirst + Ntlocal; t++)
{
sink3d=0;
for (int ivec = 0; ivec < dp.nvec; ivec++)
{
ExtractSliceLocal(evec3d,epack.evec[ivec],0,t-Ntfirst,Tdir);
sink3d += evec3d * perambulator.tensor(t, ivec, dk, inoise,dt,ds);
}
InsertSliceLocal(sink3d,phi[vecindex],0,t-Ntfirst,Tdir);
}
}
}
}
}
}
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MDistil_DistilVectors_hpp_

36
Hadrons/Modules/MDistil/LapEvec.cc Normal file
View File

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

335
Hadrons/Modules/MDistil/LapEvec.hpp Normal file
View File

@ -0,0 +1,335 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MDistil/LapEvec.hpp
Copyright (C) 2019
Author: Felix Erben <ferben@ed.ac.uk>
Author: Michael Marshall <Michael.Marshall@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_MDistil_LapEvec_hpp_
#define Hadrons_MDistil_LapEvec_hpp_
#include <Hadrons/Modules/MDistil/Distil.hpp>
BEGIN_HADRONS_NAMESPACE
BEGIN_MODULE_NAMESPACE(MDistil)
/******************************************************************************
Laplacian eigenvectors - parameters
Computes the eigenvectors of the 3D-Laplacian, built from stout-smeared
gauge links with the specified number of steps and smearing parameter rho.
The smearing is only applied to the spatial components of the gauge field,
i.e. rho_{4i} = rho_{i4} = rho_{44} = 0.
Chebyshev-preconditioning is needed for convergence of the nvec lowest
eigenvectors.
******************************************************************************/
struct StoutParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(StoutParameters,
int, steps,
double, rho)
StoutParameters() = default;
template <class ReaderClass> StoutParameters(Reader<ReaderClass>& Reader){read(Reader,"StoutSmearing",*this);}
};
struct ChebyshevParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(ChebyshevParameters,
int, PolyOrder,
double, alpha,
double, beta)
ChebyshevParameters() = default;
template <class ReaderClass> ChebyshevParameters(Reader<ReaderClass>& Reader){read(Reader,"Chebyshev",*this);}
};
struct LanczosParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(LanczosParameters,
int, Nvec,
int, Nk,
int, Np,
int, MaxIt,
double, resid,
int, IRLLog)
LanczosParameters() = default;
template <class ReaderClass> LanczosParameters(Reader<ReaderClass>& Reader){read(Reader,"Lanczos",*this);}
};
// These are the actual parameters passed to the module during construction
struct LapEvecPar: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(LapEvecPar
,std::string, gauge
,StoutParameters, Stout
,ChebyshevParameters, Cheby
,LanczosParameters, Lanczos
,std::string, FileName)
};
/******************************************************************************
Laplacian eigenvectors - Module (class) definition
******************************************************************************/
template <typename GImpl>
class TLapEvec: public Module<LapEvecPar>
{
public:
GAUGE_TYPE_ALIASES(GImpl,);
// constructor
TLapEvec(const std::string name);
// destructor
virtual ~TLapEvec(void) {};
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
protected:
std::unique_ptr<GridCartesian> gridLD; // Owned by me, so I must delete it
};
MODULE_REGISTER_TMP(LapEvec, TLapEvec<GIMPL>, MDistil);
/******************************************************************************
TLapEvec implementation
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename GImpl>
TLapEvec<GImpl>::TLapEvec(const std::string name) : Module<LapEvecPar>(name) {}
// dependencies/products ///////////////////////////////////////////////////////
template <typename GImpl>
std::vector<std::string> TLapEvec<GImpl>::getInput(void)
{
return std::vector<std::string>{par().gauge};
}
template <typename GImpl>
std::vector<std::string> TLapEvec<GImpl>::getOutput(void)
{
return {getName()}; // This is the higher dimensional eigenpack
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename GImpl>
void TLapEvec<GImpl>::setup(void)
{
GridCartesian * gridHD = env().getGrid();
MakeLowerDimGrid(gridLD,gridHD);
const int Ntlocal{gridHD->LocalDimensions()[Tdir]};
// Temporaries
envTmpLat(GaugeField, "Umu_stout");
envTmpLat(GaugeField, "Umu_smear");
envTmp(LatticeGaugeField, "UmuNoTime",1,LatticeGaugeField(gridLD.get()));
envTmp(LatticeColourVector, "src",1,LatticeColourVector(gridLD.get()));
envTmp(std::vector<LapEvecs>, "eig",1,std::vector<LapEvecs>(Ntlocal));
// Output objects
envCreate(LapEvecs, getName(), 1, par().Lanczos.Nvec, gridHD);
}
/*************************************************************************************
-Grad^2 (Peardon, 2009, pg 2, equation 3, https://arxiv.org/abs/0905.2160)
Field Type of field the operator will be applied to
GaugeField Gauge field the operator will smear using
*************************************************************************************/
template<typename Field, typename GaugeField=LatticeGaugeField>
class Laplacian3D : public LinearOperatorBase<Field>, public LinearFunction<Field> {
typedef typename GaugeField::vector_type vCoeff_t;
public:
int nd; // number of spatial dimensions
std::vector<Lattice<iColourMatrix<vCoeff_t> > > U;
// Construct this operator given a gauge field and the number of dimensions it should act on
Laplacian3D( GaugeField& gf, int dimSpatial = Tdir ) : nd{dimSpatial}
{
if (dimSpatial<1)
{
HADRONS_ERROR(Range,"Must be at least one spatial dimension");
}
for (int mu = 0 ; mu < nd ; mu++)
U.push_back(PeekIndex<LorentzIndex>(gf,mu));
}
// Apply this operator to "in", return result in "out"
void operator()(const Field& in, Field& out) {
if (nd > in.Grid()->Nd())
{
HADRONS_ERROR(Range,"nd too large");
}
conformable( in, out );
out = ( ( Real ) ( 2 * nd ) ) * in;
Field tmp_(in.Grid());
typedef typename GaugeField::vector_type vCoeff_t;
for (int mu = 0 ; mu < nd ; mu++)
{
out -= U[mu] * Cshift( in, mu, 1);
tmp_ = adj( U[mu] ) * in;
out -= Cshift(tmp_,mu,-1);
}
}
void OpDiag (const Field &in, Field &out) { HADRONS_ERROR(Definition, "OpDiag() undefined"); };
void OpDir (const Field &in, Field &out,int dir,int disp) { HADRONS_ERROR(Definition, "OpDir() undefined"); };
void Op (const Field &in, Field &out) { HADRONS_ERROR(Definition, "Op() undefined"); };
void AdjOp (const Field &in, Field &out) { HADRONS_ERROR(Definition, "AdjOp() undefined"); };
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2) { HADRONS_ERROR(Definition, "HermOpAndNorm() undefined"); };
void HermOp(const Field &in, Field &out) { operator()(in,out); };
};
template<typename Field>
class Laplacian3DHerm : public LinearFunction<Field> {
public:
OperatorFunction<Field> & poly_;
LinearOperatorBase<Field> &Linop_;
Laplacian3DHerm(OperatorFunction<Field> & poly,LinearOperatorBase<Field>& linop)
: poly_{poly}, Linop_{linop} {}
void operator()(const Field& in, Field& out)
{
poly_(Linop_,in,out);
}
};
/******************************************************************************
Calculate low-mode eigenvalues of the Laplacian
******************************************************************************/
// execution ///////////////////////////////////////////////////////////////////
template <typename GImpl>
void TLapEvec<GImpl>::execute(void)
{
const ChebyshevParameters &ChebPar{par().Cheby};
const LanczosParameters &LPar{par().Lanczos};
// Disable IRL logging if requested
LOG(Message) << "IRLLog=" << LPar.IRLLog << std::endl;
const int PreviousIRLLogState{GridLogIRL.isActive()};
GridLogIRL.Active( LPar.IRLLog == 0 ? 0 : 1 );
// Stout smearing
envGetTmp(GaugeField, Umu_smear);
Umu_smear = envGet(GaugeField, par().gauge); // The smeared field starts off as the Gauge field
LOG(Message) << "Initial plaquette: " << WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu_smear) << std::endl;
const StoutParameters &Stout{par().Stout};
if( Stout.steps )
{
envGetTmp(GaugeField, Umu_stout);
Smear_Stout<PeriodicGimplR> LS(Stout.rho, Tdir); // spatial smearing only
for (int i = 0; i < Stout.steps; i++) {
LS.smear(Umu_stout, Umu_smear);
Umu_smear = Umu_stout;
}
LOG(Message) << "Smeared plaquette: " << WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu_smear) << std::endl;
}
////////////////////////////////////////////////////////////////////////
// Invert nabla operator separately on each time-slice
////////////////////////////////////////////////////////////////////////
auto & eig4d = envGet(LapEvecs, getName() );
envGetTmp(std::vector<LapEvecs>, eig); // Eigenpack for each timeslice
envGetTmp(LatticeGaugeField, UmuNoTime); // Gauge field without time dimension
envGetTmp(LatticeColourVector, src);
GridCartesian * gridHD = env().getGrid();
const int Ntlocal{gridHD->LocalDimensions()[Tdir]};
const int Ntfirst{gridHD->LocalStarts()[Tdir]};
uint32_t ConvergenceErrors{0};
for (int t = 0; t < Ntlocal; t++ )
{
LOG(Message) << "------------------------------------------------------------" << std::endl;
LOG(Message) << " Compute eigenpack, local timeslice = " << t << " / " << Ntlocal << std::endl;
LOG(Message) << "------------------------------------------------------------" << std::endl;
eig[t].resize(LPar.Nk+LPar.Np,gridLD.get());
// Construct smearing operator
ExtractSliceLocal(UmuNoTime,Umu_smear,0,t,Tdir); // switch to 3d/4d objects
Laplacian3D<LatticeColourVector> Nabla(UmuNoTime);
LOG(Message) << "Chebyshev preconditioning to order " << ChebPar.PolyOrder
<< " with parameters (alpha,beta) = (" << ChebPar.alpha << "," << ChebPar.beta << ")" << std::endl;
Chebyshev<LatticeColourVector> Cheb(ChebPar.alpha,ChebPar.beta,ChebPar.PolyOrder);
// Construct source vector according to Test_dwf_compressed_lanczos.cc
src = 11.0; // NB: This is a dummy parameter and just needs to be non-zero
RealD nn = norm2(src);
nn = Grid::sqrt(nn);
src = src * (1.0/nn);
Laplacian3DHerm<LatticeColourVector> NablaCheby(Cheb,Nabla);
ImplicitlyRestartedLanczos<LatticeColourVector>
IRL(NablaCheby,Nabla,LPar.Nvec,LPar.Nk,LPar.Nk+LPar.Np,LPar.resid,LPar.MaxIt);
int Nconv = 0;
IRL.calc(eig[t].eval,eig[t].evec,src,Nconv);
if (Nconv < LPar.Nvec)
{
// NB: Can't assert here since we are processing local slices - i.e. not all nodes would assert
ConvergenceErrors = 1;
LOG(Error) << "MDistil::LapEvec : Not enough eigenvectors converged. If this occurs in practice, we should modify the eigensolver to iterate once more to ensure the second convergence test does not take us below the requested number of eigenvectors" << std::endl;
}
if( Nconv != LPar.Nvec )
eig[t].resize(LPar.Nvec, gridLD.get());
RotateEigen( eig[t].evec ); // Rotate the eigenvectors into our phase convention
for (int i=0;i<LPar.Nvec;i++){
InsertSliceLocal(eig[t].evec[i],eig4d.evec[i],0,t,Tdir);
if(t==0 && Ntfirst==0)
eig4d.eval[i] = eig[t].eval[i]; // TODO: Discuss: is this needed? Is there a better way?
}
}
GridLogIRL.Active( PreviousIRLLogState );
gridHD->GlobalSum(ConvergenceErrors);
if(ConvergenceErrors!=0)
{
HADRONS_ERROR(Program,"The eingensolver failed to find enough eigenvectors on at least one node");
}
// Now write out the 4d eigenvectors
std::string sEigenPackName(par().FileName);
if( !sEigenPackName.empty() )
{
eig4d.record.solverXml = parString();
ModuleBase * b{vm().getModule(par().gauge)};
std::string sOperatorXml{ "<module><id><type>" };
sOperatorXml.append( b->getRegisteredName() );
sOperatorXml.append( "</type></id><options>" );
sOperatorXml.append( b->parString() );
sOperatorXml.append( "</options></module>" );
eig4d.record.operatorXml = sOperatorXml;
sEigenPackName.append(".");
sEigenPackName.append(std::to_string(vm().getTrajectory()));
eig4d.write(sEigenPackName,false);
}
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MDistil_LapEvec_hpp_

36
Hadrons/Modules/MDistil/Noises.cc Normal file
View File

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

146
Hadrons/Modules/MDistil/Noises.hpp Normal file
View File

@ -0,0 +1,146 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MDistil/Noises.hpp
Copyright (C) 2019
Author: Felix Erben <ferben@ed.ac.uk>
Author: Michael Marshall <Michael.Marshall@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_MDistil_Noises_hpp_
#define Hadrons_MDistil_Noises_hpp_
#include <Hadrons/Modules/MDistil/Distil.hpp>
BEGIN_HADRONS_NAMESPACE
BEGIN_MODULE_NAMESPACE(MDistil)
/******************************************************************************
* Noises *
******************************************************************************/
class NoisesPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(NoisesPar,
std::string, DistilParams,
std::string, NoiseFileName)
};
template <typename FImpl>
class TNoises: public Module<NoisesPar>
{
public:
// constructor
TNoises(const std::string name);
// destructor
virtual ~TNoises(void) {};
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
};
MODULE_REGISTER_TMP(Noises, TNoises<FIMPL>, MDistil);
/******************************************************************************
* TNoises implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TNoises<FImpl>::TNoises(const std::string name) : Module<NoisesPar>(name) {}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TNoises<FImpl>::getInput(void)
{
return {par().DistilParams};
}
template <typename FImpl>
std::vector<std::string> TNoises<FImpl>::getOutput(void)
{
return {getName()};
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TNoises<FImpl>::setup(void)
{
const DistilParameters &dp{envGet(DistilParameters, par().DistilParams)};
const int Nt{env().getDim(Tdir)};
std::cout << dp.nnoise << dp.nvec << Nt << Ns << std::endl;
envCreate(NoiseTensor, getName(), 1, dp.nnoise, Nt, dp.nvec, Ns);
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TNoises<FImpl>::execute(void)
{
const DistilParameters &dp{envGet(DistilParameters, par().DistilParams)};
const int Nt{env().getDim(Tdir)};
const bool full_tdil{ dp.TI == Nt };
const bool exact_distillation{ full_tdil && dp.LI == dp.nvec };
// We use our own seeds so we can specify different noises per quark
Real rn;
auto &noise = envGet(NoiseTensor, getName());
for (int inoise = 0; inoise < dp.nnoise; inoise++)
{
for (int t = 0; t < Nt; t++)
{
for (int ivec = 0; ivec < dp.nvec; ivec++)
{
for (int is = 0; is < Ns; is++)
{
if (exact_distillation)
{
noise.tensor(inoise, t, ivec, is) = 1.;
}
else
{
random(rngSerial(),rn);
// We could use a greater number of complex roots of unity
// ... but this seems to work well
noise.tensor(inoise, t, ivec, is) = (rn > 0.5) ? -1 : 1;
}
}
}
}
}
if (env().getGrid()->IsBoss())
{
std::string sName {par().NoiseFileName};
sName.append(".");
sName.append(std::to_string(vm().getTrajectory()));
noise.write(sName.c_str());
}
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif

36
Hadrons/Modules/MDistil/PerambFromSolve.cc Normal file
View File

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

183
Hadrons/Modules/MDistil/PerambFromSolve.hpp Normal file
View File

@ -0,0 +1,183 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MDistil/PerambFromSolve.hpp
Copyright (C) 2019
Author: Felix Erben <ferben@ed.ac.uk>
Author: Michael Marshall <Michael.Marshall@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_MDistil_PerambFromSolve_hpp_
#define Hadrons_MDistil_PerambFromSolve_hpp_
#include <Hadrons/Modules/MDistil/Distil.hpp>
BEGIN_HADRONS_NAMESPACE
BEGIN_MODULE_NAMESPACE(MDistil)
/******************************************************************************
* PerambFromSolve
This module computes a perambulator from an already completed solve.
Optionally, the number of eigenvectors used in the perambulator and the
parameter LI can be chosen to be lower than the ones in the solve, allowing
for a study of the signal with different values of nvec.
LI_reduced : value of LI actually used in the computation
nvec_reduced: value of nvec actually used in the computation
LI : value of LI used to compute the 'solve'
nvec : value of nvec used to compute the 'solve'
******************************************************************************/
class PerambFromSolvePar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(PerambFromSolvePar,
std::string, lapevec,
std::string, PerambFileName,
std::string, solve,
int, nvec_reduced,
int, LI_reduced,
std::string, DistilParams);
};
template <typename FImpl>
class TPerambFromSolve: public Module<PerambFromSolvePar>
{
public:
FERM_TYPE_ALIASES(FImpl,);
// constructor
TPerambFromSolve(const std::string name);
// destructor
virtual ~TPerambFromSolve(void) {};
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
protected:
std::unique_ptr<GridCartesian> grid3d; // Owned by me, so I must delete it
};
MODULE_REGISTER_TMP(PerambFromSolve, TPerambFromSolve<FIMPL>, MDistil);
/******************************************************************************
* TPerambFromSolve implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TPerambFromSolve<FImpl>::TPerambFromSolve(const std::string name) : Module<PerambFromSolvePar>(name){}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TPerambFromSolve<FImpl>::getInput(void)
{
return {par().solve, par().lapevec, par().DistilParams};
}
template <typename FImpl>
std::vector<std::string> TPerambFromSolve<FImpl>::getOutput(void)
{
return std::vector<std::string>{getName()};
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TPerambFromSolve<FImpl>::setup(void)
{
const DistilParameters & dp{envGet(MDistil::DistilParameters, par().DistilParams)};
const int Nt{env().getDim(Tdir)};
const bool full_tdil{ dp.TI == Nt };
const int Nt_inv{ full_tdil ? 1 : dp.TI };
MakeLowerDimGrid( grid3d, env().getGrid() );
const int nvec_reduced{par().nvec_reduced};
const int LI_reduced{ par().LI_reduced};
envCreate(PerambTensor, getName(), 1, Nt,nvec_reduced,LI_reduced,dp.nnoise,Nt_inv,dp.SI);
envCreate(NoiseTensor, getName() + "_noise", 1, dp.nnoise, Nt, dp.nvec, Ns );
envTmp(LatticeColourVector, "result3d_nospin",1,LatticeColourVector(grid3d.get()));
envTmp(LatticeColourVector, "evec3d",1,LatticeColourVector(grid3d.get()));
envTmpLat(LatticeColourVector, "result4d_nospin");
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TPerambFromSolve<FImpl>::execute(void)
{
GridCartesian * grid4d = env().getGrid();
const int Ntlocal{grid4d->LocalDimensions()[3]};
const int Ntfirst{grid4d->LocalStarts()[3]};
const DistilParameters &dp{envGet(DistilParameters, par().DistilParams)};
const int Nt{env().getDim(Tdir)};
const bool full_tdil{ dp.TI == Nt };
const int Nt_inv{ full_tdil ? 1 : dp.TI };
const int nvec_reduced{par().nvec_reduced};
const int LI_reduced{ par().LI_reduced};
auto &perambulator = envGet(PerambTensor, getName());
auto &solve = envGet(std::vector<FermionField>, par().solve);
auto &epack = envGet(Grid::Hadrons::EigenPack<LatticeColourVector>, par().lapevec);
envGetTmp(LatticeColourVector, result4d_nospin);
envGetTmp(LatticeColourVector, result3d_nospin);
envGetTmp(LatticeColourVector, evec3d);
for (int inoise = 0; inoise < dp.nnoise; inoise++)
{
for (int dk = 0; dk < LI_reduced; dk++)
{
for (int dt = 0; dt < Nt_inv; dt++)
{
for (int ds = 0; ds < dp.SI; ds++)
{
for (int is = 0; is < Ns; is++)
{
result4d_nospin = peekSpin(solve[inoise+dp.nnoise*(dk+dp.LI*(dt+Nt_inv*ds))],is);
for (int t = Ntfirst; t < Ntfirst + Ntlocal; t++)
{
ExtractSliceLocal(result3d_nospin,result4d_nospin,0,t-Ntfirst,Tdir);
for (int ivec = 0; ivec < nvec_reduced; ivec++)
{
ExtractSliceLocal(evec3d,epack.evec[ivec],0,t-Ntfirst,Tdir);
pokeSpin(perambulator.tensor(t, ivec, dk, inoise,dt,ds),static_cast<Complex>(innerProduct(evec3d, result3d_nospin)),is);
LOG(Message) << "perambulator(t, ivec, dk, inoise,dt,ds)(is) = (" << t << "," << ivec << "," << dk << "," << inoise << "," << dt << "," << ds << ")(" << is << ") = " << perambulator.tensor(t, ivec, dk, inoise,dt,ds)()(is)() << std::endl;
}
}
}
}
}
}
}
if(grid4d->IsBoss())
{
std::string sPerambName{par().PerambFileName};
sPerambName.append( "." );
sPerambName.append( std::to_string(vm().getTrajectory()));
perambulator.write(sPerambName.c_str());
}
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MDistil_PerambFromSolve_hpp_

57
Hadrons/Modules/MDistil/Perambulator.cc Normal file
View File

@ -0,0 +1,57 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MDistil/Perambulator.cc
Copyright (C) 2019
Author: Felix Erben <ferben@ed.ac.uk>
Author: Michael Marshall <Michael.Marshall@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Hadrons/Modules/MDistil/Perambulator.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MDistil;
template class Grid::Hadrons::MDistil::TPerambulator<FIMPL>;
BEGIN_HADRONS_NAMESPACE
// Global constants for distillation
#ifdef HAVE_HDF5
extern const std::string NamedTensorFileExtension{".h5"};
#else
extern const std::string NamedTensorFileExtension{".dat"};
#endif
BEGIN_MODULE_NAMESPACE(MDistil)
const std::string NoiseTensor::Name__{"Noises"};
const std::array<std::string, 4> NoiseTensor::DefaultIndexNames__{"nNoise", "nT", "nVec", "nS"};
const std::string PerambTensor::Name__{"Perambulator"};
const std::array<std::string, 6> PerambTensor::DefaultIndexNames__{"nT", "nVec", "LI", "nNoise", "nT_inv", "SI"};
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE

263
Hadrons/Modules/MDistil/Perambulator.hpp Normal file
View File

@ -0,0 +1,263 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MDistil/Perambulator.hpp
Copyright (C) 2019
Author: Felix Erben <ferben@ed.ac.uk>
Author: Michael Marshall <Michael.Marshall@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_MDistil_Perambulator_hpp_
#define Hadrons_MDistil_Perambulator_hpp_
#include <Hadrons/Modules/MDistil/Distil.hpp>
BEGIN_HADRONS_NAMESPACE
BEGIN_MODULE_NAMESPACE(MDistil)
/******************************************************************************
* Perambulator *
******************************************************************************/
class PerambulatorPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(PerambulatorPar,
std::string, lapevec,
std::string, solver,
std::string, noise,
std::string, PerambFileName,
std::string, UnsmearedSinkFileName,
std::string, DistilParams);
};
template <typename FImpl>
class TPerambulator: public Module<PerambulatorPar>
{
public:
FERM_TYPE_ALIASES(FImpl,);
SOLVER_TYPE_ALIASES(FImpl,);
// constructor
TPerambulator(const std::string name);
// destructor
virtual ~TPerambulator(void) {};
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
protected:
std::unique_ptr<GridCartesian> grid3d; // Owned by me, so I must delete it
unsigned int Ls_;
};
MODULE_REGISTER_TMP(Perambulator, TPerambulator<FIMPL>, MDistil);
/******************************************************************************
* TPerambulator implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TPerambulator<FImpl>::TPerambulator(const std::string name) : Module<PerambulatorPar>(name) {}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TPerambulator<FImpl>::getInput(void)
{
return {par().lapevec, par().solver, par().noise, par().DistilParams};
}
template <typename FImpl>
std::vector<std::string> TPerambulator<FImpl>::getOutput(void)
{
return {getName(), getName() + "_unsmeared_sink"};
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TPerambulator<FImpl>::setup(void)
{
MakeLowerDimGrid(grid3d, env().getGrid());
const DistilParameters &dp = envGet(DistilParameters, par().DistilParams);
const int Nt{env().getDim(Tdir)};
const bool full_tdil{ dp.TI == Nt };
const int Nt_inv{ full_tdil ? 1 : dp.TI };
envCreate(PerambTensor, getName(), 1, Nt, dp.nvec, dp.LI, dp.nnoise, Nt_inv, dp.SI);
envCreate(std::vector<FermionField>, getName() + "_unsmeared_sink", 1,
dp.nnoise*dp.LI*Ns*Nt_inv, envGetGrid(FermionField));
envTmpLat(LatticeSpinColourVector, "dist_source");
envTmpLat(LatticeSpinColourVector, "source4d");
envTmp(LatticeSpinColourVector, "source3d",1,LatticeSpinColourVector(grid3d.get()));
envTmp(LatticeColourVector, "source3d_nospin",1,LatticeColourVector(grid3d.get()));
envTmpLat(LatticeSpinColourVector, "result4d");
envTmpLat(LatticeColourVector, "result4d_nospin");
envTmp(LatticeColourVector, "result3d_nospin",1,LatticeColourVector(grid3d.get()));
envTmp(LatticeColourVector, "evec3d",1,LatticeColourVector(grid3d.get()));
Ls_ = env().getObjectLs(par().solver);
envTmpLat(FermionField, "v4dtmp");
envTmpLat(FermionField, "v5dtmp", Ls_);
envTmpLat(FermionField, "v5dtmp_sol", Ls_);
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TPerambulator<FImpl>::execute(void)
{
const DistilParameters &dp{ envGet(DistilParameters, par().DistilParams) };
const int Nt{env().getDim(Tdir)};
const bool full_tdil{ dp.TI == Nt };
const int Nt_inv{ full_tdil ? 1 : dp.TI };
auto &solver=envGet(Solver, par().solver);
auto &mat = solver.getFMat();
envGetTmp(FermionField, v4dtmp);
envGetTmp(FermionField, v5dtmp);
envGetTmp(FermionField, v5dtmp_sol);
auto &noise = envGet(NoiseTensor, par().noise);
auto &perambulator = envGet(PerambTensor, getName());
auto &epack = envGet(LapEvecs, par().lapevec);
auto &unsmeared_sink = envGet(std::vector<FermionField>, getName() + "_unsmeared_sink");
envGetTmp(LatticeSpinColourVector, dist_source);
envGetTmp(LatticeSpinColourVector, source4d);
envGetTmp(LatticeSpinColourVector, source3d);
envGetTmp(LatticeColourVector, source3d_nospin);
envGetTmp(LatticeSpinColourVector, result4d);
envGetTmp(LatticeColourVector, result4d_nospin);
envGetTmp(LatticeColourVector, result3d_nospin);
envGetTmp(LatticeColourVector, evec3d);
GridCartesian * const grid4d{ env().getGrid() }; // Owned by environment (so I won't delete it)
const int Ntlocal{grid4d->LocalDimensions()[3]};
const int Ntfirst{grid4d->LocalStarts()[3]};
const std::string UnsmearedSinkFileName{ par().UnsmearedSinkFileName };
for (int inoise = 0; inoise < dp.nnoise; inoise++)
{
for (int dk = 0; dk < dp.LI; dk++)
{
for (int dt = 0; dt < Nt_inv; dt++)
{
for (int ds = 0; ds < dp.SI; ds++)
{
LOG(Message) << "LapH source vector from noise " << inoise << " and dilution component (d_k,d_t,d_alpha) : (" << dk << ","<< dt << "," << ds << ")" << std::endl;
dist_source = 0;
evec3d = 0;
for (int it = dt; it < Nt; it += dp.TI)
{
const int t_inv{full_tdil ? dp.tsrc : it};
if( t_inv >= Ntfirst && t_inv < Ntfirst + Ntlocal )
{
for (int ik = dk; ik < dp.nvec; ik += dp.LI)
{
for (int is = ds; is < Ns; is += dp.SI)
{
ExtractSliceLocal(evec3d,epack.evec[ik],0,t_inv-Ntfirst,Tdir);
source3d_nospin = evec3d * noise.tensor(inoise, t_inv, ik, is);
source3d=0;
pokeSpin(source3d,source3d_nospin,is);
source4d=0;
InsertSliceLocal(source3d,source4d,0,t_inv-Ntfirst,Tdir);
dist_source += source4d;
}
}
}
}
result4d=0;
v4dtmp = dist_source;
if (Ls_ == 1)
solver(result4d, v4dtmp);
else
{
mat.ImportPhysicalFermionSource(v4dtmp, v5dtmp);
solver(v5dtmp_sol, v5dtmp);
mat.ExportPhysicalFermionSolution(v5dtmp_sol, v4dtmp);
result4d = v4dtmp;
}
if (!UnsmearedSinkFileName.empty())
unsmeared_sink[inoise+dp.nnoise*(dk+dp.LI*(dt+Nt_inv*ds))] = result4d;
for (int is = 0; is < Ns; is++)
{
result4d_nospin = peekSpin(result4d,is);
for (int t = Ntfirst; t < Ntfirst + Ntlocal; t++)
{
ExtractSliceLocal(result3d_nospin,result4d_nospin,0,t-Ntfirst,Tdir);
for (int ivec = 0; ivec < dp.nvec; ivec++)
{
ExtractSliceLocal(evec3d,epack.evec[ivec],0,t-Ntfirst,Tdir);
pokeSpin(perambulator.tensor(t, ivec, dk, inoise,dt,ds),static_cast<Complex>(innerProduct(evec3d, result3d_nospin)),is);
}
}
}
}
}
}
}
// Now share my timeslice data with other members of the grid
const int NumSlices{grid4d->_processors[Tdir] / grid3d->_processors[Tdir]};
if (NumSlices > 1)
{
LOG(Debug) << "Sharing perambulator data with other nodes" << std::endl;
const int MySlice {grid4d->_processor_coor[Tdir]};
const int SliceCount {static_cast<int>(perambulator.tensor.size()/NumSlices)};
PerambTensor::Scalar * const MyData {perambulator.tensor.data()+MySlice*SliceCount};
Coordinate coor(Nd);
for (int i = 0 ; i < Tdir ; i++) coor[i] = grid4d->_processor_coor[i];
std::vector<CommsRequest_t> reqs(0);
for (int i = 1; i < NumSlices ; i++)
{
coor[Tdir] = (MySlice+i)%NumSlices;
const int SendRank { grid4d->RankFromProcessorCoor(coor) };
const int RecvSlice { ( MySlice - i + NumSlices ) % NumSlices };
coor[Tdir] = RecvSlice;
const auto RecvRank = grid4d->RankFromProcessorCoor(coor);
grid4d->SendToRecvFromBegin(reqs,MyData,SendRank, perambulator.tensor.data()
+ RecvSlice*SliceCount,RecvRank,SliceCount*sizeof(PerambTensor::Scalar));
}
grid4d->SendToRecvFromComplete(reqs);
}
// Save the perambulator to disk from the boss node
if (grid4d->IsBoss())
{
std::string sPerambName {par().PerambFileName};
sPerambName.append(".");
sPerambName.append(std::to_string(vm().getTrajectory()));
perambulator.write(sPerambName.c_str());
}
//Save the unsmeared sinks if filename specified
if (!UnsmearedSinkFileName.empty())
{
LOG(Message) << "Writing unsmeared sink to " << UnsmearedSinkFileName << std::endl;
A2AVectorsIo::write(UnsmearedSinkFileName, unsmeared_sink, false, vm().getTrajectory());
}
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MDistil_Perambulator_hpp_

7
Hadrons/Modules/MIO/LoadDistilNoise.cc Normal file
View File

@ -0,0 +1,7 @@
#include <Hadrons/Modules/MIO/LoadDistilNoise.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MIO;
template class Grid::Hadrons::MIO::TLoadDistilNoise<FIMPL>;

111
Hadrons/Modules/MIO/LoadDistilNoise.hpp Normal file
View File

@ -0,0 +1,111 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MDistil/LoadDistilNoise.hpp
Copyright (C) 2019
Author: Felix Erben <ferben@ed.ac.uk>
Author: Michael Marshall <Michael.Marshall@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_MIO_LoadDistilNoise_hpp_
#define Hadrons_MIO_LoadDistilNoise_hpp_
#include <Hadrons/Modules/MDistil/Distil.hpp>
BEGIN_HADRONS_NAMESPACE
BEGIN_MODULE_NAMESPACE(MIO)
/******************************************************************************
* LoadDistilNoise *
******************************************************************************/
class LoadDistilNoisePar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(LoadDistilNoisePar,
std::string, NoiseFileName,
std::string, DistilParams);
};
template <typename FImpl>
class TLoadDistilNoise: public Module<LoadDistilNoisePar>
{
public:
// constructor
TLoadDistilNoise(const std::string name);
// destructor
virtual ~TLoadDistilNoise(void) {};
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
};
MODULE_REGISTER_TMP(LoadDistilNoise, TLoadDistilNoise<FIMPL>, MIO);
/******************************************************************************
* TLoadDistilNoise implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TLoadDistilNoise<FImpl>::TLoadDistilNoise(const std::string name) : Module<LoadDistilNoisePar>(name) {}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TLoadDistilNoise<FImpl>::getInput(void)
{
return {par().DistilParams};
}
template <typename FImpl>
std::vector<std::string> TLoadDistilNoise<FImpl>::getOutput(void)
{
return {getName()};
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TLoadDistilNoise<FImpl>::setup(void)
{
const MDistil::DistilParameters &dp{envGet(MDistil::DistilParameters, par().DistilParams)};
const int Nt{env().getDim(Tdir)};
envCreate(MDistil::NoiseTensor, getName(), 1, dp.nnoise, Nt, dp.nvec, Ns);
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TLoadDistilNoise<FImpl>::execute(void)
{
auto &noises = envGet(MDistil::NoiseTensor, getName());
std::string sNoiseName{ par().NoiseFileName };
sNoiseName.append( 1, '.' );
sNoiseName.append( std::to_string( vm().getTrajectory() ) );
noises.read(sNoiseName.c_str());
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif

36
Hadrons/Modules/MIO/LoadPerambulator.cc Normal file
View File

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

113
Hadrons/Modules/MIO/LoadPerambulator.hpp Normal file
View File

@ -0,0 +1,113 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MDistil/LoadPerambulator.hpp
Copyright (C) 2019
Author: Felix Erben <ferben@ed.ac.uk>
Author: Michael Marshall <Michael.Marshall@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_MIO_LoadPerambulator_hpp_
#define Hadrons_MIO_LoadPerambulator_hpp_
#include <Hadrons/Modules/MDistil/Distil.hpp>
BEGIN_HADRONS_NAMESPACE
BEGIN_MODULE_NAMESPACE(MIO)
/******************************************************************************
* LoadPerambulator *
******************************************************************************/
class LoadPerambulatorPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(LoadPerambulatorPar,
std::string, PerambFileName,
std::string, DistilParams);
};
template <typename FImpl>
class TLoadPerambulator: public Module<LoadPerambulatorPar>
{
public:
// constructor
TLoadPerambulator(const std::string name);
// destructor
virtual ~TLoadPerambulator(void) {};
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
};
MODULE_REGISTER_TMP(LoadPerambulator, TLoadPerambulator<FIMPL>, MIO);
/******************************************************************************
* TLoadPerambulator implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TLoadPerambulator<FImpl>::TLoadPerambulator(const std::string name) : Module<LoadPerambulatorPar>(name) {}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TLoadPerambulator<FImpl>::getInput(void)
{
return {par().DistilParams};
}
template <typename FImpl>
std::vector<std::string> TLoadPerambulator<FImpl>::getOutput(void)
{
return {getName()};
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TLoadPerambulator<FImpl>::setup(void)
{
const MDistil::DistilParameters &dp{envGet(MDistil::DistilParameters, par().DistilParams)};
const int Nt{env().getDim(Tdir)};
const bool full_tdil{ dp.TI == Nt };
const int Nt_inv{ full_tdil ? 1 : dp.TI };
envCreate(MDistil::PerambTensor, getName(), 1, Nt,dp.nvec,dp.LI,dp.nnoise,Nt_inv,dp.SI);
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TLoadPerambulator<FImpl>::execute(void)
{
auto &perambulator = envGet(MDistil::PerambTensor, getName());
std::string sPerambName{ par().PerambFileName };
sPerambName.append( 1, '.' );
sPerambName.append( std::to_string( vm().getTrajectory() ) );
perambulator.read(sPerambName.c_str());
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif

2
Hadrons/Modules/MNPR/Amputate.hpp
View File

@ -187,7 +187,7 @@ void TAmputate<FImpl1, FImpl2>::execute(void)
result.Vamp.resize(Gamma::nGamma/2);
for( int mu=0; mu < Gamma::nGamma/2; mu++){
Gamma::Algebra gam = mu;
result.Vamp[mu] = 1/12.0*trace(adj(Gamma(mu*2+1))*g5*Sout_inv*g5*vertex[mu]*Sin_inv);
result.Vamp[mu] = static_cast<Real>( 1 / 12.0 ) * trace(adj(Gamma(mu*2+1))*g5*Sout_inv*g5*vertex[mu]*Sin_inv)()()();
LOG(Message) << "Vamp[" << mu << "] - " << result.Vamp[mu] << std::endl;
}
}

190
Hadrons/NamedTensor.hpp Normal file
View File

@ -0,0 +1,190 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/NamedTensor.hpp
Copyright (C) 2015-2019
Author: Felix Erben <ferben@ed.ac.uk>
Author: Michael Marshall <Michael.Marshall@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_NamedTensor_hpp_
#define Hadrons_NamedTensor_hpp_
#include <Hadrons/Global.hpp>
#include <Hadrons/EigenPack.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
NamedTensor contains:
1) Name of the tensor. By default, this is the tag name used for save / load
2) Eigen::Tensor of type Scalar_ and rank NumIndices_ (row-major order)
3) Name for each index
They can be persisted to / restored from disk. During restore, these validations are performed:
1) Tensor dimensionality must match
2) IndexNames are validated against current values
3) If the tensor has non-zero size, the tensor being loaded must have same extent in each dimension
******************************************************************************/
extern const std::string NamedTensorFileExtension;
template<typename Scalar_, int NumIndices_>
class NamedTensor : Serializable
{
public:
using Scalar = Scalar_;
static constexpr int NumIndices = NumIndices_;
using ET = Eigen::Tensor<Scalar_, NumIndices_, Eigen::RowMajor>;
using Index = typename ET::Index;
GRID_SERIALIZABLE_CLASS_MEMBERS(NamedTensor,
ET, tensor,
std::vector<std::string>, IndexNames );
// Name of the object and Index names as set in the constructor
const std::string &Name_;
const std::array<std::string, NumIndices_> &DefaultIndexNames_;
// Default constructor (assumes tensor will be loaded from file)
NamedTensor(const std::string &Name,
const std::array<std::string, NumIndices_> &indexNames)
: IndexNames{indexNames.begin(), indexNames.end()}, Name_{Name}, DefaultIndexNames_{indexNames} {}
// Construct a named tensor explicitly specifying size of each dimension
template<typename... IndexTypes>
NamedTensor(const std::string &Name,
const std::array<std::string, NumIndices_> &indexNames,
Eigen::Index firstDimension, IndexTypes... otherDimensions)
: tensor(firstDimension, otherDimensions...),
IndexNames{indexNames.begin(), indexNames.end()}, Name_{Name}, DefaultIndexNames_{indexNames}
{
if(sizeof...(otherDimensions) + 1 != NumIndices_)
{
HADRONS_ERROR(Argument, "NamedTensor: dimensions != tensor rank");
}
}
// Do my index names match the default for my type?
template<typename array_or_vector_of_string>
bool ValidateIndexNames( const array_or_vector_of_string &CheckNames ) const
{
return IndexNames.size() == CheckNames.size() && std::equal( IndexNames.begin(), IndexNames.end(), CheckNames.begin(),
[](const std::string &s1, const std::string &s2)
{
return s1.size() == s2.size() && std::equal( s1.begin(), s1.end(), s2.begin(),
[](const char & c1, const char & c2)
{ return c1 == c2 || std::toupper(c1) == std::toupper(c2); }); // case insensitive
});
}
bool ValidateIndexNames() const { return ValidateIndexNames(DefaultIndexNames_); }
#ifdef HAVE_HDF5
using Default_Reader = Grid::Hdf5Reader;
using Default_Writer = Grid::Hdf5Writer;
#else
using Default_Reader = Grid::BinaryReader;
using Default_Writer = Grid::BinaryWriter;
#endif
void write(const std::string &FileName, const std::string &Tag) const
{
std::string FileName_{FileName};
FileName_.append( NamedTensorFileExtension );
LOG(Message) << "Writing " << Name_ << " to file " << FileName_ << " tag " << Tag << std::endl;
Default_Writer w( FileName_ );
write( w, Tag, *this );
}
void write(const std::string &FileName) const { return write(FileName, Name_); }
// Read tensor.
// Validate:
// 1) index names (if requested)
// 2) index dimensions (if they are non-zero when called)
template<typename Reader> void read(Reader &r, bool bValidate, const std::string &Tag)
{
// Grab index names and dimensions
std::vector<std::string> OldIndexNames{std::move(IndexNames)};
const typename ET::Dimensions OldDimensions{tensor.dimensions()};
read(r, Tag, *this);
const typename ET::Dimensions & NewDimensions{tensor.dimensions()};
for (int i = 0; i < NumIndices_; i++)
if(OldDimensions[i] && OldDimensions[i] != NewDimensions[i])
{
HADRONS_ERROR(Size,"NamedTensor::read dimension size");
}
if (bValidate && !ValidateIndexNames(OldIndexNames))
{
HADRONS_ERROR(Definition,"NamedTensor::read dimension name");
}
}
template<typename Reader> void read(Reader &r, bool bValidate = true) { read(r, bValidate, Name_); }
inline void read (const std::string &FileName, bool bValidate, const std::string &Tag)
{
Default_Reader r(FileName + NamedTensorFileExtension);
read(r, bValidate, Tag);
}
inline void read (const std::string &FileName, bool bValidate= true) { return read(FileName, bValidate, Name_); }
};
/******************************************************************************
Common elements for distillation
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MDistil)
//Eigenvectors of the Laplacian
using LapEvecs = Grid::Hadrons::EigenPack<LatticeColourVector>;
// Noise vector (index order: nnoise, nt, nvec, ns)
class NoiseTensor : public NamedTensor<Complex, 4>
{
static const std::string Name__;
static const std::array<std::string, 4> DefaultIndexNames__;
public:
// Default constructor (assumes tensor will be loaded from file)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NoiseTensor() : NamedTensor{Name__, DefaultIndexNames__} {}
// Construct a named tensor explicitly specifying size of each dimension
template<typename... IndexTypes>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NoiseTensor(Eigen::Index nNoise, Eigen::Index nT, Eigen::Index nVec, Eigen::Index nS)
: NamedTensor{Name__, DefaultIndexNames__, nNoise, nT, nVec, nS} {}
};
class PerambTensor : public NamedTensor<SpinVector, 6>
{
static const std::string Name__;
static const std::array<std::string, 6> DefaultIndexNames__;
public:
// Default constructor (assumes tensor will be loaded from file)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE PerambTensor() : NamedTensor{Name__, DefaultIndexNames__} {}
// Construct a named tensor explicitly specifying size of each dimension
template<typename... IndexTypes>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE PerambTensor(Eigen::Index nT, Eigen::Index nVec, Eigen::Index LI, Eigen::Index nNoise, Eigen::Index nT_inv, Eigen::Index SI)
: NamedTensor{Name__, DefaultIndexNames__, nT, nVec, LI, nNoise, nT_inv, SI} {}
};
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_NamedTensor_hpp_

17
Hadrons/modules.inc
View File

@ -48,7 +48,9 @@ modules_cc =\
Modules/MNoise/FullVolumeSpinColorDiagonal.cc \
Modules/MIO/LoadA2AVectors.cc \
Modules/MIO/LoadBinary.cc \
Modules/MIO/LoadPerambulator.cc \
Modules/MIO/LoadA2AMatrixDiskVector.cc \
Modules/MIO/LoadDistilNoise.cc \
Modules/MIO/LoadNersc.cc \
Modules/MIO/LoadCoarseEigenPack.cc \
Modules/MIO/LoadCosmHol.cc \
@ -63,6 +65,12 @@ modules_cc =\
Modules/MScalarSUN/TrKinetic.cc \
Modules/MScalarSUN/StochFreeField.cc \
Modules/MScalarSUN/TrMag.cc \
Modules/MDistil/PerambFromSolve.cc \
Modules/MDistil/DistilVectors.cc \
Modules/MDistil/DistilPar.cc \
Modules/MDistil/LapEvec.cc \
Modules/MDistil/Perambulator.cc \
Modules/MDistil/Noises.cc \
Modules/MAction/MobiusDWF.cc \
Modules/MAction/DWF.cc \
Modules/MAction/ScaledDWF.cc \
@ -126,8 +134,10 @@ modules_hpp =\
Modules/MIO/LoadBinary.hpp \
Modules/MIO/LoadNersc.hpp \
Modules/MIO/LoadEigenPack.hpp \
Modules/MIO/LoadDistilNoise.hpp \
Modules/MIO/LoadA2AMatrixDiskVector.hpp \
Modules/MIO/LoadCoarseEigenPack.hpp \
Modules/MIO/LoadPerambulator.hpp \
Modules/MIO/LoadA2AVectors.hpp \
Modules/MIO/LoadCosmHol.hpp \
Modules/MScalarSUN/TwoPointNPR.hpp \
@ -141,6 +151,13 @@ modules_hpp =\
Modules/MScalarSUN/StochFreeField.hpp \
Modules/MScalarSUN/Div.hpp \
Modules/MScalarSUN/TrKinetic.hpp \
Modules/MDistil/PerambFromSolve.hpp \
Modules/MDistil/Distil.hpp \
Modules/MDistil/DistilVectors.hpp \
Modules/MDistil/Noises.hpp \
Modules/MDistil/LapEvec.hpp \
Modules/MDistil/Perambulator.hpp \
Modules/MDistil/DistilPar.hpp \
Modules/MAction/ZMobiusDWF.hpp \
Modules/MAction/ScaledDWF.hpp \
Modules/MAction/WilsonClover.hpp \

14
configure.ac
View File

@ -136,15 +136,15 @@ case ${ac_SFW_FP16} in
esac
############### SUMMIT JSRUN
AC_ARG_ENABLE([jsrun],
[AC_HELP_STRING([--enable-jsrun=yes|no], [enable IBMs jsrun resource manager for SUMMIT])],
[ac_JSRUN=${enable_jsrun}], [ac_JSRUN=no])
case ${ac_JSRUN} in
AC_ARG_ENABLE([summit],
[AC_HELP_STRING([--enable-summit=yes|no], [enable IBMs jsrun resource manager for SUMMIT])],
[ac_JSRUN=${enable_summit}], [ac_SUMMIT=no])
case ${ac_SUMMIT} in
no);;
yes)
AC_DEFINE([GRID_IBM_SUMMIT],[1],[Let JSRUN manage the GPU device allocation]);;
no);;
*)
AC_MSG_ERROR(["JSRUN option not supported ${ac_JSRUN}"]);;
AC_DEFINE([GRID_IBM_SUMMIT],[1],[Let JSRUN manage the GPU device allocation]);;
esac
############### Intel libraries
@ -255,7 +255,7 @@ AC_ARG_ENABLE([simd],[AC_HELP_STRING([--enable-simd=code],
AC_ARG_ENABLE([gen-simd-width],
[AS_HELP_STRING([--enable-gen-simd-width=size],
[size (in bytes) of the generic SIMD vectors (default: 32)])],
[size (in bytes) of the generic SIMD vectors (default: 64)])],
[ac_gen_simd_width=$enable_gen_simd_width],
[ac_gen_simd_width=64])

6
documentation/GridXcode/readme.md
View File

@ -119,13 +119,13 @@ Install [MacPorts][MacPorts] if you haven't done so already, and then install pa
These are the `portname`s for mandatory Grid libraries:
* git
* git-flow-avh
* gmp
* mpfr
and these are the `portname`s for optional Grid libraries:
* fftw-3
* fftw-3-single
* hdf5
* lapack
* doxygen
@ -369,7 +369,7 @@ Instead:
3. From a terminal session, locate and run your executable using `mpirun` (*the mangled name of the project build products will not be exactly the same as this example*):
`$GridPre/openmpi-3.1.3/bin/mpirun -np 2 ~/Library/Developer/Xcode/DerivedData/HelloGrid-fiyyuveptaqelbbvllomcgjyvghr/Build/Products/Debug/HelloGrid --grid 4.4.4.8 --mpi 1.1.1.2`
`$GridPre/openmpi/bin/mpirun -np 2 ~/Library/Developer/Xcode/DerivedData/HelloGrid-fiyyuveptaqelbbvllomcgjyvghr/Build/Products/Debug/HelloGrid --grid 4.4.4.8 --mpi 1.1.1.2`
The Xcode debugger will attach to the first process.

382
tests/hadrons/Test_distil.cc Normal file
View File

@ -0,0 +1,382 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Tests/Hadrons/Test_hadrons_distil.cc
Copyright (C) 2015-2019
Author: Felix Erben <ferben@ed.ac.uk>
Author: Michael Marshall <Michael.Marshall@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <typeinfo>
#include <Hadrons/Application.hpp>
#include <Hadrons/Modules.hpp>
using namespace Grid;
using namespace Hadrons;
// Very simple iterators for Eigen tensors
// The only way I could get these iterators to work is to put the begin() and end() functions in the Eigen namespace
// So if Eigen ever defines these, we'll have a conflict and have to change this
namespace Eigen {
template <typename ET>
inline typename std::enable_if<EigenIO::is_tensor<ET>::value, typename EigenIO::Traits<ET>::scalar_type *>::type
begin( ET & et ) { return reinterpret_cast<typename Grid::EigenIO::Traits<ET>::scalar_type *>(et.data()); }
template <typename ET>
inline typename std::enable_if<EigenIO::is_tensor<ET>::value, typename EigenIO::Traits<ET>::scalar_type *>::type
end( ET & et ) { return begin(et) + et.size() * EigenIO::Traits<ET>::count; }
}
/////////////////////////////////////////////////////////////
// Test creation of laplacian eigenvectors
/////////////////////////////////////////////////////////////
void test_Global(Application &application)
{
// global parameters
Application::GlobalPar globalPar;
globalPar.trajCounter.start = 1100;
globalPar.trajCounter.end = 1120;
globalPar.trajCounter.step = 20;
globalPar.runId = "test";
globalPar.graphFile = "";
globalPar.scheduleFile = "";
globalPar.saveSchedule = "false";
globalPar.parallelWriteMaxRetry = -1;
application.setPar(globalPar);
}
/////////////////////////////////////////////////////////////
// Create a random gauge with the correct name
/////////////////////////////////////////////////////////////
std::string test_Gauge(Application &application )
{
std::string sGaugeName{ "gauge" };
application.createModule<MGauge::Random>( sGaugeName );
return sGaugeName;
}
/////////////////////////////////////////////////////////////
// Test creation of laplacian eigenvectors
/////////////////////////////////////////////////////////////
void test_LapEvec(Application &application)
{
const char szModuleName[] = "LapEvec";
test_Gauge( application );
MDistil::LapEvecPar p;
p.gauge = "gauge";
p.Stout.steps = 3;
p.Stout.rho = 0.2;
p.Cheby.PolyOrder = 11;
p.Cheby.alpha = 0.55;
p.Cheby.beta = 35.5;
p.Lanczos.Nvec = 5;
p.Lanczos.Nk = 6;
p.Lanczos.Np = 2;
p.Lanczos.MaxIt = 1000;
p.Lanczos.resid = 1e-2;
p.Lanczos.IRLLog = 0;
application.createModule<MDistil::LapEvec>(szModuleName,p);
}
/////////////////////////////////////////////////////////////
// Test creation Solver
/////////////////////////////////////////////////////////////
std::string SolverName( const char * pSuffix = nullptr ) {
std::string sSolverName{ "CG" };
if( pSuffix && pSuffix[0] ) {
sSolverName.append( "_" );
sSolverName.append( pSuffix );
}
return sSolverName;
}
std::string test_Solver(Application &application, const char * pSuffix = nullptr )
{
std::string sActionName{ "DWF" };
if( pSuffix && pSuffix[0] ) {
sActionName.append( "_" );
sActionName.append( pSuffix );
}
MAction::DWF::Par actionPar;
actionPar.gauge = "gauge";
actionPar.Ls = 16;
actionPar.M5 = 1.8;
actionPar.mass = 0.005;
actionPar.boundary = "1 1 1 -1";
actionPar.twist = "0. 0. 0. 0.";
application.createModule<MAction::DWF>( sActionName, actionPar );
MSolver::RBPrecCG::Par solverPar;
solverPar.action = sActionName;
solverPar.residual = 1.0e-2;
solverPar.maxIteration = 10000;
std::string sSolverName{ SolverName( pSuffix ) };
application.createModule<MSolver::RBPrecCG>( sSolverName, solverPar );
return sSolverName;
}
/////////////////////////////////////////////////////////////
// DistilParameters
/////////////////////////////////////////////////////////////
std::string test_DPar(Application &application) {
MDistil::DistilParameters DPar;
DPar.nvec = 5;
DPar.nnoise = 1;
DPar.tsrc = 0;
DPar.LI = 5;
DPar.TI = 8;
DPar.SI = 4;
std::string sDParName{"DPar_l"};
application.createModule<MDistil::DistilPar>(sDParName,DPar);
return sDParName;
}
/////////////////////////////////////////////////////////////
// Noises
/////////////////////////////////////////////////////////////
std::string test_Noises(Application &application, const std::string &sNoiseBaseName ) {
MDistil::NoisesPar NoisePar;
NoisePar.DistilParams = "DPar_l";
NoisePar.NoiseFileName = "noise";
std::string sNoiseName{"noise"};
application.createModule<MDistil::Noises>(sNoiseName,NoisePar);
return sNoiseName;
}
/////////////////////////////////////////////////////////////
// Perambulators
/////////////////////////////////////////////////////////////
std::string PerambulatorName( const char * pszSuffix = nullptr )
{
std::string sPerambulatorName{ "Peramb" };
if( pszSuffix && pszSuffix[0] )
sPerambulatorName.append( pszSuffix );
return sPerambulatorName;
}
void test_LoadPerambulators( Application &application, const char * pszSuffix = nullptr )
{
std::string sModuleName{ "Peramb_load" };
MIO::LoadPerambulator::Par PerambPar;
PerambPar.PerambFileName = "Peramb";
PerambPar.DistilParams = "DPar_l";
test_Noises(application, sModuleName); // I want these written after solver stuff
application.createModule<MIO::LoadPerambulator>( sModuleName, PerambPar );
}
void test_Perambulators( Application &application, const char * pszSuffix = nullptr )
{
std::string sModuleName{ PerambulatorName( pszSuffix ) };
// Perambulator parameters
MDistil::Perambulator::Par PerambPar;
PerambPar.lapevec = "LapEvec";
PerambPar.PerambFileName = sModuleName;
PerambPar.solver = test_Solver( application, pszSuffix );
PerambPar.DistilParams = "DPar_l";
PerambPar.noise = "noise";
test_Noises(application, sModuleName); // I want these written after solver stuff
application.createModule<MDistil::Perambulator>( sModuleName, PerambPar );
}
/////////////////////////////////////////////////////////////
// DistilVectors
/////////////////////////////////////////////////////////////
void test_DistilVectors(Application &application, const char * pszSuffix = nullptr, const char * pszNvec = nullptr )
{
std::string sModuleName{"DistilVecs"};
if( pszSuffix )
sModuleName.append( pszSuffix );
std::string sPerambName{"Peramb"};
if( pszSuffix )
sPerambName.append( pszSuffix );
MDistil::DistilVectors::Par DistilVecPar;
DistilVecPar.noise = "noise";
DistilVecPar.rho = "rho";
DistilVecPar.phi = "phi";
DistilVecPar.perambulator = sPerambName;
DistilVecPar.lapevec = "LapEvec";
DistilVecPar.DistilParams = "DPar_l";
application.createModule<MDistil::DistilVectors>(sModuleName,DistilVecPar);
}
/////////////////////////////////////////////////////////////
// MesonSink
/////////////////////////////////////////////////////////////
void test_MesonSink(Application &application)
{
// DistilVectors parameters
MContraction::A2AMesonField::Par A2AMesonFieldPar;
//A2AMesonFieldPar.left="Peramb_unsmeared_sink";
A2AMesonFieldPar.left="g5phi";
A2AMesonFieldPar.right="Peramb_unsmeared_sink";
A2AMesonFieldPar.output="DistilFields";
A2AMesonFieldPar.gammas="Identity";
A2AMesonFieldPar.mom={"0 0 0"};
A2AMesonFieldPar.cacheBlock=2;
A2AMesonFieldPar.block=4;
application.createModule<MContraction::A2AMesonField>("DistilMesonSink",A2AMesonFieldPar);
}
/////////////////////////////////////////////////////////////
// MesonFields
/////////////////////////////////////////////////////////////
void test_MesonField(Application &application, const char * pszFileSuffix,
const char * pszObjectLeft = nullptr, const char * pszObjectRight = nullptr )
{
// DistilVectors parameters
if( pszObjectLeft == nullptr )
pszObjectLeft = pszFileSuffix;
if( pszObjectRight == nullptr )
pszObjectRight = pszObjectLeft;
MContraction::A2AMesonField::Par A2AMesonFieldPar;
A2AMesonFieldPar.left="";
A2AMesonFieldPar.right=A2AMesonFieldPar.left;
A2AMesonFieldPar.left.append( pszObjectLeft );
A2AMesonFieldPar.right.append( pszObjectRight );
A2AMesonFieldPar.output="MesonSinks";
A2AMesonFieldPar.output.append( pszFileSuffix );
A2AMesonFieldPar.gammas="Identity";
A2AMesonFieldPar.mom={"0 0 0"};
A2AMesonFieldPar.cacheBlock=2;
A2AMesonFieldPar.block=4;
std::string sObjectName{"DistilMesonField"};
sObjectName.append( pszFileSuffix );
application.createModule<MContraction::A2AMesonField>(sObjectName, A2AMesonFieldPar);
}
bool bNumber( int &ri, const char * & pstr, bool bGobbleWhiteSpace = true )
{
if( bGobbleWhiteSpace )
while( std::isspace(static_cast<unsigned char>(*pstr)) )
pstr++;
const char * p = pstr;
bool bMinus = false;
char c = * p++;
if( c == '+' )
c = * p++;
else if( c == '-' ) {
bMinus = true;
c = * p++;
}
int n = c - '0';
if( n < 0 || n > 9 )
return false;
while( * p >= '0' && * p <= '9' ) {
n = n * 10 + ( * p ) - '0';
p++;
}
if( bMinus )
n *= -1;
ri = n;
pstr = p;
return true;
}
int main(int argc, char *argv[])
{
// Decode command-line parameters. 1st one is which test to run
int iTestNum = -1;
for(int i = 1 ; i < argc ; i++ ) {
std::cout << "argv[" << i << "]=\"" << argv[i] << "\"" << std::endl;
const char * p = argv[i];
if( * p == '/' || * p == '-' ) {
p++;
char c = * p++;
switch(toupper(c)) {
case 'T':
if( bNumber( iTestNum, p ) ) {
std::cout << "Test " << iTestNum << " requested";
if( * p )
std::cout << " (ignoring trailer \"" << p << "\")";
std::cout << std::endl;
}
else
std::cout << "Invalid test \"" << &argv[i][2] << "\"" << std::endl;
break;
default:
std::cout << "Ignoring switch \"" << &argv[i][1] << "\"" << std::endl;
break;
}
}
}
// initialization //////////////////////////////////////////////////////////
Grid_init(&argc, &argv);
HadronsLogError.Active(GridLogError.isActive());
HadronsLogWarning.Active(GridLogWarning.isActive());
HadronsLogMessage.Active(GridLogMessage.isActive());
HadronsLogIterative.Active(GridLogIterative.isActive());
HadronsLogDebug.Active(GridLogDebug.isActive());
LOG(Message) << "Grid initialized" << std::endl;
// run setup ///////////////////////////////////////////////////////////////
Application application;
// For now perform free propagator test - replace this with distillation test(s)
LOG(Message) << "====== Creating xml for test " << iTestNum << " ======" << std::endl;
switch(iTestNum) {
default: // 0
LOG(Message) << "Computing Meson 2pt-function" << std::endl;
test_Global( application );
test_LapEvec( application );
test_DPar( application );
test_Perambulators( application );
test_DistilVectors( application );
test_MesonField( application, "Phi", "phi" );
test_MesonField( application, "Rho", "rho" );
break;
case 1:
LOG(Message) << "Computing Meson 2pt-function by loading perambulators" << std::endl;
test_Global( application );
test_LapEvec( application );
test_DPar( application );
test_LoadPerambulators( application );
test_DistilVectors( application, "_load" );
test_MesonField( application, "Phi", "phi" );
test_MesonField( application, "Rho", "rho" );
break;
}
// execution
static const char XmlFileName[] = "test_distil.xml";
application.saveParameterFile( XmlFileName );
const Grid::Coordinate &lat{GridDefaultLatt()};
if( lat.size() == 4 && lat[0] == 4 && lat[1] == 4 && lat[2] == 4 && lat[3] == 8 )
application.run();
else
LOG(Warning) << "The parameters in " << XmlFileName << " are designed to run on a laptop usid --grid 4.4.4.8" << std::endl;
// epilogue
LOG(Message) << "Grid is finalizing now" << std::endl;
Grid_finalize();
return EXIT_SUCCESS;
}