diff --git a/lib/Make.inc b/lib/Make.inc index f3f2edbf..dd9dce12 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/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/WilsonGaugeAction.h ./qcd/action/pseudofermion/TwoFlavour.h ./qcd/action/pseudofermion/TwoFlavourEvenOdd.h ./qcd/action/pseudofermion/TwoFlavourRatio.h ./qcd/hmc/HMC.h ./qcd/hmc/integrators/Integrator.h ./qcd/hmc/integrators/Integrator_algorithm.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/WilsonGaugeAction.h ./qcd/action/pseudofermion/EvenOddSchurDifferentiable.h ./qcd/action/pseudofermion/TwoFlavour.h ./qcd/action/pseudofermion/TwoFlavourEvenOdd.h ./qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h ./qcd/action/pseudofermion/TwoFlavourRatio.h ./qcd/hmc/HMC.h ./qcd/hmc/integrators/Integrator.h ./qcd/hmc/integrators/Integrator_algorithm.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 331400a4..bc9526ca 100644 --- a/lib/qcd/action/Actions.h +++ b/lib/qcd/action/Actions.h @@ -150,8 +150,16 @@ typedef DomainWallFermion GparityDomainWallFermionD; //////////////////////////////////////// // Pseudo fermion combinations for HMC //////////////////////////////////////// +#include #include -#include #include +#include +#include + +//Todo: RHMC +//#include +//#include +//#include +//#include #endif diff --git a/lib/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h b/lib/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h new file mode 100644 index 00000000..065d1a3b --- /dev/null +++ b/lib/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h @@ -0,0 +1,112 @@ +#ifndef QCD_EVEN_ODD_SCHUR_DIFFERENTIABLE_H +#define QCD_EVEN_ODD_SCHUR_DIFFERENTIABLE_H + +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: + INHERIT_IMPL_TYPES(Impl); + + typedef FermionOperator Matrix; + + SchurDifferentiableOperator (Matrix &Mat) : SchurDiagMooeeOperator(Mat) {}; + + void MpcDeriv(GaugeField &Force,const FermionField &U,const FermionField &V) { + + GridBase *fgrid = this->_Mat.FermionGrid(); + GridBase *fcbgrid = this->_Mat.FermionRedBlackGrid(); + GridBase *ugrid = this->_Mat.GaugeGrid(); + GridBase *ucbgrid = this->_Mat.GaugeRedBlackGrid(); + + Real coeff = 1.0; + + FermionField tmp1(fcbgrid); + FermionField tmp2(fcbgrid); + + conformable(fcbgrid,U._grid); + conformable(fcbgrid,V._grid); + + // Assert the checkerboard?? or code for either + assert(U.checkerboard==Odd); + assert(V.checkerboard==U.checkerboard); + + GaugeField ForceO(ucbgrid); + GaugeField ForceE(ucbgrid); + + // X^dag Der_oe MeeInv Meo Y + // Use Mooee as nontrivial but gauge field indept + this->_Mat.Meooe (V,tmp1); // odd->even -- implicit -0.5 factor to be applied + this->_Mat.MooeeInv(tmp1,tmp2); // even->even + this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerNo); + + // Accumulate X^dag M_oe MeeInv Der_eo Y + this->_Mat.MeooeDag (U,tmp1); // even->odd -- implicit -0.5 factor to be applied + 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; + } + + + void MpcDagDeriv(GaugeField &Force,const FermionField &U,const FermionField &V) { + + GridBase *fgrid = this->_Mat.FermionGrid(); + GridBase *fcbgrid = this->_Mat.FermionRedBlackGrid(); + GridBase *ugrid = this->_Mat.GaugeGrid(); + GridBase *ucbgrid = this->_Mat.GaugeRedBlackGrid(); + + Real coeff = 1.0; + + FermionField tmp1(fcbgrid); + FermionField tmp2(fcbgrid); + + conformable(fcbgrid,U._grid); + conformable(fcbgrid,V._grid); + + // Assert the checkerboard?? or code for either + assert(V.checkerboard==Odd); + assert(V.checkerboard==V.checkerboard); + + GaugeField ForceO(ucbgrid); + GaugeField ForceE(ucbgrid); + + // X^dag Der_oe MeeInv Meo Y + // Use Mooee as nontrivial but gauge field indept + this->_Mat.MeooeDag (V,tmp1); // odd->even -- implicit -0.5 factor to be applied + this->_Mat.MooeeInvDag(tmp1,tmp2); // even->even + this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerYes); + + // Accumulate X^dag M_oe MeeInv Der_eo Y + this->_Mat.Meooe (U,tmp1); // even->odd -- implicit -0.5 factor to be applied + 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; + } + + }; + + } +} +#endif diff --git a/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h b/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h index 9ea2de69..f07e86c4 100644 --- a/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h +++ b/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h @@ -4,108 +4,6 @@ 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: - INHERIT_IMPL_TYPES(Impl); - - typedef FermionOperator Matrix; - - SchurDifferentiableOperator (Matrix &Mat) : SchurDiagMooeeOperator(Mat) {}; - - void MpcDeriv(GaugeField &Force,const FermionField &U,const FermionField &V) { - - GridBase *fgrid = this->_Mat.FermionGrid(); - GridBase *fcbgrid = this->_Mat.FermionRedBlackGrid(); - GridBase *ugrid = this->_Mat.GaugeGrid(); - GridBase *ucbgrid = this->_Mat.GaugeRedBlackGrid(); - - Real coeff = 1.0; - - FermionField tmp1(fcbgrid); - FermionField tmp2(fcbgrid); - - conformable(fcbgrid,U._grid); - conformable(fcbgrid,V._grid); - - // Assert the checkerboard?? or code for either - assert(U.checkerboard==Odd); - assert(V.checkerboard==U.checkerboard); - - GaugeField ForceO(ucbgrid); - GaugeField ForceE(ucbgrid); - - // X^dag Der_oe MeeInv Meo Y - // Use Mooee as nontrivial but gauge field indept - this->_Mat.Meooe (V,tmp1); // odd->even -- implicit -0.5 factor to be applied - this->_Mat.MooeeInv(tmp1,tmp2); // even->even - this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerNo); - - // Accumulate X^dag M_oe MeeInv Der_eo Y - this->_Mat.MeooeDag (U,tmp1); // even->odd -- implicit -0.5 factor to be applied - 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; - } - - - void MpcDagDeriv(GaugeField &Force,const FermionField &U,const FermionField &V) { - - GridBase *fgrid = this->_Mat.FermionGrid(); - GridBase *fcbgrid = this->_Mat.FermionRedBlackGrid(); - GridBase *ugrid = this->_Mat.GaugeGrid(); - GridBase *ucbgrid = this->_Mat.GaugeRedBlackGrid(); - - Real coeff = 1.0; - - FermionField tmp1(fcbgrid); - FermionField tmp2(fcbgrid); - - conformable(fcbgrid,U._grid); - conformable(fcbgrid,V._grid); - - // Assert the checkerboard?? or code for either - assert(V.checkerboard==Odd); - assert(V.checkerboard==V.checkerboard); - - GaugeField ForceO(ucbgrid); - GaugeField ForceE(ucbgrid); - - // X^dag Der_oe MeeInv Meo Y - // Use Mooee as nontrivial but gauge field indept - this->_Mat.MeooeDag (V,tmp1); // odd->even -- implicit -0.5 factor to be applied - this->_Mat.MooeeInvDag(tmp1,tmp2); // even->even - this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerYes); - - // Accumulate X^dag M_oe MeeInv Der_eo Y - this->_Mat.Meooe (U,tmp1); // even->odd -- implicit -0.5 factor to be applied - 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; - } - - }; //////////////////////////////////////////////////////////////////////// @@ -123,7 +21,6 @@ namespace Grid{ FermionOperator & FermOp;// the basic operator OperatorFunction &DerivativeSolver; - OperatorFunction &ActionSolver; FermionField PhiOdd; // the pseudo fermion field for this trajectory @@ -154,6 +51,7 @@ namespace Grid{ // P(eta) = e^{- eta^dag eta} // // e^{x^2/2 sig^2} => sig^2 = 0.5. + RealD scale = std::sqrt(0.5); FermionField eta (FermOp.FermionGrid()); @@ -169,6 +67,7 @@ namespace Grid{ FermOp.ImportGauge(U); PCop.MpcDag(etaOdd,PhiOdd); + FermOp.MooeeDag(etaEven,PhiEven); PhiOdd =PhiOdd*scale; @@ -219,17 +118,16 @@ namespace Grid{ FermionField Y(FermOp.FermionRedBlackGrid()); GaugeField tmp(FermOp.GaugeGrid()); - SchurDifferentiableOperator PCop(FermOp); - - X=zero; - DerivativeSolver(PCop,PhiOdd,X); - PCop.Op(X,Y); + SchurDifferentiableOperator Mpc(FermOp); // Our conventions really make this UdSdU; We do not differentiate wrt Udag here. // So must take dSdU - adj(dSdU) and left multiply by mom to get dS/dt. - PCop.MpcDeriv(tmp , Y, X ); dSdU=tmp; - PCop.MpcDagDeriv(tmp , X, Y); dSdU=dSdU+tmp; + X=zero; + DerivativeSolver(Mpc,PhiOdd,X); + Mpc.Mpc(X,Y); + Mpc.MpcDeriv(tmp , Y, X ); dSdU=tmp; + Mpc.MpcDagDeriv(tmp , X, Y); dSdU=dSdU+tmp; // Treat the EE case. (MdagM)^-1 = Minv Minvdag // Deriv defaults to zero. diff --git a/lib/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h b/lib/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h new file mode 100644 index 00000000..935cec10 --- /dev/null +++ b/lib/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h @@ -0,0 +1,163 @@ +#ifndef QCD_PSEUDOFERMION_TWO_FLAVOUR_EVEN_ODD_RATIO_H +#define QCD_PSEUDOFERMION_TWO_FLAVOUR_EVEN_ODD_RATIO_H + +namespace Grid{ + namespace QCD{ + + /////////////////////////////////////// + // Two flavour ratio + /////////////////////////////////////// + template + class TwoFlavourEvenOddRatioPseudoFermionAction : public Action { + public: + INHERIT_IMPL_TYPES(Impl); + + private: + FermionOperator & NumOp;// the basic operator + FermionOperator & DenOp;// the basic operator + + OperatorFunction &DerivativeSolver; + OperatorFunction &ActionSolver; + + FermionField PhiOdd; // the pseudo fermion field for this trajectory + FermionField PhiEven; // the pseudo fermion field for this trajectory + + public: + TwoFlavourEvenOddRatioPseudoFermionAction(FermionOperator &_NumOp, + FermionOperator &_DenOp, + OperatorFunction & DS, + OperatorFunction & AS) : + NumOp(_NumOp), + DenOp(_DenOp), + DerivativeSolver(DS), + ActionSolver(AS), + PhiEven(_NumOp.FermionRedBlackGrid()), + PhiOdd(_NumOp.FermionRedBlackGrid()) + { + conformable(_NumOp.FermionGrid(), _DenOp.FermionGrid()); + conformable(_NumOp.FermionRedBlackGrid(), _DenOp.FermionRedBlackGrid()); + conformable(_NumOp.GaugeGrid(), _DenOp.GaugeGrid()); + conformable(_NumOp.GaugeRedBlackGrid(), _DenOp.GaugeRedBlackGrid()); + }; + + virtual void init(const GaugeField &U, GridParallelRNG& pRNG) { + + // P(phi) = e^{- phi^dag Vpc (MpcdagMpc)^-1 Vpcdag phi} + // + // NumOp == V + // DenOp == M + // + // Take phi_o = Vpcdag^{-1} Mpcdag eta_o ; eta_o = Mpcdag^{-1} Vpcdag Phi + // + // P(eta_o) = e^{- eta_o^dag eta_o} + // + // e^{x^2/2 sig^2} => sig^2 = 0.5. + // + RealD scale = std::sqrt(0.5); + + FermionField eta (NumOp.FermionGrid()); + FermionField etaOdd (NumOp.FermionRedBlackGrid()); + FermionField etaEven(NumOp.FermionRedBlackGrid()); + FermionField tmp (NumOp.FermionRedBlackGrid()); + + gaussian(pRNG,eta); + pickCheckerboard(Even,etaEven,eta); + pickCheckerboard(Odd,etaOdd,eta); + + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + SchurDifferentiableOperator Mpc(DenOp); + SchurDifferentiableOperator Vpc(NumOp); + + // Odd det factors + Mpc.MpcDag(etaOdd,PhiOdd); + ActionSolver(Vpc,PhiOdd,tmp); + Vpc.Mpc(tmp,PhiOdd); + + // Even det factors + DenOp.MooeeDag(etaEven,tmp); + NumOp.MooeeInvDag(tmp,PhiEven); + + PhiOdd =PhiOdd*scale; + PhiEven=PhiEven*scale; + + }; + + ////////////////////////////////////////////////////// + // S = phi^dag V (Mdag M)^-1 Vdag phi + ////////////////////////////////////////////////////// + virtual RealD S(const GaugeField &U) { + + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + SchurDifferentiableOperator Mpc(DenOp); + SchurDifferentiableOperator Vpc(NumOp); + + FermionField X(NumOp.FermionRedBlackGrid()); + FermionField Y(NumOp.FermionRedBlackGrid()); + + X=zero; + Vpc.MpcDag(PhiOdd,Y); // Y= Vdag phi + ActionSolver(Mpc,Y,X); // X= (MdagM)^-1 Vdag phi + Mpc.Mpc(X,Y); // Y= Mdag^-1 Vdag phi + + RealD action = norm2(Y); + + // The EE factorised block; normally can replace with zero if det is constant (gauge field indept) + // Only really clover term that creates this. Leave the EE portion as a future to do to make most + // rapid progresss on DWF for now. + // + // Vpc.MooeeDag(PhiEven,X); + // Mpc.MooeeInvDag(X,Y); + // action = action + norm2(Y); + + return action; + }; + + ////////////////////////////////////////////////////// + // dS/du = phi^dag dV (Mdag M)^-1 V^dag phi + // - phi^dag V (Mdag M)^-1 [ Mdag dM + dMdag M ] (Mdag M)^-1 V^dag phi + // + phi^dag V (Mdag M)^-1 dV^dag phi + ////////////////////////////////////////////////////// + virtual void deriv(const GaugeField &U,GaugeField & dSdU) { + + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + SchurDifferentiableOperator Mpc(DenOp); + SchurDifferentiableOperator Vpc(NumOp); + + FermionField X(NumOp.FermionRedBlackGrid()); + FermionField Y(NumOp.FermionRedBlackGrid()); + + GaugeField force(NumOp.GaugeGrid()); + + X=zero; + + //Y=Vdag phi + //X = (Mdag M)^-1 V^dag phi + //Y = (Mdag)^-1 V^dag phi + Vpc.MpcDag(PhiOdd,Y); // Y= Vdag phi + DerivativeSolver(Mpc,Y,X); // X= (MdagM)^-1 Vdag phi + Mpc.Mpc(X,Y); // Y= Mdag^-1 Vdag phi + + // phi^dag V (Mdag M)^-1 dV^dag phi + Vpc.MpcDagDeriv(force , X, PhiOdd ); dSdU=force; + + // phi^dag dV (Mdag M)^-1 V^dag phi + Vpc.MpcDeriv(force , PhiOdd, X ); dSdU=dSdU+force; + + // - phi^dag V (Mdag M)^-1 Mdag dM (Mdag M)^-1 V^dag phi + // - phi^dag V (Mdag M)^-1 dMdag M (Mdag M)^-1 V^dag phi + Mpc.MpcDeriv(force,Y,X); dSdU=dSdU-force; + Mpc.MpcDagDeriv(force,X,Y); dSdU=dSdU-force; + + dSdU = -Ta(dSdU); + + }; + }; + } +} +#endif diff --git a/lib/qcd/action/pseudofermion/TwoFlavourRatio.h b/lib/qcd/action/pseudofermion/TwoFlavourRatio.h index 5cfb8c38..da15e629 100644 --- a/lib/qcd/action/pseudofermion/TwoFlavourRatio.h +++ b/lib/qcd/action/pseudofermion/TwoFlavourRatio.h @@ -1,10 +1,5 @@ -<<<<<<< HEAD #ifndef QCD_PSEUDOFERMION_TWO_FLAVOUR_RATIO_H #define QCD_PSEUDOFERMION_TWO_FLAVOUR_RATIO_H -======= -#ifndef QCD_PSEUDOFERMION_TWO_FLAVOUR_H -#define QCD_PSEUDOFERMION_TWO_FLAVOUR_H ->>>>>>> ef6a9e6b07b80aea909c0a62f223fa3e66f53b3a namespace Grid{ namespace QCD{ @@ -107,7 +102,6 @@ namespace Grid{ FermionField X(NumOp.FermionGrid()); FermionField Y(NumOp.FermionGrid()); - FermionField f1(NumOp.FermionGrid()); GaugeField force(NumOp.GaugeGrid()); @@ -130,8 +124,8 @@ namespace Grid{ // - phi^dag V (Mdag M)^-1 dMdag M (Mdag M)^-1 V^dag phi DenOp.MDeriv(force,Y,X,DaggerNo); dSdU=dSdU-force; DenOp.MDeriv(force,X,Y,DaggerYes); dSdU=dSdU-force; - dSdU = - dSdU; - dSdU = Ta(dSdU); + + dSdU = - Ta(dSdU); }; }; diff --git a/tests/Make.inc b/tests/Make.inc index bfacb980..125778fe 100644 --- a/tests/Make.inc +++ b/tests/Make.inc @@ -1,5 +1,5 @@ -bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_cheby Test_contfrac_cg Test_contfrac_even_odd Test_contfrac_force Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_force Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_gparity Test_gpwilson_even_odd Test_hmc_EOWilsonFermionGauge Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonRatio Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_partfrac_force Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi +bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_cheby Test_contfrac_cg Test_contfrac_even_odd Test_contfrac_force Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_force Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_gparity Test_gpwilson_even_odd Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonRatio Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_partfrac_force Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi Test_GaugeAction_SOURCES=Test_GaugeAction.cc @@ -102,6 +102,10 @@ Test_hmc_EOWilsonFermionGauge_SOURCES=Test_hmc_EOWilsonFermionGauge.cc Test_hmc_EOWilsonFermionGauge_LDADD=-lGrid +Test_hmc_EOWilsonRatio_SOURCES=Test_hmc_EOWilsonRatio.cc +Test_hmc_EOWilsonRatio_LDADD=-lGrid + + Test_hmc_WilsonFermionGauge_SOURCES=Test_hmc_WilsonFermionGauge.cc Test_hmc_WilsonFermionGauge_LDADD=-lGrid diff --git a/tests/Test_hmc_EOWilsonRatio.cc b/tests/Test_hmc_EOWilsonRatio.cc new file mode 100644 index 00000000..6a4e2bbb --- /dev/null +++ b/tests/Test_hmc_EOWilsonRatio.cc @@ -0,0 +1,56 @@ +#include "Grid.h" + + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(4,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + + GridCartesian Fine(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout); + GridParallelRNG pRNG(&Fine); + pRNG.SeedRandomDevice(); + LatticeLorentzColourMatrix U(&Fine); + + SU3::HotConfiguration(pRNG, U); + + // simplify template declaration? Strip the lorentz from the second template + WilsonGaugeAction Waction(5.6); + + Real mass=-0.77; + Real pv =0.0; + WilsonFermionR DenOp(U,Fine,RBFine,mass); + WilsonFermionR NumOp(U,Fine,RBFine,pv); + + ConjugateGradient CG(1.0e-8,10000); + + TwoFlavourEvenOddRatioPseudoFermionAction WilsonNf2(NumOp, DenOp,CG,CG); + + //Collect actions + ActionLevel Level1; + Level1.push_back(&WilsonNf2); + Level1.push_back(&Waction); + ActionSet FullSet; + FullSet.push_back(Level1); + + + // 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); + Integrator MDynamics(&Fine,MDpar, FullSet); + + // Create HMC + HMCparameters HMCpar; + HybridMonteCarlo HMC(HMCpar, MDynamics); + + HMC.evolve(U); + +}