1
0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-14 13:57:07 +01:00

Compare commits

..

2 Commits

Author SHA1 Message Date
251b904a28 Merge branch 'release/ISC-freeze-2' 2018-06-04 21:09:48 +01:00
5a112feac3 Merge branch 'release/ISC-freeze-1' 2018-06-04 18:49:40 +01:00
196 changed files with 1420 additions and 12238 deletions

View File

@ -1,803 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./benchmarks/Benchmark_wilson.cc
Copyright (C) 2018
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: paboyle <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 */
#include <Grid/Grid.h>
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
#include "Grid/util/Profiling.h"
template<class vobj>
void sliceInnerProductMesonField(std::vector< std::vector<ComplexD> > &mat,
const std::vector<Lattice<vobj> > &lhs,
const std::vector<Lattice<vobj> > &rhs,
int orthogdim)
{
typedef typename vobj::scalar_object sobj;
typedef typename vobj::scalar_type scalar_type;
typedef typename vobj::vector_type vector_type;
int Lblock = lhs.size();
int Rblock = rhs.size();
GridBase *grid = lhs[0]._grid;
const int Nd = grid->_ndimension;
const int Nsimd = grid->Nsimd();
int Nt = grid->GlobalDimensions()[orthogdim];
assert(mat.size()==Lblock*Rblock);
for(int t=0;t<mat.size();t++){
assert(mat[t].size()==Nt);
}
int fd=grid->_fdimensions[orthogdim];
int ld=grid->_ldimensions[orthogdim];
int rd=grid->_rdimensions[orthogdim];
// will locally sum vectors first
// sum across these down to scalars
// splitting the SIMD
std::vector<vector_type,alignedAllocator<vector_type> > lvSum(rd*Lblock*Rblock);
parallel_for (int r = 0; r < rd * Lblock * Rblock; r++){
lvSum[r] = zero;
}
std::vector<scalar_type > lsSum(ld*Lblock*Rblock,scalar_type(0.0));
int e1= grid->_slice_nblock[orthogdim];
int e2= grid->_slice_block [orthogdim];
int stride=grid->_slice_stride[orthogdim];
std::cout << GridLogMessage << " Entering first parallel loop "<<std::endl;
// Parallelise over t-direction doesn't expose as much parallelism as needed for KNL
parallel_for(int r=0;r<rd;r++){
int so=r*grid->_ostride[orthogdim]; // base offset for start of plane
for(int n=0;n<e1;n++){
for(int b=0;b<e2;b++){
int ss= so+n*stride+b;
for(int i=0;i<Lblock;i++){
auto left = conjugate(lhs[i]._odata[ss]);
for(int j=0;j<Rblock;j++){
int idx = i+Lblock*j+Lblock*Rblock*r;
auto right = rhs[j]._odata[ss];
vector_type vv = left()(0)(0) * right()(0)(0)
+ left()(0)(1) * right()(0)(1)
+ left()(0)(2) * right()(0)(2)
+ left()(1)(0) * right()(1)(0)
+ left()(1)(1) * right()(1)(1)
+ left()(1)(2) * right()(1)(2)
+ left()(2)(0) * right()(2)(0)
+ left()(2)(1) * right()(2)(1)
+ left()(2)(2) * right()(2)(2)
+ left()(3)(0) * right()(3)(0)
+ left()(3)(1) * right()(3)(1)
+ left()(3)(2) * right()(3)(2);
lvSum[idx]=lvSum[idx]+vv;
}
}
}
}
}
std::cout << GridLogMessage << " Entering second parallel loop "<<std::endl;
// Sum across simd lanes in the plane, breaking out orthog dir.
parallel_for(int rt=0;rt<rd;rt++){
std::vector<int> icoor(Nd);
for(int i=0;i<Lblock;i++){
for(int j=0;j<Rblock;j++){
iScalar<vector_type> temp;
std::vector<iScalar<scalar_type> > extracted(Nsimd);
temp._internal = lvSum[i+Lblock*j+Lblock*Rblock*rt];
extract(temp,extracted);
for(int idx=0;idx<Nsimd;idx++){
grid->iCoorFromIindex(icoor,idx);
int ldx =rt+icoor[orthogdim]*rd;
int ij_dx = i+Lblock*j+Lblock*Rblock*ldx;
lsSum[ij_dx]=lsSum[ij_dx]+extracted[idx]._internal;
}
}}
}
std::cout << GridLogMessage << " Entering non parallel loop "<<std::endl;
for(int t=0;t<fd;t++)
{
int pt = t / ld; // processor plane
int lt = t % ld;
for(int i=0;i<Lblock;i++){
for(int j=0;j<Rblock;j++){
if (pt == grid->_processor_coor[orthogdim]){
int ij_dx = i + Lblock * j + Lblock * Rblock * lt;
mat[i+j*Lblock][t] = lsSum[ij_dx];
}
else{
mat[i+j*Lblock][t] = scalar_type(0.0);
}
}}
}
std::cout << GridLogMessage << " Done "<<std::endl;
// defer sum over nodes.
return;
}
template<class vobj>
void sliceInnerProductMesonFieldGamma(std::vector< std::vector<ComplexD> > &mat,
const std::vector<Lattice<vobj> > &lhs,
const std::vector<Lattice<vobj> > &rhs,
int orthogdim,
std::vector<Gamma::Algebra> gammas)
{
typedef typename vobj::scalar_object sobj;
typedef typename vobj::scalar_type scalar_type;
typedef typename vobj::vector_type vector_type;
int Lblock = lhs.size();
int Rblock = rhs.size();
GridBase *grid = lhs[0]._grid;
const int Nd = grid->_ndimension;
const int Nsimd = grid->Nsimd();
int Nt = grid->GlobalDimensions()[orthogdim];
int Ngamma = gammas.size();
assert(mat.size()==Lblock*Rblock*Ngamma);
for(int t=0;t<mat.size();t++){
assert(mat[t].size()==Nt);
}
int fd=grid->_fdimensions[orthogdim];
int ld=grid->_ldimensions[orthogdim];
int rd=grid->_rdimensions[orthogdim];
// will locally sum vectors first
// sum across these down to scalars
// splitting the SIMD
int MFrvol = rd*Lblock*Rblock*Ngamma;
int MFlvol = ld*Lblock*Rblock*Ngamma;
std::vector<vector_type,alignedAllocator<vector_type> > lvSum(MFrvol);
parallel_for (int r = 0; r < MFrvol; r++){
lvSum[r] = zero;
}
std::vector<scalar_type > lsSum(MFlvol);
parallel_for (int r = 0; r < MFlvol; r++){
lsSum[r]=scalar_type(0.0);
}
int e1= grid->_slice_nblock[orthogdim];
int e2= grid->_slice_block [orthogdim];
int stride=grid->_slice_stride[orthogdim];
std::cout << GridLogMessage << " Entering first parallel loop "<<std::endl;
// Parallelise over t-direction doesn't expose as much parallelism as needed for KNL
parallel_for(int r=0;r<rd;r++){
int so=r*grid->_ostride[orthogdim]; // base offset for start of plane
for(int n=0;n<e1;n++){
for(int b=0;b<e2;b++){
int ss= so+n*stride+b;
for(int i=0;i<Lblock;i++){
auto left = conjugate(lhs[i]._odata[ss]);
for(int j=0;j<Rblock;j++){
for(int mu=0;mu<Ngamma;mu++){
auto right = Gamma(gammas[mu])*rhs[j]._odata[ss];
vector_type vv = left()(0)(0) * right()(0)(0)
+ left()(0)(1) * right()(0)(1)
+ left()(0)(2) * right()(0)(2)
+ left()(1)(0) * right()(1)(0)
+ left()(1)(1) * right()(1)(1)
+ left()(1)(2) * right()(1)(2)
+ left()(2)(0) * right()(2)(0)
+ left()(2)(1) * right()(2)(1)
+ left()(2)(2) * right()(2)(2)
+ left()(3)(0) * right()(3)(0)
+ left()(3)(1) * right()(3)(1)
+ left()(3)(2) * right()(3)(2);
int idx = mu+i*Ngamma+Lblock*Ngamma*j+Ngamma*Lblock*Rblock*r;
lvSum[idx]=lvSum[idx]+vv;
}
}
}
}
}
}
std::cout << GridLogMessage << " Entering second parallel loop "<<std::endl;
// Sum across simd lanes in the plane, breaking out orthog dir.
parallel_for(int rt=0;rt<rd;rt++){
iScalar<vector_type> temp;
std::vector<int> icoor(Nd);
std::vector<iScalar<scalar_type> > extracted(Nsimd);
for(int i=0;i<Lblock;i++){
for(int j=0;j<Rblock;j++){
for(int mu=0;mu<Ngamma;mu++){
int ij_rdx = mu+i*Ngamma+Ngamma*Lblock*j+Ngamma*Lblock*Rblock*rt;
temp._internal = lvSum[ij_rdx];
extract(temp,extracted);
for(int idx=0;idx<Nsimd;idx++){
grid->iCoorFromIindex(icoor,idx);
int ldx =rt+icoor[orthogdim]*rd;
int ij_ldx = mu+i*Ngamma+Ngamma*Lblock*j+Ngamma*Lblock*Rblock*ldx;
lsSum[ij_ldx]=lsSum[ij_ldx]+extracted[idx]._internal;
}
}}}
}
std::cout << GridLogMessage << " Entering non parallel loop "<<std::endl;
for(int t=0;t<fd;t++)
{
int pt = t / ld; // processor plane
int lt = t % ld;
for(int i=0;i<Lblock;i++){
for(int j=0;j<Rblock;j++){
for(int mu=0;mu<Ngamma;mu++){
if (pt == grid->_processor_coor[orthogdim]){
int ij_dx = mu+i*Ngamma+Ngamma*Lblock*j+Ngamma*Lblock*Rblock* lt;
mat[mu+i*Ngamma+j*Lblock*Ngamma][t] = lsSum[ij_dx];
}
else{
mat[mu+i*Ngamma+j*Lblock*Ngamma][t] = scalar_type(0.0);
}
}}}
}
std::cout << GridLogMessage << " Done "<<std::endl;
// defer sum over nodes.
return;
}
template<class vobj>
void sliceInnerProductMesonFieldGamma1(std::vector< std::vector<ComplexD> > &mat,
const std::vector<Lattice<vobj> > &lhs,
const std::vector<Lattice<vobj> > &rhs,
int orthogdim,
std::vector<Gamma::Algebra> gammas)
{
typedef typename vobj::scalar_object sobj;
typedef typename vobj::scalar_type scalar_type;
typedef typename vobj::vector_type vector_type;
typedef iSpinMatrix<vector_type> SpinMatrix_v;
typedef iSpinMatrix<scalar_type> SpinMatrix_s;
int Lblock = lhs.size();
int Rblock = rhs.size();
GridBase *grid = lhs[0]._grid;
const int Nd = grid->_ndimension;
const int Nsimd = grid->Nsimd();
int Nt = grid->GlobalDimensions()[orthogdim];
int Ngamma = gammas.size();
assert(mat.size()==Lblock*Rblock*Ngamma);
for(int t=0;t<mat.size();t++){
assert(mat[t].size()==Nt);
}
int fd=grid->_fdimensions[orthogdim];
int ld=grid->_ldimensions[orthogdim];
int rd=grid->_rdimensions[orthogdim];
// will locally sum vectors first
// sum across these down to scalars
// splitting the SIMD
int MFrvol = rd*Lblock*Rblock;
int MFlvol = ld*Lblock*Rblock;
Vector<SpinMatrix_v > lvSum(MFrvol);
parallel_for (int r = 0; r < MFrvol; r++){
lvSum[r] = zero;
}
Vector<SpinMatrix_s > lsSum(MFlvol);
parallel_for (int r = 0; r < MFlvol; r++){
lsSum[r]=scalar_type(0.0);
}
int e1= grid->_slice_nblock[orthogdim];
int e2= grid->_slice_block [orthogdim];
int stride=grid->_slice_stride[orthogdim];
std::cout << GridLogMessage << " Entering first parallel loop "<<std::endl;
// Parallelise over t-direction doesn't expose as much parallelism as needed for KNL
parallel_for(int r=0;r<rd;r++){
int so=r*grid->_ostride[orthogdim]; // base offset for start of plane
for(int n=0;n<e1;n++){
for(int b=0;b<e2;b++){
int ss= so+n*stride+b;
for(int i=0;i<Lblock;i++){
auto left = conjugate(lhs[i]._odata[ss]);
for(int j=0;j<Rblock;j++){
SpinMatrix_v vv;
auto right = rhs[j]._odata[ss];
for(int s1=0;s1<Ns;s1++){
for(int s2=0;s2<Ns;s2++){
vv()(s2,s1)() = left()(s1)(0) * right()(s2)(0)
+ left()(s1)(1) * right()(s2)(1)
+ left()(s1)(2) * right()(s2)(2);
}}
int idx = i+Lblock*j+Lblock*Rblock*r;
lvSum[idx]=lvSum[idx]+vv;
}
}
}
}
}
std::cout << GridLogMessage << " Entering second parallel loop "<<std::endl;
// Sum across simd lanes in the plane, breaking out orthog dir.
parallel_for(int rt=0;rt<rd;rt++){
std::vector<int> icoor(Nd);
std::vector<SpinMatrix_s> extracted(Nsimd);
for(int i=0;i<Lblock;i++){
for(int j=0;j<Rblock;j++){
int ij_rdx = i+Lblock*j+Lblock*Rblock*rt;
extract(lvSum[ij_rdx],extracted);
for(int idx=0;idx<Nsimd;idx++){
grid->iCoorFromIindex(icoor,idx);
int ldx = rt+icoor[orthogdim]*rd;
int ij_ldx = i+Lblock*j+Lblock*Rblock*ldx;
lsSum[ij_ldx]=lsSum[ij_ldx]+extracted[idx];
}
}}
}
std::cout << GridLogMessage << " Entering third parallel loop "<<std::endl;
parallel_for(int t=0;t<fd;t++)
{
int pt = t / ld; // processor plane
int lt = t % ld;
for(int i=0;i<Lblock;i++){
for(int j=0;j<Rblock;j++){
if (pt == grid->_processor_coor[orthogdim]){
int ij_dx = i + Lblock * j + Lblock * Rblock * lt;
for(int mu=0;mu<Ngamma;mu++){
mat[mu+i*Ngamma+j*Lblock*Ngamma][t] = trace(lsSum[ij_dx]*Gamma(gammas[mu]));
}
}
else{
for(int mu=0;mu<Ngamma;mu++){
mat[mu+i*Ngamma+j*Lblock*Ngamma][t] = scalar_type(0.0);
}
}
}}
}
std::cout << GridLogMessage << " Done "<<std::endl;
// defer sum over nodes.
return;
}
template<class vobj>
void sliceInnerProductMesonFieldGammaMom(std::vector< std::vector<ComplexD> > &mat,
const std::vector<Lattice<vobj> > &lhs,
const std::vector<Lattice<vobj> > &rhs,
int orthogdim,
std::vector<Gamma::Algebra> gammas,
const std::vector<LatticeComplex > &mom)
{
typedef typename vobj::scalar_object sobj;
typedef typename vobj::scalar_type scalar_type;
typedef typename vobj::vector_type vector_type;
typedef iSpinMatrix<vector_type> SpinMatrix_v;
typedef iSpinMatrix<scalar_type> SpinMatrix_s;
int Lblock = lhs.size();
int Rblock = rhs.size();
GridBase *grid = lhs[0]._grid;
const int Nd = grid->_ndimension;
const int Nsimd = grid->Nsimd();
int Nt = grid->GlobalDimensions()[orthogdim];
int Ngamma = gammas.size();
int Nmom = mom.size();
assert(mat.size()==Lblock*Rblock*Ngamma*Nmom);
for(int t=0;t<mat.size();t++){
assert(mat[t].size()==Nt);
}
int fd=grid->_fdimensions[orthogdim];
int ld=grid->_ldimensions[orthogdim];
int rd=grid->_rdimensions[orthogdim];
// will locally sum vectors first
// sum across these down to scalars
// splitting the SIMD
int MFrvol = rd*Lblock*Rblock*Nmom;
int MFlvol = ld*Lblock*Rblock*Nmom;
Vector<SpinMatrix_v > lvSum(MFrvol);
parallel_for (int r = 0; r < MFrvol; r++){
lvSum[r] = zero;
}
Vector<SpinMatrix_s > lsSum(MFlvol);
parallel_for (int r = 0; r < MFlvol; r++){
lsSum[r]=scalar_type(0.0);
}
int e1= grid->_slice_nblock[orthogdim];
int e2= grid->_slice_block [orthogdim];
int stride=grid->_slice_stride[orthogdim];
std::cout << GridLogMessage << " Entering first parallel loop "<<std::endl;
// Parallelise over t-direction doesn't expose as much parallelism as needed for KNL
parallel_for(int r=0;r<rd;r++){
int so=r*grid->_ostride[orthogdim]; // base offset for start of plane
for(int n=0;n<e1;n++){
for(int b=0;b<e2;b++){
int ss= so+n*stride+b;
for(int i=0;i<Lblock;i++){
auto left = conjugate(lhs[i]._odata[ss]);
for(int j=0;j<Rblock;j++){
SpinMatrix_v vv;
auto right = rhs[j]._odata[ss];
for(int s1=0;s1<Ns;s1++){
for(int s2=0;s2<Ns;s2++){
vv()(s1,s2)() = left()(s1)(0) * right()(s2)(0)
+ left()(s1)(1) * right()(s2)(1)
+ left()(s1)(2) * right()(s2)(2);
}}
// After getting the sitewise product do the mom phase loop
int base = Nmom*i+Nmom*Lblock*j+Nmom*Lblock*Rblock*r;
// Trigger unroll
for ( int m=0;m<Nmom;m++){
int idx = m+base;
auto phase = mom[m]._odata[ss];
mac(&lvSum[idx],&vv,&phase);
}
}
}
}
}
}
std::cout << GridLogMessage << " Entering second parallel loop "<<std::endl;
// Sum across simd lanes in the plane, breaking out orthog dir.
parallel_for(int rt=0;rt<rd;rt++){
std::vector<int> icoor(Nd);
std::vector<SpinMatrix_s> extracted(Nsimd);
for(int i=0;i<Lblock;i++){
for(int j=0;j<Rblock;j++){
for(int m=0;m<Nmom;m++){
int ij_rdx = m+Nmom*i+Nmom*Lblock*j+Nmom*Lblock*Rblock*rt;
extract(lvSum[ij_rdx],extracted);
for(int idx=0;idx<Nsimd;idx++){
grid->iCoorFromIindex(icoor,idx);
int ldx = rt+icoor[orthogdim]*rd;
int ij_ldx = m+Nmom*i+Nmom*Lblock*j+Nmom*Lblock*Rblock*ldx;
lsSum[ij_ldx]=lsSum[ij_ldx]+extracted[idx];
}
}}}
}
std::cout << GridLogMessage << " Entering third parallel loop "<<std::endl;
parallel_for(int t=0;t<fd;t++)
{
int pt = t / ld; // processor plane
int lt = t % ld;
for(int i=0;i<Lblock;i++){
for(int j=0;j<Rblock;j++){
if (pt == grid->_processor_coor[orthogdim]){
for(int m=0;m<Nmom;m++){
int ij_dx = m+Nmom*i + Nmom*Lblock * j + Nmom*Lblock * Rblock * lt;
for(int mu=0;mu<Ngamma;mu++){
mat[ mu
+m*Ngamma
+i*Nmom*Ngamma
+j*Nmom*Ngamma*Lblock][t] = trace(lsSum[ij_dx]*Gamma(gammas[mu]));
}
}
}
else{
for(int mu=0;mu<Ngamma;mu++){
for(int m=0;m<Nmom;m++){
mat[mu+m*Ngamma+i*Nmom*Ngamma+j*Nmom*Lblock*Ngamma][t] = scalar_type(0.0);
}}
}
}}
}
std::cout << GridLogMessage << " Done "<<std::endl;
// defer sum over nodes.
return;
}
/*
template void sliceInnerProductMesonField<SpinColourVector>(std::vector< std::vector<ComplexD> > &mat,
const std::vector<Lattice<SpinColourVector> > &lhs,
const std::vector<Lattice<SpinColourVector> > &rhs,
int orthogdim) ;
*/
std::vector<Gamma::Algebra> Gmu4 ( {
Gamma::Algebra::GammaX,
Gamma::Algebra::GammaY,
Gamma::Algebra::GammaZ,
Gamma::Algebra::GammaT });
std::vector<Gamma::Algebra> Gmu16 ( {
Gamma::Algebra::Gamma5,
Gamma::Algebra::GammaT,
Gamma::Algebra::GammaTGamma5,
Gamma::Algebra::GammaX,
Gamma::Algebra::GammaXGamma5,
Gamma::Algebra::GammaY,
Gamma::Algebra::GammaYGamma5,
Gamma::Algebra::GammaZ,
Gamma::Algebra::GammaZGamma5,
Gamma::Algebra::Identity,
Gamma::Algebra::SigmaXT,
Gamma::Algebra::SigmaXY,
Gamma::Algebra::SigmaXZ,
Gamma::Algebra::SigmaYT,
Gamma::Algebra::SigmaYZ,
Gamma::Algebra::SigmaZT
});
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
std::vector<int> latt_size = GridDefaultLatt();
std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
std::vector<int> mpi_layout = GridDefaultMpi();
GridCartesian Grid(latt_size,simd_layout,mpi_layout);
const int Nmom=7;
int nt = latt_size[Tp];
uint64_t vol = 1;
for(int d=0;d<Nd;d++){
vol = vol*latt_size[d];
}
std::vector<int> seeds({1,2,3,4});
GridParallelRNG pRNG(&Grid);
pRNG.SeedFixedIntegers(seeds);
int Nm = atoi(argv[1]); // number of all modes (high + low)
std::vector<LatticeFermion> v(Nm,&Grid);
std::vector<LatticeFermion> w(Nm,&Grid);
std::vector<LatticeFermion> gammaV(Nm,&Grid);
std::vector<LatticeComplex> phases(Nmom,&Grid);
for(int i=0;i<Nm;i++) {
random(pRNG,v[i]);
random(pRNG,w[i]);
}
for(int i=0;i<Nmom;i++) {
phases[i] = Complex(1.0);
}
double flops = vol * (11.0 * 8.0 + 6.0) * Nm*Nm;
double byte = vol * (12.0 * sizeof(Complex) ) * Nm*Nm;
std::vector<ComplexD> ip(nt);
std::vector<std::vector<ComplexD> > MesonFields (Nm*Nm);
std::vector<std::vector<ComplexD> > MesonFields4 (Nm*Nm*4);
std::vector<std::vector<ComplexD> > MesonFields16 (Nm*Nm*16);
std::vector<std::vector<ComplexD> > MesonFields161(Nm*Nm*16);
std::vector<std::vector<ComplexD> > MesonFields16mom (Nm*Nm*16*Nmom);
std::vector<std::vector<ComplexD> > MesonFieldsRef(Nm*Nm);
for(int i=0;i<MesonFields.size();i++ ) MesonFields [i].resize(nt);
for(int i=0;i<MesonFieldsRef.size();i++) MesonFieldsRef[i].resize(nt);
for(int i=0;i<MesonFields4.size();i++ ) MesonFields4 [i].resize(nt);
for(int i=0;i<MesonFields16.size();i++ ) MesonFields16 [i].resize(nt);
for(int i=0;i<MesonFields161.size();i++ ) MesonFields161[i].resize(nt);
for(int i=0;i<MesonFields16mom.size();i++ ) MesonFields16mom [i].resize(nt);
GridLogMessage.TimingMode(1);
std::cout<<GridLogMessage << "Running loop with sliceInnerProductVector"<<std::endl;
double t0 = usecond();
for(int i=0;i<Nm;i++) {
for(int j=0;j<Nm;j++) {
sliceInnerProductVector(ip, w[i],v[j],Tp);
for(int t=0;t<nt;t++){
MesonFieldsRef[i+j*Nm][t] = ip[t];
}
}}
double t1 = usecond();
std::cout<<GridLogMessage << "Done "<< (t1-t0) <<" usecond " <<std::endl;
std::cout<<GridLogMessage << "Done "<< flops/(t1-t0) <<" mflops " <<std::endl;
std::cout<<GridLogMessage << "Done "<< byte /(t1-t0) <<" MB/s " <<std::endl;
std::cout<<GridLogMessage << "Running loop with new code for Nt="<<nt<<std::endl;
t0 = usecond();
sliceInnerProductMesonField(MesonFields,w,v,Tp);
t1 = usecond();
std::cout<<GridLogMessage << "Done "<< (t1-t0) <<" usecond " <<std::endl;
std::cout<<GridLogMessage << "Done "<< flops/(t1-t0) <<" mflops " <<std::endl;
std::cout<<GridLogMessage << "Done "<< byte /(t1-t0) <<" MB/s " <<std::endl;
std::cout<<GridLogMessage << "Running loop with Four gammas code for Nt="<<nt<<std::endl;
flops = vol * (11.0 * 8.0 + 6.0) * Nm*Nm*4;
byte = vol * (12.0 * sizeof(Complex) ) * Nm*Nm
+ vol * ( 2.0 * sizeof(Complex) ) * Nm*Nm* 4;
t0 = usecond();
sliceInnerProductMesonFieldGamma(MesonFields4,w,v,Tp,Gmu4);
t1 = usecond();
std::cout<<GridLogMessage << "Done "<< (t1-t0) <<" usecond " <<std::endl;
std::cout<<GridLogMessage << "Done "<< flops/(t1-t0) <<" mflops " <<std::endl;
std::cout<<GridLogMessage << "Done "<< byte /(t1-t0) <<" MB/s " <<std::endl;
std::cout<<GridLogMessage << "Running loop with Sixteen gammas code for Nt="<<nt<<std::endl;
flops = vol * (11.0 * 8.0 + 6.0) * Nm*Nm*16;
byte = vol * (12.0 * sizeof(Complex) ) * Nm*Nm
+ vol * ( 2.0 * sizeof(Complex) ) * Nm*Nm* 16;
t0 = usecond();
sliceInnerProductMesonFieldGamma(MesonFields16,w,v,Tp,Gmu16);
t1 = usecond();
std::cout<<GridLogMessage << "Done "<< (t1-t0) <<" usecond " <<std::endl;
std::cout<<GridLogMessage << "Done "<< flops/(t1-t0) <<" mflops " <<std::endl;
std::cout<<GridLogMessage << "Done "<< byte /(t1-t0) <<" MB/s " <<std::endl;
std::cout<<GridLogMessage << "Running loop with Sixteen gammas code1 for Nt="<<nt<<std::endl;
flops = vol * ( 2 * 8.0 + 6.0) * Nm*Nm*16;
byte = vol * (12.0 * sizeof(Complex) ) * Nm*Nm
+ vol * ( 2.0 * sizeof(Complex) ) * Nm*Nm* 16;
t0 = usecond();
sliceInnerProductMesonFieldGamma1(MesonFields161, w, v, Tp, Gmu16);
t1 = usecond();
std::cout<<GridLogMessage << "Done "<< (t1-t0) <<" usecond " <<std::endl;
std::cout<<GridLogMessage << "Done "<< flops/(t1-t0) <<" mflops " <<std::endl;
std::cout<<GridLogMessage << "Done "<< byte /(t1-t0) <<" MB/s " <<std::endl;
std::cout<<GridLogMessage << "Running loop with Sixteen gammas "<<Nmom<<" momenta "<<std::endl;
flops = vol * ( 2 * 8.0 + 6.0 + 8.0*Nmom) * Nm*Nm*16;
byte = vol * (12.0 * sizeof(Complex) ) * Nm*Nm
+ vol * ( 2.0 * sizeof(Complex) *Nmom ) * Nm*Nm* 16;
t0 = usecond();
sliceInnerProductMesonFieldGammaMom(MesonFields16mom,w,v,Tp,Gmu16,phases);
t1 = usecond();
std::cout<<GridLogMessage << "Done "<< (t1-t0) <<" usecond " <<std::endl;
std::cout<<GridLogMessage << "Done "<< flops/(t1-t0) <<" mflops " <<std::endl;
std::cout<<GridLogMessage << "Done "<< byte /(t1-t0) <<" MB/s " <<std::endl;
RealD err = 0;
RealD err2 = 0;
ComplexD diff;
ComplexD diff2;
for(int i=0;i<Nm;i++) {
for(int j=0;j<Nm;j++) {
for(int t=0;t<nt;t++){
diff = MesonFields[i+Nm*j][t] - MesonFieldsRef[i+Nm*j][t];
err += real(diff*conj(diff));
}
}}
std::cout<<GridLogMessage << "Norm error "<< err <<std::endl;
err = err*0.;
diff = diff*0.;
for (int mu = 0; mu < 16; mu++){
for (int k = 0; k < gammaV.size(); k++){
gammaV[k] = Gamma(Gmu16[mu]) * v[k];
}
for (int i = 0; i < Nm; i++){
for (int j = 0; j < Nm; j++){
sliceInnerProductVector(ip, w[i], gammaV[j], Tp);
for (int t = 0; t < nt; t++){
MesonFields[i + j * Nm][t] = ip[t];
diff = MesonFields16[mu+i*16+Nm*16*j][t] - MesonFields161[mu+i*16+Nm*16*j][t];
diff2 = MesonFields[i+j*Nm][t] - MesonFields161[mu+i*16+Nm*16*j][t];
err += real(diff*conj(diff));
err2 += real(diff2*conj(diff2));
}
}
}
}
std::cout << GridLogMessage << "Norm error 16 gamma1/16 gamma naive " << err << std::endl;
std::cout << GridLogMessage << "Norm error 16 gamma1/sliceInnerProduct " << err2 << std::endl;
Grid_finalize();
}

View File

@ -1,6 +1,6 @@
#!/usr/bin/env bash
EIGEN_URL='http://bitbucket.org/eigen/eigen/get/3.3.5.tar.bz2'
EIGEN_URL='http://bitbucket.org/eigen/eigen/get/3.3.3.tar.bz2'
echo "-- deploying Eigen source..."
wget ${EIGEN_URL} --no-check-certificate && ./scripts/update_eigen.sh `basename ${EIGEN_URL}` && rm `basename ${EIGEN_URL}`

View File

@ -480,8 +480,8 @@ GRID_LIBS=$LIBS
GRID_SHORT_SHA=`git rev-parse --short HEAD`
GRID_SHA=`git rev-parse HEAD`
GRID_BRANCH=`git rev-parse --abbrev-ref HEAD`
AM_CXXFLAGS="-I${abs_srcdir}/include -I${abs_srcdir}/Eigen/ -I${abs_srcdir}/Eigen/unsupported $AM_CXXFLAGS"
AM_CFLAGS="-I${abs_srcdir}/include -I${abs_srcdir}/Eigen/ -I${abs_srcdir}/Eigen/unsupported $AM_CFLAGS"
AM_CXXFLAGS="-I${abs_srcdir}/include $AM_CXXFLAGS"
AM_CFLAGS="-I${abs_srcdir}/include $AM_CFLAGS"
AM_LDFLAGS="-L${cwd}/lib $AM_LDFLAGS"
AC_SUBST([AM_CFLAGS])
AC_SUBST([AM_CXXFLAGS])

View File

