portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 07:55:35 +00:00
Code Releases Activity

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

Browse Source
This commit is contained in:
Peter Boyle 2019-12-09 03:53:01 -05:00
commit edd1c924eb
54 changed files with 4641 additions and 425 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

66
Grid/algorithms/CoarsenedMatrix.h
View File

@ -236,6 +236,39 @@ public:
Orthogonalise();
}
virtual void CreateSubspaceChebyshev(GridParallelRNG &RNG,LinearOperatorBase<FineField> &hermop,int nn=nbasis) {
RealD scale;
Chebyshev<FineField> Cheb(0.1,64.0,900);
FineField noise(FineGrid);
FineField Mn(FineGrid);
for(int b=0;b<nn;b++){
gaussian(RNG,noise);
scale = std::pow(norm2(noise),-0.5);
noise=noise*scale;
hermop.Op(noise,Mn); std::cout<<GridLogMessage << "noise ["<<b<<"] <n|MdagM|n> "<<norm2(Mn)<<std::endl;
Cheb(hermop,noise,Mn);
scale = std::pow(norm2(Mn),-0.5);
Mn=Mn*scale;
subspace[b] = Mn;
hermop.Op(Mn,noise); std::cout<<GridLogMessage << "filtered["<<b<<"] <f|MdagM|f> "<<norm2(noise)<<std::endl;
}
Orthogonalise();
}
};
// Fine Object == (per site) type of fine field
// nbasis == number of deflation vectors
@ -255,6 +288,7 @@ public:
////////////////////
Geometry geom;
GridBase * _grid;
int hermitian;
CartesianStencil<siteVector,siteVector,int> Stencil;
@ -271,10 +305,11 @@ public:
conformable(_grid,in.Grid());
conformable(in.Grid(),out.Grid());
RealD Nin = norm2(in);
SimpleCompressor<siteVector> compressor;
Stencil.HaloExchange(in,compressor);
auto in_v = in.View();
auto out_v = in.View();
auto out_v = out.View();
thread_for(ss,Grid()->oSites(),{
siteVector res = Zero();
siteVector nbr;
@ -296,19 +331,23 @@ public:
}
vstream(out_v[ss],res);
});
return norm2(out);
RealD Nout= norm2(out);
return Nout;
};
RealD Mdag (const CoarseVector &in, CoarseVector &out){
// // corresponds to Petrov-Galerkin coarsening
// return M(in,out);
// corresponds to Galerkin coarsening
CoarseVector tmp(Grid());
G5C(tmp, in);
M(tmp, out);
G5C(out, out);
return norm2(out);
RealD Mdag (const CoarseVector &in, CoarseVector &out)
{
if(hermitian) {
// corresponds to Petrov-Galerkin coarsening
return M(in,out);
} else {
// corresponds to Galerkin coarsening
CoarseVector tmp(Grid());
G5C(tmp, in);
M(tmp, out);
G5C(out, out);
return norm2(out);
}
};
void Mdir(const CoarseVector &in, CoarseVector &out, int dir, int disp){
@ -356,10 +395,11 @@ public:
};
CoarsenedMatrix(GridCartesian &CoarseGrid) :
CoarsenedMatrix(GridCartesian &CoarseGrid, int hermitian_=0) :
_grid(&CoarseGrid),
geom(CoarseGrid._ndimension),
hermitian(hermitian_),
Stencil(&CoarseGrid,geom.npoint,Even,geom.directions,geom.displacements,0),
A(geom.npoint,&CoarseGrid)
{

26
Grid/algorithms/LinearOperator.h
View File

@ -171,6 +171,32 @@ public:
}
};
template<class Matrix,class Field>
class NonHermitianLinearOperator : public LinearOperatorBase<Field> {
Matrix &_Mat;
public:
NonHermitianLinearOperator(Matrix &Mat): _Mat(Mat){};
// Support for coarsening to a multigrid
void OpDiag (const Field &in, Field &out) {
_Mat.Mdiag(in,out);
}
void OpDir (const Field &in, Field &out,int dir,int disp) {
_Mat.Mdir(in,out,dir,disp);
}
void Op (const Field &in, Field &out){
_Mat.M(in,out);
}
void AdjOp (const Field &in, Field &out){
_Mat.Mdag(in,out);
}
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
assert(0);
}
void HermOp(const Field &in, Field &out){
assert(0);
}
};
//////////////////////////////////////////////////////////
// Even Odd Schur decomp operators; there are several
// ways to introduce the even odd checkerboarding

38
Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h
View File

@ -86,7 +86,7 @@ public:
cp=GCRnStep(Linop,src,psi,rsq);
std::cout<<GridLogMessage<<"VPGCR("<<mmax<<","<<nstep<<") "<< steps <<" steps cp = "<<cp<<std::endl;
std::cout<<GridLogMessage<<"VPGCR("<<mmax<<","<<nstep<<") "<< steps <<" steps cp = "<<cp<<" target "<<rsq <<std::endl;
if(cp<rsq) {
@ -134,6 +134,9 @@ public:
std::vector<Field> p(mmax,grid);
std::vector<RealD> qq(mmax);
std::cout<<GridLogIterative<< " ************** "<< std::endl;
std::cout<<GridLogIterative<< " GCRnStep("<<nstep<<")"<<std::endl;
std::cout<<GridLogIterative<< " ************** "<< std::endl;
//////////////////////////////////
// initial guess x0 is taken as nonzero.
@ -143,25 +146,32 @@ public:
Linop.HermOpAndNorm(psi,Az,zAz,zAAz);
MatTimer.Stop();
LinalgTimer.Start();
r=src-Az;
LinalgTimer.Stop();
std::cout<<GridLogIterative<< " GCRnStep true residual r = src - A psi "<<norm2(r) <<std::endl;
/////////////////////
// p = Prec(r)
/////////////////////
std::cout<<GridLogIterative<< " GCRnStep apply preconditioner z= M^-1 r "<< std::endl;
std::cout<<GridLogIterative<< " --------------------------------------- "<< std::endl;
PrecTimer.Start();
Preconditioner(r,z);
PrecTimer.Stop();
std::cout<<GridLogIterative<< " --------------------------------------- "<< std::endl;
std::cout<<GridLogIterative<< " GCRnStep called Preconditioner z "<< norm2(z) <<std::endl;
MatTimer.Start();
Linop.HermOp(z,tmp);
MatTimer.Stop();
// MatTimer.Start();
// Linop.HermOp(z,tmp);
// MatTimer.Stop();
LinalgTimer.Start();
ttmp=tmp;
tmp=tmp-r;
LinalgTimer.Stop();
// LinalgTimer.Start();
// ttmp=tmp;
// tmp=tmp-r;
// LinalgTimer.Stop();
/*
std::cout<<GridLogMessage<<r<<std::endl;
@ -175,10 +185,12 @@ public:
MatTimer.Stop();
LinalgTimer.Start();
//p[0],q[0],qq[0]
p[0]= z;
q[0]= Az;
qq[0]= zAAz;
std::cout<<GridLogIterative<< " GCRnStep p0=z, q0 = A p0 " <<std::endl;
cp =norm2(r);
LinalgTimer.Stop();
@ -200,24 +212,26 @@ public:
cp = axpy_norm(r,-a,q[peri_k],r);
LinalgTimer.Stop();
std::cout<<GridLogMessage<< " VPGCR_step["<<steps<<"] resid " << cp << " target " <<rsq<<std::endl;
if((k==nstep-1)||(cp<rsq)){
return cp;
}
std::cout<<GridLogMessage<< " VPGCR_step["<<steps<<"] resid " <<sqrt(cp/rsq)<<std::endl;
std::cout<<GridLogIterative<< " GCRnStep apply preconditioner z= M^-1 r "<< std::endl;
std::cout<<GridLogIterative<< " --------------------------------------- "<< std::endl;
PrecTimer.Start();
Preconditioner(r,z);// solve Az = r
PrecTimer.Stop();
std::cout<<GridLogIterative<< " --------------------------------------- "<< std::endl;
std::cout<<GridLogIterative<< " GCRnStep called Preconditioner z "<< norm2(z) <<std::endl;
MatTimer.Start();
Linop.HermOpAndNorm(z,Az,zAz,zAAz);
Linop.HermOp(z,tmp);
MatTimer.Stop();
LinalgTimer.Start();
tmp=tmp-r;
std::cout<<GridLogMessage<< " Preconditioner resid " <<sqrt(norm2(tmp)/norm2(r))<<std::endl;
q[peri_kp]=Az;
p[peri_kp]=z;

371
Grid/algorithms/iterative/QuasiMinimalResidual.h Normal file
View File

@ -0,0 +1,371 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/algorithmsf/iterative/QuasiMinimalResidual.h
Copyright (C) 2019
Author: Peter Boyle <paboyle@ph.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 */
#pragma once
NAMESPACE_BEGIN(Grid);
template<class Field>
RealD innerG5ProductReal(Field &l, Field &r)
{
Gamma G5(Gamma::Algebra::Gamma5);
Field tmp(l.Grid());
// tmp = G5*r;
G5R5(tmp,r);
ComplexD ip =innerProduct(l,tmp);
std::cout << "innerProductRealG5R5 "<<ip<<std::endl;
return ip.real();
}
template<class Field>
class QuasiMinimalResidual : public OperatorFunction<Field> {
public:
using OperatorFunction<Field>::operator();
bool ErrorOnNoConverge;
RealD Tolerance;
Integer MaxIterations;
Integer IterationCount;
QuasiMinimalResidual(RealD tol,
Integer maxit,
bool err_on_no_conv = true)
: Tolerance(tol)
, MaxIterations(maxit)
, ErrorOnNoConverge(err_on_no_conv)
{};
#if 1
void operator()(LinearOperatorBase<Field> &LinOp, const Field &b, Field &x)
{
RealD resid;
IterationCount=0;
RealD rho, rho_1, xi, gamma, gamma_1, theta, theta_1;
RealD eta, delta, ep, beta;
GridBase *Grid = b.Grid();
Field r(Grid), d(Grid), s(Grid);
Field v(Grid), w(Grid), y(Grid), z(Grid);
Field v_tld(Grid), w_tld(Grid), y_tld(Grid), z_tld(Grid);
Field p(Grid), q(Grid), p_tld(Grid);
Real normb = norm2(b);
LinOp.Op(x,r); r = b - r;
assert(normb> 0.0);
resid = norm2(r)/normb;
if (resid <= Tolerance) {
return;
}
v_tld = r;
y = v_tld;
rho = norm2(y);
// Take Gamma5 conjugate
// Gamma G5(Gamma::Algebra::Gamma5);
// G5R5(w_tld,r);
// w_tld = G5* v_tld;
w_tld=v_tld;
z = w_tld;
xi = norm2(z);
gamma = 1.0;
eta = -1.0;
theta = 0.0;
for (int i = 1; i <= MaxIterations; i++) {
// Breakdown tests
assert( rho != 0.0);
assert( xi != 0.0);
v = (1. / rho) * v_tld;
y = (1. / rho) * y;
w = (1. / xi) * w_tld;
z = (1. / xi) * z;
ComplexD Zdelta = innerProduct(z, y); // Complex?
std::cout << "Zdelta "<<Zdelta<<std::endl;
delta = Zdelta.real();
y_tld = y;
z_tld = z;
if (i > 1) {
p = y_tld - (xi * delta / ep) * p;
q = z_tld - (rho * delta / ep) * q;
} else {
p = y_tld;
q = z_tld;
}
LinOp.Op(p,p_tld); // p_tld = A * p;
ComplexD Zep = innerProduct(q, p_tld);
ep=Zep.real();
std::cout << "Zep "<<Zep <<std::endl;
// Complex Audit
assert(abs(ep)>0);
beta = ep / delta;
assert(abs(beta)>0);
v_tld = p_tld - beta * v;
y = v_tld;
rho_1 = rho;
rho = norm2(y);
LinOp.AdjOp(q,w_tld);
w_tld = w_tld - beta * w;
z = w_tld;
xi = norm2(z);
gamma_1 = gamma;
theta_1 = theta;
theta = rho / (gamma_1 * beta);
gamma = 1.0 / sqrt(1.0 + theta * theta);
std::cout << "theta "<<theta<<std::endl;
std::cout << "gamma "<<gamma<<std::endl;
assert(abs(gamma)> 0.0);
eta = -eta * rho_1 * gamma* gamma / (beta * gamma_1 * gamma_1);
if (i > 1) {
d = eta * p + (theta_1 * theta_1 * gamma * gamma) * d;
s = eta * p_tld + (theta_1 * theta_1 * gamma * gamma) * s;
} else {
d = eta * p;
s = eta * p_tld;
}
x =x+d; // update approximation vector
r =r-s; // compute residual
if ((resid = norm2(r) / normb) <= Tolerance) {
return;
}
std::cout << "Iteration "<<i<<" resid " << resid<<std::endl;
}
assert(0);
return; // no convergence
}
#else
// QMRg5 SMP thesis
void operator()(LinearOperatorBase<Field> &LinOp, const Field &b, Field &x)
{
// Real scalars
GridBase *grid = b.Grid();
Field r(grid);
Field p_m(grid), p_m_minus_1(grid), p_m_minus_2(grid);
Field v_m(grid), v_m_minus_1(grid), v_m_plus_1(grid);
Field tmp(grid);
RealD w;
RealD z1, z2;
RealD delta_m, delta_m_minus_1;
RealD c_m_plus_1, c_m, c_m_minus_1;
RealD s_m_plus_1, s_m, s_m_minus_1;
RealD alpha, beta, gamma, epsilon;
RealD mu, nu, rho, theta, xi, chi;
RealD mod2r, mod2b;
RealD tau2, target2;
mod2b=norm2(b);
/////////////////////////
// Initial residual
/////////////////////////
LinOp.Op(x,tmp);
r = b - tmp;
/////////////////////////
// \mu = \rho = |r_0|
/////////////////////////
mod2r = norm2(r);
rho = sqrt( mod2r);
mu=rho;
std::cout << "QuasiMinimalResidual rho "<< rho<<std::endl;
/////////////////////////
// Zero negative history
/////////////////////////
v_m_plus_1 = Zero();
v_m_minus_1 = Zero();
p_m_minus_1 = Zero();
p_m_minus_2 = Zero();
// v0
v_m = (1.0/rho)*r;
/////////////////////////
// Initial coeffs
/////////////////////////
delta_m_minus_1 = 1.0;
c_m_minus_1 = 1.0;
c_m = 1.0;
s_m_minus_1 = 0.0;
s_m = 0.0;
/////////////////////////
// Set up convergence check
/////////////////////////
tau2 = mod2r;
target2 = mod2b * Tolerance*Tolerance;
for(int iter = 0 ; iter < MaxIterations; iter++){
/////////////////////////
// \delta_m = (v_m, \gamma_5 v_m)
/////////////////////////
delta_m = innerG5ProductReal(v_m,v_m);
std::cout << "QuasiMinimalResidual delta_m "<< delta_m<<std::endl;
/////////////////////////
// tmp = A v_m
/////////////////////////
LinOp.Op(v_m,tmp);
/////////////////////////
// \alpha = (v_m, \gamma_5 temp) / \delta_m
/////////////////////////
alpha = innerG5ProductReal(v_m,tmp);
alpha = alpha/delta_m ;
std::cout << "QuasiMinimalResidual alpha "<< alpha<<std::endl;
/////////////////////////
// \beta = \rho \delta_m / \delta_{m-1}
/////////////////////////
beta = rho * delta_m / delta_m_minus_1;
std::cout << "QuasiMinimalResidual beta "<< beta<<std::endl;
/////////////////////////
// \tilde{v}_{m+1} = temp - \alpha v_m - \beta v_{m-1}
/////////////////////////
v_m_plus_1 = tmp - alpha*v_m - beta*v_m_minus_1;
///////////////////////////////
// \rho = || \tilde{v}_{m+1} ||
///////////////////////////////
rho = sqrt( norm2(v_m_plus_1) );
std::cout << "QuasiMinimalResidual rho "<< rho<<std::endl;
///////////////////////////////
// v_{m+1} = \tilde{v}_{m+1}
///////////////////////////////
v_m_plus_1 = (1.0 / rho) * v_m_plus_1;
////////////////////////////////
// QMR recurrence coefficients.
////////////////////////////////
theta = s_m_minus_1 * beta;
gamma = c_m_minus_1 * beta;
epsilon = c_m * gamma + s_m * alpha;
xi = -s_m * gamma + c_m * alpha;
nu = sqrt( xi*xi + rho*rho );
c_m_plus_1 = fabs(xi) / nu;
if ( xi == 0.0 ) {
s_m_plus_1 = 1.0;
} else {
s_m_plus_1 = c_m_plus_1 * rho / xi;
}
chi = c_m_plus_1 * xi + s_m_plus_1 * rho;
std::cout << "QuasiMinimalResidual coeffs "<< theta <<" "<<gamma<<" "<< epsilon<<" "<< xi<<" "<< nu<<std::endl;
std::cout << "QuasiMinimalResidual coeffs "<< chi <<std::endl;
////////////////////////////////
//p_m=(v_m - \epsilon p_{m-1} - \theta p_{m-2}) / \chi
////////////////////////////////
p_m = (1.0/chi) * v_m - (epsilon/chi) * p_m_minus_1 - (theta/chi) * p_m_minus_2;
////////////////////////////////////////////////////////////////
// \psi = \psi + c_{m+1} \mu p_m
////////////////////////////////////////////////////////////////
x = x + ( c_m_plus_1 * mu ) * p_m;
////////////////////////////////////////
//
////////////////////////////////////////
mu = -s_m_plus_1 * mu;
delta_m_minus_1 = delta_m;
c_m_minus_1 = c_m;
c_m = c_m_plus_1;
s_m_minus_1 = s_m;
s_m = s_m_plus_1;
////////////////////////////////////
// Could use pointer swizzle games.
////////////////////////////////////
v_m_minus_1 = v_m;
v_m = v_m_plus_1;
p_m_minus_2 = p_m_minus_1;
p_m_minus_1 = p_m;
/////////////////////////////////////
// Convergence checks
/////////////////////////////////////
z1 = RealD(iter+1.0);
z2 = z1 + 1.0;
tau2 = tau2 *( z2 / z1 ) * s_m * s_m;
std::cout << " QuasiMinimumResidual iteration "<< iter<<std::endl;
std::cout << " QuasiMinimumResidual tau bound "<< tau2<<std::endl;
// Compute true residual
mod2r = tau2;
if ( 1 || (tau2 < (100.0 * target2)) ) {
LinOp.Op(x,tmp);
r = b - tmp;
mod2r = norm2(r);
std::cout << " QuasiMinimumResidual true residual is "<< mod2r<<std::endl;
}
if ( mod2r < target2 ) {
std::cout << " QuasiMinimumResidual has converged"<<std::endl;
return;
}
}
}
#endif
};
NAMESPACE_END(Grid);

9
Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h
View File

@ -383,10 +383,11 @@ void CayleyFermion5D<Impl>::MeooeDag (const FermionField &psi, FermionField &
}
template<class Impl>
void CayleyFermion5D<Impl>::Mdir (const FermionField &psi, FermionField &chi,int dir,int disp){
Meo5D(psi,this->tmp());
// Apply 4d dslash fragment
this->DhopDir(this->tmp(),chi,dir,disp);
void CayleyFermion5D<Impl>::Mdir (const FermionField &psi, FermionField &chi,int dir,int disp)
{
FermionField tmp(psi.Grid());
Meo5D(psi,tmp);
this->DhopDir(tmp,chi,dir,disp);
}
// force terms; five routines; default to Dhop on diagonal
template<class Impl>

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]))()()();
}
}
}

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

