From 62febd2823e63d88cc9c2041120a0d240cd71c8e Mon Sep 17 00:00:00 2001 From: paboyle Date: Wed, 31 Aug 2016 00:23:09 +0100 Subject: [PATCH 01/58] Wilson prop test --- tests/core/Test_fft.cc | 86 +++++++++++++++++++++++++++++++----------- 1 file changed, 64 insertions(+), 22 deletions(-) diff --git a/tests/core/Test_fft.cc b/tests/core/Test_fft.cc index 2fdaeed2..c076695d 100644 --- a/tests/core/Test_fft.cc +++ b/tests/core/Test_fft.cc @@ -46,16 +46,19 @@ int main (int argc, char ** argv) for(int d=0;d p({1,2,3,2}); @@ -72,31 +75,41 @@ int main (int argc, char ** argv) } C = exp(C*ci); - + Csav = C; S=zero; S = S+C; - FFT theFFT(&Fine); + FFT theFFT(&GRID); + + theFFT.FFT_dim(Ctilde,C,0,FFT::forward); C=Ctilde; + theFFT.FFT_dim(Ctilde,C,1,FFT::forward); C=Ctilde; std::cout << theFFT.MFlops()<<" Mflops "< seeds({1,2,3,4}); + GridParallelRNG pRNG(&GRID); + pRNG.SeedFixedIntegers(seeds); + + LatticeGaugeFieldD Umu(&GRID); + LatticeFermionD src(&GRID); gaussian(pRNG,src); + LatticeFermionD tmp(&GRID); + LatticeFermionD ref(&GRID); + + SU3::ColdConfiguration(pRNG,Umu); // Unit gauge + + RealD mass=0.1; + WilsonFermionD Dw(Umu,GRID,RBGRID,mass); + + Dw.M(src,tmp); + + std::cout << " src = " < ssource(vgrid); ssource =source; Lattice pgsource(&pencil_g); @@ -184,9 +214,10 @@ namespace Grid { istride = ostride = Ncomp*Nlow; /* distance between two elements in the same FT */ int *inembed = n, *onembed = n; - - int sign = FFTW_FORWARD; - if (inverse) sign = FFTW_BACKWARD; + scalar div; + if ( sign == backward ) div = 1.0/G; + else if ( sign == forward ) div = 1.0; + else assert(0); FFTW_plan p; { @@ -258,6 +289,7 @@ PARALLEL_FOR_LOOP sobj s; gcoor[dim] = lcoor[dim]+L*pc; peekLocalSite(s,pgresult,gcoor); + s = s * div; pokeLocalSite(s,result,lcoor); } From d15ab66aae0909685bb163c1738884e5e856f370 Mon Sep 17 00:00:00 2001 From: paboyle Date: Wed, 31 Aug 2016 00:25:22 +0100 Subject: [PATCH 05/58] FFT moves higher in include order --- lib/Grid.h | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/lib/Grid.h b/lib/Grid.h index 486ee4d3..0c5983f3 100644 --- a/lib/Grid.h +++ b/lib/Grid.h @@ -77,11 +77,10 @@ Author: paboyle #include #include #include -#include -#include - #include +#include +#include #include #include From 02e983a0cd5b13a34238598824867e50f47fece2 Mon Sep 17 00:00:00 2001 From: paboyle Date: Wed, 31 Aug 2016 00:26:02 +0100 Subject: [PATCH 06/58] Momentum space prop and free prop convolution --- lib/qcd/action/fermion/FermionOperator.h | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h index ea5583eb..daee2e5b 100644 --- a/lib/qcd/action/fermion/FermionOperator.h +++ b/lib/qcd/action/fermion/FermionOperator.h @@ -91,6 +91,19 @@ namespace Grid { virtual void Mdiag (const FermionField &in, FermionField &out) { Mooee(in,out);}; // Same as Mooee applied to both CB's virtual void Mdir (const FermionField &in, FermionField &out,int dir,int disp)=0; // case by case Wilson, Clover, Cayley, ContFrac, PartFrac + + virtual void MomentumSpacePropagator(FermionField &out,const FermionField &in) { assert(0);}; + virtual void FreePropagator(const FermionField &in,FermionField &out) { + FFT theFFT((GridCartesian *) in._grid); + + FermionField in_k(in._grid); + FermionField prop_k(in._grid); + + theFFT.FFT_all_dim(in_k,in,FFT::forward); + this->MomentumSpacePropagator(prop_k,in_k); + theFFT.FFT_all_dim(out,prop_k,FFT::backward); + }; + /////////////////////////////////////////////// // Updates gauge field during HMC /////////////////////////////////////////////// From 0c1d7e4daf4ad2a7e41a6e5f689bb07682619f58 Mon Sep 17 00:00:00 2001 From: paboyle Date: Wed, 31 Aug 2016 00:26:36 +0100 Subject: [PATCH 07/58] Mom space prop for Wilson action --- lib/qcd/action/fermion/WilsonFermion.cc | 52 +++++++++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index 59632409..caf1913d 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -136,6 +136,58 @@ namespace QCD { out.checkerboard = in.checkerboard; MooeeInv(in,out); } + template + void WilsonFermion:: MomentumSpacePropagator(FermionField &out, const FermionField &in) { + + // what type LatticeComplex + conformable(_grid,out._grid); + + typedef typename FermionField::vector_type vector_type; + typedef typename FermionField::scalar_type ScalComplex; + + typedef Lattice > LatComplex; + + Gamma::GammaMatrix Gmu [] = { + Gamma::GammaX, + Gamma::GammaY, + Gamma::GammaZ, + Gamma::GammaT + }; + + std::vector latt_size = _grid->_fdimensions; + + FermionField num (_grid); num = zero; + LatComplex wilson(_grid); wilson= zero; + LatComplex one (_grid); one = ScalComplex(1.0,0.0); + + LatComplex denom(_grid); denom= zero; + LatComplex kmu(_grid); + ScalComplex ci(0.0,1.0); + // momphase = n * 2pi / L + for(int mu=0;mu Date: Wed, 31 Aug 2016 00:27:04 +0100 Subject: [PATCH 08/58] Wilson tree level added --- lib/qcd/action/fermion/WilsonFermion.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/lib/qcd/action/fermion/WilsonFermion.h b/lib/qcd/action/fermion/WilsonFermion.h index 3de2cac4..c050c0eb 100644 --- a/lib/qcd/action/fermion/WilsonFermion.h +++ b/lib/qcd/action/fermion/WilsonFermion.h @@ -78,6 +78,8 @@ namespace Grid { virtual void MooeeInv(const FermionField &in, FermionField &out) ; virtual void MooeeInvDag(const FermionField &in, FermionField &out) ; + virtual void MomentumSpacePropagator(FermionField &out,const FermionField &in) ; + //////////////////////// // Derivative interface //////////////////////// From 8535d433a745cb50e2fd0dc5f46cf01ee3a61446 Mon Sep 17 00:00:00 2001 From: paboyle Date: Wed, 31 Aug 2016 00:27:53 +0100 Subject: [PATCH 09/58] Cold or hot must support any precisoin --- lib/qcd/utils/SUn.h | 21 ++++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/lib/qcd/utils/SUn.h b/lib/qcd/utils/SUn.h index e9403836..67da9787 100644 --- a/lib/qcd/utils/SUn.h +++ b/lib/qcd/utils/SUn.h @@ -617,22 +617,33 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g PokeIndex(out,Umu,mu); } } - static void TepidConfiguration(GridParallelRNG &pRNG,LatticeGaugeField &out){ - LatticeMatrix Umu(out._grid); + template + static void TepidConfiguration(GridParallelRNG &pRNG,GaugeField &out){ + typedef typename GaugeField::vector_type vector_type; + typedef iSUnMatrix vMatrixType; + typedef Lattice LatticeMatrixType; + + LatticeMatrixType Umu(out._grid); for(int mu=0;mu(out,Umu,mu); } } - static void ColdConfiguration(GridParallelRNG &pRNG,LatticeGaugeField &out){ - LatticeMatrix Umu(out._grid); + template + static void ColdConfiguration(GridParallelRNG &pRNG,GaugeField &out){ + typedef typename GaugeField::vector_type vector_type; + typedef iSUnMatrix vMatrixType; + typedef Lattice LatticeMatrixType; + + LatticeMatrixType Umu(out._grid); Umu=1.0; for(int mu=0;mu(out,Umu,mu); } } - static void taProj( const LatticeMatrix &in, LatticeMatrix &out){ + template + static void taProj( const LatticeMatrixType &in, LatticeMatrixType &out){ out = Ta(in); } template From 7422953e36f895f1a7499452461e0e2c4ef4716c Mon Sep 17 00:00:00 2001 From: paboyle Date: Wed, 31 Aug 2016 00:42:47 +0100 Subject: [PATCH 10/58] Poisson solver example --- tests/core/Test_poisson_fft.cc | 138 +++++++++++++++++++++++++++++++++ 1 file changed, 138 insertions(+) create mode 100644 tests/core/Test_poisson_fft.cc diff --git a/tests/core/Test_poisson_fft.cc b/tests/core/Test_poisson_fft.cc new file mode 100644 index 00000000..2252ec7a --- /dev/null +++ b/tests/core/Test_poisson_fft.cc @@ -0,0 +1,138 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_poisson_fft.cc + + Copyright (C) 2015 + +Author: Azusa Yamaguchi +Author: Peter Boyle + + 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 + +using namespace Grid; +using namespace Grid::QCD; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + int threads = GridThread::GetThreads(); + std::cout< latt_size ({N,N}); + std::vector simd_layout({vComplexD::Nsimd(),1}); + std::vector mpi_layout ({1,1}); + + int vol = 1; + int nd = latt_size.size(); + for(int d=0;dInteger(N2+W),zz,Charge); + } + + // std::cout << Charge < k(4,&GRID); + LatticeComplexD ksq(&GRID); + + ksq=zero; + for(int mu=0;mu (L/2), k[mu]-L, k[mu]); + + // std::cout << k[mu]< zero_mode(nd,0); + TComplexD Tone = ComplexD(1.0,0.0); + pokeSite(Tone,ksq,zero_mode); + + // std::cout << "Charge\n" << Charge < Date: Mon, 26 Sep 2016 09:36:14 +0100 Subject: [PATCH 11/58] FPE's on macos set up --- lib/Init.cc | 39 +++++++++++++++++++++++---------------- lib/Init.h | 1 - 2 files changed, 23 insertions(+), 17 deletions(-) diff --git a/lib/Init.cc b/lib/Init.cc index 4cb4bb50..cff43185 100644 --- a/lib/Init.cc +++ b/lib/Init.cc @@ -46,23 +46,10 @@ Author: paboyle #include #include #include -#include namespace Grid { - std::string demangle(const char* name) { - - int status = -4; // some arbitrary value to eliminate the compiler warning - - // enable c++11 by passing the flag -std=c++11 to g++ - std::unique_ptr res { - abi::__cxa_demangle(name, NULL, NULL, &status), - std::free - }; - - return (status==0) ? res.get() : name ; - } ////////////////////////////////////////////////////// // Convenience functions to access stadard command line arg @@ -343,10 +330,30 @@ void Grid_sa_signal_handler(int sig,siginfo_t *si,void * ptr) exit(0); return; }; -#ifdef GRID_FPE + + #define _GNU_SOURCE #include + +#ifdef __APPLE__ +static int +feenableexcept (unsigned int excepts) +{ + static fenv_t fenv; + unsigned int new_excepts = excepts & FE_ALL_EXCEPT, + old_excepts; // previous masks + + if ( fegetenv (&fenv) ) return -1; + old_excepts = fenv.__control & FE_ALL_EXCEPT; + + // unmask + fenv.__control &= ~new_excepts; + fenv.__mxcsr &= ~(new_excepts << 7); + + return ( fesetenv (&fenv) ? -1 : old_excepts ); +} #endif + void Grid_debug_handler_init(void) { struct sigaction sa,osa; @@ -355,9 +362,9 @@ void Grid_debug_handler_init(void) sa.sa_flags = SA_SIGINFO; sigaction(SIGSEGV,&sa,NULL); sigaction(SIGTRAP,&sa,NULL); -#ifdef GRID_FPE + feenableexcept( FE_INVALID|FE_OVERFLOW|FE_DIVBYZERO); + sigaction(SIGFPE,&sa,NULL); -#endif } } diff --git a/lib/Init.h b/lib/Init.h index 6eb07a49..8986d24b 100644 --- a/lib/Init.h +++ b/lib/Init.h @@ -52,7 +52,6 @@ namespace Grid { void GridCmdOptionCSL(std::string str,std::vector & vec); void GridCmdOptionIntVector(std::string &str,std::vector & vec); - std::string demangle(const char* name) ; void GridParseLayout(char **argv,int argc, std::vector &latt, From 567b6cf23f5b33bb101965fc7613b969373db20a Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 26 Sep 2016 09:36:51 +0100 Subject: [PATCH 12/58] demangle moves to logging --- lib/Log.cc | 15 +++++++++++++++ lib/Log.h | 3 ++- 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/lib/Log.cc b/lib/Log.cc index 2082570d..5152bd77 100644 --- a/lib/Log.cc +++ b/lib/Log.cc @@ -31,8 +31,23 @@ directory /* END LEGAL */ #include +#include + namespace Grid { + std::string demangle(const char* name) { + + int status = -4; // some arbitrary value to eliminate the compiler warning + + // enable c++11 by passing the flag -std=c++11 to g++ + std::unique_ptr res { + abi::__cxa_demangle(name, NULL, NULL, &status), + std::free + }; + + return (status==0) ? res.get() : name ; + } + GridStopWatch Logger::StopWatch; std::ostream Logger::devnull(0); diff --git a/lib/Log.h b/lib/Log.h index 156f52ee..9f9e19a7 100644 --- a/lib/Log.h +++ b/lib/Log.h @@ -143,6 +143,7 @@ extern GridLogger GridLogIterative ; extern GridLogger GridLogIntegrator ; extern Colours GridLogColours; + std::string demangle(const char* name) ; #define _NBACKTRACE (256) extern void * Grid_backtrace_buffer[_NBACKTRACE]; @@ -161,7 +162,7 @@ std::fclose(fp); \ int symbols = backtrace (Grid_backtrace_buffer,_NBACKTRACE);\ char **strings = backtrace_symbols(Grid_backtrace_buffer,symbols);\ for (int i = 0; i < symbols; i++){\ - std::fprintf (fp,"BackTrace Strings: %d %s\n",i, strings[i]); std::fflush(fp); \ + std::fprintf (fp,"BackTrace Strings: %d %s\n",i, demangle(strings[i]).c_str()); std::fflush(fp); \ }\ } #else From 16b37b956cb8018ba9acfab1424305854dde750a Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 26 Sep 2016 09:37:59 +0100 Subject: [PATCH 13/58] divide goes to ET --- lib/lattice/Lattice_base.h | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/lib/lattice/Lattice_base.h b/lib/lattice/Lattice_base.h index 3bfa5613..b97e1359 100644 --- a/lib/lattice/Lattice_base.h +++ b/lib/lattice/Lattice_base.h @@ -303,17 +303,6 @@ PARALLEL_FOR_LOOP *this = (*this)+r; return *this; } - - strong_inline friend Lattice operator / (const Lattice &lhs,const Lattice &rhs){ - conformable(lhs,rhs); - Lattice ret(lhs._grid); -PARALLEL_FOR_LOOP - for(int ss=0;ssoSites();ss++){ - ret._odata[ss] = lhs._odata[ss]*pow(rhs._odata[ss],-1.0); - } - return ret; - }; - }; // class Lattice template std::ostream& operator<< (std::ostream& stream, const Lattice &o){ From 81a7a0307698246638776d1fad7a6610b6116e83 Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 26 Sep 2016 09:38:17 +0100 Subject: [PATCH 14/58] Integer << --- lib/Simd.h | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/lib/Simd.h b/lib/Simd.h index 9568c555..adc2849d 100644 --- a/lib/Simd.h +++ b/lib/Simd.h @@ -237,6 +237,18 @@ namespace Grid { stream<<">"; return stream; } + inline std::ostream& operator<< (std::ostream& stream, const vInteger &o){ + int nn=vInteger::Nsimd(); + std::vector > buf(nn); + vstore(o,&buf[0]); + stream<<"<"; + for(int i=0;i"; + return stream; + } } From 52a39f0fcdcdf152b495a86630f377e92277985b Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 26 Sep 2016 09:38:38 +0100 Subject: [PATCH 15/58] Divide in ET --- lib/lattice/Lattice_ET.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/lib/lattice/Lattice_ET.h b/lib/lattice/Lattice_ET.h index 7ebac99d..1bb83901 100644 --- a/lib/lattice/Lattice_ET.h +++ b/lib/lattice/Lattice_ET.h @@ -261,6 +261,7 @@ GridUnopClass(UnaryExp, exp(a)); GridBinOpClass(BinaryAdd, lhs + rhs); GridBinOpClass(BinarySub, lhs - rhs); GridBinOpClass(BinaryMul, lhs *rhs); +GridBinOpClass(BinaryDiv, lhs /rhs); GridBinOpClass(BinaryAnd, lhs &rhs); GridBinOpClass(BinaryOr, lhs | rhs); @@ -385,6 +386,7 @@ GRID_DEF_UNOP(exp, UnaryExp); GRID_DEF_BINOP(operator+, BinaryAdd); GRID_DEF_BINOP(operator-, BinarySub); GRID_DEF_BINOP(operator*, BinaryMul); +GRID_DEF_BINOP(operator/, BinaryDiv); GRID_DEF_BINOP(operator&, BinaryAnd); GRID_DEF_BINOP(operator|, BinaryOr); From b6713ecb60c2a0095ae4af44911b4bd6600d03a8 Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 26 Sep 2016 09:39:09 +0100 Subject: [PATCH 16/58] Momentum space rules for Overlap, DWF untested to date --- lib/qcd/action/fermion/WilsonFermion.cc | 2 + lib/qcd/action/fermion/WilsonFermion5D.cc | 164 ++++++++++++++++++++++ lib/qcd/action/fermion/WilsonFermion5D.h | 3 + 3 files changed, 169 insertions(+) diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index caf1913d..55fbfa01 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -184,7 +184,9 @@ namespace QCD { num = num + wilson*in; // -i gmu sin k + 2 sin^2 k/2 + m denom= denom+wilson*wilson; // sin^2 k + (2 sin^2 k/2 + m)^2 + denom= one/denom; + out = num*denom; // [ -i gmu sin k + 2 sin^2 k/2 + m] / [ sin^2 k + (2 sin^2 k/2 + m)^2 ] } diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index ef7e84ab..28f02b67 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -392,6 +392,170 @@ void WilsonFermion5D::DW(const FermionField &in, FermionField &out,int dag axpy(out,4.0-M5,in,out); } +template +void WilsonFermion5D::MomentumSpacePropagatorHt(FermionField &out,const FermionField &in) +{ + // what type LatticeComplex + GridBase *_grid = _FourDimGrid; + conformable(_grid,out._grid); + + typedef typename FermionField::vector_type vector_type; + typedef typename FermionField::scalar_type ScalComplex; + typedef iSinglet Tcomplex; + typedef Lattice > LatComplex; + + Gamma::GammaMatrix Gmu [] = { + Gamma::GammaX, + Gamma::GammaY, + Gamma::GammaZ, + Gamma::GammaT + }; + + std::vector latt_size = _grid->_fdimensions; + + + LatComplex sk(_grid); sk = zero; + LatComplex sk2(_grid); sk2= zero; + + LatComplex W(_grid); W= zero; + LatComplex a(_grid); a= zero; + LatComplex cosha(_grid); a= zero; + + LatComplex one (_grid); one = ScalComplex(1.0,0.0); + + FermionField num (_grid); num = zero; + LatComplex denom(_grid); denom= zero; + LatComplex kmu(_grid); + ScalComplex ci(0.0,1.0); + + for(int mu=0;mulSites();idx++){ + Tcomplex cc; + RealD sgn; + std::vector lcoor(Nd); + _grid->LocalIndexToLocalCoor(idx,lcoor); + peekLocalSite(cc,cosha,lcoor); + sgn= ::fabs(real(cc)); + cc = ScalComplex(::acosh(sgn),0.0); + pokeLocalSite(cc,a,lcoor); + } + + std::cout << "Ht mom space num" << norm2(num)< +void WilsonFermion5D::MomentumSpacePropagatorHw(FermionField &out,const FermionField &in) +{ + // what type LatticeComplex + GridBase *_grid = _FourDimGrid; + conformable(_grid,out._grid); + + typedef typename FermionField::vector_type vector_type; + typedef typename FermionField::scalar_type ScalComplex; + + typedef Lattice > LatComplex; + + Gamma::GammaMatrix Gmu [] = { + Gamma::GammaX, + Gamma::GammaY, + Gamma::GammaZ, + Gamma::GammaT + }; + + std::vector latt_size = _grid->_fdimensions; + + + LatComplex sk(_grid); sk = zero; + LatComplex sk2(_grid); sk2= zero; + + LatComplex w_k(_grid); w_k= zero; + LatComplex b_k(_grid); b_k= zero; + + LatComplex one (_grid); one = ScalComplex(1.0,0.0); + + FermionField num (_grid); num = zero; + LatComplex denom(_grid); denom= zero; + LatComplex kmu(_grid); + ScalComplex ci(0.0,1.0); + + for(int mu=0;mu Date: Mon, 26 Sep 2016 09:42:22 +0100 Subject: [PATCH 17/58] Divide handling improved --- lib/simd/Grid_avx.h | 17 ++++++- lib/simd/Grid_avx512.h | 12 +++++ lib/simd/Grid_imci.h | 12 +++++ lib/simd/Grid_sse4.h | 15 ++++++ lib/simd/Grid_vector_types.h | 90 ++++++++++++++++++++++-------------- 5 files changed, 110 insertions(+), 36 deletions(-) diff --git a/lib/simd/Grid_avx.h b/lib/simd/Grid_avx.h index 03faabee..33496e07 100644 --- a/lib/simd/Grid_avx.h +++ b/lib/simd/Grid_avx.h @@ -365,6 +365,18 @@ namespace Optimization { } }; + struct Div{ + // Real float + inline __m256 operator()(__m256 a, __m256 b){ + return _mm256_div_ps(a,b); + } + // Real double + inline __m256d operator()(__m256d a, __m256d b){ + return _mm256_div_pd(a,b); + } + }; + + struct Conj{ // Complex single inline __m256 operator()(__m256 in){ @@ -473,7 +485,7 @@ namespace Optimization { } #endif - + /* inline std::ostream & operator << (std::ostream& stream, const __m256 a) { const float *p=(const float *)&a; @@ -486,7 +498,7 @@ namespace Optimization { stream<< "{"< > { ////////////////////////////////////// // demote a vector to real type ////////////////////////////////////// - // type alias used to simplify the syntax of std::enable_if -template -using Invoke = typename T::type; -template -using EnableIf = Invoke >; -template -using NotEnableIf = Invoke >; +template using Invoke = typename T::type; +template using EnableIf = Invoke >; +template using NotEnableIf = Invoke >; //////////////////////////////////////////////////////// // Check for complexity with type traits -template -struct is_complex : public std::false_type {}; -template <> -struct is_complex > : public std::true_type {}; -template <> -struct is_complex > : public std::true_type {}; +template struct is_complex : public std::false_type {}; +template <> struct is_complex > : public std::true_type {}; +template <> struct is_complex > : public std::true_type {}; -template -using IfReal = Invoke::value, int> >; -template -using IfComplex = Invoke::value, int> >; -template -using IfInteger = Invoke::value, int> >; +template using IfReal = Invoke::value, int> >; +template using IfComplex = Invoke::value, int> >; +template using IfInteger = Invoke::value, int> >; -template -using IfNotReal = - Invoke::value, int> >; -template -using IfNotComplex = Invoke::value, int> >; -template -using IfNotInteger = Invoke::value, int> >; +template using IfNotReal = Invoke::value, int> >; +template using IfNotComplex = Invoke::value, int> >; +template using IfNotInteger = Invoke::value, int> >; //////////////////////////////////////////////////////// // Define the operation templates functors @@ -285,6 +271,20 @@ class Grid_simd { return a * b; } + ////////////////////////////////// + // Divides + ////////////////////////////////// + friend inline Grid_simd operator/(const Scalar_type &a, Grid_simd b) { + Grid_simd va; + vsplat(va, a); + return va / b; + } + friend inline Grid_simd operator/(Grid_simd b, const Scalar_type &a) { + Grid_simd va; + vsplat(va, a); + return b / a; + } + /////////////////////// // Unary negation /////////////////////// @@ -428,7 +428,6 @@ inline void rotate(Grid_simd &ret,Grid_simd b,int nrot) ret.v = Optimization::Rotate::rotate(b.v,2*nrot); } - template inline void vbroadcast(Grid_simd &ret,const Grid_simd &src,int lane){ S* typepun =(S*) &src; @@ -512,7 +511,6 @@ template = 0> inline void vfalse(Grid_simd &ret) { vsplat(ret, 0); } - template inline void zeroit(Grid_simd &z) { vzero(z); @@ -530,7 +528,6 @@ inline void vstream(Grid_simd &out, const Grid_simd &in) { typedef typename S::value_type T; binary((T *)&out.v, in.v, VstreamSIMD()); } - template = 0> inline void vstream(Grid_simd &out, const Grid_simd &in) { out = in; @@ -569,6 +566,36 @@ inline Grid_simd operator*(Grid_simd a, Grid_simd b) { return ret; }; +// Distinguish between complex types and others +template = 0> +inline Grid_simd operator/(Grid_simd a, Grid_simd b) { + typedef Grid_simd simd; + + simd ret; + simd den; + typename simd::conv_t conv; + + ret = a * conjugate(b) ; + den = b * conjugate(b) ; + + + auto real_den = toReal(den); + + ret.v=binary(ret.v, real_den.v, DivSIMD()); + + std::cout << " Div call from complex "< + class Photon { + + public: + + INHERIT_GIMPL_TYPES(Gimpl); + + enum PhotonType { FEYNMAN_L, FEYNMAN_TL }; + + PhotonType GaugeBC; + + Photon(PhotonType _GaugeBC) : GaugeBC(_GaugeBC){}; + + void FreePropagator (const GaugeField &in,GaugeField &out) { + FFT theFFT((GridCartesian *) in._grid); + + GaugeField in_k(in._grid); + GaugeField prop_k(in._grid); + + theFFT.FFT_all_dim(in_k,in,FFT::forward); + MomentumSpacePropagator(prop_k,in_k); + theFFT.FFT_all_dim(out,prop_k,FFT::backward); + } + void MomentumSpacePropagator(GaugeField &out,const GaugeField &in) { + if ( GaugeBC == FEYNMAN_L ) { + FeynmanGaugeMomentumSpacePropagator_L(out,in); + } else if ( GaugeBC == FEYNMAN_TL ) { + FeynmanGaugeMomentumSpacePropagator_TL(out,in); + } else { // No coulomb yet + assert(0); + } + } + void FeynmanGaugeMomentumSpacePropagator_L(GaugeField &out, const GaugeField &in) { + + FeynmanGaugeMomentumSpacePropagator_TL(out,in); + + GridBase *grid = out._grid; + LatticeInteger coor(grid); + GaugeField zz(grid); zz=zero; + + // xyzt + for(int d = 0; d < grid->_ndimension-1;d++){ + LatticeCoordinate(coor,d); + out = where(coor==Integer(0),zz,out); + } + } + void FeynmanGaugeMomentumSpacePropagator_TL(GaugeField &out, const GaugeField &in) { + + + // what type LatticeComplex + GridBase *grid = out._grid; + int nd = grid->_ndimension; + + typedef typename GaugeField::vector_type vector_type; + typedef typename GaugeField::scalar_type ScalComplex; + typedef Lattice > LatComplex; + + std::vector latt_size = grid->_fdimensions; + + LatComplex denom(grid); denom= zero; + LatComplex one(grid); one = ScalComplex(1.0,0.0); + + LatComplex kmu(grid); + ScalComplex ci(0.0,1.0); + // momphase = n * 2pi / L + for(int mu=0;mu zero_mode(nd,0); + TComplexD Tone = ComplexD(1.0,0.0); + TComplexD Tzero= ComplexD(0.0,0.0); + + pokeSite(Tone,denom,zero_mode); + + denom= one/denom; + + pokeSite(Tzero,denom,zero_mode); + + out = zero; + out = in*denom; + }; + + }; + +}} +#endif From 34f887ca1c35234717568ee80c9edbab6df7eef4 Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 26 Sep 2016 09:44:36 +0100 Subject: [PATCH 20/58] Test_fft not complete; preparing for tests of momentum space DWF and Overlap feynman rules but not there yet. --- tests/core/Test_fft.cc | 98 ++++++++++++++++++++++++++++++++++++------ 1 file changed, 84 insertions(+), 14 deletions(-) diff --git a/tests/core/Test_fft.cc b/tests/core/Test_fft.cc index c076695d..a20791fb 100644 --- a/tests/core/Test_fft.cc +++ b/tests/core/Test_fft.cc @@ -27,6 +27,7 @@ Author: Peter Boyle *************************************************************************************/ /* END LEGAL */ #include +#include using namespace Grid; using namespace Grid::QCD; @@ -127,27 +128,96 @@ int main (int argc, char ** argv) pRNG.SeedFixedIntegers(seeds); LatticeGaugeFieldD Umu(&GRID); - LatticeFermionD src(&GRID); gaussian(pRNG,src); - LatticeFermionD tmp(&GRID); - LatticeFermionD ref(&GRID); SU3::ColdConfiguration(pRNG,Umu); // Unit gauge - RealD mass=0.1; - WilsonFermionD Dw(Umu,GRID,RBGRID,mass); + { + LatticeFermionD src(&GRID); gaussian(pRNG,src); + LatticeFermionD tmp(&GRID); + LatticeFermionD ref(&GRID); + + RealD mass=0.1; + WilsonFermionD Dw(Umu,GRID,RBGRID,mass); + + Dw.M(src,tmp); - Dw.M(src,tmp); + std::cout << " src = " < HermOp(Ddwf); + ConjugateGradient CG(1.0e-4,1000); + CG(HermOp,src5,result5); + result5 = zero; + + ExtractSlice(tmp,result5,0,sdir); + result4 = (tmp+G5*tmp)*0.5; + + ExtractSlice(tmp,result5,Ls-1,sdir); + result4 = result4+(tmp-G5*tmp)*0.5; + + std::cout << "src "< QEDGimplTypesD; + typedef Photon QEDGaction; + QEDGaction Maxwell(QEDGaction::FEYNMAN_L); + QEDGaction::GaugeField Prop(&GRID);Prop=zero; + QEDGaction::GaugeField Source(&GRID);Source=zero; + + Maxwell.FreePropagator (Source,Prop); + std::cout << " MaxwellFree propagator\n"; + } Grid_finalize(); } From 167cc2650ea8ad87f9920326513feea559ffc523 Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 26 Sep 2016 09:58:09 +0100 Subject: [PATCH 21/58] GNU SOURCE problem on travis --- lib/Init.cc | 2 -- 1 file changed, 2 deletions(-) diff --git a/lib/Init.cc b/lib/Init.cc index cff43185..2f9619a9 100644 --- a/lib/Init.cc +++ b/lib/Init.cc @@ -331,8 +331,6 @@ void Grid_sa_signal_handler(int sig,siginfo_t *si,void * ptr) return; }; - -#define _GNU_SOURCE #include #ifdef __APPLE__ From 9353b6edfe7fd0567f40fad0f0dfb497a5d91c22 Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 26 Sep 2016 10:09:13 +0100 Subject: [PATCH 22/58] Fenv out of grid namespace --- lib/Init.cc | 41 ++++++++++++++++++++--------------------- 1 file changed, 20 insertions(+), 21 deletions(-) diff --git a/lib/Init.cc b/lib/Init.cc index 2f9619a9..221c056d 100644 --- a/lib/Init.cc +++ b/lib/Init.cc @@ -48,6 +48,26 @@ Author: paboyle #include +#include +#ifdef __APPLE__ +static int +feenableexcept (unsigned int excepts) +{ + static fenv_t fenv; + unsigned int new_excepts = excepts & FE_ALL_EXCEPT, + old_excepts; // previous masks + + if ( fegetenv (&fenv) ) return -1; + old_excepts = fenv.__control & FE_ALL_EXCEPT; + + // unmask + fenv.__control &= ~new_excepts; + fenv.__mxcsr &= ~(new_excepts << 7); + + return ( fesetenv (&fenv) ? -1 : old_excepts ); +} +#endif + namespace Grid { @@ -331,27 +351,6 @@ void Grid_sa_signal_handler(int sig,siginfo_t *si,void * ptr) return; }; -#include - -#ifdef __APPLE__ -static int -feenableexcept (unsigned int excepts) -{ - static fenv_t fenv; - unsigned int new_excepts = excepts & FE_ALL_EXCEPT, - old_excepts; // previous masks - - if ( fegetenv (&fenv) ) return -1; - old_excepts = fenv.__control & FE_ALL_EXCEPT; - - // unmask - fenv.__control &= ~new_excepts; - fenv.__mxcsr &= ~(new_excepts << 7); - - return ( fesetenv (&fenv) ? -1 : old_excepts ); -} -#endif - void Grid_debug_handler_init(void) { struct sigaction sa,osa; From 87acd06990090acae4f19ae147ca794d243764f5 Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 26 Sep 2016 10:11:34 +0100 Subject: [PATCH 23/58] Use streaming stores --- lib/simd/Intel512common.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/simd/Intel512common.h b/lib/simd/Intel512common.h index dabbf6d8..74793638 100644 --- a/lib/simd/Intel512common.h +++ b/lib/simd/Intel512common.h @@ -139,8 +139,8 @@ Author: paboyle #define ZLOADd(OFF,PTR,ri,ir) VLOADd(OFF,PTR,ir) VSHUFd(ir,ri) -#define VSTOREf(OFF,PTR,SRC) "vmovaps " #SRC "," #OFF "*64(" #PTR ")" ";\n" -#define VSTOREd(OFF,PTR,SRC) "vmovapd " #SRC "," #OFF "*64(" #PTR ")" ";\n" +#define VSTOREf(OFF,PTR,SRC) "vmovntps " #SRC "," #OFF "*64(" #PTR ")" ";\n" +#define VSTOREd(OFF,PTR,SRC) "vmovntpd " #SRC "," #OFF "*64(" #PTR ")" ";\n" // Swaps Re/Im ; could unify this with IMCI #define VSHUFd(A,DEST) "vpshufd $0x4e," #A "," #DEST ";\n" From 082ae350c66a0f514bb5f4bd009a28c5ba9bbb35 Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 10 Oct 2016 23:42:30 +0100 Subject: [PATCH 24/58] static schedule by default --- lib/Threads.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/Threads.h b/lib/Threads.h index 9d1295e5..5abe2ca4 100644 --- a/lib/Threads.h +++ b/lib/Threads.h @@ -37,7 +37,7 @@ Author: paboyle #ifdef GRID_OMP #include -#define PARALLEL_FOR_LOOP _Pragma("omp parallel for ") +#define PARALLEL_FOR_LOOP _Pragma("omp parallel for schedule(static)") #define PARALLEL_NESTED_LOOP2 _Pragma("omp parallel for collapse(2)") #else #define PARALLEL_FOR_LOOP From 09ca32d678461ad1da759e286d7007d629584bd4 Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 10 Oct 2016 23:42:55 +0100 Subject: [PATCH 25/58] Dminus added for Cayley --- lib/qcd/action/fermion/CayleyFermion5D.cc | 24 +++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/lib/qcd/action/fermion/CayleyFermion5D.cc b/lib/qcd/action/fermion/CayleyFermion5D.cc index c4196043..57b047d4 100644 --- a/lib/qcd/action/fermion/CayleyFermion5D.cc +++ b/lib/qcd/action/fermion/CayleyFermion5D.cc @@ -50,6 +50,30 @@ namespace QCD { mass(_mass) { } +template +void CayleyFermion5D::Dminus(const FermionField &psi, FermionField &chi) +{ + int Ls=this->Ls; + FermionField tmp(psi._grid); + + this->DW(psi,tmp,DaggerNo); + + for(int s=0;s +void CayleyFermion5D::DminusDag(const FermionField &psi, FermionField &chi) +{ + int Ls=this->Ls; + FermionField tmp(psi._grid); + + this->DW(psi,tmp,DaggerYes); + + for(int s=0;s void CayleyFermion5D::M5D (const FermionField &psi, FermionField &chi) { From c0d5b99016d6ac9914c5b9eeee8a5fc6c760f54b Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 10 Oct 2016 23:43:19 +0100 Subject: [PATCH 26/58] Dminus --- lib/qcd/action/fermion/CayleyFermion5D.h | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/lib/qcd/action/fermion/CayleyFermion5D.h b/lib/qcd/action/fermion/CayleyFermion5D.h index 53a93a05..1d8c2b95 100644 --- a/lib/qcd/action/fermion/CayleyFermion5D.h +++ b/lib/qcd/action/fermion/CayleyFermion5D.h @@ -56,6 +56,9 @@ namespace Grid { virtual void M5D (const FermionField &psi, FermionField &chi); virtual void M5Ddag(const FermionField &psi, FermionField &chi); + virtual void Dminus(const FermionField &psi, FermionField &chi); + virtual void DminusDag(const FermionField &psi, FermionField &chi); + ///////////////////////////////////////////////////// // Instantiate different versions depending on Impl ///////////////////////////////////////////////////// @@ -117,6 +120,7 @@ namespace Grid { GridRedBlackCartesian &FourDimRedBlackGrid, RealD _mass,RealD _M5,const ImplParams &p= ImplParams()); + protected: void SetCoefficientsZolotarev(RealD zolohi,Approx::zolotarev_data *zdata,RealD b,RealD c); void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c); From d7ce164e6e25ef2fd2047583d5b74753b6a64769 Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 10 Oct 2016 23:43:36 +0100 Subject: [PATCH 27/58] Feynman rule for DWF --- lib/qcd/action/fermion/DomainWallFermion.h | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/lib/qcd/action/fermion/DomainWallFermion.h b/lib/qcd/action/fermion/DomainWallFermion.h index b72ca8ec..c0b6b6aa 100644 --- a/lib/qcd/action/fermion/DomainWallFermion.h +++ b/lib/qcd/action/fermion/DomainWallFermion.h @@ -42,6 +42,10 @@ namespace Grid { INHERIT_IMPL_TYPES(Impl); public: + void MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _m) { + this->MomentumSpacePropagatorHt(out,in,_m); + }; + virtual void Instantiatable(void) {}; // Constructors DomainWallFermion(GaugeField &_Umu, @@ -51,6 +55,7 @@ namespace Grid { GridRedBlackCartesian &FourDimRedBlackGrid, RealD _mass,RealD _M5,const ImplParams &p= ImplParams()) : + CayleyFermion5D(_Umu, FiveDimGrid, FiveDimRedBlackGrid, From c0145745047b2500f19b57cfa4c28d2aa7ac3ddc Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 10 Oct 2016 23:44:00 +0100 Subject: [PATCH 28/58] A "please implement me" feynman rule. If this were abstract virtual it would require/force implementation --- lib/qcd/action/fermion/FermionOperator.h | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h index daee2e5b..742c6e08 100644 --- a/lib/qcd/action/fermion/FermionOperator.h +++ b/lib/qcd/action/fermion/FermionOperator.h @@ -92,15 +92,16 @@ namespace Grid { virtual void Mdir (const FermionField &in, FermionField &out,int dir,int disp)=0; // case by case Wilson, Clover, Cayley, ContFrac, PartFrac - virtual void MomentumSpacePropagator(FermionField &out,const FermionField &in) { assert(0);}; - virtual void FreePropagator(const FermionField &in,FermionField &out) { + virtual void MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _m) { assert(0);}; + + virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass) { FFT theFFT((GridCartesian *) in._grid); FermionField in_k(in._grid); FermionField prop_k(in._grid); theFFT.FFT_all_dim(in_k,in,FFT::forward); - this->MomentumSpacePropagator(prop_k,in_k); + this->MomentumSpacePropagator(prop_k,in_k,mass); theFFT.FFT_all_dim(out,prop_k,FFT::backward); }; From 6f26d2e8d4133ce88a1ca5056b083ecdf7d545b6 Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 10 Oct 2016 23:45:18 +0100 Subject: [PATCH 29/58] Overlap tree level feynman rule --- lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h b/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h index cd23ddd5..9cab0e22 100644 --- a/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h +++ b/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h @@ -42,7 +42,11 @@ namespace Grid { INHERIT_IMPL_TYPES(Impl); public: - // Constructors + void MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _m) { + this->MomentumSpacePropagatorHw(out,in,_m); + }; + + // Constructors OverlapWilsonCayleyTanhFermion(GaugeField &_Umu, GridCartesian &FiveDimGrid, GridRedBlackCartesian &FiveDimRedBlackGrid, From 616e7cd83eb4c45b3102c6a232c5f26b61517a01 Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 10 Oct 2016 23:45:48 +0100 Subject: [PATCH 30/58] Mass parameter --- lib/qcd/action/fermion/WilsonFermion.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index 55fbfa01..36c1870c 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -137,7 +137,7 @@ namespace QCD { MooeeInv(in,out); } template - void WilsonFermion:: MomentumSpacePropagator(FermionField &out, const FermionField &in) { + void WilsonFermion:: MomentumSpacePropagator(FermionField &out, const FermionField &in,RealD _m) { // what type LatticeComplex conformable(_grid,out._grid); @@ -179,7 +179,7 @@ namespace QCD { denom=denom + sin(kmu)*sin(kmu); } - wilson = wilson + mass; // 2 sin^2 k/2 + m + wilson = wilson + _m; // 2 sin^2 k/2 + m num = num + wilson*in; // -i gmu sin k + 2 sin^2 k/2 + m From 657e0a8f4d457f0350826e2a12eff6643f031ccd Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 10 Oct 2016 23:46:10 +0100 Subject: [PATCH 31/58] Mass parameter --- lib/qcd/action/fermion/WilsonFermion.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/qcd/action/fermion/WilsonFermion.h b/lib/qcd/action/fermion/WilsonFermion.h index c050c0eb..3a1ba432 100644 --- a/lib/qcd/action/fermion/WilsonFermion.h +++ b/lib/qcd/action/fermion/WilsonFermion.h @@ -78,7 +78,7 @@ namespace Grid { virtual void MooeeInv(const FermionField &in, FermionField &out) ; virtual void MooeeInvDag(const FermionField &in, FermionField &out) ; - virtual void MomentumSpacePropagator(FermionField &out,const FermionField &in) ; + virtual void MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _mass) ; //////////////////////// // Derivative interface From 96f1d1b828f613a7d2278a441cc86ee59cfd2526 Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 10 Oct 2016 23:46:45 +0100 Subject: [PATCH 32/58] Debugged Domain wall and Overlap feynman rules (infinite Ls, finite mass). --- lib/qcd/action/fermion/WilsonFermion5D.cc | 184 ++++++++++------------ 1 file changed, 81 insertions(+), 103 deletions(-) diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 28f02b67..ea3e8535 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -393,92 +393,95 @@ void WilsonFermion5D::DW(const FermionField &in, FermionField &out,int dag } template -void WilsonFermion5D::MomentumSpacePropagatorHt(FermionField &out,const FermionField &in) +void WilsonFermion5D::MomentumSpacePropagatorHt(FermionField &out,const FermionField &in, RealD mass) { - // what type LatticeComplex - GridBase *_grid = _FourDimGrid; - conformable(_grid,out._grid); + // what type LatticeComplex + GridBase *_grid = _FourDimGrid; + conformable(_grid,out._grid); + + typedef typename FermionField::vector_type vector_type; + typedef typename FermionField::scalar_type ScalComplex; + typedef iSinglet Tcomplex; + typedef Lattice > LatComplex; + + Gamma::GammaMatrix Gmu [] = { + Gamma::GammaX, + Gamma::GammaY, + Gamma::GammaZ, + Gamma::GammaT + }; - typedef typename FermionField::vector_type vector_type; - typedef typename FermionField::scalar_type ScalComplex; - typedef iSinglet Tcomplex; - typedef Lattice > LatComplex; + std::vector latt_size = _grid->_fdimensions; - Gamma::GammaMatrix Gmu [] = { - Gamma::GammaX, - Gamma::GammaY, - Gamma::GammaZ, - Gamma::GammaT - }; + + FermionField num (_grid); num = zero; - std::vector latt_size = _grid->_fdimensions; + LatComplex sk(_grid); sk = zero; + LatComplex sk2(_grid); sk2= zero; + LatComplex W(_grid); W= zero; + LatComplex a(_grid); a= zero; + LatComplex one (_grid); one = ScalComplex(1.0,0.0); + LatComplex denom(_grid); denom= zero; + LatComplex cosha(_grid); + LatComplex kmu(_grid); + LatComplex Wea(_grid); + LatComplex Wema(_grid); - - LatComplex sk(_grid); sk = zero; - LatComplex sk2(_grid); sk2= zero; - - LatComplex W(_grid); W= zero; - LatComplex a(_grid); a= zero; - LatComplex cosha(_grid); a= zero; - - LatComplex one (_grid); one = ScalComplex(1.0,0.0); - - FermionField num (_grid); num = zero; - LatComplex denom(_grid); denom= zero; - LatComplex kmu(_grid); - ScalComplex ci(0.0,1.0); - - for(int mu=0;mulSites();idx++){ - Tcomplex cc; - RealD sgn; - std::vector lcoor(Nd); - _grid->LocalIndexToLocalCoor(idx,lcoor); - peekLocalSite(cc,cosha,lcoor); - sgn= ::fabs(real(cc)); - cc = ScalComplex(::acosh(sgn),0.0); - pokeLocalSite(cc,a,lcoor); - } + LatticeCoordinate(kmu,mu); + + RealD TwoPiL = M_PI * 2.0/ latt_size[mu]; + + kmu = TwoPiL * kmu; + + sk2 = sk2 + 2.0*sin(kmu*0.5)*sin(kmu*0.5); + sk = sk + sin(kmu) *sin(kmu); + + num = num - sin(kmu)*ci*(Gamma(Gmu[mu])*in); + + } + + W = one - M5 + sk2; - std::cout << "Ht mom space num" << norm2(num)< alpha + //////////////////////////////////////////// + cosha = (one + W*W + sk) / (W*2.0); + + // FIXME Need a Lattice acosh + for(int idx=0;idx<_grid->lSites();idx++){ + std::vector lcoor(Nd); + Tcomplex cc; + RealD sgn; + _grid->LocalIndexToLocalCoor(idx,lcoor); + peekLocalSite(cc,cosha,lcoor); + assert((double)real(cc)>=1.0); + assert(fabs((double)imag(cc))<=1.0e-15); + cc = ScalComplex(::acosh(real(cc)),0.0); + pokeLocalSite(cc,a,lcoor); + } + + Wea = ( exp( a) * W ); + Wema= ( exp(-a) * W ); + + num = num + ( one - Wema ) * mass * in; + denom= ( Wea - one ) + mass*mass * (one - Wema); + out = num/denom; } template -void WilsonFermion5D::MomentumSpacePropagatorHw(FermionField &out,const FermionField &in) +void WilsonFermion5D::MomentumSpacePropagatorHw(FermionField &out,const FermionField &in,RealD mass) { - // what type LatticeComplex + Gamma::GammaMatrix Gmu [] = { + Gamma::GammaX, + Gamma::GammaY, + Gamma::GammaZ, + Gamma::GammaT + }; + GridBase *_grid = _FourDimGrid; conformable(_grid,out._grid); @@ -487,16 +490,9 @@ void WilsonFermion5D::MomentumSpacePropagatorHw(FermionField &out,const Fe typedef Lattice > LatComplex; - Gamma::GammaMatrix Gmu [] = { - Gamma::GammaX, - Gamma::GammaY, - Gamma::GammaZ, - Gamma::GammaT - }; std::vector latt_size = _grid->_fdimensions; - LatComplex sk(_grid); sk = zero; LatComplex sk2(_grid); sk2= zero; @@ -524,35 +520,17 @@ void WilsonFermion5D::MomentumSpacePropagatorHw(FermionField &out,const Fe num = num - sin(kmu)*ci*(Gamma(Gmu[mu])*in); } + num = num + mass * in ; b_k = sk2 - M5; w_k = sqrt(sk + b_k*b_k); - std::cout << "Hw M5 " << M5< Date: Mon, 10 Oct 2016 23:47:33 +0100 Subject: [PATCH 33/58] Mass parameter --- lib/qcd/action/fermion/WilsonFermion5D.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/qcd/action/fermion/WilsonFermion5D.h b/lib/qcd/action/fermion/WilsonFermion5D.h index a56bee11..fc8b5d41 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.h +++ b/lib/qcd/action/fermion/WilsonFermion5D.h @@ -88,8 +88,8 @@ namespace Grid { virtual void DhopDerivEO(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); virtual void DhopDerivOE(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); - void MomentumSpacePropagatorHt(FermionField &out,const FermionField &in) ; - void MomentumSpacePropagatorHw(FermionField &out,const FermionField &in) ; + void MomentumSpacePropagatorHt(FermionField &out,const FermionField &in,RealD mass) ; + void MomentumSpacePropagatorHw(FermionField &out,const FermionField &in,RealD mass) ; // Implement hopping term non-hermitian hopping term; half cb or both // Implement s-diagonal DW From dc389e467c833c13ed1a4a882009351297ff1b12 Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 10 Oct 2016 23:48:05 +0100 Subject: [PATCH 34/58] axpy_ssp for any coeff type via template --- lib/qcd/utils/LinalgUtils.h | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/lib/qcd/utils/LinalgUtils.h b/lib/qcd/utils/LinalgUtils.h index c2cb95ee..e7e6a794 100644 --- a/lib/qcd/utils/LinalgUtils.h +++ b/lib/qcd/utils/LinalgUtils.h @@ -39,8 +39,8 @@ namespace QCD{ //on the 5d (rb4d) checkerboarded lattices //////////////////////////////////////////////////////////////////////// -template -void axpibg5x(Lattice &z,const Lattice &x,RealD a,RealD b) +template +void axpibg5x(Lattice &z,const Lattice &x,Coeff a,Coeff b) { z.checkerboard = x.checkerboard; conformable(x,z); @@ -57,8 +57,8 @@ PARALLEL_FOR_LOOP } } -template -void axpby_ssp(Lattice &z, RealD a,const Lattice &x,RealD b,const Lattice &y,int s,int sp) +template +void axpby_ssp(Lattice &z, Coeff a,const Lattice &x,Coeff b,const Lattice &y,int s,int sp) { z.checkerboard = x.checkerboard; conformable(x,y); @@ -72,8 +72,8 @@ PARALLEL_FOR_LOOP } } -template -void ag5xpby_ssp(Lattice &z,RealD a,const Lattice &x,RealD b,const Lattice &y,int s,int sp) +template +void ag5xpby_ssp(Lattice &z,Coeff a,const Lattice &x,Coeff b,const Lattice &y,int s,int sp) { z.checkerboard = x.checkerboard; conformable(x,y); @@ -90,8 +90,8 @@ PARALLEL_FOR_LOOP } } -template -void axpbg5y_ssp(Lattice &z,RealD a,const Lattice &x,RealD b,const Lattice &y,int s,int sp) +template +void axpbg5y_ssp(Lattice &z,Coeff a,const Lattice &x,Coeff b,const Lattice &y,int s,int sp) { z.checkerboard = x.checkerboard; conformable(x,y); @@ -108,8 +108,8 @@ PARALLEL_FOR_LOOP } } -template -void ag5xpbg5y_ssp(Lattice &z,RealD a,const Lattice &x,RealD b,const Lattice &y,int s,int sp) +template +void ag5xpbg5y_ssp(Lattice &z,Coeff a,const Lattice &x,Coeff b,const Lattice &y,int s,int sp) { z.checkerboard = x.checkerboard; conformable(x,y); @@ -127,8 +127,8 @@ PARALLEL_FOR_LOOP } } -template -void axpby_ssp_pminus(Lattice &z,RealD a,const Lattice &x,RealD b,const Lattice &y,int s,int sp) +template +void axpby_ssp_pminus(Lattice &z,Coeff a,const Lattice &x,Coeff b,const Lattice &y,int s,int sp) { z.checkerboard = x.checkerboard; conformable(x,y); @@ -144,8 +144,8 @@ PARALLEL_FOR_LOOP } } -template -void axpby_ssp_pplus(Lattice &z,RealD a,const Lattice &x,RealD b,const Lattice &y,int s,int sp) +template +void axpby_ssp_pplus(Lattice &z,Coeff a,const Lattice &x,Coeff b,const Lattice &y,int s,int sp) { z.checkerboard = x.checkerboard; conformable(x,y); From db749f103f86fbf3977cccd66b7b06f3d9687a9a Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 10 Oct 2016 23:48:35 +0100 Subject: [PATCH 35/58] Add Wilson, DWF, Overlap feynman rule tests --- tests/core/Test_fft.cc | 346 ++++++++++++++++++++++++++++++++++------- 1 file changed, 289 insertions(+), 57 deletions(-) diff --git a/tests/core/Test_fft.cc b/tests/core/Test_fft.cc index a20791fb..562a53a1 100644 --- a/tests/core/Test_fft.cc +++ b/tests/core/Test_fft.cc @@ -1,6 +1,6 @@ /************************************************************************************* - Grid physics library, www.github.com/paboyle/Grid + grid` physics library, www.github.com/paboyle/Grid Source file: ./tests/Test_cshift.cc @@ -61,18 +61,22 @@ int main (int argc, char ** argv) LatticeSpinMatrixD S(&GRID); LatticeSpinMatrixD Stilde(&GRID); - std::vector p({1,2,3,2}); + std::vector p({1,3,2,3}); one = ComplexD(1.0,0.0); zz = ComplexD(0.0,0.0); ComplexD ci(0.0,1.0); + + std::cout<<"*************************************************"< seeds({1,2,3,4}); + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds); // naughty seeding GridParallelRNG pRNG(&GRID); pRNG.SeedFixedIntegers(seeds); LatticeGaugeFieldD Umu(&GRID); SU3::ColdConfiguration(pRNG,Umu); // Unit gauge - + // Umu=zero; + //////////////////////////////////////////////////// + // Wilson test + //////////////////////////////////////////////////// { LatticeFermionD src(&GRID); gaussian(pRNG,src); LatticeFermionD tmp(&GRID); LatticeFermionD ref(&GRID); - RealD mass=0.1; + RealD mass=0.01; WilsonFermionD Dw(Umu,GRID,RBGRID,mass); Dw.M(src,tmp); - std::cout << " src = " < HermOp(Ddwf); - ConjugateGradient CG(1.0e-4,1000); - CG(HermOp,src5,result5); - result5 = zero; - - ExtractSlice(tmp,result5,0,sdir); - result4 = (tmp+G5*tmp)*0.5; + RealD TwoPiL = M_PI * 2.0/ latt_size[mu]; + + kmu = TwoPiL * kmu; + + sk2 = sk2 + 2.0*sin(kmu*0.5)*sin(kmu*0.5); + sk = sk + sin(kmu) *sin(kmu); + + // -1/2 Dw -> 1/2 gmu (eip - emip) = i sinp gmu + Kinetic = Kinetic + sin(kmu)*ci*(Gamma(Gmu[mu])*src5_p); - ExtractSlice(tmp,result5,Ls-1,sdir); - result4 = result4+(tmp-G5*tmp)*0.5; - - std::cout << "src "< point(4,0); + src=zero; + SpinColourVectorD ferm; gaussian(sRNG,ferm); + pokeSite(ferm,src,point); + + const int Ls=32; + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,&GRID); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,&GRID); + + RealD mass=0.01; + RealD M5 =0.8; + DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,GRID,RBGRID,mass,M5); + + // Momentum space prop + std::cout << " Solving by FFT and Feynman rules" < HermOp(Ddwf); + ConjugateGradient CG(1.0e-16,10000); + CG(HermOp,src5,result5); + + //////////////////////////////////////////////////////////////////////// + // Domain wall physical field propagator + //////////////////////////////////////////////////////////////////////// + /* + psi = chiralProjectMinus(psi_5[0]); + psi += chiralProjectPlus(psi_5[Ls-1]); + */ + ExtractSlice(tmp,result5,0 ,sdir); result4 = (tmp-G5*tmp)*0.5; + ExtractSlice(tmp,result5,Ls-1,sdir); result4 = result4+(tmp+G5*tmp)*0.5; + + std::cout << " Taking difference" < point(4,0); + src=zero; + SpinColourVectorD ferm; gaussian(sRNG,ferm); + pokeSite(ferm,src,point); + + const int Ls=48; + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,&GRID); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,&GRID); + + RealD mass=0.01; + RealD M5 =0.8; + + OverlapWilsonCayleyTanhFermionD Dov(Umu,*FGrid,*FrbGrid,GRID,RBGRID,mass,M5,1.0); + + // Momentum space prop + std::cout << " Solving by FFT and Feynman rules" < HermOp(Dov); + ConjugateGradient CG(1.0e-16,10000); + CG(HermOp,src5,result5); + + //////////////////////////////////////////////////////////////////////// + // Domain wall physical field propagator + //////////////////////////////////////////////////////////////////////// + /* + psi = chiralProjectMinus(psi_5[0]); + psi += chiralProjectPlus(psi_5[Ls-1]); + */ + ExtractSlice(tmp,result5,0 ,sdir); result4 = (tmp-G5*tmp)*0.5; + ExtractSlice(tmp,result5,Ls-1,sdir); result4 = result4+(tmp+G5*tmp)*0.5; + + std::cout << " Taking difference" < QEDGimplTypesD; typedef Photon QEDGaction; + QEDGaction Maxwell(QEDGaction::FEYNMAN_L); QEDGaction::GaugeField Prop(&GRID);Prop=zero; QEDGaction::GaugeField Source(&GRID);Source=zero; From 3d5c9a1ee98146322cb48eb6674b47f1fcfc6063 Mon Sep 17 00:00:00 2001 From: paboyle Date: Mon, 10 Oct 2016 23:50:13 +0100 Subject: [PATCH 36/58] No compile fix on clang++ 3.9 --- lib/simd/Grid_avx.h | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/lib/simd/Grid_avx.h b/lib/simd/Grid_avx.h index 33496e07..902142d2 100644 --- a/lib/simd/Grid_avx.h +++ b/lib/simd/Grid_avx.h @@ -450,13 +450,12 @@ namespace Optimization { }; #if defined (AVX2) || defined (AVXFMA4) -#define _mm256_alignr_epi32(ret,a,b,n) ret=(__m256) _mm256_alignr_epi8((__m256i)a,(__m256i)b,(n*4)%16) -#define _mm256_alignr_epi64(ret,a,b,n) ret=(__m256d) _mm256_alignr_epi8((__m256i)a,(__m256i)b,(n*8)%16) +#define _mm256_alignr_epi32_grid(ret,a,b,n) ret=(__m256) _mm256_alignr_epi8((__m256i)a,(__m256i)b,(n*4)%16) +#define _mm256_alignr_epi64_grid(ret,a,b,n) ret=(__m256d) _mm256_alignr_epi8((__m256i)a,(__m256i)b,(n*8)%16) #endif #if defined (AVX1) - -#define _mm256_alignr_epi32(ret,a,b,n) { \ +#define _mm256_alignr_epi32_grid(ret,a,b,n) { \ __m128 aa, bb; \ \ aa = _mm256_extractf128_ps(a,1); \ @@ -470,7 +469,7 @@ namespace Optimization { ret = _mm256_insertf128_ps(ret,aa,0); \ } -#define _mm256_alignr_epi64(ret,a,b,n) { \ +#define _mm256_alignr_epi64_grid(ret,a,b,n) { \ __m128d aa, bb; \ \ aa = _mm256_extractf128_pd(a,1); \ @@ -530,9 +529,9 @@ namespace Optimization { __m256 tmp = Permute::Permute0(in); __m256 ret; if ( n > 3 ) { - _mm256_alignr_epi32(ret,in,tmp,n); + _mm256_alignr_epi32_grid(ret,in,tmp,n); } else { - _mm256_alignr_epi32(ret,tmp,in,n); + _mm256_alignr_epi32_grid(ret,tmp,in,n); } // std::cout << " align epi32 n=" < "<< ret < 1 ) { - _mm256_alignr_epi64(ret,in,tmp,n); + _mm256_alignr_epi64_grid(ret,in,tmp,n); } else { - _mm256_alignr_epi64(ret,tmp,in,n); + _mm256_alignr_epi64_grid(ret,tmp,in,n); } // std::cout << " align epi64 n=" < "<< ret < Date: Mon, 10 Oct 2016 23:50:42 +0100 Subject: [PATCH 37/58] verbose remove --- lib/simd/Grid_vector_types.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/lib/simd/Grid_vector_types.h b/lib/simd/Grid_vector_types.h index 4c3e903d..e1592eef 100644 --- a/lib/simd/Grid_vector_types.h +++ b/lib/simd/Grid_vector_types.h @@ -583,8 +583,6 @@ inline Grid_simd operator/(Grid_simd a, Grid_simd b) { ret.v=binary(ret.v, real_den.v, DivSIMD()); - std::cout << " Div call from complex "< "<< ret < "<< ret < inline Grid::ComplexF Reduce::operator()(__m256 in){ From 6e01264bb7f5bafa2bc5b2043a93d77d8393f378 Mon Sep 17 00:00:00 2001 From: paboyle Date: Tue, 11 Oct 2016 10:03:39 +0100 Subject: [PATCH 39/58] don't use static by default --- lib/Threads.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/Threads.h b/lib/Threads.h index 5abe2ca4..9d1295e5 100644 --- a/lib/Threads.h +++ b/lib/Threads.h @@ -37,7 +37,7 @@ Author: paboyle #ifdef GRID_OMP #include -#define PARALLEL_FOR_LOOP _Pragma("omp parallel for schedule(static)") +#define PARALLEL_FOR_LOOP _Pragma("omp parallel for ") #define PARALLEL_NESTED_LOOP2 _Pragma("omp parallel for collapse(2)") #else #define PARALLEL_FOR_LOOP From f7c2aa3ba54db5edec66a7cc3af521dddd5b8dc9 Mon Sep 17 00:00:00 2001 From: paboyle Date: Wed, 12 Oct 2016 00:29:13 +0100 Subject: [PATCH 40/58] runtime by default --- lib/Threads.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/Threads.h b/lib/Threads.h index 9d1295e5..b8915d74 100644 --- a/lib/Threads.h +++ b/lib/Threads.h @@ -37,7 +37,7 @@ Author: paboyle #ifdef GRID_OMP #include -#define PARALLEL_FOR_LOOP _Pragma("omp parallel for ") +#define PARALLEL_FOR_LOOP _Pragma("omp parallel for schedule(runtime)") #define PARALLEL_NESTED_LOOP2 _Pragma("omp parallel for collapse(2)") #else #define PARALLEL_FOR_LOOP From 6b27c42dfe1d5c010e83804843a9b89bbf5480fd Mon Sep 17 00:00:00 2001 From: paboyle Date: Wed, 12 Oct 2016 00:29:39 +0100 Subject: [PATCH 41/58] Cosmetic --- lib/qcd/action/gauge/Photon.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/qcd/action/gauge/Photon.h b/lib/qcd/action/gauge/Photon.h index a7556b19..29d40f8b 100644 --- a/lib/qcd/action/gauge/Photon.h +++ b/lib/qcd/action/gauge/Photon.h @@ -77,8 +77,8 @@ namespace Grid{ out = where(coor==Integer(0),zz,out); } } - void FeynmanGaugeMomentumSpacePropagator_TL(GaugeField &out, const GaugeField &in) { + void FeynmanGaugeMomentumSpacePropagator_TL(GaugeField &out, const GaugeField &in) { // what type LatticeComplex GridBase *grid = out._grid; @@ -92,8 +92,8 @@ namespace Grid{ LatComplex denom(grid); denom= zero; LatComplex one(grid); one = ScalComplex(1.0,0.0); + LatComplex kmu(grid); - LatComplex kmu(grid); ScalComplex ci(0.0,1.0); // momphase = n * 2pi / L for(int mu=0;mu Date: Wed, 12 Oct 2016 00:29:57 +0100 Subject: [PATCH 42/58] Static required for shmem. Reading same object twice requires csum reset --- lib/parallelIO/BinaryIO.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/parallelIO/BinaryIO.h b/lib/parallelIO/BinaryIO.h index 184209dc..acd96a39 100644 --- a/lib/parallelIO/BinaryIO.h +++ b/lib/parallelIO/BinaryIO.h @@ -457,7 +457,7 @@ class BinaryIO { // available (how short sighted is that?) ////////////////////////////////////////////////////////// Umu = zero; - static uint32_t csum=0; + static uint32_t csum; csum=0; fobj fileObj; static sobj siteObj; // Static to place in symmetric region for SHMEM From 63d219498ba4a036d74c8520e7455c94543cff49 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 21 Oct 2016 13:10:13 +0100 Subject: [PATCH 43/58] first (dirty) implementation of Feynman stoctachtic EM field --- lib/qcd/action/gauge/Photon.h | 344 ++++++++++++++++++++++------------ 1 file changed, 226 insertions(+), 118 deletions(-) diff --git a/lib/qcd/action/gauge/Photon.h b/lib/qcd/action/gauge/Photon.h index 29d40f8b..ea6f21b9 100644 --- a/lib/qcd/action/gauge/Photon.h +++ b/lib/qcd/action/gauge/Photon.h @@ -1,126 +1,234 @@ - /************************************************************************************* - - Grid physics library, www.github.com/paboyle/Grid - - Source file: ./lib/qcd/action/gauge/Photon.h - - Copyright (C) 2015 - -Author: Peter Boyle - - 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 */ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/qcd/action/gauge/Photon.h + + Copyright (C) 2015 + + Author: Peter Boyle + + 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 QCD_PHOTON_ACTION_H #define QCD_PHOTON_ACTION_H namespace Grid{ - namespace QCD{ +namespace QCD{ + + template + class Photon + { + public: + INHERIT_GIMPL_TYPES(Gimpl); + enum class Gauge {Feynman, Coulomb, Landau}; + enum class ZmScheme {QedL, QedTL}; + public: + Photon(Gauge gauge, ZmScheme zmScheme); + virtual ~Photon(void) = default; + void FreePropagator(const GaugeField &in, GaugeField &out); + void MomentumSpacePropagator(const GaugeField &in, GaugeField &out); + void StochasticField(GaugeField &out, GridParallelRNG &rng); + private: + void kHatSquared(GaugeLinkField &out); + void zmSub(GaugeLinkField &out); + private: + Gauge gauge_; + ZmScheme zmScheme_; + }; - template - class Photon { - - public: - - INHERIT_GIMPL_TYPES(Gimpl); - - enum PhotonType { FEYNMAN_L, FEYNMAN_TL }; - - PhotonType GaugeBC; - - Photon(PhotonType _GaugeBC) : GaugeBC(_GaugeBC){}; - - void FreePropagator (const GaugeField &in,GaugeField &out) { - FFT theFFT((GridCartesian *) in._grid); - - GaugeField in_k(in._grid); - GaugeField prop_k(in._grid); - - theFFT.FFT_all_dim(in_k,in,FFT::forward); - MomentumSpacePropagator(prop_k,in_k); - theFFT.FFT_all_dim(out,prop_k,FFT::backward); - } - void MomentumSpacePropagator(GaugeField &out,const GaugeField &in) { - if ( GaugeBC == FEYNMAN_L ) { - FeynmanGaugeMomentumSpacePropagator_L(out,in); - } else if ( GaugeBC == FEYNMAN_TL ) { - FeynmanGaugeMomentumSpacePropagator_TL(out,in); - } else { // No coulomb yet - assert(0); - } - } - void FeynmanGaugeMomentumSpacePropagator_L(GaugeField &out, const GaugeField &in) { - - FeynmanGaugeMomentumSpacePropagator_TL(out,in); - - GridBase *grid = out._grid; - LatticeInteger coor(grid); - GaugeField zz(grid); zz=zero; - - // xyzt - for(int d = 0; d < grid->_ndimension-1;d++){ - LatticeCoordinate(coor,d); - out = where(coor==Integer(0),zz,out); - } - } - - void FeynmanGaugeMomentumSpacePropagator_TL(GaugeField &out, const GaugeField &in) { - - // what type LatticeComplex - GridBase *grid = out._grid; - int nd = grid->_ndimension; - - typedef typename GaugeField::vector_type vector_type; - typedef typename GaugeField::scalar_type ScalComplex; - typedef Lattice > LatComplex; + template + Photon::Photon(Gauge gauge, ZmScheme zmScheme) + : gauge_(gauge), zmScheme_(zmScheme) + {} + + template + void Photon::FreePropagator (const GaugeField &in,GaugeField &out) + { + FFT theFFT(in._grid); - std::vector latt_size = grid->_fdimensions; - - LatComplex denom(grid); denom= zero; - LatComplex one(grid); one = ScalComplex(1.0,0.0); - LatComplex kmu(grid); - - ScalComplex ci(0.0,1.0); - // momphase = n * 2pi / L - for(int mu=0;mu zero_mode(nd,0); - TComplexD Tone = ComplexD(1.0,0.0); - TComplexD Tzero= ComplexD(0.0,0.0); - - pokeSite(Tone,denom,zero_mode); - - denom= one/denom; - - pokeSite(Tzero,denom,zero_mode); - - out = zero; - out = in*denom; - }; - - }; - + GaugeField in_k(in._grid); + GaugeField prop_k(in._grid); + + theFFT.FFT_all_dim(in_k,in,FFT::forward); + MomentumSpacePropagator(prop_k,in_k); + theFFT.FFT_all_dim(out,prop_k,FFT::backward); + } + + template + void Photon::kHatSquared(GaugeLinkField &out) + { + GridBase *grid = out._grid; + GaugeLinkField kmu(grid); + const unsigned int nd = grid->_ndimension; + std::vector &l = grid->_fdimensions; + + out = zero; + for(int mu = 0; mu < nd; mu++) + { + RealD twoPiL = M_PI*2./l[mu]; + + LatticeCoordinate(kmu,mu); + kmu = 2.*sin(.5*twoPiL*kmu); + out = out + kmu*kmu; + } + } + + template + void Photon::zmSub(GaugeLinkField &out) + { + GridBase *grid = out._grid; + const unsigned int nd = grid->_ndimension; + + switch (zmScheme_) + { + case ZmScheme::QedTL: + { + std::vector zm(nd,0); + TComplexD Tzero = ComplexD(0.0,0.0); + + pokeSite(Tzero, out, zm); + + break; + } + case ZmScheme::QedL: + { + LatticeInteger spNrm(grid), coor(grid); + GaugeLinkField z(grid); + + spNrm = zero; + for(int d = 0; d < grid->_ndimension - 1; d++) + { + LatticeCoordinate(coor,d); + spNrm = spNrm + coor*coor; + } + out = where(spNrm == 0, 0.*out, out); + + break; + } + default: + break; + } + } + + template + void Photon::MomentumSpacePropagator(const GaugeField &in, + GaugeField &out) + { + GridBase *grid = out._grid; + const unsigned int nd = grid->_ndimension; + LatticeComplex k2Inv(grid), one(grid); + + one = ComplexD(1.0,0.0); + kHatSquared(k2Inv); + + std::vector zm(nd,0); + TComplexD Tone = ComplexD(1.0,0.0); + TComplexD Tzero= ComplexD(0.0,0.0); + + pokeSite(Tone, k2Inv, zm); + k2Inv = one/k2Inv; + pokeSite(Tzero, k2Inv,zm); + + out = in*k2Inv; + } + + template + void Photon::StochasticField(GaugeField &out, GridParallelRNG &rng) + { + auto *grid = out._grid; + const unsigned int nd = grid->_ndimension; + GaugeLinkField k2(grid), r(grid); + GaugeField aTilde(grid); + FFT fft(grid); + + kHatSquared(k2); + zmSub(k2); + for(int mu = 0; mu < nd; mu++) + { + gaussian(rng, r); + r = k2*r; + pokeLorentz(aTilde, r, mu); + } + fft.FFT_all_dim(out, aTilde, FFT::backward); + } +// template +// void Photon::FeynmanGaugeMomentumSpacePropagator_L(GaugeField &out, +// const GaugeField &in) +// { +// +// FeynmanGaugeMomentumSpacePropagator_TL(out,in); +// +// GridBase *grid = out._grid; +// LatticeInteger coor(grid); +// GaugeField zz(grid); zz=zero; +// +// // xyzt +// for(int d = 0; d < grid->_ndimension-1;d++){ +// LatticeCoordinate(coor,d); +// out = where(coor==Integer(0),zz,out); +// } +// } +// +// template +// void Photon::FeynmanGaugeMomentumSpacePropagator_TL(GaugeField &out, +// const GaugeField &in) +// { +// +// // what type LatticeComplex +// GridBase *grid = out._grid; +// int nd = grid->_ndimension; +// +// typedef typename GaugeField::vector_type vector_type; +// typedef typename GaugeField::scalar_type ScalComplex; +// typedef Lattice > LatComplex; +// +// std::vector latt_size = grid->_fdimensions; +// +// LatComplex denom(grid); denom= zero; +// LatComplex one(grid); one = ScalComplex(1.0,0.0); +// LatComplex kmu(grid); +// +// ScalComplex ci(0.0,1.0); +// // momphase = n * 2pi / L +// for(int mu=0;mu zero_mode(nd,0); +// TComplexD Tone = ComplexD(1.0,0.0); +// TComplexD Tzero= ComplexD(0.0,0.0); +// +// pokeSite(Tone,denom,zero_mode); +// +// denom= one/denom; +// +// pokeSite(Tzero,denom,zero_mode); +// +// out = zero; +// out = in*denom; +// }; + }} #endif From 462921e5491fbef99c93cfc460819890b0a63b66 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 21 Oct 2016 14:41:08 +0100 Subject: [PATCH 44/58] QED: fix stochastic field --- lib/qcd/action/gauge/Photon.h | 37 +++++++++++++++++------------------ 1 file changed, 18 insertions(+), 19 deletions(-) diff --git a/lib/qcd/action/gauge/Photon.h b/lib/qcd/action/gauge/Photon.h index ea6f21b9..aac72898 100644 --- a/lib/qcd/action/gauge/Photon.h +++ b/lib/qcd/action/gauge/Photon.h @@ -45,7 +45,7 @@ namespace QCD{ void MomentumSpacePropagator(const GaugeField &in, GaugeField &out); void StochasticField(GaugeField &out, GridParallelRNG &rng); private: - void kHatSquared(GaugeLinkField &out); + void invKHatSquared(GaugeLinkField &out); void zmSub(GaugeLinkField &out); private: Gauge gauge_; @@ -71,13 +71,17 @@ namespace QCD{ } template - void Photon::kHatSquared(GaugeLinkField &out) + void Photon::invKHatSquared(GaugeLinkField &out) { GridBase *grid = out._grid; - GaugeLinkField kmu(grid); + GaugeLinkField kmu(grid), one(grid); const unsigned int nd = grid->_ndimension; std::vector &l = grid->_fdimensions; + std::vector zm(nd,0); + TComplexD Tone = ComplexD(1.0,0.0); + TComplexD Tzero= ComplexD(0.0,0.0); + one = ComplexD(1.0,0.0); out = zero; for(int mu = 0; mu < nd; mu++) { @@ -87,6 +91,9 @@ namespace QCD{ kmu = 2.*sin(.5*twoPiL*kmu); out = out + kmu*kmu; } + pokeSite(Tone, out, zm); + out = one/out; + pokeSite(Tzero, out,zm); } template @@ -131,19 +138,10 @@ namespace QCD{ GaugeField &out) { GridBase *grid = out._grid; - const unsigned int nd = grid->_ndimension; - LatticeComplex k2Inv(grid), one(grid); + LatticeComplex k2Inv(grid); - one = ComplexD(1.0,0.0); - kHatSquared(k2Inv); - - std::vector zm(nd,0); - TComplexD Tone = ComplexD(1.0,0.0); - TComplexD Tzero= ComplexD(0.0,0.0); - - pokeSite(Tone, k2Inv, zm); - k2Inv = one/k2Inv; - pokeSite(Tzero, k2Inv,zm); + invKHatSquared(k2Inv); + zmSub(k2Inv); out = in*k2Inv; } @@ -153,16 +151,17 @@ namespace QCD{ { auto *grid = out._grid; const unsigned int nd = grid->_ndimension; - GaugeLinkField k2(grid), r(grid); + GaugeLinkField sqrtK2Inv(grid), r(grid); GaugeField aTilde(grid); FFT fft(grid); - kHatSquared(k2); - zmSub(k2); + invKHatSquared(sqrtK2Inv); + sqrtK2Inv = sqrt(real(sqrtK2Inv)); + zmSub(sqrtK2Inv); for(int mu = 0; mu < nd; mu++) { gaussian(rng, r); - r = k2*r; + r = sqrtK2Inv*r; pokeLorentz(aTilde, r, mu); } fft.FFT_all_dim(out, aTilde, FFT::backward); From 7c8f79b1479e9f15d1c82c27579f8741757a41f8 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 21 Oct 2016 15:20:12 +0100 Subject: [PATCH 45/58] more stochastic QED fixes --- lib/qcd/action/Actions.h | 1 + lib/qcd/action/gauge/Photon.h | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h index ba6e577d..a23c6495 100644 --- a/lib/qcd/action/Actions.h +++ b/lib/qcd/action/Actions.h @@ -57,6 +57,7 @@ Author: paboyle //////////////////////////////////////////// // Gauge Actions //////////////////////////////////////////// +#include #include #include diff --git a/lib/qcd/action/gauge/Photon.h b/lib/qcd/action/gauge/Photon.h index aac72898..0dd5a9dc 100644 --- a/lib/qcd/action/gauge/Photon.h +++ b/lib/qcd/action/gauge/Photon.h @@ -124,7 +124,7 @@ namespace QCD{ LatticeCoordinate(coor,d); spNrm = spNrm + coor*coor; } - out = where(spNrm == 0, 0.*out, out); + out = where(spNrm == Integer(0), 0.*out, out); break; } @@ -149,7 +149,7 @@ namespace QCD{ template void Photon::StochasticField(GaugeField &out, GridParallelRNG &rng) { - auto *grid = out._grid; + auto *grid = dynamic_cast(out._grid); const unsigned int nd = grid->_ndimension; GaugeLinkField sqrtK2Inv(grid), r(grid); GaugeField aTilde(grid); From bca861e112a0473f8215e66310489a9216b7870a Mon Sep 17 00:00:00 2001 From: Azusa Yamaguchi Date: Tue, 25 Oct 2016 14:21:48 +0100 Subject: [PATCH 46/58] Note:FFT shoud be GridFFT (Not change yet). Gauge fix with FFt is added (tests/core) --- lib/qcd/utils/SUn.h | 31 ++++ lib/qcd/utils/WilsonLoops.h | 2 +- tests/core/Test_fft.cc | 3 + tests/core/Test_fft_gfix.cc | 301 ++++++++++++++++++++++++++++++++++++ 4 files changed, 336 insertions(+), 1 deletion(-) create mode 100644 tests/core/Test_fft_gfix.cc diff --git a/lib/qcd/utils/SUn.h b/lib/qcd/utils/SUn.h index 3240098a..914c3bad 100644 --- a/lib/qcd/utils/SUn.h +++ b/lib/qcd/utils/SUn.h @@ -674,6 +674,37 @@ class SU { out += la; } } +/* + add GaugeTrans +*/ + +template + static void GaugeTransform( GaugeField &Umu, GaugeMat &g){ + GridBase *grid = Umu._grid; + conformable(grid,g._grid); + + GaugeMat U(grid); + GaugeMat ag(grid); ag = adj(g); + + for(int mu=0;mu(Umu,mu); + U = g*U*Cshift(ag, mu, 1); + PokeIndex(Umu,U,mu); + } + } + template + static void GaugeTransform( std::vector &U, GaugeMat &g){ + GridBase *grid = g._grid; + GaugeMat ag(grid); ag = adj(g); + for(int mu=0;mu + static void RandomGaugeTransform(GridParallelRNG &pRNG, GaugeField &Umu, GaugeMat &g){ + LieRandomize(pRNG,g,1.0); + GaugeTransform(Umu,g); + } // Projects the algebra components a lattice matrix (of dimension ncol*ncol -1 ) // inverse operation: FundamentalLieAlgebraMatrix diff --git a/lib/qcd/utils/WilsonLoops.h b/lib/qcd/utils/WilsonLoops.h index 10022f50..03d45c07 100644 --- a/lib/qcd/utils/WilsonLoops.h +++ b/lib/qcd/utils/WilsonLoops.h @@ -522,4 +522,4 @@ typedef WilsonLoops SU3WilsonLoops; } } -#endif \ No newline at end of file +#endif diff --git a/tests/core/Test_fft.cc b/tests/core/Test_fft.cc index 562a53a1..512d6e62 100644 --- a/tests/core/Test_fft.cc +++ b/tests/core/Test_fft.cc @@ -441,6 +441,8 @@ int main (int argc, char ** argv) } { + /* + * typedef GaugeImplTypes QEDGimplTypesD; typedef Photon QEDGaction; @@ -450,6 +452,7 @@ int main (int argc, char ** argv) Maxwell.FreePropagator (Source,Prop); std::cout << " MaxwellFree propagator\n"; + */ } Grid_finalize(); } diff --git a/tests/core/Test_fft_gfix.cc b/tests/core/Test_fft_gfix.cc new file mode 100644 index 00000000..1eda9a67 --- /dev/null +++ b/tests/core/Test_fft_gfix.cc @@ -0,0 +1,301 @@ + /************************************************************************************* + + grid` physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_cshift.cc + + Copyright (C) 2015 + +Author: Azusa Yamaguchi +Author: Peter Boyle + + 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 +#include + +using namespace Grid; +using namespace Grid::QCD; + +template +class FourierAcceleratedGaugeFixer : public Gimpl { + public: + INHERIT_GIMPL_TYPES(Gimpl); + + typedef typename Gimpl::GaugeLinkField GaugeMat; + typedef typename Gimpl::GaugeField GaugeLorentz; + + static void GaugeLinkToLieAlgebraField(const std::vector &U,std::vector &A) { + for(int mu=0;mu &A,GaugeMat &dmuAmu) { + dmuAmu=zero; + for(int mu=0;mu::avgPlaquette(Umu); + RealD org_link_trace=WilsonLoops::linkTrace(Umu); + RealD old_trace = org_link_trace; + RealD trG; + + std::vector U(Nd,grid); + GaugeMat dmuAmu(grid); + + for(int i=0;i(Umu,mu); + //trG = SteepestDescentStep(U,alpha,dmuAmu); + trG = FourierAccelSteepestDescentStep(U,alpha,dmuAmu); + for(int mu=0;mu(Umu,U[mu],mu); + // Monitor progress and convergence test + // infrequently to minimise cost overhead + if ( i %20 == 0 ) { + RealD plaq =WilsonLoops::avgPlaquette(Umu); + RealD link_trace=WilsonLoops::linkTrace(Umu); + + std::cout << GridLogMessage << " Iteration "< &U,RealD & alpha, GaugeMat & dmuAmu) { + GridBase *grid = U[0]._grid; + + std::vector A(Nd,grid); + GaugeMat g(grid); + + GaugeLinkToLieAlgebraField(U,A); + ExpiAlphaDmuAmu(A,g,alpha,dmuAmu); + + + RealD vol = grid->gSites(); + RealD trG = TensorRemove(sum(trace(g))).real()/vol/Nc; + + SU::GaugeTransform(U,g); + + return trG; + } + + static RealD FourierAccelSteepestDescentStep(std::vector &U,RealD & alpha, GaugeMat & dmuAmu) { + + GridBase *grid = U[0]._grid; + + RealD vol = grid->gSites(); + + FFT theFFT((GridCartesian *)grid); + + LatticeComplex Fp(grid); + LatticeComplex psq(grid); psq=zero; + LatticeComplex pmu(grid); + LatticeComplex one(grid); one = ComplexD(1.0,0.0); + + GaugeMat g(grid); + GaugeMat dmuAmu_p(grid); + std::vector A(Nd,grid); + + GaugeLinkToLieAlgebraField(U,A); + + DmuAmu(A,dmuAmu); + + theFFT.FFT_all_dim(dmuAmu_p,dmuAmu,FFT::forward); + + ////////////////////////////////// + // Work out Fp = psq_max/ psq... + ////////////////////////////////// + std::vector latt_size = grid->GlobalDimensions(); + std::vector coor(grid->_ndimension,0); + for(int mu=0;mu::taExp(ciadmam,g); + + RealD trG = TensorRemove(sum(trace(g))).real()/vol/Nc; + + SU::GaugeTransform(U,g); + + return trG; + } + + static void ExpiAlphaDmuAmu(const std::vector &A,GaugeMat &g,RealD & alpha, GaugeMat &dmuAmu) { + GridBase *grid = g._grid; + ComplexD cialpha(0.0,-alpha); + GaugeMat ciadmam(grid); + DmuAmu(A,dmuAmu); + ciadmam = dmuAmu*cialpha; + SU::taExp(ciadmam,g); + } +/* + //////////////////////////////////////////////////////////////// + // NB The FT for fields living on links has an extra phase in it + // Could add these to the FFT class as a later task since this code + // might be reused elsewhere ???? + //////////////////////////////////////////////////////////////// + static void InverseFourierTransformAmu(FFT &theFFT,const std::vector &Ap,std::vector &Ax) { + GridBase * grid = theFFT.Grid(); + std::vector latt_size = grid->GlobalDimensions(); + + ComplexField pmu(grid); + ComplexField pha(grid); + GaugeMat Apha(grid); + + ComplexD ci(0.0,1.0); + + for(int mu=0;mu &Ax,std::vector &Ap) { + GridBase * grid = theFFT.Grid(); + std::vector latt_size = grid->GlobalDimensions(); + + ComplexField pmu(grid); + ComplexField pha(grid); + ComplexD ci(0.0,1.0); + + // Sign convention for FFTW calls: + // A(x)= Sum_p e^ipx A(p) / V + // A(p)= Sum_p e^-ipx A(x) + + for(int mu=0;mu seeds({1,2,3,4}); + + Grid_init(&argc,&argv); + + int threads = GridThread::GetThreads(); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout( { vComplexD::Nsimd(),1,1,1}); + std::vector mpi_layout = GridDefaultMpi(); + + int vol = 1; + for(int d=0;d::avgPlaquette(Umu); + std::cout << " Initial plaquette "<::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-10, 1.0e-10); + + + plaq=WilsonLoops::avgPlaquette(Umu); + std::cout << " Final plaquette "< Date: Sat, 29 Oct 2016 11:04:02 +0100 Subject: [PATCH 47/58] Add missing volume factor in stochastic QED field --- lib/qcd/action/gauge/Photon.h | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/lib/qcd/action/gauge/Photon.h b/lib/qcd/action/gauge/Photon.h index 0dd5a9dc..409569d8 100644 --- a/lib/qcd/action/gauge/Photon.h +++ b/lib/qcd/action/gauge/Photon.h @@ -151,12 +151,19 @@ namespace QCD{ { auto *grid = dynamic_cast(out._grid); const unsigned int nd = grid->_ndimension; + std::vector latt_size = grid->_fdimensions; GaugeLinkField sqrtK2Inv(grid), r(grid); GaugeField aTilde(grid); FFT fft(grid); + Integer vol = 1; + for(int d = 0; d < nd; d++) + { + vol = vol * latt_size[d]; + } + invKHatSquared(sqrtK2Inv); - sqrtK2Inv = sqrt(real(sqrtK2Inv)); + sqrtK2Inv = sqrt(vol*real(sqrtK2Inv)); zmSub(sqrtK2Inv); for(int mu = 0; mu < nd; mu++) { From 791cb050c8e79bd0994e474044da5f055a9c05e8 Mon Sep 17 00:00:00 2001 From: paboyle Date: Tue, 1 Nov 2016 11:35:43 +0000 Subject: [PATCH 48/58] Comms improvements --- benchmarks/Benchmark_comms.cc | 86 ++- bootstrap.sh | 6 - configure.ac | 4 + lib/Cshift.h | 4 + lib/Makefile.am | 5 + lib/communicator/Communicator_base.cc | 15 +- lib/communicator/Communicator_base.h | 43 +- lib/communicator/Communicator_mpi.cc | 12 +- lib/communicator/Communicator_mpi3.cc | 12 +- lib/communicator/Communicator_mpi3_leader.cc | 733 +++++++++++++++++++ lib/communicator/Communicator_none.cc | 8 +- lib/communicator/Communicator_shmem.cc | 31 +- scripts/arm_configure.experimental | 1 - scripts/arm_configure.experimental_cortex57 | 3 - scripts/bench_wilson.sh | 9 - scripts/build-all | 46 -- scripts/configure-all | 11 - scripts/configure-commands | 89 --- scripts/configure-cray | 10 - scripts/configure-mic | 10 - scripts/copyright | 16 +- scripts/cray-modules | 2 - scripts/reconfigure_script | 4 - scripts/update_fftw.sh | 18 - scripts/wilson.gnu | 7 - 25 files changed, 898 insertions(+), 287 deletions(-) create mode 100644 lib/communicator/Communicator_mpi3_leader.cc delete mode 100644 scripts/arm_configure.experimental delete mode 100644 scripts/arm_configure.experimental_cortex57 delete mode 100755 scripts/bench_wilson.sh delete mode 100755 scripts/build-all delete mode 100755 scripts/configure-all delete mode 100755 scripts/configure-commands delete mode 100755 scripts/configure-cray delete mode 100755 scripts/configure-mic delete mode 100644 scripts/cray-modules delete mode 100755 scripts/reconfigure_script delete mode 100755 scripts/update_fftw.sh delete mode 100644 scripts/wilson.gnu diff --git a/benchmarks/Benchmark_comms.cc b/benchmarks/Benchmark_comms.cc index d0cd96c7..234d2fbb 100644 --- a/benchmarks/Benchmark_comms.cc +++ b/benchmarks/Benchmark_comms.cc @@ -42,14 +42,13 @@ int main (int argc, char ** argv) int Nloop=10; int nmu=0; - for(int mu=0;mu<4;mu++) if (mpi_layout[mu]>1) nmu++; + for(int mu=0;mu1) nmu++; + std::cout< latt_size ({lat*mpi_layout[0], + lat*mpi_layout[1], + lat*mpi_layout[2], + lat*mpi_layout[3]}); + + GridCartesian Grid(latt_size,simd_layout,mpi_layout); + + std::vector xbuf(8); + std::vector rbuf(8); + Grid.ShmBufferFreeAll(); + for(int d=0;d<8;d++){ + xbuf[d] = (HalfSpinColourVectorD *)Grid.ShmBufferMalloc(lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD)); + rbuf[d] = (HalfSpinColourVectorD *)Grid.ShmBufferMalloc(lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD)); + } + + int ncomm; + int bytes=lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD); + + double start=usecond(); + for(int i=0;i requests; + + ncomm=0; + for(int mu=0;mu<4;mu++){ + + if (mpi_layout[mu]>1 ) { + + ncomm++; + int comm_proc=1; + int xmit_to_rank; + int recv_from_rank; + + Grid.ShiftedRanks(mu,comm_proc,xmit_to_rank,recv_from_rank); + Grid.StencilSendToRecvFromBegin(requests, + (void *)&xbuf[mu][0], + xmit_to_rank, + (void *)&rbuf[mu][0], + recv_from_rank, + bytes); + + comm_proc = mpi_layout[mu]-1; + + Grid.ShiftedRanks(mu,comm_proc,xmit_to_rank,recv_from_rank); + Grid.StencilSendToRecvFromBegin(requests, + (void *)&xbuf[mu+4][0], + xmit_to_rank, + (void *)&rbuf[mu+4][0], + recv_from_rank, + bytes); + + } + } + Grid.StencilSendToRecvFromComplete(requests); + Grid.Barrier(); + + } + double stop=usecond(); + + double dbytes = bytes; + double xbytes = Nloop*dbytes*2.0*ncomm; + double rbytes = xbytes; + double bidibytes = xbytes+rbytes; + + double time = stop-start; // microseconds + + std::cout< #include #endif +#ifdef GRID_COMMS_MPI3L +#include +#endif + #ifdef GRID_COMMS_SHMEM #include // uses same implementation of communicator #endif diff --git a/lib/Makefile.am b/lib/Makefile.am index 3285a545..a779135f 100644 --- a/lib/Makefile.am +++ b/lib/Makefile.am @@ -9,6 +9,11 @@ if BUILD_COMMS_MPI3 extra_sources+=communicator/Communicator_base.cc endif +if BUILD_COMMS_MPI3L + extra_sources+=communicator/Communicator_mpi3_leader.cc + extra_sources+=communicator/Communicator_base.cc +endif + if BUILD_COMMS_SHMEM extra_sources+=communicator/Communicator_shmem.cc extra_sources+=communicator/Communicator_base.cc diff --git a/lib/communicator/Communicator_base.cc b/lib/communicator/Communicator_base.cc index 9810987d..77b4f800 100644 --- a/lib/communicator/Communicator_base.cc +++ b/lib/communicator/Communicator_base.cc @@ -31,13 +31,6 @@ namespace Grid { /////////////////////////////////////////////////////////////// // Info that is setup once and indept of cartesian layout /////////////////////////////////////////////////////////////// -int CartesianCommunicator::ShmRank; -int CartesianCommunicator::ShmSize; -int CartesianCommunicator::GroupRank; -int CartesianCommunicator::GroupSize; -int CartesianCommunicator::WorldRank; -int CartesianCommunicator::WorldSize; -int CartesianCommunicator::Slave; void * CartesianCommunicator::ShmCommBuf; ///////////////////////////////// @@ -70,12 +63,6 @@ int CartesianCommunicator::ProcessorCount(void) { return //////////////////////////////////////////////////////////////////////////////// // very VERY rarely (Log, serial RNG) we need world without a grid //////////////////////////////////////////////////////////////////////////////// -int CartesianCommunicator::RankWorld(void){ return WorldRank; }; -int CartesianCommunicator::Ranks (void) { return WorldSize; }; -int CartesianCommunicator::Nodes (void) { return GroupSize; }; -int CartesianCommunicator::Cores (void) { return ShmSize; }; -int CartesianCommunicator::NodeRank (void) { return GroupRank; }; -int CartesianCommunicator::CoreRank (void) { return ShmRank; }; void CartesianCommunicator::GlobalSum(ComplexF &c) { @@ -94,7 +81,7 @@ void CartesianCommunicator::GlobalSumVector(ComplexD *c,int N) GlobalSumVector((double *)c,2*N); } -#ifndef GRID_COMMS_MPI3 +#if !defined( GRID_COMMS_MPI3) && !defined (GRID_COMMS_MPI3L) void CartesianCommunicator::StencilSendToRecvFromBegin(std::vector &list, void *xmit, diff --git a/lib/communicator/Communicator_base.h b/lib/communicator/Communicator_base.h index 9a15d2b0..d4bb50d9 100644 --- a/lib/communicator/Communicator_base.h +++ b/lib/communicator/Communicator_base.h @@ -1,3 +1,4 @@ + /************************************************************************************* Grid physics library, www.github.com/paboyle/Grid @@ -37,6 +38,9 @@ Author: Peter Boyle #ifdef GRID_COMMS_MPI3 #include #endif +#ifdef GRID_COMMS_MPI3L +#include +#endif #ifdef GRID_COMMS_SHMEM #include #endif @@ -60,9 +64,9 @@ class CartesianCommunicator { std::vector _processor_coor; // linear processor coordinate unsigned long _ndimension; -#if defined (GRID_COMMS_MPI) || defined (GRID_COMMS_MPI3) - MPI_Comm communicator; +#if defined (GRID_COMMS_MPI) || defined (GRID_COMMS_MPI3) || defined (GRID_COMMS_MPI3L) static MPI_Comm communicator_world; + MPI_Comm communicator; typedef MPI_Request CommsRequest_t; #else typedef int CommsRequest_t; @@ -75,7 +79,15 @@ class CartesianCommunicator { // cartesian communicator on a subset of ranks, slave ranks controlled // by group leader with data xfer via shared memory //////////////////////////////////////////////////////////////////// -#ifdef GRID_COMMS_MPI3 +#ifdef GRID_COMMS_MPI3 + + static int ShmRank; + static int ShmSize; + static int GroupRank; + static int GroupSize; + static int WorldRank; + static int WorldSize; + std::vector WorldDims; std::vector GroupDims; std::vector ShmDims; @@ -83,7 +95,7 @@ class CartesianCommunicator { std::vector GroupCoor; std::vector ShmCoor; std::vector WorldCoor; - + static std::vector GroupRanks; static std::vector MyGroup; static int ShmSetup; @@ -93,13 +105,20 @@ class CartesianCommunicator { std::vector LexicographicToWorldRank; static std::vector ShmCommBufs; + #else static void ShmInitGeneric(void); static commVector ShmBufStorageVector; #endif + + ///////////////////////////////// + // Grid information and queries + // Implemented in Communicator_base.C + ///////////////////////////////// static void * ShmCommBuf; size_t heap_top; size_t heap_bytes; + void *ShmBufferSelf(void); void *ShmBuffer(int rank); void *ShmBufferTranslate(int rank,void * local_p); @@ -123,28 +142,12 @@ class CartesianCommunicator { int RankFromProcessorCoor(std::vector &coor); void ProcessorCoorFromRank(int rank,std::vector &coor); - ///////////////////////////////// - // Grid information and queries - ///////////////////////////////// - static int ShmRank; - static int ShmSize; - static int GroupSize; - static int GroupRank; - static int WorldRank; - static int WorldSize; - static int Slave; - int IsBoss(void) ; int BossRank(void) ; int ThisRank(void) ; const std::vector & ThisProcessorCoor(void) ; const std::vector & ProcessorGrid(void) ; int ProcessorCount(void) ; - static int Ranks (void); - static int Nodes (void); - static int Cores (void); - static int NodeRank (void); - static int CoreRank (void); //////////////////////////////////////////////////////////////////////////////// // very VERY rarely (Log, serial RNG) we need world without a grid diff --git a/lib/communicator/Communicator_mpi.cc b/lib/communicator/Communicator_mpi.cc index a638eebb..65ced9c7 100644 --- a/lib/communicator/Communicator_mpi.cc +++ b/lib/communicator/Communicator_mpi.cc @@ -44,13 +44,6 @@ void CartesianCommunicator::Init(int *argc, char ***argv) { MPI_Init(argc,argv); } MPI_Comm_dup (MPI_COMM_WORLD,&communicator_world); - MPI_Comm_rank(communicator_world,&WorldRank); - MPI_Comm_size(communicator_world,&WorldSize); - ShmRank=0; - ShmSize=1; - GroupRank=WorldRank; - GroupSize=WorldSize; - Slave =0; ShmInitGeneric(); } @@ -198,6 +191,11 @@ void CartesianCommunicator::Broadcast(int root,void* data, int bytes) // Should only be used prior to Grid Init finished. // Check for this? /////////////////////////////////////////////////////// +int CartesianCommunicator::RankWorld(void){ + int r; + MPI_Comm_rank(communicator_world,&r); + return r; +} void CartesianCommunicator::BroadcastWorld(int root,void* data, int bytes) { int ierr= MPI_Bcast(data, diff --git a/lib/communicator/Communicator_mpi3.cc b/lib/communicator/Communicator_mpi3.cc index 14a152ab..c707ec1f 100644 --- a/lib/communicator/Communicator_mpi3.cc +++ b/lib/communicator/Communicator_mpi3.cc @@ -30,12 +30,18 @@ Author: Peter Boyle namespace Grid { - /////////////////////////////////////////////////////////////////////////////////////////////////// // Info that is setup once and indept of cartesian layout /////////////////////////////////////////////////////////////////////////////////////////////////// int CartesianCommunicator::ShmSetup = 0; +int CartesianCommunicator::ShmRank; +int CartesianCommunicator::ShmSize; +int CartesianCommunicator::GroupRank; +int CartesianCommunicator::GroupSize; +int CartesianCommunicator::WorldRank; +int CartesianCommunicator::WorldSize; + MPI_Comm CartesianCommunicator::communicator_world; MPI_Comm CartesianCommunicator::ShmComm; MPI_Win CartesianCommunicator::ShmWindow; @@ -97,15 +103,15 @@ void CartesianCommunicator::Init(int *argc, char ***argv) { std::vector world_ranks(WorldSize); GroupRanks.resize(WorldSize); - MyGroup.resize(ShmSize); for(int r=0;r + + 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.h" +#include + +namespace Grid { + +enum { COMMAND_ISEND, COMMAND_IRECV, COMMAND_WAITALL }; + +struct Descriptor { + uint64_t buf; + size_t bytes; + int rank; + int tag; + int command; + MPI_Request request; +}; + +const int pool = 48; + +class SlaveState { +public: + volatile int head; + volatile int start; + volatile int tail; + volatile Descriptor Descrs[pool]; +}; + +class Slave { +public: + SlaveState *state; + MPI_Comm squadron; + uint64_t base; + //////////////////////////////////////////////////////////// + // Descriptor circular pointers + //////////////////////////////////////////////////////////// + Slave() {}; + + void Init(SlaveState * _state,MPI_Comm _squadron); + + void EventLoop (void) { + std::cerr<< " Entering even loop "<tail != state->head ); + std::cerr<< " FIFO drained "<< state->tail < Slaves; + + static int ShmSetup; + + static int UniverseRank; + static int UniverseSize; + + static MPI_Comm communicator_universe; + static MPI_Comm communicator_cached; + + static MPI_Comm HorizontalComm; + static int HorizontalRank; + static int HorizontalSize; + + static MPI_Comm VerticalComm; + static MPI_Win VerticalWindow; + static int VerticalSize; + static int VerticalRank; + + static std::vector VerticalShmBufs; + static std::vector > UniverseRanks; + static std::vector UserCommunicatorToWorldRanks; + + static MPI_Group WorldGroup, CachedGroup; + + static void CommunicatorInit (MPI_Comm &communicator_world, + MPI_Comm &ShmComm, + void * &ShmCommBuf); + + static void MapCommRankToWorldRank(int &hashtag, int & comm_world_peer,int tag, MPI_Comm comm,int rank); + + ///////////////////////////////////////////////////////// + // routines for master proc must handle any communicator + ///////////////////////////////////////////////////////// + + static uint64_t QueueSend(int slave,void *buf, int bytes, int tag, MPI_Comm comm,int rank) { + std::cerr<< " Queueing send "<< bytes<= units ) { + mywork = myoff = 0; + } else { + mywork = (nwork+me)/units; + myoff = basework * me; + if ( me > backfill ) + myoff+= (me-backfill); + } + return; + }; + + static void QueueMultiplexedSend(void *buf, int bytes, int tag, MPI_Comm comm,int rank) { + uint8_t * cbuf = (uint8_t *) buf; + int mywork, myoff, procs; + procs = VerticalSize-1; + for(int s=0;s MPIoffloadEngine::Slaves; + +int MPIoffloadEngine::UniverseRank; +int MPIoffloadEngine::UniverseSize; + +MPI_Comm MPIoffloadEngine::communicator_universe; +MPI_Comm MPIoffloadEngine::communicator_cached; +MPI_Group MPIoffloadEngine::WorldGroup; +MPI_Group MPIoffloadEngine::CachedGroup; + +MPI_Comm MPIoffloadEngine::HorizontalComm; +int MPIoffloadEngine::HorizontalRank; +int MPIoffloadEngine::HorizontalSize; + +MPI_Comm MPIoffloadEngine::VerticalComm; +MPI_Win MPIoffloadEngine::VerticalWindow; +int MPIoffloadEngine::VerticalSize; +int MPIoffloadEngine::VerticalRank; + +std::vector MPIoffloadEngine::VerticalShmBufs; +std::vector > MPIoffloadEngine::UniverseRanks; +std::vector MPIoffloadEngine::UserCommunicatorToWorldRanks; + +int MPIoffloadEngine::ShmSetup = 0; + +void MPIoffloadEngine::CommunicatorInit (MPI_Comm &communicator_world, + MPI_Comm &ShmComm, + void * &ShmCommBuf) +{ + int flag; + assert(ShmSetup==0); + + ////////////////////////////////////////////////////////////////////// + // Universe is all nodes prior to squadron grouping + ////////////////////////////////////////////////////////////////////// + MPI_Comm_dup (MPI_COMM_WORLD,&communicator_universe); + MPI_Comm_rank(communicator_universe,&UniverseRank); + MPI_Comm_size(communicator_universe,&UniverseSize); + + ///////////////////////////////////////////////////////////////////// + // Split into groups that can share memory (Verticals) + ///////////////////////////////////////////////////////////////////// + // MPI_Comm_split_type(communicator_universe, MPI_COMM_TYPE_SHARED, 0, MPI_INFO_NULL,&VerticalComm); + MPI_Comm_split(communicator_universe,(UniverseRank&0x1),UniverseRank,&VerticalComm); + MPI_Comm_rank(VerticalComm ,&VerticalRank); + MPI_Comm_size(VerticalComm ,&VerticalSize); + + ////////////////////////////////////////////////////////////////////// + // Split into horizontal groups by rank in squadron + ////////////////////////////////////////////////////////////////////// + MPI_Comm_split(communicator_universe,VerticalRank,UniverseRank,&HorizontalComm); + MPI_Comm_rank(HorizontalComm,&HorizontalRank); + MPI_Comm_size(HorizontalComm,&HorizontalSize); + assert(HorizontalSize*VerticalSize==UniverseSize); + + //////////////////////////////////////////////////////////////////////////////// + // What is my place in the world + //////////////////////////////////////////////////////////////////////////////// + int WorldRank=0; + if(VerticalRank==0) WorldRank = HorizontalRank; + int ierr=MPI_Allreduce(MPI_IN_PLACE,&WorldRank,1,MPI_INT,MPI_SUM,VerticalComm); + assert(ierr==0); + + //////////////////////////////////////////////////////////////////////////////// + // Where is the world in the universe? + //////////////////////////////////////////////////////////////////////////////// + UniverseRanks = std::vector >(HorizontalSize,std::vector(VerticalSize,0)); + UniverseRanks[WorldRank][VerticalRank] = UniverseRank; + for(int w=0;w cached_ranks(size); + + for(int r=0;r>0 )&0xFFFF)^((icomm>>16)&0xFFFF) + ^ ((icomm>>32)&0xFFFF)^((icomm>>48)&0xFFFF); + + hashtag = (comm_hash<<15) | tag; + +}; + +void Slave::Init(SlaveState * _state,MPI_Comm _squadron) +{ + squadron=_squadron; + state =_state; + state->head = state->tail = state->start = 0; + MPI_Barrier(squadron); + base = (uint64_t)MPIoffloadEngine::VerticalShmBufs[0]; + int rank; MPI_Comm_rank(_squadron,&rank); +} +#define PERI_PLUS(A) ( (A+1)%pool ) +void Slave::Event (void) { + + static int tail_last; + static int head_last; + static int start_last; + int ierr; + + if ( (state->tail != tail_last) + ||(state->head != head_last) + ||(state->start != start_last) + ) { + std::cerr<< " Event loop "<< state->tail << " "<< state->start<< " "<< state->head <tail; + if ( t != state->start ) { + int flag=0; + + std::cerr<< " Testing tail "<< t<<" "<< (void *)&state->Descrs[t].request + << " "<Descrs[t].request<Descrs[t].request,&flag,MPI_STATUS_IGNORE); + // ierr=MPI_Test((MPI_Request *)&state->Descrs[t].request,&flag,MPI_STATUS_IGNORE); + assert(ierr==0); + if ( flag ) { + state->tail = PERI_PLUS(t); + std::cerr<< " Tail advanced from "<< t<start; + if ( s != state->head ) { + switch ( state->Descrs[s].command ) { + case COMMAND_ISEND: + std::cerr<< " Send "<Descrs[s].buf<< "["<Descrs[s].bytes<<"]" + << " to " << state->Descrs[s].rank<< " tag" << state->Descrs[s].tag + << " Comm " << MPIoffloadEngine::communicator_universe<< std::endl; + + std::cerr<< " Request was "<Descrs[s].request<Descrs[s].buf+base), + state->Descrs[s].bytes, + MPI_CHAR, + state->Descrs[s].rank, + state->Descrs[s].tag, + MPIoffloadEngine::communicator_universe, + (MPI_Request *)&state->Descrs[s].request); + std::cerr<< " Request is "<Descrs[s].request<Descrs[0].request<start = PERI_PLUS(s); + break; + + case COMMAND_IRECV: + std::cerr<< " Recv "<Descrs[s].buf<< "["<Descrs[s].bytes<<"]" + << " to " << state->Descrs[s].rank<< " tag" << state->Descrs[s].tag + << " Comm " << MPIoffloadEngine::communicator_universe<< std::endl; + + std::cerr<< " Request was "<Descrs[s].request<Descrs[s].buf+base), + state->Descrs[s].bytes, + MPI_CHAR, + state->Descrs[s].rank, + state->Descrs[s].tag, + MPIoffloadEngine::communicator_universe, + (MPI_Request *)&state->Descrs[s].request); + std::cerr<< " Request is "<Descrs[s].request<Descrs[0].request<start = PERI_PLUS(s); + break; + + case COMMAND_WAITALL: + std::cerr<< " Wait all "<tail;t!=s; t=PERI_PLUS(t) ){ + std::cerr<< " Wait ["<Descrs[t].request <Descrs[0].request<Descrs[t].request,MPI_STATUS_IGNORE); + }; + s=PERI_PLUS(s); + state->start = s; + state->tail = s; + break; + + default: + assert(0); + break; + } + return; + } +} + ////////////////////////////////////////////////////////////////////////////// + // External interaction with the queue + ////////////////////////////////////////////////////////////////////////////// + +uint64_t Slave::QueueCommand(int command,void *buf, int bytes, int tag, MPI_Comm comm,int commrank) +{ + ///////////////////////////////////////// + // Spin; if FIFO is full until not full + ///////////////////////////////////////// + int head =state->head; + int next = PERI_PLUS(head); + + // Set up descriptor + int worldrank; + int hashtag; + MPI_Comm communicator; + MPI_Request request; + + MPIoffloadEngine::MapCommRankToWorldRank(hashtag,commrank,tag,comm,worldrank); + + int VerticalRank = MPIoffloadEngine::VerticalRank; + uint64_t relative= (uint64_t)buf - base; + state->Descrs[head].buf = relative; + state->Descrs[head].bytes = bytes; + state->Descrs[head].rank = MPIoffloadEngine::UniverseRanks[worldrank][VerticalRank]; + state->Descrs[head].tag = hashtag; + state->Descrs[head].command= command; + std::cerr<< " QueueCommand "<tail==next ); + + // Msync on weak order architectures + // Advance pointer + state->head = next; + + return 0; +} + + +/////////////////////////////////////////////////////////////////////////////////////////////////// +// Info that is setup once and indept of cartesian layout +/////////////////////////////////////////////////////////////////////////////////////////////////// + +MPI_Comm CartesianCommunicator::communicator_world; + +void CartesianCommunicator::Init(int *argc, char ***argv) +{ + int flag; + MPI_Initialized(&flag); // needed to coexist with other libs apparently + if ( !flag ) { + MPI_Init(argc,argv); + } + communicator_world = MPI_COMM_WORLD; + MPI_Comm ShmComm; + MPIoffloadEngine::CommunicatorInit (communicator_world,ShmComm,ShmCommBuf); +} +void CartesianCommunicator::ShiftedRanks(int dim,int shift,int &source,int &dest) +{ + int ierr=MPI_Cart_shift(communicator,dim,shift,&source,&dest); + assert(ierr==0); +} +int CartesianCommunicator::RankFromProcessorCoor(std::vector &coor) +{ + int rank; + int ierr=MPI_Cart_rank (communicator, &coor[0], &rank); + assert(ierr==0); + return rank; +} +void CartesianCommunicator::ProcessorCoorFromRank(int rank, std::vector &coor) +{ + coor.resize(_ndimension); + int ierr=MPI_Cart_coords (communicator, rank, _ndimension,&coor[0]); + assert(ierr==0); +} + +CartesianCommunicator::CartesianCommunicator(const std::vector &processors) +{ + _ndimension = processors.size(); + std::vector periodic(_ndimension,1); + + _Nprocessors=1; + _processors = processors; + + for(int i=0;i<_ndimension;i++){ + _Nprocessors*=_processors[i]; + } + + int Size; + MPI_Comm_size(communicator_world,&Size); + assert(Size==_Nprocessors); + + _processor_coor.resize(_ndimension); + MPI_Cart_create(communicator_world, _ndimension,&_processors[0],&periodic[0],1,&communicator); + MPI_Comm_rank (communicator,&_processor); + MPI_Cart_coords(communicator,_processor,_ndimension,&_processor_coor[0]); +}; + +void CartesianCommunicator::GlobalSum(uint32_t &u){ + int ierr=MPI_Allreduce(MPI_IN_PLACE,&u,1,MPI_UINT32_T,MPI_SUM,communicator); + assert(ierr==0); +} +void CartesianCommunicator::GlobalSum(uint64_t &u){ + int ierr=MPI_Allreduce(MPI_IN_PLACE,&u,1,MPI_UINT64_T,MPI_SUM,communicator); + assert(ierr==0); +} +void CartesianCommunicator::GlobalSum(float &f){ + int ierr=MPI_Allreduce(MPI_IN_PLACE,&f,1,MPI_FLOAT,MPI_SUM,communicator); + assert(ierr==0); +} +void CartesianCommunicator::GlobalSumVector(float *f,int N) +{ + int ierr=MPI_Allreduce(MPI_IN_PLACE,f,N,MPI_FLOAT,MPI_SUM,communicator); + assert(ierr==0); +} +void CartesianCommunicator::GlobalSum(double &d) +{ + int ierr = MPI_Allreduce(MPI_IN_PLACE,&d,1,MPI_DOUBLE,MPI_SUM,communicator); + assert(ierr==0); +} +void CartesianCommunicator::GlobalSumVector(double *d,int N) +{ + int ierr = MPI_Allreduce(MPI_IN_PLACE,d,N,MPI_DOUBLE,MPI_SUM,communicator); + assert(ierr==0); +} + +// Basic Halo comms primitive +void CartesianCommunicator::SendToRecvFrom(void *xmit, + int dest, + void *recv, + int from, + int bytes) +{ + std::vector reqs(0); + SendToRecvFromBegin(reqs,xmit,dest,recv,from,bytes); + SendToRecvFromComplete(reqs); +} + +void CartesianCommunicator::SendRecvPacket(void *xmit, + void *recv, + int sender, + int receiver, + int bytes) +{ + MPI_Status stat; + assert(sender != receiver); + int tag = sender; + if ( _processor == sender ) { + MPI_Send(xmit, bytes, MPI_CHAR,receiver,tag,communicator); + } + if ( _processor == receiver ) { + MPI_Recv(recv, bytes, MPI_CHAR,sender,tag,communicator,&stat); + } +} + +// Basic Halo comms primitive +void CartesianCommunicator::SendToRecvFromBegin(std::vector &list, + void *xmit, + int dest, + void *recv, + int from, + int bytes) +{ + MPI_Request xrq; + MPI_Request rrq; + int rank = _processor; + int ierr; + ierr =MPI_Isend(xmit, bytes, MPI_CHAR,dest,_processor,communicator,&xrq); + ierr|=MPI_Irecv(recv, bytes, MPI_CHAR,from,from,communicator,&rrq); + + assert(ierr==0); + + list.push_back(xrq); + list.push_back(rrq); +} + +void CartesianCommunicator::StencilSendToRecvFromBegin(std::vector &list, + void *xmit, + int dest, + void *recv, + int from, + int bytes) +{ + uint64_t xmit_i = (uint64_t) xmit; + uint64_t recv_i = (uint64_t) recv; + uint64_t shm = (uint64_t) ShmCommBuf; + // assert xmit and recv lie in shared memory region + assert( (xmit_i >= shm) && (xmit_i+bytes <= shm+MAX_MPI_SHM_BYTES) ); + assert( (recv_i >= shm) && (recv_i+bytes <= shm+MAX_MPI_SHM_BYTES) ); + MPIoffloadEngine::QueueMultiplexedSend(xmit,bytes,_processor,communicator,dest); + MPIoffloadEngine::QueueMultiplexedRecv(recv,bytes,from,communicator,from); +} + + +void CartesianCommunicator::StencilSendToRecvFromComplete(std::vector &list) +{ + MPIoffloadEngine::WaitAll(); +} + +void CartesianCommunicator::StencilBarrier(void) +{ +} + +void CartesianCommunicator::SendToRecvFromComplete(std::vector &list) +{ + int nreq=list.size(); + std::vector status(nreq); + int ierr = MPI_Waitall(nreq,&list[0],&status[0]); + assert(ierr==0); +} + +void CartesianCommunicator::Barrier(void) +{ + int ierr = MPI_Barrier(communicator); + assert(ierr==0); +} + +void CartesianCommunicator::Broadcast(int root,void* data, int bytes) +{ + int ierr=MPI_Bcast(data, + bytes, + MPI_BYTE, + root, + communicator); + assert(ierr==0); +} + +void CartesianCommunicator::BroadcastWorld(int root,void* data, int bytes) +{ + int ierr= MPI_Bcast(data, + bytes, + MPI_BYTE, + root, + communicator_world); + assert(ierr==0); +} + +void *CartesianCommunicator::ShmBufferSelf(void) { return ShmCommBuf; } + +void *CartesianCommunicator::ShmBuffer(int rank) { + return NULL; +} +void *CartesianCommunicator::ShmBufferTranslate(int rank,void * local_p) { + return NULL; +} + + +}; + diff --git a/lib/communicator/Communicator_none.cc b/lib/communicator/Communicator_none.cc index 198c1add..0f43f1f5 100644 --- a/lib/communicator/Communicator_none.cc +++ b/lib/communicator/Communicator_none.cc @@ -34,13 +34,6 @@ namespace Grid { void CartesianCommunicator::Init(int *argc, char *** arv) { - WorldRank = 0; - WorldSize = 1; - ShmRank=0; - ShmSize=1; - GroupRank=WorldRank; - GroupSize=WorldSize; - Slave =0; ShmInitGeneric(); } @@ -99,6 +92,7 @@ void CartesianCommunicator::SendToRecvFromComplete(std::vector & assert(0); } +int CartesianCommunicator::RankWorld(void){return 0;} void CartesianCommunicator::Barrier(void){} void CartesianCommunicator::Broadcast(int root,void* data, int bytes) {} void CartesianCommunicator::BroadcastWorld(int root,void* data, int bytes) { } diff --git a/lib/communicator/Communicator_shmem.cc b/lib/communicator/Communicator_shmem.cc index d5f3ed76..56e03224 100644 --- a/lib/communicator/Communicator_shmem.cc +++ b/lib/communicator/Communicator_shmem.cc @@ -50,11 +50,16 @@ typedef struct HandShake_t { uint64_t seq_remote; } HandShake; +std::array make_psync_init(void) { + array ret; + ret.fill(SHMEM_SYNC_VALUE); + return ret; +} +static std::array psync_init = make_psync_init(); static Vector< HandShake > XConnections; static Vector< HandShake > RConnections; - void CartesianCommunicator::Init(int *argc, char ***argv) { shmem_init(); XConnections.resize(shmem_n_pes()); @@ -65,13 +70,6 @@ void CartesianCommunicator::Init(int *argc, char ***argv) { RConnections[pe].seq_local = 0; RConnections[pe].seq_remote= 0; } - WorldSize = shmem_n_pes(); - WorldRank = shmem_my_pe(); - ShmRank=0; - ShmSize=1; - GroupRank=WorldRank; - GroupSize=WorldSize; - Slave =0; shmem_barrier_all(); ShmInitGeneric(); } @@ -103,7 +101,7 @@ void CartesianCommunicator::GlobalSum(uint32_t &u){ static long long source ; static long long dest ; static long long llwrk[_SHMEM_REDUCE_MIN_WRKDATA_SIZE]; - static long psync[_SHMEM_REDUCE_SYNC_SIZE]; + static std::array psync = psync_init; // int nreduce=1; // int pestart=0; @@ -119,7 +117,7 @@ void CartesianCommunicator::GlobalSum(uint64_t &u){ static long long source ; static long long dest ; static long long llwrk[_SHMEM_REDUCE_MIN_WRKDATA_SIZE]; - static long psync[_SHMEM_REDUCE_SYNC_SIZE]; + static std::array psync = psync_init; // int nreduce=1; // int pestart=0; @@ -135,7 +133,7 @@ void CartesianCommunicator::GlobalSum(float &f){ static float source ; static float dest ; static float llwrk[_SHMEM_REDUCE_MIN_WRKDATA_SIZE]; - static long psync[_SHMEM_REDUCE_SYNC_SIZE]; + static std::array psync = psync_init; source = f; dest =0.0; @@ -147,7 +145,7 @@ void CartesianCommunicator::GlobalSumVector(float *f,int N) static float source ; static float dest = 0 ; static float llwrk[_SHMEM_REDUCE_MIN_WRKDATA_SIZE]; - static long psync[_SHMEM_REDUCE_SYNC_SIZE]; + static std::array psync = psync_init; if ( shmem_addr_accessible(f,_processor) ){ shmem_float_sum_to_all(f,f,N,0,0,_Nprocessors,llwrk,psync); @@ -166,7 +164,7 @@ void CartesianCommunicator::GlobalSum(double &d) static double source; static double dest ; static double llwrk[_SHMEM_REDUCE_MIN_WRKDATA_SIZE]; - static long psync[_SHMEM_REDUCE_SYNC_SIZE]; + static std::array psync = psync_init; source = d; dest = 0; @@ -178,7 +176,8 @@ void CartesianCommunicator::GlobalSumVector(double *d,int N) static double source ; static double dest ; static double llwrk[_SHMEM_REDUCE_MIN_WRKDATA_SIZE]; - static long psync[_SHMEM_REDUCE_SYNC_SIZE]; + static std::array psync = psync_init; + if ( shmem_addr_accessible(d,_processor) ){ shmem_double_sum_to_all(d,d,N,0,0,_Nprocessors,llwrk,psync); @@ -295,7 +294,7 @@ void CartesianCommunicator::Barrier(void) } void CartesianCommunicator::Broadcast(int root,void* data, int bytes) { - static long psync[_SHMEM_REDUCE_SYNC_SIZE]; + static std::array psync = psync_init; static uint32_t word; uint32_t *array = (uint32_t *) data; assert( (bytes % 4)==0); @@ -318,7 +317,7 @@ void CartesianCommunicator::Broadcast(int root,void* data, int bytes) } void CartesianCommunicator::BroadcastWorld(int root,void* data, int bytes) { - static long psync[_SHMEM_REDUCE_SYNC_SIZE]; + static std::array psync = psync_init; static uint32_t word; uint32_t *array = (uint32_t *) data; assert( (bytes % 4)==0); diff --git a/scripts/arm_configure.experimental b/scripts/arm_configure.experimental deleted file mode 100644 index fbad84a5..00000000 --- a/scripts/arm_configure.experimental +++ /dev/null @@ -1 +0,0 @@ -./configure --host=arm-linux-gnueabihf CXX=clang++-3.5 CXXFLAGS='-std=c++11 -O3 -target arm-linux-gnueabihf -I/usr/arm-linux-gnueabihf/include/ -I/home/neo/Codes/gmp6.0/gmp-arm/include/ -I/usr/arm-linux-gnueabihf/include/c++/4.8.2/arm-linux-gnueabihf/ -L/home/neo/Codes/gmp6.0/gmp-arm/lib/ -I/home/neo/Codes/mpfr3.1.2/mpfr-arm/include/ -L/home/neo/Codes/mpfr3.1.2/mpfr-arm/lib/ -static -mcpu=cortex-a7' --enable-simd=NEONv7 diff --git a/scripts/arm_configure.experimental_cortex57 b/scripts/arm_configure.experimental_cortex57 deleted file mode 100644 index d229763e..00000000 --- a/scripts/arm_configure.experimental_cortex57 +++ /dev/null @@ -1,3 +0,0 @@ -#./configure --host=arm-linux-gnueabihf CXX=clang++-3.5 CXXFLAGS='-std=c++11 -O3 -target arm-linux-gnueabihf -I/usr/arm-linux-gnueabihf/include/ -I/home/neo/Codes/gmp6.0/gmp-arm/include/ -I/usr/lib/llvm-3.5/lib/clang/3.5.0/include/ -L/home/neo/Codes/gmp6.0/gmp-arm/lib/ -I/home/neo/Codes/mpfr3.1.2/mpfr-arm/include/ -L/home/neo/Codes/mpfr3.1.2/mpfr-arm/lib/ -static -mcpu=cortex-a57' --enable-simd=NEONv7 - -./configure --host=aarch64-linux-gnu CXX=clang++-3.5 CXXFLAGS='-std=c++11 -O3 -target aarch64-linux-gnu -static -I/home/neo/Codes/gmp6.0/gmp-armv8/include/ -L/home/neo/Codes/gmp6.0/gmp-armv8/lib/ -I/home/neo/Codes/mpfr3.1.2/mpfr-armv8/include/ -L/home/neo/Codes/mpfr3.1.2/mpfr-armv8/lib/ -I/usr/aarch64-linux-gnu/include/ -I/usr/aarch64-linux-gnu/include/c++/4.8.2/aarch64-linux-gnu/' --enable-simd=NEONv7 diff --git a/scripts/bench_wilson.sh b/scripts/bench_wilson.sh deleted file mode 100755 index af73d591..00000000 --- a/scripts/bench_wilson.sh +++ /dev/null @@ -1,9 +0,0 @@ -for omp in 1 2 4 -do -echo > wilson.t$omp -for vol in 4.4.4.4 4.4.4.8 4.4.8.8 4.8.8.8 8.8.8.8 8.8.8.16 8.8.16.16 8.16.16.16 -do -perf=` ./benchmarks/Grid_wilson --grid $vol --omp $omp | grep mflop | awk '{print $3}'` -echo $vol $perf >> wilson.t$omp -done -done \ No newline at end of file diff --git a/scripts/build-all b/scripts/build-all deleted file mode 100755 index b97dca19..00000000 --- a/scripts/build-all +++ /dev/null @@ -1,46 +0,0 @@ -#!/bin/bash -e - -DIRS="clang-avx clang-avx-openmp clang-avx-openmp-mpi clang-avx-mpi clang-avx2 clang-avx2-openmp clang-avx2-openmp-mpi clang-avx2-mpi clang-sse" -EXTRADIRS="g++-avx g++-sse4 icpc-avx icpc-avx2 icpc-avx512" -BLACK="\033[30m" -RED="\033[31m" -GREEN="\033[32m" -YELLOW="\033[33m" -BLUE="\033[34m" -PINK="\033[35m" -CYAN="\033[36m" -WHITE="\033[37m" -NORMAL="\033[0;39m" - -for D in $DIRS -do - -echo -echo -e $RED ============================== -echo -e $GREEN $D -echo -e $RED ============================== -echo -e $BLUE - - cd builds/$D - make clean all -j 8 - cd ../../ -echo -e $NORMAL -done - -if [ "X$1" == "Xextra" ] -then -for D in $EXTRADIRS -do - -echo -echo -e $RED ============================== -echo -e $RED $D -echo -e $RED ============================== -echo -e $BLUE - - cd builds/$D - make clean all -j 8 - cd ../../ -echo -e $NORMAL -done -fi \ No newline at end of file diff --git a/scripts/configure-all b/scripts/configure-all deleted file mode 100755 index ad91034d..00000000 --- a/scripts/configure-all +++ /dev/null @@ -1,11 +0,0 @@ -#!/bin/bash - -DIRS="clang-avx clang-avx-openmp clang-avx-openmp-mpi clang-avx-mpi clang-avx2 clang-avx2-openmp clang-avx2-openmp-mpi clang-avx2-mpi icpc-avx icpc-avx2 icpc-avx512 g++-sse4 g++-avx clang-sse icpc-avx-openmp-mpi icpc-avx-openmp" - -for D in $DIRS -do - mkdir -p builds/$D - cd builds/$D - ../../scripts/configure-commands $D - cd ../.. -done diff --git a/scripts/configure-commands b/scripts/configure-commands deleted file mode 100755 index a3599d1f..00000000 --- a/scripts/configure-commands +++ /dev/null @@ -1,89 +0,0 @@ -#!/bin/bash -WD=$1 -BLACK="\033[30m" -RED="\033[31m" -GREEN="\033[32m" -YELLOW="\033[33m" -BLUE="\033[34m" -PINK="\033[35m" -CYAN="\033[36m" -WHITE="\033[37m" -NORMAL="\033[0;39m" -echo -echo -e $RED ============================== -echo -e $GREEN $WD -echo -e $RED ============================== -echo -e $YELLOW - -case $WD in -g++-avx) - CXX=g++ ../../configure --enable-simd=AVX CXXFLAGS="-mavx -O3 -std=c++11" LIBS="-lgmp -lmpfr" --enable-comms=none - ;; -g++-avx-openmp) - CXX=g++ ../../configure --enable-simd=AVX CXXFLAGS="-mavx -fopenmp -O3 -std=c++11" LIBS="-fopenmp -lgmp -lmpfr" --enable-comms=none - ;; -g++5-sse4) - CXX=g++-5 ../../configure --enable-simd=SSE4 CXXFLAGS="-msse4 -O3 -std=c++11" LIBS="-lgmp -lmpfr" --enable-comms=none - ;; -g++5-avx) - CXX=g++-5 ../../configure --enable-simd=AVX CXXFLAGS="-mavx -O3 -std=c++11" LIBS="-lgmp -lmpfr" --enable-comms=none - ;; -icpc-avx) - CXX=icpc ../../configure --enable-simd=AVX CXXFLAGS="-mavx -O3 -std=c++11" LIBS="-lgmp -lmpfr" --enable-comms=none - ;; -icpc-avx-openmp-mpi) -CXX=icpc ../../configure --enable-simd=AVX CXXFLAGS="-mavx -fopenmp -O3 -I/opt/local/include/openmpi-mp/ -std=c++11" LDFLAGS=-L/opt/local/lib/openmpi-mp/ LIBS="-lmpi -lmpi_cxx -fopenmp -lgmp -lmpfr" --enable-comms=mpi -;; -icpc-avx-openmp) -CXX=icpc ../../configure --enable-precision=single --enable-simd=AVX CXXFLAGS="-mavx -fopenmp -O3 -std=c++11" LIBS="-fopenmp -lgmp -lmpfr" --enable-comms=mpi -;; -icpc-avx2) - CXX=icpc ../../configure --enable-simd=AVX2 CXXFLAGS="-march=core-avx2 -O3 -std=c++11" LIBS="-lgmp -lmpfr" --enable-comms=none - ;; -icpc-avx512) - CXX=icpc ../../configure --enable-simd=AVX512 CXXFLAGS="-xCOMMON-AVX512 -O3 -std=c++11" --host=none LIBS="-lgmp -lmpfr" --enable-comms=none - ;; -icpc-mic) - CXX=icpc ../../configure --host=none --enable-simd=IMCI CXXFLAGS="-mmic -O3 -std=c++11" LDFLAGS=-mmic LIBS="-lgmp -lmpfr" --enable-comms=none - ;; -icpc-mic-avx512) - CXX=icpc ../../configure --host=none --enable-simd=IMCI CXXFLAGS="-xCOMMON_AVX512 -O3 -std=c++11" LDFLAGS=-xCOMMON_AVX512 LIBS="-lgmp -lmpfr" --enable-comms=none - ;; -clang-sse) -CXX=clang++ ../../configure --enable-precision=single --enable-simd=SSE4 CXXFLAGS="-msse4 -O3 -std=c++11" LIBS="-lgmp -lmpfr" --enable-comms=none - ;; -clang-avx) -CXX=clang++ ../../configure --enable-simd=AVX CXXFLAGS="-mavx -O3 -std=c++11" LIBS="-lgmp -lmpfr" --enable-comms=none - ;; -clang-avx2) -CXX=clang++ ../../configure --enable-simd=AVX2 CXXFLAGS="-mavx2 -mfma -O3 -std=c++11" LIBS="-lgmp -lmpfr" --enable-comms=none - ;; -clang-avx-openmp) -CXX=clang-omp++ ../../configure --enable-precision=double --enable-simd=AVX CXXFLAGS="-mavx -fopenmp -O3 -std=c++11" LDFLAGS="-fopenmp" LIBS="-lgmp -lmpfr" --enable-comms=none - ;; -clang-xc30) -CXX=$HOME/Clang/install/bin/clang++ ../../configure --enable-simd=AVX CXXFLAGS="-mavx -O3 -std=c++11 -I/opt/gcc/4.9.2/snos/include/g++/x86_64-suse-linux/ -I/opt/gcc/4.9.2/snos/include/g++/ " LDFLAGS="" LIBS="-lgmp -lmpfr" --enable-comms=none - ;; -clang-xc30-openmp) -CXX=$HOME/Clang/install/bin/clang++ ../../configure --enable-simd=AVX CXXFLAGS="-mavx -fopenmp -O3 -std=c++11 -I/opt/gcc/4.9.2/snos/include/g++/x86_64-suse-linux/ -I/opt/gcc/4.9.2/snos/include/g++/ " LDFLAGS="-fopenmp" LIBS="-lgmp -lmpfr" --enable-comms=none - ;; -clang-avx2-openmp) -CXX=clang-omp++ ../../configure --enable-simd=AVX2 CXXFLAGS="-mavx2 -mfma -fopenmp -O3 -std=c++11" LDFLAGS="-fopenmp" LIBS="-lgmp -lmpfr" --enable-comms=none - ;; -clang-avx-openmp-mpi) -CXX=clang-omp++ ../../configure --enable-simd=AVX CXXFLAGS="-mavx -fopenmp -O3 -I/opt/local/include/openmpi-mp/ -std=c++11" LDFLAGS=-L/opt/local/lib/openmpi-mp/ LIBS="-lmpi -lmpi_cxx -fopenmp -lgmp -lmpfr" --enable-comms=mpi -;; -clang-avx2-openmp-mpi) -CXX=clang-omp++ ../../configure --enable-simd=AVX2 CXXFLAGS="-mavx2 -mfma -fopenmp -O3 -I/opt/local/include/openmpi-mp/ -std=c++11" LDFLAGS=-L/opt/local/lib/openmpi-mp/ LIBS="-lmpi -lmpi_cxx -fopenmp -lgmp -lmpfr" --enable-comms=mpi -;; -clang-avx-mpi) -CXX=clang++ ../../configure --enable-simd=AVX CXXFLAGS="-mavx -O3 -I/opt/local/include/openmpi-mp/ -std=c++11" LDFLAGS=-L/opt/local/lib/openmpi-mp/ LIBS="-lmpi -lmpi_cxx -lgmp -lmpfr" --enable-comms=mpi -;; -clang-avx2-mpi) -CXX=clang++ ../../configure --enable-simd=AVX2 CXXFLAGS="-mavx2 -mfma -O3 -I/opt/local/include/openmpi-mp/ -std=c++11" LDFLAGS=-L/opt/local/lib/openmpi-mp/ LIBS="-lmpi -lmpi_cxx -lgmp -lmpfr" --enable-comms=mpi -;; -clang-avx2) -CXX=clang++ ../../configure --enable-simd=AVX2 CXXFLAGS="-mavx2 -mfma -O3 -std=c++11" LDFLAGS="-L/usr/local/lib/" LIBS="-lgmp -lmpfr" --enable-comms=none -;; -esac -echo -e $NORMAL diff --git a/scripts/configure-cray b/scripts/configure-cray deleted file mode 100755 index db581493..00000000 --- a/scripts/configure-cray +++ /dev/null @@ -1,10 +0,0 @@ -#!/bin/bash -DIRS="g++-avx-openmp g++-avx clang-xc30 clang-xc30-openmp" - -for D in $DIRS -do - mkdir -p builds/$D - cd builds/$D - ../../scripts/configure-commands $D - cd ../.. -done diff --git a/scripts/configure-mic b/scripts/configure-mic deleted file mode 100755 index 668845fe..00000000 --- a/scripts/configure-mic +++ /dev/null @@ -1,10 +0,0 @@ -#!/bin/bash -DIRS="build-icpc-mic" - -for D in $DIRS -do - mkdir -p $D - cd $D - ../configure-commands - cd .. -done diff --git a/scripts/copyright b/scripts/copyright index 92772f16..cc9ed6e5 100755 --- a/scripts/copyright +++ b/scripts/copyright @@ -12,6 +12,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: $1 Copyright (C) 2015 +Copyright (C) 2016 EOF @@ -38,8 +39,21 @@ See the full license in the file "LICENSE" in the top level distribution directo /* END LEGAL */ EOF + cat message > tmp.fil -cat $1 >> tmp.fil + +NOTICE=`grep -n "END LEGAL" $1 | awk '{ print $1 }' ` + +if [ "X$NOTICE" != "X" ] +then + echo "found notice ending on line $NOTICE" + awk 'BEGIN { P=0 } { if ( P ) print } /END LEGAL/{P=1} ' $1 >> tmp.fil +else + cat $1 >> tmp.fil + +fi + + cp tmp.fil $1 shift diff --git a/scripts/cray-modules b/scripts/cray-modules deleted file mode 100644 index 95de2b0b..00000000 --- a/scripts/cray-modules +++ /dev/null @@ -1,2 +0,0 @@ -module swap PrgEnv-cray PrgEnv-intel -module swap intel/14.0.4.211 intel/15.0.2.164 diff --git a/scripts/reconfigure_script b/scripts/reconfigure_script deleted file mode 100755 index d8d7212d..00000000 --- a/scripts/reconfigure_script +++ /dev/null @@ -1,4 +0,0 @@ -aclocal -I m4 -autoheader -f -automake -f --add-missing -autoconf -f diff --git a/scripts/update_fftw.sh b/scripts/update_fftw.sh deleted file mode 100755 index 20e42080..00000000 --- a/scripts/update_fftw.sh +++ /dev/null @@ -1,18 +0,0 @@ -#!/usr/bin/env bash - -if (( $# != 1 )); then - echo "usage: `basename $0` " 1>&2 - exit 1 -fi -ARC=$1 - -INITDIR=`pwd` -rm -rf lib/fftw -mkdir lib/fftw - -ARCDIR=`tar -tf ${ARC} | head -n1 | sed -e 's@/.*@@'` -tar -xf ${ARC} -cp ${ARCDIR}/api/fftw3.h lib/fftw/ - -cd ${INITDIR} -rm -rf ${ARCDIR} diff --git a/scripts/wilson.gnu b/scripts/wilson.gnu deleted file mode 100644 index 69bca5b5..00000000 --- a/scripts/wilson.gnu +++ /dev/null @@ -1,7 +0,0 @@ -plot 'wilson.t1' u 2 w l t "AVX1-OMP=1" -replot 'wilson.t2' u 2 w l t "AVX1-OMP=2" -replot 'wilson.t4' u 2 w l t "AVX1-OMP=4" -set terminal 'pdf' -set output 'wilson_clang.pdf' -replot -quit From 7f0fc0eff5ee11e218462bab1b64a28de71c837e Mon Sep 17 00:00:00 2001 From: James Harrison Date: Tue, 1 Nov 2016 16:02:35 +0000 Subject: [PATCH 49/58] Remove explicit use of double-precision types in photon.h --- lib/qcd/action/gauge/Photon.h | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/lib/qcd/action/gauge/Photon.h b/lib/qcd/action/gauge/Photon.h index 409569d8..852ecb3e 100644 --- a/lib/qcd/action/gauge/Photon.h +++ b/lib/qcd/action/gauge/Photon.h @@ -78,14 +78,14 @@ namespace QCD{ const unsigned int nd = grid->_ndimension; std::vector &l = grid->_fdimensions; std::vector zm(nd,0); - TComplexD Tone = ComplexD(1.0,0.0); - TComplexD Tzero= ComplexD(0.0,0.0); + TComplex Tone = Complex(1.0,0.0); + TComplex Tzero= Complex(0.0,0.0); - one = ComplexD(1.0,0.0); + one = Complex(1.0,0.0); out = zero; for(int mu = 0; mu < nd; mu++) { - RealD twoPiL = M_PI*2./l[mu]; + Real twoPiL = M_PI*2./l[mu]; LatticeCoordinate(kmu,mu); kmu = 2.*sin(.5*twoPiL*kmu); @@ -93,7 +93,7 @@ namespace QCD{ } pokeSite(Tone, out, zm); out = one/out; - pokeSite(Tzero, out,zm); + pokeSite(Tzero, out, zm); } template @@ -107,7 +107,7 @@ namespace QCD{ case ZmScheme::QedTL: { std::vector zm(nd,0); - TComplexD Tzero = ComplexD(0.0,0.0); + TComplex Tzero = Complex(0.0,0.0); pokeSite(Tzero, out, zm); From bb94ddd0ebca02fb7c82b47e1befd411fafe0283 Mon Sep 17 00:00:00 2001 From: paboyle Date: Wed, 2 Nov 2016 08:07:09 +0000 Subject: [PATCH 50/58] Tidy up of mpi3; also some cleaning of the dslash controls. --- benchmarks/Benchmark_comms.cc | 134 +------------- benchmarks/Benchmark_dwf.cc | 63 +++++-- benchmarks/Benchmark_dwf_ntpf.cc | 153 ---------------- benchmarks/Benchmark_dwf_sweep.cc | 16 +- benchmarks/Benchmark_wilson_sweep.cc | 13 ++ benchmarks/Benchmark_zmm.cc | 175 ------------------- lib/Init.cc | 156 +++++++++++------ lib/communicator/Communicator_base.cc | 9 +- lib/communicator/Communicator_base.h | 2 +- lib/communicator/Communicator_mpi3_leader.cc | 147 ++++++++-------- lib/qcd/action/fermion/WilsonKernels.cc | 3 +- lib/qcd/action/fermion/WilsonKernels.h | 80 ++++++--- lib/stencil/Lebesgue.cc | 2 +- 13 files changed, 321 insertions(+), 632 deletions(-) delete mode 100644 benchmarks/Benchmark_dwf_ntpf.cc delete mode 100644 benchmarks/Benchmark_zmm.cc diff --git a/benchmarks/Benchmark_comms.cc b/benchmarks/Benchmark_comms.cc index 234d2fbb..eab9b9b9 100644 --- a/benchmarks/Benchmark_comms.cc +++ b/benchmarks/Benchmark_comms.cc @@ -48,8 +48,8 @@ int main (int argc, char ** argv) std::cout< latt_size ({lat*mpi_layout[0], @@ -124,7 +124,7 @@ int main (int argc, char ** argv) std::cout< latt_size ({lat,lat,lat,lat}); @@ -199,7 +199,7 @@ int main (int argc, char ** argv) std::cout< latt_size ({lat*mpi_layout[0], @@ -271,131 +271,5 @@ int main (int argc, char ** argv) std::cout< latt_size ({lat,lat,lat,lat}); - - GridCartesian Grid(latt_size,simd_layout,mpi_layout); - - std::vector > xbuf(8,std::vector(lat*lat*lat*Ls)); - std::vector > rbuf(8,std::vector(lat*lat*lat*Ls)); - - - int ncomm; - int bytes=lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD); - - - std::vector empty; - std::vector > requests_fwd(Nd,empty); - std::vector > requests_bwd(Nd,empty); - - for(int mu=0;mu<4;mu++){ - ncomm=0; - if (mpi_layout[mu]>1 ) { - ncomm++; - - int comm_proc; - int xmit_to_rank; - int recv_from_rank; - - comm_proc=1; - Grid.ShiftedRanks(mu,comm_proc,xmit_to_rank,recv_from_rank); - Grid.SendToRecvFromInit(requests_fwd[mu], - (void *)&xbuf[mu][0], - xmit_to_rank, - (void *)&rbuf[mu][0], - recv_from_rank, - bytes); - - comm_proc = mpi_layout[mu]-1; - Grid.ShiftedRanks(mu,comm_proc,xmit_to_rank,recv_from_rank); - Grid.SendToRecvFromInit(requests_bwd[mu], - (void *)&xbuf[mu+4][0], - xmit_to_rank, - (void *)&rbuf[mu+4][0], - recv_from_rank, - bytes); - - } - } - - { - double start=usecond(); - for(int i=0;i1 ) { - - Grid.SendToRecvFromBegin(requests_fwd[mu]); - Grid.SendToRecvFromComplete(requests_fwd[mu]); - Grid.SendToRecvFromBegin(requests_bwd[mu]); - Grid.SendToRecvFromComplete(requests_bwd[mu]); - } - } - Grid.Barrier(); - } - - double stop=usecond(); - - double dbytes = bytes; - double xbytes = Nloop*dbytes*2.0*ncomm; - double rbytes = xbytes; - double bidibytes = xbytes+rbytes; - - double time = stop-start; - - std::cout<1 ) { - - Grid.SendToRecvFromBegin(requests_fwd[mu]); - Grid.SendToRecvFromBegin(requests_bwd[mu]); - Grid.SendToRecvFromComplete(requests_fwd[mu]); - Grid.SendToRecvFromComplete(requests_bwd[mu]); - } - } - Grid.Barrier(); - } - - double stop=usecond(); - - double dbytes = bytes; - double xbytes = Nloop*dbytes*2.0*ncomm; - double rbytes = xbytes; - double bidibytes = xbytes+rbytes; - - double time = stop-start; - - std::cout< WilsonFermion5DR; typedef WilsonFermion5D WilsonFermion5DF; typedef WilsonFermion5D WilsonFermion5DD; @@ -54,10 +53,6 @@ int main (int argc, char ** argv) { Grid_init(&argc,&argv); - if( GridCmdOptionExists(argv,argv+argc,"--asynch") ){ - overlapComms = true; - } - int threads = GridThread::GetThreads(); std::cout<_Nprocessors; - for(int doasm=1;doasm<2;doasm++){ - - QCD::WilsonKernelsStatic::AsmOpt=doasm; - DomainWallFermionR Dw(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); - - std::cout<::Dhop "< WilsonFermion5DR; LatticeFermion ssrc(sFGrid); LatticeFermion sref(sFGrid); @@ -248,6 +261,16 @@ int main (int argc, char ** argv) sr_e = zero; sr_o = zero; + std::cout << GridLogMessage<< "*********************************************************" <::DhopEO "< -Author: paboyle - - 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 - -using namespace std; -using namespace Grid; -using namespace Grid::QCD; - -template -struct scal { - d internal; -}; - - Gamma::GammaMatrix Gmu [] = { - Gamma::GammaX, - Gamma::GammaY, - Gamma::GammaZ, - Gamma::GammaT - }; - -bool overlapComms = false; - - -int main (int argc, char ** argv) -{ - Grid_init(&argc,&argv); - - if( GridCmdOptionExists(argv,argv+argc,"--asynch") ){ - overlapComms = true; - } - - int threads = GridThread::GetThreads(); - std::cout< latt4 = GridDefaultLatt(); - const int Ls=16; - GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); - GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); - GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); - GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); - - std::vector seeds4({1,2,3,4}); - std::vector seeds5({5,6,7,8}); - - GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); - GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); - - LatticeFermion src (FGrid); random(RNG5,src); - LatticeFermion result(FGrid); result=zero; - LatticeFermion ref(FGrid); ref=zero; - LatticeFermion tmp(FGrid); - LatticeFermion err(FGrid); - - ColourMatrix cm = Complex(1.0,0.0); - - LatticeGaugeField Umu(UGrid); - random(RNG4,Umu); - - LatticeGaugeField Umu5d(FGrid); - - // replicate across fifth dimension - for(int ss=0;ssoSites();ss++){ - for(int s=0;s U(4,FGrid); - for(int mu=0;mu(Umu5d,mu); - } - - if (1) - { - ref = zero; - for(int mu=0;mu_Nprocessors; - - - QCD::WilsonKernelsStatic::AsmOpt=1; - - DomainWallFermionR Dw(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,params); - - std::cout< seeds({1,2,3,4}); RealD mass = 0.1; + std::cout << GridLogMessage<< "*****************************************************************" < - - 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 - - -using namespace Grid; -using namespace Grid::QCD; - - -int bench(std::ofstream &os, std::vector &latt4,int Ls); - -int main(int argc,char **argv) -{ - Grid_init(&argc,&argv); - std::ofstream os("zmm.dat"); - - os << "#V Ls Lxy Lzt C++ Asm OMP L1 " < grid({L,L,m*L,m*L}); - std::cout << GridLogMessage <<"\t"; - for(int i=0;i<4;i++) { - std::cout << grid[i]<<"x"; - } - std::cout << Ls<<"\t\t"; - bench(os,grid,Ls); - } - } - } -} - -int bench(std::ofstream &os, std::vector &latt4,int Ls) -{ - - GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(latt4, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); - GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); - GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); - GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); - - std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); - std::vector mpi_layout = GridDefaultMpi(); - int threads = GridThread::GetThreads(); - - std::vector seeds4({1,2,3,4}); - std::vector seeds5({5,6,7,8}); - - GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds4); - - LatticeFermion src (FGrid); - LatticeFermion tmp (FGrid); - LatticeFermion srce(FrbGrid); - - LatticeFermion resulto(FrbGrid); resulto=zero; - LatticeFermion resulta(FrbGrid); resulta=zero; - LatticeFermion junk(FrbGrid); junk=zero; - LatticeFermion diff(FrbGrid); - LatticeGaugeField Umu(UGrid); - - double mfc, mfa, mfo, mfl1; - - GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); - GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); - random(RNG5,src); -#if 1 - random(RNG4,Umu); -#else - int mmu=2; - std::vector U(4,UGrid); - for(int mu=0;mu(Umu,mu); - if ( mu!=mmu ) U[mu] = zero; - if ( mu==mmu ) U[mu] = 1.0; - PokeIndex(Umu,U[mu],mu); - } -#endif - pickCheckerboard(Even,srce,src); - - RealD mass=0.1; - RealD M5 =1.8; - DomainWallFermionR Dw(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); - - int ncall=50; - double t0=usecond(); - for(int i=0;i & vec) return; } +void GridCmdOptionInt(std::string &str,int & val) +{ + std::stringstream ss(str); + ss>>val; + return; +} + void GridParseLayout(char **argv,int argc, std::vector &latt, @@ -153,14 +160,12 @@ void GridParseLayout(char **argv,int argc, assert(ompthreads.size()==1); GridThread::SetThreads(ompthreads[0]); } - if( GridCmdOptionExists(argv,argv+argc,"--cores") ){ - std::vector cores(0); + int cores; arg= GridCmdOptionPayload(argv,argv+argc,"--cores"); - GridCmdOptionIntVector(arg,cores); - GridThread::SetCores(cores[0]); + GridCmdOptionInt(arg,cores); + GridThread::SetCores(cores); } - } std::string GridCmdVectorIntToString(const std::vector & vec){ @@ -169,7 +174,7 @@ std::string GridCmdVectorIntToString(const std::vector & vec){ return oss.str(); } ///////////////////////////////////////////////////////// -// +// Reinit guard ///////////////////////////////////////////////////////// static int Grid_is_initialised = 0; @@ -178,27 +183,31 @@ void Grid_init(int *argc,char ***argv) { GridLogger::StopWatch.Start(); + std::string arg; + + //////////////////////////////////// + // Shared memory block size + //////////////////////////////////// + if( GridCmdOptionExists(*argv,*argv+*argc,"--shm") ){ + int MB; + arg= GridCmdOptionPayload(*argv,*argv+*argc,"--shm"); + GridCmdOptionInt(arg,MB); + CartesianCommunicator::MAX_MPI_SHM_BYTES = MB*1024*1024; + } + CartesianCommunicator::Init(argc,argv); - // Parse command line args. + //////////////////////////////////// + // Logging + //////////////////////////////////// - std::string arg; std::vector logstreams; std::string defaultLog("Error,Warning,Message,Performance"); - GridCmdOptionCSL(defaultLog,logstreams); GridLogConfigure(logstreams); - if( GridCmdOptionExists(*argv,*argv+*argc,"--help") ){ - std::cout<= MAX_MPI_SHM_BYTES) { + std::cout<< " ShmBufferMalloc exceeded shared heap size -- try increasing with --shm flag" < /* END LEGAL */ #include "Grid.h" #include +#include namespace Grid { @@ -45,6 +46,7 @@ const int pool = 48; class SlaveState { public: + sem_t sem; volatile int head; volatile int start; volatile int tail; @@ -56,29 +58,32 @@ public: SlaveState *state; MPI_Comm squadron; uint64_t base; + int universe_rank; + int vertical_rank; //////////////////////////////////////////////////////////// // Descriptor circular pointers //////////////////////////////////////////////////////////// Slave() {}; - void Init(SlaveState * _state,MPI_Comm _squadron); + void Init(SlaveState * _state,MPI_Comm _squadron,int _universe_rank,int _vertical_rank); void EventLoop (void) { - std::cerr<< " Entering even loop "<tail != state->head ); - std::cerr<< " FIFO drained "<< state->tail <"<"<"<>0 )&0xFFFF)^((icomm>>16)&0xFFFF) ^ ((icomm>>32)&0xFFFF)^((icomm>>48)&0xFFFF); - hashtag = (comm_hash<<15) | tag; + // hashtag = (comm_hash<<15) | tag; + hashtag = tag; }; -void Slave::Init(SlaveState * _state,MPI_Comm _squadron) +void Slave::Init(SlaveState * _state,MPI_Comm _squadron,int _universe_rank,int _vertical_rank) { squadron=_squadron; + universe_rank=_universe_rank; + vertical_rank=_vertical_rank; state =_state; + std::cout << "state "<<_state<<" comm "<<_squadron<<" universe_rank"<head = state->tail = state->start = 0; - MPI_Barrier(squadron); base = (uint64_t)MPIoffloadEngine::VerticalShmBufs[0]; int rank; MPI_Comm_rank(_squadron,&rank); } #define PERI_PLUS(A) ( (A+1)%pool ) -void Slave::Event (void) { +int Slave::Event (void) { static int tail_last; static int head_last; static int start_last; int ierr; - if ( (state->tail != tail_last) - ||(state->head != head_last) - ||(state->start != start_last) - ) { - std::cerr<< " Event loop "<< state->tail << " "<< state->start<< " "<< state->head <tail; - if ( t != state->start ) { - int flag=0; - - std::cerr<< " Testing tail "<< t<<" "<< (void *)&state->Descrs[t].request - << " "<Descrs[t].request<Descrs[t].request,&flag,MPI_STATUS_IGNORE); - // ierr=MPI_Test((MPI_Request *)&state->Descrs[t].request,&flag,MPI_STATUS_IGNORE); - assert(ierr==0); - if ( flag ) { - state->tail = PERI_PLUS(t); - std::cerr<< " Tail advanced from "<< t<head ) { switch ( state->Descrs[s].command ) { case COMMAND_ISEND: - std::cerr<< " Send "<Descrs[s].buf<< "["<Descrs[s].bytes<<"]" - << " to " << state->Descrs[s].rank<< " tag" << state->Descrs[s].tag - << " Comm " << MPIoffloadEngine::communicator_universe<< std::endl; - - std::cerr<< " Request was "<Descrs[s].request<Descrs[s].buf<< "["<Descrs[s].bytes<<"]" + << " to " << state->Descrs[s].rank<< " tag" << state->Descrs[s].tag + << " Comm " << MPIoffloadEngine::communicator_universe<< " me " <Descrs[s].buf+base), state->Descrs[s].bytes, MPI_CHAR, @@ -438,18 +430,17 @@ void Slave::Event (void) { state->Descrs[s].tag, MPIoffloadEngine::communicator_universe, (MPI_Request *)&state->Descrs[s].request); - std::cerr<< " Request is "<Descrs[s].request<Descrs[0].request<start = PERI_PLUS(s); + return 1; break; case COMMAND_IRECV: - std::cerr<< " Recv "<Descrs[s].buf<< "["<Descrs[s].bytes<<"]" - << " to " << state->Descrs[s].rank<< " tag" << state->Descrs[s].tag - << " Comm " << MPIoffloadEngine::communicator_universe<< std::endl; - - std::cerr<< " Request was "<Descrs[s].request<Descrs[s].buf<< "["<Descrs[s].bytes<<"]" + << " from " << state->Descrs[s].rank<< " tag" << state->Descrs[s].tag + << " Comm " << MPIoffloadEngine::communicator_universe<< " me "<< universe_rank<< std::endl; + */ ierr=MPI_Irecv((void *)(state->Descrs[s].buf+base), state->Descrs[s].bytes, MPI_CHAR, @@ -457,30 +448,32 @@ void Slave::Event (void) { state->Descrs[s].tag, MPIoffloadEngine::communicator_universe, (MPI_Request *)&state->Descrs[s].request); - std::cerr<< " Request is "<Descrs[s].request<Descrs[0].request<Descrs[s].request<Descrs[0].request<start = PERI_PLUS(s); + return 1; break; case COMMAND_WAITALL: - std::cerr<< " Wait all "<tail;t!=s; t=PERI_PLUS(t) ){ - std::cerr<< " Wait ["<Descrs[t].request <Descrs[0].request<Descrs[t].request,MPI_STATUS_IGNORE); }; s=PERI_PLUS(s); state->start = s; state->tail = s; + MPI_Barrier(squadron); + return 1; break; default: assert(0); break; } - return; } + return 0; } ////////////////////////////////////////////////////////////////////////////// // External interaction with the queue @@ -500,17 +493,29 @@ uint64_t Slave::QueueCommand(int command,void *buf, int bytes, int tag, MPI_Comm MPI_Comm communicator; MPI_Request request; - MPIoffloadEngine::MapCommRankToWorldRank(hashtag,commrank,tag,comm,worldrank); + MPIoffloadEngine::MapCommRankToWorldRank(hashtag,worldrank,tag,comm,commrank); - int VerticalRank = MPIoffloadEngine::VerticalRank; uint64_t relative= (uint64_t)buf - base; state->Descrs[head].buf = relative; state->Descrs[head].bytes = bytes; - state->Descrs[head].rank = MPIoffloadEngine::UniverseRanks[worldrank][VerticalRank]; + state->Descrs[head].rank = MPIoffloadEngine::UniverseRanks[worldrank][vertical_rank]; state->Descrs[head].tag = hashtag; state->Descrs[head].command= command; - std::cerr<< " QueueCommand "<tail==next ); @@ -671,6 +676,8 @@ void CartesianCommunicator::StencilSendToRecvFromBegin(std::vector= shm) && (xmit_i+bytes <= shm+MAX_MPI_SHM_BYTES) ); assert( (recv_i >= shm) && (recv_i+bytes <= shm+MAX_MPI_SHM_BYTES) ); + assert(from!=_processor); + assert(dest!=_processor); MPIoffloadEngine::QueueMultiplexedSend(xmit,bytes,_processor,communicator,dest); MPIoffloadEngine::QueueMultiplexedRecv(recv,bytes,from,communicator,from); } diff --git a/lib/qcd/action/fermion/WilsonKernels.cc b/lib/qcd/action/fermion/WilsonKernels.cc index 52ee8704..43776c86 100644 --- a/lib/qcd/action/fermion/WilsonKernels.cc +++ b/lib/qcd/action/fermion/WilsonKernels.cc @@ -32,8 +32,7 @@ directory namespace Grid { namespace QCD { -int WilsonKernelsStatic::HandOpt; -int WilsonKernelsStatic::AsmOpt; +int WilsonKernelsStatic::Opt; template WilsonKernels::WilsonKernels(const ImplParams &p) : Base(p){}; diff --git a/lib/qcd/action/fermion/WilsonKernels.h b/lib/qcd/action/fermion/WilsonKernels.h index 669ee4be..919f7540 100644 --- a/lib/qcd/action/fermion/WilsonKernels.h +++ b/lib/qcd/action/fermion/WilsonKernels.h @@ -40,9 +40,9 @@ namespace QCD { //////////////////////////////////////////////////////////////////////////////////////////////////////////////// class WilsonKernelsStatic { public: + enum { OptGeneric, OptHandUnroll, OptInlineAsm }; // S-direction is INNERMOST and takes no part in the parity. - static int AsmOpt; // these are a temporary hack - static int HandOpt; // these are a temporary hack + static int Opt; // these are a temporary hack }; template class WilsonKernels : public FermionOperator , public WilsonKernelsStatic { @@ -56,24 +56,40 @@ public: template typename std::enable_if::type DiracOptDhopSite(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, SiteHalfSpinor * buf, - int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out) { + int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out) + { + switch(Opt) { #ifdef AVX512 - if (AsmOpt) { - WilsonKernels::DiracOptAsmDhopSite(st,lo,U,buf,sF,sU,Ls,Ns,in,out); - } else { -#else - { -#endif + case OptInlineAsm: for (int site = 0; site < Ns; site++) { for (int s = 0; s < Ls; s++) { - if (HandOpt) - WilsonKernels::DiracOptHandDhopSite(st,lo,U,buf,sF,sU,in,out); - else - WilsonKernels::DiracOptGenericDhopSite(st,lo,U,buf,sF,sU,in,out); + WilsonKernels::DiracOptAsmDhopSite(st,lo,U,buf,sF,sU,Ls,Ns,in,out); sF++; } sU++; } + break; +#endif + case OptHandUnroll: + for (int site = 0; site < Ns; site++) { + for (int s = 0; s < Ls; s++) { + WilsonKernels::DiracOptHandDhopSite(st,lo,U,buf,sF,sU,in,out); + sF++; + } + sU++; + } + break; + case OptGeneric: + for (int site = 0; site < Ns; site++) { + for (int s = 0; s < Ls; s++) { + WilsonKernels::DiracOptGenericDhopSite(st,lo,U,buf,sF,sU,in,out); + sF++; + } + sU++; + } + break; + default: + assert(0); } } @@ -81,7 +97,7 @@ public: typename std::enable_if<(Impl::Dimension != 3 || (Impl::Dimension == 3 && Nc != 3)) && EnableBool, void>::type DiracOptDhopSite(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, SiteHalfSpinor * buf, int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out) { - + // no kernel choice for (int site = 0; site < Ns; site++) { for (int s = 0; s < Ls; s++) { WilsonKernels::DiracOptGenericDhopSite(st, lo, U, buf, sF, sU, in, out); @@ -95,23 +111,39 @@ public: typename std::enable_if::type DiracOptDhopSiteDag(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, SiteHalfSpinor * buf, int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out) { + + switch(Opt) { #ifdef AVX512 - if (AsmOpt) { - WilsonKernels::DiracOptAsmDhopSiteDag(st,lo,U,buf,sF,sU,Ls,Ns,in,out); - } else { -#else - { -#endif + case OptInlineAsm: for (int site = 0; site < Ns; site++) { for (int s = 0; s < Ls; s++) { - if (HandOpt) - WilsonKernels::DiracOptHandDhopSiteDag(st,lo,U,buf,sF,sU,in,out); - else - WilsonKernels::DiracOptGenericDhopSiteDag(st,lo,U,buf,sF,sU,in,out); + WilsonKernels::DiracOptAsmDhopSiteDag(st,lo,U,buf,sF,sU,Ls,Ns,in,out); sF++; } sU++; } + break; +#endif + case OptHandUnroll: + for (int site = 0; site < Ns; site++) { + for (int s = 0; s < Ls; s++) { + WilsonKernels::DiracOptHandDhopSiteDag(st,lo,U,buf,sF,sU,in,out); + sF++; + } + sU++; + } + break; + case OptGeneric: + for (int site = 0; site < Ns; site++) { + for (int s = 0; s < Ls; s++) { + WilsonKernels::DiracOptGenericDhopSiteDag(st,lo,U,buf,sF,sU,in,out); + sF++; + } + sU++; + } + break; + default: + assert(0); } } diff --git a/lib/stencil/Lebesgue.cc b/lib/stencil/Lebesgue.cc index c34b5c96..c30c14c7 100644 --- a/lib/stencil/Lebesgue.cc +++ b/lib/stencil/Lebesgue.cc @@ -32,7 +32,7 @@ Author: paboyle namespace Grid { int LebesgueOrder::UseLebesgueOrder; -std::vector LebesgueOrder::Block({2,2,2,2}); +std::vector LebesgueOrder::Block({8,2,2,2}); LebesgueOrder::IndexInteger LebesgueOrder::alignup(IndexInteger n){ n--; // 1000 0011 --> 1000 0010 From 32375aca657c87674c321ce5e8fc2a82b48dcec2 Mon Sep 17 00:00:00 2001 From: paboyle Date: Wed, 2 Nov 2016 09:27:20 +0000 Subject: [PATCH 51/58] Semaphore sleep/wake up on remote processes. --- lib/communicator/Communicator_mpi3_leader.cc | 84 +++++++++++++++----- 1 file changed, 65 insertions(+), 19 deletions(-) diff --git a/lib/communicator/Communicator_mpi3_leader.cc b/lib/communicator/Communicator_mpi3_leader.cc index db0638da..f160c681 100644 --- a/lib/communicator/Communicator_mpi3_leader.cc +++ b/lib/communicator/Communicator_mpi3_leader.cc @@ -27,7 +27,24 @@ Author: Peter Boyle /* END LEGAL */ #include "Grid.h" #include + +/// bloody mac os doesn't implement unnamed semaphores since it is "optional" posix. +/// darwin dispatch semaphores don't seem to be multiprocess. +//#ifdef __APPLE__ #include +typedef sem_t *Grid_semaphore; +#define SEM_INIT(S) S = sem_open(sem_name,0,00777,0); +#define SEM_INIT_EXCL(S) sem_unlink(sem_name); S = sem_open(sem_name,O_CREAT|O_EXCL,00777,0); +#define SEM_POST(S) sem_post(S); +#define SEM_WAIT(S) sem_wait(S); +//#else +//#include +//typedef sem_t Grid_semaphore; +//#define SEM_INIT(S) +//#define SEM_INIT_EXCL(S) sem_init(&S,0); +//#define SEM_POST(S) sem_post(&S); +//#define SEM_WAIT(S) sem_wait(&S); +//#endif namespace Grid { @@ -46,7 +63,8 @@ const int pool = 48; class SlaveState { public: - sem_t sem; + Grid_semaphore sem_head; + Grid_semaphore sem_tail; volatile int head; volatile int start; volatile int tail; @@ -60,18 +78,43 @@ public: uint64_t base; int universe_rank; int vertical_rank; + char sem_name [NAME_MAX]; //////////////////////////////////////////////////////////// // Descriptor circular pointers //////////////////////////////////////////////////////////// Slave() {}; void Init(SlaveState * _state,MPI_Comm _squadron,int _universe_rank,int _vertical_rank); - + + void SemInit(void) { + sprintf(sem_name,"/Grid_mpi3_head_%d",universe_rank); + SEM_INIT(state->sem_head); + sprintf(sem_name,"/Grid_mpi3_tail_%d",universe_rank); + SEM_INIT(state->sem_tail); + } + void SemInitExcl(void) { + sprintf(sem_name,"/Grid_mpi3_head_%d",universe_rank); + SEM_INIT_EXCL(state->sem_head); + sprintf(sem_name,"/Grid_mpi3_tail_%d",universe_rank); + SEM_INIT_EXCL(state->sem_tail); + } + void WakeUpDMA(void) { + SEM_POST(state->sem_head); + }; + void WakeUpCompute(void) { + SEM_POST(state->sem_tail); + }; + void WaitForCommand(void) { + SEM_WAIT(state->sem_head); + }; + void WaitForComplete(void) { + SEM_WAIT(state->sem_tail); + }; void EventLoop (void) { std::cout<< " Entering event loop "<tail != state->head ); + QueueCommand(COMMAND_WAITALL,0,0,0,squadron,0); + WakeUpDMA(); + WaitForComplete(); + assert ( state->tail == state->head ); } }; @@ -130,25 +173,22 @@ public: // routines for master proc must handle any communicator ///////////////////////////////////////////////////////// - static uint64_t QueueSend(int slave,void *buf, int bytes, int tag, MPI_Comm comm,int rank) { + static void QueueSend(int slave,void *buf, int bytes, int tag, MPI_Comm comm,int rank) { // std::cout<< " Queueing send "<< bytes<< " slave "<< slave << " to comm "<start = s; state->tail = s; - MPI_Barrier(squadron); + + WakeUpCompute(); + return 1; break; From 757a928f9a70c0f2448cc6844f1112bab9cdf7d9 Mon Sep 17 00:00:00 2001 From: paboyle Date: Wed, 2 Nov 2016 12:37:46 +0000 Subject: [PATCH 52/58] Improvement to use own SHM_OPEN call to avoid openmpi bug. --- benchmarks/Benchmark_comms.cc | 4 +- lib/communicator/Communicator_mpi3_leader.cc | 159 ++++++++++++++----- 2 files changed, 124 insertions(+), 39 deletions(-) diff --git a/benchmarks/Benchmark_comms.cc b/benchmarks/Benchmark_comms.cc index eab9b9b9..de73bc81 100644 --- a/benchmarks/Benchmark_comms.cc +++ b/benchmarks/Benchmark_comms.cc @@ -48,7 +48,7 @@ int main (int argc, char ** argv) std::cout< #include "Grid.h" #include -/// bloody mac os doesn't implement unnamed semaphores since it is "optional" posix. -/// darwin dispatch semaphores don't seem to be multiprocess. -//#ifdef __APPLE__ +//////////////////////////////////////////////////////////////////////////////////////////////////////////////// +/// Workarounds: +/// i) bloody mac os doesn't implement unnamed semaphores since it is "optional" posix. +/// darwin dispatch semaphores don't seem to be multiprocess. +/// +/// ii) openmpi under --mca shmem posix works with two squadrons per node; +/// openmpi under default mca settings (I think --mca shmem mmap) on MacOS makes two squadrons map the SAME +/// memory as each other, despite their living on different communicators. This appears to be a bug in OpenMPI. +/// +//////////////////////////////////////////////////////////////////////////////////////////////////////////////// #include typedef sem_t *Grid_semaphore; -#define SEM_INIT(S) S = sem_open(sem_name,0,00777,0); -#define SEM_INIT_EXCL(S) sem_unlink(sem_name); S = sem_open(sem_name,O_CREAT|O_EXCL,00777,0); -#define SEM_POST(S) sem_post(S); -#define SEM_WAIT(S) sem_wait(S); -//#else -//#include -//typedef sem_t Grid_semaphore; -//#define SEM_INIT(S) -//#define SEM_INIT_EXCL(S) sem_init(&S,0); -//#define SEM_POST(S) sem_post(&S); -//#define SEM_WAIT(S) sem_wait(&S); -//#endif + +#define SEM_INIT(S) S = sem_open(sem_name,0,0600,0); assert ( S != SEM_FAILED ); +#define SEM_INIT_EXCL(S) sem_unlink(sem_name); S = sem_open(sem_name,O_CREAT|O_EXCL,0600,0); assert ( S != SEM_FAILED ); +#define SEM_POST(S) assert ( sem_post(S) == 0 ); +#define SEM_WAIT(S) assert ( sem_wait(S) == 0 ); + +#include + namespace Grid { @@ -63,8 +66,6 @@ const int pool = 48; class SlaveState { public: - Grid_semaphore sem_head; - Grid_semaphore sem_tail; volatile int head; volatile int start; volatile int tail; @@ -73,6 +74,8 @@ public: class Slave { public: + Grid_semaphore sem_head; + Grid_semaphore sem_tail; SlaveState *state; MPI_Comm squadron; uint64_t base; @@ -87,33 +90,38 @@ public: void Init(SlaveState * _state,MPI_Comm _squadron,int _universe_rank,int _vertical_rank); void SemInit(void) { - sprintf(sem_name,"/Grid_mpi3_head_%d",universe_rank); - SEM_INIT(state->sem_head); - sprintf(sem_name,"/Grid_mpi3_tail_%d",universe_rank); - SEM_INIT(state->sem_tail); + sprintf(sem_name,"/Grid_mpi3_sem_head_%d",universe_rank); + printf("SEM_NAME: %s \n",sem_name); + SEM_INIT(sem_head); + sprintf(sem_name,"/Grid_mpi3_sem_tail_%d",universe_rank); + printf("SEM_NAME: %s \n",sem_name); + SEM_INIT(sem_tail); } void SemInitExcl(void) { - sprintf(sem_name,"/Grid_mpi3_head_%d",universe_rank); - SEM_INIT_EXCL(state->sem_head); - sprintf(sem_name,"/Grid_mpi3_tail_%d",universe_rank); - SEM_INIT_EXCL(state->sem_tail); + sprintf(sem_name,"/Grid_mpi3_sem_head_%d",universe_rank); + printf("SEM_INIT_EXCL: %s \n",sem_name); + SEM_INIT_EXCL(sem_head); + sprintf(sem_name,"/Grid_mpi3_sem_tail_%d",universe_rank); + printf("SEM_INIT_EXCL: %s \n",sem_name); + SEM_INIT_EXCL(sem_tail); } void WakeUpDMA(void) { - SEM_POST(state->sem_head); + SEM_POST(sem_head); }; void WakeUpCompute(void) { - SEM_POST(state->sem_tail); + SEM_POST(sem_tail); }; void WaitForCommand(void) { - SEM_WAIT(state->sem_head); + SEM_WAIT(sem_head); }; void WaitForComplete(void) { - SEM_WAIT(state->sem_tail); + SEM_WAIT(sem_tail); }; void EventLoop (void) { std::cout<< " Entering event loop "<tail == state->head ); } }; @@ -176,19 +188,25 @@ public: static void QueueSend(int slave,void *buf, int bytes, int tag, MPI_Comm comm,int rank) { // std::cout<< " Queueing send "<< bytes<< " slave "<< slave << " to comm "< MPIoffloadEngine::VerticalShmBufs; std::vector > MPIoffloadEngine::UniverseRanks; std::vector MPIoffloadEngine::UserCommunicatorToWorldRanks; @@ -274,8 +291,12 @@ void MPIoffloadEngine::CommunicatorInit (MPI_Comm &communicator_world, ///////////////////////////////////////////////////////////////////// // Split into groups that can share memory (Verticals) ///////////////////////////////////////////////////////////////////// - // MPI_Comm_split_type(communicator_universe, MPI_COMM_TYPE_SHARED, 0, MPI_INFO_NULL,&VerticalComm); - MPI_Comm_split(communicator_universe,(UniverseRank/2),UniverseRank,&VerticalComm); +#define MPI_SHARED_MEM_DEBUG +#ifdef MPI_SHARED_MEM_DEBUG + MPI_Comm_split(communicator_universe,(UniverseRank/4),UniverseRank,&VerticalComm); +#else + MPI_Comm_split_type(communicator_universe, MPI_COMM_TYPE_SHARED, 0, MPI_INFO_NULL,&VerticalComm); +#endif MPI_Comm_rank(VerticalComm ,&VerticalRank); MPI_Comm_size(VerticalComm ,&VerticalSize); @@ -308,19 +329,83 @@ void MPIoffloadEngine::CommunicatorInit (MPI_Comm &communicator_world, ////////////////////////////////////////////////////////////////////////////////////////////////////////// // allocate the shared window for our group, pass back Shm info to CartesianCommunicator ////////////////////////////////////////////////////////////////////////////////////////////////////////// + VerticalShmBufs.resize(VerticalSize); + +#undef MPI_SHARED_MEM +#ifdef MPI_SHARED_MEM ierr = MPI_Win_allocate_shared(CartesianCommunicator::MAX_MPI_SHM_BYTES,1,MPI_INFO_NULL,VerticalComm,&ShmCommBuf,&VerticalWindow); ierr|= MPI_Win_lock_all (MPI_MODE_NOCHECK, VerticalWindow); assert(ierr==0); - - std::cout<<"SHM "<0 ) size = sizeof(SlaveState); + + sprintf(shm_name,"/Grid_mpi3_shm_%d_%d",WorldRank,r); + + shm_unlink(shm_name); + + int fd=shm_open(shm_name,O_RDWR|O_CREAT,0600); + if ( fd < 0 ) { + perror("failed shm_open"); + assert(0); + } + + ftruncate(fd, size); + + VerticalShmBufs[r] = mmap(NULL,size, PROT_READ | PROT_WRITE, MAP_SHARED, fd, 0); + + if ( VerticalShmBufs[r] == MAP_FAILED ) { + perror("failed mmap"); + assert(0); + } + + uint64_t * check = (uint64_t *) VerticalShmBufs[r]; + check[0] = WorldRank; + check[1] = r; + + // std::cout<<"SHM "<0 ) size = sizeof(SlaveState); + + sprintf(shm_name,"/Grid_mpi3_shm_%d_%d",WorldRank,r); + + int fd=shm_open(shm_name,O_RDWR|O_CREAT,0600); + if ( fd<0 ) { + perror("failed shm_open"); + assert(0); + } + VerticalShmBufs[r] = mmap(NULL,size, PROT_READ | PROT_WRITE, MAP_SHARED, fd, 0); + + uint64_t * check = (uint64_t *) VerticalShmBufs[r]; + assert(check[0]== WorldRank); + assert(check[1]== r); + std::cerr<<"SHM "< Date: Wed, 2 Nov 2016 19:54:03 +0000 Subject: [PATCH 53/58] Decrease mpi3l verbose --- lib/communicator/Communicator_mpi3_leader.cc | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/lib/communicator/Communicator_mpi3_leader.cc b/lib/communicator/Communicator_mpi3_leader.cc index b6995a44..49a3e948 100644 --- a/lib/communicator/Communicator_mpi3_leader.cc +++ b/lib/communicator/Communicator_mpi3_leader.cc @@ -48,7 +48,6 @@ typedef sem_t *Grid_semaphore; #include - namespace Grid { enum { COMMAND_ISEND, COMMAND_IRECV, COMMAND_WAITALL }; @@ -91,18 +90,18 @@ public: void SemInit(void) { sprintf(sem_name,"/Grid_mpi3_sem_head_%d",universe_rank); - printf("SEM_NAME: %s \n",sem_name); + // printf("SEM_NAME: %s \n",sem_name); SEM_INIT(sem_head); sprintf(sem_name,"/Grid_mpi3_sem_tail_%d",universe_rank); - printf("SEM_NAME: %s \n",sem_name); + // printf("SEM_NAME: %s \n",sem_name); SEM_INIT(sem_tail); } void SemInitExcl(void) { sprintf(sem_name,"/Grid_mpi3_sem_head_%d",universe_rank); - printf("SEM_INIT_EXCL: %s \n",sem_name); + // printf("SEM_INIT_EXCL: %s \n",sem_name); SEM_INIT_EXCL(sem_head); sprintf(sem_name,"/Grid_mpi3_sem_tail_%d",universe_rank); - printf("SEM_INIT_EXCL: %s \n",sem_name); + // printf("SEM_INIT_EXCL: %s \n",sem_name); SEM_INIT_EXCL(sem_tail); } void WakeUpDMA(void) { @@ -118,7 +117,7 @@ public: SEM_WAIT(sem_tail); }; void EventLoop (void) { - std::cout<< " Entering event loop "<head = state->tail = state->start = 0; base = (uint64_t)MPIoffloadEngine::VerticalShmBufs[0]; int rank; MPI_Comm_rank(_squadron,&rank); From 111bfbc6bc60b4edcb696f21c2edff4606fdb5f0 Mon Sep 17 00:00:00 2001 From: paboyle Date: Thu, 3 Nov 2016 11:40:26 +0000 Subject: [PATCH 54/58] notimestamp by default --- lib/Init.cc | 6 ++++-- tests/forces/Test_wilson_force.cc | 18 ++++++++---------- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/lib/Init.cc b/lib/Init.cc index d15e4bd1..695bb01b 100644 --- a/lib/Init.cc +++ b/lib/Init.cc @@ -234,7 +234,7 @@ void Grid_init(int *argc,char ***argv) std::cout<(mom,mommu,mu); // fourth order exponential approx - parallel_for(auto i=mom.begin();i Date: Thu, 3 Nov 2016 13:48:07 +0000 Subject: [PATCH 55/58] MPI auto configure fix --- configure.ac | 31 ++++++++++++++++++------------- 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/configure.ac b/configure.ac index a8ad40dd..9178a11d 100644 --- a/configure.ac +++ b/configure.ac @@ -253,18 +253,23 @@ AC_ARG_ENABLE([comms],[AC_HELP_STRING([--enable-comms=none|mpi|mpi-auto|mpi3|mpi case ${ac_COMMS} in none) AC_DEFINE([GRID_COMMS_NONE],[1],[GRID_COMMS_NONE] ) + comms_type='none' ;; - mpi|mpi-auto) - AC_DEFINE([GRID_COMMS_MPI],[1],[GRID_COMMS_MPI] ) - ;; - mpi3|mpi3-auto) - AC_DEFINE([GRID_COMMS_MPI3],[1],[GRID_COMMS_MPI3] ) - ;; - mpi3l) + mpi3l*) AC_DEFINE([GRID_COMMS_MPI3L],[1],[GRID_COMMS_MPI3L] ) + comms_type='mpi3l' + ;; + mpi3*) + AC_DEFINE([GRID_COMMS_MPI3],[1],[GRID_COMMS_MPI3] ) + comms_type='mpi3' + ;; + mpi*) + AC_DEFINE([GRID_COMMS_MPI],[1],[GRID_COMMS_MPI] ) + comms_type='mpi' ;; shmem) AC_DEFINE([GRID_COMMS_SHMEM],[1],[GRID_COMMS_SHMEM] ) + comms_type='shmem' ;; *) AC_MSG_ERROR([${ac_COMMS} unsupported --enable-comms option]); @@ -282,11 +287,11 @@ case ${ac_COMMS} in ;; esac -AM_CONDITIONAL(BUILD_COMMS_SHMEM,[ test "X${ac_COMMS}X" == "XshmemX" ]) -AM_CONDITIONAL(BUILD_COMMS_MPI,[ test "X${ac_COMMS}X" == "XmpiX" || test "X${ac_COMMS}X" == "Xmpi-autoX" ]) -AM_CONDITIONAL(BUILD_COMMS_MPI3,[ test "X${ac_COMMS}X" == "Xmpi3X"] ) -AM_CONDITIONAL(BUILD_COMMS_MPI3L,[ test "X${ac_COMMS}X" == "Xmpi3lX"] ) -AM_CONDITIONAL(BUILD_COMMS_NONE,[ test "X${ac_COMMS}X" == "XnoneX" ]) +AM_CONDITIONAL(BUILD_COMMS_SHMEM, [ test "${comms_type}X" == "shmemX" ]) +AM_CONDITIONAL(BUILD_COMMS_MPI, [ test "${comms_type}X" == "mpiX" ]) +AM_CONDITIONAL(BUILD_COMMS_MPI3, [ test "${comms_type}X" == "mpi3X" ] ) +AM_CONDITIONAL(BUILD_COMMS_MPI3L, [ test "${comms_type}X" == "mpi3lX" ] ) +AM_CONDITIONAL(BUILD_COMMS_NONE, [ test "${comms_type}X" == "noneX" ]) ############### RNG selection AC_ARG_ENABLE([rng],[AC_HELP_STRING([--enable-rng=ranlux48|mt19937],\ @@ -379,7 +384,7 @@ compiler version : ${ax_cv_gxx_version} ----- BUILD OPTIONS ----------------------------------- SIMD : ${ac_SIMD} Threading : ${ac_openmp} -Communications type : ${ac_COMMS} +Communications type : ${comms_type} Default precision : ${ac_PRECISION} RNG choice : ${ac_RNG} GMP : `if test "x$have_gmp" = xtrue; then echo yes; else echo no; fi` From c65d23935a4e30d8c3a501772c5afd270a27fccc Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 3 Nov 2016 13:48:20 +0000 Subject: [PATCH 56/58] README update --- README.md | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/README.md b/README.md index faf86faf..c7461368 100644 --- a/README.md +++ b/README.md @@ -20,7 +20,7 @@ License: GPL v2. Last update Nov 2016. -_Please send all pull requests to the `develop` branch._ +_Please do not send pull requests to the `master` branch which is reserved for releases._ ### Bug report @@ -29,7 +29,7 @@ _To help us tracking and solving more efficiently issues with Grid, please repor When you file an issue, please go though the following checklist: 1. Check that the code is pointing to the `HEAD` of `develop` or any commit in `master` which is tagged with a version number. -2. Give a description of the target platform (CPU, network, compiler). +2. Give a description of the target platform (CPU, network, compiler). Please give the full CPU part description, using for example `cat /proc/cpuinfo | grep 'model name' | uniq` (Linux) or `sysctl machdep.cpu.brand_string` (macOS) and the full output the `--version` option of your compiler. 3. Give the exact `configure` command used. 4. Attach `config.log`. 5. Attach `config.summary`. @@ -45,7 +45,7 @@ are provided, similar to HPF and cmfortran, and user control is given over the m array indices to both MPI tasks and SIMD processing elements. * Identically shaped arrays then be processed with perfect data parallelisation. -* Such identically shapped arrays are called conformable arrays. +* Such identically shaped arrays are called conformable arrays. The transformation is based on the observation that Cartesian array processing involves identical processing to be performed on different regions of the Cartesian array. @@ -127,14 +127,15 @@ make -C tests/ tests The following options can be use with the `--enable-simd=` option to target different communication interfaces: -| `` | Description | -| ------------- | -------------------------------------------- | -| `none` | no communications | -| `mpi[-auto]` | MPI communications | -| `mpi3[-auto]` | MPI communications using MPI 3 shared memory | -| `shmem ` | Cray SHMEM communications | +| `` | Description | +| -------------- | ------------------------------------------------------------- | +| `none` | no communications | +| `mpi[-auto]` | MPI communications | +| `mpi3[-auto]` | MPI communications using MPI 3 shared memory | +| `mpi3l[-auto]` | MPI communications using MPI 3 shared memory and leader model | +| `shmem ` | Cray SHMEM communications | -For `mpi` and `mpi3` the optional `-auto` suffix instructs the `configure` scripts to determine all the necessary compilation and linking flags. This is done by extracting the informations from the MPI wrapper specified in the environment variable `MPICXX` (if not specified `configure` will scan though a list of default names). +For the MPI interfaces the optional `-auto` suffix instructs the `configure` scripts to determine all the necessary compilation and linking flags. This is done by extracting the informations from the MPI wrapper specified in the environment variable `MPICXX` (if not specified `configure` will scan though a list of default names). ### Possible SIMD types @@ -160,7 +161,7 @@ Alternatively, some CPU codenames can be directly used: | `BGQ` | Blue Gene/Q | #### Notes: -- We currently support AVX512 only for the Intel compiler. Support for GCC and clang will appear in future versions. +- We currently support AVX512 only for the Intel compiler. Support for GCC and clang will appear in future versions of Grid when the AVX512 support within GCC and clang will be more advanced. - For BG/Q only [bgclang](http://trac.alcf.anl.gov/projects/llvm-bgq) is supported. We do not presently plan to support more compilers for this platform. - BG/Q performances are currently rather poor. This is being investigated for future versions. @@ -171,7 +172,7 @@ The following configuration is recommended for the Intel Knights Landing platfor ``` bash ../configure --enable-precision=double\ --enable-simd=KNL \ - --enable-comms=mpi3-auto \ + --enable-comms=mpi-auto \ --with-gmp= \ --with-mpfr= \ --enable-mkl \ @@ -183,10 +184,9 @@ where `` is the UNIX prefix where GMP and MPFR are installed. If you are w ``` bash ../configure --enable-precision=double\ --enable-simd=KNL \ - --enable-comms=mpi3 \ + --enable-comms=mpi \ --with-gmp= \ --with-mpfr= \ --enable-mkl \ CXX=CC CC=cc -``` - +``` \ No newline at end of file From aee44dc69472e610cf2b3041ac203e99459d5276 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 3 Nov 2016 13:54:15 +0000 Subject: [PATCH 57/58] Photon.h removed from develop branch --- lib/qcd/action/gauge/Photon.h | 240 ---------------------------------- 1 file changed, 240 deletions(-) delete mode 100644 lib/qcd/action/gauge/Photon.h diff --git a/lib/qcd/action/gauge/Photon.h b/lib/qcd/action/gauge/Photon.h deleted file mode 100644 index 852ecb3e..00000000 --- a/lib/qcd/action/gauge/Photon.h +++ /dev/null @@ -1,240 +0,0 @@ -/************************************************************************************* - - Grid physics library, www.github.com/paboyle/Grid - - Source file: ./lib/qcd/action/gauge/Photon.h - - Copyright (C) 2015 - - Author: Peter Boyle - - 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 QCD_PHOTON_ACTION_H -#define QCD_PHOTON_ACTION_H - -namespace Grid{ -namespace QCD{ - - template - class Photon - { - public: - INHERIT_GIMPL_TYPES(Gimpl); - enum class Gauge {Feynman, Coulomb, Landau}; - enum class ZmScheme {QedL, QedTL}; - public: - Photon(Gauge gauge, ZmScheme zmScheme); - virtual ~Photon(void) = default; - void FreePropagator(const GaugeField &in, GaugeField &out); - void MomentumSpacePropagator(const GaugeField &in, GaugeField &out); - void StochasticField(GaugeField &out, GridParallelRNG &rng); - private: - void invKHatSquared(GaugeLinkField &out); - void zmSub(GaugeLinkField &out); - private: - Gauge gauge_; - ZmScheme zmScheme_; - }; - - template - Photon::Photon(Gauge gauge, ZmScheme zmScheme) - : gauge_(gauge), zmScheme_(zmScheme) - {} - - template - void Photon::FreePropagator (const GaugeField &in,GaugeField &out) - { - FFT theFFT(in._grid); - - GaugeField in_k(in._grid); - GaugeField prop_k(in._grid); - - theFFT.FFT_all_dim(in_k,in,FFT::forward); - MomentumSpacePropagator(prop_k,in_k); - theFFT.FFT_all_dim(out,prop_k,FFT::backward); - } - - template - void Photon::invKHatSquared(GaugeLinkField &out) - { - GridBase *grid = out._grid; - GaugeLinkField kmu(grid), one(grid); - const unsigned int nd = grid->_ndimension; - std::vector &l = grid->_fdimensions; - std::vector zm(nd,0); - TComplex Tone = Complex(1.0,0.0); - TComplex Tzero= Complex(0.0,0.0); - - one = Complex(1.0,0.0); - out = zero; - for(int mu = 0; mu < nd; mu++) - { - Real twoPiL = M_PI*2./l[mu]; - - LatticeCoordinate(kmu,mu); - kmu = 2.*sin(.5*twoPiL*kmu); - out = out + kmu*kmu; - } - pokeSite(Tone, out, zm); - out = one/out; - pokeSite(Tzero, out, zm); - } - - template - void Photon::zmSub(GaugeLinkField &out) - { - GridBase *grid = out._grid; - const unsigned int nd = grid->_ndimension; - - switch (zmScheme_) - { - case ZmScheme::QedTL: - { - std::vector zm(nd,0); - TComplex Tzero = Complex(0.0,0.0); - - pokeSite(Tzero, out, zm); - - break; - } - case ZmScheme::QedL: - { - LatticeInteger spNrm(grid), coor(grid); - GaugeLinkField z(grid); - - spNrm = zero; - for(int d = 0; d < grid->_ndimension - 1; d++) - { - LatticeCoordinate(coor,d); - spNrm = spNrm + coor*coor; - } - out = where(spNrm == Integer(0), 0.*out, out); - - break; - } - default: - break; - } - } - - template - void Photon::MomentumSpacePropagator(const GaugeField &in, - GaugeField &out) - { - GridBase *grid = out._grid; - LatticeComplex k2Inv(grid); - - invKHatSquared(k2Inv); - zmSub(k2Inv); - - out = in*k2Inv; - } - - template - void Photon::StochasticField(GaugeField &out, GridParallelRNG &rng) - { - auto *grid = dynamic_cast(out._grid); - const unsigned int nd = grid->_ndimension; - std::vector latt_size = grid->_fdimensions; - GaugeLinkField sqrtK2Inv(grid), r(grid); - GaugeField aTilde(grid); - FFT fft(grid); - - Integer vol = 1; - for(int d = 0; d < nd; d++) - { - vol = vol * latt_size[d]; - } - - invKHatSquared(sqrtK2Inv); - sqrtK2Inv = sqrt(vol*real(sqrtK2Inv)); - zmSub(sqrtK2Inv); - for(int mu = 0; mu < nd; mu++) - { - gaussian(rng, r); - r = sqrtK2Inv*r; - pokeLorentz(aTilde, r, mu); - } - fft.FFT_all_dim(out, aTilde, FFT::backward); - } -// template -// void Photon::FeynmanGaugeMomentumSpacePropagator_L(GaugeField &out, -// const GaugeField &in) -// { -// -// FeynmanGaugeMomentumSpacePropagator_TL(out,in); -// -// GridBase *grid = out._grid; -// LatticeInteger coor(grid); -// GaugeField zz(grid); zz=zero; -// -// // xyzt -// for(int d = 0; d < grid->_ndimension-1;d++){ -// LatticeCoordinate(coor,d); -// out = where(coor==Integer(0),zz,out); -// } -// } -// -// template -// void Photon::FeynmanGaugeMomentumSpacePropagator_TL(GaugeField &out, -// const GaugeField &in) -// { -// -// // what type LatticeComplex -// GridBase *grid = out._grid; -// int nd = grid->_ndimension; -// -// typedef typename GaugeField::vector_type vector_type; -// typedef typename GaugeField::scalar_type ScalComplex; -// typedef Lattice > LatComplex; -// -// std::vector latt_size = grid->_fdimensions; -// -// LatComplex denom(grid); denom= zero; -// LatComplex one(grid); one = ScalComplex(1.0,0.0); -// LatComplex kmu(grid); -// -// ScalComplex ci(0.0,1.0); -// // momphase = n * 2pi / L -// for(int mu=0;mu zero_mode(nd,0); -// TComplexD Tone = ComplexD(1.0,0.0); -// TComplexD Tzero= ComplexD(0.0,0.0); -// -// pokeSite(Tone,denom,zero_mode); -// -// denom= one/denom; -// -// pokeSite(Tzero,denom,zero_mode); -// -// out = zero; -// out = in*denom; -// }; - -}} -#endif From 2854e601e6b1720da5b0ae3dbc9bec9906667bb4 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 3 Nov 2016 14:09:47 +0000 Subject: [PATCH 58/58] FFT test typo --- tests/core/Test_fft.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/core/Test_fft.cc b/tests/core/Test_fft.cc index 6872f0e1..93596066 100644 --- a/tests/core/Test_fft.cc +++ b/tests/core/Test_fft.cc @@ -70,7 +70,7 @@ int main (int argc, char ** argv) std::cout<<"*************************************************"<