From 155c164b0c25f35cb318ceb3b47da43307ac0ecb Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sat, 15 Aug 2015 23:25:49 +0100 Subject: [PATCH] * Finished the template/policy style introduction of gparity, except the gparity force terms. So valence sector looks ok. FermionOperatorImpl.h provides the policy classes. Expect HMC will introduce a smearing policy and a fermion representation change policy template param. Will also probably need multi-precision work. * HMC is running even-odd and non-checkerboarded (checked 4^4 wilson fermion/wilson gauge). There appears to be a bug in the multi-level integrator -- passes with single level but not with multi-level. In any case there looks to be quite a bit to clean up. This is the "const det" style implementation that is not appropriate yet for clover since it assumes that Mee is indept of the gauge fields. Easily fixed in future. --- lib/Make.inc | 2 +- lib/qcd/action/Actions.h | 4 +- lib/qcd/action/fermion/CayleyFermion5D.cc | 4 +- lib/qcd/action/fermion/CayleyFermion5D.h | 4 +- .../fermion/ContinuedFractionFermion5D.cc | 6 +- .../fermion/ContinuedFractionFermion5D.h | 4 +- lib/qcd/action/fermion/DomainWallFermion.h | 8 +- lib/qcd/action/fermion/FermionImplTypedefs.h | 11 - lib/qcd/action/fermion/FermionOperator.h | 5 +- lib/qcd/action/fermion/FermionOperatorImpl.h | 353 ++++++++++++------ lib/qcd/action/fermion/MobiusFermion.h | 12 +- .../action/fermion/MobiusZolotarevFermion.h | 12 +- .../fermion/OverlapWilsonCayleyTanhFermion.h | 6 +- .../OverlapWilsonCayleyZolotarevFermion.h | 12 +- .../OverlapWilsonContfracTanhFermion.h | 12 +- .../OverlapWilsonContfracZolotarevFermion.h | 12 +- .../OverlapWilsonPartialFractionTanhFermion.h | 6 +- ...lapWilsonPartialFractionZolotarevFermion.h | 16 +- .../fermion/PartialFractionFermion5D.cc | 15 +- .../action/fermion/PartialFractionFermion5D.h | 4 +- lib/qcd/action/fermion/ScaledShamirFermion.h | 2 +- .../action/fermion/ShamirZolotarevFermion.h | 12 +- lib/qcd/action/fermion/WilsonFermion.cc | 9 +- lib/qcd/action/fermion/WilsonFermion.h | 10 +- lib/qcd/action/fermion/WilsonFermion5D.cc | 23 +- lib/qcd/action/fermion/WilsonFermion5D.h | 8 +- lib/qcd/action/fermion/WilsonKernels.h | 54 ++- lib/qcd/action/pseudofermion/TwoFlavour.h | 3 +- .../action/pseudofermion/TwoFlavourEvenOdd.h | 41 +- lib/qcd/hmc/integrators/Integrator_base.h | 6 +- scripts/hmc.sh | 6 + scripts/substitute | 2 +- tests/Test_hmc_EOWilsonFermionGauge.cc | 10 +- tests/Test_hmc_WilsonFermionGauge.cc | 6 +- 34 files changed, 414 insertions(+), 286 deletions(-) delete mode 100644 lib/qcd/action/fermion/FermionImplTypedefs.h create mode 100755 scripts/hmc.sh diff --git a/lib/Make.inc b/lib/Make.inc index 8d595ae58..8da8d28e5 100644 --- a/lib/Make.inc +++ b/lib/Make.inc @@ -1,4 +1,4 @@ -HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/AdefGeneric.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/PrecConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/Preconditioner.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./Config.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./Init.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_unary.h ./lattice/Lattice_where.h ./Lattice.h ./Log.h ./MacroMagic.h ./Old/Tensor_peek.h ./Old/Tensor_poke.h ./parallelIO/NerscIO.h ./qcd/action/ActionBase.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionImplTypedefs.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/pseudofermion/TwoFlavour.h ./qcd/action/pseudofermion/TwoFlavourEvenOdd.h ./qcd/hmc/HMC.h ./qcd/hmc/integrators/Integrator.h ./qcd/hmc/integrators/Integrator_algorithm.h ./qcd/hmc/integrators/Integrator_base.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/WilsonLoops.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Grid_vector_unops.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_determinant.h ./tensors/Tensor_exp.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_index.h ./tensors/Tensor_inner.h ./tensors/Tensor_logical.h ./tensors/Tensor_outer.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./tensors/Tensor_unary.h ./Tensors.h ./Threads.h ./Timer.h +HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/AdefGeneric.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/PrecConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/Preconditioner.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./Config.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./Init.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_unary.h ./lattice/Lattice_where.h ./Lattice.h ./Log.h ./MacroMagic.h ./Old/Tensor_peek.h ./Old/Tensor_poke.h ./parallelIO/NerscIO.h ./qcd/action/ActionBase.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/FermionOperatorImpl.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/gauge/GaugeMap.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/pseudofermion/TwoFlavour.h ./qcd/action/pseudofermion/TwoFlavourEvenOdd.h ./qcd/hmc/HMC.h ./qcd/hmc/integrators/Integrator.h ./qcd/hmc/integrators/Integrator_algorithm.h ./qcd/hmc/integrators/Integrator_base.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/WilsonLoops.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Grid_vector_unops.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_determinant.h ./tensors/Tensor_exp.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_index.h ./tensors/Tensor_inner.h ./tensors/Tensor_logical.h ./tensors/Tensor_outer.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./tensors/Tensor_unary.h ./Tensors.h ./Threads.h ./Timer.h CCFILES=./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./Init.cc ./Log.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/hmc/HMC.cc ./qcd/hmc/integrators/Integrator.cc ./qcd/spin/Dirac.cc ./qcd/utils/SpaceTimeGrid.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h index cf44aa7a0..2c902d840 100644 --- a/lib/qcd/action/Actions.h +++ b/lib/qcd/action/Actions.h @@ -27,6 +27,8 @@ // Utility functions //////////////////////////////////////////// #include //used by all wilson type fermions +#include +#include #include //used by all wilson type fermions @@ -54,8 +56,6 @@ //////////////////////////////////////////// // Fermion operators / actions //////////////////////////////////////////// -#include -#include #include // 4d wilson like //#include diff --git a/lib/qcd/action/fermion/CayleyFermion5D.cc b/lib/qcd/action/fermion/CayleyFermion5D.cc index 76066adfc..39d04492f 100644 --- a/lib/qcd/action/fermion/CayleyFermion5D.cc +++ b/lib/qcd/action/fermion/CayleyFermion5D.cc @@ -8,12 +8,12 @@ namespace QCD { GridRedBlackCartesian &FiveDimRedBlackGrid, GridCartesian &FourDimGrid, GridRedBlackCartesian &FourDimRedBlackGrid, - RealD _mass,RealD _M5) : + RealD _mass,RealD _M5,const ImplParams &p) : WilsonFermion5D(_Umu, FiveDimGrid, FiveDimRedBlackGrid, FourDimGrid, - FourDimRedBlackGrid,_M5), + FourDimRedBlackGrid,_M5,p), mass(_mass) { } diff --git a/lib/qcd/action/fermion/CayleyFermion5D.h b/lib/qcd/action/fermion/CayleyFermion5D.h index c0d2a2922..c7eb6e807 100644 --- a/lib/qcd/action/fermion/CayleyFermion5D.h +++ b/lib/qcd/action/fermion/CayleyFermion5D.h @@ -9,7 +9,7 @@ namespace Grid { class CayleyFermion5D : public WilsonFermion5D { public: -#include + INHERIT_IMPL_TYPES(Impl); public: // override multiply @@ -64,7 +64,7 @@ namespace Grid { GridRedBlackCartesian &FiveDimRedBlackGrid, GridCartesian &FourDimGrid, GridRedBlackCartesian &FourDimRedBlackGrid, - RealD _mass,RealD _M5); + RealD _mass,RealD _M5,const ImplParams &p= ImplParams()); protected: void SetCoefficientsZolotarev(RealD zolohi,Approx::zolotarev_data *zdata,RealD b,RealD c); diff --git a/lib/qcd/action/fermion/ContinuedFractionFermion5D.cc b/lib/qcd/action/fermion/ContinuedFractionFermion5D.cc index 846e6a1b6..f562ff546 100644 --- a/lib/qcd/action/fermion/ContinuedFractionFermion5D.cc +++ b/lib/qcd/action/fermion/ContinuedFractionFermion5D.cc @@ -256,10 +256,10 @@ namespace Grid { GridRedBlackCartesian &FiveDimRedBlackGrid, GridCartesian &FourDimGrid, GridRedBlackCartesian &FourDimRedBlackGrid, - RealD _mass,RealD M5) : + RealD _mass,RealD M5,const ImplParams &p) : WilsonFermion5D(_Umu, - FiveDimGrid, FiveDimRedBlackGrid, - FourDimGrid, FourDimRedBlackGrid,M5), + FiveDimGrid, FiveDimRedBlackGrid, + FourDimGrid, FourDimRedBlackGrid,M5,p), mass(_mass) { int Ls = this->Ls; diff --git a/lib/qcd/action/fermion/ContinuedFractionFermion5D.h b/lib/qcd/action/fermion/ContinuedFractionFermion5D.h index b1a7fba27..b3f25ac57 100644 --- a/lib/qcd/action/fermion/ContinuedFractionFermion5D.h +++ b/lib/qcd/action/fermion/ContinuedFractionFermion5D.h @@ -9,7 +9,7 @@ namespace Grid { class ContinuedFractionFermion5D : public WilsonFermion5D { public: -#include + INHERIT_IMPL_TYPES(Impl); public: // override multiply @@ -41,7 +41,7 @@ namespace Grid { GridRedBlackCartesian &FiveDimRedBlackGrid, GridCartesian &FourDimGrid, GridRedBlackCartesian &FourDimRedBlackGrid, - RealD _mass,RealD M5); + RealD _mass,RealD M5,const ImplParams &p= ImplParams()); protected: diff --git a/lib/qcd/action/fermion/DomainWallFermion.h b/lib/qcd/action/fermion/DomainWallFermion.h index 00a6babab..4d088503a 100644 --- a/lib/qcd/action/fermion/DomainWallFermion.h +++ b/lib/qcd/action/fermion/DomainWallFermion.h @@ -11,7 +11,7 @@ namespace Grid { class DomainWallFermion : public CayleyFermion5D { public: -#include + INHERIT_IMPL_TYPES(Impl); public: virtual void Instantiatable(void) {}; @@ -21,13 +21,13 @@ namespace Grid { GridRedBlackCartesian &FiveDimRedBlackGrid, GridCartesian &FourDimGrid, GridRedBlackCartesian &FourDimRedBlackGrid, - RealD _mass,RealD _M5) : - + RealD _mass,RealD _M5,const ImplParams &p= ImplParams()) : + CayleyFermion5D(_Umu, FiveDimGrid, FiveDimRedBlackGrid, FourDimGrid, - FourDimRedBlackGrid,_mass,_M5) + FourDimRedBlackGrid,_mass,_M5,p) { RealD eps = 1.0; diff --git a/lib/qcd/action/fermion/FermionImplTypedefs.h b/lib/qcd/action/fermion/FermionImplTypedefs.h deleted file mode 100644 index 59f39be73..000000000 --- a/lib/qcd/action/fermion/FermionImplTypedefs.h +++ /dev/null @@ -1,11 +0,0 @@ - public: - typedef typename Impl::Simd Simd; - typedef typename Impl::FermionField FermionField; - typedef typename Impl::GaugeLinkField GaugeLinkField; - typedef typename Impl::GaugeField GaugeField; - typedef typename Impl::DoubledGaugeField DoubledGaugeField; - typedef typename Impl::SiteSpinor SiteSpinor; - typedef typename Impl::SiteHalfSpinor SiteHalfSpinor; - typedef typename Impl::Compressor Compressor; - typedef WilsonKernels Kernels; - diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h index 4903f0e10..31b6d8a9a 100644 --- a/lib/qcd/action/fermion/FermionOperator.h +++ b/lib/qcd/action/fermion/FermionOperator.h @@ -21,7 +21,10 @@ namespace Grid { class FermionOperator : public CheckerBoardedSparseMatrixBase, public Impl { public: -#include + + INHERIT_IMPL_TYPES(Impl); + + FermionOperator(const ImplParams &p= ImplParams()) : Impl(p) {}; public: diff --git a/lib/qcd/action/fermion/FermionOperatorImpl.h b/lib/qcd/action/fermion/FermionOperatorImpl.h index 38a12d323..24c64149a 100644 --- a/lib/qcd/action/fermion/FermionOperatorImpl.h +++ b/lib/qcd/action/fermion/FermionOperatorImpl.h @@ -5,35 +5,136 @@ namespace Grid { namespace QCD { + ////////////////////////////////////////////// + // Template parameter class constructs to package + // externally control Fermion implementations + // in orthogonal directions + // + // Ultimately need Impl to always define types where XXX is opaque + // + // typedef typename XXX Simd; + // typedef typename XXX GaugeLinkField; + // typedef typename XXX GaugeField; + // typedef typename XXX GaugeActField; + // typedef typename XXX FermionField; + // typedef typename XXX DoubledGaugeField; + // typedef typename XXX SiteSpinor; + // typedef typename XXX SiteHalfSpinor; + // typedef typename XXX Compressor; + // + // and Methods: + // void ImportGauge(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu) + // void DoubleStore(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu) + // void multLink(SiteHalfSpinor &phi,const SiteDoubledGaugeField &U,const SiteHalfSpinor &chi,int mu,StencilEntry *SE,CartesianStencil &St) + // void InsertForce4D(GaugeField &mat,const FermionField &Btilde,const FermionField &A,int mu) + // void InsertForce5D(GaugeField &mat,const FermionField &Btilde,const FermionField &A,int mu) + // + // + // To acquire the typedefs from "Base" (either a base class or template param) use: + // + // INHERIT_GIMPL_TYPES(Base) + // INHERIT_FIMPL_TYPES(Base) + // INHERIT_IMPL_TYPES(Base) + // + // The Fermion operators will do the following: + // + // struct MyOpParams { + // RealD mass; + // }; + // + // + // template + // class MyOp : pubic { + // public: + // + // INHERIT_ALL_IMPL_TYPES(Impl); + // + // MyOp(MyOpParams Myparm, ImplParams &ImplParam) : Impl(ImplParam) + // { + // + // }; + // + // } + ////////////////////////////////////////////// + + + //////////////////////////////////////////////////////////////////////// + // Implementation dependent gauge types + //////////////////////////////////////////////////////////////////////// + +#define INHERIT_IMPL_TYPES(Base) \ + INHERIT_GIMPL_TYPES(Base)\ + INHERIT_FIMPL_TYPES(Base) + +#define INHERIT_GIMPL_TYPES(GImpl) \ + typedef typename GImpl::Simd Simd;\ + typedef typename GImpl::GaugeLinkField GaugeLinkField;\ + typedef typename GImpl::GaugeField GaugeField; + + // Composition with smeared link, bc's etc.. probably need multiple inheritance // Variable precision "S" and variable Nc template - class WilsonImpl { + class ImplGauge { + public: + + typedef S Simd; + + template using iImplGaugeLink = iScalar > >; + template using iImplGaugeField = iVector >, Nd >; + + typedef iImplGaugeLink SiteGaugeLink; + typedef iImplGaugeField SiteGaugeField; + + typedef Lattice GaugeLinkField; // bit ugly naming; polarised gauge field, lorentz... all ugly + typedef Lattice GaugeField; + + }; + + //////////////////////////////////////////////////////////////////////// + // Implementation dependent fermion types + //////////////////////////////////////////////////////////////////////// + +#define INHERIT_FIMPL_TYPES(Impl)\ + typedef typename Impl::FermionField FermionField; \ + typedef typename Impl::DoubledGaugeField DoubledGaugeField; \ + typedef typename Impl::SiteSpinor SiteSpinor; \ + typedef typename Impl::SiteHalfSpinor SiteHalfSpinor; \ + typedef typename Impl::Compressor Compressor; \ + typedef typename Impl::ImplParams ImplParams; + + /////// + // Single flavour four spinors with colour index + /////// + template + class WilsonImpl : public ImplGauge { public: - typedef S Simd; + typedef ImplGauge Gimpl; + + + INHERIT_GIMPL_TYPES(Gimpl); template using iImplSpinor = iScalar, Ns> >; template using iImplHalfSpinor = iScalar, Nhs> >; - template using iImplGaugeLink = iScalar > >; - template using iImplGaugeField = iVector >, Nd >; template using iImplDoubledGaugeField = iVector >, Nds >; typedef iImplSpinor SiteSpinor; typedef iImplHalfSpinor SiteHalfSpinor; - typedef iImplGaugeLink SiteGaugeLink; - typedef iImplGaugeField SiteGaugeField; typedef iImplDoubledGaugeField SiteDoubledGaugeField; typedef Lattice FermionField; - typedef Lattice GaugeLinkField; // bit ugly naming; polarised gauge field, lorentz... all ugly - typedef Lattice GaugeField; typedef Lattice DoubledGaugeField; typedef WilsonCompressor Compressor; - static inline void multLink(SiteHalfSpinor &phi,const SiteDoubledGaugeField &U,const SiteHalfSpinor &chi,int mu,StencilEntry *SE,CartesianStencil &St){ + typedef struct WilsonImplParams { } ImplParams; + ImplParams Params; + WilsonImpl(const ImplParams &p= ImplParams()) : Params(p) {}; + + inline void multLink(SiteHalfSpinor &phi,const SiteDoubledGaugeField &U,const SiteHalfSpinor &chi,int mu,StencilEntry *SE,CartesianStencil &St){ mult(&phi(),&U(mu),&chi()); } + inline void DoubleStore(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu) { conformable(Uds._grid,GaugeGrid); @@ -46,52 +147,69 @@ namespace Grid { PokeIndex(Uds,U,mu+4); } } - - inline void InsertForce(GaugeField &mat,const FermionField &Btilde,const FermionField &A,int mu){ + + inline void InsertForce4D(GaugeField &mat, FermionField &Btilde, FermionField &A,int mu){ GaugeLinkField link(mat._grid); link = TraceIndex(outerProduct(Btilde,A)); PokeIndex(mat,link,mu); } + inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField Ã,int mu){ + + int Ls=Btilde._grid->_fdimensions[0]; + + GaugeLinkField tmp(mat._grid); + tmp = zero; +PARALLEL_FOR_LOOP + for(int sss=0;sssoSites();sss++){ + int sU=sss; + for(int s=0;s(outerProduct(Btilde[sF],Atilde[sF])); // ordering here + } + } + PokeIndex(mat,tmp,mu); + + } + }; - template - class GparityWilsonImpl { + //////////////////////////////////////////////////////////////////////////////////////// + // Flavour doubled spinors; is Gparity the only? what about C*? + //////////////////////////////////////////////////////////////////////////////////////// + template + class GparityWilsonImpl : public ImplGauge { public: - typedef S Simd; + typedef ImplGauge Gimpl; + + INHERIT_GIMPL_TYPES(Gimpl); template using iImplSpinor = iVector, Ns>, Ngp >; template using iImplHalfSpinor = iVector, Nhs>, Ngp >; - template using iImplGaugeField = iVector >, Nd >; - - template using iImplGaugeLink = iScalar > >; template using iImplDoubledGaugeField = iVector >, Nds >, Ngp >; typedef iImplSpinor SiteSpinor; typedef iImplHalfSpinor SiteHalfSpinor; - typedef iImplGaugeLink SiteGaugeLink; - typedef iImplGaugeField SiteGaugeField; typedef iImplDoubledGaugeField SiteDoubledGaugeField; typedef Lattice FermionField; - typedef Lattice GaugeLinkField; // bit ugly naming; polarised gauge field, lorentz... all ugly - typedef Lattice GaugeField; typedef Lattice DoubledGaugeField; - // typedef GparityWilsonCompressor Compressor; typedef WilsonCompressor Compressor; - // provide the multiply by link that is differentiated between Gparity (with flavour index) and - // non-Gparity - static inline void multLink(SiteHalfSpinor &phi,const SiteDoubledGaugeField &U,const SiteHalfSpinor &chi,int mu,StencilEntry *SE,CartesianStencil &St){ + typedef struct GparityWilsonImplParams {std::vector twists; } ImplParams; + ImplParams Params; + GparityWilsonImpl(const ImplParams &p= ImplParams()) : Params(p) {}; + + // provide the multiply by link that is differentiated between Gparity (with flavour index) and non-Gparity + inline void multLink(SiteHalfSpinor &phi,const SiteDoubledGaugeField &U,const SiteHalfSpinor &chi,int mu,StencilEntry *SE,CartesianStencil &St){ typedef SiteHalfSpinor vobj; typedef typename SiteHalfSpinor::scalar_object sobj; vobj vtmp; sobj stmp; - std::vector gpbc({0,0,0,1,0,0,0,1}); GridBase *grid = St._grid; @@ -99,125 +217,131 @@ namespace Grid { int direction = St._directions[mu]; int distance = St._distances[mu]; - int ptype = St._permute_type[mu]; - int sl = St._grid->_simd_layout[direction]; - + int ptype = St._permute_type[mu]; + int sl = St._grid->_simd_layout[direction]; + + // Fixme X.Y.Z.T hardcode in stencil + int mmu = mu % Nd; + // assert our assumptions assert((distance==1)||(distance==-1)); // nearest neighbour stencil hard code assert((sl==1)||(sl==2)); std::vector icoor; - if ( SE->_around_the_world && gpbc[mu] ) { + if ( SE->_around_the_world && Params.twists[mmu] ) { if ( sl == 2 ) { - std::vector vals(Nsimd); + std::vector vals(Nsimd); - extract(chi,vals); - for(int s=0;siCoorFromIindex(icoor,s); - - assert((icoor[direction]==0)||(icoor[direction]==1)); - - int permute_lane; - if ( distance == 1) { - permute_lane = icoor[direction]?1:0; - } else { - permute_lane = icoor[direction]?0:1; + grid->iCoorFromIindex(icoor,s); + + assert((icoor[direction]==0)||(icoor[direction]==1)); + + int permute_lane; + if ( distance == 1) { + permute_lane = icoor[direction]?1:0; + } else { + permute_lane = icoor[direction]?0:1; + } + + if ( permute_lane ) { + stmp(0) = vals[s](1); + stmp(1) = vals[s](0); + vals[s] = stmp; + } } + merge(vtmp,vals); - if ( permute_lane ) { - stmp(0) = vals[s](1); - stmp(1) = vals[s](0); - vals[s] = stmp; - } + } else { + vtmp(0) = chi(1); + vtmp(1) = chi(0); } - merge(vtmp,vals); - + mult(&phi(0),&U(0)(mu),&vtmp(0)); + mult(&phi(1),&U(1)(mu),&vtmp(1)); + } else { - vtmp(0) = chi(1); - vtmp(1) = chi(0); + mult(&phi(0),&U(0)(mu),&chi(0)); + mult(&phi(1),&U(1)(mu),&chi(1)); } - mult(&phi(0),&U(0)(mu),&vtmp(0)); - mult(&phi(1),&U(1)(mu),&vtmp(1)); - } else { - mult(&phi(0),&U(0)(mu),&chi(0)); - mult(&phi(1),&U(1)(mu),&chi(1)); } - } - - static inline void InsertForce(GaugeField &mat,const FermionField &Btilde,const FermionField &A,int mu){ - // Fixme - return; - } - static inline void DoubleStore(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu) - { - - conformable(Uds._grid,GaugeGrid); - conformable(Umu._grid,GaugeGrid); - - GaugeLinkField Utmp(GaugeGrid); - GaugeLinkField U(GaugeGrid); - GaugeLinkField Uconj(GaugeGrid); - - Lattice > coor(GaugeGrid); - - std::vector gpdirs({0,0,0,1}); - - for(int mu=0;mu(Umu,mu); - Uconj = conjugate(U); - - int neglink = GaugeGrid->GlobalDimensions()[mu]-1; + GaugeLinkField Utmp(GaugeGrid); + GaugeLinkField U(GaugeGrid); + GaugeLinkField Uconj(GaugeGrid); - if ( gpdirs[mu] ) { - Uconj = where(coor==neglink,-Uconj,Uconj); - } + Lattice > coor(GaugeGrid); + + for(int mu=0;mu(Umu,mu); + Uconj = conjugate(U); + + // This phase could come from a simple bc 1,1,-1,1 .. + int neglink = GaugeGrid->GlobalDimensions()[mu]-1; + if ( Params.twists[mu] ) { + Uconj = where(coor==neglink,-Uconj,Uconj); + } + + PARALLEL_FOR_LOOP - for(auto ss=U.begin();ss WilsonImplR; // Real.. whichever prec typedef WilsonImpl WilsonImplF; // Float typedef WilsonImpl WilsonImplD; // Double @@ -226,7 +350,6 @@ PARALLEL_FOR_LOOP typedef GparityWilsonImpl GparityWilsonImplF; // Float typedef GparityWilsonImpl GparityWilsonImplD; // Double - } } #endif diff --git a/lib/qcd/action/fermion/MobiusFermion.h b/lib/qcd/action/fermion/MobiusFermion.h index e84661dca..b8147ce56 100644 --- a/lib/qcd/action/fermion/MobiusFermion.h +++ b/lib/qcd/action/fermion/MobiusFermion.h @@ -11,7 +11,7 @@ namespace Grid { class MobiusFermion : public CayleyFermion5D { public: -#include + INHERIT_IMPL_TYPES(Impl); public: virtual void Instantiatable(void) {}; @@ -22,13 +22,13 @@ namespace Grid { GridCartesian &FourDimGrid, GridRedBlackCartesian &FourDimRedBlackGrid, RealD _mass,RealD _M5, - RealD b, RealD c) : + RealD b, RealD c,const ImplParams &p= ImplParams()) : CayleyFermion5D(_Umu, - FiveDimGrid, - FiveDimRedBlackGrid, - FourDimGrid, - FourDimRedBlackGrid,_mass,_M5) + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_mass,_M5,p) { RealD eps = 1.0; diff --git a/lib/qcd/action/fermion/MobiusZolotarevFermion.h b/lib/qcd/action/fermion/MobiusZolotarevFermion.h index 1c572326d..54428f561 100644 --- a/lib/qcd/action/fermion/MobiusZolotarevFermion.h +++ b/lib/qcd/action/fermion/MobiusZolotarevFermion.h @@ -11,7 +11,7 @@ namespace Grid { class MobiusZolotarevFermion : public CayleyFermion5D { public: -#include + INHERIT_IMPL_TYPES(Impl); public: virtual void Instantiatable(void) {}; @@ -23,13 +23,13 @@ namespace Grid { GridRedBlackCartesian &FourDimRedBlackGrid, RealD _mass,RealD _M5, RealD b, RealD c, - RealD lo, RealD hi) : + RealD lo, RealD hi,const ImplParams &p= ImplParams()) : CayleyFermion5D(_Umu, - FiveDimGrid, - FiveDimRedBlackGrid, - FourDimGrid, - FourDimRedBlackGrid,_mass,_M5) + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_mass,_M5,p) { RealD eps = lo/hi; diff --git a/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h b/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h index 0e8d4e19b..093475f32 100644 --- a/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h +++ b/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h @@ -11,7 +11,7 @@ namespace Grid { class OverlapWilsonCayleyTanhFermion : public MobiusFermion { public: -#include + INHERIT_IMPL_TYPES(Impl); public: // Constructors @@ -21,14 +21,14 @@ namespace Grid { GridCartesian &FourDimGrid, GridRedBlackCartesian &FourDimRedBlackGrid, RealD _mass,RealD _M5, - RealD scale) : + RealD scale,const ImplParams &p= ImplParams()) : // b+c=scale, b-c = 0 <=> b =c = scale/2 MobiusFermion(_Umu, FiveDimGrid, FiveDimRedBlackGrid, FourDimGrid, - FourDimRedBlackGrid,_mass,_M5,0.5*scale,0.5*scale) + FourDimRedBlackGrid,_mass,_M5,0.5*scale,0.5*scale,p) { } }; diff --git a/lib/qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h b/lib/qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h index 549292dbe..909dca02b 100644 --- a/lib/qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h +++ b/lib/qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h @@ -11,7 +11,7 @@ namespace Grid { class OverlapWilsonCayleyZolotarevFermion : public MobiusZolotarevFermion { public: -#include + INHERIT_IMPL_TYPES(Impl); public: // Constructors @@ -22,13 +22,13 @@ namespace Grid { GridCartesian &FourDimGrid, GridRedBlackCartesian &FourDimRedBlackGrid, RealD _mass,RealD _M5, - RealD lo, RealD hi) : + RealD lo, RealD hi,const ImplParams &p= ImplParams()) : // 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) + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_mass,_M5,0.5,0.5,lo,hi,p) {} diff --git a/lib/qcd/action/fermion/OverlapWilsonContfracTanhFermion.h b/lib/qcd/action/fermion/OverlapWilsonContfracTanhFermion.h index 73baffa7e..652f193ef 100644 --- a/lib/qcd/action/fermion/OverlapWilsonContfracTanhFermion.h +++ b/lib/qcd/action/fermion/OverlapWilsonContfracTanhFermion.h @@ -11,7 +11,7 @@ namespace Grid { class OverlapWilsonContFracTanhFermion : public ContinuedFractionFermion5D { public: -#include + INHERIT_IMPL_TYPES(Impl); public: virtual void Instantiatable(void){}; @@ -22,14 +22,14 @@ namespace Grid { GridCartesian &FourDimGrid, GridRedBlackCartesian &FourDimRedBlackGrid, RealD _mass,RealD _M5, - RealD scale) : + RealD scale,const ImplParams &p= ImplParams()) : // b+c=scale, b-c = 0 <=> b =c = scale/2 ContinuedFractionFermion5D(_Umu, - FiveDimGrid, - FiveDimRedBlackGrid, - FourDimGrid, - FourDimRedBlackGrid,_mass,_M5) + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_mass,_M5,p) { assert((this->Ls&0x1)==1); // Odd Ls required int nrational=this->Ls-1;// Even rational order diff --git a/lib/qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h b/lib/qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h index 3c662e920..7847df3f1 100644 --- a/lib/qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h +++ b/lib/qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h @@ -11,7 +11,7 @@ namespace Grid { class OverlapWilsonContFracZolotarevFermion : public ContinuedFractionFermion5D { public: -#include + INHERIT_IMPL_TYPES(Impl); public: virtual void Instantiatable(void){}; @@ -22,14 +22,14 @@ namespace Grid { GridCartesian &FourDimGrid, GridRedBlackCartesian &FourDimRedBlackGrid, RealD _mass,RealD _M5, - RealD lo,RealD hi): + RealD lo,RealD hi,const ImplParams &p= ImplParams()): // b+c=scale, b-c = 0 <=> b =c = scale/2 ContinuedFractionFermion5D(_Umu, - FiveDimGrid, - FiveDimRedBlackGrid, - FourDimGrid, - FourDimRedBlackGrid,_mass,_M5) + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_mass,_M5,p) { assert((this->Ls&0x1)==1); // Odd Ls required diff --git a/lib/qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h b/lib/qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h index b017a9d72..b5d89b7a5 100644 --- a/lib/qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h +++ b/lib/qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h @@ -11,7 +11,7 @@ namespace Grid { class OverlapWilsonPartialFractionTanhFermion : public PartialFractionFermion5D { public: -#include + INHERIT_IMPL_TYPES(Impl); public: virtual void Instantiatable(void){}; @@ -22,14 +22,14 @@ namespace Grid { GridCartesian &FourDimGrid, GridRedBlackCartesian &FourDimRedBlackGrid, RealD _mass,RealD _M5, - RealD scale) : + RealD scale,const ImplParams &p= ImplParams()) : // b+c=scale, b-c = 0 <=> b =c = scale/2 PartialFractionFermion5D(_Umu, FiveDimGrid, FiveDimRedBlackGrid, FourDimGrid, - FourDimRedBlackGrid,_mass,_M5) + FourDimRedBlackGrid,_mass,_M5,p) { assert((this->Ls&0x1)==1); // Odd Ls required int nrational=this->Ls-1;// Even rational order diff --git a/lib/qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h b/lib/qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h index 3294f0791..e80c97519 100644 --- a/lib/qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h +++ b/lib/qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h @@ -11,25 +11,25 @@ namespace Grid { class OverlapWilsonPartialFractionZolotarevFermion : public PartialFractionFermion5D { public: -#include + INHERIT_IMPL_TYPES(Impl); public: virtual void Instantiatable(void){}; // Constructors OverlapWilsonPartialFractionZolotarevFermion(GaugeField &_Umu, - GridCartesian &FiveDimGrid, - GridRedBlackCartesian &FiveDimRedBlackGrid, - GridCartesian &FourDimGrid, - GridRedBlackCartesian &FourDimRedBlackGrid, - RealD _mass,RealD _M5, - RealD lo,RealD hi): + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD _M5, + RealD lo,RealD hi,const ImplParams &p= ImplParams()): // b+c=scale, b-c = 0 <=> b =c = scale/2 PartialFractionFermion5D(_Umu, FiveDimGrid, FiveDimRedBlackGrid, FourDimGrid, - FourDimRedBlackGrid,_mass,_M5) + FourDimRedBlackGrid,_mass,_M5,p) { assert((this->Ls&0x1)==1); // Odd Ls required diff --git a/lib/qcd/action/fermion/PartialFractionFermion5D.cc b/lib/qcd/action/fermion/PartialFractionFermion5D.cc index 81afb653b..a008dc78a 100644 --- a/lib/qcd/action/fermion/PartialFractionFermion5D.cc +++ b/lib/qcd/action/fermion/PartialFractionFermion5D.cc @@ -369,14 +369,15 @@ namespace Grid { // Constructors template PartialFractionFermion5D::PartialFractionFermion5D(GaugeField &_Umu, - GridCartesian &FiveDimGrid, - GridRedBlackCartesian &FiveDimRedBlackGrid, - GridCartesian &FourDimGrid, - GridRedBlackCartesian &FourDimRedBlackGrid, - RealD _mass,RealD M5) : + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD M5, + const ImplParams &p) : WilsonFermion5D(_Umu, - FiveDimGrid, FiveDimRedBlackGrid, - FourDimGrid, FourDimRedBlackGrid,M5), + FiveDimGrid, FiveDimRedBlackGrid, + FourDimGrid, FourDimRedBlackGrid,M5,p), mass(_mass) { diff --git a/lib/qcd/action/fermion/PartialFractionFermion5D.h b/lib/qcd/action/fermion/PartialFractionFermion5D.h index acbc30e1d..d0f2718d9 100644 --- a/lib/qcd/action/fermion/PartialFractionFermion5D.h +++ b/lib/qcd/action/fermion/PartialFractionFermion5D.h @@ -9,7 +9,7 @@ namespace Grid { class PartialFractionFermion5D : public WilsonFermion5D { public: -#include + INHERIT_IMPL_TYPES(Impl); public: const int part_frac_chroma_convention=1; @@ -47,7 +47,7 @@ namespace Grid { GridRedBlackCartesian &FiveDimRedBlackGrid, GridCartesian &FourDimGrid, GridRedBlackCartesian &FourDimRedBlackGrid, - RealD _mass,RealD M5); + RealD _mass,RealD M5,const ImplParams &p= ImplParams()); protected: diff --git a/lib/qcd/action/fermion/ScaledShamirFermion.h b/lib/qcd/action/fermion/ScaledShamirFermion.h index de7ded0cd..2225fd4bf 100644 --- a/lib/qcd/action/fermion/ScaledShamirFermion.h +++ b/lib/qcd/action/fermion/ScaledShamirFermion.h @@ -11,7 +11,7 @@ namespace Grid { class ScaledShamirFermion : public MobiusFermion { public: -#include + INHERIT_IMPL_TYPES(Impl); public: // Constructors diff --git a/lib/qcd/action/fermion/ShamirZolotarevFermion.h b/lib/qcd/action/fermion/ShamirZolotarevFermion.h index 339a37f18..34d6fd2b4 100644 --- a/lib/qcd/action/fermion/ShamirZolotarevFermion.h +++ b/lib/qcd/action/fermion/ShamirZolotarevFermion.h @@ -11,7 +11,7 @@ namespace Grid { class ShamirZolotarevFermion : public MobiusZolotarevFermion { public: -#include + INHERIT_IMPL_TYPES(Impl); public: // Constructors @@ -23,14 +23,14 @@ namespace Grid { GridCartesian &FourDimGrid, GridRedBlackCartesian &FourDimRedBlackGrid, RealD _mass,RealD _M5, - RealD lo, RealD hi) : + RealD lo, RealD hi,const ImplParams &p= ImplParams()) : // b+c = 1; b-c = 1 => b=1, c=0 MobiusZolotarevFermion(_Umu, - FiveDimGrid, - FiveDimRedBlackGrid, - FourDimGrid, - FourDimRedBlackGrid,_mass,_M5,1.0,0.0,lo,hi) + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_mass,_M5,1.0,0.0,lo,hi,p) {} diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index ea8c2209f..caea055f4 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -15,7 +15,8 @@ namespace QCD { WilsonFermion::WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid, GridRedBlackCartesian &Hgrid, - RealD _mass) : + RealD _mass,const ImplParams &p) : + Kernels(p), _grid(&Fgrid), _cbgrid(&Hgrid), Stencil (&Fgrid,npoint,Even,directions,displacements), @@ -118,7 +119,9 @@ namespace QCD { Compressor compressor(dag); FermionField Btilde(B._grid); - + FermionField Atilde(B._grid); + Atilde = A; + st.HaloExchange(B,comm_buf,compressor); for(int mu=0;mu - class WilsonFermion : public FermionOperator, public WilsonFermionStatic + class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { -#include - public: + INHERIT_IMPL_TYPES(Impl); + typedef WilsonKernels Kernels; /////////////////////////////////////////////////////////////// // Implement the abstract base @@ -88,7 +88,9 @@ namespace Grid { WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid, GridRedBlackCartesian &Hgrid, - RealD _mass) ; + RealD _mass, + const ImplParams &p= ImplParams() + ) ; // DoubleStore impl dependent void ImportGauge(const GaugeField &_Umu); diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 10caad139..66ca67d51 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -15,7 +15,8 @@ WilsonFermion5D::WilsonFermion5D(GaugeField &_Umu, GridRedBlackCartesian &FiveDimRedBlackGrid, GridCartesian &FourDimGrid, GridRedBlackCartesian &FourDimRedBlackGrid, - RealD _M5) : + RealD _M5,const ImplParams &p) : + Kernels(p), _FiveDimGrid(&FiveDimGrid), _FiveDimRedBlackGrid(&FiveDimRedBlackGrid), _FourDimGrid(&FourDimGrid), @@ -105,11 +106,11 @@ PARALLEL_FOR_LOOP template void WilsonFermion5D::DerivInternal(CartesianStencil & st, - DoubledGaugeField & U, - GaugeField &mat, - const FermionField &A, - const FermionField &B, - int dag) + DoubledGaugeField & U, + GaugeField &mat, + const FermionField &A, + const FermionField &B, + int dag) { assert((dag==DaggerNo) ||(dag==DaggerYes)); @@ -118,7 +119,6 @@ void WilsonFermion5D::DerivInternal(CartesianStencil & st, Compressor compressor(dag); - GaugeLinkField tmp(mat._grid); FermionField Btilde(B._grid); FermionField Atilde(B._grid); @@ -157,14 +157,7 @@ PARALLEL_FOR_LOOP } - Impl::InsertForce(mat,Btilde,A,mu); - /* - tmp = zero; - for(int sss=0;sssoSites();sss++){ - tmp[sU] = tmp[sU]+ traceIndex(outerProduct(Btilde[sF],Atilde[sF])); // ordering here - } - PokeIndex(mat,tmp,mu); - */ + Impl::InsertForce5D(mat,Btilde,Atilde,mu); } } diff --git a/lib/qcd/action/fermion/WilsonFermion5D.h b/lib/qcd/action/fermion/WilsonFermion5D.h index 5159d9e1d..ba84617c2 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.h +++ b/lib/qcd/action/fermion/WilsonFermion5D.h @@ -26,9 +26,11 @@ namespace Grid { }; template - class WilsonFermion5D : public FermionOperator, public WilsonFermion5DStatic + class WilsonFermion5D : public WilsonKernels, public WilsonFermion5DStatic { -#include + INHERIT_IMPL_TYPES(Impl); + typedef WilsonKernels Kernels; + public: /////////////////////////////////////////////////////////////// // Implement the abstract base @@ -90,7 +92,7 @@ namespace Grid { GridRedBlackCartesian &FiveDimRedBlackGrid, GridCartesian &FourDimGrid, GridRedBlackCartesian &FourDimRedBlackGrid, - double _M5); + double _M5,const ImplParams &p= ImplParams()); // DoubleStore void ImportGauge(const GaugeField &_Umu); diff --git a/lib/qcd/action/fermion/WilsonKernels.h b/lib/qcd/action/fermion/WilsonKernels.h index eb0981fc4..9e14a3a33 100644 --- a/lib/qcd/action/fermion/WilsonKernels.h +++ b/lib/qcd/action/fermion/WilsonKernels.h @@ -7,43 +7,41 @@ namespace Grid { //////////////////////////////////////////////////////////////////////////////////////////////////////////////// // Helper routines that implement Wilson stencil for a single site. + // Common to both the WilsonFermion and WilsonFermion5D //////////////////////////////////////////////////////////////////////////////////////////////////////////////// - template - class WilsonKernels { + template class WilsonKernels : public FermionOperator { public: -#include - + INHERIT_IMPL_TYPES(Impl); + typedef FermionOperator Base; + public: - static void - DiracOptDhopSite(CartesianStencil &st,DoubledGaugeField &U, - std::vector > &buf, - int sF,int sU,const FermionField &in, FermionField &out); + void DiracOptDhopSite(CartesianStencil &st,DoubledGaugeField &U, + std::vector > &buf, + int sF,int sU,const FermionField &in, FermionField &out); - static void - DiracOptDhopSiteDag(CartesianStencil &st,DoubledGaugeField &U, - std::vector > &buf, - int sF,int sU,const FermionField &in,FermionField &out); + void DiracOptDhopSiteDag(CartesianStencil &st,DoubledGaugeField &U, + std::vector > &buf, + int sF,int sU,const FermionField &in,FermionField &out); - static void - DiracOptDhopDir(CartesianStencil &st,DoubledGaugeField &U, - std::vector > &buf, - int sF,int sU,const FermionField &in, FermionField &out,int dirdisp,int gamma); + void DiracOptDhopDir(CartesianStencil &st,DoubledGaugeField &U, + std::vector > &buf, + int sF,int sU,const FermionField &in, FermionField &out,int dirdisp,int gamma); - static void - DiracOptHandDhopSite(CartesianStencil &st,DoubledGaugeField &U, - std::vector > &buf, - int sF,int sU,const FermionField &in, FermionField &out){ - DiracOptDhopSite(st,U,buf,sF,sU,in,out); // will template override for Wilson Nc=3 - } + void DiracOptHandDhopSite(CartesianStencil &st,DoubledGaugeField &U, + std::vector > &buf, + int sF,int sU,const FermionField &in, FermionField &out){ + DiracOptDhopSite(st,U,buf,sF,sU,in,out); // will template override for Wilson Nc=3 + } - static void - DiracOptHandDhopSiteDag(CartesianStencil &st,DoubledGaugeField &U, - std::vector > &buf, - int sF,int sU,const FermionField &in, FermionField &out){ - DiracOptDhopSiteDag(st,U,buf,sF,sU,in,out); // will template override for Wilson Nc=3 - } + void DiracOptHandDhopSiteDag(CartesianStencil &st,DoubledGaugeField &U, + std::vector > &buf, + int sF,int sU,const FermionField &in, FermionField &out){ + DiracOptDhopSiteDag(st,U,buf,sF,sU,in,out); // will template override for Wilson Nc=3 + } + + WilsonKernels(const ImplParams &p= ImplParams()) : Base(p) {}; }; diff --git a/lib/qcd/action/pseudofermion/TwoFlavour.h b/lib/qcd/action/pseudofermion/TwoFlavour.h index 635c1ca6a..c04231a7f 100644 --- a/lib/qcd/action/pseudofermion/TwoFlavour.h +++ b/lib/qcd/action/pseudofermion/TwoFlavour.h @@ -97,9 +97,10 @@ namespace Grid{ //////////////////////////////////////////////////////////////////////// template class TwoFlavourPseudoFermionAction : public Action { + public: + INHERIT_IMPL_TYPES(Impl); private: -#include FermionOperator & FermOp;// the basic operator diff --git a/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h b/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h index 35a208c88..3cbf6ffc7 100644 --- a/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h +++ b/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h @@ -4,14 +4,21 @@ namespace Grid{ namespace QCD{ + // Base even odd HMC on the normal Mee based schur decomposition. + // + // M = (Mee Meo) = (1 0 ) (Mee 0 ) (1 Mee^{-1} Meo) + // (Moe Moo) (Moe Mee^-1 1 ) (0 Moo-Moe Mee^-1 Meo) (0 1 ) + // + // Determinant is det of middle factor + // This assumes Mee is indept of U. + // template class SchurDifferentiableOperator : public SchurDiagMooeeOperator,typename Impl::FermionField> { public: -#include - - public: - typedef FermionOperator Matrix; + INHERIT_IMPL_TYPES(Impl); + + typedef FermionOperator Matrix; SchurDifferentiableOperator (Matrix &Mat) : SchurDiagMooeeOperator(Mat) {}; @@ -19,10 +26,11 @@ namespace Grid{ GridBase *fgrid = this->_Mat.FermionGrid(); GridBase *fcbgrid = this->_Mat.FermionRedBlackGrid(); - GridBase *ugrid = this->_Mat.GaugeGrid(); + GridBase *ugrid = this->_Mat.GaugeGrid(); GridBase *ucbgrid = this->_Mat.GaugeRedBlackGrid(); Real coeff = 1.0; + FermionField tmp1(fcbgrid); FermionField tmp2(fcbgrid); @@ -31,7 +39,7 @@ namespace Grid{ // Assert the checkerboard?? or code for either assert(U.checkerboard==Odd); - assert(V.checkerboard==V.checkerboard); + assert(V.checkerboard==U.checkerboard); GaugeField ForceO(ucbgrid); GaugeField ForceE(ucbgrid); @@ -47,6 +55,9 @@ namespace Grid{ this->_Mat.MooeeInvDag(tmp1,tmp2); // even->even this->_Mat.MeoDeriv(ForceE,tmp2,V,DaggerNo); + assert(ForceE.checkerboard==Even); + assert(ForceO.checkerboard==Odd); + setCheckerboard(Force,ForceE); setCheckerboard(Force,ForceO); Force=-Force; @@ -61,6 +72,7 @@ namespace Grid{ GridBase *ucbgrid = this->_Mat.GaugeRedBlackGrid(); Real coeff = 1.0; + FermionField tmp1(fcbgrid); FermionField tmp2(fcbgrid); @@ -85,6 +97,9 @@ namespace Grid{ this->_Mat.MooeeInv(tmp1,tmp2); // even->even this->_Mat.MeoDeriv(ForceE,tmp2,V,DaggerYes); + assert(ForceE.checkerboard==Even); + assert(ForceO.checkerboard==Odd); + setCheckerboard(Force,ForceE); setCheckerboard(Force,ForceO); Force=-Force; @@ -100,7 +115,7 @@ namespace Grid{ class TwoFlavourEvenOddPseudoFermionAction : public Action { public: -#include + INHERIT_IMPL_TYPES(Impl); private: @@ -180,8 +195,8 @@ namespace Grid{ // The EE factorised block; normally can replace with zero if det is constant (gauge field indept) // Only really clover term that creates this. - FermOp.MooeeInvDag(PhiEven,Y); - action = action + norm2(Y); + // FermOp.MooeeInvDag(PhiEven,Y); + // action = action + norm2(Y); std::cout << GridLogMessage << "Pseudofermion EO action "<* ActPtr; // now force the same size as the rest of the code - typedef std::vector ActionLevel; + typedef std::vector ActionLevel; typedef std::vector ActionSet; typedef std::vector ObserverList; @@ -45,10 +45,6 @@ namespace Grid{ void generate_momenta_su3(LatticeLorentzColourMatrix&,GridParallelRNG&); } - - - - template< class IntegratorPolicy > class Integrator{ private: diff --git a/scripts/hmc.sh b/scripts/hmc.sh new file mode 100755 index 000000000..a8ebd3d0e --- /dev/null +++ b/scripts/hmc.sh @@ -0,0 +1,6 @@ +#!/bin/bash + +LOG=$1 +SWEEPS=$2 +grep Plaq $LOG | tail -n $SWEEPS | awk '{ S=S+$10} END { print S/NR} ' +grep dH log_cb | tail -n $SWEEPS | awk '{ S=S+exp(-$10)} END { print S/NR} ' diff --git a/scripts/substitute b/scripts/substitute index fab833c77..7607a115f 100755 --- a/scripts/substitute +++ b/scripts/substitute @@ -10,7 +10,7 @@ shift while (( "$#" )); do echo $1 -sed -e "s/${WAS}/${IS}/g" -i .bak $1 +sed -e "s:${WAS}:${IS}:g" -i .bak $1 shift diff --git a/tests/Test_hmc_EOWilsonFermionGauge.cc b/tests/Test_hmc_EOWilsonFermionGauge.cc index 66ec118b8..fa81fecef 100644 --- a/tests/Test_hmc_EOWilsonFermionGauge.cc +++ b/tests/Test_hmc_EOWilsonFermionGauge.cc @@ -35,16 +35,14 @@ int main (int argc, char ** argv) //Collect actions ActionLevel Level1; Level1.push_back(&WilsonNf2); - ActionLevel Level2(3); - Level2.push_back(&Waction); + Level1.push_back(&Waction); ActionSet FullSet; FullSet.push_back(Level1); - FullSet.push_back(Level2); // Create integrator - typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm - // typedef LeapFrog IntegratorAlgorithm;// change here to change the algorithm - IntegratorParameters MDpar(12,20,1.0); + // typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm + typedef LeapFrog IntegratorAlgorithm;// change here to change the algorithm + IntegratorParameters MDpar(12,40,1.0); std::vector rel ={1}; Integrator MDynamics(&Fine,MDpar, FullSet); diff --git a/tests/Test_hmc_WilsonFermionGauge.cc b/tests/Test_hmc_WilsonFermionGauge.cc index ef5980629..4e4fb3a54 100644 --- a/tests/Test_hmc_WilsonFermionGauge.cc +++ b/tests/Test_hmc_WilsonFermionGauge.cc @@ -35,17 +35,15 @@ int main (int argc, char ** argv) //Collect actions ActionLevel Level1; Level1.push_back(&WilsonNf2); - ActionLevel Level2(3); - Level2.push_back(&Waction); + Level1.push_back(&Waction); ActionSet FullSet; FullSet.push_back(Level1); - FullSet.push_back(Level2); + // FullSet.push_back(Level2); // Create integrator typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm // typedef LeapFrog IntegratorAlgorithm;// change here to change the algorithm IntegratorParameters MDpar(12,20,1.0); - std::vector rel ={1}; Integrator MDynamics(&Fine,MDpar, FullSet); // Create HMC