diff --git a/Makefile.am b/Makefile.am index fc3f6a0a..3b1d5690 100644 --- a/Makefile.am +++ b/Makefile.am @@ -1,3 +1,5 @@ # additional include paths necessary to compile the C++ library AM_CXXFLAGS = -I$(top_srcdir)/ SUBDIRS = lib tests benchmarks + +filelist: $(SUBDIRS) \ No newline at end of file diff --git a/benchmarks/Grid_comms.cc b/benchmarks/Benchmark_comms.cc similarity index 95% rename from benchmarks/Grid_comms.cc rename to benchmarks/Benchmark_comms.cc index 59d709a0..a24e7ed8 100644 --- a/benchmarks/Grid_comms.cc +++ b/benchmarks/Benchmark_comms.cc @@ -77,11 +77,12 @@ int main (int argc, char ** argv) } double stop=usecond(); - double xbytes = Nloop*bytes*2*ncomm; + double dbytes = bytes; + double xbytes = Nloop*dbytes*2.0*ncomm; 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(); + const int Ls=8; + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::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(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"<>$CONFIG_STATUS <<\_ACEOF || ac_write_fail=1 for ac_config_target in $ac_config_targets do case $ac_config_target in - "lib/Grid_config.h") CONFIG_HEADERS="$CONFIG_HEADERS lib/Grid_config.h" ;; + "lib/GridConfig.h") CONFIG_HEADERS="$CONFIG_HEADERS lib/GridConfig.h" ;; "depfiles") CONFIG_COMMANDS="$CONFIG_COMMANDS depfiles" ;; "docs/doxy.cfg") CONFIG_FILES="$CONFIG_FILES docs/doxy.cfg" ;; "Makefile") CONFIG_FILES="$CONFIG_FILES Makefile" ;; diff --git a/configure.ac b/configure.ac index 1edb3e33..74a0eebc 100644 --- a/configure.ac +++ b/configure.ac @@ -11,7 +11,7 @@ AC_CANONICAL_SYSTEM AM_INIT_AUTOMAKE(subdir-objects) AC_CONFIG_MACRO_DIR([m4]) AC_CONFIG_SRCDIR([lib/Grid.h]) -AC_CONFIG_HEADERS([lib/Grid_config.h]) +AC_CONFIG_HEADERS([lib/GridConfig.h]) m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])]) AC_MSG_NOTICE([ diff --git a/tests/InvSqrt.gnu b/lib/.dirstamp similarity index 100% rename from tests/InvSqrt.gnu rename to lib/.dirstamp diff --git a/lib/Grid_algorithms.h b/lib/Algorithms.h similarity index 100% rename from lib/Grid_algorithms.h rename to lib/Algorithms.h diff --git a/lib/Grid_aligned_allocator.h b/lib/AlignedAllocator.h similarity index 100% rename from lib/Grid_aligned_allocator.h rename to lib/AlignedAllocator.h diff --git a/lib/Cartesian.h b/lib/Cartesian.h new file mode 100644 index 00000000..db384b53 --- /dev/null +++ b/lib/Cartesian.h @@ -0,0 +1,8 @@ +#ifndef GRID_CARTESIAN_H +#define GRID_CARTESIAN_H + +#include +#include +#include + +#endif diff --git a/lib/Grid_communicator.h b/lib/Communicator.h similarity index 57% rename from lib/Grid_communicator.h rename to lib/Communicator.h index cfa6e0a7..6880adda 100644 --- a/lib/Grid_communicator.h +++ b/lib/Communicator.h @@ -1,6 +1,6 @@ #ifndef GRID_COMMUNICATOR_H #define GRID_COMMUNICATOR_H -#include +#include #endif diff --git a/lib/Grid_comparison.h b/lib/Comparison.h similarity index 98% rename from lib/Grid_comparison.h rename to lib/Comparison.h index 3f9c206d..ecd6ece0 100644 --- a/lib/Grid_comparison.h +++ b/lib/Comparison.h @@ -141,7 +141,7 @@ namespace Grid { } } -#include -#include +#include +#include #endif diff --git a/lib/Grid_cshift.h b/lib/Cshift.h similarity index 51% rename from lib/Grid_cshift.h rename to lib/Cshift.h index 10c7a3c4..3caccbf9 100644 --- a/lib/Grid_cshift.h +++ b/lib/Cshift.h @@ -1,13 +1,13 @@ #ifndef _GRID_CSHIFT_H_ #define _GRID_CSHIFT_H_ -#include +#include #ifdef GRID_COMMS_NONE -#include +#include #endif #ifdef GRID_COMMS_MPI -#include +#include #endif #endif diff --git a/lib/Grid.h b/lib/Grid.h index 1e513a46..16530434 100644 --- a/lib/Grid.h +++ b/lib/Grid.h @@ -33,7 +33,7 @@ #define strong_inline __attribute__((always_inline)) inline -#include +#include //////////////////////////////////////////////////////////// // Tunable header includes @@ -46,22 +46,22 @@ #include #endif -#include +#include -#include -#include +#include +#include -#include // subdir aggregate -#include // subdir aggregate -#include // subdir aggregate -#include -#include // subdir aggregate -#include // subdir aggregate +#include // subdir aggregate +#include // subdir aggregate +#include // subdir aggregate +#include // subdir aggregate +#include +#include // subdir aggregate +#include // subdir aggregate +#include // subdir aggregate -#include // subdir aggregate - -#include -#include +#include +#include namespace Grid { diff --git a/lib/Grid_config.h b/lib/GridConfig.h similarity index 96% rename from lib/Grid_config.h rename to lib/GridConfig.h index 914bd06c..e5b75382 100644 --- a/lib/Grid_config.h +++ b/lib/GridConfig.h @@ -1,5 +1,5 @@ -/* lib/Grid_config.h. Generated from Grid_config.h.in by configure. */ -/* lib/Grid_config.h.in. Generated from configure.ac by autoheader. */ +/* lib/GridConfig.h. Generated from GridConfig.h.in by configure. */ +/* lib/GridConfig.h.in. Generated from configure.ac by autoheader. */ /* AVX */ /* #undef AVX1 */ diff --git a/lib/Grid_config.h.in b/lib/GridConfig.h.in similarity index 98% rename from lib/Grid_config.h.in rename to lib/GridConfig.h.in index 8c958c82..afa1d1e4 100644 --- a/lib/Grid_config.h.in +++ b/lib/GridConfig.h.in @@ -1,4 +1,4 @@ -/* lib/Grid_config.h.in. Generated from configure.ac by autoheader. */ +/* lib/GridConfig.h.in. Generated from configure.ac by autoheader. */ /* AVX */ #undef AVX1 diff --git a/lib/Grid_init.cc b/lib/GridInit.cc similarity index 97% rename from lib/Grid_init.cc rename to lib/GridInit.cc index f72393cb..20371a4b 100644 --- a/lib/Grid_init.cc +++ b/lib/GridInit.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; + WilsonFermion5D::HandOptDslash=1; + } + if( GridCmdOptionExists(*argv,*argv+*argc,"--lebesgue") ){ + LebesgueOrder::UseLebesgueOrder=1; } GridParseLayout(*argv,*argc, Grid_default_latt, diff --git a/lib/Grid_cartesian.h b/lib/Grid_cartesian.h deleted file mode 100644 index c01be20a..00000000 --- a/lib/Grid_cartesian.h +++ /dev/null @@ -1,8 +0,0 @@ -#ifndef GRID_CARTESIAN_H -#define GRID_CARTESIAN_H - -#include -#include -#include - -#endif diff --git a/lib/Grid_math.h b/lib/Grid_math.h deleted file mode 100644 index 17bc09a5..00000000 --- a/lib/Grid_math.h +++ /dev/null @@ -1,16 +0,0 @@ -#ifndef GRID_MATH_H -#define GRID_MATH_H - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - - -#endif diff --git a/lib/Grid_lattice.h b/lib/Lattice.h similarity index 58% rename from lib/Grid_lattice.h rename to lib/Lattice.h index 35664aee..1f29a908 100644 --- a/lib/Grid_lattice.h +++ b/lib/Lattice.h @@ -1,6 +1,6 @@ #ifndef GRID_LATTICE_H #define GRID_LATTICE_H -#include +#include #endif diff --git a/lib/Make.inc b/lib/Make.inc new file mode 100644 index 00000000..787539b7 --- /dev/null +++ b/lib/Make.inc @@ -0,0 +1,4 @@ + +HFILES=./Cshift.h ./simd/Grid_avx.h ./simd/Grid_vector_types.h ./simd/Grid_sse4.h ./simd/Grid_avx512.h ./simd/Old/Grid_vRealD.h ./simd/Old/Grid_vComplexD.h ./simd/Old/Grid_vInteger.h ./simd/Old/Grid_vComplexF.h ./simd/Old/Grid_vRealF.h ./simd/Grid_qpx.h ./Tensors.h ./Algorithms.h ./communicator/Communicator_base.h ./lattice/Lattice_rng.h ./lattice/Lattice_reduction.h ./lattice/Lattice_transfer.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_comparison.h ./lattice/Lattice_overload.h ./lattice/Lattice_reality.h ./lattice/Lattice_local.h ./lattice/Lattice_conformable.h ./lattice/Lattice_where.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_ET.h ./lattice/Lattice_transpose.h ./lattice/Lattice_trace.h ./Stencil.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_poke.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_class.h ./tensors/Tensor_transpose.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_reality.h ./tensors/Tensor_trace.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_outer.h ./tensors/Tensor_inner.h ./tensors/Tensor_traits.h ./tensors/Tensor_Ta.h ./tensors/Tensor_peek.h ./tensors/Tensor_arith.h ./tensors/Tensor_extract_merge.h ./Communicator.h ./Cartesian.h ./parallelIO/NerscIO.h ./qcd/QCD.h ./qcd/SpaceTimeGrid.h ./qcd/LinalgUtils.h ./qcd/TwoSpinor.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/Dirac.h ./cshift/Cshift_common.h ./cshift/Cshift_none.h ./cshift/Cshift_mpi.h ./Simd.h ./GridConfig.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_red_black.h ./cartesian/Cartesian_full.h ./AlignedAllocator.h ./Lattice.h ./Threads.h ./Comparison.h ./Grid.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Zolotarev.h ./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Remez.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./stencil/Lebesgue.h + +CCFILES=./qcd/SpaceTimeGrid.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/Dirac.cc ./GridInit.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc diff --git a/lib/Makefile.am b/lib/Makefile.am index a5b89c0a..7ee52e0e 100644 --- a/lib/Makefile.am +++ b/lib/Makefile.am @@ -3,100 +3,26 @@ AM_CXXFLAGS = -I$(top_srcdir)/ extra_sources= if BUILD_COMMS_MPI - extra_sources+=communicator/Grid_communicator_mpi.cc + extra_sources+=communicator/Communicator_mpi.cc endif if BUILD_COMMS_NONE - extra_sources+=communicator/Grid_communicator_none.cc + extra_sources+=communicator/Communicator_none.cc endif # # Libraries # -lib_LIBRARIES = libGrid.a -libGrid_a_SOURCES = \ - Grid_init.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 \ - $(extra_sources) +include Make.inc + +lib_LIBRARIES = libGrid.a +libGrid_a_SOURCES = $(CCFILES) $(extra_sources) + + +# qcd/action/fermion/PartialFractionFermion5D.cc\ \ # # 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=$(HFILES) diff --git a/lib/Grid_simd.h b/lib/Simd.h similarity index 100% rename from lib/Grid_simd.h rename to lib/Simd.h diff --git a/lib/Grid_stencil.h b/lib/Stencil.h similarity index 91% rename from lib/Grid_stencil.h rename to lib/Stencil.h index fa22361b..8529e73a 100644 --- a/lib/Grid_stencil.h +++ b/lib/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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#endif diff --git a/lib/Grid_threads.h b/lib/Threads.h similarity index 97% rename from lib/Grid_threads.h rename to lib/Threads.h index 24581855..f453b860 100644 --- a/lib/Grid_threads.h +++ b/lib/Threads.h @@ -9,7 +9,7 @@ #ifdef GRID_OMP #include -#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/algorithms/LinearOperator.h b/lib/algorithms/LinearOperator.h index cb6240ba..5839f2c1 100644 --- a/lib/algorithms/LinearOperator.h +++ b/lib/algorithms/LinearOperator.h @@ -125,39 +125,7 @@ namespace Grid { }; */ - // Chroma interface defining GaugeAction - /* - template class GaugeAction - virtual const CreateGaugeState& getCreateState() const = 0; - virtual GaugeState* createState(const Q& q) const - virtual const GaugeBC& getGaugeBC() const - virtual const Set& getSet(void) const = 0; - virtual void deriv(P& result, const Handle< GaugeState >& state) const - virtual Double S(const Handle< GaugeState >& state) const = 0; - class LinearGaugeAction : public GaugeAction< multi1d, multi1d > - typedef multi1d P; - typedef multi1d Q; - virtual void staple(LatticeColorMatrix& result, - const Handle< GaugeState >& state, - int mu, int cb) const = 0; - */ - - // Chroma interface defining FermionAction - /* - template class FermAct4D : public FermionAction - virtual LinearOperator* linOp(Handle< FermState > state) const = 0; - virtual LinearOperator* lMdagM(Handle< FermState > state) const = 0; - virtual LinOpSystemSolver* invLinOp(Handle< FermState > state, - virtual MdagMSystemSolver* invMdagM(Handle< FermState > state, - virtual LinOpMultiSystemSolver* mInvLinOp(Handle< FermState > state, - virtual MdagMMultiSystemSolver* mInvMdagM(Handle< FermState > state, - virtual MdagMMultiSystemSolverAccumulate* mInvMdagMAcc(Handle< FermState > state, - virtual SystemSolver* qprop(Handle< FermState > state, - class DiffFermAct4D : public FermAct4D - virtual DiffLinearOperator* linOp(Handle< FermState > state) const = 0; - virtual DiffLinearOperator* lMdagM(Handle< FermState > state) const = 0; - */ } #endif 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/approx/Zolotarev.cc b/lib/algorithms/approx/Zolotarev.cc index 7629bbde..c73a3436 100644 --- a/lib/algorithms/approx/Zolotarev.cc +++ b/lib/algorithms/approx/Zolotarev.cc @@ -58,6 +58,8 @@ /* Compute the partial fraction expansion coefficients (alpha) from the * factored form */ +namespace Grid { +namespace Approx { static void construct_partfrac(izd *z) { int dn = z -> dn, dd = z -> dd, type = z -> type; @@ -291,7 +293,7 @@ static void sncndnFK(INTERNAL_PRECISION u, INTERNAL_PRECISION k, * Set type = 0 for the Zolotarev approximation, which is zero at x = 0, and * type = 1 for the approximation which is infinite at x = 0. */ -zolotarev_data* bfm_zolotarev(PRECISION epsilon, int n, int type) { +zolotarev_data* grid_zolotarev(PRECISION epsilon, int n, int type) { INTERNAL_PRECISION A, c, cp, kp, ksq, sn, cn, dn, Kp, Kj, z, z0, t, M, F, l, invlambda, xi, xisq, *tv, s, opl; int m, czero, ts; @@ -412,7 +414,7 @@ zolotarev_data* bfm_zolotarev(PRECISION epsilon, int n, int type) { return zd; } -zolotarev_data* bfm_higham(PRECISION epsilon, int n) { +zolotarev_data* grid_higham(PRECISION epsilon, int n) { INTERNAL_PRECISION A, M, c, cp, z, z0, t, epssq; int m, czero; zolotarev_data *zd; @@ -502,6 +504,7 @@ zolotarev_data* bfm_higham(PRECISION epsilon, int n) { free(d); return zd; } +}} #ifdef TEST @@ -707,4 +710,6 @@ int main(int argc, char** argv) { return EXIT_SUCCESS; } + + #endif /* TEST */ diff --git a/lib/algorithms/approx/Zolotarev.h b/lib/algorithms/approx/Zolotarev.h index 3f0dc58e..869e5a89 100644 --- a/lib/algorithms/approx/Zolotarev.h +++ b/lib/algorithms/approx/Zolotarev.h @@ -1,7 +1,8 @@ /* -*- Mode: C; comment-column: 22; fill-column: 79; -*- */ #ifdef __cplusplus -extern "C" { +namespace Grid { +namespace Approx { #endif #define HVERSION Header Time-stamp: <14-OCT-2004 09:26:51.00 adk@MISSCONTRARY> @@ -76,10 +77,10 @@ typedef struct { * zolotarev_data structure. The arguments must satisfy the constraints that * epsilon > 0, n > 0, and type = 0 or 1. */ -ZOLOTAREV_DATA* bfm_higham(PRECISION epsilon, int n) ; -ZOLOTAREV_DATA* bfm_zolotarev(PRECISION epsilon, int n, int type); +ZOLOTAREV_DATA* grid_higham(PRECISION epsilon, int n) ; +ZOLOTAREV_DATA* grid_zolotarev(PRECISION epsilon, int n, int type); #endif #ifdef __cplusplus -} +}} #endif 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/Cartesian_base.h similarity index 95% rename from lib/cartesian/Grid_cartesian_base.h rename to lib/cartesian/Cartesian_base.h index a74773f8..6303e38e 100644 --- a/lib/cartesian/Grid_cartesian_base.h +++ b/lib/cartesian/Cartesian_base.h @@ -2,7 +2,6 @@ #define GRID_CARTESIAN_BASE_H #include -#include namespace Grid{ @@ -21,7 +20,7 @@ public: // Give Lattice access template friend class Lattice; - GridBase(std::vector & processor_grid) : CartesianCommunicator(processor_grid) {}; + GridBase(const std::vector & processor_grid) : CartesianCommunicator(processor_grid) {}; // Physics Grid information. @@ -52,16 +51,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/Cartesian_full.h similarity index 88% rename from lib/cartesian/Grid_cartesian_full.h rename to lib/cartesian/Cartesian_full.h index 73bd08b3..2a9e0be8 100644 --- a/lib/cartesian/Grid_cartesian_full.h +++ b/lib/cartesian/Cartesian_full.h @@ -18,15 +18,18 @@ 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 CheckerBoardShift(int source_cb,int dim,int shift, int osite){ - return shift; + virtual int CheckerBoardShiftForCB(int source_cb,int dim,int shift, int ocb){ + return shift; } - GridCartesian(std::vector &dimensions, - std::vector &simd_layout, - std::vector &processor_grid + virtual int CheckerBoardShift(int source_cb,int dim,int shift, int osite){ + return shift; + } + GridCartesian(const std::vector &dimensions, + const std::vector &simd_layout, + const std::vector &processor_grid ) : GridBase(processor_grid) { /////////////////////// diff --git a/lib/cartesian/Cartesian_red_black.h b/lib/cartesian/Cartesian_red_black.h new file mode 100644 index 00000000..3a84ed49 --- /dev/null +++ b/lib/cartesian/Cartesian_red_black.h @@ -0,0 +1,196 @@ +#ifndef GRID_CARTESIAN_RED_BLACK_H +#define GRID_CARTESIAN_RED_BLACK_H + + +namespace Grid { + + static const int CbRed =0; + 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; + +// Specialise this for red black grids storing half the data like a chess board. +class GridRedBlackCartesian : public GridBase +{ +public: + std::vector _checker_dim_mask; + int _checker_dim; + + virtual int CheckerBoarded(int dim){ + if( dim==_checker_dim) return 1; + else return 0; + } + virtual int CheckerBoard(std::vector site){ + 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 CheckerBoardShiftForCB(int source_cb,int dim,int shift,int ocb){ + if(dim != _checker_dim) return shift; + + 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. + + 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){ + + 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(const GridBase *base) : GridRedBlackCartesian(base->_fdimensions,base->_simd_layout,base->_processors) {}; + + GridRedBlackCartesian(const std::vector &dimensions, + const std::vector &simd_layout, + const std::vector &processor_grid, + const std::vector &checker_dim_mask, + int checker_dim + ) : GridBase(processor_grid) + { + Init(dimensions,simd_layout,processor_grid,checker_dim_mask,checker_dim); + } + GridRedBlackCartesian(const std::vector &dimensions, + const std::vector &simd_layout, + const 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(const std::vector &dimensions, + const std::vector &simd_layout, + const std::vector &processor_grid, + const std::vector &checker_dim_mask, + int checker_dim) + { + /////////////////////// + // Grid information + /////////////////////// + _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==_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=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; + }; + +}; + +} +#endif diff --git a/lib/cartesian/Grid_cartesian_red_black.h b/lib/cartesian/Grid_cartesian_red_black.h deleted file mode 100644 index ff2d3ba8..00000000 --- a/lib/cartesian/Grid_cartesian_red_black.h +++ /dev/null @@ -1,134 +0,0 @@ -#ifndef GRID_CARTESIAN_RED_BLACK_H -#define GRID_CARTESIAN_RED_BLACK_H - - -namespace Grid { - - static const int CbRed =0; - static const int CbBlack=1; - static const int Even =CbRed; - static const int Odd =CbBlack; - static const int DaggerNo=0; - static const int DaggerYes=1; - -// Specialise this for red black grids storing half the data like a chess board. -class GridRedBlackCartesian : public GridBase -{ -public: - virtual int CheckerBoarded(int dim){ - if( dim==0) return 1; - else return 0; - } - virtual int CheckerBoard(std::vector site){ - return (site[0]+site[1]+site[2]+site[3])&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){ - - if(dim != 0) return shift; - - int fulldim =_fdimensions[0]; - 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 CheckerBoardDestination(int source_cb,int shift){ - if ((shift+_fdimensions[0])&0x1) { - return 1-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) - { - /////////////////////// - // 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; - - 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); - - 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]); - return idx; - }; - -}; - -} -#endif diff --git a/lib/communicator/Grid_communicator_base.h b/lib/communicator/Communicator_base.h similarity index 98% rename from lib/communicator/Grid_communicator_base.h rename to lib/communicator/Communicator_base.h index 47c1f525..61e19993 100644 --- a/lib/communicator/Grid_communicator_base.h +++ b/lib/communicator/Communicator_base.h @@ -27,7 +27,7 @@ class CartesianCommunicator { #endif // Constructor - CartesianCommunicator(std::vector &pdimensions_in); + CartesianCommunicator(const std::vector &pdimensions_in); // Wraps MPI_Cart routines void ShiftedRanks(int dim,int shift,int & source, int & dest); diff --git a/lib/communicator/Grid_communicator_mpi.cc b/lib/communicator/Communicator_mpi.cc similarity index 97% rename from lib/communicator/Grid_communicator_mpi.cc rename to lib/communicator/Communicator_mpi.cc index 6ef05c3d..5dd34705 100644 --- a/lib/communicator/Grid_communicator_mpi.cc +++ b/lib/communicator/Communicator_mpi.cc @@ -5,7 +5,7 @@ namespace Grid { // Should error check all MPI calls. -CartesianCommunicator::CartesianCommunicator(std::vector &processors) +CartesianCommunicator::CartesianCommunicator(const std::vector &processors) { _ndimension = processors.size(); std::vector periodic(_ndimension,1); diff --git a/lib/communicator/Grid_communicator_none.cc b/lib/communicator/Communicator_none.cc similarity index 95% rename from lib/communicator/Grid_communicator_none.cc rename to lib/communicator/Communicator_none.cc index 90eb26cc..d7eb9453 100644 --- a/lib/communicator/Grid_communicator_none.cc +++ b/lib/communicator/Communicator_none.cc @@ -1,7 +1,7 @@ #include "Grid.h" namespace Grid { -CartesianCommunicator::CartesianCommunicator(std::vector &processors) +CartesianCommunicator::CartesianCommunicator(const std::vector &processors) { _processors = processors; _ndimension = processors.size(); diff --git a/lib/cshift/Grid_cshift_common.h b/lib/cshift/Cshift_common.h similarity index 92% rename from lib/cshift/Grid_cshift_common.h rename to lib/cshift/Cshift_common.h index c369fe1c..ffa3f771 100644 --- a/lib/cshift/Grid_cshift_common.h +++ b/lib/cshift/Cshift_common.h @@ -153,7 +153,7 @@ PARALLEL_NESTED_LOOP2 ////////////////////////////////////////////////////// // local to node block strided copies ////////////////////////////////////////////////////// -template void Copy_plane(Lattice& lhs,Lattice &rhs, int dimension,int lplane,int rplane,int cbmask) +template void Copy_plane(Lattice& lhs,const Lattice &rhs, int dimension,int lplane,int rplane,int cbmask) { int rd = rhs._grid->_rdimensions[dimension]; @@ -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]); - } } @@ -183,7 +182,7 @@ PARALLEL_NESTED_LOOP2 } -template void Copy_plane_permute(Lattice& lhs,Lattice &rhs, int dimension,int lplane,int rplane,int cbmask,int permute_type) +template void Copy_plane_permute(Lattice& lhs,const Lattice &rhs, int dimension,int lplane,int rplane,int cbmask,int permute_type) { int rd = rhs._grid->_rdimensions[dimension]; @@ -213,12 +212,12 @@ PARALLEL_NESTED_LOOP2 ////////////////////////////////////////////////////// // Local to node Cshift ////////////////////////////////////////////////////// -template void Cshift_local(Lattice& ret,Lattice &rhs,int dimension,int shift) +template void Cshift_local(Lattice& ret,const Lattice &rhs,int dimension,int shift) { 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); @@ -228,7 +227,7 @@ template void Cshift_local(Lattice& ret,Lattice &rhs,int } } -template Lattice Cshift_local(Lattice &ret,Lattice &rhs,int dimension,int shift,int cbmask) +template Lattice Cshift_local(Lattice &ret,const Lattice &rhs,int dimension,int shift,int cbmask) { GridBase *grid = rhs._grid; int fd = grid->_fdimensions[dimension]; @@ -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/Cshift_mpi.h similarity index 85% rename from lib/cshift/Grid_cshift_mpi.h rename to lib/cshift/Cshift_mpi.h index 569017d8..9cdfd316 100644 --- a/lib/cshift/Grid_cshift_mpi.h +++ b/lib/cshift/Cshift_mpi.h @@ -4,7 +4,7 @@ namespace Grid { -template Lattice Cshift(Lattice &rhs,int dimension,int shift) +template Lattice Cshift(const Lattice &rhs,int dimension,int shift) { typedef typename vobj::vector_type vector_type; typedef typename vobj::scalar_type scalar_type; @@ -17,7 +17,7 @@ template Lattice Cshift(Lattice &rhs,int dimension,int s // Map to always positive shift modulo global full dimension. shift = (shift+fd)%fd; - ret.checkerboard = rhs._grid->CheckerBoardDestination(rhs.checkerboard,shift); + ret.checkerboard = rhs._grid->CheckerBoardDestination(rhs.checkerboard,shift,dimension); // the permute type int simd_layout = rhs._grid->_simd_layout[dimension]; @@ -35,12 +35,12 @@ template Lattice Cshift(Lattice &rhs,int dimension,int s return ret; } -template void Cshift_comms(Lattice& ret,Lattice &rhs,int dimension,int shift) +template void Cshift_comms(Lattice& ret,const Lattice &rhs,int dimension,int shift) { 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); @@ -50,12 +50,12 @@ template void Cshift_comms(Lattice& ret,Lattice &rhs,int } } -template void Cshift_comms_simd(Lattice& ret,Lattice &rhs,int dimension,int shift) +template void Cshift_comms_simd(Lattice& ret,const Lattice &rhs,int dimension,int shift) { 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); @@ -65,7 +65,7 @@ template void Cshift_comms_simd(Lattice& ret,Lattice &rh } } -template void Cshift_comms(Lattice &ret,Lattice &rhs,int dimension,int shift,int cbmask) +template void Cshift_comms(Lattice &ret,const Lattice &rhs,int dimension,int shift,int cbmask) { typedef typename vobj::vector_type vector_type; typedef typename vobj::scalar_type scalar_type; @@ -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(Lattice &ret,Lattice &rhs,int } } -template void Cshift_comms_simd(Lattice &ret,Lattice &rhs,int dimension,int shift,int cbmask) +template void Cshift_comms_simd(Lattice &ret,const Lattice &rhs,int dimension,int shift,int cbmask) { GridBase *grid=rhs._grid; const int Nsimd = grid->Nsimd(); @@ -162,8 +162,8 @@ template 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) +template Lattice Cshift(const 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_ET.h b/lib/lattice/Lattice_ET.h similarity index 100% rename from lib/lattice/Grid_lattice_ET.h rename to lib/lattice/Lattice_ET.h diff --git a/lib/lattice/Grid_lattice_arith.h b/lib/lattice/Lattice_arith.h similarity index 100% rename from lib/lattice/Grid_lattice_arith.h rename to lib/lattice/Lattice_arith.h diff --git a/lib/lattice/Grid_lattice_base.h b/lib/lattice/Lattice_base.h similarity index 92% rename from lib/lattice/Grid_lattice_base.h rename to lib/lattice/Lattice_base.h index 4a6d3180..6b5fe261 100644 --- a/lib/lattice/Grid_lattice_base.h +++ b/lib/lattice/Lattice_base.h @@ -176,7 +176,7 @@ PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP for(int ss=0;ss<_grid->oSites();ss++){ #ifdef STREAMING_STORES - vobj tmp = eval(tmp,ss,expr); + vobj tmp = eval(ss,expr); vstream(_odata[ss] ,tmp); #else _odata[ss]=eval(ss,expr); @@ -283,24 +283,23 @@ PARALLEL_FOR_LOOP -#include +#include #define GRID_LATTICE_EXPRESSION_TEMPLATES #ifdef GRID_LATTICE_EXPRESSION_TEMPLATES -#include +#include #else -#include +#include #endif -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include diff --git a/lib/lattice/Grid_lattice_comparison.h b/lib/lattice/Lattice_comparison.h similarity index 100% rename from lib/lattice/Grid_lattice_comparison.h rename to lib/lattice/Lattice_comparison.h diff --git a/lib/lattice/Grid_lattice_conformable.h b/lib/lattice/Lattice_conformable.h similarity index 100% rename from lib/lattice/Grid_lattice_conformable.h rename to lib/lattice/Lattice_conformable.h diff --git a/lib/lattice/Grid_lattice_coordinate.h b/lib/lattice/Lattice_coordinate.h similarity index 100% rename from lib/lattice/Grid_lattice_coordinate.h rename to lib/lattice/Lattice_coordinate.h diff --git a/lib/lattice/Grid_lattice_local.h b/lib/lattice/Lattice_local.h similarity index 100% rename from lib/lattice/Grid_lattice_local.h rename to lib/lattice/Lattice_local.h diff --git a/lib/lattice/Grid_lattice_overload.h b/lib/lattice/Lattice_overload.h similarity index 100% rename from lib/lattice/Grid_lattice_overload.h rename to lib/lattice/Lattice_overload.h diff --git a/lib/lattice/Grid_lattice_peekpoke.h b/lib/lattice/Lattice_peekpoke.h similarity index 98% rename from lib/lattice/Grid_lattice_peekpoke.h rename to lib/lattice/Lattice_peekpoke.h index 7bffdbb8..9416226a 100644 --- a/lib/lattice/Grid_lattice_peekpoke.h +++ b/lib/lattice/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_reality.h b/lib/lattice/Lattice_reality.h similarity index 100% rename from lib/lattice/Grid_lattice_reality.h rename to lib/lattice/Lattice_reality.h diff --git a/lib/lattice/Grid_lattice_reduction.h b/lib/lattice/Lattice_reduction.h similarity index 100% rename from lib/lattice/Grid_lattice_reduction.h rename to lib/lattice/Lattice_reduction.h diff --git a/lib/lattice/Grid_lattice_rng.h b/lib/lattice/Lattice_rng.h similarity index 100% rename from lib/lattice/Grid_lattice_rng.h rename to lib/lattice/Lattice_rng.h diff --git a/lib/lattice/Grid_lattice_trace.h b/lib/lattice/Lattice_trace.h similarity index 100% rename from lib/lattice/Grid_lattice_trace.h rename to lib/lattice/Lattice_trace.h diff --git a/lib/lattice/Grid_lattice_transfer.h b/lib/lattice/Lattice_transfer.h similarity index 99% rename from lib/lattice/Grid_lattice_transfer.h rename to lib/lattice/Lattice_transfer.h index 27ae27e1..f0046b0e 100644 --- a/lib/lattice/Grid_lattice_transfer.h +++ b/lib/lattice/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/lattice/Grid_lattice_transpose.h b/lib/lattice/Lattice_transpose.h similarity index 100% rename from lib/lattice/Grid_lattice_transpose.h rename to lib/lattice/Lattice_transpose.h diff --git a/lib/lattice/Grid_lattice_where.h b/lib/lattice/Lattice_where.h similarity index 100% rename from lib/lattice/Grid_lattice_where.h rename to lib/lattice/Lattice_where.h diff --git a/lib/math/Grid_math_arith.h b/lib/math/Grid_math_arith.h deleted file mode 100644 index ca90ba88..00000000 --- a/lib/math/Grid_math_arith.h +++ /dev/null @@ -1,11 +0,0 @@ -#ifndef GRID_MATH_ARITH_H -#define GRID_MATH_ARITH_H - -#include -#include -#include -#include -#include - -#endif - diff --git a/lib/parallelIO/GridNerscIO.h b/lib/parallelIO/NerscIO.h similarity index 100% rename from lib/parallelIO/GridNerscIO.h rename to lib/parallelIO/NerscIO.h 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/LinalgUtils.h b/lib/qcd/LinalgUtils.h new file mode 100644 index 00000000..2b83d115 --- /dev/null +++ b/lib/qcd/LinalgUtils.h @@ -0,0 +1,113 @@ +#ifndef GRID_QCD_LINALG_UTILS_H +#define GRID_QCD_LINALG_UTILS_H + +namespace Grid{ +namespace QCD{ +//////////////////////////////////////////////////////////////////////// +//This file brings additional linear combination assist that is helpful +//to QCD such as chiral projectors and spin matrices applied to one of the inputs. +//These routines support five-D chiral fermions and contain s-subslice indexing +//on the 5d (rb4d) checkerboarded lattices +//////////////////////////////////////////////////////////////////////// +template +void axpby_ssp(Lattice &z, RealD a,const Lattice &x,RealD b,const Lattice &y,int s,int sp) +{ + z.checkerboard = x.checkerboard; + conformable(x,y); + conformable(x,z); + GridBase *grid=x._grid; + int Ls = grid->_rdimensions[0]; +PARALLEL_FOR_LOOP + for(int ss=0;ssoSites();ss+=Ls){ // adds Ls + vobj tmp = a*x._odata[ss+s]+b*y._odata[ss+sp]; + vstream(z._odata[ss+s],tmp); + } +} + +template +void ag5xpby_ssp(Lattice &z,RealD a,const Lattice &x,RealD b,const Lattice &y,int s,int sp) +{ + z.checkerboard = x.checkerboard; + conformable(x,y); + conformable(x,z); + GridBase *grid=x._grid; + int Ls = grid->_rdimensions[0]; +PARALLEL_FOR_LOOP + for(int ss=0;ssoSites();ss+=Ls){ // adds Ls + vobj tmp; + multGamma5(tmp(),a*x._odata[ss+s]()); + tmp = tmp + b*y._odata[ss+sp]; + vstream(z._odata[ss+s],tmp); + } +} + +template +void axpbg5y_ssp(Lattice &z,RealD a,const Lattice &x,RealD b,const Lattice &y,int s,int sp) +{ + z.checkerboard = x.checkerboard; + conformable(x,y); + conformable(x,z); + GridBase *grid=x._grid; + int Ls = grid->_rdimensions[0]; +PARALLEL_FOR_LOOP + for(int ss=0;ssoSites();ss+=Ls){ // adds Ls + vobj tmp; + multGamma5(tmp(),b*y._odata[ss+sp]()); + tmp = tmp + a*x._odata[ss+s]; + vstream(z._odata[ss+s],tmp); + } +} + +template +void ag5xpbg5y_ssp(Lattice &z,RealD a,const Lattice &x,RealD b,const Lattice &y,int s,int sp) +{ + z.checkerboard = x.checkerboard; + conformable(x,y); + conformable(x,z); + GridBase *grid=x._grid; + int Ls = grid->_rdimensions[0]; +PARALLEL_FOR_LOOP + for(int ss=0;ssoSites();ss+=Ls){ // adds Ls + vobj tmp1; + vobj tmp2; + tmp1 = a*x._odata[ss+s]+b*y._odata[ss+sp]; + multGamma5(tmp2(),tmp1()); + vstream(z._odata[ss+s],tmp2); + } +} + +template +void axpby_ssp_pminus(Lattice &z,RealD a,const Lattice &x,RealD b,const Lattice &y,int s,int sp) +{ + z.checkerboard = x.checkerboard; + conformable(x,y); + conformable(x,z); + GridBase *grid=x._grid; + int Ls = grid->_rdimensions[0]; +PARALLEL_FOR_LOOP + for(int ss=0;ssoSites();ss+=Ls){ // adds Ls + vobj tmp; + spProj5m(tmp,y._odata[ss+sp]); + tmp = a*x._odata[ss+s]+b*tmp; + vstream(z._odata[ss+s],tmp); + } +} + +template +void axpby_ssp_pplus(Lattice &z,RealD a,const Lattice &x,RealD b,const Lattice &y,int s,int sp) +{ + z.checkerboard = x.checkerboard; + conformable(x,y); + conformable(x,z); + GridBase *grid=x._grid; + int Ls = grid->_rdimensions[0]; +PARALLEL_FOR_LOOP + for(int ss=0;ssoSites();ss+=Ls){ // adds Ls + vobj tmp; + spProj5p(tmp,y._odata[ss+sp]); + tmp = a*x._odata[ss+s]+b*tmp; + vstream(z._odata[ss+s],tmp); + } +} +}} +#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..4b5edc5c 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,10 @@ namespace QCD { } //namespace QCD } // Grid -#include -#include -//#include -#include +#include +#include +#include +#include +#include #endif diff --git a/lib/qcd/SpaceTimeGrid.cc b/lib/qcd/SpaceTimeGrid.cc new file mode 100644 index 00000000..284c5771 --- /dev/null +++ b/lib/qcd/SpaceTimeGrid.cc @@ -0,0 +1,52 @@ +#include + +namespace Grid { + namespace QCD { + +///////////////////////////////////////////////////////////////// +// Public interface +///////////////////////////////////////////////////////////////// +GridCartesian *SpaceTimeGrid::makeFourDimGrid(const std::vector & latt,const std::vector &simd,const std::vector &mpi) +{ + return new GridCartesian(latt,simd,mpi); +} +GridRedBlackCartesian *SpaceTimeGrid::makeFourDimRedBlackGrid(const GridCartesian *FourDimGrid) +{ + return new GridRedBlackCartesian(FourDimGrid); +} + +GridCartesian *SpaceTimeGrid::makeFiveDimGrid(int Ls,const GridCartesian *FourDimGrid) +{ + int N4=FourDimGrid->_ndimension; + + std::vector latt5(1,Ls); + std::vector simd5(1,1); + std::vector mpi5(1,1); + + for(int d=0;d_fdimensions[d]); + simd5.push_back(FourDimGrid->_simd_layout[d]); + mpi5.push_back(FourDimGrid->_processors[d]); + } + return new GridCartesian(latt5,simd5,mpi5); +} + +GridRedBlackCartesian *SpaceTimeGrid::makeFiveDimRedBlackGrid(int Ls,const GridCartesian *FourDimGrid) +{ + int N4=FourDimGrid->_ndimension; + int cbd=1; + std::vector latt5(1,Ls); + std::vector simd5(1,1); + std::vector mpi5(1,1); + std::vector cb5(1,0); + + for(int d=0;d_fdimensions[d]); + simd5.push_back(FourDimGrid->_simd_layout[d]); + mpi5.push_back(FourDimGrid->_processors[d]); + cb5.push_back( 1); + } + return new GridRedBlackCartesian(latt5,simd5,mpi5,cb5,cbd); +} + +}} diff --git a/lib/qcd/SpaceTimeGrid.h b/lib/qcd/SpaceTimeGrid.h new file mode 100644 index 00000000..0b386a0e --- /dev/null +++ b/lib/qcd/SpaceTimeGrid.h @@ -0,0 +1,18 @@ +#ifndef GRID_QCD_SPACE_TIME_GRID_H +#define GRID_QCD_SPACE_TIME_GRID_H +namespace Grid { +namespace QCD { + +class SpaceTimeGrid { + public: + + static GridCartesian *makeFourDimGrid(const std::vector & latt,const std::vector &simd,const std::vector &mpi); + static GridRedBlackCartesian *makeFourDimRedBlackGrid (const GridCartesian *FourDimGrid); + static GridCartesian *makeFiveDimGrid (int Ls,const GridCartesian *FourDimGrid); + static GridRedBlackCartesian *makeFiveDimRedBlackGrid(int Ls,const GridCartesian *FourDimGrid); + +}; + +}} + +#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..b31e9136 --- /dev/null +++ b/lib/qcd/action/Actions.h @@ -0,0 +1,98 @@ +#ifndef GRID_QCD_ACTIONS_H +#define GRID_QCD_ACTIONS_H + + +// Some reorganisation likely required as both Chroma and IroIro +// are separating the concept of the operator from that of action. +// +// The FermAction contains methods to create +// +// * Linear operators (Hermitian and non-hermitian) .. my LinearOperator +// * System solvers (Hermitian and non-hermitian) .. my OperatorFunction +// * MultiShift System solvers (Hermitian and non-hermitian) .. my OperatorFunction + + +//////////////////////////////////////////// +// Abstract base interface +//////////////////////////////////////////// +#include + +//////////////////////////////////////////// +// Utility functions +//////////////////////////////////////////// +#include //used by all wilson type fermions +#include //used by all wilson type fermions + +//////////////////////////////////////////// +// 4D formulations +//////////////////////////////////////////// +#include +//#include + +//////////////////////////////////////////// +// 5D formulations... +//////////////////////////////////////////// + +#include // used by all 5d overlap types + +////////// +// Cayley +////////// +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include + +////////////////////// +// Continued fraction +////////////////////// +#include +#include +#include + +////////////////////// +// Partial fraction +////////////////////// +#include + + + // Chroma interface defining FermionAction + /* + template class FermAct4D : public FermionAction + virtual LinearOperator* linOp(Handle< FermState > state) const = 0; + virtual LinearOperator* lMdagM(Handle< FermState > state) const = 0; + virtual LinOpSystemSolver* invLinOp(Handle< FermState > state, + virtual MdagMSystemSolver* invMdagM(Handle< FermState > state, + virtual LinOpMultiSystemSolver* mInvLinOp(Handle< FermState > state, + virtual MdagMMultiSystemSolver* mInvMdagM(Handle< FermState > state, + virtual MdagMMultiSystemSolverAccumulate* mInvMdagMAcc(Handle< FermState > state, + virtual SystemSolver* qprop(Handle< FermState > state, + class DiffFermAct4D : public FermAct4D + virtual DiffLinearOperator* linOp(Handle< FermState > state) const = 0; + virtual DiffLinearOperator* lMdagM(Handle< FermState > state) const = 0; + */ + + + // Chroma interface defining GaugeAction + /* + template class GaugeAction + virtual const CreateGaugeState& getCreateState() const = 0; + virtual GaugeState* createState(const Q& q) const + virtual const GaugeBC& getGaugeBC() const + virtual const Set& getSet(void) const = 0; + virtual void deriv(P& result, const Handle< GaugeState >& state) const + virtual Double S(const Handle< GaugeState >& state) const = 0; + + class LinearGaugeAction : public GaugeAction< multi1d, multi1d > + typedef multi1d P; + */ + +#endif diff --git a/lib/qcd/action/fermion/.dirstamp b/lib/qcd/action/fermion/.dirstamp new file mode 100644 index 00000000..e69de29b diff --git a/lib/qcd/action/fermion/CayleyFermion5D.cc b/lib/qcd/action/fermion/CayleyFermion5D.cc new file mode 100644 index 00000000..e47ff331 --- /dev/null +++ b/lib/qcd/action/fermion/CayleyFermion5D.cc @@ -0,0 +1,346 @@ +#include +namespace Grid { +namespace QCD { + + CayleyFermion5D::CayleyFermion5D(LatticeGaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD _M5) : + WilsonFermion5D(_Umu, + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_M5), + mass(_mass) + { + } + + // override multiply + RealD CayleyFermion5D::M (const LatticeFermion &psi, LatticeFermion &chi) + { + LatticeFermion Din(psi._grid); + + // Assemble Din + for(int s=0;s D1+^dag P+ D2-^dag + //D2- P+ D2+ P-D1-^dag D2+dag + + LatticeFermion Din(psi._grid); + // Apply Dw + DW(psi,Din,DaggerYes); + + for(int s=0;s=0;s--){ + axpby_ssp_pminus (chi,1.0,chi,-uee[s],chi,s,s+1); // chi[Ls] + } + } + + void CayleyFermion5D::MooeeInvDag (const LatticeFermion &psi, LatticeFermion &chi) + { + // Apply (U^{\prime})^{-dagger} + axpby_ssp (chi,1.0,psi, 0.0,psi,0,0); // chi[0]=psi[0] + for (int s=1;s=0;s--){ + axpby_ssp_pplus (chi,1.0,chi,-lee[s],chi,s,s+1); // chi[Ls] + } + } + + // Tanh + void CayleyFermion5D::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c) + { + SetCoefficientsZolotarev(1.0,zdata,b,c); + + } + //Zolo + void CayleyFermion5D::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata,RealD b,RealD c) + { + + /////////////////////////////////////////////////////////// + // The Cayley coeffs (unprec) + /////////////////////////////////////////////////////////// + omega.resize(Ls); + bs.resize(Ls); + cs.resize(Ls); + as.resize(Ls); + + // + // Ts = ( [bs+cs]Dw )^-1 ( (bs+cs) Dw ) + // -(g5 ------- -1 ) ( g5 --------- + 1 ) + // ( {2+(bs-cs)Dw} ) ( 2+(bs-cs) Dw ) + // + // bs = 1/2( (1/omega_s + 1)*b + (1/omega - 1)*c ) = 1/2( 1/omega(b+c) + (b-c) ) + // cs = 1/2( (1/omega_s - 1)*b + (1/omega + 1)*c ) = 1/2( 1/omega(b+c) - (b-c) ) + // + // bs+cs = 0.5*( 1/omega(b+c) + (b-c) + 1/omega(b+c) - (b-c) ) = 1/omega(b+c) + // bs-cs = 0.5*( 1/omega(b+c) + (b-c) - 1/omega(b+c) + (b-c) ) = b-c + // + // So + // + // Ts = ( [b+c]Dw/omega_s )^-1 ( (b+c) Dw /omega_s ) + // -(g5 ------- -1 ) ( g5 --------- + 1 ) + // ( {2+(b-c)Dw} ) ( 2+(b-c) Dw ) + // + // Ts = ( [b+c]Dw )^-1 ( (b+c) Dw ) + // -(g5 ------- -omega_s) ( g5 --------- + omega_s ) + // ( {2+(b-c)Dw} ) ( 2+(b-c) Dw ) + // + + double bpc = b+c; + double bmc = b-c; + for(int i=0; i < Ls; i++){ + as[i] = 1.0; + omega[i] = ((double)zdata->gamma[i])*zolo_hi; //NB reciprocal relative to Chroma NEF code + bs[i] = 0.5*(bpc/omega[i] + bmc); + cs[i] = 0.5*(bpc/omega[i] - bmc); + } + + //////////////////////////////////////////////////////// + // Constants for the preconditioned matrix Cayley form + //////////////////////////////////////////////////////// + bee.resize(Ls); + cee.resize(Ls); + beo.resize(Ls); + ceo.resize(Ls); + + for(int i=0;i omega; + std::vector bs; // S dependent coeffs + std::vector cs; + std::vector as; + // For preconditioning Cayley form + std::vector bee; + std::vector cee; + std::vector aee; + std::vector beo; + std::vector ceo; + std::vector aeo; + // LDU factorisation of the eeoo matrix + std::vector lee; + std::vector leem; + std::vector uee; + std::vector ueem; + std::vector dee; + + // Constructors + CayleyFermion5D(LatticeGaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD _M5); + + protected: + void SetCoefficientsZolotarev(RealD zolohi,Approx::zolotarev_data *zdata,RealD b,RealD c); + void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c); + }; + + } +} + +#endif diff --git a/lib/qcd/action/fermion/ContinuedFractionFermion5D.cc b/lib/qcd/action/fermion/ContinuedFractionFermion5D.cc new file mode 100644 index 00000000..92f6473e --- /dev/null +++ b/lib/qcd/action/fermion/ContinuedFractionFermion5D.cc @@ -0,0 +1,171 @@ +#include + +namespace Grid { + namespace QCD { + + void ContinuedFractionFermion5D::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD scale) + { + SetCoefficientsZolotarev(1.0/scale,zdata); + } + void ContinuedFractionFermion5D::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata) + { + R=(1+this->mass)/(1-this->mass); + + Beta.resize(Ls); + cc.resize(Ls); + cc_d.resize(Ls); + sqrt_cc.resize(Ls); + for(int i=0; i < Ls ; i++){ + Beta[i] = zdata -> beta[i]; + cc[i] = 1.0/Beta[i]; + cc_d[i]=sqrt(cc[i]); + } + + cc_d[Ls-1]=1.0; + for(int i=0; i < Ls-1 ; i++){ + sqrt_cc[i]= sqrt(cc[i]*cc[i+1]); + } + sqrt_cc[Ls-2]=sqrt(cc[Ls-2]); + + + ZoloHiInv =1.0/zolo_hi; + double dw_diag = (4.0-M5)*ZoloHiInv; + + See.resize(Ls); + Aee.resize(Ls); + int sign=1; + for(int s=0;s=0;s--){ + axpbg5y_ssp(chi,1.0/cc_d[s],chi,-1.0*cc_d[s+1]/See[s]/cc_d[s],chi,s,s+1); + } + } + void ContinuedFractionFermion5D::MooeeInvDag (const LatticeFermion &psi, LatticeFermion &chi) + { + MooeeInv(psi,chi); + } + + // Constructors + ContinuedFractionFermion5D::ContinuedFractionFermion5D( + LatticeGaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD M5) : + WilsonFermion5D(_Umu, + FiveDimGrid, FiveDimRedBlackGrid, + FourDimGrid, FourDimRedBlackGrid,M5), + mass(_mass) + { + assert((Ls&0x1)==1); // Odd Ls required + } + + } +} + diff --git a/lib/qcd/action/fermion/ContinuedFractionFermion5D.h b/lib/qcd/action/fermion/ContinuedFractionFermion5D.h new file mode 100644 index 00000000..f363878f --- /dev/null +++ b/lib/qcd/action/fermion/ContinuedFractionFermion5D.h @@ -0,0 +1,59 @@ +#ifndef GRID_QCD_CONTINUED_FRACTION_H +#define GRID_QCD_CONTINUED_FRACTION_H + +namespace Grid { + + namespace QCD { + + class ContinuedFractionFermion5D : public WilsonFermion5D + { + public: + + // 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); + + // virtual void Instantiatable(void)=0; + virtual void Instantiatable(void) =0; + + // Constructors + ContinuedFractionFermion5D(LatticeGaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD M5); + + protected: + + void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD scale); + void SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata);; + + Approx::zolotarev_data *zdata; + + // Cont frac + RealD mass; + RealD R; + RealD ZoloHiInv; + std::vector Beta; + std::vector cc;; + std::vector cc_d;; + std::vector sqrt_cc; + std::vector See; + std::vector Aee; + + }; + + + } +} + +#endif diff --git a/lib/qcd/action/fermion/DomainWallFermion.h b/lib/qcd/action/fermion/DomainWallFermion.h new file mode 100644 index 00000000..a25c0c3c --- /dev/null +++ b/lib/qcd/action/fermion/DomainWallFermion.h @@ -0,0 +1,46 @@ +#ifndef GRID_QCD_DOMAIN_WALL_FERMION_H +#define GRID_QCD_DOMAIN_WALL_FERMION_H + +#include + +namespace Grid { + + namespace QCD { + + class DomainWallFermion : public CayleyFermion5D + { + public: + + virtual void Instantiatable(void) {}; + // Constructors + DomainWallFermion(LatticeGaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD _M5) : + + CayleyFermion5D(_Umu, + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_mass,_M5) + + { + RealD eps = 1.0; + + Approx::zolotarev_data *zdata = Approx::grid_higham(eps,this->Ls);// eps is ignored for higham + assert(zdata->n==this->Ls); + + std::cout << "DomainWallFermion with Ls="<CayleyFermion5D::SetCoefficientsTanh(zdata,1.0,0.0); + + } + + }; + + } +} + +#endif diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h new file mode 100644 index 00000000..47c47478 --- /dev/null +++ b/lib/qcd/action/fermion/FermionOperator.h @@ -0,0 +1,48 @@ +#ifndef GRID_QCD_FERMION_OPERATOR_H +#define GRID_QCD_FERMION_OPERATOR_H + +namespace Grid { + + namespace QCD { + + ////////////////////////////////////////////////////////////////////////////// + // Four component fermions + // Should type template the vector and gauge types + // Think about multiple representations + ////////////////////////////////////////////////////////////////////////////// + template + class FermionOperator : 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/MobiusFermion.h b/lib/qcd/action/fermion/MobiusFermion.h new file mode 100644 index 00000000..33f94089 --- /dev/null +++ b/lib/qcd/action/fermion/MobiusFermion.h @@ -0,0 +1,47 @@ +#ifndef GRID_QCD_MOBIUS_FERMION_H +#define GRID_QCD_MOBIUS_FERMION_H + +#include + +namespace Grid { + + namespace QCD { + + class MobiusFermion : public CayleyFermion5D + { + public: + + virtual void Instantiatable(void) {}; + // Constructors + MobiusFermion(LatticeGaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD _M5, + RealD b, RealD c) : + + CayleyFermion5D(_Umu, + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_mass,_M5) + + { + RealD eps = 1.0; + + std::cout << "MobiusFermion (b="<Ls);// eps is ignored for higham + assert(zdata->n==this->Ls); + + // Call base setter + this->CayleyFermion5D::SetCoefficientsTanh(zdata,b,c); + + } + + }; + + } +} + +#endif diff --git a/lib/qcd/action/fermion/MobiusZolotarevFermion.h b/lib/qcd/action/fermion/MobiusZolotarevFermion.h new file mode 100644 index 00000000..1be61601 --- /dev/null +++ b/lib/qcd/action/fermion/MobiusZolotarevFermion.h @@ -0,0 +1,49 @@ +#ifndef GRID_QCD_MOBIUS_ZOLOTAREV_FERMION_H +#define GRID_QCD_MOBIUS_ZOLOTAREV_FERMION_H + +#include + +namespace Grid { + + namespace QCD { + + class MobiusZolotarevFermion : public CayleyFermion5D + { + public: + + virtual void Instantiatable(void) {}; + // Constructors + MobiusZolotarevFermion(LatticeGaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD _M5, + RealD b, RealD c, + RealD lo, RealD hi) : + + CayleyFermion5D(_Umu, + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_mass,_M5) + + { + RealD eps = lo/hi; + + Approx::zolotarev_data *zdata = Approx::grid_zolotarev(eps,this->Ls,0);// eps is ignored for higham + assert(zdata->n==this->Ls); + + std::cout << "MobiusZolotarevFermion (b="<CayleyFermion5D::SetCoefficientsZolotarev(hi,zdata,b,c); + + } + + }; + + } +} + +#endif diff --git a/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h b/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h new file mode 100644 index 00000000..e764c8ae --- /dev/null +++ b/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h @@ -0,0 +1,34 @@ +#ifndef OVERLAP_WILSON_CAYLEY_TANH_FERMION_H +#define OVERLAP_WILSON_CAYLEY_TANH_FERMION_H + +#include + +namespace Grid { + + namespace QCD { + + class OverlapWilsonCayleyTanhFermion : public MobiusFermion + { + public: + + // Constructors + OverlapWilsonCayleyTanhFermion(LatticeGaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD _M5, + RealD scale) : + + // b+c=scale, b-c = 0 <=> b =c = scale/2 + MobiusFermion(_Umu, + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_mass,_M5,0.5*scale,0.5*scale) + { + } + }; + } +} +#endif diff --git a/lib/qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h b/lib/qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h new file mode 100644 index 00000000..82c43fb7 --- /dev/null +++ b/lib/qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h @@ -0,0 +1,37 @@ +#ifndef OVERLAP_WILSON_CAYLEY_ZOLOTAREV_FERMION_H +#define OVERLAP_WILSON_CAYLEY_ZOLOTAREV_FERMION_H + +#include + +namespace Grid { + + namespace QCD { + + class OverlapWilsonCayleyZolotarevFermion : public MobiusZolotarevFermion + { + public: + + // Constructors + + OverlapWilsonCayleyZolotarevFermion(LatticeGaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD _M5, + RealD lo, RealD hi) : + // b+c=1.0, b-c = 0 <=> b =c = 1/2 + MobiusZolotarevFermion(_Umu, + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_mass,_M5,0.5,0.5,lo,hi) + + {} + + }; + + } +} + +#endif diff --git a/lib/qcd/action/fermion/OverlapWilsonContfracTanhFermion.h b/lib/qcd/action/fermion/OverlapWilsonContfracTanhFermion.h new file mode 100644 index 00000000..ed0c24dc --- /dev/null +++ b/lib/qcd/action/fermion/OverlapWilsonContfracTanhFermion.h @@ -0,0 +1,39 @@ +#ifndef OVERLAP_WILSON_CONTFRAC_TANH_FERMION_H +#define OVERLAP_WILSON_CONTFRAC_TANH_FERMION_H + +#include + +namespace Grid { + + namespace QCD { + + class OverlapWilsonContFracTanhFermion : public ContinuedFractionFermion5D + { + public: + + virtual void Instantiatable(void){}; + // Constructors + OverlapWilsonContFracTanhFermion(LatticeGaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD _M5, + RealD scale) : + + // b+c=scale, b-c = 0 <=> b =c = scale/2 + ContinuedFractionFermion5D(_Umu, + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_mass,_M5) + { + assert((Ls&0x1)==1); // Odd Ls required + int nrational=Ls-1;// Even rational order + zdata = Approx::grid_higham(1.0,nrational);// eps is ignored for higham + SetCoefficientsTanh(zdata,scale); + } + }; + } +} +#endif diff --git a/lib/qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h b/lib/qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h new file mode 100644 index 00000000..caf01133 --- /dev/null +++ b/lib/qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h @@ -0,0 +1,44 @@ +#ifndef OVERLAP_WILSON_CONTFRAC_ZOLOTAREV_FERMION_H +#define OVERLAP_WILSON_CONTFRAC_ZOLOTAREV_FERMION_H + +#include + +namespace Grid { + + namespace QCD { + + class OverlapWilsonContFracZolotarevFermion : public ContinuedFractionFermion5D + { + public: + + virtual void Instantiatable(void){}; + // Constructors + OverlapWilsonContFracZolotarevFermion(LatticeGaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD _M5, + RealD lo,RealD hi): + + // b+c=scale, b-c = 0 <=> b =c = scale/2 + ContinuedFractionFermion5D(_Umu, + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_mass,_M5) + { + assert((Ls&0x1)==1); // Odd Ls required + + int nrational=Ls-1;// Even rational order + RealD eps = lo/hi; + + Approx::zolotarev_data *zdata = Approx::grid_zolotarev(eps,nrational,0); + + SetCoefficientsZolotarev(hi,zdata); + + } + }; + } +} +#endif diff --git a/lib/qcd/action/fermion/PartialFractionFermion5D.cc b/lib/qcd/action/fermion/PartialFractionFermion5D.cc new file mode 100644 index 00000000..8b137891 --- /dev/null +++ b/lib/qcd/action/fermion/PartialFractionFermion5D.cc @@ -0,0 +1 @@ + diff --git a/lib/qcd/action/fermion/PartialFractionFermion5D.h b/lib/qcd/action/fermion/PartialFractionFermion5D.h new file mode 100644 index 00000000..95f8c0f9 --- /dev/null +++ b/lib/qcd/action/fermion/PartialFractionFermion5D.h @@ -0,0 +1,49 @@ +#ifndef GRID_QCD_PARTIAL_FRACTION_H +#define GRID_QCD_PARTIAL_FRACTION_H + +namespace Grid { + + namespace QCD { + + class PartialFractionFermion5D : public WilsonFermion5D + { + public: + + // 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); + + private: + + virtual void PartialFractionCoefficients(void); + + Approx::zolotarev_data *zdata; + + // Part frac + double R; + std::vector p; + std::vector q; + + // Constructors + PartialFractionFermion5D(LatticeGaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD M5); + + }; + + + } +} + +#endif diff --git a/lib/qcd/action/fermion/ScaledShamirFermion.h b/lib/qcd/action/fermion/ScaledShamirFermion.h new file mode 100644 index 00000000..59fb16a8 --- /dev/null +++ b/lib/qcd/action/fermion/ScaledShamirFermion.h @@ -0,0 +1,37 @@ +#ifndef GRID_QCD_SCALED_SHAMIR_FERMION_H +#define GRID_QCD_SCALED_SHAMIR_FERMION_H + +#include + +namespace Grid { + + namespace QCD { + + class ScaledShamirFermion : public MobiusFermion + { + public: + + // Constructors + ScaledShamirFermion(LatticeGaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD _M5, + RealD scale) : + + // b+c=scale, b-c = 1 <=> 2b = scale+1; 2c = scale-1 + MobiusFermion(_Umu, + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_mass,_M5,0.5*(scale+1.0),0.5*(scale-1.0)) + { + } + + }; + + } +} + +#endif diff --git a/lib/qcd/action/fermion/ShamirZolotarevFermion.h b/lib/qcd/action/fermion/ShamirZolotarevFermion.h new file mode 100644 index 00000000..6a7df439 --- /dev/null +++ b/lib/qcd/action/fermion/ShamirZolotarevFermion.h @@ -0,0 +1,39 @@ +#ifndef GRID_QCD_SHAMIR_ZOLOTAREV_FERMION_H +#define GRID_QCD_SHAMIR_ZOLOTAREV_FERMION_H + +#include + +namespace Grid { + + namespace QCD { + + class ShamirZolotarevFermion : public MobiusZolotarevFermion + { + public: + + // Constructors + + + ShamirZolotarevFermion(LatticeGaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD _M5, + RealD lo, RealD hi) : + + // b+c = 1; b-c = 1 => b=1, c=0 + MobiusZolotarevFermion(_Umu, + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_mass,_M5,1.0,0.0,lo,hi) + + {} + + }; + + } +} + +#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..9f2da251 --- /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, + RealD _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..6040d328 --- /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 FermionOperator + { + 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,RealD _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: + + RealD 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/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc new file mode 100644 index 00000000..d22701b0 --- /dev/null +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -0,0 +1,186 @@ +#include + +namespace Grid { +namespace QCD { + + // S-direction is INNERMOST and takes no part in the parity. + const std::vector WilsonFermion5D::directions ({1,2,3,4, 1, 2, 3, 4}); + const std::vector WilsonFermion5D::displacements({1,1,1,1,-1,-1,-1,-1}); + + int WilsonFermion5D::HandOptDslash; + + // 5d lattice for DWF. + WilsonFermion5D::WilsonFermion5D(LatticeGaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _M5) : + _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 + M5(_M5), + 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 WilsonFermion5D::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); + } +} +void WilsonFermion5D::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 ) { +PARALLEL_FOR_LOOP + for(int ss=0;ssoSites();ss++){ + for(int s=0;soSites();ss++){ + for(int s=0;soSites();ss++){ + for(int s=0;soSites();ss++){ + 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;} + + // full checkerboard operations; leave unimplemented as abstract for now + //virtual RealD M (const LatticeFermion &in, LatticeFermion &out)=0; + //virtual RealD Mdag (const LatticeFermion &in, LatticeFermion &out)=0; + + // half checkerboard operations; leave unimplemented as abstract for now + // virtual void Meooe (const LatticeFermion &in, LatticeFermion &out)=0; + // virtual void MeooeDag (const LatticeFermion &in, LatticeFermion &out)=0; + // virtual void Mooee (const LatticeFermion &in, LatticeFermion &out)=0; + // virtual void MooeeDag (const LatticeFermion &in, LatticeFermion &out)=0; + // virtual void MooeeInv (const LatticeFermion &in, LatticeFermion &out)=0; + // virtual void MooeeInvDag (const LatticeFermion &in, LatticeFermion &out)=0; + + // Implement hopping term non-hermitian hopping term; half cb or both + // Implement s-diagonal DW + void DW (const LatticeFermion &in, LatticeFermion &out,int dag); + 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 + WilsonFermion5D(LatticeGaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + double _M5); + + // 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; + 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/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/stamp-h1 b/lib/stamp-h1 index 753f890f..2a8e9ffe 100644 --- a/lib/stamp-h1 +++ b/lib/stamp-h1 @@ -1 +1 @@ -timestamp for lib/Grid_config.h +timestamp for lib/GridConfig.h diff --git a/lib/stencil/Lebesgue.cc b/lib/stencil/Lebesgue.cc new file mode 100644 index 00000000..977e3562 --- /dev/null +++ b/lib/stencil/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/Lebesgue.h b/lib/stencil/Lebesgue.h new file mode 100644 index 00000000..1d59b127 --- /dev/null +++ b/lib/stencil/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/Stencil_common.cc similarity index 73% rename from lib/stencil/Grid_stencil_common.cc rename to lib/stencil/Stencil_common.cc index 6cbcc890..7f894faf 100644 --- a/lib/stencil/Grid_stencil_common.cc +++ b/lib/stencil/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->CheckerBoardShiftForCB(_checkerboard,dimension,shift,Even); + sshift[1] = _grid->CheckerBoardShiftForCB(_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->CheckerBoardShiftForCB(_checkerboard,dimension,shift,Even); + sshift[1] = _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,Odd); if ( sshift[0] == sshift[1] ) { Comms(point,dimension,shift,0x3); } else { @@ -185,9 +94,9 @@ 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 sshift = _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,cb); int sx = (x+sshift)%rd; int permute_slice=0; @@ -224,8 +133,8 @@ 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 sshift= _grid->CheckerBoardShift(_checkerboard,dimension,shift,cb); + int cb= (cbmask==0x2)? Odd : Even; + int sshift= _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,cb); for(int x=0;x inline iScalar Ta(const iScalar&r) + { + iScalar ret; + ret._internal = Ta(r._internal); + return ret; + } + template inline iVector Ta(const iVector&r) + { + iVector ret; + for(int i=0;i inline iMatrix Ta(const iMatrix &arg) + { + iMatrix ret(arg); + vtype factor = (1/(double)N); + ret = (ret - adj(arg))*0.5; + ret -= trace(ret)*factor; + return ret; + } + +} +#endif diff --git a/lib/tensors/Tensor_arith.h b/lib/tensors/Tensor_arith.h new file mode 100644 index 00000000..853a19a5 --- /dev/null +++ b/lib/tensors/Tensor_arith.h @@ -0,0 +1,11 @@ +#ifndef GRID_MATH_ARITH_H +#define GRID_MATH_ARITH_H + +#include +#include +#include +#include +#include + +#endif + diff --git a/lib/math/Grid_math_arith_add.h b/lib/tensors/Tensor_arith_add.h similarity index 96% rename from lib/math/Grid_math_arith_add.h rename to lib/tensors/Tensor_arith_add.h index 19a2f407..66d5879f 100644 --- a/lib/math/Grid_math_arith_add.h +++ b/lib/tensors/Tensor_arith_add.h @@ -45,7 +45,10 @@ namespace Grid { { for(int c2=0;c2_internal[c1][c2],&lhs->_internal,&rhs->_internal[c1][c2]); + if ( c1==c2) + add(&ret->_internal[c1][c2],&lhs->_internal,&rhs->_internal[c1][c2]); + else + ret->_internal[c1][c2]=lhs->_internal[c1][c2]; }} return; } diff --git a/lib/math/Grid_math_arith_mac.h b/lib/tensors/Tensor_arith_mac.h similarity index 100% rename from lib/math/Grid_math_arith_mac.h rename to lib/tensors/Tensor_arith_mac.h diff --git a/lib/math/Grid_math_arith_mul.h b/lib/tensors/Tensor_arith_mul.h similarity index 100% rename from lib/math/Grid_math_arith_mul.h rename to lib/tensors/Tensor_arith_mul.h diff --git a/lib/math/Grid_math_arith_scalar.h b/lib/tensors/Tensor_arith_scalar.h similarity index 100% rename from lib/math/Grid_math_arith_scalar.h rename to lib/tensors/Tensor_arith_scalar.h diff --git a/lib/math/Grid_math_arith_sub.h b/lib/tensors/Tensor_arith_sub.h similarity index 99% rename from lib/math/Grid_math_arith_sub.h rename to lib/tensors/Tensor_arith_sub.h index c359349f..b7a5a54d 100644 --- a/lib/math/Grid_math_arith_sub.h +++ b/lib/tensors/Tensor_arith_sub.h @@ -44,7 +44,7 @@ template strong_inline void sub(iMat const iMatrix * __restrict__ rhs){ for(int c2=0;c2_internal[c1][c2],&lhs->_internal,&rhs->_internal[c1][c2]); } else { // Fails -- need unary minus. Catalogue other unops? @@ -60,7 +60,7 @@ template strong_inline void sub(iMat const iScalar * __restrict__ rhs){ for(int c2=0;c2_internal[c1][c2],&lhs->_internal[c1][c2],&rhs->_internal); else ret->_internal[c1][c2]=lhs->_internal[c1][c2]; diff --git a/lib/math/Grid_math_tensors.h b/lib/tensors/Tensor_class.h similarity index 100% rename from lib/math/Grid_math_tensors.h rename to lib/tensors/Tensor_class.h diff --git a/lib/Grid_extract.h b/lib/tensors/Tensor_extract_merge.h similarity index 100% rename from lib/Grid_extract.h rename to lib/tensors/Tensor_extract_merge.h diff --git a/lib/math/Grid_math_inner.h b/lib/tensors/Tensor_inner.h similarity index 100% rename from lib/math/Grid_math_inner.h rename to lib/tensors/Tensor_inner.h diff --git a/lib/math/Grid_math_outer.h b/lib/tensors/Tensor_outer.h similarity index 100% rename from lib/math/Grid_math_outer.h rename to lib/tensors/Tensor_outer.h diff --git a/lib/math/Grid_math_peek.h b/lib/tensors/Tensor_peek.h similarity index 100% rename from lib/math/Grid_math_peek.h rename to lib/tensors/Tensor_peek.h diff --git a/lib/math/Grid_math_poke.h b/lib/tensors/Tensor_poke.h similarity index 100% rename from lib/math/Grid_math_poke.h rename to lib/tensors/Tensor_poke.h diff --git a/lib/math/Grid_math_reality.h b/lib/tensors/Tensor_reality.h similarity index 94% rename from lib/math/Grid_math_reality.h rename to lib/tensors/Tensor_reality.h index 40d4c624..c35a8630 100644 --- a/lib/math/Grid_math_reality.h +++ b/lib/tensors/Tensor_reality.h @@ -2,10 +2,6 @@ #define GRID_MATH_REALITY_H namespace Grid { - /////////////////////////////////////////////////////////////////////////////////////////////////// - /////////////////////////////////////////// CONJ /////////////////////////////////////////// - /////////////////////////////////////////////////////////////////////////////////////////////////// - /////////////////////////////////////////////// // multiply by I; make recursive. /////////////////////////////////////////////// @@ -151,6 +147,9 @@ template inline iMatrix adj(const iMatrix & + + + ///////////////////////////////////////////////////////////////// // Can only take the real/imag part of scalar objects, since // lattice objects of different complex nature are non-conformable. diff --git a/lib/math/Grid_math_trace.h b/lib/tensors/Tensor_trace.h similarity index 100% rename from lib/math/Grid_math_trace.h rename to lib/tensors/Tensor_trace.h diff --git a/lib/math/Grid_math_traits.h b/lib/tensors/Tensor_traits.h similarity index 100% rename from lib/math/Grid_math_traits.h rename to lib/tensors/Tensor_traits.h diff --git a/lib/math/Grid_math_transpose.h b/lib/tensors/Tensor_transpose.h similarity index 100% rename from lib/math/Grid_math_transpose.h rename to lib/tensors/Tensor_transpose.h diff --git a/scripts/filelist b/scripts/filelist new file mode 100755 index 00000000..99db51b2 --- /dev/null +++ b/scripts/filelist @@ -0,0 +1,58 @@ +#!/bin/bash + +cd lib + +HFILES=`find . -type f -name '*.h'` +CCFILES=`find . -type f -name '*.cc' -not -name '*ommunicator*.cc'` +echo> Make.inc +echo HFILES=$HFILES >> Make.inc +echo >> Make.inc +echo CCFILES=$CCFILES >> Make.inc + +cd .. + + + +cd tests + +echo> Make.inc +TESTS=`ls T*.cc` +TESTLIST=`echo ${TESTS} | sed s/.cc//g ` + +echo > Make.inc +echo bin_PROGRAMS = ${TESTLIST} >> Make.inc +echo >> Make.inc + +for f in $TESTS +do +BNAME=`basename $f .cc` +echo >> Make.inc +echo ${BNAME}_SOURCES=$f >> Make.inc +echo ${BNAME}_LDADD=-lGrid>> Make.inc +echo >> Make.inc +done + +cd .. + + +cd benchmarks + + +echo> Make.inc +TESTS=`ls B*.cc` +TESTLIST=`echo ${TESTS} | sed s/.cc//g ` + +echo > Make.inc +echo bin_PROGRAMS = ${TESTLIST} >> Make.inc +echo >> Make.inc + +for f in $TESTS +do +BNAME=`basename $f .cc` +echo >> Make.inc +echo ${BNAME}_SOURCES=$f >> Make.inc +echo ${BNAME}_LDADD=-lGrid>> Make.inc +echo >> Make.inc +done + +cd .. diff --git a/scripts/linecount b/scripts/linecount index c0731b23..7624c15f 100755 --- a/scripts/linecount +++ b/scripts/linecount @@ -1,3 +1,3 @@ #!/bin/sh -wc -l */*.h */*/*.h */*/*/*.h */*.cc */*/*.cc */*/*/*.cc \ No newline at end of file +wc -l lib/*.h lib/*/*.h lib/*/*/*.h lib/*.cc lib/*/*.cc lib/*/*/*.cc tests/*.cc benchmarks/*.cc lib/*/*/*/*.cc lib/*/*/*/*.h \ No newline at end of file diff --git a/tests/Grid_simd_new.cc b/tests/Grid_simd_new.cc deleted file mode 100644 index 41781304..00000000 --- a/tests/Grid_simd_new.cc +++ /dev/null @@ -1,165 +0,0 @@ -#include -#include "simd/Grid_vector_types.h" -#include - -using namespace std; -using namespace Grid; -using namespace Grid::QCD; - -class funcPlus { -public: - funcPlus() {}; - template void operator()(vec &rr,vec &i1,vec &i2) const { rr = i1+i2;} - std::string name(void) const { return std::string("Plus"); } -}; -class funcMinus { -public: - funcMinus() {}; - template void operator()(vec &rr,vec &i1,vec &i2) const { rr = i1-i2;} - std::string name(void) const { return std::string("Minus"); } -}; -class funcTimes { -public: - funcTimes() {}; - template void operator()(vec &rr,vec &i1,vec &i2) const { rr = i1*i2;} - std::string name(void) const { return std::string("Times"); } -}; -class funcConj { -public: - funcConj() {}; - template void operator()(vec &rr,vec &i1,vec &i2) const { rr = conjugate(i1);} - std::string name(void) const { return std::string("Conj"); } -}; -class funcAdj { -public: - funcAdj() {}; - template void operator()(vec &rr,vec &i1,vec &i2) const { rr = adj(i1);} - std::string name(void) const { return std::string("Adj"); } -}; - -class funcTimesI { -public: - funcTimesI() {}; - template void operator()(vec &rr,vec &i1,vec &i2) const { rr = timesI(i1);} - std::string name(void) const { return std::string("timesI"); } -}; - -class funcTimesMinusI { -public: - funcTimesMinusI() {}; - template void operator()(vec &rr,vec &i1,vec &i2) const { rr = timesMinusI(i1);} - std::string name(void) const { return std::string("timesMinusI"); } -}; - -template -void Tester(const functor &func) -{ - GridSerialRNG sRNG; - sRNG.SeedRandomDevice(); - - int Nsimd = vec::Nsimd(); - - std::vector input1(Nsimd); - std::vector input2(Nsimd); - std::vector result(Nsimd); - std::vector reference(Nsimd); - - std::vector > buf(3); - vec & v_input1 = buf[0]; - vec & v_input2 = buf[1]; - vec & v_result = buf[2]; - - - for(int i=0;i(v_input1,input1); - merge(v_input2,input2); - merge(v_result,result); - - func(v_result,v_input1,v_input2); - - for(int i=0;i(v_result,result); - std::cout << " " << func.name()<0){ - std::cout<< "*****" << std::endl; - std::cout<< "["< latt_size = GridDefaultLatt(); - std::vector simd_layout = GridDefaultSimd(4,MyComplexF::Nsimd()); - std::vector mpi_layout = GridDefaultMpi(); - - GridCartesian Grid(latt_size,simd_layout,mpi_layout); - std::vector seeds({1,2,3,4}); - - // Insist that operations on random scalars gives - // identical results to on vectors. - - std::cout << "==================================="<< std::endl; - std::cout << "Testing MyComplexF "<(funcTimesI()); - Tester(funcTimesMinusI()); - Tester(funcPlus()); - Tester(funcMinus()); - Tester(funcTimes()); - Tester(funcConj()); - Tester(funcAdj()); - - std::cout << "==================================="<< std::endl; - std::cout << "Testing MyComplexD "<(funcTimesI()); - Tester(funcTimesMinusI()); - Tester(funcPlus()); - Tester(funcMinus()); - Tester(funcTimes()); - Tester(funcConj()); - Tester(funcAdj()); - - std::cout << "==================================="<< std::endl; - std::cout << "Testing MyRealF "<(funcPlus()); - Tester(funcMinus()); - Tester(funcTimes()); - Tester(funcAdj()); - - std::cout << "==================================="<< std::endl; - std::cout << "Testing MyRealD "<(funcPlus()); - Tester(funcMinus()); - Tester(funcTimes()); - Tester(funcAdj()); - - Grid_finalize(); -} diff --git a/tests/Make.inc b/tests/Make.inc new file mode 100644 index 00000000..b525874d --- /dev/null +++ b/tests/Make.inc @@ -0,0 +1,91 @@ + +bin_PROGRAMS = Test_cayley_cg Test_cayley_even_odd Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_even_odd Test_gamma Test_main Test_nersc_io Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_even_odd + + +Test_cayley_cg_SOURCES=Test_cayley_cg.cc +Test_cayley_cg_LDADD=-lGrid + + +Test_cayley_even_odd_SOURCES=Test_cayley_even_odd.cc +Test_cayley_even_odd_LDADD=-lGrid + + +Test_contfrac_cg_SOURCES=Test_contfrac_cg.cc +Test_contfrac_cg_LDADD=-lGrid + + +Test_contfrac_even_odd_SOURCES=Test_contfrac_even_odd.cc +Test_contfrac_even_odd_LDADD=-lGrid + + +Test_cshift_SOURCES=Test_cshift.cc +Test_cshift_LDADD=-lGrid + + +Test_cshift_red_black_SOURCES=Test_cshift_red_black.cc +Test_cshift_red_black_LDADD=-lGrid + + +Test_dwf_cg_prec_SOURCES=Test_dwf_cg_prec.cc +Test_dwf_cg_prec_LDADD=-lGrid + + +Test_dwf_cg_schur_SOURCES=Test_dwf_cg_schur.cc +Test_dwf_cg_schur_LDADD=-lGrid + + +Test_dwf_cg_unprec_SOURCES=Test_dwf_cg_unprec.cc +Test_dwf_cg_unprec_LDADD=-lGrid + + +Test_dwf_even_odd_SOURCES=Test_dwf_even_odd.cc +Test_dwf_even_odd_LDADD=-lGrid + + +Test_gamma_SOURCES=Test_gamma.cc +Test_gamma_LDADD=-lGrid + + +Test_main_SOURCES=Test_main.cc +Test_main_LDADD=-lGrid + + +Test_nersc_io_SOURCES=Test_nersc_io.cc +Test_nersc_io_LDADD=-lGrid + + +Test_remez_SOURCES=Test_remez.cc +Test_remez_LDADD=-lGrid + + +Test_rng_SOURCES=Test_rng.cc +Test_rng_LDADD=-lGrid + + +Test_rng_fixed_SOURCES=Test_rng_fixed.cc +Test_rng_fixed_LDADD=-lGrid + + +Test_simd_SOURCES=Test_simd.cc +Test_simd_LDADD=-lGrid + + +Test_stencil_SOURCES=Test_stencil.cc +Test_stencil_LDADD=-lGrid + + +Test_wilson_cg_prec_SOURCES=Test_wilson_cg_prec.cc +Test_wilson_cg_prec_LDADD=-lGrid + + +Test_wilson_cg_schur_SOURCES=Test_wilson_cg_schur.cc +Test_wilson_cg_schur_LDADD=-lGrid + + +Test_wilson_cg_unprec_SOURCES=Test_wilson_cg_unprec.cc +Test_wilson_cg_unprec_LDADD=-lGrid + + +Test_wilson_even_odd_SOURCES=Test_wilson_even_odd.cc +Test_wilson_even_odd_LDADD=-lGrid + diff --git a/tests/Makefile.am b/tests/Makefile.am index 3c79191c..83385001 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -2,37 +2,4 @@ AM_CXXFLAGS = -I$(top_srcdir)/lib AM_LDFLAGS = -L$(top_builddir)/lib -# -# Test code -# -bin_PROGRAMS = Grid_main Grid_stencil Grid_nersc_io Grid_cshift Grid_gamma Grid_simd Grid_rng Grid_remez Grid_rng_fixed - -Grid_main_SOURCES = Grid_main.cc -Grid_main_LDADD = -lGrid - -Grid_rng_SOURCES = Grid_rng.cc -Grid_rng_LDADD = -lGrid - -Grid_rng_fixed_SOURCES = Grid_rng_fixed.cc -Grid_rng_fixed_LDADD = -lGrid - -Grid_remez_SOURCES = Grid_remez.cc -Grid_remez_LDADD = -lGrid - -Grid_nersc_io_SOURCES = Grid_nersc_io.cc -Grid_nersc_io_LDADD = -lGrid - -Grid_cshift_SOURCES = Grid_cshift.cc -Grid_cshift_LDADD = -lGrid - -Grid_gamma_SOURCES = Grid_gamma.cc -Grid_gamma_LDADD = -lGrid - -Grid_stencil_SOURCES = Grid_stencil.cc -Grid_stencil_LDADD = -lGrid - -Grid_simd_SOURCES = Grid_simd.cc -Grid_simd_LDADD = -lGrid - -#Grid_simd_new_SOURCES = Grid_simd_new.cc -#Grid_simd_new_LDADD = -lGrid +include Make.inc diff --git a/tests/Sqrt.gnu b/tests/Sqrt.gnu deleted file mode 100644 index ae56ab97..00000000 --- a/tests/Sqrt.gnu +++ /dev/null @@ -1,2 +0,0 @@ -f(x) = 6.81384+(-2.34645e-06/(x+0.000228091))+(-1.51593e-05/(x+0.00112084))+(-6.89254e-05/(x+0.003496))+(-0.000288983/(x+0.00954309))+(-0.00119277/(x+0.024928))+(-0.0050183/(x+0.0646627))+(-0.0226449/(x+0.171576))+(-0.123767/(x+0.491792))+(-1.1705/(x+1.78667))+(-102.992/(x+18.4866)); -f(x) = 0.14676+(0.00952992/(x+5.40933e-05))+(0.0115952/(x+0.000559699))+(0.0161824/(x+0.00203338))+(0.0243252/(x+0.00582831))+(0.0379533/(x+0.0154649))+(0.060699/(x+0.0401156))+(0.100345/(x+0.104788))+(0.178335/(x+0.286042))+(0.381586/(x+0.892189))+(1.42625/(x+4.38422)); diff --git a/tests/Test_cayley_cg.cc b/tests/Test_cayley_cg.cc new file mode 100644 index 00000000..4510d4b5 --- /dev/null +++ b/tests/Test_cayley_cg.cc @@ -0,0 +1,172 @@ +#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 + }; + +template +void TestCGinversions(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5); +template +void TestCGschur(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5); + +template +void TestCGunprec(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5); + +template +void TestCGprec(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5); + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + int threads = GridThread::GetThreads(); + std::cout << "Grid is setup to use "< seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + + LatticeGaugeField Umu(UGrid); random(RNG4,Umu); + std::vector U(4,UGrid); + + RealD mass=0.1; + RealD M5 =1.8; + std::cout <<"DomainWallFermion test"<(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + RealD b=1.5;// Scale factor b+c=2, b-c=1 + RealD c=0.5; + std::cout <<"MobiusFermion test"<(Dmob,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + std::cout <<"MobiusZolotarevFermion test"<(Dzolo,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + std::cout <<"ScaledShamirFermion test"<(Dsham,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + std::cout <<"ShamirZolotarevFermion test"<(Dshamz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + std::cout <<"OverlapWilsonCayleyTanhFermion test"<(Dov,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + std::cout <<"OverlapWilsonCayleyZolotarevFermion test"<(Dovz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + Grid_finalize(); +} +template +void TestCGinversions(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5) +{ + std::cout << "Testing unpreconditioned inverter"<(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5); + std::cout << "Testing red black preconditioned inverter"<(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5); + std::cout << "Testing red black Schur inverter"<(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5); +} + +template +void TestCGunprec(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5) +{ + LatticeFermion src (FGrid); random(*RNG5,src); + LatticeFermion result(FGrid); result=zero; + + HermitianOperator HermOp(Ddwf); + ConjugateGradient CG(1.0e-8,10000); + CG(HermOp,src,result); + +} +template +void TestCGprec(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5) +{ + LatticeFermion src (FGrid); random(*RNG5,src); + LatticeFermion src_o(FrbGrid); + LatticeFermion result_o(FrbGrid); + pickCheckerboard(Odd,src_o,src); + result_o=zero; + + HermitianCheckerBoardedOperator HermOpEO(Ddwf); + ConjugateGradient CG(1.0e-8,10000); + CG(HermOpEO,src_o,result_o); +} + + +template +void TestCGschur(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5) +{ + LatticeFermion src (FGrid); random(*RNG5,src); + LatticeFermion result(FGrid); result=zero; + + ConjugateGradient CG(1.0e-8,10000); + SchurRedBlackSolve SchurSolver(CG); + SchurSolver(Ddwf,src,result); +} diff --git a/tests/Test_cayley_even_odd.cc b/tests/Test_cayley_even_odd.cc new file mode 100644 index 00000000..df28981b --- /dev/null +++ b/tests/Test_cayley_even_odd.cc @@ -0,0 +1,245 @@ +#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 + }; + + +template +void TestWhat(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, GridParallelRNG *RNG5); + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + int threads = GridThread::GetThreads(); + std::cout << "Grid is setup to use "< seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + + LatticeGaugeField Umu(UGrid); random(RNG4,Umu); + std::vector U(4,UGrid); + + RealD mass=0.1; + RealD M5 =1.8; + std::cout <<"DomainWallFermion test"<(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + RealD b=1.5;// Scale factor b+c=2, b-c=1 + RealD c=0.5; + std::cout <<"MobiusFermion test"<(Dmob,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + std::cout <<"MobiusZolotarevFermion test"<(Dzolo,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + std::cout <<"ScaledShamirFermion test"<(Dsham,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + + std::cout <<"ShamirZolotarevFermion test"<(Dshamz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + std::cout <<"OverlapWilsonCayleyTanhFermion test"<(Dov,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + std::cout <<"OverlapWilsonCayleyZolotarevFermion test"<(Dovz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + Grid_finalize(); +} + +template +void TestWhat(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5) +{ + + LatticeFermion src (FGrid); random(*RNG5,src); + LatticeFermion phi (FGrid); random(*RNG5,phi); + LatticeFermion chi (FGrid); random(*RNG5,chi); + LatticeFermion result(FGrid); result=zero; + LatticeFermion ref(FGrid); ref=zero; + LatticeFermion tmp(FGrid); tmp=zero; + LatticeFermion err(FGrid); tmp=zero; + + LatticeFermion src_e (FrbGrid); + LatticeFermion src_o (FrbGrid); + LatticeFermion r_e (FrbGrid); + LatticeFermion r_o (FrbGrid); + LatticeFermion r_eo (FGrid); + LatticeFermion r_eeoo(FGrid); + + std::cout<<"=========================================================="< * = < chi | Deo^dag| phi> "< + +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 + }; + + +template +void TestCGinversions(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5); +template +void TestCGschur(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5); + +template +void TestCGunprec(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5); + +template +void TestCGprec(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5); + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + int threads = GridThread::GetThreads(); + std::cout << "Grid is setup to use "< seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + + LatticeGaugeField Umu(UGrid); random(RNG4,Umu); + std::vector U(4,UGrid); + + RealD mass=0.1; + RealD M5 =1.8; + + + std::cout <<"OverlapWilsonContFracTanhFermion test"<(Dcf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + std::cout <<"OverlapWilsonContFracZolotarevFermion test"<(Dcfz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + Grid_finalize(); +} +template +void TestCGinversions(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5) +{ + std::cout << "Testing unpreconditioned inverter"<(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5); + std::cout << "Testing red black preconditioned inverter"<(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5); + std::cout << "Testing red black Schur inverter"<(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5); +} + +template +void TestCGunprec(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5) +{ + LatticeFermion src (FGrid); random(*RNG5,src); + LatticeFermion result(FGrid); result=zero; + + HermitianOperator HermOp(Ddwf); + ConjugateGradient CG(1.0e-8,10000); + CG(HermOp,src,result); + +} +template +void TestCGprec(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5) +{ + LatticeFermion src (FGrid); random(*RNG5,src); + LatticeFermion src_o(FrbGrid); + LatticeFermion result_o(FrbGrid); + pickCheckerboard(Odd,src_o,src); + result_o=zero; + + HermitianCheckerBoardedOperator HermOpEO(Ddwf); + ConjugateGradient CG(1.0e-8,10000); + CG(HermOpEO,src_o,result_o); +} + + +template +void TestCGschur(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5) +{ + LatticeFermion src (FGrid); random(*RNG5,src); + LatticeFermion result(FGrid); result=zero; + + ConjugateGradient CG(1.0e-8,10000); + SchurRedBlackSolve SchurSolver(CG); + SchurSolver(Ddwf,src,result); +} diff --git a/tests/Test_contfrac_even_odd.cc b/tests/Test_contfrac_even_odd.cc new file mode 100644 index 00000000..e13c1189 --- /dev/null +++ b/tests/Test_contfrac_even_odd.cc @@ -0,0 +1,223 @@ +#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 + }; + + +template +void TestWhat(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, GridParallelRNG *RNG5); + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + int threads = GridThread::GetThreads(); + std::cout << "Grid is setup to use "< seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + + LatticeGaugeField Umu(UGrid); random(RNG4,Umu); + std::vector U(4,UGrid); + + RealD mass=0.1; + RealD M5 =1.8; + + std::cout <<"OverlapWilsonContFracTanhFermion test"<(Dcf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + std::cout <<"OverlapWilsonContFracZolotarevFermion test"<(Dcfz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + Grid_finalize(); +} + +template +void TestWhat(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5) +{ + + LatticeFermion src (FGrid); random(*RNG5,src); + LatticeFermion phi (FGrid); random(*RNG5,phi); + LatticeFermion chi (FGrid); random(*RNG5,chi); + LatticeFermion result(FGrid); result=zero; + LatticeFermion ref(FGrid); ref=zero; + LatticeFermion tmp(FGrid); tmp=zero; + LatticeFermion err(FGrid); tmp=zero; + + LatticeFermion src_e (FrbGrid); + LatticeFermion src_o (FrbGrid); + LatticeFermion r_e (FrbGrid); + LatticeFermion r_o (FrbGrid); + LatticeFermion r_eo (FGrid); + LatticeFermion r_eeoo(FGrid); + + std::cout<<"=========================================================="< * = < chi | Deo^dag| phi> "< -#include using namespace Grid; using namespace Grid::QCD; diff --git a/tests/Test_cshift_red_black.cc b/tests/Test_cshift_red_black.cc new file mode 100644 index 00000000..9ffa66b1 --- /dev/null +++ b/tests/Test_cshift_red_black.cc @@ -0,0 +1,164 @@ +#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 + <<" ["< + +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); + + const int Ls=8; + + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::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 RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + + LatticeFermion src(FGrid); random(RNG5,src); + LatticeFermion result(FGrid); result=zero; + LatticeGaugeField Umu(UGrid); random(RNG4,Umu); + + std::vector U(4,UGrid); + for(int mu=0;mu(Umu,mu); + } + + RealD mass=0.1; + RealD M5=1.8; + DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); + + LatticeFermion src_o(FrbGrid); + LatticeFermion result_o(FrbGrid); + pickCheckerboard(Odd,src_o,src); + result_o=zero; + + HermitianCheckerBoardedOperator HermOpEO(Ddwf); + ConjugateGradient CG(1.0e-8,10000); + CG(HermOpEO,src_o,result_o); + + Grid_finalize(); +} diff --git a/tests/Test_dwf_cg_schur.cc b/tests/Test_dwf_cg_schur.cc new file mode 100644 index 00000000..aac4d3fd --- /dev/null +++ b/tests/Test_dwf_cg_schur.cc @@ -0,0 +1,53 @@ +#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 + }; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + const int Ls=8; + + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::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 RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + + LatticeFermion src(FGrid); random(RNG5,src); + LatticeFermion result(FGrid); result=zero; + LatticeGaugeField Umu(UGrid); random(RNG4,Umu); + + std::vector U(4,UGrid); + for(int mu=0;mu(Umu,mu); + } + + RealD mass=0.1; + RealD M5=1.8; + DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); + + ConjugateGradient CG(1.0e-8,10000); + SchurRedBlackSolve SchurSolver(CG); + SchurSolver(Ddwf,src,result); + + Grid_finalize(); +} diff --git a/tests/Test_dwf_cg_unprec.cc b/tests/Test_dwf_cg_unprec.cc new file mode 100644 index 00000000..5c9e7ad3 --- /dev/null +++ b/tests/Test_dwf_cg_unprec.cc @@ -0,0 +1,53 @@ +#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 + }; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + const int Ls=8; + + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::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 RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + + LatticeFermion src(FGrid); random(RNG5,src); + LatticeFermion result(FGrid); result=zero; + LatticeGaugeField Umu(UGrid); random(RNG4,Umu); + + std::vector U(4,UGrid); + for(int mu=0;mu(Umu,mu); + } + + RealD mass=0.1; + RealD M5=1.8; + DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); + + HermitianOperator HermOp(Ddwf); + ConjugateGradient CG(1.0e-8,10000); + CG(HermOp,src,result); + + Grid_finalize(); +} diff --git a/tests/Test_dwf_even_odd.cc b/tests/Test_dwf_even_odd.cc new file mode 100644 index 00000000..ac47bbf9 --- /dev/null +++ b/tests/Test_dwf_even_odd.cc @@ -0,0 +1,207 @@ +#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 + }; + + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + int threads = GridThread::GetThreads(); + std::cout << "Grid is setup to use "< 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 phi (FGrid); random(RNG5,phi); + LatticeFermion chi (FGrid); random(RNG5,chi); + LatticeFermion result(FGrid); result=zero; + LatticeFermion ref(FGrid); ref=zero; + LatticeFermion tmp(FGrid); tmp=zero; + LatticeFermion err(FGrid); tmp=zero; + LatticeGaugeField Umu(UGrid); random(RNG4,Umu); + std::vector U(4,UGrid); + + // Only one non-zero (y) + Umu=zero; + for(int nn=0;nn0 ) + U[nn]=zero; + pokeIndex(Umu,U[nn],nn); + } + + RealD mass=0.1; + RealD M5 =1.8; + DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); + + LatticeFermion src_e (FrbGrid); + LatticeFermion src_o (FrbGrid); + LatticeFermion r_e (FrbGrid); + LatticeFermion r_o (FrbGrid); + LatticeFermion r_eo (FGrid); + LatticeFermion r_eeoo(FGrid); + + std::cout<<"=========================================================="< * = < chi | Deo^dag| phi> "< -#include using namespace std; using namespace Grid; diff --git a/tests/Grid_main.cc b/tests/Test_main.cc similarity index 97% rename from tests/Grid_main.cc rename to tests/Test_main.cc index c526e17a..724e1563 100644 --- a/tests/Grid_main.cc +++ b/tests/Test_main.cc @@ -56,6 +56,7 @@ int main (int argc, char ** argv) GridCartesian Fine(latt_size,simd_layout,mpi_layout); GridRedBlackCartesian rbFine(latt_size,simd_layout,mpi_layout); GridParallelRNG FineRNG(&Fine); + GridSerialRNG SerialRNG; FineRNG.SeedRandomDevice(); LatticeColourMatrix Foo(&Fine); @@ -83,6 +84,9 @@ int main (int argc, char ** argv) LatticeSpinMatrix sMat(&Fine); LatticeSpinColourMatrix scMat(&Fine); + LatticeLorentzColourMatrix lcMat(&Fine); + + LatticeComplex scalar(&Fine); LatticeReal rscalar(&Fine); LatticeReal iscalar(&Fine); @@ -99,12 +103,15 @@ int main (int argc, char ** argv) random(FineRNG,cMat); random(FineRNG,sMat); random(FineRNG,scMat); + random(FineRNG,lcMat); random(FineRNG,cVec); random(FineRNG,sVec); random(FineRNG,scVec); + fflush(stdout); + TComplex tr = trace(cmat); cVec = cMat * cVec; // LatticeColourVector = LatticeColourMatrix * LatticeColourVector @@ -206,7 +213,14 @@ int main (int argc, char ** argv) scm=transpose(scm); scm=transposeIndex<1>(scm); + + random(SerialRNG, cm); + std::cout << cm << std::endl; + cm = Ta(cm); + //cm = adj(cm); + TComplex tracecm= trace(cm); + std::cout << cm << " "<< tracecm << std::endl; // Foo = Foo+scalar; // LatticeColourMatrix+Scalar @@ -219,6 +233,10 @@ int main (int argc, char ** argv) LatticeComplex trscMat(&Fine); trscMat = trace(scMat); // Trace + // LatticeComplex trlcMat(&Fine); + // trlcMat = trace(lcMat); // Trace involving iVector - now generates error + + { // Peek-ology and Poke-ology, with a little app-ology TComplex c; ColourMatrix c_m; diff --git a/tests/Grid_nersc_io.cc b/tests/Test_nersc_io.cc similarity index 98% rename from tests/Grid_nersc_io.cc rename to tests/Test_nersc_io.cc index 6fe587a6..80d78291 100644 --- a/tests/Grid_nersc_io.cc +++ b/tests/Test_nersc_io.cc @@ -1,5 +1,4 @@ #include -#include using namespace std; using namespace Grid; diff --git a/tests/Grid_remez.cc b/tests/Test_remez.cc similarity index 100% rename from tests/Grid_remez.cc rename to tests/Test_remez.cc diff --git a/tests/Grid_rng.cc b/tests/Test_rng.cc similarity index 97% rename from tests/Grid_rng.cc rename to tests/Test_rng.cc index 97f6c6b7..1731b740 100644 --- a/tests/Grid_rng.cc +++ b/tests/Test_rng.cc @@ -1,5 +1,4 @@ #include -#include using namespace std; using namespace Grid; diff --git a/tests/Grid_rng_fixed.cc b/tests/Test_rng_fixed.cc similarity index 97% rename from tests/Grid_rng_fixed.cc rename to tests/Test_rng_fixed.cc index 01e3315b..c836c93f 100644 --- a/tests/Grid_rng_fixed.cc +++ b/tests/Test_rng_fixed.cc @@ -1,5 +1,4 @@ #include -#include using namespace std; using namespace Grid; diff --git a/tests/Grid_simd.cc b/tests/Test_simd.cc similarity index 99% rename from tests/Grid_simd.cc rename to tests/Test_simd.cc index db600fe4..f4858d06 100644 --- a/tests/Grid_simd.cc +++ b/tests/Test_simd.cc @@ -1,5 +1,4 @@ #include -#include using namespace std; using namespace Grid; diff --git a/tests/Grid_stencil.cc b/tests/Test_stencil.cc similarity index 100% rename from tests/Grid_stencil.cc rename to tests/Test_stencil.cc diff --git a/benchmarks/Grid_wilson_cg_prec.cc b/tests/Test_wilson_cg_prec.cc similarity index 88% rename from benchmarks/Grid_wilson_cg_prec.cc rename to tests/Test_wilson_cg_prec.cc index e86ec820..e376349c 100644 --- a/benchmarks/Grid_wilson_cg_prec.cc +++ b/tests/Test_wilson_cg_prec.cc @@ -42,9 +42,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); @@ -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/tests/Test_wilson_cg_schur.cc similarity index 96% rename from benchmarks/Grid_wilson_cg_schur.cc rename to tests/Test_wilson_cg_schur.cc index af630ae1..28db1d4b 100644 --- a/benchmarks/Grid_wilson_cg_schur.cc +++ b/tests/Test_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/tests/Test_wilson_cg_unprec.cc similarity index 92% rename from benchmarks/Grid_wilson_cg_unprec.cc rename to tests/Test_wilson_cg_unprec.cc index 15302aab..905dfde5 100644 --- a/benchmarks/Grid_wilson_cg_unprec.cc +++ b/tests/Test_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/tests/Test_wilson_even_odd.cc similarity index 97% rename from benchmarks/Grid_wilson_evenodd.cc rename to tests/Test_wilson_even_odd.cc index 4bc8c357..3ebc4709 100644 --- a/benchmarks/Grid_wilson_evenodd.cc +++ b/tests/Test_wilson_even_odd.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); @@ -68,8 +68,6 @@ int main (int argc, char ** argv) LatticeFermion r_o (&RBGrid); LatticeFermion r_eo (&Grid); - const int Even=0; - const int Odd=1; std::cout<<"=========================================================="<