From 66d997e03102bf52f6310c0f285be99820d34868 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sun, 31 May 2015 15:09:02 +0100 Subject: [PATCH] Large scale change to support 5d fermion formulations. Have 5d replicated wilson with 4d gauge working and matrix regressing to Ls copies of wilson. --- benchmarks/Grid_comms.cc | 2 +- benchmarks/Grid_dwf.cc | 145 +++++++++++ benchmarks/Grid_wilson.cc | 9 +- benchmarks/Grid_wilson_cg_prec.cc | 6 +- benchmarks/Grid_wilson_cg_schur.cc | 2 +- benchmarks/Grid_wilson_cg_unprec.cc | 4 +- benchmarks/Grid_wilson_evenodd.cc | 2 +- benchmarks/Makefile.am | 5 +- lib/Grid.h | 2 +- lib/Grid_init.cc | 6 +- lib/Grid_stencil.h | 33 +-- lib/Grid_threads.h | 2 +- lib/Makefile.am | 166 +++++++------ lib/algorithms/SparseMatrix.h | 7 +- lib/algorithms/iterative/SchurRedBlack.h | 4 +- lib/cartesian/Grid_cartesian_base.h | 11 +- lib/cartesian/Grid_cartesian_full.h | 7 +- lib/cartesian/Grid_cartesian_red_black.h | 196 ++++++++++----- lib/cshift/Grid_cshift_common.h | 14 +- lib/cshift/Grid_cshift_mpi.h | 16 +- lib/cshift/Grid_cshift_none.h | 2 +- lib/lattice/Grid_lattice_peekpoke.h | 2 +- lib/lattice/Grid_lattice_transfer.h | 3 +- lib/qcd/{Grid_qcd_dirac.cc => Dirac.cc} | 0 lib/qcd/{Grid_qcd_dirac.h => Dirac.h} | 0 lib/qcd/Grid_qcd_wilson_dop.cc | 217 ----------------- lib/qcd/Grid_qcd_wilson_dop.h | 105 -------- lib/qcd/{Grid_qcd.h => QCD.h} | 17 +- lib/qcd/{Grid_qcd_2spinor.h => TwoSpinor.h} | 0 lib/qcd/action/Actions.h | 10 + lib/qcd/action/fermion/FermionAction.h | 47 ++++ .../action/fermion/FiveDimWilsonFermion.cc | 228 ++++++++++++++++++ lib/qcd/action/fermion/FiveDimWilsonFermion.h | 108 +++++++++ lib/qcd/action/fermion/WilsonCompressor.h | 61 +++++ lib/qcd/action/fermion/WilsonFermion.cc | 163 +++++++++++++ lib/qcd/action/fermion/WilsonFermion.h | 87 +++++++ .../fermion/WilsonKernels.cc} | 46 ++-- lib/qcd/action/fermion/WilsonKernels.h | 42 ++++ .../fermion/WilsonKernelsHand.cc} | 57 ++--- lib/stencil/Grid_lebesgue.cc | 103 ++++++++ lib/stencil/Grid_lebesgue.h | 29 +++ lib/stencil/Grid_stencil_common.cc | 103 +------- scripts/linecount | 2 +- tests/Grid_cshift_red_black.cc | 165 +++++++++++++ tests/Makefile.am | 5 +- 45 files changed, 1549 insertions(+), 692 deletions(-) create mode 100644 benchmarks/Grid_dwf.cc rename lib/qcd/{Grid_qcd_dirac.cc => Dirac.cc} (100%) rename lib/qcd/{Grid_qcd_dirac.h => Dirac.h} (100%) delete mode 100644 lib/qcd/Grid_qcd_wilson_dop.cc delete mode 100644 lib/qcd/Grid_qcd_wilson_dop.h rename lib/qcd/{Grid_qcd.h => QCD.h} (97%) rename lib/qcd/{Grid_qcd_2spinor.h => TwoSpinor.h} (100%) create mode 100644 lib/qcd/action/Actions.h create mode 100644 lib/qcd/action/fermion/FermionAction.h create mode 100644 lib/qcd/action/fermion/FiveDimWilsonFermion.cc create mode 100644 lib/qcd/action/fermion/FiveDimWilsonFermion.h create mode 100644 lib/qcd/action/fermion/WilsonCompressor.h create mode 100644 lib/qcd/action/fermion/WilsonFermion.cc create mode 100644 lib/qcd/action/fermion/WilsonFermion.h rename lib/qcd/{Grid_qcd_dhop.cc => action/fermion/WilsonKernels.cc} (87%) create mode 100644 lib/qcd/action/fermion/WilsonKernels.h rename lib/qcd/{Grid_qcd_dhop_hand.cc => action/fermion/WilsonKernelsHand.cc} (92%) create mode 100644 lib/stencil/Grid_lebesgue.cc create mode 100644 lib/stencil/Grid_lebesgue.h create mode 100644 tests/Grid_cshift_red_black.cc diff --git a/benchmarks/Grid_comms.cc b/benchmarks/Grid_comms.cc index c218070b..a24e7ed8 100644 --- a/benchmarks/Grid_comms.cc +++ b/benchmarks/Grid_comms.cc @@ -82,7 +82,7 @@ int main (int argc, char ** argv) double rbytes = xbytes; double bidibytes = xbytes+rbytes; - double time = stop-start; + double time = stop-start; // microseconds std::cout << lat<<"\t\t"< + +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 + }; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + int threads = GridThread::GetThreads(); + std::cout << "Grid is setup to use "< latt4 = GridDefaultLatt(); + std::vector simd4 = GridDefaultSimd(Nd,vComplexF::Nsimd()); + std::vector mpi4 = GridDefaultMpi(); + + assert(latt4.size()==4 ); + assert(simd4.size()==4 ); + assert(mpi4.size() ==4 ); + + const int Ls=1; + std::vector latt5({Ls,latt4[0],latt4[1],latt4[2],latt4[3]}); + std::vector simd5({1 ,simd4[0],simd4[1],simd4[2],simd4[3]}); + std::vector mpi5({1 , mpi4[0], mpi4[1], mpi4[2], mpi4[3]}); + std::vector cb5({0,1,1,1,1}); // Checkerboard 4d only + int cbd=1; // use dim-1 to reduce + + // Four dim grid for gauge field U + GridCartesian UGrid(latt4,simd4,mpi4); + GridRedBlackCartesian UrbGrid(&UGrid); + + // Five dim grid for fermions F + GridCartesian FGrid(latt5,simd5,mpi5); + GridRedBlackCartesian FrbGrid(latt5,simd5,mpi5,cb5,cbd); + + std::vector seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + + 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); + + GridParallelRNG RNG4(&UGrid); RNG4.SeedFixedIntegers(seeds4); + 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(Umu,U[nn],nn); } +#endif for(int mu=0;mu(Umu,mu); @@ -80,9 +82,9 @@ int main (int argc, char ** argv) } } } - + ref = -0.5*ref; RealD mass=0.1; - WilsonMatrix Dw(Umu,Grid,RBGrid,mass); + WilsonFermion Dw(Umu,Grid,RBGrid,mass); std::cout << "Calling Dw"< HermOp(Dw); + // HermitianOperator HermOp(Dw); // ConjugateGradient CG(1.0e-8,10000); // CG(HermOp,src,result); @@ -53,7 +53,7 @@ int main (int argc, char ** argv) pickCheckerboard(Odd,src_o,src); result_o=zero; - HermitianCheckerBoardedOperator HermOpEO(Dw); + HermitianCheckerBoardedOperator HermOpEO(Dw); ConjugateGradient CG(1.0e-8,10000); CG(HermOpEO,src_o,result_o); diff --git a/benchmarks/Grid_wilson_cg_schur.cc b/benchmarks/Grid_wilson_cg_schur.cc index af630ae1..28db1d4b 100644 --- a/benchmarks/Grid_wilson_cg_schur.cc +++ b/benchmarks/Grid_wilson_cg_schur.cc @@ -37,7 +37,7 @@ int main (int argc, char ** argv) LatticeFermion resid(&Grid); RealD mass=0.5; - WilsonMatrix Dw(Umu,Grid,RBGrid,mass); + WilsonFermion Dw(Umu,Grid,RBGrid,mass); ConjugateGradient CG(1.0e-8,10000); SchurRedBlackSolve SchurSolver(CG); diff --git a/benchmarks/Grid_wilson_cg_unprec.cc b/benchmarks/Grid_wilson_cg_unprec.cc index 15302aab..905dfde5 100644 --- a/benchmarks/Grid_wilson_cg_unprec.cc +++ b/benchmarks/Grid_wilson_cg_unprec.cc @@ -47,9 +47,9 @@ int main (int argc, char ** argv) } RealD mass=0.5; - WilsonMatrix Dw(Umu,Grid,RBGrid,mass); + WilsonFermion Dw(Umu,Grid,RBGrid,mass); - HermitianOperator HermOp(Dw); + HermitianOperator HermOp(Dw); ConjugateGradient CG(1.0e-8,10000); CG(HermOp,src,result); diff --git a/benchmarks/Grid_wilson_evenodd.cc b/benchmarks/Grid_wilson_evenodd.cc index 4bc8c357..a073139d 100644 --- a/benchmarks/Grid_wilson_evenodd.cc +++ b/benchmarks/Grid_wilson_evenodd.cc @@ -60,7 +60,7 @@ int main (int argc, char ** argv) RealD mass=0.1; - WilsonMatrix Dw(Umu,Grid,RBGrid,mass); + WilsonFermion Dw(Umu,Grid,RBGrid,mass); LatticeFermion src_e (&RBGrid); LatticeFermion src_o (&RBGrid); diff --git a/benchmarks/Makefile.am b/benchmarks/Makefile.am index 95ae5eca..b2649669 100644 --- a/benchmarks/Makefile.am +++ b/benchmarks/Makefile.am @@ -5,11 +5,14 @@ AM_LDFLAGS = -L$(top_builddir)/lib # # Test code # -bin_PROGRAMS = Grid_wilson Grid_comms Grid_memory_bandwidth Grid_su3 Grid_wilson_cg_unprec Grid_wilson_evenodd Grid_wilson_cg_prec Grid_wilson_cg_schur +bin_PROGRAMS = Grid_wilson Grid_comms Grid_memory_bandwidth Grid_su3 Grid_wilson_cg_unprec Grid_wilson_evenodd Grid_wilson_cg_prec Grid_wilson_cg_schur Grid_dwf Grid_wilson_SOURCES = Grid_wilson.cc Grid_wilson_LDADD = -lGrid +Grid_dwf_SOURCES = Grid_dwf.cc +Grid_dwf_LDADD = -lGrid + Grid_wilson_evenodd_SOURCES = Grid_wilson_evenodd.cc Grid_wilson_evenodd_LDADD = -lGrid diff --git a/lib/Grid.h b/lib/Grid.h index 1e513a46..7fa56892 100644 --- a/lib/Grid.h +++ b/lib/Grid.h @@ -60,7 +60,7 @@ #include // subdir aggregate -#include +#include #include namespace Grid { diff --git a/lib/Grid_init.cc b/lib/Grid_init.cc index f72393cb..580e602a 100644 --- a/lib/Grid_init.cc +++ b/lib/Grid_init.cc @@ -142,7 +142,11 @@ void Grid_init(int *argc,char ***argv) Grid_quiesce_nodes(); } if( GridCmdOptionExists(*argv,*argv+*argc,"--dslash-opt") ){ - WilsonMatrix::HandOptDslash=1; + WilsonFermion::HandOptDslash=1; + FiveDimWilsonFermion::HandOptDslash=1; + } + if( GridCmdOptionExists(*argv,*argv+*argc,"--lebesgue") ){ + LebesgueOrder::UseLebesgueOrder=1; } GridParseLayout(*argv,*argc, Grid_default_latt, diff --git a/lib/Grid_stencil.h b/lib/Grid_stencil.h index fa22361b..50d22453 100644 --- a/lib/Grid_stencil.h +++ b/lib/Grid_stencil.h @@ -1,6 +1,8 @@ #ifndef GRID_STENCIL_H #define GRID_STENCIL_H +#include // subdir aggregate + ////////////////////////////////////////////////////////////////////////////////////////// // Must not lose sight that goal is to be able to construct really efficient // gather to a point stencil code. CSHIFT is not the best way, so need @@ -38,29 +40,12 @@ ////////////////////////////////////////////////////////////////////////////////////////// namespace Grid { - - + class CartesianStencil { // Stencil runs along coordinate axes only; NO diagonal fill in. public: typedef uint32_t StencilInteger; - - - - StencilInteger alignup(StencilInteger n){ - n--; // 1000 0011 --> 1000 0010 - n |= n >> 1; // 1000 0010 | 0100 0001 = 1100 0011 - n |= n >> 2; // 1100 0011 | 0011 0000 = 1111 0011 - n |= n >> 4; // 1111 0011 | 0000 1111 = 1111 1111 - n |= n >> 8; // ... (At this point all bits are 1, so further bitwise-or - n |= n >> 16; // operations produce no effect.) - n++; // 1111 1111 --> 1 0000 0000 - return n; - }; - void LebesgueOrder(void); - - std::vector _LebesgueReorder; int _checkerboard; int _npoints; // Move to template param? @@ -131,8 +116,8 @@ namespace Grid { // Gather phase int sshift [2]; if ( comm_dim ) { - sshift[0] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,0); - sshift[1] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,1); + sshift[0] = _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,Even); + sshift[1] = _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,Odd); if ( sshift[0] == sshift[1] ) { if (splice_dim) { GatherStartCommsSimd(source,dimension,shift,0x3,u_comm_buf,u_comm_offset,compress); @@ -179,8 +164,8 @@ namespace Grid { std::vector > send_buf(buffer_size); // hmm... std::vector > recv_buf(buffer_size); - int cb= (cbmask==0x2)? 1 : 0; - int sshift= _grid->CheckerBoardShift(rhs.checkerboard,dimension,shift,cb); + int cb= (cbmask==0x2)? Odd : Even; + int sshift= _grid->CheckerBoardShiftForCB(rhs.checkerboard,dimension,shift,cb); for(int x=0;xCheckerBoardShift(rhs.checkerboard,dimension,shift,cb); + int cb = (cbmask==0x2)? Odd : Even; + int sshift= _grid->CheckerBoardShiftForCB(rhs.checkerboard,dimension,shift,cb); // loop over outer coord planes orthog to dim for(int x=0;x -#define PARALLEL_FOR_LOOP _Pragma("omp parallel for") +#define PARALLEL_FOR_LOOP _Pragma("omp parallel for ") #define PARALLEL_NESTED_LOOP2 _Pragma("omp parallel for collapse(2)") #else #define PARALLEL_FOR_LOOP diff --git a/lib/Makefile.am b/lib/Makefile.am index a5b89c0a..aa531df7 100644 --- a/lib/Makefile.am +++ b/lib/Makefile.am @@ -16,87 +16,103 @@ endif lib_LIBRARIES = libGrid.a libGrid_a_SOURCES = \ Grid_init.cc \ + stencil/Grid_lebesgue.cc \ stencil/Grid_stencil_common.cc \ - qcd/Grid_qcd_dirac.cc \ - qcd/Grid_qcd_dhop.cc \ - qcd/Grid_qcd_dhop_hand.cc \ - qcd/Grid_qcd_wilson_dop.cc \ algorithms/approx/Zolotarev.cc \ algorithms/approx/Remez.cc \ + qcd/action/fermion/FiveDimWilsonFermion.cc\ + qcd/action/fermion/WilsonFermion.cc\ + qcd/action/fermion/WilsonKernels.cc\ + qcd/action/fermion/WilsonKernelsHand.cc\ + qcd/Dirac.cc\ $(extra_sources) # # Include files # -nobase_include_HEADERS = algorithms/approx/bigfloat.h \ - algorithms/approx/Chebyshev.h \ - algorithms/approx/Remez.h \ - algorithms/approx/Zolotarev.h \ - algorithms/iterative/ConjugateGradient.h \ - algorithms/iterative/NormalEquations.h \ - algorithms/iterative/SchurRedBlack.h \ - algorithms/LinearOperator.h \ - algorithms/SparseMatrix.h \ - cartesian/Grid_cartesian_base.h \ - cartesian/Grid_cartesian_full.h \ - cartesian/Grid_cartesian_red_black.h \ - communicator/Grid_communicator_base.h \ - cshift/Grid_cshift_common.h \ - cshift/Grid_cshift_mpi.h \ - cshift/Grid_cshift_none.h \ - Grid.h \ - Grid_algorithms.h \ - Grid_aligned_allocator.h \ - Grid_cartesian.h \ - Grid_communicator.h \ - Grid_comparison.h \ - Grid_cshift.h \ - Grid_extract.h \ - Grid_lattice.h \ - Grid_math.h \ - Grid_simd.h \ - Grid_stencil.h \ - Grid_threads.h \ - lattice/Grid_lattice_arith.h \ - lattice/Grid_lattice_base.h \ - lattice/Grid_lattice_comparison.h \ - lattice/Grid_lattice_conformable.h \ - lattice/Grid_lattice_coordinate.h \ - lattice/Grid_lattice_ET.h \ - lattice/Grid_lattice_local.h \ - lattice/Grid_lattice_overload.h \ - lattice/Grid_lattice_peekpoke.h \ - lattice/Grid_lattice_reality.h \ - lattice/Grid_lattice_reduction.h \ - lattice/Grid_lattice_rng.h \ - lattice/Grid_lattice_trace.h \ - lattice/Grid_lattice_transfer.h \ - lattice/Grid_lattice_transpose.h \ - lattice/Grid_lattice_where.h \ - math/Grid_math_arith.h \ - math/Grid_math_arith_add.h \ - math/Grid_math_arith_mac.h \ - math/Grid_math_arith_mul.h \ - math/Grid_math_arith_scalar.h \ - math/Grid_math_arith_sub.h \ - math/Grid_math_inner.h \ - math/Grid_math_outer.h \ - math/Grid_math_peek.h \ - math/Grid_math_poke.h \ - math/Grid_math_reality.h \ - math/Grid_math_tensors.h \ - math/Grid_math_trace.h \ - math/Grid_math_traits.h \ - math/Grid_math_transpose.h \ - parallelIO/GridNerscIO.h \ - qcd/Grid_qcd.h \ - qcd/Grid_qcd_2spinor.h \ - qcd/Grid_qcd_dirac.h \ - qcd/Grid_qcd_wilson_dop.h \ - simd/Grid_vector_types.h \ - simd/Grid_sse4.h \ - simd/Grid_avx.h \ - simd/Grid_avx512.h - - +nobase_include_HEADERS=\ + ./algorithms/approx/bigfloat.h\ + ./algorithms/approx/bigfloat_double.h\ + ./algorithms/approx/Chebyshev.h\ + ./algorithms/approx/Remez.h\ + ./algorithms/approx/Zolotarev.h\ + ./algorithms/iterative/ConjugateGradient.h\ + ./algorithms/iterative/NormalEquations.h\ + ./algorithms/iterative/SchurRedBlack.h\ + ./algorithms/LinearOperator.h\ + ./algorithms/SparseMatrix.h\ + ./cartesian/Grid_cartesian_base.h\ + ./cartesian/Grid_cartesian_full.h\ + ./cartesian/Grid_cartesian_red_black.h\ + ./communicator/Grid_communicator_base.h\ + ./cshift/Grid_cshift_common.h\ + ./cshift/Grid_cshift_mpi.h\ + ./cshift/Grid_cshift_none.h\ + ./Grid.h\ + ./Grid_algorithms.h\ + ./Grid_aligned_allocator.h\ + ./Grid_cartesian.h\ + ./Grid_communicator.h\ + ./Grid_comparison.h\ + ./Grid_config.h\ + ./Grid_cshift.h\ + ./Grid_extract.h\ + ./Grid_lattice.h\ + ./Grid_math.h\ + ./Grid_simd.h\ + ./Grid_stencil.h\ + ./Grid_threads.h\ + ./lattice/Grid_lattice_arith.h\ + ./lattice/Grid_lattice_base.h\ + ./lattice/Grid_lattice_comparison.h\ + ./lattice/Grid_lattice_conformable.h\ + ./lattice/Grid_lattice_coordinate.h\ + ./lattice/Grid_lattice_ET.h\ + ./lattice/Grid_lattice_local.h\ + ./lattice/Grid_lattice_overload.h\ + ./lattice/Grid_lattice_peekpoke.h\ + ./lattice/Grid_lattice_reality.h\ + ./lattice/Grid_lattice_reduction.h\ + ./lattice/Grid_lattice_rng.h\ + ./lattice/Grid_lattice_trace.h\ + ./lattice/Grid_lattice_transfer.h\ + ./lattice/Grid_lattice_transpose.h\ + ./lattice/Grid_lattice_where.h\ + ./math/Grid_math_arith.h\ + ./math/Grid_math_arith_add.h\ + ./math/Grid_math_arith_mac.h\ + ./math/Grid_math_arith_mul.h\ + ./math/Grid_math_arith_scalar.h\ + ./math/Grid_math_arith_sub.h\ + ./math/Grid_math_inner.h\ + ./math/Grid_math_outer.h\ + ./math/Grid_math_peek.h\ + ./math/Grid_math_poke.h\ + ./math/Grid_math_reality.h\ + ./math/Grid_math_tensors.h\ + ./math/Grid_math_trace.h\ + ./math/Grid_math_traits.h\ + ./math/Grid_math_transpose.h\ + ./parallelIO/GridNerscIO.h\ + ./qcd/action/Actions.h\ + ./qcd/action/fermion/FermionAction.h\ + ./qcd/action/fermion/FiveDimWilsonFermion.h\ + ./qcd/action/fermion/WilsonCompressor.h\ + ./qcd/action/fermion/WilsonFermion.h\ + ./qcd/action/fermion/WilsonKernels.h\ + ./qcd/Dirac.h\ + ./qcd/QCD.h\ + ./qcd/TwoSpinor.h\ + ./qcd/FermionAction.h\ + ./simd/Grid_avx.h\ + ./simd/Grid_avx512.h\ + ./simd/Grid_qpx.h\ + ./simd/Grid_sse4.h\ + ./simd/Grid_vector_types.h\ + ./simd/Old/Grid_vComplexD.h\ + ./simd/Old/Grid_vComplexF.h\ + ./simd/Old/Grid_vInteger.h\ + ./simd/Old/Grid_vRealD.h\ + ./simd/Old/Grid_vRealF.h\ + ./stencil/Grid_lebesgue.h diff --git a/lib/algorithms/SparseMatrix.h b/lib/algorithms/SparseMatrix.h index 7bfe959b..9c955e9a 100644 --- a/lib/algorithms/SparseMatrix.h +++ b/lib/algorithms/SparseMatrix.h @@ -10,7 +10,7 @@ namespace Grid { ///////////////////////////////////////////////////////////////////////////////////////////// template class SparseMatrixBase { public: - GridBase *_grid; + virtual GridBase *Grid(void) =0; // Full checkerboar operations virtual RealD M (const Field &in, Field &out)=0; virtual RealD Mdag (const Field &in, Field &out)=0; @@ -19,7 +19,6 @@ namespace Grid { ni=M(in,tmp); no=Mdag(tmp,out); } - SparseMatrixBase(GridBase *grid) : _grid(grid) {}; }; ///////////////////////////////////////////////////////////////////////////////////////////// @@ -27,7 +26,7 @@ namespace Grid { ///////////////////////////////////////////////////////////////////////////////////////////// template class CheckerBoardedSparseMatrixBase : public SparseMatrixBase { public: - GridBase *_cbgrid; + virtual GridBase *RedBlackGrid(void)=0; // half checkerboard operaions virtual void Meooe (const Field &in, Field &out)=0; virtual void Mooee (const Field &in, Field &out)=0; @@ -62,9 +61,7 @@ namespace Grid { Field tmp(in._grid); ni=Mpc(in,tmp); no=MpcDag(tmp,out); - // std::cout<<"MpcDagMpc "<(grid), _cbgrid(cbgrid) {}; }; } diff --git a/lib/algorithms/iterative/SchurRedBlack.h b/lib/algorithms/iterative/SchurRedBlack.h index d3a76d7f..109f554b 100644 --- a/lib/algorithms/iterative/SchurRedBlack.h +++ b/lib/algorithms/iterative/SchurRedBlack.h @@ -60,8 +60,8 @@ namespace Grid { // FIXME CGdiagonalMee not implemented virtual function // FIXME use CBfactorise to control schur decomp - GridBase *grid = _Matrix._cbgrid; - GridBase *fgrid= _Matrix._grid; + GridBase *grid = _Matrix.RedBlackGrid(); + GridBase *fgrid= _Matrix.Grid(); Field src_e(grid); Field src_o(grid); diff --git a/lib/cartesian/Grid_cartesian_base.h b/lib/cartesian/Grid_cartesian_base.h index a74773f8..e93125c1 100644 --- a/lib/cartesian/Grid_cartesian_base.h +++ b/lib/cartesian/Grid_cartesian_base.h @@ -52,16 +52,13 @@ public: //////////////////////////////////////////////////////////////// virtual int CheckerBoarded(int dim)=0; virtual int CheckerBoard(std::vector site)=0; - virtual int CheckerBoardDestination(int source_cb,int shift)=0; + virtual int CheckerBoardDestination(int source_cb,int shift,int dim)=0; virtual int CheckerBoardShift(int source_cb,int dim,int shift,int osite)=0; - inline int CheckerBoardFromOindex (int Oindex){ + virtual int CheckerBoardShiftForCB(int source_cb,int dim,int shift,int cb)=0; + int CheckerBoardFromOindex (int Oindex){ std::vector ocoor; oCoorFromOindex(ocoor,Oindex); - int ss=0; - for(int d=0;d<_ndimension;d++){ - ss=ss+ocoor[d]; - } - return ss&0x1; + return CheckerBoard(ocoor); } ////////////////////////////////////////////////////////////////////////////////////////////// diff --git a/lib/cartesian/Grid_cartesian_full.h b/lib/cartesian/Grid_cartesian_full.h index 73bd08b3..330bbfaf 100644 --- a/lib/cartesian/Grid_cartesian_full.h +++ b/lib/cartesian/Grid_cartesian_full.h @@ -18,11 +18,14 @@ public: virtual int CheckerBoard(std::vector site){ return 0; } - virtual int CheckerBoardDestination(int cb,int shift){ + virtual int CheckerBoardDestination(int cb,int shift,int dim){ return 0; } + virtual int CheckerBoardShiftForCB(int source_cb,int dim,int shift, int ocb){ + return shift; + } virtual int CheckerBoardShift(int source_cb,int dim,int shift, int osite){ - return shift; + return shift; } GridCartesian(std::vector &dimensions, std::vector &simd_layout, diff --git a/lib/cartesian/Grid_cartesian_red_black.h b/lib/cartesian/Grid_cartesian_red_black.h index ff2d3ba8..ace36edb 100644 --- a/lib/cartesian/Grid_cartesian_red_black.h +++ b/lib/cartesian/Grid_cartesian_red_black.h @@ -8,6 +8,10 @@ namespace Grid { static const int CbBlack=1; static const int Even =CbRed; static const int Odd =CbBlack; + + // Perhaps these are misplaced and + // should be in sparse matrix. + // Also should make these a named enum type static const int DaggerNo=0; static const int DaggerYes=1; @@ -15,116 +19,174 @@ namespace Grid { class GridRedBlackCartesian : public GridBase { public: + std::vector _checker_dim_mask; + int _checker_dim; + virtual int CheckerBoarded(int dim){ - if( dim==0) return 1; + if( dim==_checker_dim) return 1; else return 0; } virtual int CheckerBoard(std::vector site){ - return (site[0]+site[1]+site[2]+site[3])&0x1; + int linear=0; + assert(site.size()==_ndimension); + for(int d=0;d<_ndimension;d++){ + if(_checker_dim_mask[d]) + linear=linear+site[d]; + } + return (linear&0x1); } + // Depending on the cb of site, we toggle source cb. // for block #b, element #e = (b, e) // we need - virtual int CheckerBoardShift(int source_cb,int dim,int shift,int osite){ + virtual int CheckerBoardShiftForCB(int source_cb,int dim,int shift,int ocb){ + if(dim != _checker_dim) return shift; - if(dim != 0) return shift; - - int fulldim =_fdimensions[0]; + int fulldim =_fdimensions[dim]; shift = (shift+fulldim)%fulldim; // Probably faster with table lookup; // or by looping over x,y,z and multiply rather than computing checkerboard. - int ocb=CheckerBoardFromOindex(osite); if ( (source_cb+ocb)&1 ) { + return (shift)/2; } else { return (shift+1)/2; } } + virtual int CheckerBoardShift(int source_cb,int dim,int shift,int osite){ - virtual int CheckerBoardDestination(int source_cb,int shift){ - if ((shift+_fdimensions[0])&0x1) { + if(dim != _checker_dim) return shift; + + int ocb=CheckerBoardFromOindex(osite); + + return CheckerBoardShiftForCB(source_cb,dim,shift,ocb); + } + + virtual int CheckerBoardDestination(int source_cb,int shift,int dim){ + if ( _checker_dim_mask[dim] ) { + // If _fdimensions[checker_dim] is odd, then shifting by 1 in other dims + // does NOT cause a parity hop. + int add=(dim==_checker_dim) ? 0 : _fdimensions[_checker_dim]; + if ( (shift+add) &0x1) { return 1-source_cb; } else { return source_cb; } + } else { + return source_cb; + + } }; GridRedBlackCartesian(GridBase *base) : GridRedBlackCartesian(base->_fdimensions,base->_simd_layout,base->_processors) {}; GridRedBlackCartesian(std::vector &dimensions, std::vector &simd_layout, - std::vector &processor_grid ) : GridBase(processor_grid) + std::vector &processor_grid, + std::vector &checker_dim_mask, + int checker_dim + ) : GridBase(processor_grid) + { + Init(dimensions,simd_layout,processor_grid,checker_dim_mask,checker_dim); + } + GridRedBlackCartesian(std::vector &dimensions, + std::vector &simd_layout, + std::vector &processor_grid) : GridBase(processor_grid) + { + std::vector checker_dim_mask(dimensions.size(),1); + Init(dimensions,simd_layout,processor_grid,checker_dim_mask,0); + } + void Init(std::vector &dimensions, + std::vector &simd_layout, + std::vector &processor_grid, + std::vector &checker_dim_mask, + int checker_dim) { /////////////////////// // Grid information /////////////////////// - _ndimension = dimensions.size(); - - _fdimensions.resize(_ndimension); - _gdimensions.resize(_ndimension); - _ldimensions.resize(_ndimension); - _rdimensions.resize(_ndimension); - _simd_layout.resize(_ndimension); - - _ostride.resize(_ndimension); - _istride.resize(_ndimension); - - _fsites = _gsites = _osites = _isites = 1; + _checker_dim = checker_dim; + assert(checker_dim_mask[checker_dim]==1); + _ndimension = dimensions.size(); + assert(checker_dim_mask.size()==_ndimension); + assert(processor_grid.size()==_ndimension); + assert(simd_layout.size()==_ndimension); + + _fdimensions.resize(_ndimension); + _gdimensions.resize(_ndimension); + _ldimensions.resize(_ndimension); + _rdimensions.resize(_ndimension); + _simd_layout.resize(_ndimension); + + _ostride.resize(_ndimension); + _istride.resize(_ndimension); + + _fsites = _gsites = _osites = _isites = 1; + + _checker_dim_mask=checker_dim_mask; - for(int d=0;d<_ndimension;d++){ - _fdimensions[d] = dimensions[d]; - _gdimensions[d] = _fdimensions[d]; - _fsites = _fsites * _fdimensions[d]; - _gsites = _gsites * _gdimensions[d]; - - if (d==0) _gdimensions[0] = _gdimensions[0]/2; // Remove a checkerboard - _ldimensions[d] = _gdimensions[d]/_processors[d]; - - // Use a reduced simd grid - _simd_layout[d] = simd_layout[d]; - _rdimensions[d]= _ldimensions[d]/_simd_layout[d]; - - _osites *= _rdimensions[d]; - _isites *= _simd_layout[d]; - - // Addressing support - if ( d==0 ) { - _ostride[d] = 1; - _istride[d] = 1; - } else { - _ostride[d] = _ostride[d-1]*_rdimensions[d-1]; - _istride[d] = _istride[d-1]*_simd_layout[d-1]; - } - } - - //////////////////////////////////////////////////////////////////////////////////////////// - // subplane information - //////////////////////////////////////////////////////////////////////////////////////////// - _slice_block.resize(_ndimension); - _slice_stride.resize(_ndimension); - _slice_nblock.resize(_ndimension); + for(int d=0;d<_ndimension;d++){ + _fdimensions[d] = dimensions[d]; + _gdimensions[d] = _fdimensions[d]; + _fsites = _fsites * _fdimensions[d]; + _gsites = _gsites * _gdimensions[d]; - int block =1; - int nblock=1; - for(int d=0;d<_ndimension;d++) nblock*=_rdimensions[d]; - - for(int d=0;d<_ndimension;d++){ - nblock/=_rdimensions[d]; - _slice_block[d] =block; - _slice_stride[d]=_ostride[d]*_rdimensions[d]; - _slice_nblock[d]=nblock; - block = block*_rdimensions[d]; - } + if (d==_checker_dim) { + _gdimensions[d] = _gdimensions[d]/2; // Remove a checkerboard + } + _ldimensions[d] = _gdimensions[d]/_processors[d]; + + // Use a reduced simd grid + _simd_layout[d] = simd_layout[d]; + _rdimensions[d]= _ldimensions[d]/_simd_layout[d]; + + _osites *= _rdimensions[d]; + _isites *= _simd_layout[d]; + + // Addressing support + if ( d==0 ) { + _ostride[d] = 1; + _istride[d] = 1; + } else { + _ostride[d] = _ostride[d-1]*_rdimensions[d-1]; + _istride[d] = _istride[d-1]*_simd_layout[d-1]; + } + } + //////////////////////////////////////////////////////////////////////////////////////////// + // subplane information + //////////////////////////////////////////////////////////////////////////////////////////// + _slice_block.resize(_ndimension); + _slice_stride.resize(_ndimension); + _slice_nblock.resize(_ndimension); + + int block =1; + int nblock=1; + for(int d=0;d<_ndimension;d++) nblock*=_rdimensions[d]; + + for(int d=0;d<_ndimension;d++){ + nblock/=_rdimensions[d]; + _slice_block[d] =block; + _slice_stride[d]=_ostride[d]*_rdimensions[d]; + _slice_nblock[d]=nblock; + block = block*_rdimensions[d]; + } + }; protected: virtual int oIndex(std::vector &coor) { - int idx=_ostride[0]*((coor[0]/2)%_rdimensions[0]); - for(int d=1;d<_ndimension;d++) idx+=_ostride[d]*(coor[d]%_rdimensions[d]); + int idx=0; + for(int d=0;d<_ndimension;d++) { + if( d==_checker_dim ) { + idx+=_ostride[d]*((coor[d]/2)%_rdimensions[d]); + } else { + idx+=_ostride[d]*(coor[d]%_rdimensions[d]); + } + } return idx; }; diff --git a/lib/cshift/Grid_cshift_common.h b/lib/cshift/Grid_cshift_common.h index c369fe1c..06e812d9 100644 --- a/lib/cshift/Grid_cshift_common.h +++ b/lib/cshift/Grid_cshift_common.h @@ -175,7 +175,6 @@ PARALLEL_NESTED_LOOP2 if ( ocb&cbmask ) { //lhs._odata[lo+o]=rhs._odata[ro+o]; vstream(lhs._odata[lo+o],rhs._odata[ro+o]); - } } @@ -217,8 +216,8 @@ template void Cshift_local(Lattice& ret,Lattice &rhs,int { int sshift[2]; - sshift[0] = rhs._grid->CheckerBoardShift(rhs.checkerboard,dimension,shift,0); - sshift[1] = rhs._grid->CheckerBoardShift(rhs.checkerboard,dimension,shift,1); + sshift[0] = rhs._grid->CheckerBoardShiftForCB(rhs.checkerboard,dimension,shift,Even); + sshift[1] = rhs._grid->CheckerBoardShiftForCB(rhs.checkerboard,dimension,shift,Odd); if ( sshift[0] == sshift[1] ) { Cshift_local(ret,rhs,dimension,shift,0x3); @@ -239,8 +238,7 @@ template Lattice Cshift_local(Lattice &ret,Lattice // Map to always positive shift modulo global full dimension. shift = (shift+fd)%fd; - ret.checkerboard = grid->CheckerBoardDestination(rhs.checkerboard,shift); - + ret.checkerboard = grid->CheckerBoardDestination(rhs.checkerboard,shift,dimension); // the permute type int permute_dim =grid->PermuteDim(dimension); int permute_type=grid->PermuteType(dimension); @@ -250,11 +248,11 @@ template Lattice Cshift_local(Lattice &ret,Lattice int o = 0; int bo = x * grid->_ostride[dimension]; - int cb= (cbmask==0x2)? 1 : 0; + int cb= (cbmask==0x2)? Odd : Even; - int sshift = grid->CheckerBoardShift(rhs.checkerboard,dimension,shift,cb); + int sshift = grid->CheckerBoardShiftForCB(rhs.checkerboard,dimension,shift,cb); int sx = (x+sshift)%rd; - + int permute_slice=0; if(permute_dim){ int wrap = sshift/rd; diff --git a/lib/cshift/Grid_cshift_mpi.h b/lib/cshift/Grid_cshift_mpi.h index 569017d8..8c0badcd 100644 --- a/lib/cshift/Grid_cshift_mpi.h +++ b/lib/cshift/Grid_cshift_mpi.h @@ -39,8 +39,8 @@ template void Cshift_comms(Lattice& ret,Lattice &rhs,int { int sshift[2]; - sshift[0] = rhs._grid->CheckerBoardShift(rhs.checkerboard,dimension,shift,0); - sshift[1] = rhs._grid->CheckerBoardShift(rhs.checkerboard,dimension,shift,1); + sshift[0] = rhs._grid->CheckerBoardShiftForCB(rhs.checkerboard,dimension,shift,Even); + sshift[1] = rhs._grid->CheckerBoardShiftForCB(rhs.checkerboard,dimension,shift,Odd); if ( sshift[0] == sshift[1] ) { Cshift_comms(ret,rhs,dimension,shift,0x3); @@ -54,8 +54,8 @@ template void Cshift_comms_simd(Lattice& ret,Lattice &rh { int sshift[2]; - sshift[0] = rhs._grid->CheckerBoardShift(rhs.checkerboard,dimension,shift,0); - sshift[1] = rhs._grid->CheckerBoardShift(rhs.checkerboard,dimension,shift,1); + sshift[0] = rhs._grid->CheckerBoardShiftForCB(rhs.checkerboard,dimension,shift,Even); + sshift[1] = rhs._grid->CheckerBoardShiftForCB(rhs.checkerboard,dimension,shift,Odd); if ( sshift[0] == sshift[1] ) { Cshift_comms_simd(ret,rhs,dimension,shift,0x3); @@ -87,8 +87,8 @@ template void Cshift_comms(Lattice &ret,Lattice &rhs,int std::vector > send_buf(buffer_size); std::vector > recv_buf(buffer_size); - int cb= (cbmask==0x2)? 1 : 0; - int sshift= rhs._grid->CheckerBoardShift(rhs.checkerboard,dimension,shift,cb); + int cb= (cbmask==0x2)? Odd : Even; + int sshift= rhs._grid->CheckerBoardShiftForCB(rhs.checkerboard,dimension,shift,cb); for(int x=0;x void Cshift_comms_simd(Lattice &ret,Lattice &r /////////////////////////////////////////// // Work out what to send where /////////////////////////////////////////// - int cb = (cbmask==0x2)? 1 : 0; - int sshift= grid->CheckerBoardShift(rhs.checkerboard,dimension,shift,cb); + int cb = (cbmask==0x2)? Odd : Even; + int sshift= grid->CheckerBoardShiftForCB(rhs.checkerboard,dimension,shift,cb); // loop over outer coord planes orthog to dim for(int x=0;x Lattice Cshift(Lattice &rhs,int dimension,int shift) { Lattice ret(rhs._grid); - ret.checkerboard = rhs._grid->CheckerBoardDestination(rhs.checkerboard,shift); + ret.checkerboard = rhs._grid->CheckerBoardDestination(rhs.checkerboard,shift,dimension); Cshift_local(ret,rhs,dimension,shift); return ret; } diff --git a/lib/lattice/Grid_lattice_peekpoke.h b/lib/lattice/Grid_lattice_peekpoke.h index 7bffdbb8..9416226a 100644 --- a/lib/lattice/Grid_lattice_peekpoke.h +++ b/lib/lattice/Grid_lattice_peekpoke.h @@ -116,7 +116,7 @@ PARALLEL_FOR_LOOP int Nsimd = grid->Nsimd(); - assert( l.checkerboard== l._grid->CheckerBoard(site)); + assert( l.checkerboard == l._grid->CheckerBoard(site)); assert( sizeof(sobj)*Nsimd == sizeof(vobj)); int rank,odx,idx; diff --git a/lib/lattice/Grid_lattice_transfer.h b/lib/lattice/Grid_lattice_transfer.h index 27ae27e1..f0046b0e 100644 --- a/lib/lattice/Grid_lattice_transfer.h +++ b/lib/lattice/Grid_lattice_transfer.h @@ -32,7 +32,6 @@ PARALLEL_FOR_LOOP cbos=half._grid->CheckerBoard(coor); if (cbos==cb) { - half._odata[ssh] = full._odata[ss]; ssh++; } @@ -45,7 +44,7 @@ PARALLEL_FOR_LOOP for(int ss=0;ssoSites();ss++){ std::vector coor; int cbos; - + full._grid->oCoorFromOindex(coor,ss); cbos=half._grid->CheckerBoard(coor); diff --git a/lib/qcd/Grid_qcd_dirac.cc b/lib/qcd/Dirac.cc similarity index 100% rename from lib/qcd/Grid_qcd_dirac.cc rename to lib/qcd/Dirac.cc diff --git a/lib/qcd/Grid_qcd_dirac.h b/lib/qcd/Dirac.h similarity index 100% rename from lib/qcd/Grid_qcd_dirac.h rename to lib/qcd/Dirac.h diff --git a/lib/qcd/Grid_qcd_wilson_dop.cc b/lib/qcd/Grid_qcd_wilson_dop.cc deleted file mode 100644 index 9a3f5f6a..00000000 --- a/lib/qcd/Grid_qcd_wilson_dop.cc +++ /dev/null @@ -1,217 +0,0 @@ -#include - -namespace Grid { -namespace QCD { - -const std::vector WilsonMatrix::directions ({0,1,2,3, 0, 1, 2, 3}); -const std::vector WilsonMatrix::displacements({1,1,1,1,-1,-1,-1,-1}); - - int WilsonMatrix::HandOptDslash; - - class WilsonCompressor { - public: - int mu; - int dag; - - WilsonCompressor(int _dag){ - mu=0; - dag=_dag; - assert((dag==0)||(dag==1)); - } - void Point(int p) { - mu=p; - }; - - vHalfSpinColourVector operator () (const vSpinColourVector &in) - { - vHalfSpinColourVector ret; - int mudag=mu; - if (dag) { - mudag=(mu+Nd)%(2*Nd); - } - switch(mudag) { - case Xp: - spProjXp(ret,in); - break; - case Yp: - spProjYp(ret,in); - break; - case Zp: - spProjZp(ret,in); - break; - case Tp: - spProjTp(ret,in); - break; - case Xm: - spProjXm(ret,in); - break; - case Ym: - spProjYm(ret,in); - break; - case Zm: - spProjZm(ret,in); - break; - case Tm: - spProjTm(ret,in); - break; - default: - assert(0); - break; - } - return ret; - } - }; - - WilsonMatrix::WilsonMatrix(LatticeGaugeField &_Umu,GridCartesian &Fgrid,GridRedBlackCartesian &Hgrid, double _mass) : - CheckerBoardedSparseMatrixBase(&Fgrid,&Hgrid), - Stencil ( _grid,npoint,Even,directions,displacements), - StencilEven(_cbgrid,npoint,Even,directions,displacements), // source is Even - StencilOdd (_cbgrid,npoint,Odd ,directions,displacements), // source is Odd - mass(_mass), - Umu(_grid), - UmuEven(_cbgrid), - UmuOdd (_cbgrid) - { - // Allocate the required comms buffer - comm_buf.resize(Stencil._unified_buffer_size); // this is always big enough to contain EO - - DoubleStore(Umu,_Umu); - pickCheckerboard(Even,UmuEven,Umu); - pickCheckerboard(Odd ,UmuOdd,Umu); - } - -void WilsonMatrix::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu) -{ - LatticeColourMatrix U(_grid); - - for(int mu=0;mu(Umu,mu); - pokeIndex(Uds,U,mu); - U = adj(Cshift(U,mu,-1)); - pokeIndex(Uds,U,mu+4); - } -} - -RealD WilsonMatrix::M(const LatticeFermion &in, LatticeFermion &out) -{ - out.checkerboard=in.checkerboard; - Dhop(in,out,DaggerNo); - out = (4+mass)*in - 0.5*out ; // FIXME : axpby_norm! fusion fun - return norm2(out); -} -RealD WilsonMatrix::Mdag(const LatticeFermion &in, LatticeFermion &out) -{ - out.checkerboard=in.checkerboard; - Dhop(in,out,DaggerYes); - out = (4+mass)*in - 0.5*out ; // FIXME : axpby_norm! fusion fun - return norm2(out); -} - -void WilsonMatrix::Meooe(const LatticeFermion &in, LatticeFermion &out) -{ - if ( in.checkerboard == Odd ) { - DhopEO(in,out,DaggerNo); - } else { - DhopOE(in,out,DaggerNo); - } - out = (-0.5)*out; // FIXME : scale factor in Dhop -} -void WilsonMatrix::MeooeDag(const LatticeFermion &in, LatticeFermion &out) -{ - if ( in.checkerboard == Odd ) { - DhopEO(in,out,DaggerYes); - } else { - DhopOE(in,out,DaggerYes); - } - out = (-0.5)*out; // FIXME : scale factor in Dhop -} -void WilsonMatrix::Mooee(const LatticeFermion &in, LatticeFermion &out) -{ - out.checkerboard = in.checkerboard; - out = (4.0+mass)*in; - return ; -} -void WilsonMatrix::MooeeDag(const LatticeFermion &in, LatticeFermion &out) -{ - out.checkerboard = in.checkerboard; - Mooee(in,out); -} -void WilsonMatrix::MooeeInv(const LatticeFermion &in, LatticeFermion &out) -{ - out.checkerboard = in.checkerboard; - out = (1.0/(4.0+mass))*in; - return ; -} -void WilsonMatrix::MooeeInvDag(const LatticeFermion &in, LatticeFermion &out) -{ - out.checkerboard = in.checkerboard; - MooeeInv(in,out); -} - -void WilsonMatrix::DhopInternal(CartesianStencil & st,LatticeDoubledGaugeField & U, - const LatticeFermion &in, LatticeFermion &out,int dag) -{ - assert((dag==DaggerNo) ||(dag==DaggerYes)); - WilsonCompressor compressor(dag); - st.HaloExchange(in,comm_buf,compressor); - - if ( dag == DaggerYes ) { - if( HandOptDslash ) { -PARALLEL_FOR_LOOP - for(int sss=0;sssoSites();sss++){ - DiracOptHand::DhopSiteDag(st,U,comm_buf,sss,in,out); - } - } else { -PARALLEL_FOR_LOOP - for(int sss=0;sssoSites();sss++){ - DiracOpt::DhopSiteDag(st,U,comm_buf,sss,in,out); - } - } - } else { - if( HandOptDslash ) { -PARALLEL_FOR_LOOP - for(int sss=0;sssoSites();sss++){ - DiracOptHand::DhopSite(st,U,comm_buf,sss,in,out); - } - } else { -PARALLEL_FOR_LOOP - for(int sss=0;sssoSites();sss++){ - DiracOpt::DhopSite(st,U,comm_buf,sss,in,out); - } - } - } -} -void WilsonMatrix::DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag) -{ - conformable(in._grid,_cbgrid); // verifies half grid - conformable(in._grid,out._grid); // drops the cb check - - assert(in.checkerboard==Even); - out.checkerboard = Odd; - - DhopInternal(StencilEven,UmuOdd,in,out,dag); -} -void WilsonMatrix::DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag) -{ - conformable(in._grid,_cbgrid); // verifies half grid - conformable(in._grid,out._grid); // drops the cb check - - assert(in.checkerboard==Odd); - out.checkerboard = Even; - - DhopInternal(StencilOdd,UmuEven,in,out,dag); -} -void WilsonMatrix::Dhop(const LatticeFermion &in, LatticeFermion &out,int dag) -{ - conformable(in._grid,_grid); // verifies full grid - conformable(in._grid,out._grid); - - out.checkerboard = in.checkerboard; - - DhopInternal(Stencil,Umu,in,out,dag); -} - -}} - - - diff --git a/lib/qcd/Grid_qcd_wilson_dop.h b/lib/qcd/Grid_qcd_wilson_dop.h deleted file mode 100644 index 87418603..00000000 --- a/lib/qcd/Grid_qcd_wilson_dop.h +++ /dev/null @@ -1,105 +0,0 @@ -#ifndef GRID_QCD_WILSON_DOP_H -#define GRID_QCD_WILSON_DOP_H - - -namespace Grid { - - namespace QCD { - - // Should be in header? - const int Xp = 0; - const int Yp = 1; - const int Zp = 2; - const int Tp = 3; - const int Xm = 4; - const int Ym = 5; - const int Zm = 6; - const int Tm = 7; - - class WilsonMatrix : public CheckerBoardedSparseMatrixBase - { - //NB r=1; - public: - static int HandOptDslash; - - double mass; - // GridBase * grid; // Inherited - // GridBase * cbgrid; - - //Defines the stencils for even and odd - CartesianStencil Stencil; - CartesianStencil StencilEven; - CartesianStencil StencilOdd; - - // Copy of the gauge field , with even and odd subsets - LatticeDoubledGaugeField Umu; - LatticeDoubledGaugeField UmuEven; - LatticeDoubledGaugeField UmuOdd; - - static const int npoint=8; - static const std::vector directions ; - static const std::vector displacements; - static const int Xp,Xm,Yp,Ym,Zp,Zm,Tp,Tm; - - // Comms buffer - std::vector > comm_buf; - - // Constructor - WilsonMatrix(LatticeGaugeField &_Umu,GridCartesian &Fgrid,GridRedBlackCartesian &Hgrid,double _mass); - - // DoubleStore - void DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu); - - // override multiply - virtual RealD M (const LatticeFermion &in, LatticeFermion &out); - virtual RealD Mdag (const LatticeFermion &in, LatticeFermion &out); - - // half checkerboard operaions - virtual void Meooe (const LatticeFermion &in, LatticeFermion &out); - virtual void MeooeDag (const LatticeFermion &in, LatticeFermion &out); - virtual void Mooee (const LatticeFermion &in, LatticeFermion &out); - virtual void MooeeDag (const LatticeFermion &in, LatticeFermion &out); - virtual void MooeeInv (const LatticeFermion &in, LatticeFermion &out); - virtual void MooeeInvDag (const LatticeFermion &in, LatticeFermion &out); - - // non-hermitian hopping term; half cb or both - void Dhop (const LatticeFermion &in, LatticeFermion &out,int dag); - void DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag); - void DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag); - void DhopInternal(CartesianStencil & st,LatticeDoubledGaugeField &U, - const LatticeFermion &in, LatticeFermion &out,int dag); - - typedef iScalar > matrix; - - - }; - - - class DiracOpt { - public: - // These ones will need to be package intelligently. WilsonType base class - // for use by DWF etc.. - static void DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, - std::vector > &buf, - int ss,const LatticeFermion &in, LatticeFermion &out); - static void DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U, - std::vector > &buf, - int ss,const LatticeFermion &in, LatticeFermion &out); - - }; - class DiracOptHand { - public: - // These ones will need to be package intelligently. WilsonType base class - // for use by DWF etc.. - static void DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, - std::vector > &buf, - int ss,const LatticeFermion &in, LatticeFermion &out); - static void DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U, - std::vector > &buf, - int ss,const LatticeFermion &in, LatticeFermion &out); - - }; - - } -} -#endif diff --git a/lib/qcd/Grid_qcd.h b/lib/qcd/QCD.h similarity index 97% rename from lib/qcd/Grid_qcd.h rename to lib/qcd/QCD.h index 959f3529..7c45eb23 100644 --- a/lib/qcd/Grid_qcd.h +++ b/lib/qcd/QCD.h @@ -4,6 +4,16 @@ namespace Grid{ namespace QCD { + + static const int Xp = 0; + static const int Yp = 1; + static const int Zp = 2; + static const int Tp = 3; + static const int Xm = 4; + static const int Ym = 5; + static const int Zm = 6; + static const int Tm = 7; + static const int Nc=3; static const int Ns=4; static const int Nd=4; @@ -297,9 +307,8 @@ namespace QCD { } //namespace QCD } // Grid -#include -#include -//#include -#include +#include +#include +#include #endif diff --git a/lib/qcd/Grid_qcd_2spinor.h b/lib/qcd/TwoSpinor.h similarity index 100% rename from lib/qcd/Grid_qcd_2spinor.h rename to lib/qcd/TwoSpinor.h diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h new file mode 100644 index 00000000..c4e8a2f0 --- /dev/null +++ b/lib/qcd/action/Actions.h @@ -0,0 +1,10 @@ +#ifndef GRID_QCD_ACTIONS_H +#define GRID_QCD_ACTIONS_H + +#include +#include +#include +#include +#include + +#endif diff --git a/lib/qcd/action/fermion/FermionAction.h b/lib/qcd/action/fermion/FermionAction.h new file mode 100644 index 00000000..1b05174b --- /dev/null +++ b/lib/qcd/action/fermion/FermionAction.h @@ -0,0 +1,47 @@ +#ifndef GRID_QCD_WILSON_DOP_H +#define GRID_QCD_WILSON_DOP_H + +namespace Grid { + + namespace QCD { + + ////////////////////////////////////////////////////////////////////////////// + // Four component fermions + // Should type template the vector and gauge types + // Think about multiple representations + ////////////////////////////////////////////////////////////////////////////// + template + class FermionAction : public CheckerBoardedSparseMatrixBase + { + public: + + GridBase * Grid(void) { return FermionGrid(); }; // this is all the linalg routines need to know + GridBase * RedBlackGrid(void) { return FermionRedBlackGrid(); }; + + virtual GridBase *FermionGrid(void) =0; + virtual GridBase *FermionRedBlackGrid(void) =0; + virtual GridBase *GaugeGrid(void) =0; + virtual GridBase *GaugeRedBlackGrid(void) =0; + + // override multiply + virtual RealD M (const FermionField &in, FermionField &out)=0; + virtual RealD Mdag (const FermionField &in, FermionField &out)=0; + + // half checkerboard operaions + virtual void Meooe (const FermionField &in, FermionField &out)=0; + virtual void MeooeDag (const FermionField &in, FermionField &out)=0; + virtual void Mooee (const FermionField &in, FermionField &out)=0; + virtual void MooeeDag (const FermionField &in, FermionField &out)=0; + virtual void MooeeInv (const FermionField &in, FermionField &out)=0; + virtual void MooeeInvDag (const FermionField &in, FermionField &out)=0; + + // non-hermitian hopping term; half cb or both + virtual void Dhop (const FermionField &in, FermionField &out,int dag)=0; + virtual void DhopOE(const FermionField &in, FermionField &out,int dag)=0; + virtual void DhopEO(const FermionField &in, FermionField &out,int dag)=0; + + }; + + } +} +#endif diff --git a/lib/qcd/action/fermion/FiveDimWilsonFermion.cc b/lib/qcd/action/fermion/FiveDimWilsonFermion.cc new file mode 100644 index 00000000..43645899 --- /dev/null +++ b/lib/qcd/action/fermion/FiveDimWilsonFermion.cc @@ -0,0 +1,228 @@ +#include + +namespace Grid { +namespace QCD { + + // S-direction is INNERMOST and takes no part in the parity. + const std::vector FiveDimWilsonFermion::directions ({1,2,3,4, 1, 2, 3, 4}); + const std::vector FiveDimWilsonFermion::displacements({1,1,1,1,-1,-1,-1,-1}); + + int FiveDimWilsonFermion::HandOptDslash; + + // 5d lattice for DWF. + FiveDimWilsonFermion::FiveDimWilsonFermion(LatticeGaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + double _mass) : + _FiveDimGrid(&FiveDimGrid), + _FiveDimRedBlackGrid(&FiveDimRedBlackGrid), + _FourDimGrid(&FourDimGrid), + _FourDimRedBlackGrid(&FourDimRedBlackGrid), + Stencil (_FiveDimGrid,npoint,Even,directions,displacements), + StencilEven(_FiveDimRedBlackGrid,npoint,Even,directions,displacements), // source is Even + StencilOdd (_FiveDimRedBlackGrid,npoint,Odd ,directions,displacements), // source is Odd + mass(_mass), + Umu(_FourDimGrid), + UmuEven(_FourDimRedBlackGrid), + UmuOdd (_FourDimRedBlackGrid), + Lebesgue(_FourDimGrid), + LebesgueEvenOdd(_FourDimRedBlackGrid) +{ + // some assertions + assert(FiveDimGrid._ndimension==5); + assert(FourDimGrid._ndimension==4); + + assert(FiveDimRedBlackGrid._ndimension==5); + assert(FourDimRedBlackGrid._ndimension==4); + + assert(FiveDimRedBlackGrid._checker_dim==1); + + // Dimension zero of the five-d is the Ls direction + Ls=FiveDimGrid._fdimensions[0]; + assert(FiveDimRedBlackGrid._fdimensions[0]==Ls); + assert(FiveDimRedBlackGrid._processors[0] ==1); + assert(FiveDimRedBlackGrid._simd_layout[0]==1); + assert(FiveDimGrid._processors[0] ==1); + assert(FiveDimGrid._simd_layout[0] ==1); + + // Other dimensions must match the decomposition of the four-D fields + for(int d=0;d<4;d++){ + assert(FourDimRedBlackGrid._fdimensions[d] ==FourDimGrid._fdimensions[d]); + assert(FiveDimRedBlackGrid._fdimensions[d+1]==FourDimGrid._fdimensions[d]); + + assert(FourDimRedBlackGrid._processors[d] ==FourDimGrid._processors[d]); + assert(FiveDimRedBlackGrid._processors[d+1] ==FourDimGrid._processors[d]); + + assert(FourDimRedBlackGrid._simd_layout[d] ==FourDimGrid._simd_layout[d]); + assert(FiveDimRedBlackGrid._simd_layout[d+1]==FourDimGrid._simd_layout[d]); + + assert(FiveDimGrid._fdimensions[d+1] ==FourDimGrid._fdimensions[d]); + assert(FiveDimGrid._processors[d+1] ==FourDimGrid._processors[d]); + assert(FiveDimGrid._simd_layout[d+1] ==FourDimGrid._simd_layout[d]); + } + + // Allocate the required comms buffer + comm_buf.resize(Stencil._unified_buffer_size); // this is always big enough to contain EO + + DoubleStore(Umu,_Umu); + pickCheckerboard(Even,UmuEven,Umu); + pickCheckerboard(Odd ,UmuOdd,Umu); +} +void FiveDimWilsonFermion::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu) +{ + conformable(Uds._grid,GaugeGrid()); + conformable(Umu._grid,GaugeGrid()); + LatticeColourMatrix U(GaugeGrid()); + for(int mu=0;mu(Umu,mu); + pokeIndex(Uds,U,mu); + U = adj(Cshift(U,mu,-1)); + pokeIndex(Uds,U,mu+4); + } +} + +RealD FiveDimWilsonFermion::M(const LatticeFermion &in, LatticeFermion &out) +{ + out.checkerboard=in.checkerboard; + Dhop(in,out,DaggerNo); + return axpy_norm(out,5.0-M5,in,out); +} +RealD FiveDimWilsonFermion::Mdag(const LatticeFermion &in, LatticeFermion &out) +{ + out.checkerboard=in.checkerboard; + Dhop(in,out,DaggerYes); + return axpy_norm(out,5.0-M5,in,out); +} +void FiveDimWilsonFermion::Meooe(const LatticeFermion &in, LatticeFermion &out) +{ + if ( in.checkerboard == Odd ) { + DhopEO(in,out,DaggerNo); + } else { + DhopOE(in,out,DaggerNo); + } +} +void FiveDimWilsonFermion::MeooeDag(const LatticeFermion &in, LatticeFermion &out) +{ + if ( in.checkerboard == Odd ) { + DhopEO(in,out,DaggerYes); + } else { + DhopOE(in,out,DaggerYes); + } +} +void FiveDimWilsonFermion::Mooee(const LatticeFermion &in, LatticeFermion &out) +{ + out.checkerboard = in.checkerboard; + out = (5.0-M5)*in; + return ; +} +void FiveDimWilsonFermion::MooeeDag(const LatticeFermion &in, LatticeFermion &out) +{ + out.checkerboard = in.checkerboard; + Mooee(in,out); +} +void FiveDimWilsonFermion::MooeeInv(const LatticeFermion &in, LatticeFermion &out) +{ + out.checkerboard = in.checkerboard; + out = (1.0/(5.0-M5))*in; + return ; +} +void FiveDimWilsonFermion::MooeeInvDag(const LatticeFermion &in, LatticeFermion &out) +{ + out.checkerboard = in.checkerboard; + MooeeInv(in,out); +} +void FiveDimWilsonFermion::DhopInternal(CartesianStencil & st, LebesgueOrder &lo, + LatticeDoubledGaugeField & U, + const LatticeFermion &in, LatticeFermion &out,int dag) +{ + assert((dag==DaggerNo) ||(dag==DaggerYes)); + + WilsonCompressor compressor(dag); + + st.HaloExchange(in,comm_buf,compressor); + + // Dhop takes the 4d grid from U, and makes a 5d index for fermion + // Not loop ordering and data layout. + // Designed to create + // - per thread reuse in L1 cache for U + // - 8 linear access unit stride streams per thread for Fermion for hw prefetchable. + if ( dag == DaggerYes ) { + if( HandOptDslash ) { + for(int ss=0;ssoSites();ss++){ + int sU=lo.Reorder(ss); +PARALLEL_FOR_LOOP + for(int s=0;soSites();ss++){ + int sU=lo.Reorder(ss); +PARALLEL_FOR_LOOP + for(int s=0;soSites();ss++){ + int sU=lo.Reorder(ss); + for(int s=0;soSites();ss++){ + int sU=lo.Reorder(ss); +PARALLEL_FOR_LOOP + for(int s=0;s + { + public: + /////////////////////////////////////////////////////////////// + // Implement the abstract base + /////////////////////////////////////////////////////////////// + GridBase *GaugeGrid(void) { return _FourDimGrid ;} + GridBase *GaugeRedBlackGrid(void) { return _FourDimRedBlackGrid ;} + GridBase *FermionGrid(void) { return _FiveDimGrid;} + GridBase *FermionRedBlackGrid(void) { return _FiveDimRedBlackGrid;} + + // override multiply + virtual RealD M (const LatticeFermion &in, LatticeFermion &out); + virtual RealD Mdag (const LatticeFermion &in, LatticeFermion &out); + + // half checkerboard operaions + virtual void Meooe (const LatticeFermion &in, LatticeFermion &out); + virtual void MeooeDag (const LatticeFermion &in, LatticeFermion &out); + virtual void Mooee (const LatticeFermion &in, LatticeFermion &out); + virtual void MooeeDag (const LatticeFermion &in, LatticeFermion &out); + virtual void MooeeInv (const LatticeFermion &in, LatticeFermion &out); + virtual void MooeeInvDag (const LatticeFermion &in, LatticeFermion &out); + + // non-hermitian hopping term; half cb or both + void Dhop (const LatticeFermion &in, LatticeFermion &out,int dag); + void DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag); + void DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag); + + /////////////////////////////////////////////////////////////// + // New methods added + /////////////////////////////////////////////////////////////// + void DhopInternal(CartesianStencil & st, + LebesgueOrder &lo, + LatticeDoubledGaugeField &U, + const LatticeFermion &in, + LatticeFermion &out, + int dag); + + // Constructors + FiveDimWilsonFermion(LatticeGaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + double _mass); + + // DoubleStore + void DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu); + + /////////////////////////////////////////////////////////////// + // Data members require to support the functionality + /////////////////////////////////////////////////////////////// + static int HandOptDslash; // these are a temporary hack + + protected: + + // Add these to the support from Wilson + GridBase *_FourDimGrid; + GridBase *_FourDimRedBlackGrid; + GridBase *_FiveDimGrid; + GridBase *_FiveDimRedBlackGrid; + + static const int npoint=8; + static const std::vector directions ; + static const std::vector displacements; + + double M5; + double mass; + int Ls; + + //Defines the stencils for even and odd + CartesianStencil Stencil; + CartesianStencil StencilEven; + CartesianStencil StencilOdd; + + // Copy of the gauge field , with even and odd subsets + LatticeDoubledGaugeField Umu; + LatticeDoubledGaugeField UmuEven; + LatticeDoubledGaugeField UmuOdd; + + LebesgueOrder Lebesgue; + LebesgueOrder LebesgueEvenOdd; + + // Comms buffer + std::vector > comm_buf; + + }; + } +} + +#endif diff --git a/lib/qcd/action/fermion/WilsonCompressor.h b/lib/qcd/action/fermion/WilsonCompressor.h new file mode 100644 index 00000000..d9fe977a --- /dev/null +++ b/lib/qcd/action/fermion/WilsonCompressor.h @@ -0,0 +1,61 @@ +#ifndef GRID_QCD_WILSON_COMPRESSOR_H +#define GRID_QCD_WILSON_COMPRESSOR_H + +namespace Grid { +namespace QCD { + + class WilsonCompressor { + public: + int mu; + int dag; + + WilsonCompressor(int _dag){ + mu=0; + dag=_dag; + assert((dag==0)||(dag==1)); + } + void Point(int p) { + mu=p; + }; + + vHalfSpinColourVector operator () (const vSpinColourVector &in) + { + vHalfSpinColourVector ret; + int mudag=mu; + if (dag) { + mudag=(mu+Nd)%(2*Nd); + } + switch(mudag) { + case Xp: + spProjXp(ret,in); + break; + case Yp: + spProjYp(ret,in); + break; + case Zp: + spProjZp(ret,in); + break; + case Tp: + spProjTp(ret,in); + break; + case Xm: + spProjXm(ret,in); + break; + case Ym: + spProjYm(ret,in); + break; + case Zm: + spProjZm(ret,in); + break; + case Tm: + spProjTm(ret,in); + break; + default: + assert(0); + break; + } + return ret; + } + }; +}} // namespace close +#endif diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc new file mode 100644 index 00000000..aa30a7fa --- /dev/null +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -0,0 +1,163 @@ +#include + +namespace Grid { +namespace QCD { + +const std::vector WilsonFermion::directions ({0,1,2,3, 0, 1, 2, 3}); +const std::vector WilsonFermion::displacements({1,1,1,1,-1,-1,-1,-1}); + +int WilsonFermion::HandOptDslash; + +WilsonFermion::WilsonFermion(LatticeGaugeField &_Umu, + GridCartesian &Fgrid, + GridRedBlackCartesian &Hgrid, + double _mass) : + _grid(&Fgrid), + _cbgrid(&Hgrid), + Stencil (&Fgrid,npoint,Even,directions,displacements), + StencilEven(&Hgrid,npoint,Even,directions,displacements), // source is Even + StencilOdd (&Hgrid,npoint,Odd ,directions,displacements), // source is Odd + mass(_mass), + Umu(&Fgrid), + UmuEven(&Hgrid), + UmuOdd (&Hgrid) +{ + // Allocate the required comms buffer + comm_buf.resize(Stencil._unified_buffer_size); // this is always big enough to contain EO + DoubleStore(Umu,_Umu); + pickCheckerboard(Even,UmuEven,Umu); + pickCheckerboard(Odd ,UmuOdd,Umu); +} + +void WilsonFermion::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu) +{ + conformable(Uds._grid,GaugeGrid()); + conformable(Umu._grid,GaugeGrid()); + LatticeColourMatrix U(GaugeGrid()); + for(int mu=0;mu(Umu,mu); + pokeIndex(Uds,U,mu); + U = adj(Cshift(U,mu,-1)); + pokeIndex(Uds,U,mu+4); + } +} + +RealD WilsonFermion::M(const LatticeFermion &in, LatticeFermion &out) +{ + out.checkerboard=in.checkerboard; + Dhop(in,out,DaggerNo); + return axpy_norm(out,4+mass,in,out); +} +RealD WilsonFermion::Mdag(const LatticeFermion &in, LatticeFermion &out) +{ + out.checkerboard=in.checkerboard; + Dhop(in,out,DaggerYes); + return axpy_norm(out,4+mass,in,out); +} + +void WilsonFermion::Meooe(const LatticeFermion &in, LatticeFermion &out) +{ + if ( in.checkerboard == Odd ) { + DhopEO(in,out,DaggerNo); + } else { + DhopOE(in,out,DaggerNo); + } +} +void WilsonFermion::MeooeDag(const LatticeFermion &in, LatticeFermion &out) +{ + if ( in.checkerboard == Odd ) { + DhopEO(in,out,DaggerYes); + } else { + DhopOE(in,out,DaggerYes); + } +} +void WilsonFermion::Mooee(const LatticeFermion &in, LatticeFermion &out) +{ + out.checkerboard = in.checkerboard; + out = (4.0+mass)*in; + return ; +} +void WilsonFermion::MooeeDag(const LatticeFermion &in, LatticeFermion &out) +{ + out.checkerboard = in.checkerboard; + Mooee(in,out); +} +void WilsonFermion::MooeeInv(const LatticeFermion &in, LatticeFermion &out) +{ + out.checkerboard = in.checkerboard; + out = (1.0/(4.0+mass))*in; + return ; +} +void WilsonFermion::MooeeInvDag(const LatticeFermion &in, LatticeFermion &out) +{ + out.checkerboard = in.checkerboard; + MooeeInv(in,out); +} + +void WilsonFermion::DhopInternal(CartesianStencil & st,LatticeDoubledGaugeField & U, + const LatticeFermion &in, LatticeFermion &out,int dag) +{ + assert((dag==DaggerNo) ||(dag==DaggerYes)); + WilsonCompressor compressor(dag); + st.HaloExchange(in,comm_buf,compressor); + + if ( dag == DaggerYes ) { + if( HandOptDslash ) { +PARALLEL_FOR_LOOP + for(int sss=0;sssoSites();sss++){ + DiracOptHand::DhopSiteDag(st,U,comm_buf,sss,sss,in,out); + } + } else { +PARALLEL_FOR_LOOP + for(int sss=0;sssoSites();sss++){ + DiracOpt::DhopSiteDag(st,U,comm_buf,sss,sss,in,out); + } + } + } else { + if( HandOptDslash ) { +PARALLEL_FOR_LOOP + for(int sss=0;sssoSites();sss++){ + DiracOptHand::DhopSite(st,U,comm_buf,sss,sss,in,out); + } + } else { +PARALLEL_FOR_LOOP + for(int sss=0;sssoSites();sss++){ + DiracOpt::DhopSite(st,U,comm_buf,sss,sss,in,out); + } + } + } +} +void WilsonFermion::DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag) +{ + conformable(in._grid,_cbgrid); // verifies half grid + conformable(in._grid,out._grid); // drops the cb check + + assert(in.checkerboard==Even); + out.checkerboard = Odd; + + DhopInternal(StencilEven,UmuOdd,in,out,dag); +} +void WilsonFermion::DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag) +{ + conformable(in._grid,_cbgrid); // verifies half grid + conformable(in._grid,out._grid); // drops the cb check + + assert(in.checkerboard==Odd); + out.checkerboard = Even; + + DhopInternal(StencilOdd,UmuEven,in,out,dag); +} +void WilsonFermion::Dhop(const LatticeFermion &in, LatticeFermion &out,int dag) +{ + conformable(in._grid,_grid); // verifies full grid + conformable(in._grid,out._grid); + + out.checkerboard = in.checkerboard; + + DhopInternal(Stencil,Umu,in,out,dag); +} + +}} + + + diff --git a/lib/qcd/action/fermion/WilsonFermion.h b/lib/qcd/action/fermion/WilsonFermion.h new file mode 100644 index 00000000..5c208131 --- /dev/null +++ b/lib/qcd/action/fermion/WilsonFermion.h @@ -0,0 +1,87 @@ +#ifndef GRID_QCD_WILSON_FERMION_H +#define GRID_QCD_WILSON_FERMION_H + +namespace Grid { + + namespace QCD { + + class WilsonFermion : public FermionAction + { + public: + + /////////////////////////////////////////////////////////////// + // Implement the abstract base + /////////////////////////////////////////////////////////////// + GridBase *GaugeGrid(void) { return _grid ;} + GridBase *GaugeRedBlackGrid(void) { return _cbgrid ;} + GridBase *FermionGrid(void) { return _grid;} + GridBase *FermionRedBlackGrid(void) { return _cbgrid;} + + // override multiply + virtual RealD M (const LatticeFermion &in, LatticeFermion &out); + virtual RealD Mdag (const LatticeFermion &in, LatticeFermion &out); + + // half checkerboard operaions + void Meooe (const LatticeFermion &in, LatticeFermion &out); + void MeooeDag (const LatticeFermion &in, LatticeFermion &out); + virtual void Mooee (const LatticeFermion &in, LatticeFermion &out); // remain virtual so we + virtual void MooeeDag (const LatticeFermion &in, LatticeFermion &out); // can derive Clover + virtual void MooeeInv (const LatticeFermion &in, LatticeFermion &out); // from Wilson bas + virtual void MooeeInvDag (const LatticeFermion &in, LatticeFermion &out); + + // non-hermitian hopping term; half cb or both + void Dhop (const LatticeFermion &in, LatticeFermion &out,int dag); + void DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag); + void DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag); + + /////////////////////////////////////////////////////////////// + // Extra methods added by derived + /////////////////////////////////////////////////////////////// + void DhopInternal(CartesianStencil & st, + LatticeDoubledGaugeField &U, + const LatticeFermion &in, + LatticeFermion &out, + int dag); + + // Constructor + WilsonFermion(LatticeGaugeField &_Umu,GridCartesian &Fgrid,GridRedBlackCartesian &Hgrid,double _mass); + + // DoubleStore + void DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu); + + /////////////////////////////////////////////////////////////// + // Data members require to support the functionality + /////////////////////////////////////////////////////////////// + static int HandOptDslash; // these are a temporary hack + static int MortonOrder; + + protected: + + double mass; + + GridBase * _grid; + GridBase * _cbgrid; + + static const int npoint=8; + static const std::vector directions ; + static const std::vector displacements; + + //Defines the stencils for even and odd + CartesianStencil Stencil; + CartesianStencil StencilEven; + CartesianStencil StencilOdd; + + // Copy of the gauge field , with even and odd subsets + LatticeDoubledGaugeField Umu; + LatticeDoubledGaugeField UmuEven; + LatticeDoubledGaugeField UmuOdd; + + // Comms buffer + std::vector > comm_buf; + + + }; + + } +} +#endif diff --git a/lib/qcd/Grid_qcd_dhop.cc b/lib/qcd/action/fermion/WilsonKernels.cc similarity index 87% rename from lib/qcd/Grid_qcd_dhop.cc rename to lib/qcd/action/fermion/WilsonKernels.cc index 1e5dcd16..879bbb7a 100644 --- a/lib/qcd/Grid_qcd_dhop.cc +++ b/lib/qcd/action/fermion/WilsonKernels.cc @@ -4,8 +4,8 @@ namespace Grid { namespace QCD { void DiracOpt::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, - std::vector > &buf, - int ss,const LatticeFermion &in, LatticeFermion &out) + std::vector > &buf, + int sF,int sU,const LatticeFermion &in, LatticeFermion &out) { vHalfSpinColourVector tmp; vHalfSpinColourVector chi; @@ -16,6 +16,7 @@ void DiracOpt::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, //#define VERBOSE( A) if ( ss<10 ) { std::cout << "site " < > &buf, - int ss,const LatticeFermion &in, LatticeFermion &out) + std::vector > &buf, + int sF,int sU,const LatticeFermion &in, LatticeFermion &out) { vHalfSpinColourVector tmp; vHalfSpinColourVector chi; @@ -173,6 +174,7 @@ void DiracOpt::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U, int offset,local,perm, ptype; // Xp + int ss=sF; offset = st._offsets [Xm][ss]; local = st._is_local[Xm][ss]; perm = st._permute[Xm][ss]; @@ -186,7 +188,7 @@ void DiracOpt::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[ss](Xm),&chi()); + mult(&Uchi(),&U._odata[sU](Xm),&chi()); spReconXp(result,Uchi); // Yp @@ -202,7 +204,7 @@ void DiracOpt::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[ss](Ym),&chi()); + mult(&Uchi(),&U._odata[sU](Ym),&chi()); accumReconYp(result,Uchi); // Zp @@ -218,7 +220,7 @@ void DiracOpt::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[ss](Zm),&chi()); + mult(&Uchi(),&U._odata[sU](Zm),&chi()); accumReconZp(result,Uchi); // Tp @@ -234,7 +236,7 @@ void DiracOpt::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[ss](Tm),&chi()); + mult(&Uchi(),&U._odata[sU](Tm),&chi()); accumReconTp(result,Uchi); // Xm @@ -252,7 +254,7 @@ void DiracOpt::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[ss](Xp),&chi()); + mult(&Uchi(),&U._odata[sU](Xp),&chi()); accumReconXm(result,Uchi); // Ym @@ -269,7 +271,7 @@ void DiracOpt::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[ss](Yp),&chi()); + mult(&Uchi(),&U._odata[sU](Yp),&chi()); accumReconYm(result,Uchi); // Zm @@ -285,7 +287,7 @@ void DiracOpt::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[ss](Zp),&chi()); + mult(&Uchi(),&U._odata[sU](Zp),&chi()); accumReconZm(result,Uchi); // Tm @@ -301,9 +303,9 @@ void DiracOpt::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[ss](Tp),&chi()); + mult(&Uchi(),&U._odata[sU](Tp),&chi()); accumReconTm(result,Uchi); - vstream(out._odata[ss],result); + vstream(out._odata[ss],result*(-0.5)); } }} diff --git a/lib/qcd/action/fermion/WilsonKernels.h b/lib/qcd/action/fermion/WilsonKernels.h new file mode 100644 index 00000000..cc535022 --- /dev/null +++ b/lib/qcd/action/fermion/WilsonKernels.h @@ -0,0 +1,42 @@ +#ifndef GRID_QCD_DHOP_H +#define GRID_QCD_DHOP_H + +namespace Grid { + + namespace QCD { + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // Helper classes that implement Wilson stencil for a single site. + //////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + // Generic version works for any Nc and with extra flavour indices + class DiracOpt { + public: + // These ones will need to be package intelligently. WilsonType base class + // for use by DWF etc.. + static void DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, + std::vector > &buf, + int sF,int sU,const LatticeFermion &in, LatticeFermion &out); + static void DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U, + std::vector > &buf, + int sF,int sU,const LatticeFermion &in, LatticeFermion &out); + + }; + + // Hand unrolled for Nc=3, one flavour + class DiracOptHand { + public: + // These ones will need to be package intelligently. WilsonType base class + // for use by DWF etc.. + static void DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, + std::vector > &buf, + int sF,int sU,const LatticeFermion &in, LatticeFermion &out); + static void DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U, + std::vector > &buf, + int sF,int sU,const LatticeFermion &in, LatticeFermion &out); + + }; + + } +} +#endif diff --git a/lib/qcd/Grid_qcd_dhop_hand.cc b/lib/qcd/action/fermion/WilsonKernelsHand.cc similarity index 92% rename from lib/qcd/Grid_qcd_dhop_hand.cc rename to lib/qcd/action/fermion/WilsonKernelsHand.cc index f8d464fb..019b668b 100644 --- a/lib/qcd/Grid_qcd_dhop_hand.cc +++ b/lib/qcd/action/fermion/WilsonKernelsHand.cc @@ -27,7 +27,7 @@ Chi_12 = ref()(1)(2); #define MULT_2SPIN(A)\ - auto & ref(U._odata[ss](A)); \ + auto & ref(U._odata[sU](A)); \ U_00 = ref()(0,0);\ U_10 = ref()(1,0);\ U_20 = ref()(2,0);\ @@ -282,7 +282,7 @@ namespace QCD { void DiracOptHand::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, std::vector > &buf, - int ss,const LatticeFermion &in, LatticeFermion &out) + int sF,int sU,const LatticeFermion &in, LatticeFermion &out) { REGISTER vComplex result_00; // 12 regs on knc REGISTER vComplex result_01; @@ -338,7 +338,8 @@ void DiracOptHand::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, int offset,local,perm, ptype; - + int ss=sF; + // Xp offset = st._offsets [Xp][ss]; local = st._is_local[Xp][ss]; @@ -514,24 +515,24 @@ void DiracOptHand::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, { vSpinColourVector & ref (out._odata[ss]); - vstream(ref()(0)(0),result_00); - vstream(ref()(0)(1),result_01); - vstream(ref()(0)(2),result_02); - vstream(ref()(1)(0),result_10); - vstream(ref()(1)(1),result_11); - vstream(ref()(1)(2),result_12); - vstream(ref()(2)(0),result_20); - vstream(ref()(2)(1),result_21); - vstream(ref()(2)(2),result_22); - vstream(ref()(3)(0),result_30); - vstream(ref()(3)(1),result_31); - vstream(ref()(3)(2),result_32); + vstream(ref()(0)(0),result_00*(-0.5)); + vstream(ref()(0)(1),result_01*(-0.5)); + vstream(ref()(0)(2),result_02*(-0.5)); + vstream(ref()(1)(0),result_10*(-0.5)); + vstream(ref()(1)(1),result_11*(-0.5)); + vstream(ref()(1)(2),result_12*(-0.5)); + vstream(ref()(2)(0),result_20*(-0.5)); + vstream(ref()(2)(1),result_21*(-0.5)); + vstream(ref()(2)(2),result_22*(-0.5)); + vstream(ref()(3)(0),result_30*(-0.5)); + vstream(ref()(3)(1),result_31*(-0.5)); + vstream(ref()(3)(2),result_32*(-0.5)); } } void DiracOptHand::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U, std::vector > &buf, - int ss,const LatticeFermion &in, LatticeFermion &out) + int ss,int sU,const LatticeFermion &in, LatticeFermion &out) { REGISTER vComplex result_00; // 12 regs on knc REGISTER vComplex result_01; @@ -752,18 +753,18 @@ void DiracOptHand::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U, { vSpinColourVector & ref (out._odata[ss]); - vstream(ref()(0)(0),result_00); - vstream(ref()(0)(1),result_01); - vstream(ref()(0)(2),result_02); - vstream(ref()(1)(0),result_10); - vstream(ref()(1)(1),result_11); - vstream(ref()(1)(2),result_12); - vstream(ref()(2)(0),result_20); - vstream(ref()(2)(1),result_21); - vstream(ref()(2)(2),result_22); - vstream(ref()(3)(0),result_30); - vstream(ref()(3)(1),result_31); - vstream(ref()(3)(2),result_32); + vstream(ref()(0)(0),result_00*(-0.5)); + vstream(ref()(0)(1),result_01*(-0.5)); + vstream(ref()(0)(2),result_02*(-0.5)); + vstream(ref()(1)(0),result_10*(-0.5)); + vstream(ref()(1)(1),result_11*(-0.5)); + vstream(ref()(1)(2),result_12*(-0.5)); + vstream(ref()(2)(0),result_20*(-0.5)); + vstream(ref()(2)(1),result_21*(-0.5)); + vstream(ref()(2)(2),result_22*(-0.5)); + vstream(ref()(3)(0),result_30*(-0.5)); + vstream(ref()(3)(1),result_31*(-0.5)); + vstream(ref()(3)(2),result_32*(-0.5)); } } }} diff --git a/lib/stencil/Grid_lebesgue.cc b/lib/stencil/Grid_lebesgue.cc new file mode 100644 index 00000000..977e3562 --- /dev/null +++ b/lib/stencil/Grid_lebesgue.cc @@ -0,0 +1,103 @@ +#include + +namespace Grid { + +int LebesgueOrder::UseLebesgueOrder; + +LebesgueOrder::IndexInteger LebesgueOrder::alignup(IndexInteger n){ + n--; // 1000 0011 --> 1000 0010 + n |= n >> 1; // 1000 0010 | 0100 0001 = 1100 0011 + n |= n >> 2; // 1100 0011 | 0011 0000 = 1111 0011 + n |= n >> 4; // 1111 0011 | 0000 1111 = 1111 1111 + n |= n >> 8; // ... (At this point all bits are 1, so further bitwise-or + n |= n >> 16; // operations produce no effect.) + n++; // 1111 1111 --> 1 0000 0000 + return n; +}; + +LebesgueOrder::LebesgueOrder(GridBase *grid) +{ + _LebesgueReorder.resize(0); + + // Align up dimensions to power of two. + const IndexInteger one=1; + IndexInteger ND = grid->_ndimension; + std::vector dims(ND); + std::vector adims(ND); + std::vector > bitlist(ND); + + for(IndexInteger mu=0;mu_rdimensions[mu]; + assert ( dims[mu] != 0 ); + adims[mu] = alignup(dims[mu]); + } + + // List which bits of padded volume coordinate contribute; this strategy + // i) avoids recursion + // ii) has loop lengths at most the width of a 32 bit word. + int sitebit=0; + int split=2; + for(int mu=0;mu ax(ND); + + for(IndexInteger asite=0;asitedims[mu]-1 ) contained = 0; + + } + + if ( contained ) { + int site = ax[0] + + dims[0]*ax[1] + +dims[0]*dims[1]*ax[2] + +dims[0]*dims[1]*dims[2]*ax[3]; + + _LebesgueReorder.push_back(site); + } + } + assert( _LebesgueReorder.size() == vol ); +} +} diff --git a/lib/stencil/Grid_lebesgue.h b/lib/stencil/Grid_lebesgue.h new file mode 100644 index 00000000..1d59b127 --- /dev/null +++ b/lib/stencil/Grid_lebesgue.h @@ -0,0 +1,29 @@ +#ifndef GRID_LEBESGUE_H +#define GRID_LEBESGUE_H + +#include + +// Lebesgue, Morton, Z-graph ordering assistance +namespace Grid { + + class LebesgueOrder { + public: + + static int UseLebesgueOrder; + + typedef uint32_t IndexInteger; + + inline IndexInteger Reorder(IndexInteger ss) { + return UseLebesgueOrder ? _LebesgueReorder[ss] : ss; + }; + + IndexInteger alignup(IndexInteger n); + + LebesgueOrder(GridBase *grid); + + private: + std::vector _LebesgueReorder; + + }; +} +#endif diff --git a/lib/stencil/Grid_stencil_common.cc b/lib/stencil/Grid_stencil_common.cc index 6cbcc890..f0f8c581 100644 --- a/lib/stencil/Grid_stencil_common.cc +++ b/lib/stencil/Grid_stencil_common.cc @@ -3,95 +3,6 @@ namespace Grid { - -void CartesianStencil::LebesgueOrder(void) -{ - _LebesgueReorder.resize(0); - - // Align up dimensions to power of two. - const StencilInteger one=1; - StencilInteger ND = _grid->_ndimension; - std::vector dims(ND); - std::vector adims(ND); - std::vector > bitlist(ND); - - - for(StencilInteger mu=0;mu_rdimensions[mu]; - assert ( dims[mu] != 0 ); - adims[mu] = alignup(dims[mu]); - } - - // List which bits of padded volume coordinate contribute; this strategy - // i) avoids recursion - // ii) has loop lengths at most the width of a 32 bit word. - int sitebit=0; - int split=24; - for(int mu=0;mu ax(ND); - - for(StencilInteger asite=0;asitedims[mu]-1 ) contained = 0; - - } - - if ( contained ) { - int site = ax[0] - + dims[0]*ax[1] - +dims[0]*dims[1]*ax[2] - +dims[0]*dims[1]*dims[2]*ax[3]; - - _LebesgueReorder.push_back(site); - } - } - - assert( _LebesgueReorder.size() == vol ); -} - CartesianStencil::CartesianStencil(GridBase *grid, int npoints, int checkerboard, @@ -110,8 +21,6 @@ void CartesianStencil::LebesgueOrder(void) _unified_buffer_size=0; _request_count =0; - LebesgueOrder(); - int osites = _grid->oSites(); for(int i=0;iCheckerBoardShift(_checkerboard,dimension,shift,0); - sshift[1] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,1); + sshift[0] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,Even); + sshift[1] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,Odd); if ( sshift[0] == sshift[1] ) { Local(point,dimension,shift,0x3); @@ -154,8 +63,8 @@ void CartesianStencil::LebesgueOrder(void) } } else { // All permute extract done in comms phase prior to Stencil application // So tables are the same whether comm_dim or splice_dim - sshift[0] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,0); - sshift[1] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,1); + sshift[0] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,Even); + sshift[1] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,Odd); if ( sshift[0] == sshift[1] ) { Comms(point,dimension,shift,0x3); } else { @@ -185,7 +94,7 @@ void CartesianStencil::LebesgueOrder(void) int o = 0; int bo = x * _grid->_ostride[dimension]; - int cb= (cbmask==0x2)? 1 : 0; + int cb= (cbmask==0x2)? Odd : Even; int sshift = _grid->CheckerBoardShift(_checkerboard,dimension,shift,cb); int sx = (x+sshift)%rd; @@ -224,7 +133,7 @@ void CartesianStencil::LebesgueOrder(void) _comm_buf_size[point] = buffer_size; // Size of _one_ plane. Multiple planes may be gathered and // send to one or more remote nodes. - int cb= (cbmask==0x2)? 1 : 0; + int cb= (cbmask==0x2)? Odd : Even; int sshift= _grid->CheckerBoardShift(_checkerboard,dimension,shift,cb); for(int x=0;x +#include + +using namespace Grid; +using namespace Grid::QCD; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + int Nd = latt_size.size(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplexF::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + + std::vector mask(Nd,1); + mask[0]=0; + + GridCartesian Fine (latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout,mask,1); + + GridParallelRNG FineRNG(&Fine); FineRNG.SeedRandomDevice(); + + LatticeComplex U(&Fine); + LatticeComplex ShiftU(&Fine); + LatticeComplex rbShiftU(&Fine); + LatticeComplex Ue(&RBFine); + LatticeComplex Uo(&RBFine); + LatticeComplex ShiftUe(&RBFine); + LatticeComplex ShiftUo(&RBFine); + LatticeComplex lex(&Fine); + lex=zero; + Integer stride =1; + { + double nrm; + LatticeComplex coor(&Fine); + + for(int d=0;d coor(4); + + std::cout << "Checking the non-checkerboard shift"< scoor(coor); + scoor[dir] = (scoor[dir]+shift)%latt_size[dir]; + + Integer slex = scoor[0] + + latt_size[0]*scoor[1] + + latt_size[0]*latt_size[1]*scoor[2] + + latt_size[0]*latt_size[1]*latt_size[2]*scoor[3]; + + Complex scm(slex); + + double nrm = abs(scm-cm()()()); + std::vector peer(4); + int index=real(cm); + Fine.CoorFromIndex(peer,index,latt_size); + + if (nrm > 0){ + std::cerr<<"FAIL shift "<< shift<<" in dir "<< dir + <<" ["< scoor(coor); + scoor[dir] = (scoor[dir]+shift)%latt_size[dir]; + + Integer slex = scoor[0] + + latt_size[0]*scoor[1] + + latt_size[0]*latt_size[1]*scoor[2] + + latt_size[0]*latt_size[1]*latt_size[2]*scoor[3]; + + Complex scm(slex); + + nrm = abs(scm-cm()()()); + std::vector peer(4); + int index=real(cm); + Fine.CoorFromIndex(peer,index,latt_size); + + if (nrm > 0){ + std::cerr<<"FAIL shift "<< shift<<" in dir "<< dir + <<" ["<