portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-09-20 09:15:38 +01:00
Code Releases Activity

Merge pull request #174 from fionnoh/a2a_basics

Browse Source 
A2A basics
This commit is contained in:
Antonin Portelli 2018-08-03 16:32:14 +01:00 committed by GitHub
commit 2cfa0b0e6b
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
23 changed files with 2127 additions and 215 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

600
benchmarks/Benchmark_meson_field.cc
View File

@ -66,7 +66,11 @@ void sliceInnerProductMesonField(std::vector< std::vector<ComplexD> > &mat,
// 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);
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];
@ -87,7 +91,6 @@ void sliceInnerProductMesonField(std::vector< std::vector<ComplexD> > &mat,
for(int j=0;j<Rblock;j++){
int idx = i+Lblock*j+Lblock*Rblock*r;
auto right = rhs[j]._odata[ss];
#if 1
vector_type vv = left()(0)(0) * right()(0)(0)
+ left()(0)(1) * right()(0)(1)
+ left()(0)(2) * right()(0)(2)
@ -100,9 +103,6 @@ void sliceInnerProductMesonField(std::vector< std::vector<ComplexD> > &mat,
+ left()(3)(0) * right()(3)(0)
+ left()(3)(1) * right()(3)(1)
+ left()(3)(2) * right()(3)(2);
#else
vector_type vv = TensorRemove(innerProduct(left,right));
#endif
lvSum[idx]=lvSum[idx]+vv;
}
}
@ -140,11 +140,19 @@ void sliceInnerProductMesonField(std::vector< std::vector<ComplexD> > &mat,
}
std::cout << GridLogMessage << " Entering non parallel loop "<<std::endl;
for(int t=0;t<ld;t++){
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++){
int ij_dx = i+Lblock*j+Lblock*Rblock*t;
mat[i+j*Lblock][t] = lsSum[ij_dx];
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;
@ -152,6 +160,447 @@ void sliceInnerProductMesonField(std::vector< std::vector<ComplexD> > &mat,
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,
@ -159,6 +608,31 @@ template void sliceInnerProductMesonField<SpinColourVector>(std::vector< std::ve
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);
@ -167,7 +641,8 @@ int main (int argc, char ** argv)
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++){
@ -179,28 +654,40 @@ int main (int argc, char ** argv)
pRNG.SeedFixedIntegers(seeds);
const int Nm = 32; // number of all modes (high + low)
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<Nm;i++) {
for(int j=0;j<Nm;j++) {
MesonFields [i+j*Nm].resize(nt);
MesonFieldsRef[i+j*Nm].resize(nt);
}}
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);
@ -214,18 +701,70 @@ int main (int argc, char ** argv)
}
}}
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;
double t2 = usecond();
t0 = usecond();
sliceInnerProductMesonField(MesonFields,w,v,Tp);
double t3 = usecond();
std::cout<<GridLogMessage << "Done "<< flops/(t3-t2) <<" mflops " <<std::endl;
std::cout<<GridLogMessage << "Done "<< byte /(t3-t2) <<" MB/s " <<std::endl;
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++) {
@ -236,6 +775,29 @@ int main (int argc, char ** argv)
}}
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();
}

2
bootstrap.sh
View File

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

4
configure.ac
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 $AM_CXXFLAGS"
AM_CFLAGS="-I${abs_srcdir}/include $AM_CFLAGS"
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_LDFLAGS="-L${cwd}/lib $AM_LDFLAGS"
AC_SUBST([AM_CFLAGS])
AC_SUBST([AM_CXXFLAGS])

146
extras/Hadrons/AllToAllReduction.hpp Normal file
View File

@ -0,0 +1,146 @@
#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_

217
extras/Hadrons/AllToAllVectors.hpp Normal file
View File

@ -0,0 +1,217 @@
#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;
const int Nl, Nh;
const bool return_5d;
std::vector<Field> w_high_5d, v_high_5d, w_high_4d, v_high_4d;
public:
A2AModesSchurDiagTwo(const std::vector<Field> *_evec, const std::vector<RealD> *_eval,
Matrix &_action,
Solver &_solver,
const int _Nl, const int _Nh,
const bool _return_5d)
: evec(_evec), eval(_eval),
action(_action),
solver(_solver),
Nl(_Nl), Nh(_Nh),
return_5d(_return_5d)
{
init_resize(1, Nh);
if (return_5d) init_resize(Nh, Nh);
};
void init_resize(const size_t size_5d, const size_t size_4d)
{
GridBase *grid_5d = action.Grid();
GridBase *grid_4d = action.GaugeGrid();
w_high_5d.resize(size_5d, grid_5d);
v_high_5d.resize(size_5d, grid_5d);
w_high_4d.resize(size_4d, grid_4d);
v_high_4d.resize(size_4d, grid_4d);
}
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_

1
extras/Hadrons/Makefile.am
View File

@ -14,6 +14,7 @@ libHadrons_a_SOURCES = \
libHadrons_adir = $(pkgincludedir)/Hadrons
nobase_libHadrons_a_HEADERS = \
$(modules_hpp) \
AllToAllVectors.hpp \
Application.hpp \
EigenPack.hpp \
Environment.hpp \

96
extras/Hadrons/Modules.hpp
View File

@ -1,57 +1,59 @@
#include <Grid/Hadrons/Modules/MContraction/Baryon.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/WardIdentity.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>
#include <Grid/Hadrons/Modules/MFermion/FreeProp.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/LocalCoherenceLanczos.hpp>
#include <Grid/Hadrons/Modules/MSolver/RBPrecCG.hpp>
#include <Grid/Hadrons/Modules/MGauge/UnitEm.hpp>
#include <Grid/Hadrons/Modules/MGauge/StoutSmearing.hpp>
#include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
#include <Grid/Hadrons/Modules/MGauge/Random.hpp>
#include <Grid/Hadrons/Modules/MGauge/FundtoHirep.hpp>
#include <Grid/Hadrons/Modules/MGauge/StochEm.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/VPCounterTerms.hpp>
#include <Grid/Hadrons/Modules/MScalar/ScalarVP.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/MobiusDWF.hpp>
#include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
#include <Grid/Hadrons/Modules/MAction/WilsonClover.hpp>
#include <Grid/Hadrons/Modules/MAction/ZMobiusDWF.hpp>
#include <Grid/Hadrons/Modules/MAction/ScaledDWF.hpp>
#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/ShiftProbe.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/TimeMomProbe.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/MScalarSUN/Utils.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TransProj.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/Grad.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TrKinetic.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/LoadNersc.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/WeakHamiltonianEye.hpp>
#include <Grid/Hadrons/Modules/MContraction/Baryon.hpp>
#include <Grid/Hadrons/Modules/MContraction/Meson.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.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/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/MSolver/RBPrecCG.hpp>
#include <Grid/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp>
#include <Grid/Hadrons/Modules/MSolver/A2AVectors.hpp>