@ -1,146 +0,0 @@
#ifndef A2A_Reduction_hpp_
#define A2A_Reduction_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Environment.hpp>
#include <Grid/Hadrons/Solver.hpp>
BEGIN_HADRONS_NAMESPACE
////////////////////////////////////////////
// A2A Meson Field Inner Product
////////////////////////////////////////////
template <class FermionField>
void sliceInnerProductMesonField(std::vector<std::vector<ComplexD>> &mat,
const std::vector<Lattice<FermionField>> &lhs,
const std::vector<Lattice<FermionField>> &rhs,
int orthogdim)
{
typedef typename FermionField::scalar_type scalar_type;
typedef typename FermionField::vector_type vector_type;
int Lblock = lhs.size();
int Rblock = rhs.size();
GridBase *grid = lhs[0]._grid;
const int Nd = grid->_ndimension;
const int Nsimd = grid->Nsimd();
int Nt = grid->GlobalDimensions()[orthogdim];
assert(mat.size() == Lblock * Rblock);
for (int t = 0; t < mat.size(); t++)
{
assert(mat[t].size() == Nt);
}
int fd = grid->_fdimensions[orthogdim];
int ld = grid->_ldimensions[orthogdim];
int rd = grid->_rdimensions[orthogdim];
// will locally sum vectors first
// sum across these down to scalars
// splitting the SIMD
std::vector<vector_type, alignedAllocator<vector_type>> lvSum(rd * Lblock * Rblock);
for(int r=0;r<rd * Lblock * Rblock;r++)
{
lvSum[r]=zero;
}
std::vector<scalar_type> lsSum(ld * Lblock * Rblock, scalar_type(0.0));
int e1 = grid->_slice_nblock[orthogdim];
int e2 = grid->_slice_block[orthogdim];
int stride = grid->_slice_stride[orthogdim];
// std::cout << GridLogMessage << " Entering first parallel loop " << std::endl;
// Parallelise over t-direction doesn't expose as much parallelism as needed for KNL
parallel_for(int r = 0; r < rd; r++)
{
int so = r * grid->_ostride[orthogdim]; // base offset for start of plane
for (int n = 0; n < e1; n++)
{
for (int b = 0; b < e2; b++)
{
int ss = so + n * stride + b;
for (int i = 0; i < Lblock; i++)
{
auto left = conjugate(lhs[i]._odata[ss]);
for (int j = 0; j < Rblock; j++)
{
int idx = i + Lblock * j + Lblock * Rblock * r;
auto right = rhs[j]._odata[ss];
vector_type vv = left()(0)(0) * right()(0)(0)
+ left()(0)(1) * right()(0)(1)
+ left()(0)(2) * right()(0)(2)
+ left()(1)(0) * right()(1)(0)
+ left()(1)(1) * right()(1)(1)
+ left()(1)(2) * right()(1)(2)
+ left()(2)(0) * right()(2)(0)
+ left()(2)(1) * right()(2)(1)
+ left()(2)(2) * right()(2)(2)
+ left()(3)(0) * right()(3)(0)
+ left()(3)(1) * right()(3)(1)
+ left()(3)(2) * right()(3)(2);
lvSum[idx] = lvSum[idx] + vv;
}
}
}
}
}
// std::cout << GridLogMessage << " Entering second parallel loop " << std::endl;
// Sum across simd lanes in the plane, breaking out orthog dir.
parallel_for(int rt = 0; rt < rd; rt++)
{
std::vector<int> icoor(Nd);
for (int i = 0; i < Lblock; i++)
{
for (int j = 0; j < Rblock; j++)
{
iScalar<vector_type> temp;
std::vector<iScalar<scalar_type>> extracted(Nsimd);
temp._internal = lvSum[i + Lblock * j + Lblock * Rblock * rt];
extract(temp, extracted);
for (int idx = 0; idx < Nsimd; idx++)
{
grid->iCoorFromIindex(icoor, idx);
int ldx = rt + icoor[orthogdim] * rd;
int ij_dx = i + Lblock * j + Lblock * Rblock * ldx;
lsSum[ij_dx] = lsSum[ij_dx] + extracted[idx]._internal;
}
}
}
}
// std::cout << GridLogMessage << " Entering non parallel loop " << std::endl;
for (int t = 0; t < fd; t++)
{
int pt = t/ld; // processor plane
int lt = t%ld;
for (int i = 0; i < Lblock; i++)
{
for (int j = 0; j < Rblock; j++)
{
if (pt == grid->_processor_coor[orthogdim])
{
int ij_dx = i + Lblock * j + Lblock * Rblock * lt;
mat[i + j * Lblock][t] = lsSum[ij_dx];
}
else
{
mat[i + j * Lblock][t] = scalar_type(0.0);
}
}
}
}
// std::cout << GridLogMessage << " Done " << std::endl;
// defer sum over nodes.
return;
}
END_HADRONS_NAMESPACE
#endif // A2A_Reduction_hpp_

View File

@ -1,210 +0,0 @@
#ifndef A2A_Vectors_hpp_
#define A2A_Vectors_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Environment.hpp>
#include <Grid/Hadrons/Solver.hpp>
BEGIN_HADRONS_NAMESPACE
////////////////////////////////
// A2A Modes
////////////////////////////////
template <class Field, class Matrix, class Solver>
class A2AModesSchurDiagTwo
{
private:
const std::vector<Field> *evec;
const std::vector<RealD> *eval;
Matrix &action;
Solver &solver;
std::vector<Field> w_high_5d, v_high_5d, w_high_4d, v_high_4d;
const int Nl, Nh;
const bool return_5d;
public:
int getNl (void ) {return Nl;}
int getNh (void ) {return Nh;}
int getN (void ) {return Nh+Nl;}
A2AModesSchurDiagTwo(const std::vector<Field> *_evec, const std::vector<RealD> *_eval,
Matrix &_action,
Solver &_solver,
std::vector<Field> _w_high_5d, std::vector<Field> _v_high_5d,
std::vector<Field> _w_high_4d, std::vector<Field> _v_high_4d,
const int _Nl, const int _Nh,
const bool _return_5d)
: evec(_evec), eval(_eval),
action(_action),
solver(_solver),
w_high_5d(_w_high_5d), v_high_5d(_v_high_5d),
w_high_4d(_w_high_4d), v_high_4d(_v_high_4d),
Nl(_Nl), Nh(_Nh),
return_5d(_return_5d){};
void high_modes(Field &source_5d, Field &w_source_5d, Field &source_4d, int i)
{
int i5d;
LOG(Message) << "A2A high modes for i = " << i << std::endl;
i5d = 0;
if (return_5d) i5d = i;
this->high_mode_v(action, solver, source_5d, v_high_5d[i5d], v_high_4d[i]);
this->high_mode_w(w_source_5d, source_4d, w_high_5d[i5d], w_high_4d[i]);
}
void return_v(int i, Field &vout_5d, Field &vout_4d)
{
if (i < Nl)
{
this->low_mode_v(action, evec->at(i), eval->at(i), vout_5d, vout_4d);
}
else
{
vout_4d = v_high_4d[i - Nl];
if (!(return_5d)) i = Nl;
vout_5d = v_high_5d[i - Nl];
}
}
void return_w(int i, Field &wout_5d, Field &wout_4d)
{
if (i < Nl)
{
this->low_mode_w(action, evec->at(i), eval->at(i), wout_5d, wout_4d);
}
else
{
wout_4d = w_high_4d[i - Nl];
if (!(return_5d)) i = Nl;
wout_5d = w_high_5d[i - Nl];
}
}
void low_mode_v(Matrix &action, const Field &evec, const RealD &eval, Field &vout_5d, Field &vout_4d)
{
GridBase *grid = action.RedBlackGrid();
Field src_o(grid);
Field sol_e(grid);
Field sol_o(grid);
Field tmp(grid);
src_o = evec;
src_o.checkerboard = Odd;
pickCheckerboard(Even, sol_e, vout_5d);
pickCheckerboard(Odd, sol_o, vout_5d);
/////////////////////////////////////////////////////
// v_ie = -(1/eval_i) * MeeInv Meo MooInv evec_i
/////////////////////////////////////////////////////
action.MooeeInv(src_o, tmp);
assert(tmp.checkerboard == Odd);
action.Meooe(tmp, sol_e);
assert(sol_e.checkerboard == Even);
action.MooeeInv(sol_e, tmp);
assert(tmp.checkerboard == Even);
sol_e = (-1.0 / eval) * tmp;
assert(sol_e.checkerboard == Even);
/////////////////////////////////////////////////////
// v_io = (1/eval_i) * MooInv evec_i
/////////////////////////////////////////////////////
action.MooeeInv(src_o, tmp);
assert(tmp.checkerboard == Odd);
sol_o = (1.0 / eval) * tmp;
assert(sol_o.checkerboard == Odd);
setCheckerboard(vout_5d, sol_e);
assert(sol_e.checkerboard == Even);
setCheckerboard(vout_5d, sol_o);
assert(sol_o.checkerboard == Odd);
action.ExportPhysicalFermionSolution(vout_5d, vout_4d);
}
void low_mode_w(Matrix &action, const Field &evec, const RealD &eval, Field &wout_5d, Field &wout_4d)
{
GridBase *grid = action.RedBlackGrid();
SchurDiagTwoOperator<Matrix, Field> _HermOpEO(action);
Field src_o(grid);
Field sol_e(grid);
Field sol_o(grid);
Field tmp(grid);
GridBase *fgrid = action.Grid();
Field tmp_wout(fgrid);
src_o = evec;
src_o.checkerboard = Odd;
pickCheckerboard(Even, sol_e, tmp_wout);
pickCheckerboard(Odd, sol_o, tmp_wout);
/////////////////////////////////////////////////////
// w_ie = - MeeInvDag MoeDag Doo evec_i
/////////////////////////////////////////////////////
_HermOpEO.Mpc(src_o, tmp);
assert(tmp.checkerboard == Odd);
action.MeooeDag(tmp, sol_e);
assert(sol_e.checkerboard == Even);
action.MooeeInvDag(sol_e, tmp);
assert(tmp.checkerboard == Even);
sol_e = (-1.0) * tmp;
/////////////////////////////////////////////////////
// w_io = Doo evec_i
/////////////////////////////////////////////////////
_HermOpEO.Mpc(src_o, sol_o);
assert(sol_o.checkerboard == Odd);
setCheckerboard(tmp_wout, sol_e);
assert(sol_e.checkerboard == Even);
setCheckerboard(tmp_wout, sol_o);
assert(sol_o.checkerboard == Odd);
action.DminusDag(tmp_wout, wout_5d);
action.ExportPhysicalFermionSource(wout_5d, wout_4d);
}
void high_mode_v(Matrix &action, Solver &solver, const Field &source, Field &vout_5d, Field &vout_4d)
{
GridBase *fgrid = action.Grid();
solver(vout_5d, source); // Note: solver is solver(out, in)
action.ExportPhysicalFermionSolution(vout_5d, vout_4d);
}
void high_mode_w(const Field &w_source_5d, const Field &source_4d, Field &wout_5d, Field &wout_4d)
{
wout_5d = w_source_5d;
wout_4d = source_4d;
}
};
// TODO: A2A for coarse eigenvectors
// template <class FineField, class CoarseField, class Matrix, class Solver>
// class A2ALMSchurDiagTwoCoarse : public A2AModesSchurDiagTwo<FineField, Matrix, Solver>
// {
// private:
// const std::vector<FineField> &subspace;
// const std::vector<CoarseField> &evec_coarse;
// const std::vector<RealD> &eval_coarse;
// Matrix &action;
// public:
// A2ALMSchurDiagTwoCoarse(const std::vector<FineField> &_subspace, const std::vector<CoarseField> &_evec_coarse, const std::vector<RealD> &_eval_coarse, Matrix &_action)
// : subspace(_subspace), evec_coarse(_evec_coarse), eval_coarse(_eval_coarse), action(_action){};
// void operator()(int i, FineField &vout, FineField &wout)
// {
// FineField prom_evec(subspace[0]._grid);
// blockPromote(evec_coarse[i], prom_evec, subspace);
// this->low_mode_v(action, prom_evec, eval_coarse[i], vout);
// this->low_mode_w(action, prom_evec, eval_coarse[i], wout);
// }
// };
END_HADRONS_NAMESPACE
#endif // A2A_Vectors_hpp_

View File

@ -28,7 +28,6 @@ See the full license in the file "LICENSE" in the top level distribution directo
#include <Grid/Hadrons/Application.hpp>
#include <Grid/Hadrons/GeneticScheduler.hpp>
#include <Grid/Hadrons/Modules.hpp>
using namespace Grid;
using namespace QCD;
@ -41,9 +40,6 @@ using namespace Hadrons;
* Application implementation *
******************************************************************************/
// constructors ////////////////////////////////////////////////////////////////
#define MACOUT(macro) macro << " (" << #macro << ")"
#define MACOUTS(macro) HADRONS_STR(macro) << " (" << #macro << ")"
Application::Application(void)
{
initLogger();
@ -54,22 +50,9 @@ Application::Application(void)
loc[d] /= mpi[d];
locVol_ *= loc[d];
}
LOG(Message) << "====== HADRONS APPLICATION STARTING ======" << std::endl;
LOG(Message) << "** Dimensions" << std::endl;
LOG(Message) << "Global lattice : " << dim << std::endl;
LOG(Message) << "MPI partition : " << mpi << std::endl;
LOG(Message) << "Local lattice : " << loc << std::endl;
LOG(Message) << std::endl;
LOG(Message) << "** Default parameters (and associated C macro)" << std::endl;
LOG(Message) << "ASCII output precision : " << MACOUT(DEFAULT_ASCII_PREC) << std::endl;
LOG(Message) << "Fermion implementation : " << MACOUTS(FIMPL) << std::endl;
LOG(Message) << "z-Fermion implementation: " << MACOUTS(ZFIMPL) << std::endl;
LOG(Message) << "Scalar implementation : " << MACOUTS(SIMPL) << std::endl;
LOG(Message) << "Gauge implementation : " << MACOUTS(GIMPL) << std::endl;
LOG(Message) << "Eigenvector base size : "
<< MACOUT(HADRONS_DEFAULT_LANCZOS_NBASIS) << std::endl;
LOG(Message) << "Schur decomposition : " << MACOUTS(HADRONS_DEFAULT_SCHUR) << std::endl;
LOG(Message) << std::endl;
LOG(Message) << "Global lattice: " << dim << std::endl;
LOG(Message) << "MPI partition : " << mpi << std::endl;
LOG(Message) << "Local lattice : " << loc << std::endl;
}
Application::Application(const Application::GlobalPar &par)
@ -130,12 +113,12 @@ void Application::parseParameterFile(const std::string parameterFileName)
setPar(par);
if (!push(reader, "modules"))
{
HADRONS_ERROR(Parsing, "Cannot open node 'modules' in parameter file '"
HADRON_ERROR(Parsing, "Cannot open node 'modules' in parameter file '"
+ parameterFileName + "'");
}
if (!push(reader, "module"))
{
HADRONS_ERROR(Parsing, "Cannot open node 'modules/module' in parameter file '"
HADRON_ERROR(Parsing, "Cannot open node 'modules/module' in parameter file '"
+ parameterFileName + "'");
}
do
@ -149,27 +132,24 @@ void Application::parseParameterFile(const std::string parameterFileName)
void Application::saveParameterFile(const std::string parameterFileName)
{
XmlWriter writer(parameterFileName);
ObjectId id;
const unsigned int nMod = vm().getNModule();
LOG(Message) << "Saving application to '" << parameterFileName << "'..." << std::endl;
if (env().getGrid()->IsBoss())
write(writer, "parameters", getPar());
push(writer, "modules");
for (unsigned int i = 0; i < nMod; ++i)
{
XmlWriter writer(parameterFileName);
ObjectId id;
const unsigned int nMod = vm().getNModule();
write(writer, "parameters", getPar());
push(writer, "modules");
for (unsigned int i = 0; i < nMod; ++i)
{
push(writer, "module");
id.name = vm().getModuleName(i);
id.type = vm().getModule(i)->getRegisteredName();
write(writer, "id", id);
vm().getModule(i)->saveParameters(writer, "options");
pop(writer);
}
pop(writer);
push(writer, "module");
id.name = vm().getModuleName(i);
id.type = vm().getModule(i)->getRegisteredName();
write(writer, "id", id);
vm().getModule(i)->saveParameters(writer, "options");
pop(writer);
}
pop(writer);
pop(writer);
}
// schedule computation ////////////////////////////////////////////////////////
@ -184,24 +164,20 @@ void Application::schedule(void)
void Application::saveSchedule(const std::string filename)
{
TextWriter writer(filename);
std::vector<std::string> program;
if (!scheduled_)
{
HADRON_ERROR(Definition, "Computation not scheduled");
}
LOG(Message) << "Saving current schedule to '" << filename << "'..."
<< std::endl;
if (env().getGrid()->IsBoss())
for (auto address: program_)
{
TextWriter writer(filename);
std::vector<std::string> program;
if (!scheduled_)
{
HADRONS_ERROR(Definition, "Computation not scheduled");
}
for (auto address: program_)
{
program.push_back(vm().getModuleName(address));
}
write(writer, "schedule", program);
program.push_back(vm().getModuleName(address));
}
write(writer, "schedule", program);
}
void Application::loadSchedule(const std::string filename)
@ -224,7 +200,7 @@ void Application::printSchedule(void)
{
if (!scheduled_)
{
HADRONS_ERROR(Definition, "Computation not scheduled");
HADRON_ERROR(Definition, "Computation not scheduled");
}
auto peak = vm().memoryNeeded(program_);
LOG(Message) << "Schedule (memory needed: " << sizeString(peak) << "):"

View File

@ -31,7 +31,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/VirtualMachine.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/Modules.hpp>
BEGIN_HADRONS_NAMESPACE

View File

@ -1,323 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/EigenPack.hpp
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_EigenPack_hpp_
#define Hadrons_EigenPack_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/algorithms/iterative/Deflation.h>
#include <Grid/algorithms/iterative/LocalCoherenceLanczos.h>
BEGIN_HADRONS_NAMESPACE
// Lanczos type
#ifndef HADRONS_DEFAULT_LANCZOS_NBASIS
#define HADRONS_DEFAULT_LANCZOS_NBASIS 60
#endif
template <typename F>
class EigenPack
{
public:
typedef F Field;
struct PackRecord
{
std::string operatorXml, solverXml;
};
struct VecRecord: Serializable
{
GRID_SERIALIZABLE_CLASS_MEMBERS(VecRecord,
unsigned int, index,
double, eval);
VecRecord(void): index(0), eval(0.) {}
};
public:
std::vector<RealD> eval;
std::vector<F> evec;
PackRecord record;
public:
EigenPack(void) = default;
virtual ~EigenPack(void) = default;
EigenPack(const size_t size, GridBase *grid)
{
resize(size, grid);
}
void resize(const size_t size, GridBase *grid)
{
eval.resize(size);
evec.resize(size, grid);
}
virtual void read(const std::string fileStem, const bool multiFile, const int traj = -1)
{
if (multiFile)
{
for(int k = 0; k < evec.size(); ++k)
{
basicReadSingle(evec[k], eval[k], evecFilename(fileStem, k, traj), k);
}
}
else
{
basicRead(evec, eval, evecFilename(fileStem, -1, traj), evec.size());
}
}
virtual void write(const std::string fileStem, const bool multiFile, const int traj = -1)
{
if (multiFile)
{
for(int k = 0; k < evec.size(); ++k)
{
basicWriteSingle(evecFilename(fileStem, k, traj), evec[k], eval[k], k);
}
}
else
{
basicWrite(evecFilename(fileStem, -1, traj), evec, eval, evec.size());
}
}
protected:
std::string evecFilename(const std::string stem, const int vec, const int traj)
{
std::string t = (traj < 0) ? "" : ("." + std::to_string(traj));
if (vec == -1)
{
return stem + t + ".bin";
}
else
{
return stem + t + "/v" + std::to_string(vec) + ".bin";
}
}
template <typename T>
void basicRead(std::vector<T> &evec, std::vector<double> &eval,
const std::string filename, const unsigned int size)
{
ScidacReader binReader;
binReader.open(filename);
binReader.skipPastObjectRecord(SCIDAC_FILE_XML);
for(int k = 0; k < size; ++k)
{
VecRecord vecRecord;
LOG(Message) << "Reading eigenvector " << k << std::endl;
binReader.readScidacFieldRecord(evec[k], vecRecord);
if (vecRecord.index != k)
{
HADRONS_ERROR(Io, "Eigenvector " + std::to_string(k) + " has a"
+ " wrong index (expected " + std::to_string(vecRecord.index)
+ ") in file '" + filename + "'");
}
eval[k] = vecRecord.eval;
}
binReader.close();
}
template <typename T>
void basicReadSingle(T &evec, double &eval, const std::string filename,
const unsigned int index)
{
ScidacReader binReader;
VecRecord vecRecord;
binReader.open(filename);
binReader.skipPastObjectRecord(SCIDAC_FILE_XML);
LOG(Message) << "Reading eigenvector " << index << std::endl;
binReader.readScidacFieldRecord(evec, vecRecord);
if (vecRecord.index != index)
{
HADRONS_ERROR(Io, "Eigenvector " + std::to_string(index) + " has a"
+ " wrong index (expected " + std::to_string(vecRecord.index)
+ ") in file '" + filename + "'");
}
eval = vecRecord.eval;
binReader.close();
}
template <typename T>
void basicWrite(const std::string filename, std::vector<T> &evec,
const std::vector<double> &eval, const unsigned int size)
{
ScidacWriter binWriter(evec[0]._grid->IsBoss());
XmlWriter xmlWriter("", "eigenPackPar");
makeFileDir(filename, evec[0]._grid);
xmlWriter.pushXmlString(record.operatorXml);
xmlWriter.pushXmlString(record.solverXml);
binWriter.open(filename);
binWriter.writeLimeObject(1, 1, xmlWriter, "parameters", SCIDAC_FILE_XML);
for(int k = 0; k < size; ++k)
{
VecRecord vecRecord;
vecRecord.index = k;
vecRecord.eval = eval[k];
LOG(Message) << "Writing eigenvector " << k << std::endl;
binWriter.writeScidacFieldRecord(evec[k], vecRecord, DEFAULT_ASCII_PREC);
}
binWriter.close();
}
template <typename T>
void basicWriteSingle(const std::string filename, T &evec,
const double eval, const unsigned int index)
{
ScidacWriter binWriter(evec._grid->IsBoss());
XmlWriter xmlWriter("", "eigenPackPar");
VecRecord vecRecord;
makeFileDir(filename, evec._grid);
xmlWriter.pushXmlString(record.operatorXml);
xmlWriter.pushXmlString(record.solverXml);
binWriter.open(filename);
binWriter.writeLimeObject(1, 1, xmlWriter, "parameters", SCIDAC_FILE_XML);
vecRecord.index = index;
vecRecord.eval = eval;
LOG(Message) << "Writing eigenvector " << index << std::endl;
binWriter.writeScidacFieldRecord(evec, vecRecord, DEFAULT_ASCII_PREC);
binWriter.close();
}
};
template <typename FineF, typename CoarseF>
class CoarseEigenPack: public EigenPack<FineF>
{
public:
typedef CoarseF CoarseField;
public:
std::vector<RealD> evalCoarse;
std::vector<CoarseF> evecCoarse;
public:
CoarseEigenPack(void) = default;
virtual ~CoarseEigenPack(void) = default;
CoarseEigenPack(const size_t sizeFine, const size_t sizeCoarse,
GridBase *gridFine, GridBase *gridCoarse)
{
resize(sizeFine, sizeCoarse, gridFine, gridCoarse);
}
void resize(const size_t sizeFine, const size_t sizeCoarse,
GridBase *gridFine, GridBase *gridCoarse)
{
EigenPack<FineF>::resize(sizeFine, gridFine);
evalCoarse.resize(sizeCoarse);
evecCoarse.resize(sizeCoarse, gridCoarse);
}
void readFine(const std::string fileStem, const bool multiFile, const int traj = -1)
{
if (multiFile)
{
for(int k = 0; k < this->evec.size(); ++k)
{
this->basicReadSingle(this->evec[k], this->eval[k], this->evecFilename(fileStem + "_fine", k, traj), k);
}
}
else
{
this->basicRead(this->evec, this->eval, this->evecFilename(fileStem + "_fine", -1, traj), this->evec.size());
}
}
void readCoarse(const std::string fileStem, const bool multiFile, const int traj = -1)
{
if (multiFile)
{
for(int k = 0; k < evecCoarse.size(); ++k)
{
this->basicReadSingle(evecCoarse[k], evalCoarse[k], this->evecFilename(fileStem + "_coarse", k, traj), k);
}
}
else
{
this->basicRead(evecCoarse, evalCoarse, this->evecFilename(fileStem + "_coarse", -1, traj), evecCoarse.size());
}
}
virtual void read(const std::string fileStem, const bool multiFile, const int traj = -1)
{
readFine(fileStem, multiFile, traj);
readCoarse(fileStem, multiFile, traj);
}
void writeFine(const std::string fileStem, const bool multiFile, const int traj = -1)
{
if (multiFile)
{
for(int k = 0; k < this->evec.size(); ++k)
{
this->basicWriteSingle(this->evecFilename(fileStem + "_fine", k, traj), this->evec[k], this->eval[k], k);
}
}
else
{
this->basicWrite(this->evecFilename(fileStem + "_fine", -1, traj), this->evec, this->eval, this->evec.size());
}
}
void writeCoarse(const std::string fileStem, const bool multiFile, const int traj = -1)
{
if (multiFile)
{
for(int k = 0; k < evecCoarse.size(); ++k)
{
this->basicWriteSingle(this->evecFilename(fileStem + "_coarse", k, traj), evecCoarse[k], evalCoarse[k], k);
}
}
else
{
this->basicWrite(this->evecFilename(fileStem + "_coarse", -1, traj), evecCoarse, evalCoarse, evecCoarse.size());
}
}
virtual void write(const std::string fileStem, const bool multiFile, const int traj = -1)
{
writeFine(fileStem, multiFile, traj);
writeCoarse(fileStem, multiFile, traj);
}
};
template <typename FImpl>
using FermionEigenPack = EigenPack<typename FImpl::FermionField>;
template <typename FImpl, int nBasis>
using CoarseFermionEigenPack = CoarseEigenPack<
typename FImpl::FermionField,
typename LocalCoherenceLanczos<typename FImpl::SiteSpinor,
typename FImpl::SiteComplex,
nBasis>::CoarseField>;
END_HADRONS_NAMESPACE
#endif // Hadrons_EigenPack_hpp_

View File

@ -35,7 +35,7 @@ using namespace QCD;
using namespace Hadrons;
#define ERROR_NO_ADDRESS(address)\
HADRONS_ERROR(Definition, "no object with address " + std::to_string(address));
HADRON_ERROR(Definition, "no object with address " + std::to_string(address));
/******************************************************************************
* Environment implementation *
@ -49,10 +49,11 @@ Environment::Environment(void)
dim_, GridDefaultSimd(nd_, vComplex::Nsimd()),
GridDefaultMpi()));
gridRb4d_.reset(SpaceTimeGrid::makeFourDimRedBlackGrid(grid4d_.get()));
vol_ = 1.;
for (auto d: dim_)
auto loc = getGrid()->LocalDimensions();
locVol_ = 1;
for (unsigned int d = 0; d < loc.size(); ++d)
{
vol_ *= d;
locVol_ *= loc[d];
}
rng4d_.reset(new GridParallelRNG(grid4d_.get()));
}
@ -60,7 +61,7 @@ Environment::Environment(void)
// grids ///////////////////////////////////////////////////////////////////////
void Environment::createGrid(const unsigned int Ls)
{
if ((Ls > 1) and (grid5d_.find(Ls) == grid5d_.end()))
if (grid5d_.find(Ls) == grid5d_.end())
{
auto g = getGrid();
@ -69,49 +70,6 @@ void Environment::createGrid(const unsigned int Ls)
}
}
void Environment::createCoarseGrid(const std::vector<int> &blockSize,
const unsigned int Ls)
{
int nd = getNd();
std::vector<int> fineDim = getDim(), coarseDim;
unsigned int cLs;
auto key4d = blockSize, key5d = blockSize;
createGrid(Ls);
coarseDim.resize(nd);
for (int d = 0; d < coarseDim.size(); d++)
{
coarseDim[d] = fineDim[d]/blockSize[d];
if (coarseDim[d]*blockSize[d] != fineDim[d])
{
HADRONS_ERROR(Size, "Fine dimension " + std::to_string(d)
+ " (" + std::to_string(fineDim[d])
+ ") not divisible by coarse dimension ("
+ std::to_string(coarseDim[d]) + ")");
}
}
if (blockSize.size() > nd)
{
cLs = Ls/blockSize[nd];
if (cLs*blockSize[nd] != Ls)
{
HADRONS_ERROR(Size, "Fine Ls (" + std::to_string(Ls)
+ ") not divisible by coarse Ls ("
+ std::to_string(cLs) + ")");
}
key4d.resize(nd);
key5d.push_back(Ls);
}
gridCoarse4d_[key4d].reset(
SpaceTimeGrid::makeFourDimGrid(coarseDim,
GridDefaultSimd(nd, vComplex::Nsimd()), GridDefaultMpi()));
if (Ls > 1)
{
gridCoarse5d_[key5d].reset(
SpaceTimeGrid::makeFiveDimGrid(cLs, gridCoarse4d_[key4d].get()));
}
}
GridCartesian * Environment::getGrid(const unsigned int Ls) const
{
try
@ -127,7 +85,7 @@ GridCartesian * Environment::getGrid(const unsigned int Ls) const
}
catch(std::out_of_range &)
{
HADRONS_ERROR(Definition, "no grid with Ls= " + std::to_string(Ls));
HADRON_ERROR(Definition, "no grid with Ls= " + std::to_string(Ls));
}
}
@ -146,31 +104,7 @@ GridRedBlackCartesian * Environment::getRbGrid(const unsigned int Ls) const
}
catch(std::out_of_range &)
{
HADRONS_ERROR(Definition, "no red-black grid with Ls= " + std::to_string(Ls));
}
}
GridCartesian * Environment::getCoarseGrid(
const std::vector<int> &blockSize, const unsigned int Ls) const
{
auto key = blockSize;
try
{
if (Ls == 1)
{
key.resize(getNd());
return gridCoarse4d_.at(key).get();
}
else
{
key.push_back(Ls);
return gridCoarse5d_.at(key).get();
}
}
catch(std::out_of_range &)
{
HADRONS_ERROR(Definition, "no coarse grid with Ls= " + std::to_string(Ls));
HADRON_ERROR(Definition, "no red-black 5D grid with Ls= " + std::to_string(Ls));
}
}
@ -189,9 +123,9 @@ int Environment::getDim(const unsigned int mu) const
return dim_[mu];
}
double Environment::getVolume(void) const
unsigned long int Environment::getLocalVolume(void) const
{
return vol_;
return locVol_;
}
// random number generator /////////////////////////////////////////////////////
@ -220,7 +154,7 @@ void Environment::addObject(const std::string name, const int moduleAddress)
}
else
{
HADRONS_ERROR(Definition, "object '" + name + "' already exists");
HADRON_ERROR(Definition, "object '" + name + "' already exists");
}
}
@ -243,7 +177,7 @@ unsigned int Environment::getObjectAddress(const std::string name) const
}
else
{
HADRONS_ERROR(Definition, "no object with name '" + name + "'");
HADRON_ERROR(Definition, "no object with name '" + name + "'");
}
}
@ -336,7 +270,7 @@ int Environment::getObjectModule(const std::string name) const
unsigned int Environment::getObjectLs(const unsigned int address) const
{
if (hasCreatedObject(address))
if (hasObject(address))
{
return object_[address].Ls;
}

View File

@ -78,7 +78,7 @@ private:
Size size{0};
Storage storage{Storage::object};
unsigned int Ls{0};
const std::type_info *type{nullptr}, *derivedType{nullptr};
const std::type_info *type{nullptr};
std::string name;
int module{-1};
std::unique_ptr<Object> data{nullptr};
@ -86,16 +86,12 @@ private:
public:
// grids
void createGrid(const unsigned int Ls);
void createCoarseGrid(const std::vector<int> &blockSize,
const unsigned int Ls = 1);
GridCartesian * getGrid(const unsigned int Ls = 1) const;
GridRedBlackCartesian * getRbGrid(const unsigned int Ls = 1) const;
GridCartesian * getCoarseGrid(const std::vector<int> &blockSize,
const unsigned int Ls = 1) const;
std::vector<int> getDim(void) const;
int getDim(const unsigned int mu) const;
unsigned long int getLocalVolume(void) const;
unsigned int getNd(void) const;
double getVolume(void) const;
// random number generator
void setSeed(const std::vector<int> &seed);
GridParallelRNG * get4dRng(void) const;
@ -114,10 +110,6 @@ public:
Ts && ... args);
void setObjectModule(const unsigned int objAddress,
const int modAddress);
template <typename B, typename T>
T * getDerivedObject(const unsigned int address) const;
template <typename B, typename T>
T * getDerivedObject(const std::string name) const;
template <typename T>
T * getObject(const unsigned int address) const;
template <typename T>
@ -155,7 +147,7 @@ public:
void printContent(void) const;
private:
// general
double vol_;
unsigned long int locVol_;
bool protect_{true};
// grids
std::vector<int> dim_;
@ -163,8 +155,6 @@ private:
std::map<unsigned int, GridPt> grid5d_;
GridRbPt gridRb4d_;
std::map<unsigned int, GridRbPt> gridRb5d_;
std::map<std::vector<int>, GridPt> gridCoarse4d_;
std::map<std::vector<int>, GridPt> gridCoarse5d_;
unsigned int nd_;
// random number generator
RngPt rng4d_;
@ -186,7 +176,7 @@ Holder<T>::Holder(T *pt)
template <typename T>
T & Holder<T>::get(void) const
{
return *objPt_.get();
return &objPt_.get();
}
template <typename T>
@ -226,26 +216,24 @@ void Environment::createDerivedObject(const std::string name,
{
MemoryProfiler::stats = &memStats;
}
size_t initMem = MemoryProfiler::stats->currentlyAllocated;
object_[address].storage = storage;
object_[address].Ls = Ls;
size_t initMem = MemoryProfiler::stats->currentlyAllocated;
object_[address].storage = storage;
object_[address].Ls = Ls;
object_[address].data.reset(new Holder<B>(new T(std::forward<Ts>(args)...)));
object_[address].size = MemoryProfiler::stats->maxAllocated - initMem;
object_[address].type = &typeid(B);
object_[address].derivedType = &typeid(T);
object_[address].size = MemoryProfiler::stats->maxAllocated - initMem;
object_[address].type = &typeid(T);
if (MemoryProfiler::stats == &memStats)
{
MemoryProfiler::stats = nullptr;
}
}
// object already exists, no error if it is a cache, error otherwise
else if ((object_[address].storage != Storage::cache) or
(object_[address].storage != storage) or
(object_[address].name != name) or
(object_[address].type != &typeid(B)) or
(object_[address].derivedType != &typeid(T)))
else if ((object_[address].storage != Storage::cache) or
(object_[address].storage != storage) or
(object_[address].name != name) or
(object_[address].type != &typeid(T)))
{
HADRONS_ERROR(Definition, "object '" + name + "' already allocated");
HADRON_ERROR(Definition, "object '" + name + "' already allocated");
}
}
@ -258,64 +246,36 @@ void Environment::createObject(const std::string name,
createDerivedObject<T, T>(name, storage, Ls, std::forward<Ts>(args)...);
}
template <typename B, typename T>
T * Environment::getDerivedObject(const unsigned int address) const
template <typename T>
T * Environment::getObject(const unsigned int address) const
{
if (hasObject(address))
{
if (hasCreatedObject(address))
{
if (auto h = dynamic_cast<Holder<B> *>(object_[address].data.get()))
if (auto h = dynamic_cast<Holder<T> *>(object_[address].data.get()))
{
if (&typeid(T) == &typeid(B))
{
return dynamic_cast<T *>(h->getPt());
}
else
{
if (auto hder = dynamic_cast<T *>(h->getPt()))
{
return hder;
}
else
{
HADRONS_ERROR(Definition, "object with address " + std::to_string(address) +
" cannot be casted to '" + typeName(&typeid(T)) +
"' (has type '" + typeName(&typeid(h->get())) + "')");
}
}
return h->getPt();
}
else
{
HADRONS_ERROR(Definition, "object with address " + std::to_string(address) +
" does not have type '" + typeName(&typeid(B)) +
HADRON_ERROR(Definition, "object with address " + std::to_string(address) +
" does not have type '" + typeName(&typeid(T)) +
"' (has type '" + getObjectType(address) + "')");
}
}
else
{
HADRONS_ERROR(Definition, "object with address " + std::to_string(address) +
HADRON_ERROR(Definition, "object with address " + std::to_string(address) +
" is empty");
}
}
else
{
HADRONS_ERROR(Definition, "no object with address " + std::to_string(address));
HADRON_ERROR(Definition, "no object with address " + std::to_string(address));
}
}
template <typename B, typename T>
T * Environment::getDerivedObject(const std::string name) const
{
return getDerivedObject<B, T>(getObjectAddress(name));
}
template <typename T>
T * Environment::getObject(const unsigned int address) const
{
return getDerivedObject<T, T>(address);
}
template <typename T>
T * Environment::getObject(const std::string name) const
{
@ -338,7 +298,7 @@ bool Environment::isObjectOfType(const unsigned int address) const
}
else
{
HADRONS_ERROR(Definition, "no object with address " + std::to_string(address));
HADRON_ERROR(Definition, "no object with address " + std::to_string(address));
}
}

View File

@ -27,8 +27,6 @@ See the full license in the file "LICENSE" in the top level distribution directo
/* END LEGAL */
#include <Grid/Hadrons/Exceptions.hpp>
#include <Grid/Hadrons/VirtualMachine.hpp>
#include <Grid/Hadrons/Module.hpp>
#ifndef ERR_SUFF
#define ERR_SUFF " (" + loc + ")"
@ -49,7 +47,6 @@ CONST_EXC(Definition, Logic("definition error: " + msg, loc))
CONST_EXC(Implementation, Logic("implementation error: " + msg, loc))
CONST_EXC(Range, Logic("range error: " + msg, loc))
CONST_EXC(Size, Logic("size error: " + msg, loc))
// runtime errors
CONST_EXC(Runtime, runtime_error(msg + ERR_SUFF))
CONST_EXC(Argument, Runtime("argument error: " + msg, loc))
@ -58,24 +55,3 @@ CONST_EXC(Memory, Runtime("memory error: " + msg, loc))
CONST_EXC(Parsing, Runtime("parsing error: " + msg, loc))
CONST_EXC(Program, Runtime("program error: " + msg, loc))
CONST_EXC(System, Runtime("system error: " + msg, loc))
// abort functions
void Grid::Hadrons::Exceptions::abort(const std::exception& e)
{
auto &vm = VirtualMachine::getInstance();
int mod = vm.getCurrentModule();
LOG(Error) << "FATAL ERROR -- Exception " << typeName(&typeid(e))
<< std::endl;
if (mod >= 0)
{
LOG(Error) << "During execution of module '"
<< vm.getModuleName(mod) << "' (address " << mod << ")"
<< std::endl;
}
LOG(Error) << e.what() << std::endl;
LOG(Error) << "Aborting program" << std::endl;
Grid_finalize();
exit(EXIT_FAILURE);
}

View File

@ -34,10 +34,11 @@ See the full license in the file "LICENSE" in the top level distribution directo
#include <Grid/Hadrons/Global.hpp>
#endif
#define HADRONS_SRC_LOC std::string(__FUNCTION__) + " at " \
+ std::string(__FILE__) + ":" + std::to_string(__LINE__)
#define HADRONS_ERROR(exc, msg)\
throw(Exceptions::exc(msg, HADRONS_SRC_LOC));
#define SRC_LOC std::string(__FUNCTION__) + " at " + std::string(__FILE__) + ":"\
+ std::to_string(__LINE__)
#define HADRON_ERROR(exc, msg)\
LOG(Error) << msg << std::endl;\
throw(Exceptions::exc(msg, SRC_LOC));
#define DECL_EXC(name, base) \
class name: public base\
@ -56,7 +57,6 @@ namespace Exceptions
DECL_EXC(Implementation, Logic);
DECL_EXC(Range, Logic);
DECL_EXC(Size, Logic);
// runtime errors
DECL_EXC(Runtime, std::runtime_error);
DECL_EXC(Argument, Runtime);
@ -65,9 +65,6 @@ namespace Exceptions
DECL_EXC(Parsing, Runtime);
DECL_EXC(Program, Runtime);
DECL_EXC(System, Runtime);
// abort functions
void abort(const std::exception& e);
}
END_HADRONS_NAMESPACE

View File

@ -94,7 +94,7 @@ std::unique_ptr<T> Factory<T>::create(const std::string type,
}
catch (std::out_of_range &)
{
HADRONS_ERROR(Argument, "object of type '" + type + "' unknown");
HADRON_ERROR(Argument, "object of type '" + type + "' unknown");
}
return func(name);

View File

@ -57,9 +57,7 @@ public:
virtual ~GeneticScheduler(void) = default;
// access
const Gene & getMinSchedule(void);
V getMinValue(void);
// reset population
void initPopulation(void);
int getMinValue(void);
// breed a new generation
void nextGeneration(void);
// heuristic benchmarks
@ -78,6 +76,8 @@ public:
return out;
}
private:
// evolution steps
void initPopulation(void);
void doCrossover(void);
void doMutation(void);
// genetic operators
@ -116,7 +116,7 @@ GeneticScheduler<V, T>::getMinSchedule(void)
}
template <typename V, typename T>
V GeneticScheduler<V, T>::getMinValue(void)
int GeneticScheduler<V, T>::getMinValue(void)
{
return population_.begin()->first;
}
@ -130,7 +130,7 @@ void GeneticScheduler<V, T>::nextGeneration(void)
{
initPopulation();
}
//LOG(Debug) << "Starting population:\n" << *this << std::endl;
LOG(Debug) << "Starting population:\n" << *this << std::endl;
// random mutations
//PARALLEL_FOR_LOOP
@ -138,7 +138,7 @@ void GeneticScheduler<V, T>::nextGeneration(void)
{
doMutation();
}
//LOG(Debug) << "After mutations:\n" << *this << std::endl;
LOG(Debug) << "After mutations:\n" << *this << std::endl;
// mating
//PARALLEL_FOR_LOOP
@ -146,14 +146,14 @@ void GeneticScheduler<V, T>::nextGeneration(void)
{
doCrossover();
}
//LOG(Debug) << "After mating:\n" << *this << std::endl;
LOG(Debug) << "After mating:\n" << *this << std::endl;
// grim reaper
auto it = population_.begin();
std::advance(it, par_.popSize);
population_.erase(it, population_.end());
//LOG(Debug) << "After grim reaper:\n" << *this << std::endl;
LOG(Debug) << "After grim reaper:\n" << *this << std::endl;
}
// evolution steps /////////////////////////////////////////////////////////////