@ -41,7 +41,10 @@ public:
typedef typename FImpl::SitePropagator pobj;
typedef typename ComplexField::vector_object vobj;
typedef Lattice<iSpinMatrix<typename FImpl::Simd>> SpinMatrixField;
typedef typename SpinMatrixField::vector_object sobj;
static const int epsilon[6][3] ;
static const Complex epsilon_sgn[6];
@ -81,6 +84,69 @@ public:
const char * quarks_right,
const int parity,
robj &result);
private:
template <class mobj, class mobj2, class robj>
static void Sigma_to_Nucleon_Q1_Eye_site(const mobj &Dq_loop,
const mobj2 &Du_spec,
const mobj &Dd_tf,
const mobj &Ds_ti,
const Gamma Gamma_H,
const Gamma GammaB_sigma,
const Gamma GammaB_nucl,
robj &result);
template <class mobj, class mobj2, class robj>
static void Sigma_to_Nucleon_Q1_NonEye_site(const mobj &Du_ti,
const mobj &Du_tf,
const mobj2 &Du_spec,
const mobj &Dd_tf,
const mobj &Ds_ti,
const Gamma Gamma_H,
const Gamma GammaB_sigma,
const Gamma GammaB_nucl,
robj &result);
template <class mobj, class mobj2, class robj>
static void Sigma_to_Nucleon_Q2_Eye_site(const mobj &Dq_loop,
const mobj2 &Du_spec,
const mobj &Dd_tf,
const mobj &Ds_ti,
const Gamma Gamma_H,
const Gamma GammaB_sigma,
const Gamma GammaB_nucl,
robj &result);
template <class mobj, class mobj2, class robj>
static void Sigma_to_Nucleon_Q2_NonEye_site(const mobj &Du_ti,
const mobj &Du_tf,
const mobj2 &Du_spec,
const mobj &Dd_tf,
const mobj &Ds_ti,
const Gamma Gamma_H,
const Gamma GammaB_sigma,
const Gamma GammaB_nucl,
robj &result);
public:
template <class mobj>
static void Sigma_to_Nucleon_Eye(const PropagatorField &qq_loop,
const mobj &Du_spec,
const PropagatorField &qd_tf,
const PropagatorField &qs_ti,
const Gamma Gamma_H,
const Gamma GammaB_sigma,
const Gamma GammaB_nucl,
const std::string op,
SpinMatrixField &stn_corr);
template <class mobj>
static void Sigma_to_Nucleon_NonEye(const PropagatorField &qq_ti,
const PropagatorField &qq_tf,
const mobj &Du_spec,
const PropagatorField &qd_tf,
const PropagatorField &qs_ti,
const Gamma Gamma_H,
const Gamma GammaB_sigma,
const Gamma GammaB_nucl,
const std::string op,
SpinMatrixField &stn_corr);
};
template <class FImpl>
@ -254,4 +320,305 @@ void BaryonUtils<FImpl>::ContractBaryons_Sliced(const mobj &D1,
result=Zero();
baryon_site(D1,D2,D3,GammaA_left,GammaB_left,GammaA_right,GammaB_right,parity,wick_contraction,result);
}
/***********************************************************************
* End of Baryon 2pt-function code. *
* *
* The following code is for Sigma -> N rare hypeon decays *
**********************************************************************/
/* Dq_loop is a quark line from t_H to t_H
* Du_spec is a quark line from t_i to t_f
* Dd_tf is a quark line from t_f to t_H
* Ds_ti is a quark line from t_i to t_H */
template <class FImpl>
template <class mobj, class mobj2, class robj>
void BaryonUtils<FImpl>::Sigma_to_Nucleon_Q1_Eye_site(const mobj &Dq_loop,
const mobj2 &Du_spec,
const mobj &Dd_tf,
const mobj &Ds_ti,
const Gamma Gamma_H,
const Gamma GammaB_sigma,
const Gamma GammaB_nucl,
robj &result)
{
Gamma g5(Gamma::Algebra::Gamma5);
auto DuG = Du_spec * GammaB_nucl;
// Gamma^B * Ds * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5)
auto GDsGDd = GammaB_sigma * Ds_ti * Gamma_H * g5 * adj(Dd_tf) * g5;
// Dq_loop * \gamma_\mu^L
auto DqG = Dq_loop * Gamma_H;
for (int ie_n=0; ie_n < 6 ; ie_n++){
int a_n = epsilon[ie_n][0]; //a
int b_n = epsilon[ie_n][1]; //b
int c_n = epsilon[ie_n][2]; //c
for (int ie_s=0; ie_s < 6 ; ie_s++){
int a_s = epsilon[ie_s][0]; //a'
int b_s = epsilon[ie_s][1]; //b'
int c_s = epsilon[ie_s][2]; //c'
for (int alpha_s=0; alpha_s<Ns; alpha_s++){
for (int beta_n=0; beta_n<Ns; beta_n++){
auto GDsGDd_ab_bb = GDsGDd()(alpha_s,beta_n)(b_s,b_n);
for (int tau2=0; tau2<Ns; tau2++){
for (int j=0; j<Nc; j++){
auto DqG_tt_jj = DqG()(tau2,tau2)(j,j);
auto ee_GDGDDG = epsilon_sgn[ie_n] * epsilon_sgn[ie_s] * GDsGDd_ab_bb * DqG_tt_jj;
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_GDGDDG * DuG()(alpha_s, beta_n)(a_s,a_n) * Du_spec()(gamma_s,gamma_n)(c_s,c_n);
result()(gamma_s,gamma_n)() -= ee_GDGDDG * DuG()(gamma_s, beta_n)(c_s,a_n) * Du_spec()(alpha_s,gamma_n)(a_s,c_n);
}}
}}
}}
}
}
}
/* Du_ti is a quark line from t_i to t_H
* Du_tf is a quark line from t_f to t_H
* Du_spec is a quark line from t_i to t_f
* Dd_tf is a quark line from t_f to t_H
* Ds_ti is a quark line from t_i to t_H */
template <class FImpl>
template <class mobj, class mobj2, class robj>
void BaryonUtils<FImpl>::Sigma_to_Nucleon_Q1_NonEye_site(const mobj &Du_ti,
const mobj &Du_tf,
const mobj2 &Du_spec,
const mobj &Dd_tf,
const mobj &Ds_ti,
const Gamma Gamma_H,
const Gamma GammaB_sigma,
const Gamma GammaB_nucl,
robj &result)
{
Gamma g5(Gamma::Algebra::Gamma5);
auto DuG = Du_spec * GammaB_nucl;
auto adjDu = g5 * adj(Du_tf) * g5;
auto adjDuG = adjDu * GammaB_nucl;
// Gamma^B * Ds * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5)
auto GDsGDd = GammaB_sigma * Ds_ti * Gamma_H * g5 * adj(Dd_tf) * g5;
// Dq_loop * \gamma_\mu^L
auto DuGH = Du_ti * Gamma_H;
for (int ie_n=0; ie_n < 6 ; ie_n++){
int a_n = epsilon[ie_n][0]; //a
int b_n = epsilon[ie_n][1]; //b
int c_n = epsilon[ie_n][2]; //c
for (int ie_s=0; ie_s < 6 ; ie_s++){
int a_s = epsilon[ie_s][0]; //a'
int b_s = epsilon[ie_s][1]; //b'
int c_s = epsilon[ie_s][2]; //c'
for (int alpha_s=0; alpha_s<Ns; alpha_s++){
for (int beta_n=0; beta_n<Ns; beta_n++){
auto GDsGDd_ab_bb = GDsGDd()(alpha_s,beta_n)(b_s,b_n);
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 = 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 = 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);
result()(gamma_s,gamma_n)() -= ee_GDGDDG_a * adjDuG()(tau2, beta_n)(j,a_n) * Du_spec()(gamma_s,gamma_n)(c_s,c_n);
result()(gamma_s,gamma_n)() -= ee_GDGDDG_c * DuG()(alpha_s, beta_n)(a_s,a_n) * adjDu()(tau2,gamma_n)(j,c_n);
}
}
}}
}}
}
}
}
//Equivalent to "One-trace"
/* Dq_loop is a quark line from t_H to t_H
* Du_spec is a quark line from t_i to t_f
* Dd_tf is a quark line from t_f to t_H
* Ds_ti is a quark line from t_i to t_H */
template <class FImpl>
template <class mobj, class mobj2, class robj>
void BaryonUtils<FImpl>::Sigma_to_Nucleon_Q2_Eye_site(const mobj &Dq_loop,
const mobj2 &Du_spec,
const mobj &Dd_tf,
const mobj &Ds_ti,
const Gamma Gamma_H,
const Gamma GammaB_sigma,
const Gamma GammaB_nucl,
robj &result)
{
Gamma g5(Gamma::Algebra::Gamma5);
auto DuG = Du_spec * GammaB_nucl;
// Gamma^B * Ds * \gamma_\mu^L
auto GDsG = GammaB_sigma * Ds_ti * Gamma_H;
// Dq_loop * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5)
auto DqGDd = Dq_loop * Gamma_H * g5 * adj(Dd_tf) * g5;
for (int ie_n=0; ie_n < 6 ; ie_n++){
int a_n = epsilon[ie_n][0]; //a
int b_n = epsilon[ie_n][1]; //b
int c_n = epsilon[ie_n][2]; //c
for (int ie_s=0; ie_s < 6 ; ie_s++){
int a_s = epsilon[ie_s][0]; //a'
int b_s = epsilon[ie_s][1]; //b'
int c_s = epsilon[ie_s][2]; //c'
for (int alpha_s=0; alpha_s<Ns; alpha_s++){
for (int tau=0; tau<Ns; tau++){
for (int i=0; i<Nc; i++){
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 = 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);
result()(gamma_s,gamma_n)() += ee_GDGDGD * DuG()(gamma_s, beta_n)(c_s,a_n) * Du_spec()(alpha_s,gamma_n)(a_s,c_n);
}}
}
}}}
}
}
}
/* Du_ti is a quark line from t_i to t_H
* Du_tf is a quark line from t_f to t_H
* Du_spec is a quark line from t_i to t_f
* Dd_tf is a quark line from t_f to t_H
* Ds_ti is a quark line from t_i to t_H */
template <class FImpl>
template <class mobj, class mobj2, class robj>
void BaryonUtils<FImpl>::Sigma_to_Nucleon_Q2_NonEye_site(const mobj &Du_ti,
const mobj &Du_tf,
const mobj2 &Du_spec,
const mobj &Dd_tf,
const mobj &Ds_ti,
const Gamma Gamma_H,
const Gamma GammaB_sigma,
const Gamma GammaB_nucl,
robj &result)
{
Gamma g5(Gamma::Algebra::Gamma5);
auto DuG = Du_spec * GammaB_nucl;
auto adjDu = g5 * adj(Du_tf) * g5;
auto adjDuG = adjDu * GammaB_nucl;
// Gamma^B * Ds * \gamma_\mu^L
auto GDsG = GammaB_sigma * Ds_ti * Gamma_H;
// Du * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5)
auto DuGDd = Du_ti * Gamma_H * g5 * adj(Dd_tf) * g5;
for (int ie_n=0; ie_n < 6 ; ie_n++){
int a_n = epsilon[ie_n][0]; //a
int b_n = epsilon[ie_n][1]; //b
int c_n = epsilon[ie_n][2]; //c
for (int ie_s=0; ie_s < 6 ; ie_s++){
int a_s = epsilon[ie_s][0]; //a'
int b_s = epsilon[ie_s][1]; //b'
int c_s = epsilon[ie_s][2]; //c'
for (int alpha_s=0; alpha_s<Ns; alpha_s++){
for (int tau=0; tau<Ns; tau++){
for (int i=0; i<Nc; i++){
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 = 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 = 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);
result()(gamma_s,gamma_n)() += ee_GDGDGD_a * adjDuG()(tau, beta_n)(i,a_n) * Du_spec()(gamma_s,gamma_n)(c_s,c_n);
result()(gamma_s,gamma_n)() += ee_GDGDGD_c * DuG()(alpha_s, beta_n)(a_s,a_n) * adjDu()(tau,gamma_n)(i,c_n);
}
}
}
}}}
}
}
}
template<class FImpl>
template <class mobj>
void BaryonUtils<FImpl>::Sigma_to_Nucleon_Eye(const PropagatorField &qq_loop,
const mobj &Du_spec,
const PropagatorField &qd_tf,
const PropagatorField &qs_ti,
const Gamma Gamma_H,
const Gamma GammaB_sigma,
const Gamma GammaB_nucl,
const std::string op,
SpinMatrixField &stn_corr)
{
GridBase *grid = qs_ti.Grid();
auto vcorr= stn_corr.View();
auto vq_loop = qq_loop.View();
auto vd_tf = qd_tf.View();
auto vs_ti = qs_ti.View();
// accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
thread_for(ss,grid->oSites(),{
auto Dq_loop = vq_loop[ss];
auto Dd_tf = vd_tf[ss];
auto Ds_ti = vs_ti[ss];
sobj result=Zero();
if(op == "Q1"){
Sigma_to_Nucleon_Q1_Eye_site(Dq_loop,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result);
} else if(op == "Q2"){
Sigma_to_Nucleon_Q2_Eye_site(Dq_loop,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result);
} else {
assert(0 && "Weak Operator not correctly specified");
}
vcorr[ss] = result;
} );//end loop over lattice sites
}
template<class FImpl>
template <class mobj>
void BaryonUtils<FImpl>::Sigma_to_Nucleon_NonEye(const PropagatorField &qq_ti,
const PropagatorField &qq_tf,
const mobj &Du_spec,
const PropagatorField &qd_tf,
const PropagatorField &qs_ti,
const Gamma Gamma_H,
const Gamma GammaB_sigma,
const Gamma GammaB_nucl,
const std::string op,
SpinMatrixField &stn_corr)
{
GridBase *grid = qs_ti.Grid();
auto vcorr= stn_corr.View();
auto vq_ti = qq_ti.View();
auto vq_tf = qq_tf.View();
auto vd_tf = qd_tf.View();
auto vs_ti = qs_ti.View();
// accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
thread_for(ss,grid->oSites(),{
auto Dq_ti = vq_ti[ss];
auto Dq_tf = vq_tf[ss];
auto Dd_tf = vd_tf[ss];
auto Ds_ti = vs_ti[ss];
sobj result=Zero();
if(op == "Q1"){
Sigma_to_Nucleon_Q1_NonEye_site(Dq_ti,Dq_tf,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result);
} else if(op == "Q2"){
Sigma_to_Nucleon_Q2_NonEye_site(Dq_ti,Dq_tf,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result);
} else {
assert(0 && "Weak Operator not correctly specified");
}
vcorr[ss] = result;
} );//end loop over lattice sites
}
NAMESPACE_END(Grid);