8
extras/Hadrons/Modules/MContraction/A2AMesonField.cc Normal file
View File

@ -0,0 +1,8 @@
#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>;

480
extras/Hadrons/Modules/MContraction/A2AMesonField.hpp Normal file
View File

@ -0,0 +1,480 @@
#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 <unsupported/Eigen/CXX11/Tensor>
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);
// Arithmetic help. Move to Grid??
virtual void MesonField(Eigen::Tensor<ComplexD,5> &mat,
const LatticeFermion *lhs,
const LatticeFermion *rhs,
std::vector<Gamma::Algebra> gammas,
const std::vector<LatticeComplex > &mom,
int orthogdim,
double &t0,
double &t1,
double &t2,
double &t3);
};
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_);
}
//////////////////////////////////////////////////////////////////////////////////
// Cache blocked arithmetic routine
// Could move to Grid ???
//////////////////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TA2AMesonField<FImpl>::MesonField(Eigen::Tensor<ComplexD,5> &mat,
const LatticeFermion *lhs_wi,
const LatticeFermion *rhs_vj,
std::vector<Gamma::Algebra> gammas,
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 = gammas.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();
// 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.
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??
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(gammas[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();
grid->GlobalSumVector(&mat(0,0,0,0,0),Nmom*Ngamma*Nt*Lblock*Rblock);
t3+=usecond();
}
// 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();
MesonField(mesonFieldBlocked, &w[ii], &v[jj], gammas, phases,Tp,
t_int_0,t_int_1,t_int_2,t_int_3);
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;
LOG(Message) << " Intern0 "<<(t_int_0)/1.0e6<< " seconds " << std::endl;
LOG(Message) << " Intern1 "<<(t_int_1)/1.0e6<< " seconds " << std::endl;
LOG(Message) << " Intern2 "<<(t_int_2)/1.0e6<< " seconds " << std::endl;
LOG(Message) << " Intern3 "<<(t_int_3)/1.0e6<< " seconds " << std::endl;
double t_kernel = t_int_0 + t_int_1;
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> 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_

8
extras/Hadrons/Modules/MSolver/A2AVectors.cc Normal file
View File

@ -0,0 +1,8 @@
#include <Grid/Hadrons/Modules/MSolver/A2AVectors.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MSolver;
template class Grid::Hadrons::MSolver::TA2AVectors<FIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>;
template class Grid::Hadrons::MSolver::TA2AVectors<ZFIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>;

235
extras/Hadrons/Modules/MSolver/A2AVectors.hpp Normal file
View File

@ -0,0 +1,235 @@
#ifndef Hadrons_MSolver_A2AVectors_hpp_
#define Hadrons_MSolver_A2AVectors_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
#include <Grid/Hadrons/Solver.hpp>
#include <Grid/Hadrons/EigenPack.hpp>
#include <Grid/Hadrons/AllToAllVectors.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* A2AVectors *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MSolver)
class A2AVectorsPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(A2AVectorsPar,
bool, return_5d,
int, Nl,
int, N,
std::vector<std::string>, sources,
std::string, action,
std::string, eigenPack,
std::string, solver);
};
template <typename FImpl, int nBasis>
class TA2AVectors : public Module<A2AVectorsPar>
{
public:
FERM_TYPE_ALIASES(FImpl,);
SOLVER_TYPE_ALIASES(FImpl,);
typedef FermionEigenPack<FImpl> EPack;
typedef CoarseFermionEigenPack<FImpl, nBasis> CoarseEPack;
typedef A2AModesSchurDiagTwo<typename FImpl::FermionField, FMat, Solver> A2ABase;
public:
// constructor
TA2AVectors(const std::string name);
// destructor
virtual ~TA2AVectors(void) {};
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getReference(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
private:
unsigned int Ls_;
std::string className_;
};
MODULE_REGISTER_TMP(A2AVectors, ARG(TA2AVectors<FIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>), MSolver);
MODULE_REGISTER_TMP(ZA2AVectors, ARG(TA2AVectors<ZFIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>), MSolver);
/******************************************************************************
* TA2AVectors implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl, int nBasis>
TA2AVectors<FImpl, nBasis>::TA2AVectors(const std::string name)
: Module<A2AVectorsPar>(name)
, className_ (name + "_class")
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl, int nBasis>
std::vector<std::string> TA2AVectors<FImpl, nBasis>::getInput(void)
{
int Nl = par().Nl;
std::string sub_string = "";
if (Nl > 0) sub_string = "_subtract";
std::vector<std::string> in = {par().solver + sub_string};
int n = par().sources.size();
for (unsigned int t = 0; t < n; t += 1)
{
in.push_back(par().sources[t]);
}
return in;
}
template <typename FImpl, int nBasis>
std::vector<std::string> TA2AVectors<FImpl, nBasis>::getReference(void)
{
std::vector<std::string> ref = {par().action};
if (!par().eigenPack.empty())
{
ref.push_back(par().eigenPack);
}
return ref;
}
template <typename FImpl, int nBasis>
std::vector<std::string> TA2AVectors<FImpl, nBasis>::getOutput(void)
{
std::vector<std::string> out = {getName(), className_};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl, int nBasis>
void TA2AVectors<FImpl, nBasis>::setup(void)
{
int N = par().N;
int Nl = par().Nl;
int Nh = N - Nl;
bool return_5d = par().return_5d;
int Ls;
std::string sub_string = "";
if (Nl > 0) sub_string = "_subtract";
auto &solver = envGet(Solver, par().solver + sub_string);
Ls = env().getObjectLs(par().solver + sub_string);
auto &action = envGet(FMat, par().action);
envTmpLat(FermionField, "ferm_src", Ls);
envTmpLat(FermionField, "unphys_ferm", Ls);
envTmpLat(FermionField, "tmp");
std::vector<FermionField> *evec;
const std::vector<RealD> *eval;
if (Nl > 0)
{
// Low modes
auto &epack = envGet(EPack, par().eigenPack);
LOG(Message) << "Creating a2a vectors " << getName() <<
" using eigenpack '" << par().eigenPack << "' ("
<< epack.evec.size() << " modes)" <<
" and " << Nh << " high modes." << std::endl;
evec = &epack.evec;
eval = &epack.eval;
}
else
{
LOG(Message) << "Creating a2a vectors " << getName() <<
" using " << Nh << " high modes only." << std::endl;
}
envCreate(A2ABase, className_, Ls,
evec, eval,
action,
solver,
Nl, Nh,
return_5d);
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl, int nBasis>
void TA2AVectors<FImpl, nBasis>::execute(void)
{
auto &action = envGet(FMat, par().action);
int Nt = env().getDim(Tp);
int Nc = FImpl::Dimension;
int Ls_;
int Nl = par().Nl;
std::string sub_string = "";
if (Nl > 0) sub_string = "_subtract";
Ls_ = env().getObjectLs(par().solver + sub_string);
auto &a2areturn = envGet(A2ABase, className_);
// High modes
auto sources = par().sources;
int Nsrc = par().sources.size();
envGetTmp(FermionField, ferm_src);
envGetTmp(FermionField, unphys_ferm);
envGetTmp(FermionField, tmp);
int N_count = 0;
for (unsigned int s = 0; s < Ns; ++s)
for (unsigned int c = 0; c < Nc; ++c)
for (unsigned int T = 0; T < Nsrc; T++)
{
auto &prop_src = envGet(PropagatorField, sources[T]);
LOG(Message) << "A2A src for s = " << s << " , c = " << c << ", T = " << T << std::endl;
// source conversion for 4D sources
if (!env().isObject5d(sources[T]))
{
if (Ls_ == 1)
{
PropToFerm<FImpl>(ferm_src, prop_src, s, c);
tmp = ferm_src;
}
else
{
PropToFerm<FImpl>(tmp, prop_src, s, c);
action.ImportPhysicalFermionSource(tmp, ferm_src);
action.ImportUnphysicalFermion(tmp, unphys_ferm);
}
}
// source conversion for 5D sources
else
{
if (Ls_ != env().getObjectLs(sources[T]))
{
HADRONS_ERROR(Size, "Ls mismatch between quark action and source");
}
else
{
PropToFerm<FImpl>(ferm_src, prop_src, s, c);
action.ExportPhysicalFermionSolution(ferm_src, tmp);
unphys_ferm = ferm_src;
}
}
LOG(Message) << "a2areturn.high_modes Ncount = " << N_count << std::endl;
a2areturn.high_modes(ferm_src, unphys_ferm, tmp, N_count);
N_count++;
}
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MSolver_A2AVectors_hpp_

25
extras/Hadrons/Modules/MSolver/RBPrecCG.hpp
View File

@ -118,7 +118,7 @@ std::vector<std::string> TRBPrecCG<FImpl, nBasis>::getReference(void)
template <typename FImpl, int nBasis>
std::vector<std::string> TRBPrecCG<FImpl, nBasis>::getOutput(void)
{
std::vector<std::string> out = {getName()};
std::vector<std::string> out = {getName(), getName() + "_subtract"};
return out;
}
@ -158,7 +158,7 @@ void TRBPrecCG<FImpl, nBasis>::setup(void)
guesser.reset(new CoarseGuesser(epack.evec, epack.evecCoarse,
epack.evalCoarse));
}
catch (Exceptions::ObjectDefinition &)
catch (Exceptions::Definition &e)
{
auto &epack = envGet(EPack, par().eigenPack);
@ -168,19 +168,22 @@ void TRBPrecCG<FImpl, nBasis>::setup(void)
guesser.reset(new FineGuesser(epack.evec, epack.eval));
}
}
auto solver = [&mat, guesser, this](FermionField &sol,
const FermionField &source)
{
ConjugateGradient<FermionField> cg(par().residual,
par().maxIteration);
HADRONS_DEFAULT_SCHUR_SOLVE<FermionField> schurSolver(cg);
schurSolver(mat, source, sol, *guesser);
auto makeSolver = [&mat, guesser, this](bool subGuess) {
return [&mat, guesser, subGuess, this](FermionField &sol,
const FermionField &source) {
ConjugateGradient<FermionField> cg(par().residual,
par().maxIteration);
HADRONS_DEFAULT_SCHUR_SOLVE<FermionField> schurSolver(cg);
schurSolver.subtractGuess(subGuess);
schurSolver(mat, source, sol, *guesser);
};
};
auto solver = makeSolver(false);
envCreate(Solver, getName(), Ls, solver, mat);
auto solver_subtract = makeSolver(true);
envCreate(Solver, getName() + "_subtract", Ls, solver_subtract, mat);
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl, int nBasis>
void TRBPrecCG<FImpl, nBasis>::execute(void)

193
extras/Hadrons/modules.inc
View File

@ -1,115 +1,120 @@
modules_cc =\
Modules/MContraction/WeakHamiltonianEye.cc \
Modules/MContraction/Baryon.cc \
Modules/MContraction/Meson.cc \
Modules/MContraction/WeakNeutral4ptDisc.cc \
Modules/MContraction/WeakHamiltonianNonEye.cc \
Modules/MContraction/WardIdentity.cc \
Modules/MContraction/DiscLoop.cc \
Modules/MContraction/Gamma3pt.cc \
Modules/MFermion/FreeProp.cc \
Modules/MFermion/GaugeProp.cc \
Modules/MSource/Point.cc \
Modules/MSource/Wall.cc \
Modules/MSource/SeqConserved.cc \
Modules/MSource/SeqGamma.cc \
Modules/MSource/Z2.cc \
Modules/MSink/Point.cc \
Modules/MSink/Smear.cc \
Modules/MSolver/RBPrecCG.cc \
Modules/MSolver/LocalCoherenceLanczos.cc \
Modules/MGauge/StoutSmearing.cc \
Modules/MGauge/Unit.cc \
Modules/MGauge/UnitEm.cc \
Modules/MGauge/StochEm.cc \
Modules/MGauge/Random.cc \
Modules/MGauge/FundtoHirep.cc \
Modules/MUtilities/TestSeqGamma.cc \
Modules/MUtilities/TestSeqConserved.cc \
Modules/MLoop/NoiseLoop.cc \
Modules/MScalarSUN/ShiftProbe.cc \
Modules/MScalarSUN/Grad.cc \
Modules/MScalarSUN/TwoPointNPR.cc \
Modules/MScalarSUN/Div.cc \
Modules/MScalarSUN/TrMag.cc \
Modules/MScalarSUN/TransProj.cc \
Modules/MScalarSUN/TwoPoint.cc \
Modules/MScalarSUN/TrKinetic.cc \
Modules/MScalarSUN/TrPhi.cc \
Modules/MScalarSUN/EMT.cc \
Modules/MScalarSUN/TimeMomProbe.cc \
Modules/MScalarSUN/StochFreeField.cc \
Modules/MScalar/FreeProp.cc \
Modules/MScalar/VPCounterTerms.cc \
Modules/MScalar/ChargedProp.cc \
Modules/MScalar/ScalarVP.cc \
Modules/MAction/Wilson.cc \
Modules/MLoop/NoiseLoop.cc \
Modules/MIO/LoadBinary.cc \
Modules/MIO/LoadCoarseEigenPack.cc \
Modules/MIO/LoadNersc.cc \
Modules/MIO/LoadEigenPack.cc \
Modules/MSink/Smear.cc \
Modules/MSink/Point.cc \
Modules/MFermion/FreeProp.cc \
Modules/MFermion/GaugeProp.cc \
Modules/MGauge/Random.cc \
Modules/MGauge/StochEm.cc \
Modules/MGauge/StoutSmearing.cc \
Modules/MGauge/Unit.cc \
Modules/MGauge/Random.cc \
Modules/MGauge/UnitEm.cc \
Modules/MGauge/FundtoHirep.cc \
Modules/MUtilities/TestSeqGamma.cc \
Modules/MUtilities/TestSeqConserved.cc \
Modules/MSource/Z2.cc \
Modules/MSource/Point.cc \
Modules/MSource/SeqGamma.cc \
Modules/MSource/Wall.cc \
Modules/MSource/SeqConserved.cc \
Modules/MContraction/Meson.cc \
Modules/MContraction/WardIdentity.cc \
Modules/MContraction/WeakHamiltonianNonEye.cc \
Modules/MContraction/Baryon.cc \
Modules/MContraction/DiscLoop.cc \
Modules/MContraction/WeakHamiltonianEye.cc \
Modules/MContraction/A2AMesonField.cc \
Modules/MContraction/WeakNeutral4ptDisc.cc \
Modules/MContraction/Gamma3pt.cc \
Modules/MAction/MobiusDWF.cc \
Modules/MAction/ZMobiusDWF.cc \
Modules/MAction/WilsonClover.cc \
Modules/MAction/Wilson.cc \
Modules/MAction/DWF.cc \
Modules/MAction/ScaledDWF.cc \
Modules/MScalarSUN/TrPhi.cc \
Modules/MScalarSUN/Grad.cc \
Modules/MScalarSUN/TimeMomProbe.cc \
Modules/MScalarSUN/TrMag.cc \
Modules/MScalarSUN/TrKinetic.cc \
Modules/MScalarSUN/EMT.cc \
Modules/MScalarSUN/ShiftProbe.cc \
Modules/MScalarSUN/TransProj.cc \
Modules/MScalarSUN/StochFreeField.cc \
Modules/MScalarSUN/TwoPoint.cc \
Modules/MScalarSUN/TwoPointNPR.cc \
Modules/MScalarSUN/Div.cc \
Modules/MIO/LoadEigenPack.cc \
Modules/MIO/LoadBinary.cc \
Modules/MIO/LoadNersc.cc \
Modules/MIO/LoadCoarseEigenPack.cc
Modules/MAction/ZMobiusDWF.cc \
Modules/MSolver/A2AVectors.cc \
Modules/MSolver/RBPrecCG.cc \
Modules/MSolver/LocalCoherenceLanczos.cc
modules_hpp =\
Modules/MContraction/Baryon.hpp \
Modules/MContraction/Meson.hpp \
Modules/MContraction/WeakHamiltonian.hpp \
Modules/MContraction/WeakHamiltonianNonEye.hpp \
Modules/MContraction/DiscLoop.hpp \
Modules/MContraction/WeakNeutral4ptDisc.hpp \
Modules/MContraction/Gamma3pt.hpp \
Modules/MContraction/WardIdentity.hpp \
Modules/MContraction/WeakHamiltonianEye.hpp \
Modules/MFermion/FreeProp.hpp \
Modules/MFermion/GaugeProp.hpp \
Modules/MSource/SeqGamma.hpp \
Modules/MSource/Point.hpp \
Modules/MSource/Wall.hpp \
Modules/MSource/Z2.hpp \
Modules/MSource/SeqConserved.hpp \
Modules/MSink/Smear.hpp \
Modules/MSink/Point.hpp \
Modules/MSolver/LocalCoherenceLanczos.hpp \
Modules/MSolver/RBPrecCG.hpp \
Modules/MGauge/UnitEm.hpp \
Modules/MGauge/StoutSmearing.hpp \
Modules/MGauge/Unit.hpp \
Modules/MGauge/Random.hpp \
Modules/MGauge/FundtoHirep.hpp \
Modules/MGauge/StochEm.hpp \
Modules/MUtilities/TestSeqGamma.hpp \
Modules/MUtilities/TestSeqConserved.hpp \
Modules/MLoop/NoiseLoop.hpp \
Modules/MScalar/FreeProp.hpp \
Modules/MScalar/VPCounterTerms.hpp \
Modules/MScalar/ScalarVP.hpp \
Modules/MScalar/Scalar.hpp \
Modules/MScalar/ChargedProp.hpp \
Modules/MAction/DWF.hpp \
Modules/MAction/MobiusDWF.hpp \
Modules/MAction/Wilson.hpp \
Modules/MAction/WilsonClover.hpp \
Modules/MAction/ZMobiusDWF.hpp \
Modules/MAction/ScaledDWF.hpp \
Modules/MScalarSUN/TrKinetic.hpp \
Modules/MScalarSUN/TimeMomProbe.hpp \
Modules/MScalarSUN/StochFreeField.hpp \
Modules/MScalarSUN/TwoPointNPR.hpp \
Modules/MScalarSUN/ShiftProbe.hpp \
Modules/MScalarSUN/Grad.hpp \
Modules/MScalarSUN/TransProj.hpp \
Modules/MScalarSUN/Div.hpp \
Modules/MScalarSUN/TimeMomProbe.hpp \
Modules/MScalarSUN/TrMag.hpp \
Modules/MScalarSUN/ShiftProbe.hpp \
Modules/MScalarSUN/Utils.hpp \
Modules/MScalarSUN/EMT.hpp \
Modules/MScalarSUN/TwoPoint.hpp \
Modules/MScalarSUN/TrPhi.hpp \
Modules/MScalarSUN/Utils.hpp \
Modules/MScalarSUN/TransProj.hpp \
Modules/MScalarSUN/Grad.hpp \
Modules/MScalarSUN/TrKinetic.hpp \
Modules/MScalar/FreeProp.hpp \
Modules/MScalar/Scalar.hpp \
Modules/MScalar/ScalarVP.hpp \
Modules/MScalar/ChargedProp.hpp \
Modules/MScalar/VPCounterTerms.hpp \
Modules/MLoop/NoiseLoop.hpp \
Modules/MIO/LoadEigenPack.hpp \
Modules/MIO/LoadNersc.hpp \
Modules/MIO/LoadCoarseEigenPack.hpp \
Modules/MIO/LoadBinary.hpp
Modules/MIO/LoadBinary.hpp \
Modules/MIO/LoadNersc.hpp \
Modules/MSink/Smear.hpp \
Modules/MSink/Point.hpp \
Modules/MFermion/FreeProp.hpp \
Modules/MFermion/GaugeProp.hpp \
Modules/MGauge/FundtoHirep.hpp \
Modules/MGauge/Random.hpp \
Modules/MGauge/StoutSmearing.hpp \
Modules/MGauge/Unit.hpp \
Modules/MGauge/StochEm.hpp \
Modules/MGauge/UnitEm.hpp \
Modules/MUtilities/TestSeqGamma.hpp \
Modules/MUtilities/TestSeqConserved.hpp \
Modules/MSource/SeqConserved.hpp \
Modules/MSource/Z2.hpp \
Modules/MSource/Wall.hpp \
Modules/MSource/SeqGamma.hpp \
Modules/MSource/Point.hpp \
Modules/MContraction/WeakHamiltonianEye.hpp \
Modules/MContraction/Baryon.hpp \
Modules/MContraction/Meson.hpp \
Modules/MContraction/WeakHamiltonian.hpp \
Modules/MContraction/WeakNeutral4ptDisc.hpp \
Modules/MContraction/Gamma3pt.hpp \
Modules/MContraction/DiscLoop.hpp \
Modules/MContraction/WeakHamiltonianNonEye.hpp \
Modules/MContraction/WardIdentity.hpp \
Modules/MContraction/A2AMesonField.hpp \
Modules/MAction/WilsonClover.hpp \
Modules/MAction/ScaledDWF.hpp \
Modules/MAction/MobiusDWF.hpp \
Modules/MAction/Wilson.hpp \
Modules/MAction/DWF.hpp \
Modules/MAction/ZMobiusDWF.hpp \
Modules/MSolver/RBPrecCG.hpp \
Modules/MSolver/LocalCoherenceLanczos.hpp \
Modules/MSolver/A2AVectors.hpp

2
lib/Grid_Eigen_Dense.h
View File

@ -3,7 +3,7 @@
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wdeprecated-declarations"
#endif
#include <Grid/Eigen/Dense>
#include <Eigen/Dense>
#if defined __GNUC__
#pragma GCC diagnostic pop
#endif

2
lib/Makefile.am
View File

@ -34,4 +34,4 @@ HFILES += $(extra_headers)
libGrid_a_SOURCES = $(CCFILES)
libGrid_adir = $(pkgincludedir)
nobase_dist_pkginclude_HEADERS = $(HFILES) $(eigen_files) Config.h
nobase_dist_pkginclude_HEADERS = $(HFILES) $(eigen_files) $(eigen_unsupp_files) Config.h

4
lib/algorithms/iterative/Deflation.h
View File

@ -63,7 +63,7 @@ public:
DeflatedGuesser(const std::vector<Field> & _evec,const std::vector<RealD> & _eval) : evec(_evec), eval(_eval) {};
virtual void operator()(const Field &src,Field &guess) {
virtual void operator()(const Field &src,Field &guess) {
guess = zero;
assert(evec.size()==eval.size());
auto N = evec.size();
@ -71,6 +71,7 @@ public:
const Field& tmp = evec[i];
axpy(guess,TensorRemove(innerProduct(tmp,src)) / eval[i],tmp,guess);
}
guess.checkerboard = src.checkerboard;
}
};
@ -101,6 +102,7 @@ public:
axpy(guess_coarse,TensorRemove(innerProduct(tmp,src_coarse)) / eval_coarse[i],tmp,guess_coarse);
}
blockPromote(guess_coarse,guess,subspace);
guess.checkerboard = src.checkerboard;
};
};

124
lib/algorithms/iterative/SchurRedBlack.h
View File

@ -95,16 +95,26 @@ namespace Grid {
private:
OperatorFunction<Field> & _HermitianRBSolver;
int CBfactorise;
bool subGuess;
public:
/////////////////////////////////////////////////////
// Wrap the usual normal equations Schur trick
/////////////////////////////////////////////////////
SchurRedBlackStaggeredSolve(OperatorFunction<Field> &HermitianRBSolver) :
SchurRedBlackStaggeredSolve(OperatorFunction<Field> &HermitianRBSolver, const bool initSubGuess = false) :
_HermitianRBSolver(HermitianRBSolver)
{
CBfactorise=0;
subtractGuess(initSubGuess);
};
void subtractGuess(const bool initSubGuess)
{
subGuess = initSubGuess;
}
bool isSubtractGuess(void)
{
return subGuess;
}
template<class Matrix>
void operator() (Matrix & _Matrix,const Field &in, Field &out){
@ -150,9 +160,12 @@ namespace Grid {
// Call the red-black solver
//////////////////////////////////////////////////////////////
std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver calling the Mpc solver" <<std::endl;
guess(src_o,sol_o);
guess(src_o, sol_o);
Mtmp = sol_o;
_HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Odd);
std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver called the Mpc solver" <<std::endl;
// Fionn A2A boolean behavioural control
if (subGuess) sol_o = sol_o-Mtmp;
///////////////////////////////////////////////////
// sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
@ -167,11 +180,15 @@ namespace Grid {
std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver inserted solution" <<std::endl;
// Verify the unprec residual
_Matrix.M(out,resid);
resid = resid-in;
RealD ns = norm2(in);
RealD nr = norm2(resid);
std::cout<<GridLogMessage << "SchurRedBlackStaggered solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
if ( ! subGuess ) {
_Matrix.M(out,resid);
resid = resid-in;
RealD ns = norm2(in);
RealD nr = norm2(resid);
std::cout<<GridLogMessage << "SchurRedBlackStaggered solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
} else {
std::cout << GridLogMessage << "Guess subtracted after solve." << std::endl;
}
}
};
template<class Field> using SchurRedBlackStagSolve = SchurRedBlackStaggeredSolve<Field>;
@ -184,15 +201,25 @@ namespace Grid {
private:
OperatorFunction<Field> & _HermitianRBSolver;
int CBfactorise;
bool subGuess;
public:
/////////////////////////////////////////////////////
// Wrap the usual normal equations Schur trick
/////////////////////////////////////////////////////
SchurRedBlackDiagMooeeSolve(OperatorFunction<Field> &HermitianRBSolver,int cb=0) : _HermitianRBSolver(HermitianRBSolver)
SchurRedBlackDiagMooeeSolve(OperatorFunction<Field> &HermitianRBSolver,int cb=0, const bool initSubGuess = false) : _HermitianRBSolver(HermitianRBSolver)
{
CBfactorise=cb;
subtractGuess(initSubGuess);
};
void subtractGuess(const bool initSubGuess)
{
subGuess = initSubGuess;
}
bool isSubtractGuess(void)
{
return subGuess;
}
template<class Matrix>
void operator() (Matrix & _Matrix,const Field &in, Field &out){
ZeroGuesser<Field> guess;
@ -236,7 +263,10 @@ namespace Grid {
//////////////////////////////////////////////////////////////
std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl;
guess(src_o,sol_o);
Mtmp = sol_o;
_HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Odd);
// Fionn A2A boolean behavioural control
if (subGuess) sol_o = sol_o-Mtmp;
///////////////////////////////////////////////////
// sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
@ -249,12 +279,16 @@ namespace Grid {
setCheckerboard(out,sol_o); assert( sol_o.checkerboard ==Odd );
// Verify the unprec residual
_Matrix.M(out,resid);
resid = resid-in;
RealD ns = norm2(in);
RealD nr = norm2(resid);
if ( ! subGuess ) {
_Matrix.M(out,resid);
resid = resid-in;
RealD ns = norm2(in);
RealD nr = norm2(resid);
std::cout<<GridLogMessage << "SchurRedBlackDiagMooee solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
std::cout<<GridLogMessage << "SchurRedBlackDiagMooee solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
} else {
std::cout << GridLogMessage << "Guess subtracted after solve." << std::endl;
}
}
};
@ -267,16 +301,26 @@ namespace Grid {
private:
OperatorFunction<Field> & _HermitianRBSolver;
int CBfactorise;
bool subGuess;
public:
/////////////////////////////////////////////////////
// Wrap the usual normal equations Schur trick
/////////////////////////////////////////////////////
SchurRedBlackDiagTwoSolve(OperatorFunction<Field> &HermitianRBSolver) :
SchurRedBlackDiagTwoSolve(OperatorFunction<Field> &HermitianRBSolver, const bool initSubGuess = false) :
_HermitianRBSolver(HermitianRBSolver)
{
CBfactorise=0;
CBfactorise = 0;
subtractGuess(initSubGuess);
};
void subtractGuess(const bool initSubGuess)
{
subGuess = initSubGuess;
}
bool isSubtractGuess(void)
{
return subGuess;
}
template<class Matrix>
void operator() (Matrix & _Matrix,const Field &in, Field &out){
@ -322,8 +366,11 @@ namespace Grid {
std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl;
// _HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Odd);
guess(src_o,tmp);
Mtmp = tmp;
_HermitianRBSolver(_HermOpEO,src_o,tmp); assert(tmp.checkerboard==Odd);
_Matrix.MooeeInv(tmp,sol_o); assert( sol_o.checkerboard ==Odd);
// Fionn A2A boolean behavioural control
if (subGuess) tmp = tmp-Mtmp;
_Matrix.MooeeInv(tmp,sol_o); assert( sol_o.checkerboard ==Odd);
///////////////////////////////////////////////////
// sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
@ -336,12 +383,16 @@ namespace Grid {
setCheckerboard(out,sol_o); assert( sol_o.checkerboard ==Odd );
// Verify the unprec residual
_Matrix.M(out,resid);
resid = resid-in;
RealD ns = norm2(in);
RealD nr = norm2(resid);
if ( ! subGuess ) {
_Matrix.M(out,resid);
resid = resid-in;
RealD ns = norm2(in);
RealD nr = norm2(resid);
std::cout<<GridLogMessage << "SchurRedBlackDiagTwo solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
std::cout<<GridLogMessage << "SchurRedBlackDiagTwo solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
} else {
std::cout << GridLogMessage << "Guess subtracted after solve." << std::endl;
}
}
};
///////////////////////////////////////////////////////////////////////////////////////////////////////
@ -352,16 +403,26 @@ namespace Grid {
private:
LinearFunction<Field> & _HermitianRBSolver;
int CBfactorise;
bool subGuess;
public:
/////////////////////////////////////////////////////
// Wrap the usual normal equations Schur trick
/////////////////////////////////////////////////////
SchurRedBlackDiagTwoMixed(LinearFunction<Field> &HermitianRBSolver) :
SchurRedBlackDiagTwoMixed(LinearFunction<Field> &HermitianRBSolver, const bool initSubGuess = false) :
_HermitianRBSolver(HermitianRBSolver)
{
CBfactorise=0;
subtractGuess(initSubGuess);
};
void subtractGuess(const bool initSubGuess)
{
subGuess = initSubGuess;
}
bool isSubtractGuess(void)
{
return subGuess;
}
template<class Matrix>
void operator() (Matrix & _Matrix,const Field &in, Field &out){
@ -408,7 +469,10 @@ namespace Grid {
// _HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Odd);
// _HermitianRBSolver(_HermOpEO,src_o,tmp); assert(tmp.checkerboard==Odd);
guess(src_o,tmp);
_HermitianRBSolver(src_o,tmp); assert(tmp.checkerboard==Odd);
Mtmp = tmp;
_HermitianRBSolver(_HermOpEO,src_o,tmp); assert(tmp.checkerboard==Odd);
// Fionn A2A boolean behavioural control
if (subGuess) tmp = tmp-Mtmp;
_Matrix.MooeeInv(tmp,sol_o); assert( sol_o.checkerboard ==Odd);
///////////////////////////////////////////////////
@ -422,12 +486,16 @@ namespace Grid {
setCheckerboard(out,sol_o); assert( sol_o.checkerboard ==Odd );
// Verify the unprec residual
_Matrix.M(out,resid);
resid = resid-in;
RealD ns = norm2(in);
RealD nr = norm2(resid);
if ( ! subGuess ) {
_Matrix.M(out,resid);
resid = resid-in;
RealD ns = norm2(in);
RealD nr = norm2(resid);
std::cout<<GridLogMessage << "SchurRedBlackDiagTwo solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
std::cout << GridLogMessage << "SchurRedBlackDiagTwo solver true unprec resid " << std::sqrt(nr / ns) << " nr " << nr << " ns " << ns << std::endl;
} else {
std::cout << GridLogMessage << "Guess subtracted after solve." << std::endl;
}
}
};

111
lib/lattice/Lattice_reduction.h
View File

@ -39,7 +39,7 @@ template<class vobj> inline RealD norm2(const Lattice<vobj> &arg){
// Double inner product
template<class vobj>
inline ComplexD innerProduct(const Lattice<vobj> &left,const Lattice<vobj> &right)
inline ComplexD innerProduct(const Lattice<vobj> &left,const Lattice<vobj> &right)
{
typedef typename vobj::scalar_type scalar_type;
typedef typename vobj::vector_typeD vector_type;
@ -274,6 +274,115 @@ template<class vobj> inline void sliceSum(const Lattice<vobj> &Data,std::vector<
}
}
template<class vobj>
static void mySliceInnerProductVector( std::vector<ComplexD> & result, const Lattice<vobj> &lhs,const Lattice<vobj> &rhs,int orthogdim)
{
// std::cout << GridLogMessage << "Start mySliceInnerProductVector" << std::endl;
typedef typename vobj::scalar_type scalar_type;
std::vector<scalar_type> lsSum;
localSliceInnerProductVector(result, lhs, rhs, lsSum, orthogdim);
globalSliceInnerProductVector(result, lhs, lsSum, orthogdim);
// std::cout << GridLogMessage << "End mySliceInnerProductVector" << std::endl;
}
template <class vobj>
static void localSliceInnerProductVector(std::vector<ComplexD> &result, const Lattice<vobj> &lhs, const Lattice<vobj> &rhs, std::vector<typename vobj::scalar_type> &lsSum, int orthogdim)
{
// std::cout << GridLogMessage << "Start prep" << std::endl;
typedef typename vobj::vector_type vector_type;
typedef typename vobj::scalar_type scalar_type;
GridBase *grid = lhs._grid;
assert(grid!=NULL);
conformable(grid,rhs._grid);
const int Nd = grid->_ndimension;
const int Nsimd = grid->Nsimd();
assert(orthogdim >= 0);
assert(orthogdim < Nd);
int fd=grid->_fdimensions[orthogdim];
int ld=grid->_ldimensions[orthogdim];
int rd=grid->_rdimensions[orthogdim];
// std::cout << GridLogMessage << "Start alloc" << std::endl;
std::vector<vector_type,alignedAllocator<vector_type> > lvSum(rd); // will locally sum vectors first
lsSum.resize(ld,scalar_type(0.0)); // sum across these down to scalars
std::vector<iScalar<scalar_type>> extracted(Nsimd); // splitting the SIMD
// std::cout << GridLogMessage << "End alloc" << std::endl;
result.resize(fd); // And then global sum to return the same vector to every node for IO to file
for(int r=0;r<rd;r++){
lvSum[r]=zero;
}
int e1= grid->_slice_nblock[orthogdim];
int e2= grid->_slice_block [orthogdim];
int stride=grid->_slice_stride[orthogdim];
// std::cout << GridLogMessage << "End prep" << std::endl;
// std::cout << GridLogMessage << "Start parallel inner product, _rd = " << rd << std::endl;
vector_type vv;
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;
vv = TensorRemove(innerProduct(lhs._odata[ss], rhs._odata[ss]));
lvSum[r] = lvSum[r] + vv;
}
}
}
// std::cout << GridLogMessage << "End parallel inner product" << std::endl;
// Sum across simd lanes in the plane, breaking out orthog dir.
std::vector<int> icoor(Nd);
for(int rt=0;rt<rd;rt++){
iScalar<vector_type> temp;
temp._internal = lvSum[rt];
extract(temp,extracted);
for(int idx=0;idx<Nsimd;idx++){
grid->iCoorFromIindex(icoor,idx);
int ldx =rt+icoor[orthogdim]*rd;
lsSum[ldx]=lsSum[ldx]+extracted[idx]._internal;
}
}
// std::cout << GridLogMessage << "End sum over simd lanes" << std::endl;
}
template <class vobj>
static void globalSliceInnerProductVector(std::vector<ComplexD> &result, const Lattice<vobj> &lhs, std::vector<typename vobj::scalar_type> &lsSum, int orthogdim)
{
typedef typename vobj::scalar_type scalar_type;
GridBase *grid = lhs._grid;
int fd = result.size();
int ld = lsSum.size();
// sum over nodes.
std::vector<scalar_type> gsum;
gsum.resize(fd, scalar_type(0.0));
// std::cout << GridLogMessage << "Start of gsum[t] creation:" << std::endl;
for(int t=0;t<fd;t++){
int pt = t/ld; // processor plane
int lt = t%ld;
if ( pt == grid->_processor_coor[orthogdim] ) {
gsum[t]=lsSum[lt];
}
}
// std::cout << GridLogMessage << "End of gsum[t] creation:" << std::endl;
// std::cout << GridLogMessage << "Start of GlobalSumVector:" << std::endl;
grid->GlobalSumVector(&gsum[0], fd);
// std::cout << GridLogMessage << "End of GlobalSumVector:" << std::endl;
result = gsum;
}
template<class vobj>
static void sliceInnerProductVector( std::vector<ComplexD> & result, const Lattice<vobj> &lhs,const Lattice<vobj> &rhs,int orthogdim)
{

27
lib/qcd/action/fermion/CayleyFermion5D.cc
View File

@ -67,6 +67,33 @@ void CayleyFermion5D<Impl>::ExportPhysicalFermionSolution(const FermionField &so
axpby_ssp_pplus (tmp, 1., tmp , 1., solution5d, 0, Ls-1);
ExtractSlice(exported4d, tmp, 0, 0);
}
template<class Impl>
void CayleyFermion5D<Impl>::ExportPhysicalFermionSource(const FermionField &solution5d,FermionField &exported4d)
{
int Ls = this->Ls;
FermionField tmp(this->FermionGrid());
tmp = solution5d;
conformable(solution5d._grid,this->FermionGrid());
conformable(exported4d._grid,this->GaugeGrid());
axpby_ssp_pplus (tmp, 0., solution5d, 1., solution5d, 0, 0);
axpby_ssp_pminus(tmp, 1., tmp , 1., solution5d, 0, Ls-1);
ExtractSlice(exported4d, tmp, 0, 0);
}
template<class Impl>
void CayleyFermion5D<Impl>::ImportUnphysicalFermion(const FermionField &input4d,FermionField &imported5d)
{
int Ls = this->Ls;
FermionField tmp(this->FermionGrid());
conformable(imported5d._grid,this->FermionGrid());
conformable(input4d._grid ,this->GaugeGrid());
tmp = zero;
InsertSlice(input4d, tmp, 0 , 0);
InsertSlice(input4d, tmp, Ls-1, 0);
axpby_ssp_pplus (tmp, 0., tmp, 1., tmp, 0, 0);
axpby_ssp_pminus(tmp, 0., tmp, 1., tmp, Ls-1, Ls-1);
imported5d=tmp;
}
template<class Impl>
void CayleyFermion5D<Impl>::ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d)
{

2
lib/qcd/action/fermion/CayleyFermion5D.h
View File

@ -89,7 +89,9 @@ namespace Grid {
virtual void Dminus(const FermionField &psi, FermionField &chi);
virtual void DminusDag(const FermionField &psi, FermionField &chi);
virtual void ExportPhysicalFermionSolution(const FermionField &solution5d,FermionField &exported4d);
virtual void ExportPhysicalFermionSource(const FermionField &solution5d, FermionField &exported4d);
virtual void ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d);
virtual void ImportUnphysicalFermion(const FermionField &solution5d, FermionField &exported4d);
/////////////////////////////////////////////////////
// Instantiate different versions depending on Impl

8
lib/qcd/action/fermion/FermionOperator.h
View File

@ -162,10 +162,18 @@ namespace Grid {
{
imported = input;
};
virtual void ImportUnphysicalFermion(const FermionField &input,FermionField &imported)
{
imported=input;
};
virtual void ExportPhysicalFermionSolution(const FermionField &solution,FermionField &exported)
{
exported=solution;
};
virtual void ExportPhysicalFermionSource(const FermionField &solution,FermionField &exported)
{
exported=solution;
};
};
}

2
lib/qcd/action/fermion/WilsonCloverFermion.cc
View File

@ -27,7 +27,7 @@
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
#include <Grid/Eigen/Dense>
//#include <Grid/Eigen/Dense>
#include <Grid/qcd/spin/Dirac.h>
namespace Grid

45
scripts/update_eigen.sh
View File

@ -7,13 +7,42 @@ fi
ARC=$1
INITDIR=`pwd`
rm -rf lib/Eigen
ARCDIR=`tar -tf ${ARC} | head -n1 | sed -e 's@/.*@@'`
##################
#untar
##################
tar -xf ${ARC}
cd ${ARCDIR}
(tar -cf - Eigen --exclude='*.txt' 2>/dev/null) | tar -xf - -C ../lib/
cd ../lib
echo 'eigen_files =\' > Eigen.inc
find Eigen -type f -print | sed 's/^/ /;$q;s/$/ \\/' >> Eigen.inc
ARCDIR=`tar -tf ${ARC} | head -n1 | sed -e 's@/.*@@'`
rm -f ${ARC}
###############################
# Link to a deterministic name
###############################
mv ${ARCDIR} Eigen
# Eigen source headers
cd ${INITDIR}/Eigen
echo 'eigen_files =\' > ${INITDIR}/lib/Eigen.inc
find Eigen -name "*.h" -print | sed 's/^/ /;$q;s/$/ \\/' >> ${INITDIR}/lib/Eigen.inc
cd ${INITDIR}
rm -rf ${ARCDIR}
echo 'eigen_unsupp_files =\' >> ${INITDIR}/lib/Eigen.inc
find Eigen/unsupported/Eigen -name "*.h" -print | sed 's/^/ /;$q;s/$/ \\/' >> ${INITDIR}/lib/Eigen.inc
###################################
# back to home
###################################
cd ${INITDIR}
#########################################
# Make grid includes happy
#########################################
mkdir ${INITDIR}/lib/Eigen/
ln -s ${INITDIR}/Eigen/Eigen/* ${INITDIR}/lib/Eigen/
ln -s ${INITDIR}/Eigen/unsupported ${INITDIR}/lib/Eigen/