View File

@ -37,38 +37,20 @@ HadronsLogger Hadrons::HadronsLogWarning(1,"Warning");
HadronsLogger Hadrons::HadronsLogMessage(1,"Message");
HadronsLogger Hadrons::HadronsLogIterative(1,"Iterative");
HadronsLogger Hadrons::HadronsLogDebug(1,"Debug");
HadronsLogger Hadrons::HadronsLogIRL(1,"IRL");
void Hadrons::initLogger(void)
{
auto w = std::string("Hadrons").length();
int cw = 8;
auto w = std::string("Hadrons").length();
GridLogError.setTopWidth(w);
GridLogWarning.setTopWidth(w);
GridLogMessage.setTopWidth(w);
GridLogIterative.setTopWidth(w);
GridLogDebug.setTopWidth(w);
GridLogIRL.setTopWidth(w);
GridLogError.setChanWidth(cw);
GridLogWarning.setChanWidth(cw);
GridLogMessage.setChanWidth(cw);
GridLogIterative.setChanWidth(cw);
GridLogDebug.setChanWidth(cw);
GridLogIRL.setChanWidth(cw);
HadronsLogError.Active(true);
HadronsLogWarning.Active(true);
HadronsLogError.Active(GridLogError.isActive());
HadronsLogWarning.Active(GridLogWarning.isActive());
HadronsLogMessage.Active(GridLogMessage.isActive());
HadronsLogIterative.Active(GridLogIterative.isActive());
HadronsLogDebug.Active(GridLogDebug.isActive());
HadronsLogIRL.Active(GridLogIRL.isActive());
HadronsLogError.setChanWidth(cw);
HadronsLogWarning.setChanWidth(cw);
HadronsLogMessage.setChanWidth(cw);
HadronsLogIterative.setChanWidth(cw);
HadronsLogDebug.setChanWidth(cw);
HadronsLogIRL.setChanWidth(cw);
}
// type utilities //////////////////////////////////////////////////////////////
@ -92,84 +74,3 @@ const std::string Hadrons::resultFileExt = "h5";
#else
const std::string Hadrons::resultFileExt = "xml";
#endif
// recursive mkdir /////////////////////////////////////////////////////////////
int Hadrons::mkdir(const std::string dirName)
{
if (!dirName.empty() and access(dirName.c_str(), R_OK|W_OK|X_OK))
{
mode_t mode755;
char tmp[MAX_PATH_LENGTH];
char *p = NULL;
size_t len;
mode755 = S_IRWXU|S_IRGRP|S_IXGRP|S_IROTH|S_IXOTH;
snprintf(tmp, sizeof(tmp), "%s", dirName.c_str());
len = strlen(tmp);
if(tmp[len - 1] == '/')
{
tmp[len - 1] = 0;
}
for(p = tmp + 1; *p; p++)
{
if(*p == '/')
{
*p = 0;
::mkdir(tmp, mode755);
*p = '/';
}
}
return ::mkdir(tmp, mode755);
}
else
{
return 0;
}
}
std::string Hadrons::basename(const std::string &s)
{
constexpr char sep = '/';
size_t i = s.rfind(sep, s.length());
if (i != std::string::npos)
{
return s.substr(i+1, s.length() - i);
}
else
{
return s;
}
}
std::string Hadrons::dirname(const std::string &s)
{
constexpr char sep = '/';
size_t i = s.rfind(sep, s.length());
if (i != std::string::npos)
{
return s.substr(0, i);
}
else
{
return "";
}
}
void Hadrons::makeFileDir(const std::string filename, GridBase *g)
{
if (g->IsBoss())
{
std::string dir = dirname(filename);
int status = mkdir(dir);
if (status)
{
HADRONS_ERROR(Io, "cannot create directory '" + dir
+ "' ( " + std::strerror(errno) + ")");
}
}
}

View File