106
Hadrons/DilutedNoise.hpp
View File

@ -52,6 +52,7 @@ public:
const std::vector<FermionField> & getNoise(void) const;
const FermionField & operator[](const unsigned int i) const;
FermionField & operator[](const unsigned int i);
void normalise(Real norm);
void resize(const unsigned int nNoise);
unsigned int size(void) const;
GridCartesian *getGrid(void) const;
@ -93,6 +94,21 @@ private:
unsigned int nSrc_;
};
template <typename FImpl>
class SparseSpinColorDiagonalNoise: public DilutedNoise<FImpl>
{
public:
typedef typename FImpl::FermionField FermionField;
public:
// constructor/destructor
SparseSpinColorDiagonalNoise(GridCartesian *g, unsigned int n_src, unsigned int n_sparse);
virtual ~SparseSpinColorDiagonalNoise(void) = default;
// generate noise
virtual void generateNoise(GridParallelRNG &rng);
private:
unsigned int nSrc_;
unsigned int nSparse_;
};
/******************************************************************************
* DilutedNoise template implementation *
@ -138,6 +154,15 @@ DilutedNoise<FImpl>::operator[](const unsigned int i)
return noise_[i];
}
template <typename FImpl>
void DilutedNoise<FImpl>::normalise(Real norm)
{
for(int i=0;i<noise_.size();i++)
{
noise_[i] = norm*noise_[i];
}
}
template <typename FImpl>
void DilutedNoise<FImpl>::resize(const unsigned int nNoise)
{
@ -245,6 +270,87 @@ void FullVolumeSpinColorDiagonalNoise<FImpl>::generateNoise(GridParallelRNG &rng
}
}
/******************************************************************************
* SparseSpinColorDiagonalNoise template implementation *
******************************************************************************/
template <typename FImpl>
SparseSpinColorDiagonalNoise<FImpl>::
SparseSpinColorDiagonalNoise(GridCartesian *g, unsigned int nSrc, unsigned int nSparse)
: DilutedNoise<FImpl>(g, nSrc*Ns*FImpl::Dimension), nSrc_(nSrc), nSparse_(nSparse)
{}
template <typename FImpl>
void SparseSpinColorDiagonalNoise<FImpl>::generateNoise(GridParallelRNG &rng)
{
typedef decltype(peekColour((*this)[0], 0)) SpinField;
auto &noise = *this;
auto g = this->getGrid();
auto nd = g->GlobalDimensions().size();
auto nc = FImpl::Dimension;
LatticeInteger coor(g), coorTot(g); coorTot = 0.;
Complex shift(1., 1.);
LatticeComplex eta(g), etaSparse(g);
SpinField etas(g);
unsigned int i = 0;
unsigned int j = 0;
unsigned int nSrc_ec;
if(nSrc_%nSparse_==0)
{
nSrc_ec = nSrc_/nSparse_;
}
else
{
nSrc_ec = (nSrc_ - nSrc_%nSparse_)/nSparse_;
}
for (unsigned int n = 0; n < nSrc_; ++n)
{
bernoulli(rng, eta);
eta = (2.*eta - shift)*(1./::sqrt(2.));
if(nSparse_ != 1)
{
assert(g->GlobalDimensions()[1]%nSparse_ == 0);
// # 0 # 0
// 0 # 0 #
// # 0 # 0
// 0 # 0 #
coorTot = 0;
for(unsigned int d = 0; d < nd; ++d)
{
LatticeCoordinate(coor, d);
coorTot = coorTot + coor;
}
coorTot = coorTot + j;
eta = where(mod(coorTot,nSparse_), 0.*eta, eta);
}
for (unsigned int s = 0; s < Ns; ++s)
{
etas = Zero();
pokeSpin(etas, eta, s);
for (unsigned int c = 0; c < nc; ++c)
{
noise[i] = Zero();
pokeColour(noise[i], etas, c);
i++;
/**/
}
}
((n+1)%nSrc_ec == 0) ? j++: 0;
}
Real norm = sqrt(1./nSrc_ec);
this->normalise(norm);
}
END_HADRONS_NAMESPACE
#endif // Hadrons_DilutedNoise_hpp_

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:

