diff --git a/lib/Make.inc b/lib/Make.inc index 8d595ae5..8da8d28e 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 cf44aa7a..2c902d84 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 76066adf..39d04492 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 c0d2a292..c7eb6e80 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 846e6a1b..f562ff54 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 b1a7fba2..b3f25ac5 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 00a6baba..4d088503 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 59f39be7..00000000 --- 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 4903f0e1..31b6d8a9 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 38a12d32..24c64149 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 e84661dc..b8147ce5 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 1c572326..54428f56 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 0e8d4e19..093475f3 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 549292db..909dca02 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 73baffa7..652f193e 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 3c662e92..7847df3f 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 b017a9d7..b5d89b7a 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 3294f079..e80c9751 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 81afb653..a008dc78 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 acbc30e1..d0f2718d 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 de7ded0c..2225fd4b 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 339a37f1..34d6fd2b 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 ea8c2209..caea055f 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 10caad13..66ca67d5 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 5159d9e1..ba84617c 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 eb0981fc..9e14a3a3 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 635c1ca6..c04231a7 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 35a208c8..3cbf6ffc 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 00000000..a8ebd3d0 --- /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 fab833c7..7607a115 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 66ec118b..fa81fece 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 ef598062..4e4fb3a5 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