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

Hadrons: copyright update

Browse Source
This commit is contained in:
Antonin Portelli 2018-08-14 12:15:24 +01:00
parent f2d3e41cf2
commit 89d2fac92e
20 changed files with 401 additions and 250 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
extras/Hadrons/A2AVectors.hpp
View File

@ -219,4 +219,4 @@ void A2AVectorsSchurDiagTwo<FImpl>::makeHighModeW5D(FermionField &wout_4d,
END_HADRONS_NAMESPACE
#endif // A2A_Vectors_hpp_
#endif // A2A_Vectors_hpp_

2
extras/Hadrons/DilutedNoise.hpp
View File

@ -188,4 +188,4 @@ void TimeDilutedSpinColorDiagonalNoise<FImpl>::generateNoise(GridParallelRNG &rn
END_HADRONS_NAMESPACE
#endif // Hadrons_DilutedNoise_hpp_
#endif // Hadrons_DilutedNoise_hpp_

2
extras/Hadrons/Global.cc
View File

@ -192,7 +192,7 @@ void Hadrons::printTimeProfile(const std::map<std::string, GridTime> &timing,
}
for (auto &rt: rtiming)
{
LOG(Message) << std::right << std::setw(width) << rt.second << ": "
LOG(Message) << std::setw(width) << rt.second << ": "
<< rt.first << " us (" << std::fixed
<< std::setprecision(1)
<< static_cast<double>(rt.first)/dtotal*100 << "%)"

22
extras/Hadrons/Module.cc
View File

@ -78,12 +78,21 @@ void ModuleBase::startTimer(const std::string &name)
GridTime ModuleBase::getTimer(const std::string &name)
{
GridTime t;
if (!name.empty())
{
stopTimer(name);
t = timer_.at(name).Elapsed();
startTimer(name);
try
{
bool running = timer_.at(name).isRunning();
if (running) stopTimer(name);
t = timer_.at(name).Elapsed();
if (running) startTimer(name);
}
catch (std::out_of_range &)
{
t = GridTime::zero();
}
}
else
{
@ -93,6 +102,11 @@ GridTime ModuleBase::getTimer(const std::string &name)
return t;
}
double ModuleBase::getDTimer(const std::string &name)
{
return static_cast<double>(getTimer(name).count());
}
void ModuleBase::startCurrentTimer(const std::string &name)
{
if (!name.empty())

3
extras/Hadrons/Module.hpp
View File

@ -129,8 +129,6 @@ if (env().getGrid()->IsBoss() and !ioStem.empty())\
}\
}
#define MODULE_TIMER(name) this->startCurrentTimer(name)
/******************************************************************************
* Module class *
******************************************************************************/
@ -166,6 +164,7 @@ public:
// timers
void startTimer(const std::string &name);
GridTime getTimer(const std::string &name);
double getDTimer(const std::string &name);
void startCurrentTimer(const std::string &name);
void stopTimer(const std::string &name);
void stopCurrentTimer(void);

35
extras/Hadrons/Modules.hpp
View File

@ -1,5 +1,40 @@
/*************************************************************************************
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: Guido Cossu <guido.cossu@ed.ac.uk>
Author: James Harrison <jch1g10@soton.ac.uk>
Author: Lanny91 <andrew.lawson@gmail.com>
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: Vera Guelpers <vmg1n14@soton.ac.uk>
Author: fionnoh <fionnoh@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
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/A2AMesonField.hpp>
#include <Grid/Hadrons/Modules/MContraction/A2AMesonFieldKernels.hpp>
#include <Grid/Hadrons/Modules/MContraction/Meson.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>

1
extras/Hadrons/Modules/MContraction/A2AMesonField.cc
View File

@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MContraction/A2AMesonField.cc
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
Author: paboyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify

282
extras/Hadrons/Modules/MContraction/A2AMesonField.hpp
View File

@ -34,8 +34,8 @@ See the full license in the file "LICENSE" in the top level distribution directo
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
#include <Grid/Hadrons/A2AVectors.hpp>
#include <Grid/Eigen/unsupported/CXX11/Tensor>
#include <Grid/Hadrons/Modules/MSolver/A2AVectors.hpp>
#include <Grid/Hadrons/Modules/MContraction/A2AMesonFieldKernels.hpp>
BEGIN_HADRONS_NAMESPACE
@ -44,9 +44,6 @@ BEGIN_HADRONS_NAMESPACE
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MContraction)
typedef std::pair<Gamma::Algebra, Gamma::Algebra> GammaPair;
class A2AMesonFieldPar: Serializable
{
public:
@ -66,7 +63,6 @@ class TA2AMesonField : public Module<A2AMesonFieldPar>
public:
FERM_TYPE_ALIASES(FImpl,);
SOLVER_TYPE_ALIASES(FImpl,);
typedef Eigen::TensorMap<Eigen::Tensor<Complex, 5, Eigen::RowMajor>> MesonField;
public:
// constructor
TA2AMesonField(const std::string name);
@ -144,24 +140,24 @@ void TA2AMesonField<FImpl>::setup(void)
mom_.clear();
if (par().gammas == "all")
{
gamma_ = {
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
};
gamma_ = {
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
};
}
else
{
@ -227,7 +223,7 @@ void TA2AMesonField<FImpl>::execute(void)
if (!hasPhase_)
{
MODULE_TIMER("Momentum phases");
startTimer("Momentum phases");
for (unsigned int j = 0; j < nmom; ++j)
{
Complex i(0.0,1.0);
@ -243,9 +239,9 @@ void TA2AMesonField<FImpl>::execute(void)
ph[j] = exp((Real)(2*M_PI)*i*ph[j]);
}
hasPhase_ = true;
stopTimer("Momentum phases");
}
//////////////////////////////////////////////////////////////////////////
// i,j is first loop over SchurBlock factors reusing 5D matrices
// ii,jj is second loop over cacheBlock factors for high perf contractoin
@ -256,12 +252,8 @@ void TA2AMesonField<FImpl>::execute(void)
double flops = 0.0;
double bytes = 0.0;
double vol = env().getVolume();
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;
envGetTmp(Vector<Complex>, mfBuf);
envGetTmp(Vector<Complex>, mfCache);
@ -277,9 +269,6 @@ void TA2AMesonField<FImpl>::execute(void)
int N_ii = MIN(N_i-i,block);
int N_jj = MIN(N_j-j,block);
t_schur-=usecond();
t_schur+=usecond();
LOG(Message) << "Meson field block "
<< j/block + NBlock_j*i/block + 1
<< "/" << NBlock_i*NBlock_j << " [" << i <<" .. "
@ -296,17 +285,17 @@ void TA2AMesonField<FImpl>::execute(void)
int N_jjj = MIN(N_jj-jj,cacheBlock);
MesonField mfCacheBlock(mfCache.data(),nmom,ngamma,nt,N_iii,N_jjj);
t_contr-=usecond();
makeBlock(mfCacheBlock, &w[i+ii], &v[j+jj], gamma_, ph,
env().getNd() - 1, t_int_0, t_int_1, t_int_2, t_int_3);
t_contr+=usecond();
startTimer("contraction: total");
makeMesonFieldBlock(mfCacheBlock, &w[i+ii], &v[j+jj], gamma_, ph,
env().getNd() - 1, this);
stopTimer("contraction: total");
// flops for general N_c & N_s
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;
MODULE_TIMER("Cache copy");
startTimer("cache copy");
for(int iii=0;iii< N_iii;iii++)
for(int jjj=0;jjj< N_jjj;jjj++)
for(int m =0;m< nmom;m++)
@ -315,223 +304,44 @@ void TA2AMesonField<FImpl>::execute(void)
{
mfBlock(m,g,t,ii+iii,jj+jjj) = mfCacheBlock(m,g,t,iii,jjj);
}
stopTimer("cache copy");
}
// IO
if (!par().output.empty())
{
double blockSize, ioTime;
MODULE_TIMER("IO");
LOG(Message) << "Writing block to disk" << std::endl;
ioTime = -getDTimer("IO: write block");
startTimer("IO: total");
for(int m = 0; m < nmom; m++)
for(int g = 0; g < ngamma; g++)
{
if ((i == 0) and (j == 0))
{
startTimer("IO: file creation");
initFile(m, g);
stopTimer("IO: file creation");
}
startTimer("IO: write block");
saveBlock(mfBlock, m, g, i, j);
stopTimer("IO: write block");
}
blockSize = static_cast<double>(nmom*ngamma*nt*N_ii*N_jj*sizeof(Complex));
ioTime = static_cast<double>(this->getTimer("IO").count());
LOG(Message) << "HDF5 IO " << blockSize/ioTime*1.0e6/1024/1024
stopTimer("IO: total");
blockSize = static_cast<double>(nmom*ngamma*nt*N_ii*N_jj*sizeof(Complex));
ioTime += getDTimer("IO: write block");
LOG(Message) << "HDF5 IO done " << blockSize/ioTime*1.0e6/1024/1024
<< " MB/s" << std::endl;
}
}
double nodes = env().getGrid()->NodeCount();
double t_kernel = t_int_0 + t_int_1;
double t_kernel = getDTimer("contraction: colour trace & mom.")
+ getDTimer("contraction: local space sum");
LOG(Message) << "Perf " << flops/(t_kernel)/1.0e3/nodes << " Gflop/s/node " << std::endl;
LOG(Message) << "Perf " << bytes/(t_kernel)/1.0e3/nodes << " GB/s/node " << std::endl;
}
//////////////////////////////////////////////////////////////////////////////////
// Cache blocked arithmetic routine
// Could move to Grid ???
//////////////////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TA2AMesonField<FImpl>::makeBlock(MesonField &mat,
const FermionField *lhs_wi,
const FermionField *rhs_vj,
std::vector<Gamma::Algebra> gamma,
const std::vector<LatticeComplex> &mom,
int orthogdim,
double &t0,
double &t1,
double &t2,
double &t3)
{
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;
int Lblock = mat.dimension(3);
int Rblock = mat.dimension(4);
GridBase *grid = lhs_wi[0]._grid;
const int Nd = grid->_ndimension;
const int Nsimd = grid->Nsimd();
int Nt = grid->GlobalDimensions()[orthogdim];
int Ngamma = gamma.size();
int Nmom = mom.size();
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];
t0-=usecond();
MODULE_TIMER("Colour trace * mom.");
// Nested parallelism would be ok
// Wasting cores here. Test case r
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_wi[i]._odata[ss]);
for(int j=0;j<Rblock;j++)
{
SpinMatrix_v vv;
auto right = rhs_vj[j]._odata[ss];
for(int s1=0;s1<Ns;s1++)
for(int s2=0;s2<Ns;s2++)
{
vv()(s1,s2)() = left()(s2)(0) * right()(s1)(0)
+ left()(s2)(1) * right()(s1)(1)
+ left()(s2)(2) * right()(s1)(2);
}
// After getting the sitewise product do the mom phase loop
int base = Nmom*i+Nmom*Lblock*j+Nmom*Lblock*Rblock*r;
for ( int m=0;m<Nmom;m++)
{
int idx = m+base;
auto phase = mom[m]._odata[ss];
mac(&lvSum[idx],&vv,&phase);
}
}
}
}
}
t0+=usecond();
// Sum across simd lanes in the plane, breaking out orthog dir.
MODULE_TIMER("Local space sum");
t1-=usecond();
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];
}
}
}
t1+=usecond();
assert(mat.dimension(0) == Nmom);
assert(mat.dimension(1) == Ngamma);
assert(mat.dimension(2) == Nt);
t2-=usecond();
// ld loop and local only??
MODULE_TIMER("Spin trace");
int pd = grid->_processors[orthogdim];
int pc = grid->_processor_coor[orthogdim];
parallel_for_nest2(int lt=0;lt<ld;lt++)
{
for(int pt=0;pt<pd;pt++)
{
int t = lt + pt*ld;
if (pt == pc)
{
for(int i=0;i<Lblock;i++)
for(int j=0;j<Rblock;j++)
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++)
{
// this is a bit slow
mat(m,mu,t,i,j) = trace(lsSum[ij_dx]*Gamma(gamma[mu]));
}
}
}
else
{
const scalar_type zz(0.0);
for(int i=0;i<Lblock;i++)
for(int j=0;j<Rblock;j++)
for(int mu=0;mu<Ngamma;mu++)
for(int m=0;m<Nmom;m++)
{
mat(m,mu,t,i,j) =zz;
}
}
}
}
t2+=usecond();
////////////////////////////////////////////////////////////////////
// This global sum is taking as much as 50% of time on 16 nodes
// Vector size is 7 x 16 x 32 x 16 x 16 x sizeof(complex) = 2MB - 60MB depending on volume
// Healthy size that should suffice
////////////////////////////////////////////////////////////////////
t3-=usecond();
MODULE_TIMER("Global sum");
grid->GlobalSumVector(&mat(0,0,0,0,0),Nmom*Ngamma*Nt*Lblock*Rblock);
t3+=usecond();
LOG(Message) << "Perf " << flops/t_kernel/1.0e3/nodes << " Gflop/s/node " << std::endl;
LOG(Message) << "Perf " << bytes/t_kernel/1.0e3/nodes << " GB/s/node " << std::endl;
}
// IO

226
extras/Hadrons/Modules/MContraction/A2AMesonFieldKernels.hpp Normal file
View File

@ -0,0 +1,226 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MContraction/A2AMesonFieldKernels.hpp
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_MContraction_A2AMesonFieldKernels_hpp_
#define Hadrons_MContraction_A2AMesonFieldKernels_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Eigen/unsupported/CXX11/Tensor>
BEGIN_HADRONS_NAMESPACE
BEGIN_MODULE_NAMESPACE(MContraction)
typedef Eigen::TensorMap<Eigen::Tensor<Complex, 5, Eigen::RowMajor>> MesonField;
////////////////////////////////////////////////////////////////////////////////
// Cache blocked arithmetic routine
// Could move to Grid ???
////////////////////////////////////////////////////////////////////////////////
template <typename Field>
void makeMesonFieldBlock(MesonField &mat,
const Field *lhs_wi,
const Field *rhs_vj,
std::vector<Gamma::Algebra> gamma,
const std::vector<LatticeComplex> &mom,
int orthogdim,
ModuleBase *caller = nullptr)
{
typedef typename Field::vector_object 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;
int Lblock = mat.dimension(3);
int Rblock = mat.dimension(4);
GridBase *grid = lhs_wi[0]._grid;
const int Nd = grid->_ndimension;
const int Nsimd = grid->Nsimd();
int Nt = grid->GlobalDimensions()[orthogdim];
int Ngamma = gamma.size();
int Nmom = mom.size();
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];
if (caller) caller->startTimer("contraction: colour trace & mom.");
// Nested parallelism would be ok
// Wasting cores here. Test case r
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_wi[i]._odata[ss]);
for(int j=0;j<Rblock;j++)
{
SpinMatrix_v vv;
auto right = rhs_vj[j]._odata[ss];
for(int s1=0;s1<Ns;s1++)
for(int s2=0;s2<Ns;s2++)
{
vv()(s1,s2)() = left()(s2)(0) * right()(s1)(0)
+ left()(s2)(1) * right()(s1)(1)
+ left()(s2)(2) * right()(s1)(2);
}
// After getting the sitewise product do the mom phase loop
int base = Nmom*i+Nmom*Lblock*j+Nmom*Lblock*Rblock*r;
for ( int m=0;m<Nmom;m++)
{
int idx = m+base;
auto phase = mom[m]._odata[ss];
mac(&lvSum[idx],&vv,&phase);
}
}
}
}
}
if (caller) caller->stopTimer("contraction: colour trace & mom.");
// Sum across simd lanes in the plane, breaking out orthog dir.
if (caller) caller->startTimer("contraction: local space sum");
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];
}
}
}
if (caller) caller->stopTimer("contraction: local space sum");
// ld loop and local only??
if (caller) caller->startTimer("contraction: spin trace");
int pd = grid->_processors[orthogdim];
int pc = grid->_processor_coor[orthogdim];
parallel_for_nest2(int lt=0;lt<ld;lt++)
{
for(int pt=0;pt<pd;pt++)
{
int t = lt + pt*ld;
if (pt == pc)
{
for(int i=0;i<Lblock;i++)
for(int j=0;j<Rblock;j++)
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++)
{
// this is a bit slow
mat(m,mu,t,i,j) = trace(lsSum[ij_dx]*Gamma(gamma[mu]));
}
}
}
else
{
const scalar_type zz(0.0);
for(int i=0;i<Lblock;i++)
for(int j=0;j<Rblock;j++)
for(int mu=0;mu<Ngamma;mu++)
for(int m=0;m<Nmom;m++)
{
mat(m,mu,t,i,j) =zz;
}
}
}
}
if (caller) caller->stopTimer("contraction: spin trace");
////////////////////////////////////////////////////////////////////
// This global sum is taking as much as 50% of time on 16 nodes
// Vector size is 7 x 16 x 32 x 16 x 16 x sizeof(complex) = 2MB - 60MB depending on volume
// Healthy size that should suffice
////////////////////////////////////////////////////////////////////
if (caller) caller->startTimer("contraction: global sum");
grid->GlobalSumVector(&mat(0,0,0,0,0),Nmom*Ngamma*Nt*Lblock*Rblock);
if (caller) caller->stopTimer("contraction: global sum");
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif //Hadrons_MContraction_A2AMesonField_hpp_

1
extras/Hadrons/Modules/MFermion/FreeProp.cc
View File

@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MFermion/FreeProp.cc
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
Author: Vera Guelpers <vmg1n14@soton.ac.uk>
This program is free software; you can redistribute it and/or modify

1
extras/Hadrons/Modules/MGauge/UnitEm.cc
View File

@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MGauge/UnitEm.cc
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
Author: James Harrison <jch1g10@soton.ac.uk>
This program is free software; you can redistribute it and/or modify

1
extras/Hadrons/Modules/MGauge/UnitEm.hpp
View File

@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MGauge/UnitEm.hpp
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
Author: James Harrison <jch1g10@soton.ac.uk>
This program is free software; you can redistribute it and/or modify

1
extras/Hadrons/Modules/MScalar/ScalarVP.cc
View File

@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MScalar/ScalarVP.cc
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
Author: James Harrison <jch1g10@soton.ac.uk>
This program is free software; you can redistribute it and/or modify

1
extras/Hadrons/Modules/MScalar/ScalarVP.hpp
View File

@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MScalar/ScalarVP.hpp
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
Author: James Harrison <jch1g10@soton.ac.uk>
This program is free software; you can redistribute it and/or modify

1
extras/Hadrons/Modules/MScalar/VPCounterTerms.cc
View File

@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MScalar/VPCounterTerms.cc
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
Author: James Harrison <jch1g10@soton.ac.uk>
This program is free software; you can redistribute it and/or modify

1
extras/Hadrons/Modules/MScalar/VPCounterTerms.hpp
View File

@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MScalar/VPCounterTerms.hpp
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
Author: James Harrison <jch1g10@soton.ac.uk>
This program is free software; you can redistribute it and/or modify

14
extras/Hadrons/Modules/MSolver/A2AVectors.hpp
View File

@ -178,7 +178,7 @@ void TA2AVectors<FImpl, Pack>::execute(void)
{
auto &epack = envGet(Pack, par().eigenPack);
MODULE_TIMER("V low mode");
startTimer("V low mode");
LOG(Message) << "V vector i = " << il << " (low mode)" << std::endl;
if (Ls == 1)
{
@ -189,7 +189,8 @@ void TA2AVectors<FImpl, Pack>::execute(void)
envGetTmp(FermionField, f5);
a2a.makeLowModeV5D(v[il], f5, epack.evec[il], epack.eval[il]);
}
MODULE_TIMER("W low mode");
stopTimer("V low mode");
startTimer("W low mode");
LOG(Message) << "W vector i = " << il << " (low mode)" << std::endl;
if (Ls == 1)
{
@ -199,13 +200,14 @@ void TA2AVectors<FImpl, Pack>::execute(void)
{
envGetTmp(FermionField, f5);
a2a.makeLowModeW5D(w[il], f5, epack.evec[il], epack.eval[il]);
}
}
stopTimer("W low mode");
}
// High modes
for (unsigned int ih = 0; ih < noise.size(); ih++)
{
MODULE_TIMER("V high mode");
startTimer("V high mode");
LOG(Message) << "V vector i = " << Nl_ + ih
<< " (" << ((Nl_ > 0) ? "high " : "")
<< "stochastic mode)" << std::endl;
@ -218,7 +220,8 @@ void TA2AVectors<FImpl, Pack>::execute(void)
envGetTmp(FermionField, f5);
a2a.makeHighModeV5D(v[Nl_ + ih], f5, noise[ih]);
}
MODULE_TIMER("W high mode");
stopTimer("V high mode");
startTimer("W high mode");
LOG(Message) << "W vector i = " << Nl_ + ih
<< " (" << ((Nl_ > 0) ? "high " : "")
<< "stochastic mode)" << std::endl;
@ -231,6 +234,7 @@ void TA2AVectors<FImpl, Pack>::execute(void)
envGetTmp(FermionField, f5);
a2a.makeHighModeW5D(w[Nl_ + ih], f5, noise[ih]);
}
stopTimer("W high mode");
}
}

27
extras/Hadrons/Modules/MUtilities/RandomVectors.cc
View File

@ -1,3 +1,30 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MUtilities/RandomVectors.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/MUtilities/RandomVectors.hpp>
using namespace Grid;

27
extras/Hadrons/Modules/MUtilities/RandomVectors.hpp
View File

@ -1,3 +1,30 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MUtilities/RandomVectors.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_MUtilities_RandomVectors_hpp_
#define Hadrons_MUtilities_RandomVectors_hpp_

1
extras/Hadrons/modules.inc
View File

@ -61,6 +61,7 @@ modules_cc =\
modules_hpp =\
Modules/MContraction/Baryon.hpp \
Modules/MContraction/A2AMesonField.hpp \
Modules/MContraction/A2AMesonFieldKernels.hpp \
Modules/MContraction/Meson.hpp \
Modules/MContraction/WeakHamiltonian.hpp \
Modules/MContraction/WeakHamiltonianNonEye.hpp \