137
Hadrons/Modules.hpp
View File

@ -1,73 +1,86 @@
#include <Hadrons/Modules/MFermion/GaugeProp.hpp>
#include <Hadrons/Modules/MFermion/EMLepton.hpp>
#include <Hadrons/Modules/MFermion/FreeProp.hpp>
#include <Hadrons/Modules/MIO/LoadCosmHol.hpp>
#include <Hadrons/Modules/MIO/LoadEigenPack.hpp>
#include <Hadrons/Modules/MIO/LoadA2AVectors.hpp>
#include <Hadrons/Modules/MIO/LoadA2AMatrixDiskVector.hpp>
#include <Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp>
#include <Hadrons/Modules/MIO/LoadNersc.hpp>
#include <Hadrons/Modules/MIO/LoadBinary.hpp>
#include <Hadrons/Modules/MAction/ZMobiusDWF.hpp>
#include <Hadrons/Modules/MAction/MobiusDWF.hpp>
#include <Hadrons/Modules/MAction/Wilson.hpp>
#include <Hadrons/Modules/MAction/DWF.hpp>
#include <Hadrons/Modules/MAction/WilsonClover.hpp>
#include <Hadrons/Modules/MAction/MobiusDWF.hpp>
#include <Hadrons/Modules/MAction/ScaledDWF.hpp>
#include <Hadrons/Modules/MUtilities/RandomVectors.hpp>
#include <Hadrons/Modules/MUtilities/PrecisionCast.hpp>
#include <Hadrons/Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp>
#include <Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp>
#include <Hadrons/Modules/MContraction/DiscLoop.hpp>
#include <Hadrons/Modules/MContraction/Meson.hpp>
#include <Hadrons/Modules/MContraction/WeakMesonDecayKl2.hpp>
#include <Hadrons/Modules/MContraction/A2AMesonField.hpp>
#include <Hadrons/Modules/MContraction/WeakNonEye3pt.hpp>
#include <Hadrons/Modules/MContraction/Gamma3pt.hpp>
#include <Hadrons/Modules/MAction/Wilson.hpp>
#include <Hadrons/Modules/MAction/WilsonClover.hpp>
#include <Hadrons/Modules/MAction/ZMobiusDWF.hpp>
#include <Hadrons/Modules/MContraction/A2AAslashField.hpp>
#include <Hadrons/Modules/MContraction/A2AFourQuarkContraction.hpp>
#include <Hadrons/Modules/MContraction/Baryon.hpp>
#include <Hadrons/Modules/MContraction/WeakEye3pt.hpp>
#include <Hadrons/Modules/MContraction/A2ALoop.hpp>
#include <Hadrons/Modules/MSink/Point.hpp>
#include <Hadrons/Modules/MSink/Smear.hpp>
#include <Hadrons/Modules/MScalarSUN/StochFreeField.hpp>
#include <Hadrons/Modules/MScalarSUN/EMT.hpp>
#include <Hadrons/Modules/MScalarSUN/Utils.hpp>
#include <Hadrons/Modules/MScalarSUN/TwoPoint.hpp>
#include <Hadrons/Modules/MScalarSUN/TransProj.hpp>
#include <Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp>
#include <Hadrons/Modules/MScalarSUN/TrPhi.hpp>
#include <Hadrons/Modules/MScalarSUN/Grad.hpp>
#include <Hadrons/Modules/MScalarSUN/TrMag.hpp>
#include <Hadrons/Modules/MScalarSUN/Div.hpp>
#include <Hadrons/Modules/MScalarSUN/TrKinetic.hpp>
#include <Hadrons/Modules/MNPR/Amputate.hpp>
#include <Hadrons/Modules/MNPR/FourQuark.hpp>
#include <Hadrons/Modules/MNPR/Bilinear.hpp>
#include <Hadrons/Modules/MSource/SeqAslash.hpp>
#include <Hadrons/Modules/MSource/Momentum.hpp>
#include <Hadrons/Modules/MSource/Z2.hpp>
#include <Hadrons/Modules/MSource/Point.hpp>
#include <Hadrons/Modules/MSource/Gauss.hpp>
#include <Hadrons/Modules/MSource/SeqConserved.hpp>
#include <Hadrons/Modules/MSource/Wall.hpp>
#include <Hadrons/Modules/MSource/SeqGamma.hpp>
#include <Hadrons/Modules/MSource/Convolution.hpp>
#include <Hadrons/Modules/MGauge/Random.hpp>
#include <Hadrons/Modules/MContraction/A2AMesonField.hpp>
#include <Hadrons/Modules/MContraction/Baryon.hpp>
#include <Hadrons/Modules/MContraction/DiscLoop.hpp>
#include <Hadrons/Modules/MContraction/Gamma3pt.hpp>
#include <Hadrons/Modules/MContraction/Meson.hpp>
#include <Hadrons/Modules/MContraction/SigmaToNucleonEye.hpp>
#include <Hadrons/Modules/MContraction/SigmaToNucleonNonEye.hpp>
#include <Hadrons/Modules/MContraction/WeakEye3pt.hpp>
#include <Hadrons/Modules/MContraction/WeakMesonDecayKl2.hpp>
#include <Hadrons/Modules/MContraction/WeakNonEye3pt.hpp>
#include <Hadrons/Modules/MDistil/Distil.hpp>
#include <Hadrons/Modules/MDistil/DistilPar.hpp>
#include <Hadrons/Modules/MDistil/DistilVectors.hpp>
#include <Hadrons/Modules/MDistil/LapEvec.hpp>
#include <Hadrons/Modules/MDistil/Noises.hpp>
#include <Hadrons/Modules/MDistil/PerambFromSolve.hpp>
#include <Hadrons/Modules/MDistil/Perambulator.hpp>
#include <Hadrons/Modules/MFermion/EMLepton.hpp>
#include <Hadrons/Modules/MFermion/FreeProp.hpp>
#include <Hadrons/Modules/MFermion/GaugeProp.hpp>
#include <Hadrons/Modules/MGauge/Electrify.hpp>
#include <Hadrons/Modules/MGauge/FundtoHirep.hpp>
#include <Hadrons/Modules/MGauge/StochEm.hpp>
#include <Hadrons/Modules/MGauge/UnitEm.hpp>
#include <Hadrons/Modules/MGauge/GaugeFix.hpp>
#include <Hadrons/Modules/MGauge/Random.hpp>
#include <Hadrons/Modules/MGauge/StochEm.hpp>
#include <Hadrons/Modules/MGauge/StoutSmearing.hpp>
#include <Hadrons/Modules/MGauge/Unit.hpp>
#include <Hadrons/Modules/MGauge/Electrify.hpp>
#include <Hadrons/Modules/MSolver/A2AVectors.hpp>
#include <Hadrons/Modules/MSolver/RBPrecCG.hpp>
#include <Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp>
#include <Hadrons/Modules/MSolver/Guesser.hpp>
#include <Hadrons/Modules/MSolver/MixedPrecisionRBPrecCG.hpp>
#include <Hadrons/Modules/MSolver/A2AAslashVectors.hpp>
#include <Hadrons/Modules/MGauge/UnitEm.hpp>
#include <Hadrons/Modules/MIO/LoadA2AMatrixDiskVector.hpp>
#include <Hadrons/Modules/MIO/LoadA2AVectors.hpp>
#include <Hadrons/Modules/MIO/LoadBinary.hpp>
#include <Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp>
#include <Hadrons/Modules/MIO/LoadCosmHol.hpp>
#include <Hadrons/Modules/MIO/LoadDistilNoise.hpp>
#include <Hadrons/Modules/MIO/LoadEigenPack.hpp>
#include <Hadrons/Modules/MIO/LoadNersc.hpp>
#include <Hadrons/Modules/MIO/LoadPerambulator.hpp>
#include <Hadrons/Modules/MNPR/Amputate.hpp>
#include <Hadrons/Modules/MNPR/Bilinear.hpp>
#include <Hadrons/Modules/MNPR/FourQuark.hpp>
#include <Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp>
#include <Hadrons/Modules/MNoise/SparseSpinColorDiagonal.hpp>
#include <Hadrons/Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp>
#include <Hadrons/Modules/MScalar/ChargedProp.hpp>
#include <Hadrons/Modules/MScalar/Scalar.hpp>
#include <Hadrons/Modules/MScalar/FreeProp.hpp>
#include <Hadrons/Modules/MScalar/Scalar.hpp>
#include <Hadrons/Modules/MScalarSUN/Div.hpp>
#include <Hadrons/Modules/MScalarSUN/EMT.hpp>
#include <Hadrons/Modules/MScalarSUN/Grad.hpp>
#include <Hadrons/Modules/MScalarSUN/StochFreeField.hpp>
#include <Hadrons/Modules/MScalarSUN/TrKinetic.hpp>
#include <Hadrons/Modules/MScalarSUN/TrMag.hpp>
#include <Hadrons/Modules/MScalarSUN/TrPhi.hpp>
#include <Hadrons/Modules/MScalarSUN/TransProj.hpp>
#include <Hadrons/Modules/MScalarSUN/TwoPoint.hpp>
#include <Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp>
#include <Hadrons/Modules/MScalarSUN/Utils.hpp>
#include <Hadrons/Modules/MSink/Point.hpp>
#include <Hadrons/Modules/MSink/Smear.hpp>
#include <Hadrons/Modules/MSolver/A2AAslashVectors.hpp>
#include <Hadrons/Modules/MSolver/A2AVectors.hpp>
#include <Hadrons/Modules/MSolver/Guesser.hpp>
#include <Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp>
#include <Hadrons/Modules/MSolver/MixedPrecisionRBPrecCG.hpp>
#include <Hadrons/Modules/MSolver/RBPrecCG.hpp>
#include <Hadrons/Modules/MSource/Convolution.hpp>
#include <Hadrons/Modules/MSource/Gauss.hpp>
#include <Hadrons/Modules/MSource/Momentum.hpp>
#include <Hadrons/Modules/MSource/MomentumPhase.hpp>
#include <Hadrons/Modules/MSource/Point.hpp>
#include <Hadrons/Modules/MSource/SeqAslash.hpp>
#include <Hadrons/Modules/MSource/SeqConserved.hpp>
#include <Hadrons/Modules/MSource/SeqGamma.hpp>
#include <Hadrons/Modules/MSource/Wall.hpp>
#include <Hadrons/Modules/MSource/Z2.hpp>
#include <Hadrons/Modules/MUtilities/PrecisionCast.hpp>
#include <Hadrons/Modules/MUtilities/RandomVectors.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)));
}

7
Hadrons/Modules/MContraction/SigmaToNucleonEye.cc Normal file
View File

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

218
Hadrons/Modules/MContraction/SigmaToNucleonEye.hpp Normal file
View File