@ -39,40 +39,30 @@ See the full license in the file "LICENSE" in the top level distribution directo
#define SITE_SIZE_TYPE size_t
#endif
#ifndef DEFAULT_ASCII_PREC
#define DEFAULT_ASCII_PREC 16
#endif
#define BEGIN_HADRONS_NAMESPACE \
namespace Grid {\
using namespace QCD;\
namespace Hadrons {\
using Grid::operator<<;
#define END_HADRONS_NAMESPACE }}
#define BEGIN_MODULE_NAMESPACE(name)\
namespace name {\
using Grid::operator<<;
#define END_MODULE_NAMESPACE }
/* the 'using Grid::operator<<;' statement prevents a very nasty compilation
* error with GCC 5 (clang & GCC 6 compile fine without it).
*/
#define BEGIN_HADRONS_NAMESPACE \
namespace Grid {\
using namespace QCD;\
namespace Hadrons {\
using Grid::operator<<;\
using Grid::operator>>;
#define END_HADRONS_NAMESPACE }}
#define BEGIN_MODULE_NAMESPACE(name)\
namespace name {\
using Grid::operator<<;\
using Grid::operator>>;
#define END_MODULE_NAMESPACE }
#ifndef FIMPL
#define FIMPL WilsonImplR
#endif
#ifndef ZFIMPL
#define ZFIMPL ZWilsonImplR
#endif
#ifndef SIMPL
#define SIMPL ScalarImplCR
#endif
#ifndef GIMPL
#define GIMPL PeriodicGimplR
#define GIMPL GimplTypesR
#endif
BEGIN_HADRONS_NAMESPACE
@ -93,15 +83,17 @@ typedef typename SImpl::Field ScalarField##suffix;\
typedef typename SImpl::Field PropagatorField##suffix;
#define SOLVER_TYPE_ALIASES(FImpl, suffix)\
typedef Solver<FImpl> Solver##suffix;
typedef std::function<void(FermionField##suffix &,\
const FermionField##suffix &)> SolverFn##suffix;
#define SINK_TYPE_ALIASES(suffix)\
typedef std::function<SlicedPropagator##suffix\
(const PropagatorField##suffix &)> SinkFn##suffix;
#define FG_TYPE_ALIASES(FImpl, suffix)\
#define FGS_TYPE_ALIASES(FImpl, suffix)\
FERM_TYPE_ALIASES(FImpl, suffix)\
GAUGE_TYPE_ALIASES(FImpl, suffix)
GAUGE_TYPE_ALIASES(FImpl, suffix)\
SOLVER_TYPE_ALIASES(FImpl, suffix)
// logger
class HadronsLogger: public Logger
@ -112,14 +104,13 @@ public:
};
#define LOG(channel) std::cout << HadronsLog##channel
#define HADRONS_DEBUG_VAR(var) LOG(Debug) << #var << "= " << (var) << std::endl;
#define DEBUG_VAR(var) LOG(Debug) << #var << "= " << (var) << std::endl;
extern HadronsLogger HadronsLogError;
extern HadronsLogger HadronsLogWarning;
extern HadronsLogger HadronsLogMessage;
extern HadronsLogger HadronsLogIterative;
extern HadronsLogger HadronsLogDebug;
extern HadronsLogger HadronsLogIRL;
void initLogger(void);
@ -189,28 +180,6 @@ typedef XmlWriter ResultWriter;
#define RESULT_FILE_NAME(name) \
name + "." + std::to_string(vm().getTrajectory()) + "." + resultFileExt
// recursive mkdir
#define MAX_PATH_LENGTH 512u
int mkdir(const std::string dirName);
std::string basename(const std::string &s);
std::string dirname(const std::string &s);
void makeFileDir(const std::string filename, GridBase *g);
// default Schur convention
#ifndef HADRONS_DEFAULT_SCHUR
#define HADRONS_DEFAULT_SCHUR DiagTwo
#endif
#define _HADRONS_SCHUR_OP_(conv) Schur##conv##Operator
#define HADRONS_SCHUR_OP(conv) _HADRONS_SCHUR_OP_(conv)
#define HADRONS_DEFAULT_SCHUR_OP HADRONS_SCHUR_OP(HADRONS_DEFAULT_SCHUR)
#define _HADRONS_SCHUR_SOLVE_(conv) SchurRedBlack##conv##Solve
#define HADRONS_SCHUR_SOLVE(conv) _HADRONS_SCHUR_SOLVE_(conv)
#define HADRONS_DEFAULT_SCHUR_SOLVE HADRONS_SCHUR_SOLVE(HADRONS_DEFAULT_SCHUR)
// stringify macro
#define _HADRONS_STR(x) #x
#define HADRONS_STR(x) _HADRONS_STR(x)
END_HADRONS_NAMESPACE
#include <Grid/Hadrons/Exceptions.hpp>

View File

@ -184,7 +184,7 @@ void Graph<T>::removeVertex(const T &value)
}
else
{
HADRONS_ERROR(Range, "vertex does not exists");
HADRON_ERROR(Range, "vertex does not exists");
}
// remove all edges containing the vertex
@ -213,7 +213,7 @@ void Graph<T>::removeEdge(const Edge &e)
}
else
{
HADRONS_ERROR(Range, "edge does not exists");
HADRON_ERROR(Range, "edge does not exists");
}
}
@ -259,7 +259,7 @@ void Graph<T>::mark(const T &value, const bool doMark)
}
else
{
HADRONS_ERROR(Range, "vertex does not exists");
HADRON_ERROR(Range, "vertex does not exists");
}
}
@ -297,7 +297,7 @@ bool Graph<T>::isMarked(const T &value) const
}
else
{
HADRONS_ERROR(Range, "vertex does not exists");
HADRON_ERROR(Range, "vertex does not exists");
return false;
}
@ -543,7 +543,7 @@ std::vector<T> Graph<T>::topoSort(void)
{
if (tmpMarked.at(v))
{
HADRONS_ERROR(Range, "cannot topologically sort a cyclic graph");
HADRON_ERROR(Range, "cannot topologically sort a cyclic graph");
}
if (!isMarked(v))
{
@ -602,7 +602,7 @@ std::vector<T> Graph<T>::topoSort(Gen &gen)
{
if (tmpMarked.at(v))
{
HADRONS_ERROR(Range, "cannot topologically sort a cyclic graph");
HADRON_ERROR(Range, "cannot topologically sort a cyclic graph");
}
if (!isMarked(v))
{

View File

@ -56,21 +56,14 @@ int main(int argc, char *argv[])
Grid_init(&argc, &argv);
// execution
try
Application application(parameterFileName);
application.parseParameterFile(parameterFileName);
if (!scheduleFileName.empty())
{
Application application(parameterFileName);
application.parseParameterFile(parameterFileName);
if (!scheduleFileName.empty())
{
application.loadSchedule(scheduleFileName);
}
application.run();
}
catch (const std::exception& e)
{
Exceptions::abort(e);
application.loadSchedule(scheduleFileName);
}
application.run();
// epilogue
LOG(Message) << "Grid is finalizing now" << std::endl;

View File

@ -2,7 +2,7 @@
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MAction/DWF.cc
Source file: extras/Hadrons/HadronsXmlSchedule.cc
Copyright (C) 2015-2018
@ -25,11 +25,41 @@ with this program; if not, write to the Free Software Foundation, Inc.,
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Hadrons/Modules/MAction/DWF.hpp>
#include <Grid/Hadrons/Application.hpp>
using namespace Grid;
using namespace QCD;
using namespace Hadrons;
using namespace MAction;
template class Grid::Hadrons::MAction::TDWF<FIMPL>;
int main(int argc, char *argv[])
{
// parse command line
std::string parameterFileName, scheduleFileName;
if (argc < 3)
{
std::cerr << "usage: " << argv[0] << " <parameter file> <schedule output> [Grid options]";
std::cerr << std::endl;
std::exit(EXIT_FAILURE);
}
parameterFileName = argv[1];
scheduleFileName = argv[2];
// initialization
Grid_init(&argc, &argv);
// execution
Application application;
application.parseParameterFile(parameterFileName);
application.schedule();
application.printSchedule();
application.saveSchedule(scheduleFileName);
// epilogue
LOG(Message) << "Grid is finalizing now" << std::endl;
Grid_finalize();
return EXIT_SUCCESS;
}

View File

@ -1,5 +1,5 @@
lib_LIBRARIES = libHadrons.a
bin_PROGRAMS = HadronsXmlRun
bin_PROGRAMS = HadronsXmlRun HadronsXmlSchedule
include modules.inc
@ -14,10 +14,7 @@ libHadrons_a_SOURCES = \
libHadrons_adir = $(pkgincludedir)/Hadrons
nobase_libHadrons_a_HEADERS = \
$(modules_hpp) \
AllToAllVectors.hpp \
AllToAllReduction.hpp \
Application.hpp \
EigenPack.hpp \
Environment.hpp \
Exceptions.hpp \
Factory.hpp \
@ -27,8 +24,10 @@ nobase_libHadrons_a_HEADERS = \
Module.hpp \
Modules.hpp \
ModuleFactory.hpp \
Solver.hpp \
VirtualMachine.hpp
HadronsXmlRun_SOURCES = HadronsXmlRun.cc
HadronsXmlRun_LDADD = libHadrons.a -lGrid
HadronsXmlSchedule_SOURCES = HadronsXmlSchedule.cc
HadronsXmlSchedule_LDADD = libHadrons.a -lGrid

View File

@ -49,7 +49,7 @@ std::string ModuleBase::getName(void) const
// get factory registration name if available
std::string ModuleBase::getRegisteredName(void)
{
HADRONS_ERROR(Definition, "module '" + getName() + "' has no registered type"
HADRON_ERROR(Definition, "module '" + getName() + "' has no registered type"
+ " in the factory");
}

View File

@ -35,7 +35,32 @@ See the full license in the file "LICENSE" in the top level distribution directo
BEGIN_HADRONS_NAMESPACE
// module registration macros
#define MODULE_REGISTER(mod, base, ns)\
#define MODULE_REGISTER(mod, base)\
class mod: public base\
{\
public:\
typedef base Base;\
using Base::Base;\
virtual std::string getRegisteredName(void)\
{\
return std::string(#mod);\
}\
};\
class mod##ModuleRegistrar\
{\
public:\
mod##ModuleRegistrar(void)\
{\
ModuleFactory &modFac = ModuleFactory::getInstance();\
modFac.registerBuilder(#mod, [&](const std::string name)\
{\
return std::unique_ptr<mod>(new mod(name));\
});\
}\
};\
static mod##ModuleRegistrar mod##ModuleRegistrarInstance;
#define MODULE_REGISTER_NS(mod, base, ns)\
class mod: public base\
{\
public:\
@ -60,19 +85,12 @@ public:\
};\
static ns##mod##ModuleRegistrar ns##mod##ModuleRegistrarInstance;
#define MODULE_REGISTER_TMP(mod, base, ns)\
extern template class base;\
MODULE_REGISTER(mod, ARG(base), ns);
#define ARG(...) __VA_ARGS__
#define MACRO_REDIRECT(arg1, arg2, arg3, macro, ...) macro
#define envGet(type, name)\
*env().template getObject<type>(name)
#define envGetDerived(base, type, name)\
*env().template getDerivedObject<base, type>(name)
#define envGetTmp(type, var)\
type &var = *env().template getObject<type>(getName() + "_tmp_" + #var)
@ -119,16 +137,6 @@ envTmp(type, name, Ls, env().getGrid(Ls))
#define envTmpLat(...)\
MACRO_REDIRECT(__VA_ARGS__, envTmpLat5, envTmpLat4)(__VA_ARGS__)
#define saveResult(ioStem, name, result)\
if (env().getGrid()->IsBoss() and !ioStem.empty())\
{\
makeFileDir(ioStem, env().getGrid());\
{\
ResultWriter _writer(RESULT_FILE_NAME(ioStem));\
write(_writer, name, result);\
}\
}
/******************************************************************************
* Module class *
******************************************************************************/
@ -154,8 +162,6 @@ public:
// parse parameters
virtual void parseParameters(XmlReader &reader, const std::string name) = 0;
virtual void saveParameters(XmlWriter &writer, const std::string name) = 0;
// parameter string
virtual std::string parString(void) const = 0;
// setup
virtual void setup(void) {};
virtual void execute(void) = 0;
@ -184,11 +190,9 @@ public:
// parse parameters
virtual void parseParameters(XmlReader &reader, const std::string name);
virtual void saveParameters(XmlWriter &writer, const std::string name);
// parameter string
virtual std::string parString(void) const;
// parameter access
const P & par(void) const;
void setPar(const P &par);
const P & par(void) const;
void setPar(const P &par);
private:
P par_;
};
@ -211,8 +215,6 @@ public:
push(writer, "options");
pop(writer);
};
// parameter string (empty)
virtual std::string parString(void) const {return "";};
};
/******************************************************************************
@ -235,16 +237,6 @@ void Module<P>::saveParameters(XmlWriter &writer, const std::string name)
write(writer, name, par_);
}
template <typename P>
std::string Module<P>::parString(void) const
{
XmlWriter writer("", "");
write(writer, par_.SerialisableClassName(), par_);
return writer.string();
}
template <typename P>
const P & Module<P>::par(void) const
{

View File

@ -1,62 +1,65 @@
#include <Grid/Hadrons/Modules/MScalarSUN/TrKinetic.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/StochFreeField.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/Grad.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TransProj.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/Div.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TrMag.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/EMT.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TwoPoint.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TrPhi.hpp>
#include <Grid/Hadrons/Modules/MScalar/FreeProp.hpp>
#include <Grid/Hadrons/Modules/MScalar/Scalar.hpp>
#include <Grid/Hadrons/Modules/MScalar/ScalarVP.hpp>
#include <Grid/Hadrons/Modules/MScalar/ChargedProp.hpp>
#include <Grid/Hadrons/Modules/MScalar/VPCounterTerms.hpp>
#include <Grid/Hadrons/Modules/MLoop/NoiseLoop.hpp>
#include <Grid/Hadrons/Modules/MIO/LoadEigenPack.hpp>
#include <Grid/Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp>
#include <Grid/Hadrons/Modules/MIO/LoadBinary.hpp>
#include <Grid/Hadrons/Modules/MIO/LoadNersc.hpp>
#include <Grid/Hadrons/Modules/MSink/Smear.hpp>
#include <Grid/Hadrons/Modules/MSink/Point.hpp>
#include <Grid/Hadrons/Modules/MFermion/FreeProp.hpp>
#include <Grid/Hadrons/Modules/MFermion/GaugeProp.hpp>
#include <Grid/Hadrons/Modules/MGauge/FundtoHirep.hpp>
#include <Grid/Hadrons/Modules/MGauge/Random.hpp>
#include <Grid/Hadrons/Modules/MGauge/StoutSmearing.hpp>
#include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
#include <Grid/Hadrons/Modules/MGauge/StochEm.hpp>
#include <Grid/Hadrons/Modules/MGauge/UnitEm.hpp>
#include <Grid/Hadrons/Modules/MUtilities/TestSeqGamma.hpp>
#include <Grid/Hadrons/Modules/MUtilities/TestSeqConserved.hpp>
#include <Grid/Hadrons/Modules/MSource/SeqConserved.hpp>
#include <Grid/Hadrons/Modules/MSource/Z2.hpp>
#include <Grid/Hadrons/Modules/MSource/Wall.hpp>
#include <Grid/Hadrons/Modules/MSource/SeqGamma.hpp>
#include <Grid/Hadrons/Modules/MSource/Point.hpp>
#include <Grid/Hadrons/Modules/MContraction/MesonFieldGamma.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules.hpp
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
Author: Lanny91 <andrew.lawson@gmail.com>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Hadrons/Modules/MContraction/Baryon.hpp>
#include <Grid/Hadrons/Modules/MContraction/A2APionField.hpp>
#include <Grid/Hadrons/Modules/MContraction/Meson.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>
#include <Grid/Hadrons/Modules/MContraction/DiscLoop.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp>
#include <Grid/Hadrons/Modules/MContraction/Gamma3pt.hpp>
#include <Grid/Hadrons/Modules/MContraction/DiscLoop.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>
#include <Grid/Hadrons/Modules/MContraction/A2AMeson.hpp>
#include <Grid/Hadrons/Modules/MContraction/WardIdentity.hpp>
#include <Grid/Hadrons/Modules/MContraction/A2AMesonField.hpp>
#include <Grid/Hadrons/Modules/MAction/WilsonClover.hpp>
#include <Grid/Hadrons/Modules/MAction/ScaledDWF.hpp>
#include <Grid/Hadrons/Modules/MAction/MobiusDWF.hpp>
#include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
#include <Grid/Hadrons/Modules/MAction/DWF.hpp>
#include <Grid/Hadrons/Modules/MAction/ZMobiusDWF.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>
#include <Grid/Hadrons/Modules/MFermion/GaugeProp.hpp>
#include <Grid/Hadrons/Modules/MSource/SeqGamma.hpp>
#include <Grid/Hadrons/Modules/MSource/Point.hpp>
#include <Grid/Hadrons/Modules/MSource/Wall.hpp>
#include <Grid/Hadrons/Modules/MSource/Z2.hpp>
#include <Grid/Hadrons/Modules/MSource/SeqConserved.hpp>
#include <Grid/Hadrons/Modules/MSink/Smear.hpp>
#include <Grid/Hadrons/Modules/MSink/Point.hpp>
#include <Grid/Hadrons/Modules/MSolver/RBPrecCG.hpp>
#include <Grid/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp>
#include <Grid/Hadrons/Modules/MSolver/A2AVectors.hpp>
#include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
#include <Grid/Hadrons/Modules/MGauge/Random.hpp>
#include <Grid/Hadrons/Modules/MGauge/StochEm.hpp>
#include <Grid/Hadrons/Modules/MGauge/FundtoHirep.hpp>
#include <Grid/Hadrons/Modules/MUtilities/TestSeqGamma.hpp>
#include <Grid/Hadrons/Modules/MUtilities/TestSeqConserved.hpp>
#include <Grid/Hadrons/Modules/MLoop/NoiseLoop.hpp>
#include <Grid/Hadrons/Modules/MScalar/FreeProp.hpp>
#include <Grid/Hadrons/Modules/MScalar/Scalar.hpp>
#include <Grid/Hadrons/Modules/MScalar/ChargedProp.hpp>
#include <Grid/Hadrons/Modules/MAction/DWF.hpp>
#include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
#include <Grid/Hadrons/Modules/MAction/WilsonClover.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/Div.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TrMag.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TwoPoint.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TrPhi.hpp>
#include <Grid/Hadrons/Modules/MIO/LoadNersc.hpp>
#include <Grid/Hadrons/Modules/MIO/LoadBinary.hpp>

View File

@ -56,12 +56,12 @@ template <typename FImpl>
class TDWF: public Module<DWFPar>
{
public:
FG_TYPE_ALIASES(FImpl,);
FGS_TYPE_ALIASES(FImpl,);
public:
// constructor
TDWF(const std::string name);
// destructor
virtual ~TDWF(void) {};
virtual ~TDWF(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
@ -72,8 +72,7 @@ protected:
virtual void execute(void);
};
extern template class TDWF<FIMPL>;
MODULE_REGISTER_TMP(DWF, TDWF<FIMPL>, MAction);
MODULE_REGISTER_NS(DWF, TDWF<FIMPL>, MAction);
/******************************************************************************
* DWF template implementation *

View File

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

View File

@ -1,109 +0,0 @@
#ifndef Hadrons_MAction_MobiusDWF_hpp_
#define Hadrons_MAction_MobiusDWF_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* Mobius domain-wall fermion action *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MAction)
class MobiusDWFPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(MobiusDWFPar,
std::string , gauge,
unsigned int, Ls,
double , mass,
double , M5,
double , b,
double , c,
std::string , boundary);
};
template <typename FImpl>
class TMobiusDWF: public Module<MobiusDWFPar>
{
public:
FG_TYPE_ALIASES(FImpl,);
public:
// constructor
TMobiusDWF(const std::string name);
// destructor
virtual ~TMobiusDWF(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(MobiusDWF, TMobiusDWF<FIMPL>, MAction);
/******************************************************************************
* TMobiusDWF implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TMobiusDWF<FImpl>::TMobiusDWF(const std::string name)
: Module<MobiusDWFPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TMobiusDWF<FImpl>::getInput(void)
{
std::vector<std::string> in = {par().gauge};
return in;
}
template <typename FImpl>
std::vector<std::string> TMobiusDWF<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TMobiusDWF<FImpl>::setup(void)
{
LOG(Message) << "Setting up Mobius domain wall fermion matrix with m= "
<< par().mass << ", M5= " << par().M5 << ", Ls= " << par().Ls
<< ", b= " << par().b << ", c= " << par().c
<< " using gauge field '" << par().gauge << "'"
<< std::endl;
LOG(Message) << "Fermion boundary conditions: " << par().boundary
<< std::endl;
env().createGrid(par().Ls);
auto &U = envGet(LatticeGaugeField, par().gauge);
auto &g4 = *env().getGrid();
auto &grb4 = *env().getRbGrid();
auto &g5 = *env().getGrid(par().Ls);
auto &grb5 = *env().getRbGrid(par().Ls);
std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
typename MobiusFermion<FImpl>::ImplParams implParams(boundary);
envCreateDerived(FMat, MobiusFermion<FImpl>, getName(), par().Ls, U, g5,
grb5, g4, grb4, par().mass, par().M5, par().b, par().c,
implParams);
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TMobiusDWF<FImpl>::execute(void)
{}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MAction_MobiusDWF_hpp_

View File

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

View File

@ -1,108 +0,0 @@
#ifndef Hadrons_MAction_ScaledDWF_hpp_
#define Hadrons_MAction_ScaledDWF_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* Scaled domain wall fermion *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MAction)
class ScaledDWFPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(ScaledDWFPar,
std::string , gauge,
unsigned int, Ls,
double , mass,
double , M5,
double , scale,
std::string , boundary);
};
template <typename FImpl>
class TScaledDWF: public Module<ScaledDWFPar>
{
public:
FG_TYPE_ALIASES(FImpl,);
public:
// constructor
TScaledDWF(const std::string name);
// destructor
virtual ~TScaledDWF(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(ScaledDWF, TScaledDWF<FIMPL>, MAction);
/******************************************************************************
* TScaledDWF implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TScaledDWF<FImpl>::TScaledDWF(const std::string name)
: Module<ScaledDWFPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TScaledDWF<FImpl>::getInput(void)
{
std::vector<std::string> in = {par().gauge};
return in;
}
template <typename FImpl>
std::vector<std::string> TScaledDWF<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TScaledDWF<FImpl>::setup(void)
{
LOG(Message) << "Setting up scaled domain wall fermion matrix with m= "
<< par().mass << ", M5= " << par().M5 << ", Ls= " << par().Ls
<< ", scale= " << par().scale
<< " using gauge field '" << par().gauge << "'"
<< std::endl;
LOG(Message) << "Fermion boundary conditions: " << par().boundary
<< std::endl;
env().createGrid(par().Ls);
auto &U = envGet(LatticeGaugeField, par().gauge);
auto &g4 = *env().getGrid();
auto &grb4 = *env().getRbGrid();
auto &g5 = *env().getGrid(par().Ls);
auto &grb5 = *env().getRbGrid(par().Ls);
std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
typename MobiusFermion<FImpl>::ImplParams implParams(boundary);
envCreateDerived(FMat, ScaledShamirFermion<FImpl>, getName(), par().Ls, U, g5,
grb5, g4, grb4, par().mass, par().M5, par().scale,
implParams);
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TScaledDWF<FImpl>::execute(void)
{}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MAction_ScaledDWF_hpp_

View File

@ -1,35 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MAction/Wilson.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 <Grid/Hadrons/Modules/MAction/Wilson.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MAction;
template class Grid::Hadrons::MAction::TWilson<FIMPL>;

View File

@ -54,12 +54,12 @@ template <typename FImpl>
class TWilson: public Module<WilsonPar>
{
public:
FG_TYPE_ALIASES(FImpl,);
FGS_TYPE_ALIASES(FImpl,);
public:
// constructor
TWilson(const std::string name);
// destructor
virtual ~TWilson(void) {};
virtual ~TWilson(void) = default;
// dependencies/products
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
@ -70,7 +70,7 @@ protected:
virtual void execute(void);
};
MODULE_REGISTER_TMP(Wilson, TWilson<FIMPL>, MAction);
MODULE_REGISTER_NS(Wilson, TWilson<FIMPL>, MAction);
/******************************************************************************
* TWilson template implementation *
@ -102,7 +102,7 @@ std::vector<std::string> TWilson<FImpl>::getOutput(void)
template <typename FImpl>
void TWilson<FImpl>::setup(void)
{
LOG(Message) << "Setting up Wilson fermion matrix with m= " << par().mass
LOG(Message) << "Setting up TWilson fermion matrix with m= " << par().mass
<< " using gauge field '" << par().gauge << "'" << std::endl;
LOG(Message) << "Fermion boundary conditions: " << par().boundary
<< std::endl;

View File

@ -1,35 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MAction/WilsonClover.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 <Grid/Hadrons/Modules/MAction/WilsonClover.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MAction;
template class Grid::Hadrons::MAction::TWilsonClover<FIMPL>;

View File

@ -2,13 +2,12 @@
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MAction/WilsonClover.hpp
Source file: extras/Hadrons/Modules/MAction/Wilson.hpp
Copyright (C) 2015-2018
Copyright (C) 2015
Copyright (C) 2016
Author: Antonin Portelli <antonin.portelli@me.com>
Author: Guido Cossu <guido.cossu@ed.ac.uk>
Author: pretidav <david.preti@csic.es>
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
@ -38,7 +37,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* Wilson clover quark action *
* TWilson quark action *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MAction)
@ -59,12 +58,12 @@ template <typename FImpl>
class TWilsonClover: public Module<WilsonCloverPar>
{
public:
FG_TYPE_ALIASES(FImpl,);
FGS_TYPE_ALIASES(FImpl,);
public:
// constructor
TWilsonClover(const std::string name);
// destructor
virtual ~TWilsonClover(void) {};
virtual ~TWilsonClover(void) = default;
// dependencies/products
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
@ -74,7 +73,7 @@ public:
virtual void execute(void);
};
MODULE_REGISTER_TMP(WilsonClover, TWilsonClover<FIMPL>, MAction);
MODULE_REGISTER_NS(WilsonClover, TWilsonClover<FIMPL>, MAction);
/******************************************************************************
* TWilsonClover template implementation *
@ -106,7 +105,13 @@ std::vector<std::string> TWilsonClover<FImpl>::getOutput(void)
template <typename FImpl>
void TWilsonClover<FImpl>::setup(void)
{
LOG(Message) << "Setting up Wilson clover fermion matrix with m= " << par().mass
//unsigned int size;
// size = 2*env().template lattice4dSize<typename FImpl::DoubledGaugeField>();
// env().registerObject(getName(), size);
LOG(Message) << "Setting up TWilsonClover fermion matrix with m= " << par().mass
<< " using gauge field '" << par().gauge << "'" << std::endl;
LOG(Message) << "Fermion boundary conditions: " << par().boundary
<< std::endl;
@ -123,12 +128,23 @@ void TWilsonClover<FImpl>::setup(void)
par().csw_t,
par().clover_anisotropy,
implParams);
//FMat *fMatPt = new WilsonCloverFermion<FImpl>(U, grid, gridRb, par().mass,
// par().csw_r,
// par().csw_t,
// par().clover_anisotropy,
// implParams);
//env().setObject(getName(), fMatPt);
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TWilsonClover<FImpl>::execute()
{}
{
}
END_MODULE_NAMESPACE

View File

@ -1,35 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MAction/ZMobiusDWF.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 <Grid/Hadrons/Modules/MAction/ZMobiusDWF.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MAction;
template class Grid::Hadrons::MAction::TZMobiusDWF<ZFIMPL>;

View File

@ -1,143 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_MAction_ZMobiusDWF_hpp_
#define Hadrons_MAction_ZMobiusDWF_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* z-Mobius domain-wall fermion action *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MAction)
class ZMobiusDWFPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(ZMobiusDWFPar,
std::string , gauge,
unsigned int , Ls,
double , mass,
double , M5,
double , b,
double , c,
std::vector<std::complex<double>>, omega,
std::string , boundary);
};
template <typename FImpl>
class TZMobiusDWF: public Module<ZMobiusDWFPar>
{
public:
FG_TYPE_ALIASES(FImpl,);
public:
// constructor
TZMobiusDWF(const std::string name);
// destructor
virtual ~TZMobiusDWF(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(ZMobiusDWF, TZMobiusDWF<ZFIMPL>, MAction);
/******************************************************************************
* TZMobiusDWF implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TZMobiusDWF<FImpl>::TZMobiusDWF(const std::string name)
: Module<ZMobiusDWFPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TZMobiusDWF<FImpl>::getInput(void)
{
std::vector<std::string> in = {par().gauge};
return in;
}
template <typename FImpl>
std::vector<std::string> TZMobiusDWF<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TZMobiusDWF<FImpl>::setup(void)
{
LOG(Message) << "Setting up z-Mobius domain wall fermion matrix with m= "
<< par().mass << ", M5= " << par().M5 << ", Ls= " << par().Ls
<< ", b= " << par().b << ", c= " << par().c
<< " using gauge field '" << par().gauge << "'"
<< std::endl;
LOG(Message) << "Omegas: " << std::endl;
for (unsigned int i = 0; i < par().omega.size(); ++i)
{
LOG(Message) << " omega[" << i << "]= " << par().omega[i] << std::endl;
}
LOG(Message) << "Fermion boundary conditions: " << par().boundary
<< std::endl;
env().createGrid(par().Ls);
auto &U = envGet(LatticeGaugeField, par().gauge);
auto &g4 = *env().getGrid();
auto &grb4 = *env().getRbGrid();
auto &g5 = *env().getGrid(par().Ls);
auto &grb5 = *env().getRbGrid(par().Ls);
auto omega = par().omega;
std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
typename ZMobiusFermion<FImpl>::ImplParams implParams(boundary);
envCreateDerived(FMat, ZMobiusFermion<FImpl>, getName(), par().Ls, U, g5,
grb5, g4, grb4, par().mass, par().M5, omega,
par().b, par().c, implParams);
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TZMobiusDWF<FImpl>::execute(void)
{}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MAction_ZMobiusDWF_hpp_

View File

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

View File

@ -1,207 +0,0 @@
#ifndef Hadrons_MContraction_A2AMeson_hpp_
#define Hadrons_MContraction_A2AMeson_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
#include <Grid/Hadrons/AllToAllVectors.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* A2AMeson *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MContraction)
typedef std::pair<Gamma::Algebra, Gamma::Algebra> GammaPair;
class A2AMesonPar : Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(A2AMesonPar,
int, Nl,
int, N,
std::string, A2A1,
std::string, A2A2,
std::string, gammas,
std::string, output);
};
template <typename FImpl>
class TA2AMeson : public Module<A2AMesonPar>
{
public:
FERM_TYPE_ALIASES(FImpl, );
SOLVER_TYPE_ALIASES(FImpl, );
typedef A2AModesSchurDiagTwo<typename FImpl::FermionField, FMat, Solver> A2ABase;
class Result : Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
Gamma::Algebra, gamma_snk,
Gamma::Algebra, gamma_src,
std::vector<Complex>, corr);
};
public:
// constructor
TA2AMeson(const std::string name);
// destructor
virtual ~TA2AMeson(void){};
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
virtual void parseGammaString(std::vector<GammaPair> &gammaList);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
};
MODULE_REGISTER(A2AMeson, ARG(TA2AMeson<FIMPL>), MContraction);
MODULE_REGISTER(ZA2AMeson, ARG(TA2AMeson<ZFIMPL>), MContraction);
/******************************************************************************
* TA2AMeson implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TA2AMeson<FImpl>::TA2AMeson(const std::string name)
: Module<A2AMesonPar>(name)
{
}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TA2AMeson<FImpl>::getInput(void)
{
std::vector<std::string> in = {par().A2A1 + "_class", par().A2A2 + "_class"};
in.push_back(par().A2A1 + "_w_high_4d");
in.push_back(par().A2A2 + "_v_high_4d");
return in;
}
template <typename FImpl>
std::vector<std::string> TA2AMeson<FImpl>::getOutput(void)
{
std::vector<std::string> out = {};
return out;
}
template <typename FImpl>
void TA2AMeson<FImpl>::parseGammaString(std::vector<GammaPair> &gammaList)
{
gammaList.clear();
// Parse individual contractions from input string.
gammaList = strToVec<GammaPair>(par().gammas);
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TA2AMeson<FImpl>::setup(void)
{
int nt = env().getDim(Tp);
int N = par().N;
int Ls_ = env().getObjectLs(par().A2A1 + "_class");
envTmp(std::vector<FermionField>, "w1", 1, N, FermionField(env().getGrid(1)));
envTmp(std::vector<FermionField>, "v1", 1, N, FermionField(env().getGrid(1)));
envTmpLat(FermionField, "tmpv_5d", Ls_);
envTmpLat(FermionField, "tmpw_5d", Ls_);
envTmp(std::vector<ComplexD>, "MF_x", 1, nt);
envTmp(std::vector<ComplexD>, "MF_y", 1, nt);
envTmp(std::vector<ComplexD>, "tmp", 1, nt);
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TA2AMeson<FImpl>::execute(void)
{
LOG(Message) << "Computing A2A meson contractions" << std::endl;
Result result;
Gamma g5(Gamma::Algebra::Gamma5);
std::vector<GammaPair> gammaList;
int nt = env().getDim(Tp);
parseGammaString(gammaList);
result.gamma_snk = gammaList[0].first;
result.gamma_src = gammaList[0].second;
result.corr.resize(nt);
int Nl = par().Nl;
int N = par().N;
LOG(Message) << "N for A2A cont: " << N << std::endl;
envGetTmp(std::vector<ComplexD>, MF_x);
envGetTmp(std::vector<ComplexD>, MF_y);
envGetTmp(std::vector<ComplexD>, tmp);
for (unsigned int t = 0; t < nt; ++t)
{
tmp[t] = TensorRemove(MF_x[t] * MF_y[t] * 0.0);
}
Gamma gSnk(gammaList[0].first);
Gamma gSrc(gammaList[0].second);
auto &a2a1_fn = envGet(A2ABase, par().A2A1 + "_class");
envGetTmp(std::vector<FermionField>, w1);
envGetTmp(std::vector<FermionField>, v1);
envGetTmp(FermionField, tmpv_5d);
envGetTmp(FermionField, tmpw_5d);
LOG(Message) << "Finding v and w vectors for N = " << N << std::endl;
for (int i = 0; i < N; i++)
{
a2a1_fn.return_v(i, tmpv_5d, v1[i]);
a2a1_fn.return_w(i, tmpw_5d, w1[i]);
}
LOG(Message) << "Found v and w vectors for N = " << N << std::endl;
for (unsigned int i = 0; i < N; i++)
{
v1[i] = gSnk * v1[i];
}
int ty;
for (unsigned int i = 0; i < N; i++)
{
for (unsigned int j = 0; j < N; j++)
{
mySliceInnerProductVector(MF_x, w1[i], v1[j], Tp);
mySliceInnerProductVector(MF_y, w1[j], v1[i], Tp);
for (unsigned int t = 0; t < nt; ++t)
{
for (unsigned int tx = 0; tx < nt; tx++)
{
ty = (tx + t) % nt;
tmp[t] += TensorRemove((MF_x[tx]) * (MF_y[ty]));
}
}
}
if (i % 10 == 0)
{
LOG(Message) << "MF for i = " << i << " of " << N << std::endl;
}
}
double NTinv = 1.0 / static_cast<double>(nt);
for (unsigned int t = 0; t < nt; ++t)
{
result.corr[t] = NTinv * tmp[t];
}
saveResult(par().output, "meson", result);
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MContraction_A2AMeson_hpp_

View File

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

View File

@ -1,279 +0,0 @@
#ifndef Hadrons_MContraction_A2AMesonField_hpp_
#define Hadrons_MContraction_A2AMesonField_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
#include <Grid/Hadrons/AllToAllVectors.hpp>
#include <Grid/Hadrons/Modules/MContraction/A2Autils.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* A2AMesonField *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MContraction)
typedef std::pair<Gamma::Algebra, Gamma::Algebra> GammaPair;
class A2AMesonFieldPar : Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(A2AMesonFieldPar,
int, cacheBlock,
int, schurBlock,
int, Nmom,
int, N,
int, Nl,
std::string, A2A,
std::string, output);
};
template <typename FImpl>
class TA2AMesonField : public Module<A2AMesonFieldPar>
{
public:
FERM_TYPE_ALIASES(FImpl, );
SOLVER_TYPE_ALIASES(FImpl, );
typedef A2AModesSchurDiagTwo<typename FImpl::FermionField, FMat, Solver> A2ABase;
public:
// constructor
TA2AMesonField(const std::string name);
// destructor
virtual ~TA2AMesonField(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(A2AMesonField, ARG(TA2AMesonField<FIMPL>), MContraction);
MODULE_REGISTER(ZA2AMesonField, ARG(TA2AMesonField<ZFIMPL>), MContraction);
/******************************************************************************
* TA2AMesonField implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TA2AMesonField<FImpl>::TA2AMesonField(const std::string name)
: Module<A2AMesonFieldPar>(name)
{
}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TA2AMesonField<FImpl>::getInput(void)
{
std::vector<std::string> in = {par().A2A + "_class"};
in.push_back(par().A2A + "_w_high_4d");
in.push_back(par().A2A + "_v_high_4d");
return in;
}
template <typename FImpl>
std::vector<std::string> TA2AMesonField<FImpl>::getOutput(void)
{
std::vector<std::string> out = {};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TA2AMesonField<FImpl>::setup(void)
{
auto &a2a = envGet(A2ABase, par().A2A + "_class");
int nt = env().getDim(Tp);
int Nl = par().Nl;
int N = par().N;
int Ls_ = env().getObjectLs(par().A2A + "_class");
// Four D fields
envTmp(std::vector<FermionField>, "w", 1, par().schurBlock, FermionField(env().getGrid(1)));
envTmp(std::vector<FermionField>, "v", 1, par().schurBlock, FermionField(env().getGrid(1)));
// 5D tmp
envTmpLat(FermionField, "tmp_5d", Ls_);
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TA2AMesonField<FImpl>::execute(void)
{
LOG(Message) << "Computing A2A meson field" << std::endl;
auto &a2a = envGet(A2ABase, par().A2A + "_class");
// 2+6+4+4 = 16 gammas
// Ordering defined here
std::vector<Gamma::Algebra> gammas ( {
Gamma::Algebra::Gamma5,
Gamma::Algebra::Identity,
Gamma::Algebra::GammaX,
Gamma::Algebra::GammaY,
Gamma::Algebra::GammaZ,
Gamma::Algebra::GammaT,
Gamma::Algebra::GammaXGamma5,
Gamma::Algebra::GammaYGamma5,
Gamma::Algebra::GammaZGamma5,
Gamma::Algebra::GammaTGamma5,
Gamma::Algebra::SigmaXY,
Gamma::Algebra::SigmaXZ,
Gamma::Algebra::SigmaXT,
Gamma::Algebra::SigmaYZ,
Gamma::Algebra::SigmaYT,
Gamma::Algebra::SigmaZT
});
///////////////////////////////////////////////
// Square assumption for now Nl = Nr = N
///////////////////////////////////////////////
int nt = env().getDim(Tp);
int nx = env().getDim(Xp);
int ny = env().getDim(Yp);
int nz = env().getDim(Zp);
int N = par().N;
int Nl = par().Nl;
int ngamma = gammas.size();
int schurBlock = par().schurBlock;
int cacheBlock = par().cacheBlock;
int nmom = par().Nmom;
///////////////////////////////////////////////
// Momentum setup
///////////////////////////////////////////////
GridBase *grid = env().getGrid(1);
std::vector<LatticeComplex> phases(nmom,grid);
for(int m=0;m<nmom;m++){
phases[m] = Complex(1.0); // All zero momentum for now
}
Eigen::Tensor<ComplexD,5> mesonField (nmom,ngamma,nt,N,N);
LOG(Message) << "N = Nh+Nl for A2A MesonField is " << N << std::endl;
envGetTmp(std::vector<FermionField>, w);
envGetTmp(std::vector<FermionField>, v);
envGetTmp(FermionField, tmp_5d);
LOG(Message) << "Finding v and w vectors for N = " << N << std::endl;
//////////////////////////////////////////////////////////////////////////
// i,j is first loop over SchurBlock factors reusing 5D matrices
// ii,jj is second loop over cacheBlock factors for high perf contractoin
// iii,jjj are loops within cacheBlock
// Total index is sum of these i+ii+iii etc...
//////////////////////////////////////////////////////////////////////////
double flops = 0.0;
double bytes = 0.0;
double vol = nx*ny*nz*nt;
double t_schur=0;
double t_contr=0;
double t_int_0=0;
double t_int_1=0;
double t_int_2=0;
double t_int_3=0;
double t0 = usecond();
int N_i = N;
int N_j = N;
for(int i=0;i<N_i;i+=schurBlock){ //loop over SchurBlocking to suppress 5D matrix overhead
for(int j=0;j<N_j;j+=schurBlock){
///////////////////////////////////////////////////////////////
// Get the W and V vectors for this schurBlock^2 set of terms
///////////////////////////////////////////////////////////////
int N_ii = MIN(N_i-i,schurBlock);
int N_jj = MIN(N_j-j,schurBlock);
t_schur-=usecond();
for(int ii =0;ii < N_ii;ii++) a2a.return_w(i+ii, tmp_5d, w[ii]);
for(int jj =0;jj < N_jj;jj++) a2a.return_v(j+jj, tmp_5d, v[jj]);
t_schur+=usecond();
LOG(Message) << "Found w vectors " << i <<" .. " << i+N_ii-1 << std::endl;
LOG(Message) << "Found v vectors " << j <<" .. " << j+N_jj-1 << std::endl;
///////////////////////////////////////////////////////////////
// Series of cache blocked chunks of the contractions within this SchurBlock
///////////////////////////////////////////////////////////////
for(int ii=0;ii<N_ii;ii+=cacheBlock){
for(int jj=0;jj<N_jj;jj+=cacheBlock){
int N_iii = MIN(N_ii-ii,cacheBlock);
int N_jjj = MIN(N_jj-jj,cacheBlock);
Eigen::Tensor<ComplexD,5> mesonFieldBlocked(nmom,ngamma,nt,N_iii,N_jjj);
t_contr-=usecond();
A2Autils<FImpl>::MesonField(mesonFieldBlocked,
&w[ii],
&v[jj], gammas, phases,Tp);
t_contr+=usecond();
flops += vol * ( 2 * 8.0 + 6.0 + 8.0*nmom) * N_iii*N_jjj*ngamma;
bytes += vol * (12.0 * sizeof(Complex) ) * N_iii*N_jjj
+ vol * ( 2.0 * sizeof(Complex) *nmom ) * N_iii*N_jjj* ngamma;
///////////////////////////////////////////////////////////////
// Copy back to full meson field tensor
///////////////////////////////////////////////////////////////
parallel_for_nest2(int iii=0;iii< N_iii;iii++) {
for(int jjj=0;jjj< N_jjj;jjj++) {
for(int m =0;m< nmom;m++) {
for(int g =0;g< ngamma;g++) {
for(int t =0;t< nt;t++) {
mesonField(m,g,t,i+ii+iii,j+jj+jjj) = mesonFieldBlocked(m,g,t,iii,jjj);
}}}
}}
}}
}}
double nodes=grid->NodeCount();
double t1 = usecond();
LOG(Message) << " Contraction of MesonFields took "<<(t1-t0)/1.0e6<< " seconds " << std::endl;
LOG(Message) << " Schur "<<(t_schur)/1.0e6<< " seconds " << std::endl;
LOG(Message) << " Contr "<<(t_contr)/1.0e6<< " seconds " << std::endl;
/////////////////////////////////////////////////////////////////////////
// Test: Build the pion correlator (two end)
// < PI_ij(t0) PI_ji (t0+t) >
/////////////////////////////////////////////////////////////////////////
std::vector<ComplexD> corr(nt,ComplexD(0.0));
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
int m=0; // first momentum
int g=0; // first gamma in above ordering is gamma5 for pion
for(int t0=0;t0<nt;t0++){
for(int t=0;t<nt;t++){
int tt = (t0+t)%nt;
corr[t] += mesonField(m,g,t0,i,j)* mesonField(m,g,tt,j,i);
}}
}}
for(int t=0;t<nt;t++) corr[t] = corr[t]/ (double)nt;
for(int t=0;t<nt;t++) LOG(Message) << " " << t << " " << corr[t]<<std::endl;
// saveResult(par().output, "meson", result);
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MContraction_A2AMesonField_hpp_

View File

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

View File

@ -1,502 +0,0 @@
#ifndef Hadrons_MContraction_A2APionField_hpp_
#define Hadrons_MContraction_A2APionField_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
#include <Grid/Hadrons/AllToAllVectors.hpp>
#include <Grid/Hadrons/Modules/MContraction/A2Autils.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* A2APionField *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MContraction)
typedef std::pair<Gamma::Algebra, Gamma::Algebra> GammaPair;
class A2APionFieldPar : Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(A2APionFieldPar,
int, cacheBlock,
int, schurBlock,
int, Nmom,
std::string, A2A_i,
std::string, A2A_j,
std::string, output);
};
template <typename FImpl>
class TA2APionField : public Module<A2APionFieldPar>
{
public:
FERM_TYPE_ALIASES(FImpl, );
SOLVER_TYPE_ALIASES(FImpl, );
typedef typename FImpl::SiteSpinor vobj;
typedef typename vobj::scalar_object sobj;
typedef typename vobj::scalar_type scalar_type;
typedef typename vobj::vector_type vector_type;
typedef iSpinMatrix<vector_type> SpinMatrix_v;
typedef iSpinMatrix<scalar_type> SpinMatrix_s;
typedef iSinglet<vector_type> Scalar_v;
typedef iSinglet<scalar_type> Scalar_s;
typedef A2AModesSchurDiagTwo<typename FImpl::FermionField, FMat, Solver> A2ABase;
public:
// constructor
TA2APionField(const std::string name);
// destructor
virtual ~TA2APionField(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(A2APionField, ARG(TA2APionField<FIMPL>), MContraction);
MODULE_REGISTER(ZA2APionField, ARG(TA2APionField<ZFIMPL>), MContraction);
/******************************************************************************
* TA2APionField implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TA2APionField<FImpl>::TA2APionField(const std::string name)
: Module<A2APionFieldPar>(name)
{
}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TA2APionField<FImpl>::getInput(void)
{
std::vector<std::string> in;
in.push_back(par().A2A_i + "_class");
in.push_back(par().A2A_i + "_w_high_4d");
in.push_back(par().A2A_i + "_v_high_4d");
in.push_back(par().A2A_j + "_class");
in.push_back(par().A2A_j + "_w_high_4d");
in.push_back(par().A2A_j + "_v_high_4d");
return in;
}
template <typename FImpl>
std::vector<std::string> TA2APionField<FImpl>::getOutput(void)
{
std::vector<std::string> out = {};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TA2APionField<FImpl>::setup(void)
{
// Four D fields
envTmp(std::vector<FermionField>, "wi", 1, par().schurBlock, FermionField(env().getGrid(1)));
envTmp(std::vector<FermionField>, "vi", 1, par().schurBlock, FermionField(env().getGrid(1)));
envTmp(std::vector<FermionField>, "wj", 1, par().schurBlock, FermionField(env().getGrid(1)));
envTmp(std::vector<FermionField>, "vj", 1, par().schurBlock, FermionField(env().getGrid(1)));
// 5D tmp
int Ls_i = env().getObjectLs(par().A2A_i + "_class");
envTmpLat(FermionField, "tmp_5d", Ls_i);
int Ls_j= env().getObjectLs(par().A2A_j + "_class");
assert ( Ls_i == Ls_j );
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TA2APionField<FImpl>::execute(void)
{
LOG(Message) << "Computing A2A Pion fields" << std::endl;
auto &a2a_i = envGet(A2ABase, par().A2A_i + "_class");
auto &a2a_j = envGet(A2ABase, par().A2A_j + "_class");
///////////////////////////////////////////////
// Square assumption for now Nl = Nr = N
///////////////////////////////////////////////
int nt = env().getDim(Tp);
int nx = env().getDim(Xp);
int ny = env().getDim(Yp);
int nz = env().getDim(Zp);
// int N_i = a2a_i.par().N;
// int N_j = a2a_j.par().N;
int N_i = a2a_i.getN();
int N_j = a2a_j.getN();
int nmom=par().Nmom;
int schurBlock = par().schurBlock;
int cacheBlock = par().cacheBlock;
///////////////////////////////////////////////
// Momentum setup
///////////////////////////////////////////////
GridBase *grid = env().getGrid(1);
std::vector<LatticeComplex> phases(nmom,grid);
for(int m=0;m<nmom;m++){
phases[m] = Complex(1.0); // All zero momentum for now
}
///////////////////////////////////////////////////////////////////////
// i and j represent different flavours, hits, with different ranks.
// in general non-square case.
///////////////////////////////////////////////////////////////////////
Eigen::Tensor<ComplexD,4> pionFieldWVmom_ij (nmom,nt,N_i,N_j);
Eigen::Tensor<ComplexD,3> pionFieldWV_ij (nt,N_i,N_j);
Eigen::Tensor<ComplexD,4> pionFieldWVmom_ji (nmom,nt,N_j,N_i);
Eigen::Tensor<ComplexD,3> pionFieldWV_ji (nt,N_j,N_i);
LOG(Message) << "Rank for A2A PionField is " << N_i << " x "<<N_j << std::endl;
envGetTmp(std::vector<FermionField>, wi);
envGetTmp(std::vector<FermionField>, vi);
envGetTmp(std::vector<FermionField>, wj);
envGetTmp(std::vector<FermionField>, vj);
envGetTmp(FermionField, tmp_5d);
LOG(Message) << "Finding v and w vectors " << std::endl;
//////////////////////////////////////////////////////////////////////////
// i,j is first loop over SchurBlock factors reusing 5D matrices
// ii,jj is second loop over cacheBlock factors for high perf contractoin
// iii,jjj are loops within cacheBlock
// Total index is sum of these i+ii+iii etc...
//////////////////////////////////////////////////////////////////////////
double flops = 0.0;
double bytes = 0.0;
double vol = nx*ny*nz*nt;
double vol3 = nx*ny*nz;
double t_schur=0;
double t_contr_vwm=0;
double t_contr_vw=0;
double t_contr_ww=0;
double t_contr_vv=0;
double tt0 = usecond();
for(int i=0;i<N_i;i+=schurBlock){ //loop over SchurBlocking to suppress 5D matrix overhead
for(int j=0;j<N_j;j+=schurBlock){
///////////////////////////////////////////////////////////////
// Get the W and V vectors for this schurBlock^2 set of terms
///////////////////////////////////////////////////////////////
int N_ii = MIN(N_i-i,schurBlock);
int N_jj = MIN(N_j-j,schurBlock);
t_schur-=usecond();
for(int ii =0;ii < N_ii;ii++) a2a_i.return_w(i+ii, tmp_5d, wi[ii]);
for(int jj =0;jj < N_jj;jj++) a2a_j.return_w(j+jj, tmp_5d, wj[jj]);
for(int ii =0;ii < N_ii;ii++) a2a_i.return_v(i+ii, tmp_5d, vi[ii]);
for(int jj =0;jj < N_jj;jj++) a2a_j.return_v(j+jj, tmp_5d, vj[jj]);
t_schur+=usecond();
LOG(Message) << "Found i w&v vectors " << i <<" .. " << i+N_ii-1 << std::endl;
LOG(Message) << "Found j w&v vectors " << j <<" .. " << j+N_jj-1 << std::endl;
///////////////////////////////////////////////////////////////
// Series of cache blocked chunks of the contractions within this SchurBlock
///////////////////////////////////////////////////////////////
for(int ii=0;ii<N_ii;ii+=cacheBlock){
for(int jj=0;jj<N_jj;jj+=cacheBlock){
int N_iii = MIN(N_ii-ii,cacheBlock);
int N_jjj = MIN(N_jj-jj,cacheBlock);
Eigen::Tensor<ComplexD,4> pionFieldWVmomB_ij(nmom,nt,N_iii,N_jjj);
Eigen::Tensor<ComplexD,4> pionFieldWVmomB_ji(nmom,nt,N_jjj,N_iii);
Eigen::Tensor<ComplexD,3> pionFieldWVB_ij(nt,N_iii,N_jjj);
Eigen::Tensor<ComplexD,3> pionFieldWVB_ji(nt,N_jjj,N_iii);
t_contr_vwm-=usecond();
A2Autils<FImpl>::PionFieldWVmom(pionFieldWVmomB_ij, &wi[ii], &vj[jj], phases,Tp);
A2Autils<FImpl>::PionFieldWVmom(pionFieldWVmomB_ji, &wj[jj], &vi[ii], phases,Tp);
t_contr_vwm+=usecond();
t_contr_vw-=usecond();
A2Autils<FImpl>::PionFieldWV(pionFieldWVB_ij, &wi[ii], &vj[jj],Tp);
A2Autils<FImpl>::PionFieldWV(pionFieldWVB_ji, &wj[jj], &vi[ii],Tp);
t_contr_vw+=usecond();
flops += vol * ( 2 * 8.0 + 6.0 + 8.0*nmom) * N_iii*N_jjj;
bytes += vol * (12.0 * sizeof(Complex) ) * N_iii*N_jjj
+ vol * ( 2.0 * sizeof(Complex) *nmom ) * N_iii*N_jjj;
///////////////////////////////////////////////////////////////
// Copy back to full meson field tensor
///////////////////////////////////////////////////////////////
parallel_for_nest2(int iii=0;iii< N_iii;iii++) {
for(int jjj=0;jjj< N_jjj;jjj++) {
for(int m =0;m< nmom;m++) {
for(int t =0;t< nt;t++) {
pionFieldWVmom_ij(m,t,i+ii+iii,j+jj+jjj) = pionFieldWVmomB_ij(m,t,iii,jjj);
pionFieldWVmom_ji(m,t,j+jj+jjj,i+ii+iii) = pionFieldWVmomB_ji(m,t,jjj,iii);
}}
for(int t =0;t< nt;t++) {
pionFieldWV_ij(t,i+ii+iii,j+jj+jjj) = pionFieldWVB_ij(t,iii,jjj);
pionFieldWV_ji(t,j+jj+jjj,i+ii+iii) = pionFieldWVB_ji(t,jjj,iii);
}
}}
}}
}}
double nodes=grid->NodeCount();
double tt1 = usecond();
LOG(Message) << " Contraction of PionFields took "<<(tt1-tt0)/1.0e6<< " seconds " << std::endl;
LOG(Message) << " Schur "<<(t_schur)/1.0e6<< " seconds " << std::endl;
LOG(Message) << " Contr WVmom "<<(t_contr_vwm)/1.0e6<< " seconds " << std::endl;
LOG(Message) << " Contr WV "<<(t_contr_vw)/1.0e6<< " seconds " << std::endl;
double t_kernel = t_contr_vwm;
LOG(Message) << " Arith "<<flops/(t_kernel)/1.0e3/nodes<< " Gflop/s / node " << std::endl;
LOG(Message) << " Arith "<<bytes/(t_kernel)/1.0e3/nodes<< " GB/s /node " << std::endl;
/////////////////////////////////////////////////////////////////////////
// Test: Build the pion correlator (two end)
// < PI_ij(t0) PI_ji (t0+t) >
/////////////////////////////////////////////////////////////////////////
std::vector<ComplexD> corrMom(nt,ComplexD(0.0));
for(int i=0;i<N_i;i++){
for(int j=0;j<N_j;j++){
int m=0; // first momentum
for(int t0=0;t0<nt;t0++){
for(int t=0;t<nt;t++){
int tt = (t0+t)%nt;
corrMom[t] += pionFieldWVmom_ij(m,t0,i,j)* pionFieldWVmom_ji(m,tt,j,i);
}}
}}
for(int t=0;t<nt;t++) corrMom[t] = corrMom[t]/ (double)nt;
for(int t=0;t<nt;t++) LOG(Message) << " C_vwm " << t << " " << corrMom[t]<<std::endl;
/////////////////////////////////////////////////////////////////////////
// Test: Build the pion correlator (two end) from zero mom contraction
// < PI_ij(t0) PI_ji (t0+t) >
/////////////////////////////////////////////////////////////////////////
std::vector<ComplexD> corr(nt,ComplexD(0.0));
for(int i=0;i<N_i;i++){
for(int j=0;j<N_j;j++){
for(int t0=0;t0<nt;t0++){
for(int t=0;t<nt;t++){
int tt = (t0+t)%nt;
corr[t] += pionFieldWV_ij(t0,i,j)* pionFieldWV_ji(tt,j,i);
}}
}}
for(int t=0;t<nt;t++) corr[t] = corr[t]/ (double)nt;
for(int t=0;t<nt;t++) LOG(Message) << " C_vw " << t << " " << corr[t]<<std::endl;
/////////////////////////////////////////////////////////////////////////
// Test: Build the pion correlator from zero mom contraction with revers
// charge flow
/////////////////////////////////////////////////////////////////////////
std::vector<ComplexD> corr_wwvv(nt,ComplexD(0.0));
wi.resize(N_i,grid);
vi.resize(N_i,grid);
wj.resize(N_j,grid);
vj.resize(N_j,grid);
for(int i =0;i < N_i;i++) a2a_i.return_v(i, tmp_5d, vi[i]);
for(int i =0;i < N_i;i++) a2a_i.return_w(i, tmp_5d, wi[i]);
for(int j =0;j < N_j;j++) a2a_j.return_v(j, tmp_5d, vj[j]);
for(int j =0;j < N_j;j++) a2a_j.return_w(j, tmp_5d, wj[j]);
Eigen::Tensor<ComplexD,3> pionFieldWW_ij (nt,N_i,N_j);
Eigen::Tensor<ComplexD,3> pionFieldVV_ji (nt,N_j,N_i);
Eigen::Tensor<ComplexD,3> pionFieldWW_ji (nt,N_j,N_i);
Eigen::Tensor<ComplexD,3> pionFieldVV_ij (nt,N_i,N_j);
A2Autils<FImpl>::PionFieldWW(pionFieldWW_ij, &wi[0], &wj[0],Tp);
A2Autils<FImpl>::PionFieldVV(pionFieldVV_ji, &vj[0], &vi[0],Tp);
A2Autils<FImpl>::PionFieldWW(pionFieldWW_ji, &wj[0], &wi[0],Tp);
A2Autils<FImpl>::PionFieldVV(pionFieldVV_ij, &vi[0], &vj[0],Tp);
for(int i=0;i<N_i;i++){
for(int j=0;j<N_j;j++){
for(int t0=0;t0<nt;t0++){
for(int t=0;t<nt;t++){
int tt = (t0+t)%nt;
corr_wwvv[t] += pionFieldWW_ij(t0,i,j)* pionFieldVV_ji(tt,j,i);
corr_wwvv[t] += pionFieldWW_ji(t0,j,i)* pionFieldVV_ij(tt,i,j);
}}
}}
for(int t=0;t<nt;t++) corr_wwvv[t] = corr_wwvv[t] / vol /2.0 ; // (ij+ji noise contribs if i!=j ).
for(int t=0;t<nt;t++) LOG(Message) << " C_wwvv " << t << " " << corr_wwvv[t]<<std::endl;
/////////////////////////////////////////////////////////////////////////
// This is only correct if there are NO low modes
// Use the "ii" case to construct possible Z wall source one end trick
/////////////////////////////////////////////////////////////////////////
std::vector<ComplexD> corr_z2(nt,ComplexD(0.0));
Eigen::Tensor<ComplexD,3> pionFieldWW (nt,N_i,N_i);
Eigen::Tensor<ComplexD,3> pionFieldVV (nt,N_i,N_i);
A2Autils<FImpl>::PionFieldWW(pionFieldWW, &wi[0], &wi[0],Tp);
A2Autils<FImpl>::PionFieldVV(pionFieldVV, &vi[0], &vi[0],Tp);
for(int i=0;i<N_i;i++){
for(int t0=0;t0<nt;t0++){
for(int t=0;t<nt;t++){
int tt = (t0+t)%nt;
corr_z2[t] += pionFieldWW(t0,i,i) * pionFieldVV(tt,i,i) /vol ;
}}
}
LOG(Message) << " C_z2 WARNING only correct if Nl == 0 "<<std::endl;
for(int t=0;t<nt;t++) LOG(Message) << " C_z2 " << t << " " << corr_z2[t]<<std::endl;
/////////////////////////////////////////////////////////////////////////
// Test: Build a bag contraction
/////////////////////////////////////////////////////////////////////////
Eigen::Tensor<ComplexD,2> DeltaF2_fig8 (nt,16);
Eigen::Tensor<ComplexD,2> DeltaF2_trtr (nt,16);
Eigen::Tensor<ComplexD,1> denom0 (nt);
Eigen::Tensor<ComplexD,1> denom1 (nt);
const int dT=16;
A2Autils<FImpl>::DeltaFeq2 (dT,dT,DeltaF2_fig8,DeltaF2_trtr,
denom0,denom1,
pionFieldWW_ij,&vi[0],&vj[0],Tp);
{
int g=0; // O_{VV+AA}
for(int t=0;t<nt;t++)
LOG(Message) << " Bag [" << t << ","<<g<<"] "
<< (DeltaF2_fig8(t,g)+DeltaF2_trtr(t,g))
/ ( 8.0/3.0 * denom0[t]*denom1[t])
<<std::endl;
}
/////////////////////////////////////////////////////////////////////////
// Test: Build a bag contraction the Z2 way
// Build a wall bag comparison assuming no low modes
/////////////////////////////////////////////////////////////////////////
LOG(Message) << " Bag_z2 WARNING only correct if Nl == 0 "<<std::endl;
int t0=0;
int t1=dT;
int Nl=0;
LatticePropagator Qd0(grid);
LatticePropagator Qd1(grid);
LatticePropagator Qs0(grid);
LatticePropagator Qs1(grid);
for(int s=0;s<4;s++){
for(int c=0;c<3;c++){
int idx0 = Nl+t0*12+s*3+c;
int idx1 = Nl+t1*12+s*3+c;
FermToProp<FImpl>(Qd0, vi[idx0], s, c);
FermToProp<FImpl>(Qd1, vi[idx1], s, c);
FermToProp<FImpl>(Qs0, vj[idx0], s, c);
FermToProp<FImpl>(Qs1, vj[idx1], s, c);
}
}
std::vector<Gamma::Algebra> gammas ( {
Gamma::Algebra::GammaX,
Gamma::Algebra::GammaY,
Gamma::Algebra::GammaZ,
Gamma::Algebra::GammaT,
Gamma::Algebra::GammaXGamma5,
Gamma::Algebra::GammaYGamma5,
Gamma::Algebra::GammaZGamma5,
Gamma::Algebra::GammaTGamma5,
Gamma::Algebra::Identity,
Gamma::Algebra::Gamma5,
Gamma::Algebra::SigmaXY,
Gamma::Algebra::SigmaXZ,
Gamma::Algebra::SigmaXT,
Gamma::Algebra::SigmaYZ,
Gamma::Algebra::SigmaYT,
Gamma::Algebra::SigmaZT
});
auto G5 = Gamma::Algebra::Gamma5;
LatticePropagator anti_d0 = adj( Gamma(G5) * Qd0 * Gamma(G5));
LatticePropagator anti_d1 = adj( Gamma(G5) * Qd1 * Gamma(G5));
LatticeComplex TR1(grid);
LatticeComplex TR2(grid);
LatticeComplex Wick1(grid);
LatticeComplex Wick2(grid);
LatticePropagator PR1(grid);
LatticePropagator PR2(grid);
PR1 = Qs0 * Gamma(G5) * anti_d0;
PR2 = Qs1 * Gamma(G5) * anti_d1;
for(int g=0;g<Nd*Nd;g++){
auto g1 = gammas[g];
Gamma G1 (g1);
TR1 = trace( PR1 * G1 );
TR2 = trace( PR2 * G1 );
Wick1 = TR1*TR2;
Wick2 = trace( PR1* G1 * PR2 * G1 );
std::vector<TComplex> C1;
std::vector<TComplex> C2;
std::vector<TComplex> C3;
sliceSum(Wick1,C1, Tp);
sliceSum(Wick2,C2, Tp);
sliceSum(TR1 ,C3, Tp);
/*
if(g<5){
for(int t=0;t<C1.size();t++){
LOG(Message) << " Wick1["<<g<<","<<t<< "] "<< C1[t]<<std::endl;
}
for(int t=0;t<C2.size();t++){
LOG(Message) << " Wick2["<<g<<","<<t<< "] "<< C2[t]<<std::endl;
}
}
if( (g==9) || (g==7) ){ // P and At in above ordering
for(int t=0;t<C3.size();t++){
LOG(Message) << " <G|P>["<<g<<","<<t<< "] "<< C3[t]<<std::endl;
}
}
*/
}
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MContraction_A2APionField_hpp_

File diff suppressed because it is too large Load Diff

View File

@ -1,35 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MContraction/Baryon.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 <Grid/Hadrons/Modules/MContraction/Baryon.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MContraction;
template class Grid::Hadrons::MContraction::TBaryon<FIMPL,FIMPL,FIMPL>;

View File

@ -68,7 +68,7 @@ public:
// constructor
TBaryon(const std::string name);
// destructor
virtual ~TBaryon(void) {};
virtual ~TBaryon(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
@ -79,7 +79,7 @@ protected:
virtual void execute(void);
};
MODULE_REGISTER_TMP(Baryon, ARG(TBaryon<FIMPL, FIMPL, FIMPL>), MContraction);
MODULE_REGISTER_NS(Baryon, ARG(TBaryon<FIMPL, FIMPL, FIMPL>), MContraction);
/******************************************************************************
* TBaryon implementation *
@ -122,6 +122,7 @@ void TBaryon<FImpl1, FImpl2, FImpl3>::execute(void)
<< " quarks '" << par().q1 << "', '" << par().q2 << "', and '"
<< par().q3 << "'" << std::endl;
ResultWriter writer(RESULT_FILE_NAME(par().output));
auto &q1 = envGet(PropagatorField1, par().q1);
auto &q2 = envGet(PropagatorField2, par().q2);
auto &q3 = envGet(PropagatorField3, par().q2);
@ -130,7 +131,7 @@ void TBaryon<FImpl1, FImpl2, FImpl3>::execute(void)
// FIXME: do contractions
// saveResult(par().output, "meson", result);
// write(writer, "meson", result);
}
END_MODULE_NAMESPACE

View File

@ -1,35 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MContraction/DiscLoop.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 <Grid/Hadrons/Modules/MContraction/DiscLoop.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MContraction;
template class Grid::Hadrons::MContraction::TDiscLoop<FIMPL>;

View File

@ -65,7 +65,7 @@ public:
// constructor
TDiscLoop(const std::string name);
// destructor
virtual ~TDiscLoop(void) {};
virtual ~TDiscLoop(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
@ -76,7 +76,7 @@ protected:
virtual void execute(void);
};
MODULE_REGISTER_TMP(DiscLoop, TDiscLoop<FIMPL>, MContraction);
MODULE_REGISTER_NS(DiscLoop, TDiscLoop<FIMPL>, MContraction);
/******************************************************************************
* TDiscLoop implementation *
@ -119,6 +119,7 @@ void TDiscLoop<FImpl>::execute(void)
<< "' using '" << par().q_loop << "' with " << par().gamma
<< " insertion." << std::endl;
ResultWriter writer(RESULT_FILE_NAME(par().output));
auto &q_loop = envGet(PropagatorField, par().q_loop);
Gamma gamma(par().gamma);
std::vector<TComplex> buf;
@ -127,13 +128,15 @@ void TDiscLoop<FImpl>::execute(void)
envGetTmp(LatticeComplex, c);
c = trace(gamma*q_loop);
sliceSum(c, buf, Tp);
result.gamma = par().gamma;
result.corr.resize(buf.size());
for (unsigned int t = 0; t < buf.size(); ++t)
{
result.corr[t] = TensorRemove(buf[t]);
}
saveResult(par().output, "disc", result);
write(writer, "disc", result);
}
END_MODULE_NAMESPACE

View File

@ -1,35 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MContraction/Gamma3pt.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 <Grid/Hadrons/Modules/MContraction/Gamma3pt.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MContraction;
template class Grid::Hadrons::MContraction::TGamma3pt<FIMPL,FIMPL,FIMPL>;

View File

@ -96,7 +96,7 @@ public:
// constructor
TGamma3pt(const std::string name);
// destructor
virtual ~TGamma3pt(void) {};
virtual ~TGamma3pt(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
@ -107,7 +107,7 @@ protected:
virtual void execute(void);
};
MODULE_REGISTER_TMP(Gamma3pt, ARG(TGamma3pt<FIMPL, FIMPL, FIMPL>), MContraction);
MODULE_REGISTER_NS(Gamma3pt, ARG(TGamma3pt<FIMPL, FIMPL, FIMPL>), MContraction);
/******************************************************************************
* TGamma3pt implementation *
@ -153,6 +153,7 @@ void TGamma3pt<FImpl1, FImpl2, FImpl3>::execute(void)
// Initialise variables. q2 and q3 are normal propagators, q1 may be
// sink smeared.
ResultWriter writer(RESULT_FILE_NAME(par().output));
auto &q1 = envGet(SlicedPropagator1, par().q1);
auto &q2 = envGet(PropagatorField2, par().q2);
auto &q3 = envGet(PropagatorField2, par().q3);
@ -174,7 +175,8 @@ void TGamma3pt<FImpl1, FImpl2, FImpl3>::execute(void)
{
result.corr[t] = TensorRemove(buf[t]);
}
saveResult(par().output, "gamma3pt", result);
write(writer, "gamma3pt", result);
}
END_MODULE_NAMESPACE

View File

@ -1,35 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MContraction/Meson.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 <Grid/Hadrons/Modules/MContraction/Meson.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MContraction;
template class Grid::Hadrons::MContraction::TMeson<FIMPL,FIMPL>;

View File

@ -45,8 +45,8 @@ BEGIN_HADRONS_NAMESPACE
- q1: input propagator 1 (string)
- q2: input propagator 2 (string)
- gammas: gamma products to insert at sink & source, pairs of gamma matrices
(space-separated strings) in round brackets (i.e. (g_sink g_src)),
in a sequence (e.g. "(Gamma5 Gamma5)(Gamma5 GammaT)").
(space-separated strings) in angled brackets (i.e. <g_sink g_src>),
in a sequence (e.g. "<Gamma5 Gamma5><Gamma5 GammaT>").
Special values: "all" - perform all possible contractions.
- sink: module to compute the sink to use in contraction (string).
@ -90,7 +90,7 @@ public:
// constructor
TMeson(const std::string name);
// destructor
virtual ~TMeson(void) {};
virtual ~TMeson(void) = default;
// dependencies/products
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
@ -102,7 +102,7 @@ protected:
virtual void execute(void);
};
MODULE_REGISTER_TMP(Meson, ARG(TMeson<FIMPL, FIMPL>), MContraction);
MODULE_REGISTER_NS(Meson, ARG(TMeson<FIMPL, FIMPL>), MContraction);
/******************************************************************************
* TMeson implementation *
@ -172,6 +172,7 @@ void TMeson<FImpl1, FImpl2>::execute(void)
<< " quarks '" << par().q1 << "' and '" << par().q2 << "'"
<< std::endl;
ResultWriter writer(RESULT_FILE_NAME(par().output));
std::vector<TComplex> buf;
std::vector<Result> result;
Gamma g5(Gamma::Algebra::Gamma5);
@ -238,7 +239,7 @@ void TMeson<FImpl1, FImpl2>::execute(void)
}
}
}
saveResult(par().output, "meson", result);
write(writer, "meson", result);
}
END_MODULE_NAMESPACE

View File

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

View File

@ -1,269 +0,0 @@
#ifndef Hadrons_MContraction_MesonFieldGamma_hpp_
#define Hadrons_MContraction_MesonFieldGamma_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
#include <Grid/Hadrons/AllToAllVectors.hpp>
#include <Grid/Hadrons/AllToAllReduction.hpp>
#include <Grid/Grid_Eigen_Dense.h>
#include <fstream>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* MesonFieldGamma *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MContraction)
class MesonFieldPar : Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(MesonFieldPar,
int, Nl,
int, N,
int, Nblock,
std::string, A2A1,
std::string, A2A2,
std::string, gammas,
std::string, output);
};
template <typename FImpl>
class TMesonFieldGamma : public Module<MesonFieldPar>
{
public:
FERM_TYPE_ALIASES(FImpl, );
SOLVER_TYPE_ALIASES(FImpl, );
typedef A2AModesSchurDiagTwo<typename FImpl::FermionField, FMat, Solver> A2ABase;
class Result : Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
Gamma::Algebra, gamma,
std::vector<std::vector<std::vector<ComplexD>>>, MesonField);
};
public:
// constructor
TMesonFieldGamma(const std::string name);
// destructor
virtual ~TMesonFieldGamma(void){};
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
virtual void parseGammaString(std::vector<Gamma::Algebra> &gammaList);
virtual void vectorOfWs(std::vector<FermionField> &w, int i, int Nblock, FermionField &tmpw_5d, std::vector<FermionField> &vec_w);
virtual void vectorOfVs(std::vector<FermionField> &v, int j, int Nblock, FermionField &tmpv_5d, std::vector<FermionField> &vec_v);
virtual void gammaMult(std::vector<FermionField> &v, Gamma gamma);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
};
MODULE_REGISTER(MesonFieldGamma, ARG(TMesonFieldGamma<FIMPL>), MContraction);
MODULE_REGISTER(ZMesonFieldGamma, ARG(TMesonFieldGamma<ZFIMPL>), MContraction);
/******************************************************************************
* TMesonFieldGamma implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TMesonFieldGamma<FImpl>::TMesonFieldGamma(const std::string name)
: Module<MesonFieldPar>(name)
{
}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TMesonFieldGamma<FImpl>::getInput(void)
{
std::vector<std::string> in = {par().A2A1 + "_class", par().A2A2 + "_class"};
in.push_back(par().A2A1 + "_w_high_4d");
in.push_back(par().A2A2 + "_v_high_4d");
return in;
}
template <typename FImpl>
std::vector<std::string> TMesonFieldGamma<FImpl>::getOutput(void)
{
std::vector<std::string> out = {};
return out;
}
template <typename FImpl>
void TMesonFieldGamma<FImpl>::parseGammaString(std::vector<Gamma::Algebra> &gammaList)
{
gammaList.clear();
// Determine gamma matrices to insert at source/sink.
if (par().gammas.compare("all") == 0)
{
// Do all contractions.
for (unsigned int i = 1; i < Gamma::nGamma; i += 2)
{
gammaList.push_back(((Gamma::Algebra)i));
}
}
else
{
// Parse individual contractions from input string.
gammaList = strToVec<Gamma::Algebra>(par().gammas);
}
}
template <typename FImpl>
void TMesonFieldGamma<FImpl>::vectorOfWs(std::vector<FermionField> &w, int i, int Nblock, FermionField &tmpw_5d, std::vector<FermionField> &vec_w)
{
for (unsigned int ni = 0; ni < Nblock; ni++)
{
vec_w[ni] = w[i + ni];
}
}
template <typename FImpl>
void TMesonFieldGamma<FImpl>::vectorOfVs(std::vector<FermionField> &v, int j, int Nblock, FermionField &tmpv_5d, std::vector<FermionField> &vec_v)
{
for (unsigned int nj = 0; nj < Nblock; nj++)
{
vec_v[nj] = v[j+nj];
}
}
template <typename FImpl>
void TMesonFieldGamma<FImpl>::gammaMult(std::vector<FermionField> &v, Gamma gamma)
{
int Nblock = v.size();
for (unsigned int nj = 0; nj < Nblock; nj++)
{
v[nj] = gamma * v[nj];
}
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TMesonFieldGamma<FImpl>::setup(void)
{
int nt = env().getDim(Tp);
int N = par().N;
int Nblock = par().Nblock;
int Ls_ = env().getObjectLs(par().A2A1 + "_class");
envTmpLat(FermionField, "tmpv_5d", Ls_);
envTmpLat(FermionField, "tmpw_5d", Ls_);
envTmp(std::vector<FermionField>, "w", 1, N, FermionField(env().getGrid(1)));
envTmp(std::vector<FermionField>, "v", 1, N, FermionField(env().getGrid(1)));
envTmp(Eigen::MatrixXcd, "MF", 1, Eigen::MatrixXcd::Zero(nt, N * N));
envTmp(std::vector<FermionField>, "w_block", 1, Nblock, FermionField(env().getGrid(1)));
envTmp(std::vector<FermionField>, "v_block", 1, Nblock, FermionField(env().getGrid(1)));
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TMesonFieldGamma<FImpl>::execute(void)
{
LOG(Message) << "Computing A2A meson field for gamma = " << par().gammas << ", taking w from " << par().A2A1 << " and v from " << par().A2A2 << std::endl;
int N = par().N;
int nt = env().getDim(Tp);
int Nblock = par().Nblock;
std::vector<Result> result;
std::vector<Gamma::Algebra> gammaResultList;
std::vector<Gamma> gammaList;
parseGammaString(gammaResultList);
result.resize(gammaResultList.size());
Gamma g5(Gamma::Algebra::Gamma5);
gammaList.resize(gammaResultList.size(), g5);
for (unsigned int i = 0; i < result.size(); ++i)
{
result[i].gamma = gammaResultList[i];
result[i].MesonField.resize(N, std::vector<std::vector<ComplexD>>(N, std::vector<ComplexD>(nt)));
Gamma gamma(gammaResultList[i]);
gammaList[i] = gamma;
}
auto &a2a1 = envGet(A2ABase, par().A2A1 + "_class");
auto &a2a2 = envGet(A2ABase, par().A2A2 + "_class");
envGetTmp(FermionField, tmpv_5d);
envGetTmp(FermionField, tmpw_5d);
envGetTmp(std::vector<FermionField>, v);
envGetTmp(std::vector<FermionField>, w);
LOG(Message) << "Finding v and w vectors for N = " << N << std::endl;
for (int i = 0; i < N; i++)
{
a2a2.return_v(i, tmpv_5d, v[i]);
a2a1.return_w(i, tmpw_5d, w[i]);
}
LOG(Message) << "Found v and w vectors for N = " << N << std::endl;
std::vector<std::vector<ComplexD>> MesonField_ij;
LOG(Message) << "Before blocked MFs, Nblock = " << Nblock << std::endl;
envGetTmp(std::vector<FermionField>, v_block);
envGetTmp(std::vector<FermionField>, w_block);
MesonField_ij.resize(Nblock * Nblock, std::vector<ComplexD>(nt));
envGetTmp(Eigen::MatrixXcd, MF);
LOG(Message) << "Before blocked MFs, Nblock = " << Nblock << std::endl;
for (unsigned int i = 0; i < N; i += Nblock)
{
vectorOfWs(w, i, Nblock, tmpw_5d, w_block);
for (unsigned int j = 0; j < N; j += Nblock)
{
vectorOfVs(v, j, Nblock, tmpv_5d, v_block);
for (unsigned int k = 0; k < result.size(); k++)
{
gammaMult(v_block, gammaList[k]);
sliceInnerProductMesonField(MesonField_ij, w_block, v_block, Tp);
for (unsigned int nj = 0; nj < Nblock; nj++)
{
for (unsigned int ni = 0; ni < Nblock; ni++)
{
MF.col((i + ni) + (j + nj) * N) = Eigen::VectorXcd::Map(&MesonField_ij[nj * Nblock + ni][0], MesonField_ij[nj * Nblock + ni].size());
}
}
}
}
if (i % 10 == 0)
{
LOG(Message) << "MF for i = " << i << " of " << N << std::endl;
}
}
LOG(Message) << "Before Global sum, Nblock = " << Nblock << std::endl;
v_block[0]._grid->GlobalSumVector(MF.data(), MF.size());
LOG(Message) << "After Global sum, Nblock = " << Nblock << std::endl;
for (unsigned int i = 0; i < N; i++)
{
for (unsigned int j = 0; j < N; j++)
{
for (unsigned int k = 0; k < result.size(); k++)
{
for (unsigned int t = 0; t < nt; t++)
{
result[k].MesonField[i][j][t] = MF.col(i + N * j)[t];
}
}
}
}
saveResult(par().output, "meson", result);
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MContraction_MesonFieldGm_hpp_

View File

@ -1,35 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MContraction/WardIdentity.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 <Grid/Hadrons/Modules/MContraction/WardIdentity.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MContraction;
template class Grid::Hadrons::MContraction::TWardIdentity<FIMPL>;

View File

@ -71,7 +71,7 @@ public:
// constructor
TWardIdentity(const std::string name);
// destructor
virtual ~TWardIdentity(void) {};
virtual ~TWardIdentity(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
@ -84,7 +84,7 @@ private:
unsigned int Ls_;
};
MODULE_REGISTER_TMP(WardIdentity, TWardIdentity<FIMPL>, MContraction);
MODULE_REGISTER_NS(WardIdentity, TWardIdentity<FIMPL>, MContraction);
/******************************************************************************
* TWardIdentity implementation *
@ -119,7 +119,7 @@ void TWardIdentity<FImpl>::setup(void)
Ls_ = env().getObjectLs(par().q);
if (Ls_ != env().getObjectLs(par().action))
{
HADRONS_ERROR(Size, "Ls mismatch between quark action and propagator");
HADRON_ERROR(Size, "Ls mismatch between quark action and propagator");
}
envTmpLat(PropagatorField, "tmp");
envTmpLat(PropagatorField, "vector_WI");

View File

@ -97,7 +97,7 @@ public:\
/* constructor */ \
T##modname(const std::string name);\
/* destructor */ \
virtual ~T##modname(void) {};\
virtual ~T##modname(void) = default;\
/* dependency relation */ \
virtual std::vector<std::string> getInput(void);\
virtual std::vector<std::string> getOutput(void);\
@ -109,7 +109,7 @@ protected:\
/* execution */ \
virtual void execute(void);\
};\
MODULE_REGISTER(modname, T##modname, MContraction);
MODULE_REGISTER_NS(modname, T##modname, MContraction);
END_MODULE_NAMESPACE

View File

@ -104,6 +104,7 @@ void TWeakHamiltonianEye::execute(void)
<< par().q2 << ", '" << par().q3 << "' and '" << par().q4
<< "'." << std::endl;
ResultWriter writer(RESULT_FILE_NAME(par().output));
auto &q1 = envGet(SlicedPropagator, par().q1);
auto &q2 = envGet(PropagatorField, par().q2);
auto &q3 = envGet(PropagatorField, par().q3);
@ -146,6 +147,5 @@ void TWeakHamiltonianEye::execute(void)
SUM_MU(expbuf, E_body[mu]*E_loop[mu])
MAKE_DIAG(expbuf, corrbuf, result[E_diag], "HW_E")
// IO
saveResult(par().output, "HW_Eye", result);
write(writer, "HW_Eye", result);
}

View File

@ -104,6 +104,7 @@ void TWeakHamiltonianNonEye::execute(void)
<< par().q2 << ", '" << par().q3 << "' and '" << par().q4
<< "'." << std::endl;
ResultWriter writer(RESULT_FILE_NAME(par().output));
auto &q1 = envGet(PropagatorField, par().q1);
auto &q2 = envGet(PropagatorField, par().q2);
auto &q3 = envGet(PropagatorField, par().q3);
@ -143,6 +144,5 @@ void TWeakHamiltonianNonEye::execute(void)
SUM_MU(expbuf, W_i_side_loop[mu]*W_f_side_loop[mu])
MAKE_DIAG(expbuf, corrbuf, result[W_diag], "HW_W")
// IO
saveResult(par().output, "HW_NonEye", result);
write(writer, "HW_NonEye", result);
}

View File

@ -104,6 +104,7 @@ void TWeakNeutral4ptDisc::execute(void)
<< par().q2 << ", '" << par().q3 << "' and '" << par().q4
<< "'." << std::endl;
ResultWriter writer(RESULT_FILE_NAME(par().output));
auto &q1 = envGet(PropagatorField, par().q1);
auto &q2 = envGet(PropagatorField, par().q2);
auto &q3 = envGet(PropagatorField, par().q3);
@ -137,6 +138,5 @@ void TWeakNeutral4ptDisc::execute(void)
expbuf *= curr;
MAKE_DIAG(expbuf, corrbuf, result[neut_disc_2_diag], "HW_disc0_2")
// IO
saveResult(par().output, "HW_disc0", result);
write(writer, "HW_disc0", result);
}

View File

@ -1,36 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MFermion/FreeProp.cc
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
Author: Vera Guelpers <V.M.Guelpers@soton.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 <Grid/Hadrons/Modules/MFermion/FreeProp.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MFermion;
template class Grid::Hadrons::MFermion::TFreeProp<FIMPL>;

View File

@ -1,187 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MFermion/FreeProp.hpp
Copyright (C) 2015-2018
Author: Vera Guelpers <V.M.Guelpers@soton.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_MFermion_FreeProp_hpp_
#define Hadrons_MFermion_FreeProp_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* FreeProp *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MFermion)
class FreePropPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(FreePropPar,
std::string, source,
std::string, action,
double, mass,
std::string, twist);
};
template <typename FImpl>
class TFreeProp: public Module<FreePropPar>
{
public:
FG_TYPE_ALIASES(FImpl,);
public:
// constructor
TFreeProp(const std::string name);
// destructor
virtual ~TFreeProp(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);
private:
unsigned int Ls_;
};
MODULE_REGISTER_TMP(FreeProp, TFreeProp<FIMPL>, MFermion);
/******************************************************************************
* TFreeProp implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TFreeProp<FImpl>::TFreeProp(const std::string name)
: Module<FreePropPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TFreeProp<FImpl>::getInput(void)
{
std::vector<std::string> in = {par().source, par().action};
return in;
}
template <typename FImpl>
std::vector<std::string> TFreeProp<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName(), getName() + "_5d"};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TFreeProp<FImpl>::setup(void)
{
Ls_ = env().getObjectLs(par().action);
envCreateLat(PropagatorField, getName());
envTmpLat(FermionField, "source", Ls_);
envTmpLat(FermionField, "sol", Ls_);
envTmpLat(FermionField, "tmp");
if (Ls_ > 1)
{
envCreateLat(PropagatorField, getName() + "_5d", Ls_);
}
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TFreeProp<FImpl>::execute(void)
{
LOG(Message) << "Computing free fermion propagator '" << getName() << "'"
<< std::endl;
std::string propName = (Ls_ == 1) ? getName() : (getName() + "_5d");
auto &prop = envGet(PropagatorField, propName);
auto &fullSrc = envGet(PropagatorField, par().source);
auto &mat = envGet(FMat, par().action);
RealD mass = par().mass;
envGetTmp(FermionField, source);
envGetTmp(FermionField, sol);
envGetTmp(FermionField, tmp);
LOG(Message) << "Calculating a free Propagator with mass " << mass
<< " using the action '" << par().action
<< "' on source '" << par().source << "'" << std::endl;
for (unsigned int s = 0; s < Ns; ++s)
for (unsigned int c = 0; c < FImpl::Dimension; ++c)
{
LOG(Message) << "Calculation for spin= " << s << ", color= " << c
<< std::endl;
// source conversion for 4D sources
if (!env().isObject5d(par().source))
{
if (Ls_ == 1)
{
PropToFerm<FImpl>(source, fullSrc, s, c);
}
else
{
PropToFerm<FImpl>(tmp, fullSrc, s, c);
mat.ImportPhysicalFermionSource(tmp, source);
}
}
// source conversion for 5D sources
else
{
if (Ls_ != env().getObjectLs(par().source))
{
HADRONS_ERROR(Size, "Ls mismatch between quark action and source");
}
else
{
PropToFerm<FImpl>(source, fullSrc, s, c);
}
}
sol = zero;
std::vector<Real> twist = strToVec<Real>(par().twist);
if(twist.size() != Nd) HADRONS_ERROR(Size, "number of twist angles does not match number of dimensions");
mat.FreePropagator(source,sol,mass,twist);
FermToProp<FImpl>(prop, sol, s, c);
// create 4D propagators from 5D one if necessary
if (Ls_ > 1)
{
PropagatorField &p4d = envGet(PropagatorField, getName());
mat.ExportPhysicalFermionSolution(sol, tmp);
FermToProp<FImpl>(p4d, tmp, s, c);
}
}
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MFermion_FreeProp_hpp_

View File

@ -1,35 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MFermion/GaugeProp.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 <Grid/Hadrons/Modules/MFermion/GaugeProp.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MFermion;
template class Grid::Hadrons::MFermion::TGaugeProp<FIMPL>;
template class Grid::Hadrons::MFermion::TGaugeProp<ZFIMPL>;

View File

@ -7,9 +7,7 @@ Source file: extras/Hadrons/Modules/MFermion/GaugeProp.hpp
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
Author: Guido Cossu <guido.cossu@ed.ac.uk>
Author: Lanny91 <andrew.lawson@gmail.com>
Author: pretidav <david.preti@csic.es>
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
@ -35,10 +33,30 @@ See the full license in the file "LICENSE" in the top level distribution directo
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
#include <Grid/Hadrons/Solver.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* 5D -> 4D and 4D -> 5D conversions. *
******************************************************************************/
template<class vobj> // Note that 5D object is modified.
inline void make_4D(Lattice<vobj> &in_5d, Lattice<vobj> &out_4d, int Ls)
{
axpby_ssp_pminus(in_5d, 0., in_5d, 1., in_5d, 0, 0);
axpby_ssp_pplus(in_5d, 1., in_5d, 1., in_5d, 0, Ls-1);
ExtractSlice(out_4d, in_5d, 0, 0);
}
template<class vobj>
inline void make_5D(Lattice<vobj> &in_4d, Lattice<vobj> &out_5d, int Ls)
{
out_5d = zero;
InsertSlice(in_4d, out_5d, 0, 0);
InsertSlice(in_4d, out_5d, Ls-1, 0);
axpby_ssp_pplus(out_5d, 0., out_5d, 1., out_5d, 0, 0);
axpby_ssp_pminus(out_5d, 0., out_5d, 1., out_5d, Ls-1, Ls-1);
}
/******************************************************************************
* GaugeProp *
******************************************************************************/
@ -56,13 +74,12 @@ template <typename FImpl>
class TGaugeProp: public Module<GaugePropPar>
{
public:
FG_TYPE_ALIASES(FImpl,);
SOLVER_TYPE_ALIASES(FImpl,);
FGS_TYPE_ALIASES(FImpl,);
public:
// constructor
TGaugeProp(const std::string name);
// destructor
virtual ~TGaugeProp(void) {};
virtual ~TGaugeProp(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
@ -73,12 +90,10 @@ protected:
virtual void execute(void);
private:
unsigned int Ls_;
Solver *solver_{nullptr};
SolverFn *solver_{nullptr};
};
MODULE_REGISTER_TMP(GaugeProp, TGaugeProp<FIMPL>, MFermion);
MODULE_REGISTER_TMP(ZGaugeProp, TGaugeProp<ZFIMPL>, MFermion);
MODULE_REGISTER_NS(GaugeProp, TGaugeProp<FIMPL>, MFermion);
/******************************************************************************
* TGaugeProp implementation *
******************************************************************************/
@ -130,8 +145,7 @@ void TGaugeProp<FImpl>::execute(void)
std::string propName = (Ls_ == 1) ? getName() : (getName() + "_5d");
auto &prop = envGet(PropagatorField, propName);
auto &fullSrc = envGet(PropagatorField, par().source);
auto &solver = envGet(Solver, par().solver);
auto &mat = solver.getFMat();
auto &solver = envGet(SolverFn, par().solver);
envGetTmp(FermionField, source);
envGetTmp(FermionField, sol);
@ -139,12 +153,11 @@ void TGaugeProp<FImpl>::execute(void)
LOG(Message) << "Inverting using solver '" << par().solver
<< "' on source '" << par().source << "'" << std::endl;
for (unsigned int s = 0; s < Ns; ++s)
for (unsigned int c = 0; c < FImpl::Dimension; ++c)
for (unsigned int c = 0; c < FImpl::Dimension; ++c)
{
LOG(Message) << "Inversion for spin= " << s << ", color= " << c
<< std::endl;
// source conversion for 4D sources
LOG(Message) << "Import source" << std::endl;
if (!env().isObject5d(par().source))
{
if (Ls_ == 1)
@ -154,7 +167,7 @@ void TGaugeProp<FImpl>::execute(void)
else
{
PropToFerm<FImpl>(tmp, fullSrc, s, c);
mat.ImportPhysicalFermionSource(tmp, source);
make_5D(tmp, source, Ls_);
}
}
// source conversion for 5D sources
@ -162,23 +175,21 @@ void TGaugeProp<FImpl>::execute(void)
{
if (Ls_ != env().getObjectLs(par().source))
{
HADRONS_ERROR(Size, "Ls mismatch between quark action and source");
HADRON_ERROR(Size, "Ls mismatch between quark action and source");
}
else
{
PropToFerm<FImpl>(source, fullSrc, s, c);
}
}
LOG(Message) << "Solve" << std::endl;
sol = zero;
solver(sol, source);
LOG(Message) << "Export solution" << std::endl;
FermToProp<FImpl>(prop, sol, s, c);
// create 4D propagators from 5D one if necessary
if (Ls_ > 1)
{
PropagatorField &p4d = envGet(PropagatorField, getName());
mat.ExportPhysicalFermionSolution(sol, tmp);
make_4D(sol, tmp, Ls_);
FermToProp<FImpl>(p4d, tmp, s, c);
}
}

View File

@ -4,11 +4,9 @@ Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MGauge/FundtoHirep.cc
Copyright (C) 2015-2018
Copyright (C) 2015
Copyright (C) 2016
Author: Antonin Portelli <antonin.portelli@me.com>
Author: Guido Cossu <guido.cossu@ed.ac.uk>
Author: pretidav <david.preti@csic.es>
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
@ -44,8 +42,7 @@ TFundtoHirep<Rep>::TFundtoHirep(const std::string name)
template <class Rep>
std::vector<std::string> TFundtoHirep<Rep>::getInput(void)
{
std::vector<std::string> in = {par().gaugeconf};
std::vector<std::string> in;
return in;
}
@ -53,7 +50,6 @@ template <class Rep>
std::vector<std::string> TFundtoHirep<Rep>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
@ -61,19 +57,19 @@ std::vector<std::string> TFundtoHirep<Rep>::getOutput(void)
template <typename Rep>
void TFundtoHirep<Rep>::setup(void)
{
envCreateLat(Rep::LatticeField, getName());
envCreateLat(typename Rep::LatticeField, getName());
}
// execution ///////////////////////////////////////////////////////////////////
template <class Rep>
void TFundtoHirep<Rep>::execute(void)
{
auto &U = *env().template getObject<LatticeGaugeField>(par().gaugeconf);
LOG(Message) << "Transforming Representation" << std::endl;
auto &U = envGet(LatticeGaugeField, par().gaugeconf);
auto &URep = envGet(Rep::LatticeField, getName());
Rep TargetRepresentation(U._grid);
TargetRepresentation.update_representation(U);
auto &URep = envGet(typename Rep::LatticeField, getName());
URep = TargetRepresentation.U;
}

View File

@ -4,10 +4,11 @@ Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MGauge/FundtoHirep.hpp
Copyright (C) 2015-2018
Copyright (C) 2015
Copyright (C) 2016
Author: Antonin Portelli <antonin.portelli@me.com>
Author: pretidav <david.preti@csic.es>
Author: David Preti <david.preti@to.infn.it>
Guido Cossu <guido.cossu@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
@ -55,7 +56,7 @@ public:
// constructor
TFundtoHirep(const std::string name);
// destructor
virtual ~TFundtoHirep(void) {};
virtual ~TFundtoHirep(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
@ -65,9 +66,9 @@ public:
void execute(void);
};
//MODULE_REGISTER_TMP(FundtoAdjoint, TFundtoHirep<AdjointRepresentation>, MGauge);
//MODULE_REGISTER_TMP(FundtoTwoIndexSym, TFundtoHirep<TwoIndexSymmetricRepresentation>, MGauge);
//MODULE_REGISTER_TMP(FundtoTwoIndexAsym, TFundtoHirep<TwoIndexAntiSymmetricRepresentation>, MGauge);
//MODULE_REGISTER_NS(FundtoAdjoint, TFundtoHirep<AdjointRepresentation>, MGauge);
//MODULE_REGISTER_NS(FundtoTwoIndexSym, TFundtoHirep<TwoIndexSymmetricRepresentation>, MGauge);
//MODULE_REGISTER_NS(FundtoTwoIndexAsym, TFundtoHirep<TwoIndexAntiSymmetricRepresentation>, MGauge);
END_MODULE_NAMESPACE

View File

@ -46,7 +46,7 @@ public:
// constructor
TRandom(const std::string name);
// destructor
virtual ~TRandom(void) {};
virtual ~TRandom(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
@ -57,7 +57,7 @@ protected:
virtual void execute(void);
};
MODULE_REGISTER(Random, TRandom, MGauge);
MODULE_REGISTER_NS(Random, TRandom, MGauge);
END_MODULE_NAMESPACE

View File

@ -7,8 +7,6 @@ Source file: extras/Hadrons/Modules/MGauge/StochEm.cc
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
Author: James Harrison <j.harrison@soton.ac.uk>
Author: Vera Guelpers <vmg1n14@soton.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
@ -59,24 +57,25 @@ std::vector<std::string> TStochEm::getOutput(void)
// setup ///////////////////////////////////////////////////////////////////////
void TStochEm::setup(void)
{
weightDone_ = env().hasCreatedObject("_" + getName() + "_weight");
envCacheLat(EmComp, "_" + getName() + "_weight");
if (!env().hasCreatedObject("_" + getName() + "_weight"))
{
envCacheLat(EmComp, "_" + getName() + "_weight");
}
envCreateLat(EmField, getName());
}
// execution ///////////////////////////////////////////////////////////////////
void TStochEm::execute(void)
{
LOG(Message) << "Generating stochastic EM potential..." << std::endl;
LOG(Message) << "Generating stochatic EM potential..." << std::endl;
std::vector<Real> improvements = strToVec<Real>(par().improvement);
PhotonR photon(par().gauge, par().zmScheme, improvements, par().G0_qedInf);
PhotonR photon(par().gauge, par().zmScheme);
auto &a = envGet(EmField, getName());
auto &w = envGet(EmComp, "_" + getName() + "_weight");
if (!weightDone_)
if (!env().hasCreatedObject("_" + getName() + "_weight"))
{
LOG(Message) << "Caching stochastic EM potential weight (gauge: "
LOG(Message) << "Caching stochatic EM potential weight (gauge: "
<< par().gauge << ", zero-mode scheme: "
<< par().zmScheme << ")..." << std::endl;
photon.StochasticWeight(w);

View File

@ -7,8 +7,6 @@ Source file: extras/Hadrons/Modules/MGauge/StochEm.hpp
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
Author: James Harrison <j.harrison@soton.ac.uk>
Author: Vera Guelpers <vmg1n14@soton.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
@ -46,9 +44,7 @@ class StochEmPar: Serializable
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(StochEmPar,
PhotonR::Gauge, gauge,
PhotonR::ZmScheme, zmScheme,
std::string, improvement,
Real, G0_qedInf);
PhotonR::ZmScheme, zmScheme);
};
class TStochEm: public Module<StochEmPar>
@ -60,7 +56,7 @@ public:
// constructor
TStochEm(const std::string name);
// destructor
virtual ~TStochEm(void) {};
virtual ~TStochEm(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
@ -69,11 +65,9 @@ protected:
virtual void setup(void);
// execution
virtual void execute(void);
private:
bool weightDone_;
};
MODULE_REGISTER(StochEm, TStochEm, MGauge);
MODULE_REGISTER_NS(StochEm, TStochEm, MGauge);
END_MODULE_NAMESPACE

View File

@ -1,7 +0,0 @@
#include <Grid/Hadrons/Modules/MGauge/StoutSmearing.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MGauge;
template class Grid::Hadrons::MGauge::TStoutSmearing<GIMPL>;

View File

@ -1,104 +0,0 @@
#ifndef Hadrons_MGauge_StoutSmearing_hpp_
#define Hadrons_MGauge_StoutSmearing_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* Stout smearing *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MGauge)
class StoutSmearingPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(StoutSmearingPar,
std::string, gauge,
unsigned int, steps,
double, rho);
};
template <typename GImpl>
class TStoutSmearing: public Module<StoutSmearingPar>
{
public:
typedef typename GImpl::Field GaugeField;
public:
// constructor
TStoutSmearing(const std::string name);
// destructor
virtual ~TStoutSmearing(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(StoutSmearing, TStoutSmearing<GIMPL>, MGauge);
/******************************************************************************
* TStoutSmearing implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename GImpl>
TStoutSmearing<GImpl>::TStoutSmearing(const std::string name)
: Module<StoutSmearingPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename GImpl>
std::vector<std::string> TStoutSmearing<GImpl>::getInput(void)
{
std::vector<std::string> in = {par().gauge};
return in;
}
template <typename GImpl>
std::vector<std::string> TStoutSmearing<GImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename GImpl>
void TStoutSmearing<GImpl>::setup(void)
{
envCreateLat(GaugeField, getName());
envTmpLat(GaugeField, "buf");
}
// execution ///////////////////////////////////////////////////////////////////
template <typename GImpl>
void TStoutSmearing<GImpl>::execute(void)
{
LOG(Message) << "Smearing '" << par().gauge << "' with " << par().steps
<< " step" << ((par().steps > 1) ? "s" : "")
<< " of stout smearing and rho= " << par().rho << std::endl;
Smear_Stout<GImpl> smearer(par().rho);
auto &U = envGet(GaugeField, par().gauge);
auto &Usmr = envGet(GaugeField, getName());
envGetTmp(GaugeField, buf);
buf = U;
for (unsigned int n = 0; n < par().steps; ++n)
{
smearer.smear(Usmr, buf);
buf = Usmr;
}
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MGauge_StoutSmearing_hpp_

View File

@ -46,7 +46,7 @@ public:
// constructor
TUnit(const std::string name);
// destructor
virtual ~TUnit(void) {};
virtual ~TUnit(void) = default;
// dependencies/products
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
@ -57,7 +57,7 @@ protected:
virtual void execute(void);
};
MODULE_REGISTER(Unit, TUnit, MGauge);
MODULE_REGISTER_NS(Unit, TUnit, MGauge);
END_MODULE_NAMESPACE

View File

@ -1,69 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MGauge/StochEm.cc
Copyright (C) 2015
Copyright (C) 2016
Author: James Harrison <j.harrison@soton.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 <Grid/Hadrons/Modules/MGauge/UnitEm.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MGauge;
/******************************************************************************
* TStochEm implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
TUnitEm::TUnitEm(const std::string name)
: Module<NoPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
std::vector<std::string> TUnitEm::getInput(void)
{
return std::vector<std::string>();
}
std::vector<std::string> TUnitEm::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
void TUnitEm::setup(void)
{
envCreateLat(EmField, getName());
}
// execution ///////////////////////////////////////////////////////////////////
void TUnitEm::execute(void)
{
PhotonR photon(0, 0); // Just chose arbitrary input values here
auto &a = envGet(EmField, getName());
LOG(Message) << "Generating unit EM potential..." << std::endl;
photon.UnitField(a);
}

View File

@ -1,69 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MGauge/StochEm.hpp
Copyright (C) 2015
Copyright (C) 2016
Author: James Harrison <j.harrison@soton.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_MGauge_UnitEm_hpp_
#define Hadrons_MGauge_UnitEm_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* StochEm *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MGauge)
class TUnitEm: public Module<NoPar>
{
public:
typedef PhotonR::GaugeField EmField;
typedef PhotonR::GaugeLinkField EmComp;
public:
// constructor
TUnitEm(const std::string name);
// destructor
virtual ~TUnitEm(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);
};
MODULE_REGISTER(UnitEm, TUnitEm, MGauge);
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MGauge_UnitEm_hpp_

View File

@ -1,40 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MIO/LoadBinary.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 <Grid/Hadrons/Modules/MIO/LoadBinary.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MIO;
template class Grid::Hadrons::MIO::TLoadBinary<GIMPL>;
template class Grid::Hadrons::MIO::TLoadBinary<ScalarNxNAdjImplR<2>>;
template class Grid::Hadrons::MIO::TLoadBinary<ScalarNxNAdjImplR<3>>;
template class Grid::Hadrons::MIO::TLoadBinary<ScalarNxNAdjImplR<4>>;
template class Grid::Hadrons::MIO::TLoadBinary<ScalarNxNAdjImplR<5>>;
template class Grid::Hadrons::MIO::TLoadBinary<ScalarNxNAdjImplR<6>>;

View File

@ -61,7 +61,7 @@ public:
// constructor
TLoadBinary(const std::string name);
// destructor
virtual ~TLoadBinary(void) {};
virtual ~TLoadBinary(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
@ -71,12 +71,12 @@ public:
virtual void execute(void);
};
MODULE_REGISTER_TMP(LoadBinary, TLoadBinary<GIMPL>, MIO);
MODULE_REGISTER_TMP(LoadBinaryScalarSU2, TLoadBinary<ScalarNxNAdjImplR<2>>, MIO);
MODULE_REGISTER_TMP(LoadBinaryScalarSU3, TLoadBinary<ScalarNxNAdjImplR<3>>, MIO);
MODULE_REGISTER_TMP(LoadBinaryScalarSU4, TLoadBinary<ScalarNxNAdjImplR<4>>, MIO);
MODULE_REGISTER_TMP(LoadBinaryScalarSU5, TLoadBinary<ScalarNxNAdjImplR<5>>, MIO);
MODULE_REGISTER_TMP(LoadBinaryScalarSU6, TLoadBinary<ScalarNxNAdjImplR<6>>, MIO);
MODULE_REGISTER_NS(LoadBinary, TLoadBinary<GIMPL>, MIO);
MODULE_REGISTER_NS(LoadBinaryScalarSU2, TLoadBinary<ScalarNxNAdjImplR<2>>, MIO);
MODULE_REGISTER_NS(LoadBinaryScalarSU3, TLoadBinary<ScalarNxNAdjImplR<3>>, MIO);
MODULE_REGISTER_NS(LoadBinaryScalarSU4, TLoadBinary<ScalarNxNAdjImplR<4>>, MIO);
MODULE_REGISTER_NS(LoadBinaryScalarSU5, TLoadBinary<ScalarNxNAdjImplR<5>>, MIO);
MODULE_REGISTER_NS(LoadBinaryScalarSU6, TLoadBinary<ScalarNxNAdjImplR<6>>, MIO);
/******************************************************************************
* TLoadBinary implementation *

View File

@ -1,35 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MIO/LoadCoarseEigenPack.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 <Grid/Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MIO;
template class Grid::Hadrons::MIO::TLoadCoarseEigenPack<CoarseFermionEigenPack<FIMPL,HADRONS_DEFAULT_LANCZOS_NBASIS>>;

View File

@ -1,135 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_MIO_LoadCoarseEigenPack_hpp_
#define Hadrons_MIO_LoadCoarseEigenPack_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
#include <Grid/Hadrons/EigenPack.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* Load local coherence eigen vectors/values package *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MIO)
class LoadCoarseEigenPackPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(LoadCoarseEigenPackPar,
std::string, filestem,
bool, multiFile,
unsigned int, sizeFine,
unsigned int, sizeCoarse,
unsigned int, Ls,
std::vector<int>, blockSize);
};
template <typename Pack>
class TLoadCoarseEigenPack: public Module<LoadCoarseEigenPackPar>
{
public:
typedef CoarseEigenPack<typename Pack::Field, typename Pack::CoarseField> BasePack;
template <typename vtype>
using iImplScalar = iScalar<iScalar<iScalar<vtype>>>;
typedef iImplScalar<typename Pack::Field::vector_type> SiteComplex;
public:
// constructor
TLoadCoarseEigenPack(const std::string name);
// destructor
virtual ~TLoadCoarseEigenPack(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(LoadCoarseFermionEigenPack, ARG(TLoadCoarseEigenPack<CoarseFermionEigenPack<FIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>>), MIO);
/******************************************************************************
* TLoadCoarseEigenPack implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename Pack>
TLoadCoarseEigenPack<Pack>::TLoadCoarseEigenPack(const std::string name)
: Module<LoadCoarseEigenPackPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename Pack>
std::vector<std::string> TLoadCoarseEigenPack<Pack>::getInput(void)
{
std::vector<std::string> in;
return in;
}
template <typename Pack>
std::vector<std::string> TLoadCoarseEigenPack<Pack>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename Pack>
void TLoadCoarseEigenPack<Pack>::setup(void)
{
env().createGrid(par().Ls);
env().createCoarseGrid(par().blockSize, par().Ls);
envCreateDerived(BasePack, Pack, getName(), par().Ls, par().sizeFine,
par().sizeCoarse, env().getRbGrid(par().Ls),
env().getCoarseGrid(par().blockSize, par().Ls));
}
// execution ///////////////////////////////////////////////////////////////////
template <typename Pack>
void TLoadCoarseEigenPack<Pack>::execute(void)
{
auto cg = env().getCoarseGrid(par().blockSize, par().Ls);
auto &epack = envGetDerived(BasePack, Pack, getName());
Lattice<SiteComplex> dummy(cg);
epack.read(par().filestem, par().multiFile, vm().getTrajectory());
LOG(Message) << "Block Gramm-Schmidt pass 1"<< std::endl;
blockOrthogonalise(dummy, epack.evec);
LOG(Message) << "Block Gramm-Schmidt pass 2"<< std::endl;
blockOrthogonalise(dummy, epack.evec);
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MIO_LoadCoarseEigenPack_hpp_

View File

@ -1,35 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MIO/LoadEigenPack.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 <Grid/Hadrons/Modules/MIO/LoadEigenPack.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MIO;
template class Grid::Hadrons::MIO::TLoadEigenPack<FermionEigenPack<FIMPL>>;

View File

@ -1,123 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MIO/LoadEigenPack.hpp
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_MIO_LoadEigenPack_hpp_
#define Hadrons_MIO_LoadEigenPack_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
#include <Grid/Hadrons/EigenPack.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* Load eigen vectors/values package *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MIO)
class LoadEigenPackPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(LoadEigenPackPar,
std::string, filestem,
bool, multiFile,
unsigned int, size,
unsigned int, Ls);
};
template <typename Pack>
class TLoadEigenPack: public Module<LoadEigenPackPar>
{
public:
typedef EigenPack<typename Pack::Field> BasePack;
public:
// constructor
TLoadEigenPack(const std::string name);
// destructor
virtual ~TLoadEigenPack(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(LoadFermionEigenPack, TLoadEigenPack<FermionEigenPack<FIMPL>>, MIO);
/******************************************************************************
* TLoadEigenPack implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename Pack>
TLoadEigenPack<Pack>::TLoadEigenPack(const std::string name)
: Module<LoadEigenPackPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename Pack>
std::vector<std::string> TLoadEigenPack<Pack>::getInput(void)
{
std::vector<std::string> in;
return in;
}
template <typename Pack>
std::vector<std::string> TLoadEigenPack<Pack>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename Pack>
void TLoadEigenPack<Pack>::setup(void)
{
env().createGrid(par().Ls);
envCreateDerived(BasePack, Pack, getName(), par().Ls, par().size,
env().getRbGrid(par().Ls));
}
// execution ///////////////////////////////////////////////////////////////////
template <typename Pack>
void TLoadEigenPack<Pack>::execute(void)
{
auto &epack = envGetDerived(BasePack, Pack, getName());
epack.read(par().filestem, par().multiFile, vm().getTrajectory());
epack.eval.resize(par().size);
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MIO_LoadEigenPack_hpp_

View File

@ -71,4 +71,6 @@ void TLoadNersc::execute(void)
auto &U = envGet(LatticeGaugeField, getName());
NerscIO::readConfiguration(U, header, fileName);
LOG(Message) << "NERSC header:" << std::endl;
dump_meta_data(header, LOG(Message));
}

View File

@ -52,7 +52,7 @@ public:
// constructor
TLoadNersc(const std::string name);
// destructor
virtual ~TLoadNersc(void) {};
virtual ~TLoadNersc(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
@ -62,7 +62,7 @@ public:
virtual void execute(void);
};
MODULE_REGISTER(LoadNersc, TLoadNersc, MIO);
MODULE_REGISTER_NS(LoadNersc, TLoadNersc, MIO);
END_MODULE_NAMESPACE

View File

@ -1,35 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MLoop/NoiseLoop.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 <Grid/Hadrons/Modules/MLoop/NoiseLoop.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MLoop;
template class Grid::Hadrons::MLoop::TNoiseLoop<FIMPL>;

View File

@ -71,7 +71,7 @@ public:
// constructor
TNoiseLoop(const std::string name);
// destructor
virtual ~TNoiseLoop(void) {};
virtual ~TNoiseLoop(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
@ -82,7 +82,7 @@ protected:
virtual void execute(void);
};
MODULE_REGISTER_TMP(NoiseLoop, TNoiseLoop<FIMPL>, MLoop);
MODULE_REGISTER_NS(NoiseLoop, TNoiseLoop<FIMPL>, MLoop);
/******************************************************************************
* TNoiseLoop implementation *

View File

@ -51,8 +51,7 @@ std::vector<std::string> TChargedProp::getInput(void)
std::vector<std::string> TChargedProp::getOutput(void)
{
std::vector<std::string> out = {getName(), getName()+"_0", getName()+"_Q",
getName()+"_Sun", getName()+"_Tad"};
std::vector<std::string> out = {getName()};
return out;
}
@ -67,27 +66,18 @@ void TChargedProp::setup(void)
phaseName_.push_back("_shiftphase_" + std::to_string(mu));
}
GFSrcName_ = getName() + "_DinvSrc";
prop0Name_ = getName() + "_0";
propQName_ = getName() + "_Q";
propSunName_ = getName() + "_Sun";
propTadName_ = getName() + "_Tad";
fftName_ = getName() + "_fft";
freeMomPropDone_ = env().hasCreatedObject(freeMomPropName_);
GFSrcDone_ = env().hasCreatedObject(GFSrcName_);
phasesDone_ = env().hasCreatedObject(phaseName_[0]);
prop0Done_ = env().hasCreatedObject(prop0Name_);
envCacheLat(ScalarField, freeMomPropName_);
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
envCacheLat(ScalarField, phaseName_[mu]);
}
envCacheLat(ScalarField, GFSrcName_);
envCacheLat(ScalarField, prop0Name_);
envCreateLat(ScalarField, getName());
envCreateLat(ScalarField, propQName_);
envCreateLat(ScalarField, propSunName_);
envCreateLat(ScalarField, propTadName_);
envTmpLat(ScalarField, "buf");
envTmpLat(ScalarField, "result");
envTmpLat(ScalarField, "Amu");
@ -105,125 +95,80 @@ void TChargedProp::execute(void)
<< " (mass= " << par().mass
<< ", charge= " << par().charge << ")..." << std::endl;
auto &prop = envGet(ScalarField, getName());
auto &prop0 = envGet(ScalarField, prop0Name_);
auto &propQ = envGet(ScalarField, propQName_);
auto &propSun = envGet(ScalarField, propSunName_);
auto &propTad = envGet(ScalarField, propTadName_);
auto &GFSrc = envGet(ScalarField, GFSrcName_);
auto &G = envGet(ScalarField, freeMomPropName_);
auto &fft = envGet(FFT, fftName_);
double q = par().charge;
envGetTmp(ScalarField, buf);
auto &prop = envGet(ScalarField, getName());
auto &GFSrc = envGet(ScalarField, GFSrcName_);
auto &G = envGet(ScalarField, freeMomPropName_);
auto &fft = envGet(FFT, fftName_);
double q = par().charge;
envGetTmp(ScalarField, result);
envGetTmp(ScalarField, buf);
// -G*momD1*G*F*Src (momD1 = F*D1*Finv)
propQ = GFSrc;
momD1(propQ, fft);
propQ = -G*propQ;
propSun = -propQ;
fft.FFT_dim(propQ, propQ, env().getNd()-1, FFT::backward);
// G*F*Src
prop = GFSrc;
// G*momD1*G*momD1*G*F*Src (here buf = G*momD1*G*F*Src)
momD1(propSun, fft);
propSun = G*propSun;
fft.FFT_dim(propSun, propSun, env().getNd()-1, FFT::backward);
// - q*G*momD1*G*F*Src (momD1 = F*D1*Finv)
buf = GFSrc;
momD1(buf, fft);
buf = G*buf;
prop = prop - q*buf;
// -G*momD2*G*F*Src (momD2 = F*D2*Finv)
propTad = GFSrc;
momD2(propTad, fft);
propTad = -G*propTad;
fft.FFT_dim(propTad, propTad, env().getNd()-1, FFT::backward);
// + q^2*G*momD1*G*momD1*G*F*Src (here buf = G*momD1*G*F*Src)
momD1(buf, fft);
prop = prop + q*q*G*buf;
// - q^2*G*momD2*G*F*Src (momD2 = F*D2*Finv)
buf = GFSrc;
momD2(buf, fft);
prop = prop - q*q*G*buf;
// final FT
fft.FFT_all_dim(prop, prop, FFT::backward);
// full charged scalar propagator
fft.FFT_dim(buf, GFSrc, env().getNd()-1, FFT::backward);
prop = buf + q*propQ + q*q*propSun + q*q*propTad;
// OUTPUT IF NECESSARY
if (!par().output.empty())
{
Result result;
TComplex site;
std::vector<int> siteCoor;
LOG(Message) << "Saving momentum-projected propagator to '"
<< RESULT_FILE_NAME(par().output) << "'..."
<< std::endl;
result.projection.resize(par().outputMom.size());
result.lattice_size = env().getGrid()->_fdimensions;
result.mass = par().mass;
result.charge = q;
siteCoor.resize(env().getNd());
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
std::string filename = par().output + "." +
std::to_string(vm().getTrajectory());
LOG(Message) << "Saving zero-momentum projection to '"
<< filename << "'..." << std::endl;
ResultWriter writer(RESULT_FILE_NAME(par().output));
std::vector<TComplex> vecBuf;
std::vector<Complex> result;
sliceSum(prop, vecBuf, Tp);
result.resize(vecBuf.size());
for (unsigned int t = 0; t < vecBuf.size(); ++t)
{
result.projection[i_p].momentum = strToVec<int>(par().outputMom[i_p]);
LOG(Message) << "Calculating (" << par().outputMom[i_p]
<< ") momentum projection" << std::endl;
result.projection[i_p].corr_0.resize(env().getGrid()->_fdimensions[env().getNd()-1]);
result.projection[i_p].corr.resize(env().getGrid()->_fdimensions[env().getNd()-1]);
result.projection[i_p].corr_Q.resize(env().getGrid()->_fdimensions[env().getNd()-1]);
result.projection[i_p].corr_Sun.resize(env().getGrid()->_fdimensions[env().getNd()-1]);
result.projection[i_p].corr_Tad.resize(env().getGrid()->_fdimensions[env().getNd()-1]);
for (unsigned int j = 0; j < env().getNd()-1; ++j)
{
siteCoor[j] = result.projection[i_p].momentum[j];
}
for (unsigned int t = 0; t < result.projection[i_p].corr.size(); ++t)
{
siteCoor[env().getNd()-1] = t;
peekSite(site, prop, siteCoor);
result.projection[i_p].corr[t]=TensorRemove(site);
peekSite(site, buf, siteCoor);
result.projection[i_p].corr_0[t]=TensorRemove(site);
peekSite(site, propQ, siteCoor);
result.projection[i_p].corr_Q[t]=TensorRemove(site);
peekSite(site, propSun, siteCoor);
result.projection[i_p].corr_Sun[t]=TensorRemove(site);
peekSite(site, propTad, siteCoor);
result.projection[i_p].corr_Tad[t]=TensorRemove(site);
}
result[t] = TensorRemove(vecBuf[t]);
}
saveResult(par().output, "prop", result);
write(writer, "charge", q);
write(writer, "prop", result);
}
std::vector<int> mask(env().getNd(),1);
mask[env().getNd()-1] = 0;
fft.FFT_dim_mask(prop, prop, mask, FFT::backward);
fft.FFT_dim_mask(propQ, propQ, mask, FFT::backward);
fft.FFT_dim_mask(propSun, propSun, mask, FFT::backward);
fft.FFT_dim_mask(propTad, propTad, mask, FFT::backward);
}
void TChargedProp::makeCaches(void)
{
auto &freeMomProp = envGet(ScalarField, freeMomPropName_);
auto &GFSrc = envGet(ScalarField, GFSrcName_);
auto &prop0 = envGet(ScalarField, prop0Name_);
auto &fft = envGet(FFT, fftName_);
if (!freeMomPropDone_)
{
LOG(Message) << "Caching momentum-space free scalar propagator"
LOG(Message) << "Caching momentum space free scalar propagator"
<< " (mass= " << par().mass << ")..." << std::endl;
SIMPL::MomentumSpacePropagator(freeMomProp, par().mass);
}
if (!GFSrcDone_)
{
FFT fft(env().getGrid());
auto &source = envGet(ScalarField, par().source);
LOG(Message) << "Caching G*F*src..." << std::endl;
fft.FFT_all_dim(GFSrc, source, FFT::forward);
GFSrc = freeMomProp*GFSrc;
}
if (!prop0Done_)
{
LOG(Message) << "Caching position-space free scalar propagator..."
<< std::endl;
fft.FFT_all_dim(prop0, GFSrc, FFT::backward);
}
if (!phasesDone_)
{
std::vector<int> &l = env().getGrid()->_fdimensions;
@ -240,14 +185,6 @@ void TChargedProp::makeCaches(void)
phase_.push_back(&phmu);
}
}
else
{
phase_.clear();
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
phase_.push_back(env().getObject<ScalarField>(phaseName_[mu]));
}
}
}
void TChargedProp::momD1(ScalarField &s, FFT &fft)

View File

@ -7,7 +7,6 @@ Source file: extras/Hadrons/Modules/MScalar/ChargedProp.hpp
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
Author: James Harrison <j.harrison@soton.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
@ -48,8 +47,7 @@ public:
std::string, source,
double, mass,
double, charge,
std::string, output,
std::vector<std::string>, outputMom);
std::string, output);
};
class TChargedProp: public Module<ChargedPropPar>
@ -58,31 +56,11 @@ public:
SCALAR_TYPE_ALIASES(SIMPL,);
typedef PhotonR::GaugeField EmField;
typedef PhotonR::GaugeLinkField EmComp;
class Result: Serializable
{
public:
class Projection: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(Projection,
std::vector<int>, momentum,
std::vector<Complex>, corr,
std::vector<Complex>, corr_0,
std::vector<Complex>, corr_Q,
std::vector<Complex>, corr_Sun,
std::vector<Complex>, corr_Tad);
};
GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
std::vector<int>, lattice_size,
double, mass,
double, charge,
std::vector<Projection>, projection);
};
public:
// constructor
TChargedProp(const std::string name);
// destructor
virtual ~TChargedProp(void) {};
virtual ~TChargedProp(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
@ -96,15 +74,13 @@ private:
void momD1(ScalarField &s, FFT &fft);
void momD2(ScalarField &s, FFT &fft);
private:
bool freeMomPropDone_, GFSrcDone_, prop0Done_,
phasesDone_;
std::string freeMomPropName_, GFSrcName_, prop0Name_,
propQName_, propSunName_, propTadName_, fftName_;
bool freeMomPropDone_, GFSrcDone_, phasesDone_;
std::string freeMomPropName_, GFSrcName_, fftName_;
std::vector<std::string> phaseName_;
std::vector<ScalarField *> phase_;
};
MODULE_REGISTER(ChargedProp, TChargedProp, MScalar);
MODULE_REGISTER_NS(ChargedProp, TChargedProp, MScalar);
END_MODULE_NAMESPACE

View File

@ -83,6 +83,8 @@ void TFreeProp::execute(void)
if (!par().output.empty())
{
TextWriter writer(par().output + "." +
std::to_string(vm().getTrajectory()));
std::vector<TComplex> buf;
std::vector<Complex> result;
@ -92,6 +94,6 @@ void TFreeProp::execute(void)
{
result[t] = TensorRemove(buf[t]);
}
saveResult(par().output, "freeprop", result);
write(writer, "prop", result);
}
}

View File

@ -56,7 +56,7 @@ public:
// constructor
TFreeProp(const std::string name);
// destructor
virtual ~TFreeProp(void) {};
virtual ~TFreeProp(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
@ -70,7 +70,7 @@ private:
bool freePropDone_;
};
MODULE_REGISTER(FreeProp, TFreeProp, MScalar);
MODULE_REGISTER_NS(FreeProp, TFreeProp, MScalar);
END_MODULE_NAMESPACE

View File

@ -1,536 +0,0 @@
#include <Grid/Hadrons/Modules/MScalar/ChargedProp.hpp>
#include <Grid/Hadrons/Modules/MScalar/ScalarVP.hpp>
#include <Grid/Hadrons/Modules/MScalar/Scalar.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MScalar;
/*
* Scalar QED vacuum polarisation up to O(alpha)
*
* Conserved vector 2-point function diagram notation:
* _______
* / \
* U_nu * * U_mu
* \_______/
*
* ( adj(S(a\hat{nu}|x)) U_mu(x) S(0|x+a\hat{mu}) U_nu(0) )
* = 2 Re( - )
* ( adj(S(a\hat{nu}|x+a\hat{mu})) adj(U_mu(x)) S(0|x) U_nu(0) )
*
*
* _______
* / \
* free = 1 * * 1
* \_______/
*
*
*
* _______
* / \
* S = iA_nu * * iA_mu
* \_______/
*
*
* Delta_1
* ___*___
* / \
* X = 1 * * 1
* \___*___/
* Delta_1
*
* Delta_1 Delta_1
* ___*___ ___*___
* / \ / \
* 1 * * iA_mu + iA_nu * * 1
* \_______/ \_______/
* 4C = _______ _______
* / \ / \
* + 1 * * iA_mu + iA_nu * * 1
* \___*___/ \___*___/
* Delta_1 Delta_1
*
* Delta_1 Delta_1
* _*___*_ _______
* / \ / \
* 2E = 1 * * 1 + 1 * * 1
* \_______/ \_*___*_/
* Delta_1 Delta_1
*
* Delta_2
* ___*___ _______
* / \ / \
* 2T = 1 * * 1 + 1 * * 1
* \_______/ \___*___/
* Delta_2
*
*
* _______
* / \
* srcT = -A_nu^2/2 * * 1
* \_______/
*
*
*
* _______
* / \
* snkT = 1 * * -A_mu^2/2
* \_______/
*
* Full VP to O(alpha) = free + q^2*(S+X+4C+2E+2T+srcT+snkT)
*/
/******************************************************************************
* TScalarVP implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
TScalarVP::TScalarVP(const std::string name)
: Module<ScalarVPPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
std::vector<std::string> TScalarVP::getInput(void)
{
prop0Name_ = par().scalarProp + "_0";
propQName_ = par().scalarProp + "_Q";
propSunName_ = par().scalarProp + "_Sun";
propTadName_ = par().scalarProp + "_Tad";
std::vector<std::string> in = {par().emField, prop0Name_, propQName_,
propSunName_, propTadName_};
return in;
}
std::vector<std::string> TScalarVP::getOutput(void)
{
std::vector<std::string> out;
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
// out.push_back(getName() + "_propQ_" + std::to_string(mu));
for (unsigned int nu = 0; nu < env().getNd(); ++nu)
{
out.push_back(getName() + "_" + std::to_string(mu)
+ "_" + std::to_string(nu));
}
}
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
void TScalarVP::setup(void)
{
freeMomPropName_ = FREEMOMPROP(static_cast<TChargedProp *>(vm().getModule(par().scalarProp))->par().mass);
GFSrcName_ = par().scalarProp + "_DinvSrc";
fftName_ = par().scalarProp + "_fft";
phaseName_.clear();
muPropQName_.clear();
vpTensorName_.clear();
momPhaseName_.clear();
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
phaseName_.push_back("_shiftphase_" + std::to_string(mu));
muPropQName_.push_back(getName() + "_propQ_" + std::to_string(mu));
std::vector<std::string> vpTensorName_mu;
for (unsigned int nu = 0; nu < env().getNd(); ++nu)
{
vpTensorName_mu.push_back(getName() + "_" + std::to_string(mu)
+ "_" + std::to_string(nu));
}
vpTensorName_.push_back(vpTensorName_mu);
}
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
momPhaseName_.push_back("_momentumphase_" + std::to_string(i_p));
}
}
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
envCreateLat(ScalarField, muPropQName_[mu]);
for (unsigned int nu = 0; nu < env().getNd(); ++nu)
{
envCreateLat(ScalarField, vpTensorName_[mu][nu]);
}
}
if (!par().output.empty())
{
momPhasesDone_ = env().hasCreatedObject(momPhaseName_[0]);
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
envCacheLat(ScalarField, momPhaseName_[i_p]);
}
}
envTmpLat(ScalarField, "buf");
envTmpLat(ScalarField, "result");
envTmpLat(ScalarField, "Amu");
envTmpLat(ScalarField, "Usnk");
envTmpLat(ScalarField, "tmpProp");
}
// execution ///////////////////////////////////////////////////////////////////
void TScalarVP::execute(void)
{
// CACHING ANALYTIC EXPRESSIONS
makeCaches();
Complex ci(0.0,1.0);
Real q = static_cast<TChargedProp *>(vm().getModule(par().scalarProp))->par().charge;
auto &prop0 = envGet(ScalarField, prop0Name_);
auto &propQ = envGet(ScalarField, propQName_);
auto &propSun = envGet(ScalarField, propSunName_);
auto &propTad = envGet(ScalarField, propTadName_);
auto &GFSrc = envGet(ScalarField, GFSrcName_);
auto &G = envGet(ScalarField, freeMomPropName_);
auto &fft = envGet(FFT, fftName_);
phase_.clear();
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
auto &phmu = envGet(ScalarField, phaseName_[mu]);
phase_.push_back(&phmu);
}
// PROPAGATORS FROM SHIFTED SOURCES
LOG(Message) << "Computing O(q) charged scalar propagators..."
<< std::endl;
std::vector<ScalarField *> muPropQ;
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
auto &propmu = envGet(ScalarField, muPropQName_[mu]);
// -G*momD1*G*F*tau_mu*Src (momD1 = F*D1*Finv)
propmu = adj(*phase_[mu])*GFSrc;
momD1(propmu, fft);
propmu = -G*propmu;
fft.FFT_all_dim(propmu, propmu, FFT::backward);
muPropQ.push_back(&propmu);
}
// CONTRACTIONS
auto &A = envGet(EmField, par().emField);
envGetTmp(ScalarField, buf);
envGetTmp(ScalarField, result);
envGetTmp(ScalarField, Amu);
envGetTmp(ScalarField, Usnk);
envGetTmp(ScalarField, tmpProp);
TComplex Anu0, Usrc;
std::vector<int> coor0 = {0, 0, 0, 0};
std::vector<std::vector<ScalarField *> > vpTensor;
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
std::vector<ScalarField *> vpTensor_mu;
for (unsigned int nu = 0; nu < env().getNd(); ++nu)
{
auto &vpmunu = envGet(ScalarField, vpTensorName_[mu][nu]);
vpTensor_mu.push_back(&vpmunu);
}
vpTensor.push_back(vpTensor_mu);
}
// Prepare output data structure if necessary
Result outputData;
if (!par().output.empty())
{
outputData.projection.resize(par().outputMom.size());
outputData.lattice_size = env().getGrid()->_fdimensions;
outputData.mass = static_cast<TChargedProp *>(vm().getModule(par().scalarProp))->par().mass;
outputData.charge = q;
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
outputData.projection[i_p].momentum = strToVec<int>(par().outputMom[i_p]);
outputData.projection[i_p].pi.resize(env().getNd());
outputData.projection[i_p].pi_free.resize(env().getNd());
outputData.projection[i_p].pi_2E.resize(env().getNd());
outputData.projection[i_p].pi_2T.resize(env().getNd());
outputData.projection[i_p].pi_S.resize(env().getNd());
outputData.projection[i_p].pi_4C.resize(env().getNd());
outputData.projection[i_p].pi_X.resize(env().getNd());
outputData.projection[i_p].pi_srcT.resize(env().getNd());
outputData.projection[i_p].pi_snkT.resize(env().getNd());
for (unsigned int nu = 0; nu < env().getNd(); ++nu)
{
outputData.projection[i_p].pi[nu].resize(env().getNd());
outputData.projection[i_p].pi_free[nu].resize(env().getNd());
outputData.projection[i_p].pi_2E[nu].resize(env().getNd());
outputData.projection[i_p].pi_2T[nu].resize(env().getNd());
outputData.projection[i_p].pi_S[nu].resize(env().getNd());
outputData.projection[i_p].pi_4C[nu].resize(env().getNd());
outputData.projection[i_p].pi_X[nu].resize(env().getNd());
outputData.projection[i_p].pi_srcT[nu].resize(env().getNd());
outputData.projection[i_p].pi_snkT[nu].resize(env().getNd());
}
}
}
// Do contractions
for (unsigned int nu = 0; nu < env().getNd(); ++nu)
{
peekSite(Anu0, peekLorentz(A, nu), coor0);
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
LOG(Message) << "Computing Pi[" << mu << "][" << nu << "]..."
<< std::endl;
Amu = peekLorentz(A, mu);
// free
tmpProp = Cshift(prop0, nu, -1); // S_0(0|x-a\hat{\nu})
// = S_0(a\hat{\nu}|x)
Usrc = Complex(1.0,0.0);
vpContraction(result, prop0, tmpProp, Usrc, mu);
*vpTensor[mu][nu] = result;
// Do momentum projections if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].pi_free[mu][nu], result,
i_p);
}
}
tmpProp = result; // Just using tmpProp as a temporary ScalarField
// here (buf is modified by calls to writeVP())
// srcT
result = tmpProp * (-0.5)*Anu0*Anu0;
*vpTensor[mu][nu] += q*q*result;
// Do momentum projections if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].pi_srcT[mu][nu], result,
i_p);
}
}
// snkT
result = tmpProp * (-0.5)*Amu*Amu;
*vpTensor[mu][nu] += q*q*result;
// Do momentum projections if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].pi_snkT[mu][nu], result,
i_p);
}
}
// S
tmpProp = Cshift(prop0, nu, -1); // S_0(a\hat{\nu}|x)
Usrc = ci*Anu0;
Usnk = ci*Amu;
vpContraction(result, prop0, tmpProp, Usrc, Usnk, mu);
*vpTensor[mu][nu] += q*q*result;
// Do momentum projections if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].pi_S[mu][nu], result,
i_p);
}
}
// 4C
tmpProp = Cshift(prop0, nu, -1); // S_0(a\hat{\nu}|x)
Usrc = Complex(1.0,0.0);
Usnk = ci*Amu;
vpContraction(result, propQ, tmpProp, Usrc, Usnk, mu);
Usrc = ci*Anu0;
vpContraction(buf, propQ, tmpProp, Usrc, mu);
result += buf;
vpContraction(buf, prop0, *muPropQ[nu], Usrc, mu);
result += buf;
Usrc = Complex(1.0,0.0);
Usnk = ci*Amu;
vpContraction(buf, prop0, *muPropQ[nu], Usrc, Usnk, mu);
result += buf;
*vpTensor[mu][nu] += q*q*result;
// Do momentum projections if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].pi_4C[mu][nu], result,
i_p);
}
}
// X
Usrc = Complex(1.0,0.0);
vpContraction(result, propQ, *muPropQ[nu], Usrc, mu);
*vpTensor[mu][nu] += q*q*result;
// Do momentum projections if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].pi_X[mu][nu], result,
i_p);
}
}
// 2E
tmpProp = Cshift(prop0, nu, -1); // S_0(a\hat{\nu}|x)
Usrc = Complex(1.0,0.0);
vpContraction(result, propSun, tmpProp, Usrc, mu);
tmpProp = Cshift(propSun, nu, -1); // S_\Sigma(0|x-a\hat{\nu})
//(Note: <S(0|x-a\hat{\nu})> = <S(a\hat{\nu}|x)>)
vpContraction(buf, prop0, tmpProp, Usrc, mu);
result += buf;
*vpTensor[mu][nu] += q*q*result;
// Do momentum projections if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].pi_2E[mu][nu], result,
i_p);
}
}
// 2T
tmpProp = Cshift(prop0, nu, -1); // S_0(a\hat{\nu}|x)
Usrc = Complex(1.0,0.0);
vpContraction(result, propTad, tmpProp, Usrc, mu);
tmpProp = Cshift(propTad, nu, -1); // S_T(0|x-a\hat{\nu})
vpContraction(buf, prop0, tmpProp, Usrc, mu);
result += buf;
*vpTensor[mu][nu] += q*q*result;
// Do momentum projections if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].pi_2T[mu][nu], result,
i_p);
}
}
// Do momentum projections of full VP if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].pi[mu][nu],
*vpTensor[mu][nu], i_p);
}
}
}
}
// OUTPUT IF NECESSARY
if (!par().output.empty())
{
LOG(Message) << "Saving momentum-projected HVP to '"
<< RESULT_FILE_NAME(par().output) << "'..."
<< std::endl;
saveResult(par().output, "HVP", outputData);
}
}
void TScalarVP::makeCaches(void)
{
envGetTmp(ScalarField, buf);
if ( (!par().output.empty()) && (!momPhasesDone_) )
{
LOG(Message) << "Caching phases for momentum projections..."
<< std::endl;
std::vector<int> &l = env().getGrid()->_fdimensions;
Complex ci(0.0,1.0);
// Calculate phase factors
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
std::vector<int> mom = strToVec<int>(par().outputMom[i_p]);
auto &momph_ip = envGet(ScalarField, momPhaseName_[i_p]);
momph_ip = zero;
for (unsigned int j = 0; j < env().getNd()-1; ++j)
{
Real twoPiL = M_PI*2./l[j];
LatticeCoordinate(buf, j);
buf = mom[j]*twoPiL*buf;
momph_ip = momph_ip + buf;
}
momph_ip = exp(-ci*momph_ip);
momPhase_.push_back(&momph_ip);
}
}
}
void TScalarVP::vpContraction(ScalarField &vp,
ScalarField &prop_0_x, ScalarField &prop_nu_x,
TComplex u_src, ScalarField &u_snk, int mu)
{
// Note: this function assumes a point source is used.
vp = adj(prop_nu_x) * u_snk * Cshift(prop_0_x, mu, 1) * u_src;
vp -= Cshift(adj(prop_nu_x), mu, 1) * adj(u_snk) * prop_0_x * u_src;
vp = 2.0*real(vp);
}
void TScalarVP::vpContraction(ScalarField &vp,
ScalarField &prop_0_x, ScalarField &prop_nu_x,
TComplex u_src, int mu)
{
// Note: this function assumes a point source is used.
vp = adj(prop_nu_x) * Cshift(prop_0_x, mu, 1) * u_src;
vp -= Cshift(adj(prop_nu_x), mu, 1) * prop_0_x * u_src;
vp = 2.0*real(vp);
}
void TScalarVP::project(std::vector<Complex> &projection, const ScalarField &vp, int i_p)
{
std::vector<TComplex> vecBuf;
envGetTmp(ScalarField, buf);
buf = vp*(*momPhase_[i_p]);
sliceSum(buf, vecBuf, Tp);
projection.resize(vecBuf.size());
for (unsigned int t = 0; t < vecBuf.size(); ++t)
{
projection[t] = TensorRemove(vecBuf[t]);
}
}
void TScalarVP::momD1(ScalarField &s, FFT &fft)
{
auto &A = envGet(EmField, par().emField);
Complex ci(0.0,1.0);
envGetTmp(ScalarField, buf);
envGetTmp(ScalarField, result);
envGetTmp(ScalarField, Amu);
result = zero;
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
Amu = peekLorentz(A, mu);
buf = (*phase_[mu])*s;
fft.FFT_all_dim(buf, buf, FFT::backward);
buf = Amu*buf;
fft.FFT_all_dim(buf, buf, FFT::forward);
result = result - ci*buf;
}
fft.FFT_all_dim(s, s, FFT::backward);
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
Amu = peekLorentz(A, mu);
buf = Amu*s;
fft.FFT_all_dim(buf, buf, FFT::forward);
result = result + ci*adj(*phase_[mu])*buf;
}
s = result;
}

View File

@ -1,101 +0,0 @@
#ifndef Hadrons_MScalar_ScalarVP_hpp_
#define Hadrons_MScalar_ScalarVP_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* Scalar vacuum polarisation *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MScalar)
class ScalarVPPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(ScalarVPPar,
std::string, emField,
std::string, scalarProp,
std::string, output,
std::vector<std::string>, outputMom);
};
class TScalarVP: public Module<ScalarVPPar>
{
public:
SCALAR_TYPE_ALIASES(SIMPL,);
typedef PhotonR::GaugeField EmField;
typedef PhotonR::GaugeLinkField EmComp;
class Result: Serializable
{
public:
class Projection: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(Projection,
std::vector<int>, momentum,
std::vector<std::vector<std::vector<Complex>>>, pi,
std::vector<std::vector<std::vector<Complex>>>, pi_free,
std::vector<std::vector<std::vector<Complex>>>, pi_2E,
std::vector<std::vector<std::vector<Complex>>>, pi_2T,
std::vector<std::vector<std::vector<Complex>>>, pi_S,
std::vector<std::vector<std::vector<Complex>>>, pi_4C,
std::vector<std::vector<std::vector<Complex>>>, pi_X,
std::vector<std::vector<std::vector<Complex>>>, pi_srcT,
std::vector<std::vector<std::vector<Complex>>>, pi_snkT);
};
GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
std::vector<int>, lattice_size,
double, mass,
double, charge,
std::vector<Projection>, projection);
};
public:
// constructor
TScalarVP(const std::string name);
// destructor
virtual ~TScalarVP(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);
private:
void makeCaches(void);
// conserved vector two-point contraction
void vpContraction(ScalarField &vp,
ScalarField &prop_0_x, ScalarField &prop_nu_x,
TComplex u_src, ScalarField &u_snk, int mu);
// conserved vector two-point contraction with unit gauge link at sink
void vpContraction(ScalarField &vp,
ScalarField &prop_0_x, ScalarField &prop_nu_x,
TComplex u_src, int mu);
// write momentum-projected vacuum polarisation to file(s)
void project(std::vector<Complex> &projection, const ScalarField &vp,
int i_p);
// momentum-space Delta_1 insertion
void momD1(ScalarField &s, FFT &fft);
private:
bool momPhasesDone_;
std::string freeMomPropName_, GFSrcName_,
prop0Name_, propQName_,
propSunName_, propTadName_,
fftName_;
std::vector<std::string> phaseName_, muPropQName_,
momPhaseName_;
std::vector<std::vector<std::string> > vpTensorName_;
std::vector<ScalarField *> phase_, momPhase_;
};
MODULE_REGISTER(ScalarVP, TScalarVP, MScalar);
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MScalar_ScalarVP_hpp_

View File

@ -1,232 +0,0 @@
#include <Grid/Hadrons/Modules/MScalar/VPCounterTerms.hpp>
#include <Grid/Hadrons/Modules/MScalar/Scalar.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MScalar;
/******************************************************************************
* TVPCounterTerms implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
TVPCounterTerms::TVPCounterTerms(const std::string name)
: Module<VPCounterTermsPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
std::vector<std::string> TVPCounterTerms::getInput(void)
{
std::vector<std::string> in = {par().source};
return in;
}
std::vector<std::string> TVPCounterTerms::getOutput(void)
{
std::vector<std::string> out;
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
void TVPCounterTerms::setup(void)
{
freeMomPropName_ = FREEMOMPROP(par().mass);
phaseName_.clear();
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
phaseName_.push_back("_shiftphase_" + std::to_string(mu));
}
GFSrcName_ = getName() + "_DinvSrc";
phatsqName_ = getName() + "_pHatSquared";
prop0Name_ = getName() + "_freeProp";
twoscalarName_ = getName() + "_2scalarProp";
psquaredName_ = getName() + "_psquaredProp";
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
momPhaseName_.push_back("_momentumphase_" + std::to_string(i_p));
}
}
envCreateLat(ScalarField, freeMomPropName_);
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
envCreateLat(ScalarField, phaseName_[mu]);
}
envCreateLat(ScalarField, phatsqName_);
envCreateLat(ScalarField, GFSrcName_);
envCreateLat(ScalarField, prop0Name_);
envCreateLat(ScalarField, twoscalarName_);
envCreateLat(ScalarField, psquaredName_);
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
envCacheLat(ScalarField, momPhaseName_[i_p]);
}
}
envTmpLat(ScalarField, "buf");
envTmpLat(ScalarField, "tmp_vp");
envTmpLat(ScalarField, "vpPhase");
}
// execution ///////////////////////////////////////////////////////////////////
void TVPCounterTerms::execute(void)
{
auto &source = envGet(ScalarField, par().source);
Complex ci(0.0,1.0);
FFT fft(env().getGrid());
envGetTmp(ScalarField, buf);
envGetTmp(ScalarField, tmp_vp);
// Momentum-space free scalar propagator
auto &G = envGet(ScalarField, freeMomPropName_);
SIMPL::MomentumSpacePropagator(G, par().mass);
// Phases and hat{p}^2
auto &phatsq = envGet(ScalarField, phatsqName_);
std::vector<int> &l = env().getGrid()->_fdimensions;
LOG(Message) << "Calculating shift phases..." << std::endl;
phatsq = zero;
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
Real twoPiL = M_PI*2./l[mu];
auto &phmu = envGet(ScalarField, phaseName_[mu]);
LatticeCoordinate(buf, mu);
phmu = exp(ci*twoPiL*buf);
phase_.push_back(&phmu);
buf = 2.*sin(.5*twoPiL*buf);
phatsq = phatsq + buf*buf;
}
// G*F*src
auto &GFSrc = envGet(ScalarField, GFSrcName_);
fft.FFT_all_dim(GFSrc, source, FFT::forward);
GFSrc = G*GFSrc;
// Position-space free scalar propagator
auto &prop0 = envGet(ScalarField, prop0Name_);
prop0 = GFSrc;
fft.FFT_all_dim(prop0, prop0, FFT::backward);
// Propagators for counter-terms
auto &twoscalarProp = envGet(ScalarField, twoscalarName_);
auto &psquaredProp = envGet(ScalarField, psquaredName_);
twoscalarProp = G*GFSrc;
fft.FFT_all_dim(twoscalarProp, twoscalarProp, FFT::backward);
psquaredProp = G*phatsq*GFSrc;
fft.FFT_all_dim(psquaredProp, psquaredProp, FFT::backward);
// Prepare output data structure if necessary
Result outputData;
if (!par().output.empty())
{
outputData.projection.resize(par().outputMom.size());
outputData.lattice_size = env().getGrid()->_fdimensions;
outputData.mass = par().mass;
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
outputData.projection[i_p].momentum = strToVec<int>(par().outputMom[i_p]);
outputData.projection[i_p].twoScalar.resize(env().getNd());
outputData.projection[i_p].threeScalar.resize(env().getNd());
outputData.projection[i_p].pSquaredInsertion.resize(env().getNd());
for (unsigned int nu = 0; nu < env().getNd(); ++nu)
{
outputData.projection[i_p].twoScalar[nu].resize(env().getNd());
outputData.projection[i_p].threeScalar[nu].resize(env().getNd());
outputData.projection[i_p].pSquaredInsertion[nu].resize(env().getNd());
}
// Calculate phase factors
auto &momph_ip = envGet(ScalarField, momPhaseName_[i_p]);
momph_ip = zero;
for (unsigned int j = 0; j < env().getNd()-1; ++j)
{
Real twoPiL = M_PI*2./l[j];
LatticeCoordinate(buf, j);
buf = outputData.projection[i_p].momentum[j]*twoPiL*buf;
momph_ip = momph_ip + buf;
}
momph_ip = exp(-ci*momph_ip);
momPhase_.push_back(&momph_ip);
}
}
// Contractions
for (unsigned int nu = 0; nu < env().getNd(); ++nu)
{
buf = adj(Cshift(prop0, nu, -1));
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
// Two-scalar loop
tmp_vp = buf * Cshift(prop0, mu, 1);
tmp_vp -= Cshift(buf, mu, 1) * prop0;
tmp_vp = 2.0*real(tmp_vp);
// Output if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].twoScalar[mu][nu],
tmp_vp, i_p);
}
}
// Three-scalar loop (no vertex)
tmp_vp = buf * Cshift(twoscalarProp, mu, 1);
tmp_vp -= Cshift(buf, mu, 1) * twoscalarProp;
tmp_vp = 2.0*real(tmp_vp);
// Output if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].threeScalar[mu][nu],
tmp_vp, i_p);
}
}
// Three-scalar loop (hat{p}^2 insertion)
tmp_vp = buf * Cshift(psquaredProp, mu, 1);
tmp_vp -= Cshift(buf, mu, 1) * psquaredProp;
tmp_vp = 2.0*real(tmp_vp);
// Output if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].pSquaredInsertion[mu][nu],
tmp_vp, i_p);
}
}
}
}
// OUTPUT IF NECESSARY
if (!par().output.empty())
{
LOG(Message) << "Saving momentum-projected correlators to '"
<< RESULT_FILE_NAME(par().output) << "'..."
<< std::endl;
saveResult(par().output, "scalar_loops", outputData);
}
}
void TVPCounterTerms::project(std::vector<Complex> &projection, const ScalarField &vp, int i_p)
{
std::vector<TComplex> vecBuf;
envGetTmp(ScalarField, vpPhase);
vpPhase = vp*(*momPhase_[i_p]);
sliceSum(vpPhase, vecBuf, Tp);
projection.resize(vecBuf.size());
for (unsigned int t = 0; t < vecBuf.size(); ++t)
{
projection[t] = TensorRemove(vecBuf[t]);
}
}

View File

@ -1,75 +0,0 @@
#ifndef Hadrons_MScalar_VPCounterTerms_hpp_
#define Hadrons_MScalar_VPCounterTerms_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* VPCounterTerms *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MScalar)
class VPCounterTermsPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(VPCounterTermsPar,
std::string, source,
double, mass,
std::string, output,
std::vector<std::string>, outputMom);
};
class TVPCounterTerms: public Module<VPCounterTermsPar>
{
public:
SCALAR_TYPE_ALIASES(SIMPL,);
class Result: Serializable
{
public:
class Projection: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(Projection,
std::vector<int>, momentum,
std::vector<std::vector<std::vector<Complex>>>, twoScalar,
std::vector<std::vector<std::vector<Complex>>>, threeScalar,
std::vector<std::vector<std::vector<Complex>>>, pSquaredInsertion);
};
GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
std::vector<int>, lattice_size,
double, mass,
std::vector<Projection>, projection);
};
public:
// constructor
TVPCounterTerms(const std::string name);
// destructor
virtual ~TVPCounterTerms(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);
private:
void project(std::vector<Complex> &projection, const ScalarField &vp, int i_p);
private:
std::string freeMomPropName_, GFSrcName_, phatsqName_, prop0Name_,
twoscalarName_, twoscalarVertexName_,
psquaredName_, psquaredVertexName_;
std::vector<std::string> phaseName_, momPhaseName_;
std::vector<ScalarField *> phase_, momPhase_;
};
MODULE_REGISTER(VPCounterTerms, TVPCounterTerms, MScalar);
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MScalar_VPCounterTerms_hpp_

View File

@ -1,39 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MScalarSUN/Div.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 <Grid/Hadrons/Modules/MScalarSUN/Div.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MScalarSUN;
template class Grid::Hadrons::MScalarSUN::TDiv<ScalarNxNAdjImplR<2>>;
template class Grid::Hadrons::MScalarSUN::TDiv<ScalarNxNAdjImplR<3>>;
template class Grid::Hadrons::MScalarSUN::TDiv<ScalarNxNAdjImplR<4>>;
template class Grid::Hadrons::MScalarSUN::TDiv<ScalarNxNAdjImplR<5>>;
template class Grid::Hadrons::MScalarSUN::TDiv<ScalarNxNAdjImplR<6>>;

View File

@ -31,43 +31,42 @@ See the full license in the file "LICENSE" in the top level distribution directo
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* Divergence of a vector field *
* Div *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MScalarSUN)
class DivPar: Serializable
{
public:
GRID_SERIALIZABLE_ENUM(DiffType, undef, forward, 1, backward, 2, central, 3);
GRID_SERIALIZABLE_CLASS_MEMBERS(DivPar,
std::vector<std::string>, op,
DiffType, type,
std::string, output);
};
class DivResult: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(DivResult,
DiffType, type,
Complex, value);
};
template <typename SImpl>
class TDiv: public Module<DivPar>
{
public:
typedef typename SImpl::Field Field;
typedef typename SImpl::ComplexField ComplexField;
class Result: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
DivPar::DiffType, type,
Complex, value);
};
public:
// constructor
TDiv(const std::string name);
// destructor
virtual ~TDiv(void) {};
virtual ~TDiv(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
@ -77,14 +76,14 @@ public:
virtual void execute(void);
};
MODULE_REGISTER_TMP(DivSU2, TDiv<ScalarNxNAdjImplR<2>>, MScalarSUN);
MODULE_REGISTER_TMP(DivSU3, TDiv<ScalarNxNAdjImplR<3>>, MScalarSUN);
MODULE_REGISTER_TMP(DivSU4, TDiv<ScalarNxNAdjImplR<4>>, MScalarSUN);
MODULE_REGISTER_TMP(DivSU5, TDiv<ScalarNxNAdjImplR<5>>, MScalarSUN);
MODULE_REGISTER_TMP(DivSU6, TDiv<ScalarNxNAdjImplR<6>>, MScalarSUN);
MODULE_REGISTER_NS(DivSU2, TDiv<ScalarNxNAdjImplR<2>>, MScalarSUN);
MODULE_REGISTER_NS(DivSU3, TDiv<ScalarNxNAdjImplR<3>>, MScalarSUN);
MODULE_REGISTER_NS(DivSU4, TDiv<ScalarNxNAdjImplR<4>>, MScalarSUN);
MODULE_REGISTER_NS(DivSU5, TDiv<ScalarNxNAdjImplR<5>>, MScalarSUN);
MODULE_REGISTER_NS(DivSU6, TDiv<ScalarNxNAdjImplR<6>>, MScalarSUN);
/******************************************************************************
* TDiv implementation *
* TDiv implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename SImpl>
@ -113,7 +112,7 @@ void TDiv<SImpl>::setup(void)
{
if (par().op.size() != env().getNd())
{
HADRONS_ERROR(Size, "the number of components differs from number of dimensions");
HADRON_ERROR(Size, "the number of components differs from number of dimensions");
}
envCreateLat(ComplexField, getName());
}
@ -127,7 +126,7 @@ void TDiv<SImpl>::execute(void)
LOG(Message) << "Computing the " << par().type << " divergence of [";
for (unsigned int mu = 0; mu < nd; ++mu)
{
std::cout << "'" << par().op[mu] << ((mu == nd - 1) ? "']" : "', ");
std::cout << par().op[mu] << ((mu == nd - 1) ? "]" : ", ");
}
std::cout << std::endl;
@ -136,15 +135,27 @@ void TDiv<SImpl>::execute(void)
for (unsigned int mu = 0; mu < nd; ++mu)
{
auto &op = envGet(ComplexField, par().op[mu]);
dmuAcc(div, op, mu, par().type);
switch(par().type)
{
case DivPar::DiffType::backward:
div += op - Cshift(op, mu, -1);
break;
case DivPar::DiffType::forward:
div += Cshift(op, mu, 1) - op;
break;
case DivPar::DiffType::central:
div += 0.5*(Cshift(op, mu, 1) - Cshift(op, mu, -1));
break;
}
}
if (!par().output.empty())
{
DivResult r;
Result r;
ResultWriter writer(RESULT_FILE_NAME(par().output));
r.type = par().type;
r.value = TensorRemove(sum(div));
saveResult(par().output, "div", r);
write(writer, "div", r);
}
}

View File

@ -1,39 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MScalarSUN/EMT.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 <Grid/Hadrons/Modules/MScalarSUN/EMT.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MScalarSUN;
template class Grid::Hadrons::MScalarSUN::TEMT<ScalarNxNAdjImplR<2>>;
template class Grid::Hadrons::MScalarSUN::TEMT<ScalarNxNAdjImplR<3>>;
template class Grid::Hadrons::MScalarSUN::TEMT<ScalarNxNAdjImplR<4>>;
template class Grid::Hadrons::MScalarSUN::TEMT<ScalarNxNAdjImplR<5>>;
template class Grid::Hadrons::MScalarSUN::TEMT<ScalarNxNAdjImplR<6>>;

View File

@ -1,217 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MScalarSUN/EMT.hpp
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_MScalarSUN_EMT_hpp_
#define Hadrons_MScalarSUN_EMT_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* Energy-momentum tensor *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MScalarSUN)
class EMTPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(EMTPar,
std::string, kinetic,
std::string, phiPow,
std::string, improvement,
double , m2,
double , lambda,
double , g,
double , xi,
std::string, output);
};
class EMTResult: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(EMTResult,
std::vector<std::vector<Complex>>, value,
double, m2,
double, lambda,
double, g,
double, xi);
};
template <typename SImpl>
class TEMT: public Module<EMTPar>
{
public:
typedef typename SImpl::Field Field;
typedef typename SImpl::ComplexField ComplexField;
public:
// constructor
TEMT(const std::string name);
// destructor
virtual ~TEMT(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(EMTSU2, TEMT<ScalarNxNAdjImplR<2>>, MScalarSUN);
MODULE_REGISTER_TMP(EMTSU3, TEMT<ScalarNxNAdjImplR<3>>, MScalarSUN);
MODULE_REGISTER_TMP(EMTSU4, TEMT<ScalarNxNAdjImplR<4>>, MScalarSUN);
MODULE_REGISTER_TMP(EMTSU5, TEMT<ScalarNxNAdjImplR<5>>, MScalarSUN);
MODULE_REGISTER_TMP(EMTSU6, TEMT<ScalarNxNAdjImplR<6>>, MScalarSUN);
/******************************************************************************
* TEMT implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename SImpl>
TEMT<SImpl>::TEMT(const std::string name)
: Module<EMTPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename SImpl>
std::vector<std::string> TEMT<SImpl>::getInput(void)
{
std::vector<std::string> in;
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
for (unsigned int nu = mu; nu < env().getNd(); ++nu)
{
in.push_back(varName(par().kinetic, mu, nu));
if (!par().improvement.empty())
{
in.push_back(varName(par().improvement, mu, nu));
}
}
in.push_back(varName(par().kinetic, "sum"));
in.push_back(varName(par().phiPow, 2));
in.push_back(varName(par().phiPow, 4));
return in;
}
template <typename SImpl>
std::vector<std::string> TEMT<SImpl>::getOutput(void)
{
std::vector<std::string> out;
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
for (unsigned int nu = mu; nu < env().getNd(); ++nu)
{
out.push_back(varName(getName(), mu, nu));
}
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename SImpl>
void TEMT<SImpl>::setup(void)
{
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
for (unsigned int nu = mu; nu < env().getNd(); ++nu)
{
envCreateLat(ComplexField, varName(getName(), mu, nu));
}
}
// execution ///////////////////////////////////////////////////////////////////
template <typename SImpl>
void TEMT<SImpl>::execute(void)
{
LOG(Message) << "Computing energy-momentum tensor" << std::endl;
LOG(Message) << " kinetic terms: '" << par().kinetic << "'" << std::endl;
LOG(Message) << " tr(phi^n): '" << par().phiPow << "'" << std::endl;
if (!par().improvement.empty())
{
LOG(Message) << " improvement: '" << par().improvement << "'" << std::endl;
}
LOG(Message) << " m^2= " << par().m2 << std::endl;
LOG(Message) << " lambda= " << par().lambda << std::endl;
LOG(Message) << " g= " << par().g << std::endl;
if (!par().improvement.empty())
{
LOG(Message) << " xi= " << par().xi << std::endl;
}
const unsigned int N = SImpl::Group::Dimension, nd = env().getNd();
auto &trphi2 = envGet(ComplexField, varName(par().phiPow, 2));
auto &trphi4 = envGet(ComplexField, varName(par().phiPow, 4));
auto &sumkin = envGet(ComplexField, varName(par().kinetic, "sum"));
EMTResult result;
if (!par().output.empty())
{
result.m2 = par().m2;
result.g = par().g;
result.lambda = par().lambda;
result.xi = par().xi;
result.value.resize(nd, std::vector<Complex>(nd));
}
for (unsigned int mu = 0; mu < nd; ++mu)
for (unsigned int nu = mu; nu < nd; ++nu)
{
auto &out = envGet(ComplexField, varName(getName(), mu, nu));
auto &trkin = envGet(ComplexField, varName(par().kinetic, mu, nu));
out = 2.*trkin;
if (!par().improvement.empty())
{
auto &imp = envGet(ComplexField, varName(par().improvement, mu, nu));
out += par().xi*imp;
}
if (mu == nu)
{
out -= sumkin + par().m2*trphi2 + par().lambda*trphi4;
}
out *= N/par().g;
if (!par().output.empty())
{
result.value[mu][nu] = TensorRemove(sum(out));
result.value[mu][nu] = result.value[nu][mu];
}
}
if (!par().output.empty())
{
saveResult(par().output, "emt", result);
}
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MScalarSUN_EMT_hpp_

View File

@ -1,38 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MScalarSUN/Grad.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 <Grid/Hadrons/Modules/MScalarSUN/Grad.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MScalarSUN;
template class Grid::Hadrons::MScalarSUN::TGrad<ScalarNxNAdjImplR<2>>;
template class Grid::Hadrons::MScalarSUN::TGrad<ScalarNxNAdjImplR<3>>;
template class Grid::Hadrons::MScalarSUN::TGrad<ScalarNxNAdjImplR<4>>;
template class Grid::Hadrons::MScalarSUN::TGrad<ScalarNxNAdjImplR<5>>;
template class Grid::Hadrons::MScalarSUN::TGrad<ScalarNxNAdjImplR<6>>;

View File

@ -1,166 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MScalarSUN/Grad.hpp
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_MScalarSUN_Grad_hpp_
#define Hadrons_MScalarSUN_Grad_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* Gradient of a complex field *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MScalarSUN)
class GradPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(GradPar,
std::string, op,
DiffType, type,
std::string, output);
};
class GradResult: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(GradResult,
DiffType, type,
std::vector<Complex>, value);
};
template <typename SImpl>
class TGrad: public Module<GradPar>
{
public:
typedef typename SImpl::Field Field;
typedef typename SImpl::ComplexField ComplexField;
public:
// constructor
TGrad(const std::string name);
// destructor
virtual ~TGrad(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(GradSU2, TGrad<ScalarNxNAdjImplR<2>>, MScalarSUN);
MODULE_REGISTER_TMP(GradSU3, TGrad<ScalarNxNAdjImplR<3>>, MScalarSUN);
MODULE_REGISTER_TMP(GradSU4, TGrad<ScalarNxNAdjImplR<4>>, MScalarSUN);
MODULE_REGISTER_TMP(GradSU5, TGrad<ScalarNxNAdjImplR<5>>, MScalarSUN);
MODULE_REGISTER_TMP(GradSU6, TGrad<ScalarNxNAdjImplR<6>>, MScalarSUN);
/******************************************************************************
* TGrad implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename SImpl>
TGrad<SImpl>::TGrad(const std::string name)
: Module<GradPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename SImpl>
std::vector<std::string> TGrad<SImpl>::getInput(void)
{
std::vector<std::string> in = {par().op};
return in;
}
template <typename SImpl>
std::vector<std::string> TGrad<SImpl>::getOutput(void)
{
std::vector<std::string> out;
const auto nd = env().getNd();
for (unsigned int mu = 0; mu < nd; ++mu)
{
out.push_back(varName(getName(), mu));
}
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename SImpl>
void TGrad<SImpl>::setup(void)
{
const auto nd = env().getNd();
for (unsigned int mu = 0; mu < nd; ++mu)
{
envCreateLat(ComplexField, varName(getName(), mu));
}
}
// execution ///////////////////////////////////////////////////////////////////
template <typename SImpl>
void TGrad<SImpl>::execute(void)
{
LOG(Message) << "Computing the " << par().type << " gradient of '"
<< par().op << "'" << std::endl;
const unsigned int nd = env().getNd();
GradResult result;
auto &op = envGet(ComplexField, par().op);
if (!par().output.empty())
{
result.type = par().type;
result.value.resize(nd);
}
for (unsigned int mu = 0; mu < nd; ++mu)
{
auto &der = envGet(ComplexField, varName(getName(), mu));
dmu(der, op, mu, par().type);
if (!par().output.empty())
{
result.value[mu] = TensorRemove(sum(der));
}
}
if (!par().output.empty())
{
saveResult(par().output, "grad", result);
}
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MScalarSUN_Grad_hpp_

View File

@ -1,39 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MScalarSUN/ShiftProbe.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 <Grid/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MScalarSUN;
template class Grid::Hadrons::MScalarSUN::TShiftProbe<ScalarNxNAdjImplR<2>>;
template class Grid::Hadrons::MScalarSUN::TShiftProbe<ScalarNxNAdjImplR<3>>;
template class Grid::Hadrons::MScalarSUN::TShiftProbe<ScalarNxNAdjImplR<4>>;
template class Grid::Hadrons::MScalarSUN::TShiftProbe<ScalarNxNAdjImplR<5>>;
template class Grid::Hadrons::MScalarSUN::TShiftProbe<ScalarNxNAdjImplR<6>>;

View File

@ -1,177 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_MScalarSUN_ShiftProbe_hpp_
#define Hadrons_MScalarSUN_ShiftProbe_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* Ward identity phi^n probe with fields at different positions *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MScalarSUN)
typedef std::pair<int, int> ShiftPair;
class ShiftProbePar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(ShiftProbePar,
std::string, field,
std::string, shifts,
std::string, output);
};
class ShiftProbeResult: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(ShiftProbeResult,
std::string, shifts,
Complex, value);
};
template <typename SImpl>
class TShiftProbe: public Module<ShiftProbePar>
{
public:
typedef typename SImpl::Field Field;
typedef typename SImpl::ComplexField ComplexField;
public:
// constructor
TShiftProbe(const std::string name);
// destructor
virtual ~TShiftProbe(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(ShiftProbeSU2, TShiftProbe<ScalarNxNAdjImplR<2>>, MScalarSUN);
MODULE_REGISTER_TMP(ShiftProbeSU3, TShiftProbe<ScalarNxNAdjImplR<3>>, MScalarSUN);
MODULE_REGISTER_TMP(ShiftProbeSU4, TShiftProbe<ScalarNxNAdjImplR<4>>, MScalarSUN);
MODULE_REGISTER_TMP(ShiftProbeSU5, TShiftProbe<ScalarNxNAdjImplR<5>>, MScalarSUN);
MODULE_REGISTER_TMP(ShiftProbeSU6, TShiftProbe<ScalarNxNAdjImplR<6>>, MScalarSUN);
/******************************************************************************
* TShiftProbe implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename SImpl>
TShiftProbe<SImpl>::TShiftProbe(const std::string name)
: Module<ShiftProbePar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename SImpl>
std::vector<std::string> TShiftProbe<SImpl>::getInput(void)
{
std::vector<std::string> in = {par().field};
return in;
}
template <typename SImpl>
std::vector<std::string> TShiftProbe<SImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename SImpl>
void TShiftProbe<SImpl>::setup(void)
{
envTmpLat(Field, "acc");
envCreateLat(ComplexField, getName());
}
// execution ///////////////////////////////////////////////////////////////////
template <typename SImpl>
void TShiftProbe<SImpl>::execute(void)
{
LOG(Message) << "Creating shift probe for shifts " << par().shifts
<< std::endl;
std::vector<ShiftPair> shift;
double sign;
auto &phi = envGet(Field, par().field);
auto &probe = envGet(ComplexField, getName());
shift = strToVec<ShiftPair>(par().shifts);
if (shift.size() % 2 != 0)
{
HADRONS_ERROR(Size, "the number of shifts is odd");
}
sign = (shift.size() % 4 == 0) ? 1. : -1.;
for (auto &s: shift)
{
if (s.first >= env().getNd())
{
HADRONS_ERROR(Size, "dimension to large for shift <"
+ std::to_string(s.first) + " "
+ std::to_string(s.second) + ">" );
}
}
envGetTmp(Field, acc);
acc = 1.;
for (unsigned int i = 0; i < shift.size(); ++i)
{
if (shift[i].second == 0)
{
acc *= phi;
}
else
{
acc *= Cshift(phi, shift[i].first, shift[i].second);
}
}
probe = sign*trace(acc);
if (!par().output.empty())
{
ShiftProbeResult r;
r.shifts = par().shifts;
r.value = TensorRemove(sum(probe));
saveResult(par().output, "probe", r);
}
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MScalarSUN_ShiftProbe_hpp_

View File

@ -1,11 +0,0 @@
#include <Grid/Hadrons/Modules/MScalarSUN/StochFreeField.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MScalarSUN;
template class Grid::Hadrons::MScalarSUN::TStochFreeField<ScalarNxNAdjImplR<2>>;
template class Grid::Hadrons::MScalarSUN::TStochFreeField<ScalarNxNAdjImplR<3>>;
template class Grid::Hadrons::MScalarSUN::TStochFreeField<ScalarNxNAdjImplR<4>>;
template class Grid::Hadrons::MScalarSUN::TStochFreeField<ScalarNxNAdjImplR<5>>;
template class Grid::Hadrons::MScalarSUN::TStochFreeField<ScalarNxNAdjImplR<6>>;

View File

@ -1,151 +0,0 @@
#ifndef Hadrons_MScalarSUN_StochFreeField_hpp_
#define Hadrons_MScalarSUN_StochFreeField_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* stochastic free SU(N) scalar field *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MScalarSUN)
class StochFreeFieldPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(StochFreeFieldPar,
double, m2,
double, g);
};
template <typename SImpl>
class TStochFreeField: public Module<StochFreeFieldPar>
{
public:
typedef typename SImpl::Field Field;
typedef typename SImpl::ComplexField ComplexField;
typedef typename SImpl::Group Group;
typedef typename SImpl::SiteField::scalar_object Site;
public:
// constructor
TStochFreeField(const std::string name);
// destructor
virtual ~TStochFreeField(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);
private:
bool create_weight;
};
MODULE_REGISTER_TMP(StochFreeFieldSU2, TStochFreeField<ScalarNxNAdjImplR<2>>, MScalarSUN);
MODULE_REGISTER_TMP(StochFreeFieldSU3, TStochFreeField<ScalarNxNAdjImplR<3>>, MScalarSUN);
MODULE_REGISTER_TMP(StochFreeFieldSU4, TStochFreeField<ScalarNxNAdjImplR<4>>, MScalarSUN);
MODULE_REGISTER_TMP(StochFreeFieldSU5, TStochFreeField<ScalarNxNAdjImplR<5>>, MScalarSUN);
MODULE_REGISTER_TMP(StochFreeFieldSU6, TStochFreeField<ScalarNxNAdjImplR<6>>, MScalarSUN);
/******************************************************************************
* TStochFreeField implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename SImpl>
TStochFreeField<SImpl>::TStochFreeField(const std::string name)
: Module<StochFreeFieldPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename SImpl>
std::vector<std::string> TStochFreeField<SImpl>::getInput(void)
{
std::vector<std::string> in;
return in;
}
template <typename SImpl>
std::vector<std::string> TStochFreeField<SImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename SImpl>
void TStochFreeField<SImpl>::setup(void)
{
create_weight = false;
if (!env().hasCreatedObject("_" + getName() + "_weight"))
{
envCacheLat(ComplexField, "_" + getName() + "_weight");
create_weight = true;
}
envTmpLat(Field, "phift");
envTmpLat(ComplexField, "ca");
envCreateLat(Field, getName());
}
// execution ///////////////////////////////////////////////////////////////////
template <typename SImpl>
void TStochFreeField<SImpl>::execute(void)
{
LOG(Message) << "Generating stochastic scalar field" << std::endl;
const unsigned int N = Group::Dimension;
const unsigned int Nadj = Group::AdjointDimension;
auto &phi = envGet(Field, getName());
auto &w = envGet(ComplexField, "_" + getName() + "_weight");
auto &rng = *env().get4dRng();
double trphi2;
FFT fft(env().getGrid());
Integer vol;
vol = 1;
for(int d = 0; d < env().getNd(); d++)
{
vol = vol*env().getDim(d);
}
if (create_weight)
{
LOG(Message) << "Caching momentum-space scalar action" << std::endl;
SImpl::MomentumSpacePropagator(w, sqrt(par().m2));
w *= par().g/N;
w = sqrt(vol)*sqrt(w);
}
LOG(Message) << "Generating random momentum-space field" << std::endl;
envGetTmp(Field, phift);
envGetTmp(ComplexField, ca);
phift = zero;
for (int a = 0; a < Nadj; ++a)
{
Site ta;
gaussian(rng, ca);
Group::generator(a, ta);
phift += ca*ta;
}
phift *= w;
LOG(Message) << "Field Fourier transform" << std::endl;
fft.FFT_all_dim(phi, phift, FFT::backward);
phi = 0.5*(phi - adj(phi));
trphi2 = -TensorRemove(sum(trace(phi*phi))).real()/vol;
LOG(Message) << "tr(phi^2)= " << trphi2 << std::endl;
// ComplexField phi2(env().getGrid());
// phi2=trace(phi*phi);
// std::cout << phi2 << std::endl;
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MScalarSUN_StochFreeField_hpp_

View File

@ -1,11 +0,0 @@
#include <Grid/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MScalarSUN;
template class Grid::Hadrons::MScalarSUN::TTimeMomProbe<ScalarNxNAdjImplR<2>>;
template class Grid::Hadrons::MScalarSUN::TTimeMomProbe<ScalarNxNAdjImplR<3>>;
template class Grid::Hadrons::MScalarSUN::TTimeMomProbe<ScalarNxNAdjImplR<4>>;
template class Grid::Hadrons::MScalarSUN::TTimeMomProbe<ScalarNxNAdjImplR<5>>;
template class Grid::Hadrons::MScalarSUN::TTimeMomProbe<ScalarNxNAdjImplR<6>>;

Some files were not shown because too many files have changed in this diff Show More