@ -0,0 +1,218 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MContraction/SigmaToNucleonEye.hpp
Copyright (C) 2015-2019
Author: Antonin Portelli <antonin.portelli@me.com>
Author: Felix Erben <felix.erben@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_MContraction_SigmaToNucleonEye_hpp_
#define Hadrons_MContraction_SigmaToNucleonEye_hpp_
#include <Hadrons/Global.hpp>
#include <Hadrons/Module.hpp>
#include <Hadrons/ModuleFactory.hpp>
#include <Grid/qcd/utils/BaryonUtils.h>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* SigmaToNucleonEye *
******************************************************************************/
/*
* Sigma-to-nucleon 3-pt diagrams, eye topologies.
*
* Schematics: qqLoop |
* /->-¬ |
* / \ | qsTi G qdTf
* \ / | /---->------*------>----¬
* qsTi \ / qdTf | / /-*-¬ \
* /----->-----* *----->----¬ | / / G \ \
* * G G * | * \ / qqLoop *
* |\ /| | |\ \-<-/ /|
* | \ / | | | \ / |
* | \---------->---------/ | | | \----------->----------/ |
* \ quSpec / | \ quSpec /
* \ / | \ /
* \---------->---------/ | \----------->----------/
* quSpec | quSpec
*
* analogously to the rare-kaon naming, the left diagram is named 'one-trace' and
* the diagram on the right 'two-trace'
*
* Propagators:
* * qqLoop
* * quSpec, source at ti
* * qdTf, source at tf
* * qsTi, source at ti
*/
BEGIN_MODULE_NAMESPACE(MContraction)
class SigmaToNucleonEyePar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(SigmaToNucleonEyePar,
std::string, qqLoop,
std::string, quSpec,
std::string, qdTf,
std::string, qsTi,
unsigned int, tf,
std::string, sink,
std::string, output);
};
template <typename FImpl>
class TSigmaToNucleonEye: public Module<SigmaToNucleonEyePar>
{
public:
FERM_TYPE_ALIASES(FImpl,);
BASIC_TYPE_ALIASES(ScalarImplCR, Scalar);
SINK_TYPE_ALIASES(Scalar);
typedef typename SpinMatrixField::vector_object::scalar_object SpinMatrix;
class Metadata: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(Metadata,
Gamma::Algebra, gammaH,
Gamma::Algebra, gammaASigma,
Gamma::Algebra, gammaBSigma,
Gamma::Algebra, gammaANucl,
Gamma::Algebra, gammaBNucl,
int, trace);
};
typedef Correlator<Metadata, SpinMatrix> Result;
public:
// constructor
TSigmaToNucleonEye(const std::string name);
// destructor
virtual ~TSigmaToNucleonEye(void) {};
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
protected:
// setup
virtual void setup(void);
// execution
virtual void execute(void);
// Which gamma algebra was specified
Gamma::Algebra al;
};
MODULE_REGISTER_TMP(SigmaToNucleonEye, ARG(TSigmaToNucleonEye<FIMPL>), MContraction);
/******************************************************************************
* TSigmaToNucleonEye implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TSigmaToNucleonEye<FImpl>::TSigmaToNucleonEye(const std::string name)
: Module<SigmaToNucleonEyePar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TSigmaToNucleonEye<FImpl>::getInput(void)
{
std::vector<std::string> input = {par().qqLoop, par().quSpec, par().qdTf, par().qsTi, par().sink};
return input;
}
template <typename FImpl>
std::vector<std::string> TSigmaToNucleonEye<FImpl>::getOutput(void)
{
std::vector<std::string> out = {};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TSigmaToNucleonEye<FImpl>::setup(void)
{
envTmpLat(SpinMatrixField, "c");
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TSigmaToNucleonEye<FImpl>::execute(void)
{
const Gamma GammaB(Gamma::Algebra::SigmaXZ); // C*gamma_5
const Gamma Id(Gamma::Algebra::Identity); // C*gamma_5
LOG(Message) << "Computing sigma-to-nucleon contractions '" << getName() << "'" << std::endl;
LOG(Message) << "' with (Gamma^A,Gamma^B)_sigma = ( Identity, C*gamma_5 ) and (Gamma^A,Gamma^B)_nucl = ( Identity, C*gamma_5 )" << std::endl;
LOG(Message) << " using sink " << par().sink << "." << std::endl;
envGetTmp(SpinMatrixField, c);
std::vector<SpinMatrix> buf;
std::vector<Result> result;
Result r;
r.info.gammaASigma = Id.g;
r.info.gammaBSigma = GammaB.g;
r.info.gammaANucl = Id.g;
r.info.gammaBNucl = GammaB.g;
auto &qqLoop = envGet(PropagatorField, par().qqLoop);
auto &quSpec = envGet(SlicedPropagator, par().quSpec);
auto &qdTf = envGet(PropagatorField, par().qdTf);
auto &qsTi = envGet(PropagatorField, par().qsTi);
auto qut = quSpec[par().tf];
for (auto &G: Gamma::gall)
{
r.info.gammaH = G.g;
//Operator Q1, equivalent to the two-trace case in the rare-kaons module
c=Zero();
BaryonUtils<FIMPL>::Sigma_to_Nucleon_Eye(qqLoop,qut,qdTf,qsTi,G,GammaB,GammaB,"Q1",c);
sliceSum(c,buf,Tp);
r.corr.clear();
for (unsigned int t = 0; t < buf.size(); ++t)
{
r.corr.push_back(buf[t]);
}
r.info.trace = 2;
result.push_back(r);
//Operator Q2, equivalent to the one-trace case in the rare-kaons module
c=Zero();
BaryonUtils<FIMPL>::Sigma_to_Nucleon_Eye(qqLoop,qut,qdTf,qsTi,G,GammaB,GammaB,"Q2",c);
sliceSum(c,buf,Tp);
r.corr.clear();
for (unsigned int t = 0; t < buf.size(); ++t)
{
r.corr.push_back(buf[t]);
}
r.info.trace = 1;
result.push_back(r);
}
saveResult(par().output, "stnEye", result);
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MContraction_SigmaToNucleonEye_hpp_

7
Hadrons/Modules/MContraction/SigmaToNucleonNonEye.cc Normal file
View File

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

224
Hadrons/Modules/MContraction/SigmaToNucleonNonEye.hpp Normal file
View File

@ -0,0 +1,224 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MContraction/SigmaToNucleonNonEye.hpp
Copyright (C) 2015-2019
Author: Antonin Portelli <antonin.portelli@me.com>
Author: Felix Erben <felix.erben@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_MContraction_SigmaToNucleonNonEye_hpp_
#define Hadrons_MContraction_SigmaToNucleonNonEye_hpp_
#include <Hadrons/Global.hpp>
#include <Hadrons/Module.hpp>
#include <Hadrons/ModuleFactory.hpp>
#include <Grid/qcd/utils/BaryonUtils.h>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* SigmaToNucleonNonEye *
******************************************************************************/
/*
* Sigma-to-Nucleon 3-pt diagrams, non-eye topologies.
*
* Schematic:
* qsTi quTf | qsTi qdTf
* /-->--¬ /-->--¬ | /-->--¬ /-->--¬
* / \ / \ | / \ / \
* / \ / \ | / \ / \
* / \ / \ | / \ / \
* * * G * | * G * * G *
* |\ * G | | |\ / \ /|
* | \ / \ /| | | \ / \ / |
* | \ / \ / | | | \ / \ / |
* | \ / \ / | | | \-->--/ \-->--/ |
* \ \-->--/ \-->--/ / | \ quTi quTf /
* \ quTi qdTf / | \ /
* \ / | \ /
* \--------->----------/ | \--------->-----------/
* quSpec | quSpec
*
*
* analogously to the rare-kaon naming, the left diagram is named 'one-trace' and
* the diagram on the right 'two-trace'
*
* Propagators:
* * quTi, source at ti
* * quTf, source at tf
* * quSpec, source at ti
* * qdTf, source at tf
* * qsTi, source at ti
*/
BEGIN_MODULE_NAMESPACE(MContraction)
class SigmaToNucleonNonEyePar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(SigmaToNucleonNonEyePar,
std::string, quTi,
std::string, quTf,
std::string, quSpec,
std::string, qdTf,
std::string, qsTi,
unsigned int, tf,
std::string, sink,
std::string, output);
};
template <typename FImpl>
class TSigmaToNucleonNonEye: public Module<SigmaToNucleonNonEyePar>
{
public:
FERM_TYPE_ALIASES(FImpl,);
BASIC_TYPE_ALIASES(ScalarImplCR, Scalar);
SINK_TYPE_ALIASES(Scalar);
typedef typename SpinMatrixField::vector_object::scalar_object SpinMatrix;
class Metadata: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(Metadata,
Gamma::Algebra, gammaH,
Gamma::Algebra, gammaASigma,
Gamma::Algebra, gammaBSigma,
Gamma::Algebra, gammaANucl,
Gamma::Algebra, gammaBNucl,
int, trace);
};
typedef Correlator<Metadata, SpinMatrix> Result;
public:
// constructor
TSigmaToNucleonNonEye(const std::string name);
// destructor
virtual ~TSigmaToNucleonNonEye(void) {};
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
protected:
// setup
virtual void setup(void);
// execution
virtual void execute(void);
// Which gamma algebra was specified
Gamma::Algebra al;
};
MODULE_REGISTER_TMP(SigmaToNucleonNonEye, ARG(TSigmaToNucleonNonEye<FIMPL>), MContraction);
/******************************************************************************
* TSigmaToNucleonNonEye implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TSigmaToNucleonNonEye<FImpl>::TSigmaToNucleonNonEye(const std::string name)
: Module<SigmaToNucleonNonEyePar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TSigmaToNucleonNonEye<FImpl>::getInput(void)
{
std::vector<std::string> input = {par().quTi, par().quTf, par().quSpec, par().qdTf, par().qsTi, par().sink};
return input;
}
template <typename FImpl>
std::vector<std::string> TSigmaToNucleonNonEye<FImpl>::getOutput(void)
{
std::vector<std::string> out = {};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TSigmaToNucleonNonEye<FImpl>::setup(void)
{
envTmpLat(SpinMatrixField, "c");
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TSigmaToNucleonNonEye<FImpl>::execute(void)
{
const Gamma GammaB(Gamma::Algebra::SigmaXZ); // C*gamma_5
const Gamma Id(Gamma::Algebra::Identity); // C*gamma_5
LOG(Message) << "Computing sigma-to-nucleon contractions '" << getName() << "'" << std::endl;
LOG(Message) << "' with (Gamma^A,Gamma^B)_sigma = ( Identity, C*gamma_5 ) and (Gamma^A,Gamma^B)_nucl = ( Identity, C*gamma_5 )" << std::endl;
LOG(Message) << " using sink " << par().sink << "." << std::endl;
envGetTmp(SpinMatrixField, c);
std::vector<SpinMatrix> buf;
std::vector<Result> result;
Result r;
r.info.gammaASigma = Id.g;
r.info.gammaBSigma = GammaB.g;
r.info.gammaANucl = Id.g;
r.info.gammaBNucl = GammaB.g;
auto &quTi = envGet(PropagatorField, par().quTi);
auto &quTf = envGet(PropagatorField, par().quTf);
auto &quSpec = envGet(SlicedPropagator, par().quSpec);
auto &qdTf = envGet(PropagatorField, par().qdTf);
auto &qsTi = envGet(PropagatorField, par().qsTi);
auto qut = quSpec[par().tf];
for (auto &G: Gamma::gall)
{
r.info.gammaH = G.g;
//Operator Q1, equivalent to the two-trace case in the rare-kaons module
c=Zero();
BaryonUtils<FIMPL>::Sigma_to_Nucleon_NonEye(quTi,quTf,qut,qdTf,qsTi,G,GammaB,GammaB,"Q1",c);
sliceSum(c,buf,Tp);
r.corr.clear();
for (unsigned int t = 0; t < buf.size(); ++t)
{
r.corr.push_back(buf[t]);
}
r.info.trace = 2;
result.push_back(r);
//Operator Q2, equivalent to the one-trace case in the rare-kaons module
c=Zero();
BaryonUtils<FIMPL>::Sigma_to_Nucleon_NonEye(quTi,quTf,qut,qdTf,qsTi,G,GammaB,GammaB,"Q2",c);
sliceSum(c,buf,Tp);
r.corr.clear();
for (unsigned int t = 0; t < buf.size(); ++t)
{
r.corr.push_back(buf[t]);
}
r.info.trace = 1;
result.push_back(r);
}
saveResult(par().output, "stnNonEye", result);
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MContraction_SigmaToNucleonNonEye_hpp_

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;
}
}

8
Hadrons/Modules/MNoise/SparseSpinColorDiagonal.cc Normal file
View File

@ -0,0 +1,8 @@
#include <Hadrons/Modules/MNoise/SparseSpinColorDiagonal.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MNoise;
template class Grid::Hadrons::MNoise::TSparseSpinColorDiagonal<FIMPL>;
template class Grid::Hadrons::MNoise::TSparseSpinColorDiagonal<ZFIMPL>;

96
Hadrons/Modules/MNoise/SparseSpinColorDiagonal.hpp Normal file
View File

@ -0,0 +1,96 @@
#ifndef Hadrons_MNoise_SparseSpinColorDiagonal_hpp_
#define Hadrons_MNoise_SparseSpinColorDiagonal_hpp_
#include <Hadrons/Global.hpp>
#include <Hadrons/Module.hpp>
#include <Hadrons/ModuleFactory.hpp>
#include <Hadrons/DilutedNoise.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* SparseSpinColorDiagonal *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MNoise)
class SparseSpinColorDiagonalPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(SparseSpinColorDiagonalPar,
unsigned int, nsrc,
unsigned int, nsparse);
};
template <typename FImpl>
class TSparseSpinColorDiagonal: public Module<SparseSpinColorDiagonalPar>
{
public:
FERM_TYPE_ALIASES(FImpl,);
public:
// constructor
TSparseSpinColorDiagonal(const std::string name);
// destructor
virtual ~TSparseSpinColorDiagonal(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(SparseSpinColorDiagonal, TSparseSpinColorDiagonal<FIMPL>, MNoise);
MODULE_REGISTER_TMP(ZSparseSpinColorDiagonal, TSparseSpinColorDiagonal<ZFIMPL>, MNoise);
/******************************************************************************
* TSparseSpinColorDiagonal implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TSparseSpinColorDiagonal<FImpl>::TSparseSpinColorDiagonal(const std::string name)
: Module<SparseSpinColorDiagonalPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TSparseSpinColorDiagonal<FImpl>::getInput(void)
{
std::vector<std::string> in;
return in;
}
template <typename FImpl>
std::vector<std::string> TSparseSpinColorDiagonal<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TSparseSpinColorDiagonal<FImpl>::setup(void)
{
envCreateDerived(DilutedNoise<FImpl>,
SparseSpinColorDiagonalNoise<FImpl>,
getName(), 1, envGetGrid(FermionField), par().nsrc, par().nsparse);
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TSparseSpinColorDiagonal<FImpl>::execute(void)
{
auto &noise = envGet(DilutedNoise<FImpl>, getName());
LOG(Message) << "Generating sparse, spin-color diagonal noise with nSparse = "
<< par().nsparse << std::endl;
noise.generateNoise(rng4d());
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MNoise_SparseSpinColorDiagonal_hpp_

2
Hadrons/Modules/MSource/MomentumPhase.hpp
View File

@ -140,7 +140,7 @@ void TMomentumPhase<FImpl>::execute(void)
envGetTmp(LatticeComplex, coor);
p = strToVec<Real>(par().mom);
ph = zero;
ph = Zero();
for(unsigned int mu = 0; mu < env().getNd(); mu++)
{
LatticeCoordinate(coor, mu);

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_

6
Hadrons/make_module_list.sh
View File

@ -1,12 +1,12 @@
#!/usr/bin/env bash
echo 'modules_cc =\' > modules.inc
find Modules -name '*.cc' -type f -print | sed 's/^/ /;$q;s/$/ \\/' >> modules.inc
find Modules -name '*.cc' -type f -print | sort | sed 's/^/ /;$q;s/$/ \\/' >> modules.inc
echo '' >> modules.inc
echo 'modules_hpp =\' >> modules.inc
find Modules -name '*.hpp' -type f -print | sed 's/^/ /;$q;s/$/ \\/' >> modules.inc
find Modules -name '*.hpp' -type f -print | sort | sed 's/^/ /;$q;s/$/ \\/' >> modules.inc
echo '' >> modules.inc
rm -f Modules.hpp
for f in `find Modules -name '*.hpp'`; do
for f in `find Modules -name '*.hpp' | sort`; do
echo "#include <Hadrons/${f}>" >> Modules.hpp
done

261
Hadrons/modules.inc
View File

@ -1,147 +1,172 @@
modules_cc =\
Modules/MFermion/FreeProp.cc \
Modules/MFermion/GaugeProp.cc \
Modules/MFermion/EMLepton.cc \
Modules/MIO/LoadA2AVectors.cc \
Modules/MIO/LoadEigenPack.cc \
Modules/MIO/LoadCosmHol.cc \
Modules/MIO/LoadBinary.cc \
Modules/MIO/LoadA2AMatrixDiskVector.cc \
Modules/MIO/LoadCoarseEigenPack.cc \
Modules/MIO/LoadNersc.cc \
Modules/MAction/ZMobiusDWF.cc \
Modules/MAction/DWF.cc \
Modules/MAction/MobiusDWF.cc \
Modules/MAction/ScaledDWF.cc \
Modules/MAction/Wilson.cc \
Modules/MAction/DWF.cc \
Modules/MAction/WilsonClover.cc \
Modules/MAction/MobiusDWF.cc \
Modules/MUtilities/RandomVectors.cc \
Modules/MUtilities/PrecisionCast.cc \
Modules/MNoise/FullVolumeSpinColorDiagonal.cc \
Modules/MNoise/TimeDilutedSpinColorDiagonal.cc \
Modules/MContraction/Gamma3pt.cc \
Modules/MContraction/A2AFourQuarkContraction.cc \
Modules/MAction/ZMobiusDWF.cc \
Modules/MContraction/A2AAslashField.cc \
Modules/MContraction/A2AFourQuarkContraction.cc \
Modules/MContraction/A2ALoop.cc \
Modules/MContraction/WeakEye3pt.cc \
Modules/MContraction/WeakNonEye3pt.cc \
Modules/MContraction/A2AMesonField.cc \
Modules/MContraction/DiscLoop.cc \
Modules/MContraction/Baryon.cc \
Modules/MContraction/WeakMesonDecayKl2.cc \
Modules/MContraction/DiscLoop.cc \
Modules/MContraction/Gamma3pt.cc \
Modules/MContraction/Meson.cc \
Modules/MSink/Point.cc \
Modules/MSink/Smear.cc \
Modules/MScalarSUN/TrPhi.cc \
Modules/MScalarSUN/TrMag.cc \
Modules/MScalarSUN/TrKinetic.cc \
Modules/MScalarSUN/TwoPoint.cc \
Modules/MScalarSUN/Grad.cc \
Modules/MScalarSUN/TwoPointNPR.cc \
Modules/MScalarSUN/StochFreeField.cc \
Modules/MScalarSUN/TransProj.cc \
Modules/MScalarSUN/EMT.cc \
Modules/MScalarSUN/Div.cc \
Modules/MContraction/SigmaToNucleonEye.cc \
Modules/MContraction/SigmaToNucleonNonEye.cc \
Modules/MContraction/WeakEye3pt.cc \
Modules/MContraction/WeakMesonDecayKl2.cc \
Modules/MContraction/WeakNonEye3pt.cc \
Modules/MDistil/DistilPar.cc \
Modules/MDistil/DistilVectors.cc \
Modules/MDistil/LapEvec.cc \
Modules/MDistil/Noises.cc \
Modules/MDistil/PerambFromSolve.cc \
Modules/MDistil/Perambulator.cc \
Modules/MFermion/EMLepton.cc \
Modules/MFermion/FreeProp.cc \
Modules/MFermion/GaugeProp.cc \
Modules/MGauge/Electrify.cc \
Modules/MGauge/FundtoHirep.cc \
Modules/MGauge/GaugeFix.cc \
Modules/MGauge/Random.cc \
Modules/MGauge/StochEm.cc \
Modules/MGauge/StoutSmearing.cc \
Modules/MGauge/Unit.cc \
Modules/MGauge/UnitEm.cc \
Modules/MIO/LoadA2AMatrixDiskVector.cc \
Modules/MIO/LoadA2AVectors.cc \
Modules/MIO/LoadBinary.cc \
Modules/MIO/LoadCoarseEigenPack.cc \
Modules/MIO/LoadCosmHol.cc \
Modules/MIO/LoadDistilNoise.cc \
Modules/MIO/LoadEigenPack.cc \
Modules/MIO/LoadNersc.cc \
Modules/MIO/LoadPerambulator.cc \
Modules/MNPR/Amputate.cc \
Modules/MNPR/Bilinear.cc \
Modules/MNPR/FourQuark.cc \
Modules/MSource/SeqGamma.cc \
Modules/MSource/Z2.cc \
Modules/MSource/Convolution.cc \
Modules/MSource/Momentum.cc \
Modules/MSource/Wall.cc \
Modules/MSource/Point.cc \
Modules/MSource/SeqAslash.cc \
Modules/MSource/Gauss.cc \
Modules/MSource/SeqConserved.cc \
Modules/MGauge/StoutSmearing.cc \
Modules/MGauge/GaugeFix.cc \
Modules/MGauge/Electrify.cc \
Modules/MGauge/Random.cc \
Modules/MGauge/Unit.cc \
Modules/MGauge/UnitEm.cc \
Modules/MGauge/FundtoHirep.cc \
Modules/MGauge/StochEm.cc \
Modules/MSolver/RBPrecCG.cc \
Modules/MNoise/FullVolumeSpinColorDiagonal.cc \
Modules/MNoise/SparseSpinColorDiagonal.cc \
Modules/MNoise/TimeDilutedSpinColorDiagonal.cc \
Modules/MScalar/ChargedProp.cc \
Modules/MScalar/FreeProp.cc \
Modules/MScalarSUN/Div.cc \
Modules/MScalarSUN/EMT.cc \
Modules/MScalarSUN/Grad.cc \
Modules/MScalarSUN/StochFreeField.cc \
Modules/MScalarSUN/TrKinetic.cc \
Modules/MScalarSUN/TrMag.cc \
Modules/MScalarSUN/TrPhi.cc \
Modules/MScalarSUN/TransProj.cc \
Modules/MScalarSUN/TwoPoint.cc \
Modules/MScalarSUN/TwoPointNPR.cc \
Modules/MSink/Point.cc \
Modules/MSink/Smear.cc \
Modules/MSolver/A2AAslashVectors.cc \
Modules/MSolver/MixedPrecisionRBPrecCG.cc \
Modules/MSolver/A2AVectors.cc \
Modules/MSolver/LocalCoherenceLanczos.cc \
Modules/MScalar/ChargedProp.cc \
Modules/MScalar/FreeProp.cc
Modules/MSolver/MixedPrecisionRBPrecCG.cc \
Modules/MSolver/RBPrecCG.cc \
Modules/MSource/Convolution.cc \
Modules/MSource/Gauss.cc \
Modules/MSource/Momentum.cc \
Modules/MSource/MomentumPhase.cc \
Modules/MSource/Point.cc \
Modules/MSource/SeqAslash.cc \
Modules/MSource/SeqConserved.cc \
Modules/MSource/SeqGamma.cc \
Modules/MSource/Wall.cc \
Modules/MSource/Z2.cc \
Modules/MUtilities/PrecisionCast.cc \
Modules/MUtilities/RandomVectors.cc
modules_hpp =\
Modules/MFermion/GaugeProp.hpp \
Modules/MFermion/EMLepton.hpp \
Modules/MFermion/FreeProp.hpp \
Modules/MIO/LoadCosmHol.hpp \
Modules/MIO/LoadEigenPack.hpp \
Modules/MIO/LoadA2AVectors.hpp \
Modules/MIO/LoadA2AMatrixDiskVector.hpp \
Modules/MIO/LoadCoarseEigenPack.hpp \
Modules/MIO/LoadNersc.hpp \
Modules/MIO/LoadBinary.hpp \
Modules/MAction/ZMobiusDWF.hpp \
Modules/MAction/MobiusDWF.hpp \
Modules/MAction/Wilson.hpp \
Modules/MAction/DWF.hpp \
Modules/MAction/WilsonClover.hpp \
Modules/MAction/MobiusDWF.hpp \
Modules/MAction/ScaledDWF.hpp \
Modules/MUtilities/RandomVectors.hpp \
Modules/MUtilities/PrecisionCast.hpp \
Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp \
Modules/MNoise/FullVolumeSpinColorDiagonal.hpp \
Modules/MContraction/DiscLoop.hpp \
Modules/MContraction/Meson.hpp \
Modules/MContraction/WeakMesonDecayKl2.hpp \
Modules/MContraction/A2AMesonField.hpp \
Modules/MContraction/WeakNonEye3pt.hpp \
Modules/MContraction/Gamma3pt.hpp \
Modules/MAction/Wilson.hpp \
Modules/MAction/WilsonClover.hpp \
Modules/MAction/ZMobiusDWF.hpp \
Modules/MContraction/A2AAslashField.hpp \
Modules/MContraction/A2AFourQuarkContraction.hpp \
Modules/MContraction/Baryon.hpp \
Modules/MContraction/WeakEye3pt.hpp \
Modules/MContraction/A2ALoop.hpp \
Modules/MSink/Point.hpp \
Modules/MSink/Smear.hpp \
Modules/MScalarSUN/StochFreeField.hpp \
Modules/MScalarSUN/EMT.hpp \
Modules/MScalarSUN/Utils.hpp \
Modules/MScalarSUN/TwoPoint.hpp \
Modules/MScalarSUN/TransProj.hpp \
Modules/MScalarSUN/TwoPointNPR.hpp \
Modules/MScalarSUN/TrPhi.hpp \
Modules/MScalarSUN/Grad.hpp \
Modules/MScalarSUN/TrMag.hpp \
Modules/MScalarSUN/Div.hpp \
Modules/MScalarSUN/TrKinetic.hpp \
Modules/MNPR/Amputate.hpp \
Modules/MNPR/FourQuark.hpp \
Modules/MNPR/Bilinear.hpp \
Modules/MSource/SeqAslash.hpp \
Modules/MSource/Momentum.hpp \
Modules/MSource/Z2.hpp \
Modules/MSource/Point.hpp \
Modules/MSource/Gauss.hpp \
Modules/MSource/SeqConserved.hpp \
Modules/MSource/Wall.hpp \
Modules/MSource/SeqGamma.hpp \
Modules/MSource/Convolution.hpp \
Modules/MGauge/Random.hpp \
Modules/MContraction/A2AMesonField.hpp \
Modules/MContraction/Baryon.hpp \
Modules/MContraction/DiscLoop.hpp \
Modules/MContraction/Gamma3pt.hpp \
Modules/MContraction/Meson.hpp \
Modules/MContraction/SigmaToNucleonEye.hpp \
Modules/MContraction/SigmaToNucleonNonEye.hpp \
Modules/MContraction/WeakEye3pt.hpp \
Modules/MContraction/WeakMesonDecayKl2.hpp \
Modules/MContraction/WeakNonEye3pt.hpp \
Modules/MDistil/Distil.hpp \
Modules/MDistil/DistilPar.hpp \
Modules/MDistil/DistilVectors.hpp \
Modules/MDistil/LapEvec.hpp \
Modules/MDistil/Noises.hpp \
Modules/MDistil/PerambFromSolve.hpp \
Modules/MDistil/Perambulator.hpp \
Modules/MFermion/EMLepton.hpp \
Modules/MFermion/FreeProp.hpp \
Modules/MFermion/GaugeProp.hpp \
Modules/MGauge/Electrify.hpp \
Modules/MGauge/FundtoHirep.hpp \
Modules/MGauge/StochEm.hpp \
Modules/MGauge/UnitEm.hpp \
Modules/MGauge/GaugeFix.hpp \
Modules/MGauge/Random.hpp \
Modules/MGauge/StochEm.hpp \
Modules/MGauge/StoutSmearing.hpp \
Modules/MGauge/Unit.hpp \
Modules/MGauge/Electrify.hpp \
Modules/MSolver/A2AVectors.hpp \
Modules/MSolver/RBPrecCG.hpp \
Modules/MSolver/LocalCoherenceLanczos.hpp \
Modules/MSolver/Guesser.hpp \
Modules/MSolver/MixedPrecisionRBPrecCG.hpp \
Modules/MSolver/A2AAslashVectors.hpp \
Modules/MGauge/UnitEm.hpp \
Modules/MIO/LoadA2AMatrixDiskVector.hpp \
Modules/MIO/LoadA2AVectors.hpp \
Modules/MIO/LoadBinary.hpp \
Modules/MIO/LoadCoarseEigenPack.hpp \
Modules/MIO/LoadCosmHol.hpp \
Modules/MIO/LoadDistilNoise.hpp \
Modules/MIO/LoadEigenPack.hpp \
Modules/MIO/LoadNersc.hpp \
Modules/MIO/LoadPerambulator.hpp \
Modules/MNPR/Amputate.hpp \
Modules/MNPR/Bilinear.hpp \
Modules/MNPR/FourQuark.hpp \
Modules/MNoise/FullVolumeSpinColorDiagonal.hpp \
Modules/MNoise/SparseSpinColorDiagonal.hpp \
Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp \
Modules/MScalar/ChargedProp.hpp \
Modules/MScalar/FreeProp.hpp \
Modules/MScalar/Scalar.hpp \
Modules/MScalar/FreeProp.hpp
Modules/MScalarSUN/Div.hpp \
Modules/MScalarSUN/EMT.hpp \
Modules/MScalarSUN/Grad.hpp \
Modules/MScalarSUN/StochFreeField.hpp \
Modules/MScalarSUN/TrKinetic.hpp \
Modules/MScalarSUN/TrMag.hpp \
Modules/MScalarSUN/TrPhi.hpp \
Modules/MScalarSUN/TransProj.hpp \
Modules/MScalarSUN/TwoPoint.hpp \
Modules/MScalarSUN/TwoPointNPR.hpp \
Modules/MScalarSUN/Utils.hpp \
Modules/MSink/Point.hpp \
Modules/MSink/Smear.hpp \
Modules/MSolver/A2AAslashVectors.hpp \
Modules/MSolver/A2AVectors.hpp \
Modules/MSolver/Guesser.hpp \
Modules/MSolver/LocalCoherenceLanczos.hpp \
Modules/MSolver/MixedPrecisionRBPrecCG.hpp \
Modules/MSolver/RBPrecCG.hpp \
Modules/MSource/Convolution.hpp \
Modules/MSource/Gauss.hpp \
Modules/MSource/Momentum.hpp \
Modules/MSource/MomentumPhase.hpp \
Modules/MSource/Point.hpp \
Modules/MSource/SeqAslash.hpp \
Modules/MSource/SeqConserved.hpp \
Modules/MSource/SeqGamma.hpp \
Modules/MSource/Wall.hpp \
Modules/MSource/Z2.hpp \
Modules/MUtilities/PrecisionCast.hpp \
Modules/MUtilities/RandomVectors.hpp

2
bootstrap.sh
View File

@ -1,6 +1,6 @@
#!/usr/bin/env bash
EIGEN_URL='http://bitbucket.org/eigen/eigen/get/3.3.7.tar.bz2'
EIGEN_URL='https://gitlab.com/libeigen/eigen/-/archive/3.3.7/eigen-3.3.7.tar.bz2'
echo "-- deploying Eigen source..."
ARC=`basename ${EIGEN_URL}`

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;
}

154
tests/hadrons/Test_sigma_to_nucleon.cc Normal file
View File

@ -0,0 +1,154 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Tests/Hadrons/Test_hadrons_spectrum.cc
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Hadrons/Application.hpp>
#include <Hadrons/Modules.hpp>
using namespace Grid;
using namespace Hadrons;
int main(int argc, char *argv[])
{
// 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;
std::vector<std::string> flavour = {"l", "s", "c"};
std::vector<double> mass = {.01, .04, .2 };
// global parameters
Application::GlobalPar globalPar;
globalPar.trajCounter.start = 1500;
globalPar.trajCounter.end = 1520;
globalPar.trajCounter.step = 20;
globalPar.runId = "test";
application.setPar(globalPar);
// gauge field
application.createModule<MGauge::Unit>("gauge");
// sources
MSource::Point::Par ptPar;
ptPar.position = "0 0 0 0";
application.createModule<MSource::Point>("pt_0", ptPar);
ptPar.position = "0 0 0 4";
application.createModule<MSource::Point>("pt_4", ptPar);
// sink
MSink::Point::Par sinkPar;
sinkPar.mom = "0 0 0";
application.createModule<MSink::ScalarPoint>("sink", sinkPar);
application.createModule<MSink::Point>("sink_spec", sinkPar);
// set fermion boundary conditions to be periodic space, antiperiodic time.
std::string boundary = "1 1 1 -1";
std::string twist = "0. 0. 0. 0.";
for (unsigned int i = 0; i < flavour.size(); ++i)
{
// actions
MAction::DWF::Par actionPar;
actionPar.gauge = "gauge";
actionPar.Ls = 12;
actionPar.M5 = 1.8;
actionPar.mass = mass[i];
actionPar.boundary = boundary;
actionPar.twist = twist;
application.createModule<MAction::DWF>("DWF_" + flavour[i], actionPar);
// solvers
MSolver::RBPrecCG::Par solverPar;
solverPar.action = "DWF_" + flavour[i];
solverPar.residual = 1.0e-8;
solverPar.maxIteration = 10000;
application.createModule<MSolver::RBPrecCG>("CG_" + flavour[i],
solverPar);
}
// propagators
MFermion::GaugeProp::Par quarkPar;
quarkPar.solver = "CG_l";
quarkPar.source = "pt_0";
application.createModule<MFermion::GaugeProp>("Qpt_l_0", quarkPar);
quarkPar.source = "pt_4";
application.createModule<MFermion::GaugeProp>("Qpt_l_4", quarkPar);
quarkPar.solver = "CG_s";
quarkPar.source = "pt_0";
application.createModule<MFermion::GaugeProp>("Qpt_s_0", quarkPar);
//This should be a loop - how do I make this?
quarkPar.solver = "CG_c";
quarkPar.source = "pt_0";
application.createModule<MFermion::GaugeProp>("Qpt_c_loop", quarkPar);
quarkPar.solver = "CG_l";
quarkPar.source = "pt_0";
application.createModule<MFermion::GaugeProp>("Qpt_l_loop", quarkPar);
MSink::Smear::Par smearPar;
smearPar.q="Qpt_l_0";
smearPar.sink = "sink_spec";
application.createModule<MSink::Smear>("Qpt_u_spec",smearPar);
MContraction::SigmaToNucleonEye::Par EyePar;
EyePar.output = "SigmaToNucleon/Eye_u";
EyePar.qqLoop = "Qpt_l_loop";
EyePar.quSpec = "Qpt_u_spec";
EyePar.qdTf = "Qpt_l_4";
EyePar.qsTi = "Qpt_s_0";
EyePar.tf = 4;
EyePar.sink = "sink";
application.createModule<MContraction::SigmaToNucleonEye>("SigmaToNucleonEye_u", EyePar);
EyePar.output = "SigmaToNucleon/Eye_c";
EyePar.qqLoop = "Qpt_c_loop";
application.createModule<MContraction::SigmaToNucleonEye>("SigmaToNucleonEye_c", EyePar);
MContraction::SigmaToNucleonNonEye::Par NonEyePar;
NonEyePar.output = "SigmaToNucleon/NonEye";
NonEyePar.quTi = "Qpt_l_0";
NonEyePar.quTf = "Qpt_l_4";
NonEyePar.quSpec = "Qpt_u_spec";
NonEyePar.qdTf = "Qpt_l_4";
NonEyePar.qsTi = "Qpt_s_0";
NonEyePar.tf = 4;
NonEyePar.sink = "sink";
application.createModule<MContraction::SigmaToNucleonNonEye>("SigmaToNucleonNonEye", NonEyePar);
// execution
application.saveParameterFile("stn.xml");
application.run();
// epilogue
LOG(Message) << "Grid is finalizing now" << std::endl;
Grid_finalize();
return EXIT_SUCCESS;
}

229
tests/solver/Test_dwf_hdcr.cc
View File

@ -1,4 +1,4 @@
/*************************************************************************************
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
@ -111,11 +111,7 @@ public:
}
void operator()(const FineField &in, FineField & out) {
if ( params.domaindecompose ) {
operatorSAP(in,out);
} else {
operatorCheby(in,out);
}
operatorCheby(in,out);
}
////////////////////////////////////////////////////////////////////////
@ -232,66 +228,6 @@ public:
}
#endif
void SAP (const FineField & src,FineField & psi){
Lattice<iScalar<vInteger> > coor(src.Grid());
Lattice<iScalar<vInteger> > subset(src.Grid());
FineField r(src.Grid());
FineField zz(src.Grid()); zz=Zero();
FineField vec1(src.Grid());
FineField vec2(src.Grid());
const RealD block=params.domainsize;
subset=Zero();
for(int mu=0;mu<Nd;mu++){
LatticeCoordinate(coor,mu+1);
coor = div(coor,block);
subset = subset+coor;
}
subset = mod(subset,(Integer)2);
ShiftedMdagMLinearOperator<Matrix,FineField> fMdagMOp(_SmootherMatrix,0.0);
Chebyshev<FineField> Cheby (params.lo,params.hi,params.order,InverseApproximation);
RealD resid;
for(int i=0;i<params.steps;i++){
// Even domain residual
_FineOperator.Op(psi,vec1);// this is the G5 herm bit
r= src - vec1 ;
resid = norm2(r) /norm2(src);
std::cout << "SAP "<<i<<" resid "<<resid<<std::endl;
// Even domain solve
r= where(subset==(Integer)0,r,zz);
_SmootherOperator.AdjOp(r,vec1);
Cheby(fMdagMOp,vec1,vec2); // solves MdagM = g5 M g5M
psi = psi + vec2;
// Odd domain residual
_FineOperator.Op(psi,vec1);// this is the G5 herm bit
r= src - vec1 ;
r= where(subset==(Integer)1,r,zz);
resid = norm2(r) /norm2(src);
std::cout << "SAP "<<i<<" resid "<<resid<<std::endl;
// Odd domain solve
_SmootherOperator.AdjOp(r,vec1);
Cheby(fMdagMOp,vec1,vec2); // solves MdagM = g5 M g5M
psi = psi + vec2;
_FineOperator.Op(psi,vec1);// this is the G5 herm bit
r= src - vec1 ;
resid = norm2(r) /norm2(src);
std::cout << "SAP "<<i<<" resid "<<resid<<std::endl;
}
};
void SmootherTest (const FineField & in){
FineField vec1(in.Grid());
@ -308,6 +244,8 @@ public:
for(int ilo=0;ilo<3;ilo++){
for(int ord=5;ord<50;ord*=2){
std::cout << " lo "<<lo[ilo]<<" order "<<ord<<std::endl;
_SmootherOperator.AdjOp(in,vec1);
Chebyshev<FineField> Cheby (lo[ilo],70.0,ord,InverseApproximation);
@ -329,7 +267,6 @@ public:
CoarseVector Csol(_CoarseOperator.Grid()); Csol=Zero();
ConjugateGradient<CoarseVector> CG(3.0e-3,100000);
// ConjugateGradient<FineField> fCG(3.0e-2,1000);
HermitianLinearOperator<CoarseOperator,CoarseVector> HermOp(_CoarseOperator);
MdagMLinearOperator<CoarseOperator,CoarseVector> MdagMOp(_CoarseOperator);
@ -339,14 +276,11 @@ public:
FineField vec1(in.Grid());
FineField vec2(in.Grid());
// Chebyshev<FineField> Cheby (0.5,70.0,30,InverseApproximation);
// Chebyshev<FineField> ChebyAccu(0.5,70.0,30,InverseApproximation);
Chebyshev<FineField> Cheby (params.lo,params.hi,params.order,InverseApproximation);
Chebyshev<FineField> ChebyAccu(params.lo,params.hi,params.order,InverseApproximation);
// Cheby.JacksonSmooth();
// ChebyAccu.JacksonSmooth();
// _Aggregates.ProjectToSubspace (Csrc,in);
_Aggregates.ProjectToSubspace (Csrc,in);
// _Aggregates.PromoteFromSubspace(Csrc,out);
// std::cout<<GridLogMessage<<"Completeness: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
@ -366,8 +300,6 @@ public:
_SmootherOperator.AdjOp(in,vec1);// this is the G5 herm bit
ChebyAccu(fMdagMOp,vec1,out); // solves MdagM = g5 M g5M
std::cout<<GridLogMessage << "Smoother norm "<<norm2(out)<<std::endl;
// Update with residual for out
_FineOperator.Op(out,vec1);// this is the G5 herm bit
vec1 = in - vec1; // tmp = in - A Min
@ -377,10 +309,11 @@ public:
std::cout<<GridLogMessage << "Smoother resid "<<std::sqrt(r/Ni)<< " " << r << " " << Ni <<std::endl;
_Aggregates.ProjectToSubspace (Csrc,vec1);
HermOp.AdjOp(Csrc,Ctmp);// Normal equations
HermOp.AdjOp(Csrc,Ctmp);// Normal equations // This appears to be zero.
CG(MdagMOp,Ctmp,Csol);
_Aggregates.PromoteFromSubspace(Csol,vec1); // Ass^{-1} [in - A Min]_s
// Q = Q[in - A Min]
// Q = Q[in - A Min]
out = out+vec1;
// Three preconditioner smoothing -- hermitian if C3 = C1
@ -402,69 +335,6 @@ public:
}
void operatorSAP(const FineField &in, FineField & out) {
CoarseVector Csrc(_CoarseOperator.Grid());
CoarseVector Ctmp(_CoarseOperator.Grid());
CoarseVector Csol(_CoarseOperator.Grid()); Csol=Zero();
ConjugateGradient<CoarseVector> CG(1.0e-3,100000);
HermitianLinearOperator<CoarseOperator,CoarseVector> HermOp(_CoarseOperator);
MdagMLinearOperator<CoarseOperator,CoarseVector> MdagMOp(_CoarseOperator);
FineField vec1(in.Grid());
FineField vec2(in.Grid());
_Aggregates.ProjectToSubspace (Csrc,in);
_Aggregates.PromoteFromSubspace(Csrc,out);
std::cout<<GridLogMessage<<"Completeness: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
// To make a working smoother for indefinite operator
// must multiply by "Mdag" (ouch loses all low mode content)
// and apply to poly approx of (mdagm)^-1.
// so that we end up with an odd polynomial.
SAP(in,out);
// Update with residual for out
_FineOperator.Op(out,vec1);// this is the G5 herm bit
vec1 = in - vec1; // tmp = in - A Min
RealD r = norm2(vec1);
RealD Ni = norm2(in);
std::cout<<GridLogMessage << "SAP resid "<<std::sqrt(r/Ni)<< " " << r << " " << Ni <<std::endl;
_Aggregates.ProjectToSubspace (Csrc,vec1);
HermOp.AdjOp(Csrc,Ctmp);// Normal equations
CG(MdagMOp,Ctmp,Csol);
_Aggregates.PromoteFromSubspace(Csol,vec1); // Ass^{-1} [in - A Min]_s
// Q = Q[in - A Min]
out = out+vec1;
// Three preconditioner smoothing -- hermitian if C3 = C1
// Recompute error
_FineOperator.Op(out,vec1);// this is the G5 herm bit
vec1 = in - vec1; // tmp = in - A Min
r=norm2(vec1);
std::cout<<GridLogMessage << "Coarse resid "<<std::sqrt(r/Ni)<<std::endl;
// Reapply smoother
SAP(vec1,vec2);
out =out+vec2;
// Update with residual for out
_FineOperator.Op(out,vec1);// this is the G5 herm bit
vec1 = in - vec1; // tmp = in - A Min
r = norm2(vec1);
Ni = norm2(in);
std::cout<<GridLogMessage << "SAP resid(post) "<<std::sqrt(r/Ni)<< " " << r << " " << Ni <<std::endl;
}
};
int main (int argc, char ** argv)
@ -559,8 +429,9 @@ int main (int argc, char ** argv)
assert ( (nbasis & 0x1)==0);
int nb=nbasis/2;
std::cout<<GridLogMessage << " nbasis/2 = "<<nb<<std::endl;
Aggregates.CreateSubspace(RNG5,HermDefOp,nb);
// Aggregates.CreateSubspace(RNG5,HermDefOp,nb);
// Aggregates.CreateSubspaceLanczos(RNG5,HermDefOp,nb);
Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb);
for(int n=0;n<nb;n++){
G5R5(Aggregates.subspace[n+nb],Aggregates.subspace[n]);
std::cout<<GridLogMessage<<n<<" subspace "<<norm2(Aggregates.subspace[n+nb])<<" "<<norm2(Aggregates.subspace[n]) <<std::endl;
@ -580,7 +451,7 @@ int main (int argc, char ** argv)
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
Gamma5R5HermitianLinearOperator<DomainWallFermionR,LatticeFermion> HermIndefOp(Ddwf);
Gamma5R5HermitianLinearOperator<DomainWallFermionR,LatticeFermion> HermIndefOpDD(DdwfDD);
CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LDOp(*Coarse5d);
CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LDOp(*Coarse5d,1); // Hermitian matrix
LDOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
@ -596,14 +467,14 @@ int main (int argc, char ** argv)
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
MdagMLinearOperator<CoarseOperator,CoarseVector> PosdefLdop(LDOp);
ConjugateGradient<CoarseVector> CG(1.0e-6,100000);
// CG(PosdefLdop,c_src,c_res);
CG(PosdefLdop,c_src,c_res);
// std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// std::cout<<GridLogMessage << "Solving indef-MCR on coarse space "<< std::endl;
// std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// HermitianLinearOperator<CoarseOperator,CoarseVector> HermIndefLdop(LDOp);
// ConjugateResidual<CoarseVector> MCR(1.0e-6,100000);
//MCR(HermIndefLdop,c_src,c_res);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Solving indef-MCR on coarse space "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// HermitianLinearOperator<CoarseOperator,CoarseVector> HermIndefLdop(LDOp);
// ConjugateResidual<CoarseVector> MCR(1.0e-6,100000);
// MCR(HermIndefLdop,c_src,c_res);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Building deflation preconditioner "<< std::endl;
@ -613,38 +484,48 @@ int main (int argc, char ** argv)
HermIndefOp,Ddwf,
HermIndefOp,Ddwf);
MultiGridPreconditioner <vSpinColourVector,vTComplex,nbasis,DomainWallFermionR> PreconDD(Aggregates, LDOp,
HermIndefOp,Ddwf,
HermIndefOpDD,DdwfDD);
TrivialPrecon<LatticeFermion> simple;
// MultiGridPreconditioner <vSpinColourVector,vTComplex,nbasis,DomainWallFermionR> PreconDD(Aggregates, LDOp,
// HermIndefOp,Ddwf,
// HermIndefOpDD,DdwfDD);
// TrivialPrecon<LatticeFermion> simple;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Testing smoother efficacy"<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// Precon.SmootherTest(src);
Precon.SmootherTest(src);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Testing DD smoother efficacy"<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// std::cout<<GridLogMessage << "Testing DD smoother efficacy"<< std::endl;
// std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// PreconDD.SmootherTest(src);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Testing SAP smoother efficacy"<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// std::cout<<GridLogMessage << "Testing SAP smoother efficacy"<< std::endl;
// std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// PreconDD.SAP(src,result);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Unprec CG "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// std::cout<<GridLogMessage << "Unprec CG "<< std::endl;
// std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// TrivialPrecon<LatticeFermion> simple;
// ConjugateGradient<LatticeFermion> fCG(1.0e-8,100000);
// fCG(HermDefOp,src,result);
// exit(0);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Testing GCR on indef matrix "<< std::endl;
std::cout<<GridLogMessage << "Red Black Prec CG "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
LatticeFermion src_o(FrbGrid);
LatticeFermion result_o(FrbGrid);
pickCheckerboard(Odd,src_o,src);
result_o=Zero();
SchurDiagMooeeOperator<DomainWallFermionR,LatticeFermion> HermOpEO(Ddwf);
ConjugateGradient<LatticeFermion> pCG(1.0e-8,10000);
// pCG(HermOpEO,src_o,result_o);
// std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// std::cout<<GridLogMessage << "Testing GCR on indef matrix "<< std::endl;
// std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// PrecGeneralisedConjugateResidual<LatticeFermion> UPGCR(1.0e-8,100000,simple,8,128);
// UPGCR(HermIndefOp,src,result);
@ -656,9 +537,9 @@ int main (int argc, char ** argv)
Precon.PowerMethod(src);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Building a two level DDPGCR "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// std::cout<<GridLogMessage << "Building a two level DDPGCR "<< std::endl;
// std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// PrecGeneralisedConjugateResidual<LatticeFermion> PGCRDD(1.0e-8,100000,PreconDD,8,128);
// result=Zero();
// std::cout<<GridLogMessage<<"checking norm src "<<norm2(src)<<std::endl;
@ -672,20 +553,6 @@ int main (int argc, char ** argv)
result=Zero();
PGCR(HermIndefOp,src,result);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Red Black Prec CG "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
SchurDiagMooeeOperator<DomainWallFermionR,LatticeFermion> HermOpEO(Ddwf);
ConjugateGradient<LatticeFermion> pCG(1.0e-8,10000);
LatticeFermion src_o(FrbGrid);
LatticeFermion result_o(FrbGrid);
pickCheckerboard(Odd,src_o,src);
result_o=Zero();
pCG(HermOpEO,src_o,result_o);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Done "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;

84
tests/solver/Test_dwf_qmr_unprec.cc Normal file
View File

@ -0,0 +1,84 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_dwf_cr_unprec.cc
Copyright (C) 2019
Author: Peter Boyle <pboyle@bnl.gov>
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 <Grid/Grid.h>
#include <Grid/algorithms/iterative/QuasiMinimalResidual.h>
using namespace std;
using namespace Grid;
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
const int Ls=8;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
LatticeFermion src(FGrid); random(RNG5,src);
LatticeFermion result(FGrid); result=Zero();
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
for(int mu=0;mu<Nd;mu++){
U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
}
RealD tol = 1.0e-8;
RealD maxit=2000;
QuasiMinimalResidual<LatticeFermion> QMR(tol,maxit);
GeneralisedMinimalResidual<LatticeFermion> GMR(tol, maxit, 32,false);
RealD mass=0.0;
RealD M5=-1.8;
DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
Gamma5R5HermitianLinearOperator<DomainWallFermionR,LatticeFermion> g5HermOp(Ddwf);
QMR(g5HermOp,src,result);
GMR(g5HermOp,src,result);
NonHermitianLinearOperator<DomainWallFermionR,LatticeFermion> NonHermOp(Ddwf);
QMR(NonHermOp,src,result);
GMR(NonHermOp,src,result);
MdagMLinearOperator<DomainWallFermionR,LatticeFermion> HermOp(Ddwf);
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
CG(HermOp,src,result);
Grid_finalize();
}

66
tests/solver/Test_wilson_qmr_unprec.cc Normal file
View File

@ -0,0 +1,66 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_dwf_cr_unprec.cc
Copyright (C) 2019
Author: Peter Boyle <pboyle@bnl.gov>
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 <Grid/Grid.h>
#include <Grid/algorithms/iterative/QuasiMinimalResidual.h>
using namespace std;
using namespace Grid;
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
GridCartesian * Grid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
GridRedBlackCartesian * rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(Grid);
std::vector<int> seeds4({1,2,3,4});
GridParallelRNG RNG4(Grid); RNG4.SeedFixedIntegers(seeds4);
LatticeFermion src(Grid); random(RNG4,src);
LatticeFermion result(Grid); result=Zero();
LatticeGaugeField Umu(Grid); SU3::HotConfiguration(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,Grid);
for(int mu=0;mu<Nd;mu++){
U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
}
QuasiMinimalResidual<LatticeFermion> QMR(1.0e-8,10000);
RealD mass=0.0;
RealD M5=1.8;
WilsonFermionR Dw(Umu,*Grid,*rbGrid,mass);
NonHermitianLinearOperator<WilsonFermionR,LatticeFermion> NonHermOp(Dw);
QMR(NonHermOp,src,result);
Grid_finalize();
}