From afeabe0d23865517ee8ed9a84659f33b194070d5 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sun, 16 Aug 2015 00:14:10 +0100 Subject: [PATCH 01/11] Tidying --- lib/qcd/action/Actions.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h index 0abec0ec..cb8ef664 100644 --- a/lib/qcd/action/Actions.h +++ b/lib/qcd/action/Actions.h @@ -6,12 +6,10 @@ // are separating the concept of the operator from that of action. // // The FermAction contains methods to create -// // * Linear operators (Hermitian and non-hermitian) .. my LinearOperator // * System solvers (Hermitian and non-hermitian) .. my OperatorFunction // * MultiShift System solvers (Hermitian and non-hermitian) .. my OperatorFunction - //////////////////////////////////////////// // Abstract base interface //////////////////////////////////////////// From 353d66def1d725b075ea0962d66d7a68b450cfbe Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sun, 16 Aug 2015 01:41:05 +0100 Subject: [PATCH 02/11] Unused apparently --- lib/qcd/hmc/integrators/Integrator_base.h | 220 ---------------------- 1 file changed, 220 deletions(-) delete mode 100644 lib/qcd/hmc/integrators/Integrator_base.h diff --git a/lib/qcd/hmc/integrators/Integrator_base.h b/lib/qcd/hmc/integrators/Integrator_base.h deleted file mode 100644 index ea088f8d..00000000 --- a/lib/qcd/hmc/integrators/Integrator_base.h +++ /dev/null @@ -1,220 +0,0 @@ -//-------------------------------------------------------------------- -/*! @file Integrator_base.h - * @brief Declaration of classes for the abstract Molecular Dynamics integrator - * - * @author Guido Cossu - */ -//-------------------------------------------------------------------- - -#ifndef INTEGRATOR_INCLUDED -#define INTEGRATOR_INCLUDED - -#include - -class Observer; - - -/*! @brief Abstract base class for Molecular Dynamics management */ - -namespace Grid{ - namespace QCD{ - - typedef Action* ActPtr; // now force the same size as the rest of the code - typedef std::vector ActionLevel; - typedef std::vector ActionSet; - typedef std::vector ObserverList; - - class LeapFrog; - - - struct IntegratorParameters{ - int Nexp; - int MDsteps; // number of outer steps - RealD trajL; // trajectory length - RealD stepsize; - - IntegratorParameters(int Nexp_, - int MDsteps_, - RealD trajL_): - Nexp(Nexp_),MDsteps(MDsteps_),trajL(trajL_),stepsize(trajL/MDsteps){}; - }; - - - namespace MDutils{ - void generate_momenta(LatticeLorentzColourMatrix&,GridParallelRNG&); - void generate_momenta_su3(LatticeLorentzColourMatrix&,GridParallelRNG&); - } - - template< class IntegratorPolicy > - class Integrator{ - private: - IntegratorParameters Params; - const ActionSet as; - const std::vector Nrel; //relative step size per level - //ObserverList observers; // not yet - std::unique_ptr P; - - IntegratorPolicy TheIntegrator;// contains parameters too - - void update_P(LatticeLorentzColourMatrix&U, int level,double ep){ - for(int a=0; aderiv(U,force); - - Complex dSdt=0.0; - for(int mu=0;mu(force,mu); - mommu=PeekIndex(*P,mu); - - dSdt += sum(trace(forcemu*(*P))); - - } - std::cout << GridLogMessage << " action "<(U, mu); - Pmu=PeekIndex(*P, mu); - Umu = expMat(Pmu, Complex(ep, 0.0))*Umu; - } - - } - - void register_observers(); - void notify_observers(); - - friend void IntegratorPolicy::step (LatticeLorentzColourMatrix& U, - int level, std::vector& clock, - Integrator* Integ); - public: - Integrator(IntegratorParameters Par, - ActionSet& Aset, std::vector Nrel_): - Params(Par),as(Aset),Nrel(Nrel_){ - assert(as.size() == Nrel.size()); - }; - - ~Integrator(){} - - - //Initialization of momenta and actions - void init(LatticeLorentzColourMatrix& U, - GridParallelRNG& pRNG){ - std::cout<init(U, pRNG); - } - } - } - - - - RealD S(LatticeLorentzColourMatrix& U){ - // Momenta - LatticeComplex Hloc = - trace((*P)*adj(*P)); - Complex Hsum = sum(Hloc); - - RealD H = Hsum.real(); - - // Actions - for(int level=0; levelS(U); - - return H; - } - - - - void integrate(LatticeLorentzColourMatrix& U, int level){ - std::vector clock; - clock.resize(as.size(),0); - for(int step=0; step< Params.MDsteps; ++step) // MD step - TheIntegrator.step(U,0,clock, *(this)); - } - }; - - - class MinimumNorm2{ - const double lambda = 0.1931833275037836; - public: - void step (LatticeLorentzColourMatrix& U, int level, std::vector& clock); - - }; - - class LeapFrog{ - public: - void step (LatticeLorentzColourMatrix& U, - int level, std::vector& clock, - Integrator* Integ){ - // cl : current level - // fl : final level - // eps : current step size - - int fl = Integ->as.size() -1; - double eps = Integ->Params.stepsize; - - // Get current level step size - for(int l=0; l<=level; ++l) eps/= Integ->Nrel[l]; - - int fin = 1; - for(int l=0; l<=level; ++l) fin*= Integ->Nrel[l]; - fin = 2*Integ->Params.MDsteps*fin - 1; - - for(int e=0; eNrel[level]; ++e){ - - if(clock[level] == 0){ // initial half step - Integ->update_P(U, level,eps/2); - ++clock[level]; - for(int l=0; lupdate_U(U, eps); - for(int l=0; lupdate_P(U, level,eps/2); - - ++clock[level]; - for(int l=0; lupdate_P(U, level,eps); - - clock[level]+=2; - for(int l=0; l Date: Mon, 17 Aug 2015 23:14:48 +0100 Subject: [PATCH 03/11] Adding PV pseudofermion in prep for DWF HMC. Not compiled this yet, but cloned in from BFM. --- .../action/pseudofermion/TwoFlavourRatio.h | 132 ++++++++++++++++++ 1 file changed, 132 insertions(+) create mode 100644 lib/qcd/action/pseudofermion/TwoFlavourRatio.h diff --git a/lib/qcd/action/pseudofermion/TwoFlavourRatio.h b/lib/qcd/action/pseudofermion/TwoFlavourRatio.h new file mode 100644 index 00000000..1098c85a --- /dev/null +++ b/lib/qcd/action/pseudofermion/TwoFlavourRatio.h @@ -0,0 +1,132 @@ +#ifndef QCD_PSEUDOFERMION_TWO_FLAVOUR_H +#define QCD_PSEUDOFERMION_TWO_FLAVOUR_H + +namespace Grid{ + namespace QCD{ + + /////////////////////////////////////// + // Two flavour ratio + /////////////////////////////////////// + template + class TwoFlavourRatioPseudoFermionAction : public Action { + public: + INHERIT_IMPL_TYPES(Impl); + + private: + FermionOperator & NumOp;// the basic operator + FermionOperator & DenOp;// the basic operator + + OperatorFunction &DerivativeSolver; + OperatorFunction &ActionSolver; + + FermionField Phi; // the pseudo fermion field for this trajectory + + public: + TwoFlavourRatioPseudoFermionAction(FermionOperator &_NumOp, + FermionOperator &_DenOp, + OperatorFunction & DS, + OperatorFunction & AS + ) : NumOp(_NumOp), DenOp(_DenOp), DerivativeSolver(DS), ActionSolver(AS), Phi(Op.FermionGrid()) {}; + + virtual void init(const GaugeField &U, GridParallelRNG& pRNG) { + + // P(phi) = e^{- phi^dag V (MdagM)^-1 Vdag phi} + // + // phi = Vdag^{-1} Mdag eta + // + // P(eta) = e^{- eta^dag eta} + // + // e^{x^2/2 sig^2} => sig^2 = 0.5. + // + // So eta should be of width sig = 1/sqrt(2). + // and must multiply by 0.707.... + // + RealD scale = std::sqrt(0.5); + + FermionField eta(NumOp.FermionGrid()); + + gaussian(pRNG,eta); + + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + MdagMLinearOperator ,FermionField> MdagMOp(NumOp); + + DenOp.Mdag(eta,Phi); // Mdag eta + ActionSolver(MdagMOp,Phi,tmp); // (VdagV)^-1 Mdag eta = V^-1 Vdag^-1 Mdag eta + NumOp.M(tmp,Phi); // Vdag^-1 Mdag eta + + Phi=Phi*scale; + + }; + + ////////////////////////////////////////////////////// + // S = phi^dag V (Mdag M)^-1 Vdag phi + ////////////////////////////////////////////////////// + virtual RealD S(const GaugeField &U) { + + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + FermionField X(NumOp.FermionGrid()); + FermionField Y(NumOp.FermionGrid()); + + MdagMLinearOperator ,FermionField> MdagMOp(DenOp); + + X=zero; + NumOp.Mdag(Phi,Y); // Vdag phi + ActionSolver(MdagMOp,Y,X); // MdagMinv Vdag phi + MdagMOp.Op(X,Y); // Y=Mdaginv Vdag phi + + RealD 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); + + MdagMLinearOperator ,FermionField> MdagMOp(DenOp); + + FermionField X(NumOp.FermionGrid()); + FermionField Y(NumOp.FermionGrid()); + FermionField f1(NumOp.FermionGrid()); + + GaugeField force(FermOp.GaugeGrid()); + + X=zero; + + //f1=Vdag phi + NumOp.Mdag(phi,f1); + + //X = (Mdag M)^-1 V^dag phi + DerivativeSolver(MdagMOp,f1,X); + + //Y = (Mdag)^-1 V^dag phi + DenOp.M(X,Y); + + // phi^dag V (Mdag M)^-1 dV^dag phi + NumOp.MDeriv(force , X, Phi, DaggerYes ); dSdU=force; + + // phi^dag dV (Mdag M)^-1 V^dag phi + NumOp.MDeriv(force , Phi, X ,DaggerNo ); 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 + DenOp.MDeriv(force,Y,X,DaggerNo); dSdU=dSdU-force; + DenOp.MDeriv(force,X,Y,DaggerYes); dSdU=dSdU-force; + + dSdU = Ta(dSdU); + + }; + }; + } +} +#endif From 0bc38a69ce70f961ac56bb8d67816c3daced970c Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Tue, 18 Aug 2015 09:19:42 +0100 Subject: [PATCH 04/11] Adding PV pseudofermion in prep for DWF HMC. Not compiled this yet, but cloned in from BFM. --- .../action/pseudofermion/TwoFlavourRatio.h | 132 ++++++++++++++++++ 1 file changed, 132 insertions(+) create mode 100644 lib/qcd/action/pseudofermion/TwoFlavourRatio.h diff --git a/lib/qcd/action/pseudofermion/TwoFlavourRatio.h b/lib/qcd/action/pseudofermion/TwoFlavourRatio.h new file mode 100644 index 00000000..1098c85a --- /dev/null +++ b/lib/qcd/action/pseudofermion/TwoFlavourRatio.h @@ -0,0 +1,132 @@ +#ifndef QCD_PSEUDOFERMION_TWO_FLAVOUR_H +#define QCD_PSEUDOFERMION_TWO_FLAVOUR_H + +namespace Grid{ + namespace QCD{ + + /////////////////////////////////////// + // Two flavour ratio + /////////////////////////////////////// + template + class TwoFlavourRatioPseudoFermionAction : public Action { + public: + INHERIT_IMPL_TYPES(Impl); + + private: + FermionOperator & NumOp;// the basic operator + FermionOperator & DenOp;// the basic operator + + OperatorFunction &DerivativeSolver; + OperatorFunction &ActionSolver; + + FermionField Phi; // the pseudo fermion field for this trajectory + + public: + TwoFlavourRatioPseudoFermionAction(FermionOperator &_NumOp, + FermionOperator &_DenOp, + OperatorFunction & DS, + OperatorFunction & AS + ) : NumOp(_NumOp), DenOp(_DenOp), DerivativeSolver(DS), ActionSolver(AS), Phi(Op.FermionGrid()) {}; + + virtual void init(const GaugeField &U, GridParallelRNG& pRNG) { + + // P(phi) = e^{- phi^dag V (MdagM)^-1 Vdag phi} + // + // phi = Vdag^{-1} Mdag eta + // + // P(eta) = e^{- eta^dag eta} + // + // e^{x^2/2 sig^2} => sig^2 = 0.5. + // + // So eta should be of width sig = 1/sqrt(2). + // and must multiply by 0.707.... + // + RealD scale = std::sqrt(0.5); + + FermionField eta(NumOp.FermionGrid()); + + gaussian(pRNG,eta); + + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + MdagMLinearOperator ,FermionField> MdagMOp(NumOp); + + DenOp.Mdag(eta,Phi); // Mdag eta + ActionSolver(MdagMOp,Phi,tmp); // (VdagV)^-1 Mdag eta = V^-1 Vdag^-1 Mdag eta + NumOp.M(tmp,Phi); // Vdag^-1 Mdag eta + + Phi=Phi*scale; + + }; + + ////////////////////////////////////////////////////// + // S = phi^dag V (Mdag M)^-1 Vdag phi + ////////////////////////////////////////////////////// + virtual RealD S(const GaugeField &U) { + + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + FermionField X(NumOp.FermionGrid()); + FermionField Y(NumOp.FermionGrid()); + + MdagMLinearOperator ,FermionField> MdagMOp(DenOp); + + X=zero; + NumOp.Mdag(Phi,Y); // Vdag phi + ActionSolver(MdagMOp,Y,X); // MdagMinv Vdag phi + MdagMOp.Op(X,Y); // Y=Mdaginv Vdag phi + + RealD 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); + + MdagMLinearOperator ,FermionField> MdagMOp(DenOp); + + FermionField X(NumOp.FermionGrid()); + FermionField Y(NumOp.FermionGrid()); + FermionField f1(NumOp.FermionGrid()); + + GaugeField force(FermOp.GaugeGrid()); + + X=zero; + + //f1=Vdag phi + NumOp.Mdag(phi,f1); + + //X = (Mdag M)^-1 V^dag phi + DerivativeSolver(MdagMOp,f1,X); + + //Y = (Mdag)^-1 V^dag phi + DenOp.M(X,Y); + + // phi^dag V (Mdag M)^-1 dV^dag phi + NumOp.MDeriv(force , X, Phi, DaggerYes ); dSdU=force; + + // phi^dag dV (Mdag M)^-1 V^dag phi + NumOp.MDeriv(force , Phi, X ,DaggerNo ); 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 + DenOp.MDeriv(force,Y,X,DaggerNo); dSdU=dSdU-force; + DenOp.MDeriv(force,X,Y,DaggerYes); dSdU=dSdU-force; + + dSdU = Ta(dSdU); + + }; + }; + } +} +#endif From 6212807a7742f4378e289300c28b33581d41fad5 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Tue, 18 Aug 2015 09:21:29 +0100 Subject: [PATCH 05/11] Small dh obtained in two flavour ratio so looks ok. --- lib/Make.inc | 2 +- lib/qcd/action/Actions.h | 1 + lib/qcd/action/pseudofermion/TwoFlavour.h | 9 --- .../action/pseudofermion/TwoFlavourRatio.h | 37 ++++++------ tests/Make.inc | 6 +- tests/Test_hmc_WilsonRatio.cc | 56 +++++++++++++++++++ 6 files changed, 83 insertions(+), 28 deletions(-) create mode 100644 tests/Test_hmc_WilsonRatio.cc diff --git a/lib/Make.inc b/lib/Make.inc index 8da8d28e..f3f2edbf 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/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 +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 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 cb8ef664..331400a4 100644 --- a/lib/qcd/action/Actions.h +++ b/lib/qcd/action/Actions.h @@ -152,5 +152,6 @@ typedef DomainWallFermion GparityDomainWallFermionD; //////////////////////////////////////// #include #include +#include #endif diff --git a/lib/qcd/action/pseudofermion/TwoFlavour.h b/lib/qcd/action/pseudofermion/TwoFlavour.h index c04231a7..98714c98 100644 --- a/lib/qcd/action/pseudofermion/TwoFlavour.h +++ b/lib/qcd/action/pseudofermion/TwoFlavour.h @@ -4,15 +4,6 @@ namespace Grid{ namespace QCD{ - // Placeholder comments: - - /////////////////////////////////////// - // Two flavour ratio - /////////////////////////////////////// - // S = phi^dag V (Mdag M)^-1 V^dag phi - // 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 /////////////////////////////////////// // One flavour rational diff --git a/lib/qcd/action/pseudofermion/TwoFlavourRatio.h b/lib/qcd/action/pseudofermion/TwoFlavourRatio.h index 1098c85a..7279a051 100644 --- a/lib/qcd/action/pseudofermion/TwoFlavourRatio.h +++ b/lib/qcd/action/pseudofermion/TwoFlavourRatio.h @@ -1,5 +1,5 @@ -#ifndef QCD_PSEUDOFERMION_TWO_FLAVOUR_H -#define QCD_PSEUDOFERMION_TWO_FLAVOUR_H +#ifndef QCD_PSEUDOFERMION_TWO_FLAVOUR_RATIO_H +#define QCD_PSEUDOFERMION_TWO_FLAVOUR_RATIO_H namespace Grid{ namespace QCD{ @@ -26,30 +26,35 @@ namespace Grid{ FermionOperator &_DenOp, OperatorFunction & DS, OperatorFunction & AS - ) : NumOp(_NumOp), DenOp(_DenOp), DerivativeSolver(DS), ActionSolver(AS), Phi(Op.FermionGrid()) {}; + ) : NumOp(_NumOp), DenOp(_DenOp), DerivativeSolver(DS), ActionSolver(AS), Phi(_NumOp.FermionGrid()) {}; virtual void init(const GaugeField &U, GridParallelRNG& pRNG) { // P(phi) = e^{- phi^dag V (MdagM)^-1 Vdag phi} // - // phi = Vdag^{-1} Mdag eta + // NumOp == V + // DenOp == M + // + // Take phi = Vdag^{-1} Mdag eta ; eta = Mdag^{-1} Vdag Phi // // P(eta) = e^{- eta^dag eta} // // e^{x^2/2 sig^2} => sig^2 = 0.5. // - // So eta should be of width sig = 1/sqrt(2). - // and must multiply by 0.707.... + // So eta should be of width sig = 1/sqrt(2) and must multiply by 0.707.... // RealD scale = std::sqrt(0.5); FermionField eta(NumOp.FermionGrid()); + FermionField tmp(NumOp.FermionGrid()); gaussian(pRNG,eta); NumOp.ImportGauge(U); DenOp.ImportGauge(U); + // Note: this hard codes normal equations type solvers; alternate implementation needed for + // non-herm style solvers. MdagMLinearOperator ,FermionField> MdagMOp(NumOp); DenOp.Mdag(eta,Phi); // Mdag eta @@ -74,9 +79,9 @@ namespace Grid{ MdagMLinearOperator ,FermionField> MdagMOp(DenOp); X=zero; - NumOp.Mdag(Phi,Y); // Vdag phi - ActionSolver(MdagMOp,Y,X); // MdagMinv Vdag phi - MdagMOp.Op(X,Y); // Y=Mdaginv Vdag phi + NumOp.Mdag(Phi,Y); // Y= Vdag phi + ActionSolver(MdagMOp,Y,X); // X= (MdagM)^-1 Vdag phi + DenOp.M(X,Y); // Y= Mdag^-1 Vdag phi RealD action = norm2(Y); @@ -99,18 +104,16 @@ namespace Grid{ FermionField Y(NumOp.FermionGrid()); FermionField f1(NumOp.FermionGrid()); - GaugeField force(FermOp.GaugeGrid()); + GaugeField force(NumOp.GaugeGrid()); X=zero; - //f1=Vdag phi - NumOp.Mdag(phi,f1); - + //Y=Vdag phi //X = (Mdag M)^-1 V^dag phi - DerivativeSolver(MdagMOp,f1,X); - //Y = (Mdag)^-1 V^dag phi - DenOp.M(X,Y); + NumOp.Mdag(Phi,Y); // Y= Vdag phi + DerivativeSolver(MdagMOp,Y,X); // X= (MdagM)^-1 Vdag phi + DenOp.M(X,Y); // Y= Mdag^-1 Vdag phi // phi^dag V (Mdag M)^-1 dV^dag phi NumOp.MDeriv(force , X, Phi, DaggerYes ); dSdU=force; @@ -122,7 +125,7 @@ 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); }; diff --git a/tests/Make.inc b/tests/Make.inc index c765cdf7..bfacb980 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_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_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 @@ -110,6 +110,10 @@ Test_hmc_WilsonGauge_SOURCES=Test_hmc_WilsonGauge.cc Test_hmc_WilsonGauge_LDADD=-lGrid +Test_hmc_WilsonRatio_SOURCES=Test_hmc_WilsonRatio.cc +Test_hmc_WilsonRatio_LDADD=-lGrid + + Test_lie_generators_SOURCES=Test_lie_generators.cc Test_lie_generators_LDADD=-lGrid diff --git a/tests/Test_hmc_WilsonRatio.cc b/tests/Test_hmc_WilsonRatio.cc new file mode 100644 index 00000000..c7cdec39 --- /dev/null +++ b/tests/Test_hmc_WilsonRatio.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 FermOp(U,Fine,RBFine,mass); + WilsonFermionR NumOp(U,Fine,RBFine,pv); + + ConjugateGradient CG(1.0e-8,10000); + + TwoFlavourRatioPseudoFermionAction WilsonNf2(NumOp,FermOp,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); + +} From bdcbfe9310eef5cdb18f48d57e4bbd015e9d7557 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Tue, 18 Aug 2015 10:37:08 +0100 Subject: [PATCH 06/11] Even Odd two flavour ratio added and dH == small --- lib/Make.inc | 2 +- lib/qcd/action/Actions.h | 10 +- .../EvenOddSchurDifferentiable.h | 112 ++++++++++++ .../action/pseudofermion/TwoFlavourEvenOdd.h | 118 +------------ .../pseudofermion/TwoFlavourEvenOddRatio.h | 163 ++++++++++++++++++ .../action/pseudofermion/TwoFlavourRatio.h | 10 +- tests/Make.inc | 6 +- tests/Test_hmc_EOWilsonRatio.cc | 56 ++++++ 8 files changed, 356 insertions(+), 121 deletions(-) create mode 100644 lib/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h create mode 100644 lib/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h create mode 100644 tests/Test_hmc_EOWilsonRatio.cc 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); + +} From 2dd9ad7b0f7fd3e46a2baff46ad9b41167ab4f23 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Tue, 18 Aug 2015 10:43:32 +0100 Subject: [PATCH 07/11] Update TODO list --- TODO | 55 ++++++++++++++++++++++++++++++------------------------- 1 file changed, 30 insertions(+), 25 deletions(-) diff --git a/TODO b/TODO index 47067c70..432b9028 100644 --- a/TODO +++ b/TODO @@ -1,25 +1,30 @@ - PseudoFermions -TwoFlavour -- done. Fill in rest: - - TwoFlavourEvenOddConstEE - - TwoFlavourEvenOdd - - TwoFlavourRatioEvenOddConstEE - - TwoFlavourRatioEvenOdd +Done: Cayley, Partial , ContFrac force terms. + +Done: + - TwoFlavour + - TwoFlavourEvenOdd - TwoFlavourRatio + - TwoFlavourRatioEvenOdd + => generalise to non-const EE +fill in: - OneFlavourRationalEvenOddConstEE - OneFlavourRationalEvenOdd - OneFlavourRationalRatioEvenOddConstEE - OneFlavourRationalRatioEvenOdd - OneFlavourRationalRatio +- Clean up HMC +- Integrators + - Force Gradient + - Multi-timescale looks broken and operating on single timescale for now. + Fix/debug/rewrite this + - MacroMagic -> readers -- Link smearing/boundary conds - -- Clean up HMC - -- Test Cayley, Partial , ContFrac force terms. +- Link smearing/boundary conds; Policy class based implementation - Test DWF HMC @@ -38,13 +43,16 @@ TwoFlavour -- done. Fill in rest: ================================================================ * Extract/merge/set cleanup ; too many variants; rationalise and call simpler ones -* Used #define repetitive sequences to minimise code. + * Rewrite core tensor arithmetic support to be more systematic + = Use #define repetitive sequences to minimise code, decrease line count by thousands possible, + with more robust and maintainable implementation. + * Ensure we ET as much as possible; move unop functions into ET framework. - tests with expression args to all functions - * FIXME audit + * const audit Insert/Extract @@ -57,10 +65,10 @@ Insert/Extract * Thread scaling tests Xeon, XeonPhi -** Make the Tensor types and Complex etc... play more nicely. +* Make the Tensor types and Complex etc... play more nicely. - TensorRemove is a hack, come up with a long term rationalised approach to Complex vs. Scalar > > - QDP forces use of "toDouble" to get back to non tensor scalar. This role is presently taken TensorRemove, but I - want to introduce a syntax that does not require this. + QDP forces use of "toDouble" to get back to non tensor scalar. This role is presently taken TensorRemove, but I + want to introduce a syntax that does not require this. - Reductions that contract indices on a site should always demote the tensor structure. norm2(), innerProduct. @@ -93,10 +101,6 @@ Insert/Extract *** New Functionality ================================================================ -* - BinaryWriter, TextWriter etc... - - use protocol buffers? replace xmlReader/Writer ec.. - - Binary use htonll, htonl - * CovariantShift support -----Use a class to store gauge field? (parallel transport?) * Parallel io improvements @@ -138,14 +142,15 @@ Algorithms (lots of reuse/port from BFM) * Gauge - Wilson, symanzik, iwasaki -* rb4d support for 5th dimension in Mobius. - * Flavour matrices? -* Pauli, SU subgroup, etc.. -* su3 exponentiation & log etc.. [Jamie's code?] -* TaProj -* FFTnD ? +* Pauli, SU subgroup, etc.. + +* su3 exponentiation & log etc.. [Jamie's code?] + +* TaProj + +* FFTnD ? ====================================================================================================== FUNCTIONALITY: it pleases me to keep track of things I have done (keeps me arguably sane) From 5c364f80823fe58d2d9edb78c2cb7a9a0776adc6 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Tue, 18 Aug 2015 14:40:08 +0100 Subject: [PATCH 08/11] One flavour rational unprec added; untested but does compile. Moving param structs into a single header for later connection to file I/O using macromagic.h --- TODO | 5 +- lib/algorithms/approx/MultiShiftFunction.h | 19 +- lib/qcd/action/ActionParams.h | 26 +++ lib/qcd/action/Actions.h | 9 +- lib/qcd/action/fermion/FermionOperator.h | 4 +- lib/qcd/action/fermion/FermionOperatorImpl.h | 7 +- .../action/pseudofermion/OneFlavourRational.h | 170 ++++++++++++++++++ lib/qcd/action/pseudofermion/TwoFlavour.h | 79 -------- .../action/pseudofermion/TwoFlavourEvenOdd.h | 7 +- .../pseudofermion/TwoFlavourEvenOddRatio.h | 11 +- 10 files changed, 238 insertions(+), 99 deletions(-) create mode 100644 lib/qcd/action/ActionParams.h create mode 100644 lib/qcd/action/pseudofermion/OneFlavourRational.h diff --git a/TODO b/TODO index 432b9028..dc7065d1 100644 --- a/TODO +++ b/TODO @@ -21,8 +21,11 @@ fill in: - Force Gradient - Multi-timescale looks broken and operating on single timescale for now. Fix/debug/rewrite this + - Sign of force term. + - Prefer "RefreshInternal" or such like to "init" in naming + - Rename "Ta" as too unclear -- MacroMagic -> readers +- MacroMagic -> virtual reader class. - Link smearing/boundary conds; Policy class based implementation diff --git a/lib/algorithms/approx/MultiShiftFunction.h b/lib/algorithms/approx/MultiShiftFunction.h index fceed2eb..c9ed39e3 100644 --- a/lib/algorithms/approx/MultiShiftFunction.h +++ b/lib/algorithms/approx/MultiShiftFunction.h @@ -11,20 +11,27 @@ public: std::vector tolerances; RealD norm; RealD lo,hi; + MultiShiftFunction(int n,RealD _lo,RealD _hi): poles(n), residues(n), lo(_lo), hi(_hi) {;}; RealD approx(RealD x); void csv(std::ostream &out); void gnuplot(std::ostream &out); - MultiShiftFunction(AlgRemez & remez,double tol,bool inverse) : - order(remez.getDegree()), - tolerances(remez.getDegree(),tol), - poles(remez.getDegree()), - residues(remez.getDegree()) + + void Init(AlgRemez & remez,double tol,bool inverse) { + order=remez.getDegree(); + tolerances.resize(remez.getDegree(),tol); + poles.resize(remez.getDegree()); + residues.resize(remez.getDegree()); remez.getBounds(lo,hi); if ( inverse ) remez.getIPFE (&residues[0],&poles[0],&norm); - else remez.getPFE (&residues[0],&poles[0],&norm); + else remez.getPFE (&residues[0],&poles[0],&norm); } + MultiShiftFunction(AlgRemez & remez,double tol,bool inverse) + { + Init(remez,tol,inverse); + } + }; } #endif diff --git a/lib/qcd/action/ActionParams.h b/lib/qcd/action/ActionParams.h new file mode 100644 index 00000000..d866339c --- /dev/null +++ b/lib/qcd/action/ActionParams.h @@ -0,0 +1,26 @@ +#ifndef GRID_QCD_ACTION_PARAMS_H +#define GRID_QCD_ACTION_PARAMS_H + +namespace Grid { +namespace QCD { + + // These can move into a params header and be given MacroMagic serialisation + struct GparityWilsonImplParams { + std::vector twists; + }; + + struct WilsonImplParams { }; + + struct OneFlavourRationalParams { + RealD lo; + RealD hi; + int precision=64; + int degree=10; + RealD tolerance; // Vector? + RealD MaxIter; // Vector? + OneFlavourRationalParams (RealD lo,RealD hi,int precision=64,int degree = 10); + }; + +}} + +#endif diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h index bc9526ca..3e62a509 100644 --- a/lib/qcd/action/Actions.h +++ b/lib/qcd/action/Actions.h @@ -14,6 +14,7 @@ // Abstract base interface //////////////////////////////////////////// #include +#include //////////////////////////////////////////// // Gauge Actions @@ -157,9 +158,9 @@ typedef DomainWallFermion GparityDomainWallFermionD; #include //Todo: RHMC -//#include -//#include -//#include -//#include +#include +//#include +//#include +//#include #endif diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h index e728909e..7cdfc0d1 100644 --- a/lib/qcd/action/fermion/FermionOperator.h +++ b/lib/qcd/action/fermion/FermionOperator.h @@ -32,6 +32,8 @@ namespace Grid { virtual RealD Mdag (const FermionField &in, FermionField &out)=0; // half checkerboard operaions + virtual int ConstEE(void) { return 1; }; // clover returns zero as EE depends on gauge field + virtual void Meooe (const FermionField &in, FermionField &out)=0; virtual void MeooeDag (const FermionField &in, FermionField &out)=0; virtual void Mooee (const FermionField &in, FermionField &out)=0; @@ -49,7 +51,7 @@ namespace Grid { virtual void MDeriv (GaugeField &mat,const FermionField &U,const FermionField &V,int dag){DhopDeriv(mat,U,V,dag);}; virtual void MoeDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag){DhopDerivOE(mat,U,V,dag);}; virtual void MeoDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag){DhopDerivEO(mat,U,V,dag);}; - virtual void MooDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag){mat=zero;}; + virtual void MooDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag){mat=zero;}; // Clover can override these virtual void MeeDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag){mat=zero;}; virtual void DhopDeriv (GaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0; diff --git a/lib/qcd/action/fermion/FermionOperatorImpl.h b/lib/qcd/action/fermion/FermionOperatorImpl.h index 24c64149..1d90cace 100644 --- a/lib/qcd/action/fermion/FermionOperatorImpl.h +++ b/lib/qcd/action/fermion/FermionOperatorImpl.h @@ -5,6 +5,7 @@ namespace Grid { namespace QCD { + ////////////////////////////////////////////// // Template parameter class constructs to package // externally control Fermion implementations @@ -126,8 +127,7 @@ namespace Grid { typedef Lattice DoubledGaugeField; typedef WilsonCompressor Compressor; - - typedef struct WilsonImplParams { } ImplParams; + typedef WilsonImplParams ImplParams; ImplParams Params; WilsonImpl(const ImplParams &p= ImplParams()) : Params(p) {}; @@ -177,6 +177,7 @@ PARALLEL_FOR_LOOP //////////////////////////////////////////////////////////////////////////////////////// // Flavour doubled spinors; is Gparity the only? what about C*? //////////////////////////////////////////////////////////////////////////////////////// + template class GparityWilsonImpl : public ImplGauge { public: @@ -198,7 +199,7 @@ PARALLEL_FOR_LOOP typedef WilsonCompressor Compressor; - typedef struct GparityWilsonImplParams {std::vector twists; } ImplParams; + typedef GparityWilsonImplParams ImplParams; ImplParams Params; GparityWilsonImpl(const ImplParams &p= ImplParams()) : Params(p) {}; diff --git a/lib/qcd/action/pseudofermion/OneFlavourRational.h b/lib/qcd/action/pseudofermion/OneFlavourRational.h new file mode 100644 index 00000000..e7596065 --- /dev/null +++ b/lib/qcd/action/pseudofermion/OneFlavourRational.h @@ -0,0 +1,170 @@ +#ifndef QCD_PSEUDOFERMION_ONE_FLAVOUR_RATIONAL_H +#define QCD_PSEUDOFERMION_ONE_FLAVOUR_RATIONAL_H + +namespace Grid{ + namespace QCD{ + + /////////////////////////////////////// + // One flavour rational + /////////////////////////////////////// + + // S_f = chi^dag * N(M^dag*M)/D(M^dag*M) * chi + // + // Here, M is some operator + // N and D makeup the rat. poly + // + + template + class OneFlavourRationalPseudoFermionAction : public Action { + public: + INHERIT_IMPL_TYPES(Impl); + + typedef OneFlavourRationalParams Params; + Params param; + + MultiShiftFunction PowerHalf ; + MultiShiftFunction PowerNegHalf; + MultiShiftFunction PowerQuarter; + MultiShiftFunction PowerNegQuarter; + + private: + + FermionOperator & FermOp;// the basic operator + + // NOT using "Nroots"; IroIro is -- perhaps later, but this wasn't good for us historically + // and hasenbusch works better + + FermionField Phi; // the pseudo fermion field for this trajectory + + public: + + OneFlavourRationalPseudoFermionAction(FermionOperator &Op, + Params & p + ) : FermOp(Op), Phi(Op.FermionGrid()), param(p) + { + AlgRemez remez(param.lo,param.hi,param.precision); + + // MdagM^(+- 1/2) + std::cout< sig^2 = 0.5. + // + // So eta should be of width sig = 1/sqrt(2). + + RealD scale = std::sqrt(0.5); + + FermionField eta(FermOp.FermionGrid()); + + gaussian(pRNG,eta); + + FermOp.ImportGauge(U); + + // mutishift CG + MdagMLinearOperator ,FermionField> MdagMOp(FermOp); + ConjugateGradientMultiShift msCG(param.MaxIter,PowerQuarter); + msCG(MdagMOp,eta,Phi); + + Phi=Phi*scale; + + }; + + ////////////////////////////////////////////////////// + // S = phi^dag (Mdag M)^-1/2 phi + ////////////////////////////////////////////////////// + virtual RealD S(const GaugeField &U) { + + FermOp.ImportGauge(U); + + FermionField Y(FermOp.FermionGrid()); + + MdagMLinearOperator ,FermionField> MdagMOp(FermOp); + + ConjugateGradientMultiShift msCG(param.MaxIter,PowerNegQuarter); + + msCG(MdagMOp,Phi,Y); + + RealD action = norm2(Y); + std::cout << GridLogMessage << "Pseudofermion action FIXME -- is -1/4 solve or -1/2 solve faster??? "< MPhi_k (Npole,FermOp.FermionGrid()); + + FermionField X(FermOp.FermionGrid()); + FermionField Y(FermOp.FermionGrid()); + + GaugeField tmp(FermOp.GaugeGrid()); + + FermOp.ImportGauge(U); + + MdagMLinearOperator ,FermionField> MdagMOp(FermOp); + + ConjugateGradientMultiShift msCG(param.MaxIter,PowerNegHalf); + + msCG(MdagMOp,Phi,MPhi_k); + + dSdU = zero; + for(int k=0;k Date: Tue, 18 Aug 2015 16:58:56 +0100 Subject: [PATCH 09/11] rhmc for 1+1 wilson is conserving dH~0. A good days work ;) --- tests/Test_rhmc_Wilson1p1.cc | 56 ++++++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) create mode 100644 tests/Test_rhmc_Wilson1p1.cc diff --git a/tests/Test_rhmc_Wilson1p1.cc b/tests/Test_rhmc_Wilson1p1.cc new file mode 100644 index 00000000..f44e121a --- /dev/null +++ b/tests/Test_rhmc_Wilson1p1.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; + WilsonFermionR FermOp(U,Fine,RBFine,mass); + + // 1+1 flavour + OneFlavourRationalParams Params(1.0e-4,64.0,1000,1.0e-6); + OneFlavourRationalPseudoFermionAction WilsonNf1a(FermOp,Params); + OneFlavourRationalPseudoFermionAction WilsonNf1b(FermOp,Params); + + //Collect actions + ActionLevel Level1; + Level1.push_back(&WilsonNf1a); + Level1.push_back(&WilsonNf1b); + Level1.push_back(&Waction); + + ActionSet FullSet; + FullSet.push_back(Level1); + + // Create integrator + typedef MinimumNorm2 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); + +} From 48db72259e66d6486d35e1e7ff2df6da74057648 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Tue, 18 Aug 2015 18:37:39 +0100 Subject: [PATCH 10/11] EvenOdd schur decomposed mpcdagmpc version of rhmc determinant. dH is also small and plaquette looks right. --- .../pseudofermion/OneFlavourEvenOddRational.h | 185 ++++++++++++++++++ tests/Test_rhmc_EOWilson1p1.cc | 55 ++++++ 2 files changed, 240 insertions(+) create mode 100644 lib/qcd/action/pseudofermion/OneFlavourEvenOddRational.h create mode 100644 tests/Test_rhmc_EOWilson1p1.cc diff --git a/lib/qcd/action/pseudofermion/OneFlavourEvenOddRational.h b/lib/qcd/action/pseudofermion/OneFlavourEvenOddRational.h new file mode 100644 index 00000000..beb3abc1 --- /dev/null +++ b/lib/qcd/action/pseudofermion/OneFlavourEvenOddRational.h @@ -0,0 +1,185 @@ +#ifndef QCD_PSEUDOFERMION_ONE_FLAVOUR_EVEN_ODD_RATIONAL_H +#define QCD_PSEUDOFERMION_ONE_FLAVOUR_EVEN_ODD_RATIONAL_H + +namespace Grid{ + namespace QCD{ + + /////////////////////////////////////// + // One flavour rational + /////////////////////////////////////// + + // S_f = chi^dag * N(Mpc^dag*Mpc)/D(Mpc^dag*Mpc) * chi + // + // Here, M is some operator + // N and D makeup the rat. poly + // + + template + class OneFlavourEvenOddRationalPseudoFermionAction : public Action { + public: + INHERIT_IMPL_TYPES(Impl); + + typedef OneFlavourRationalParams Params; + Params param; + + MultiShiftFunction PowerHalf ; + MultiShiftFunction PowerNegHalf; + MultiShiftFunction PowerQuarter; + MultiShiftFunction PowerNegQuarter; + + private: + + FermionOperator & FermOp;// the basic operator + + // NOT using "Nroots"; IroIro is -- perhaps later, but this wasn't good for us historically + // and hasenbusch works better + + FermionField PhiEven; // the pseudo fermion field for this trajectory + FermionField PhiOdd; // the pseudo fermion field for this trajectory + + + public: + + OneFlavourEvenOddRationalPseudoFermionAction(FermionOperator &Op, + Params & p ) : FermOp(Op), + PhiEven(Op.FermionRedBlackGrid()), + PhiOdd (Op.FermionRedBlackGrid()), + param(p) + { + AlgRemez remez(param.lo,param.hi,param.precision); + + // MdagM^(+- 1/2) + std::cout< sig^2 = 0.5. + // + // So eta should be of width sig = 1/sqrt(2). + + RealD scale = std::sqrt(0.5); + + FermionField eta (FermOp.FermionGrid()); + FermionField etaOdd (FermOp.FermionRedBlackGrid()); + FermionField etaEven(FermOp.FermionRedBlackGrid()); + + gaussian(pRNG,eta); eta=eta*scale; + + pickCheckerboard(Even,etaEven,eta); + pickCheckerboard(Odd,etaOdd,eta); + + FermOp.ImportGauge(U); + + // mutishift CG + SchurDifferentiableOperator Mpc(FermOp); + ConjugateGradientMultiShift msCG(param.MaxIter,PowerQuarter); + msCG(Mpc,etaOdd,PhiOdd); + + ////////////////////////////////////////////////////// + // FIXME : Clover term not yet.. + ////////////////////////////////////////////////////// + + assert(FermOp.ConstEE() == 1); + PhiEven = zero; + + }; + + ////////////////////////////////////////////////////// + // S = phi^dag (Mdag M)^-1/2 phi + ////////////////////////////////////////////////////// + virtual RealD S(const GaugeField &U) { + + FermOp.ImportGauge(U); + + FermionField Y(FermOp.FermionRedBlackGrid()); + + SchurDifferentiableOperator Mpc(FermOp); + + ConjugateGradientMultiShift msCG(param.MaxIter,PowerNegQuarter); + + msCG(Mpc,PhiOdd,Y); + + RealD action = norm2(Y); + std::cout << GridLogMessage << "Pseudofermion action FIXME -- is -1/4 solve or -1/2 solve faster??? "< MPhi_k (Npole,FermOp.FermionRedBlackGrid()); + + FermionField X(FermOp.FermionRedBlackGrid()); + FermionField Y(FermOp.FermionRedBlackGrid()); + + GaugeField tmp(FermOp.GaugeGrid()); + + FermOp.ImportGauge(U); + + SchurDifferentiableOperator Mpc(FermOp); + + ConjugateGradientMultiShift msCG(param.MaxIter,PowerNegHalf); + + msCG(Mpc,PhiOdd,MPhi_k); + + dSdU = zero; + for(int k=0;k 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; + WilsonFermionR FermOp(U,Fine,RBFine,mass); + + // 1+1 flavour + OneFlavourRationalParams Params(1.0e-4,64.0,1000,1.0e-6); + OneFlavourEvenOddRationalPseudoFermionAction WilsonNf1a(FermOp,Params); + OneFlavourEvenOddRationalPseudoFermionAction WilsonNf1b(FermOp,Params); + + //Collect actions + ActionLevel Level1; + Level1.push_back(&WilsonNf1a); + Level1.push_back(&WilsonNf1b); + Level1.push_back(&Waction); + + ActionSet FullSet; + FullSet.push_back(Level1); + + // Create integrator + typedef MinimumNorm2 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); + +} From dd6bb73ee01b721a94b67e8570b77ec179cd391e Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Wed, 19 Aug 2015 04:58:40 +0100 Subject: [PATCH 11/11] Added one flavour rational ratios (unprec) --- lib/Make.inc | 2 +- lib/algorithms/approx/MultiShiftFunction.h | 2 + .../iterative/ConjugateGradientMultiShift.h | 6 +- lib/qcd/action/ActionParams.h | 11 +- lib/qcd/action/Actions.h | 4 +- .../pseudofermion/OneFlavourRationalRatio.h | 226 ++++++++++++++++++ .../action/pseudofermion/TwoFlavourEvenOdd.h | 2 +- tests/Make.inc | 14 +- tests/Test_rhmc_WilsonRatio.cc | 58 +++++ 9 files changed, 315 insertions(+), 10 deletions(-) create mode 100644 lib/qcd/action/pseudofermion/OneFlavourRationalRatio.h create mode 100644 tests/Test_rhmc_WilsonRatio.cc diff --git a/lib/Make.inc b/lib/Make.inc index dd9dce12..d2320011 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/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 +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/ActionParams.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/OneFlavourEvenOddRational.h ./qcd/action/pseudofermion/OneFlavourRational.h ./qcd/action/pseudofermion/OneFlavourRationalRatio.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/algorithms/approx/MultiShiftFunction.h b/lib/algorithms/approx/MultiShiftFunction.h index c9ed39e3..1bd6faae 100644 --- a/lib/algorithms/approx/MultiShiftFunction.h +++ b/lib/algorithms/approx/MultiShiftFunction.h @@ -27,6 +27,8 @@ public: if ( inverse ) remez.getIPFE (&residues[0],&poles[0],&norm); else remez.getPFE (&residues[0],&poles[0],&norm); } + // Allow deferred initialisation + MultiShiftFunction(void){}; MultiShiftFunction(AlgRemez & remez,double tol,bool inverse) { Init(remez,tol,inverse); diff --git a/lib/algorithms/iterative/ConjugateGradientMultiShift.h b/lib/algorithms/iterative/ConjugateGradientMultiShift.h index c6e7650d..14177abc 100644 --- a/lib/algorithms/iterative/ConjugateGradientMultiShift.h +++ b/lib/algorithms/iterative/ConjugateGradientMultiShift.h @@ -27,10 +27,14 @@ public: void operator() (LinearOperatorBase &Linop, const Field &src, Field &psi) { - GridBase *grid = src._grid; int nshift = shifts.order; std::vector results(nshift,grid); + (*this)(Linop,src,results,psi); +} +void operator() (LinearOperatorBase &Linop, const Field &src, std::vector &results, Field &psi) +{ + int nshift = shifts.order; (*this)(Linop,src,results); diff --git a/lib/qcd/action/ActionParams.h b/lib/qcd/action/ActionParams.h index d866339c..55481bc2 100644 --- a/lib/qcd/action/ActionParams.h +++ b/lib/qcd/action/ActionParams.h @@ -14,11 +14,14 @@ namespace QCD { struct OneFlavourRationalParams { RealD lo; RealD hi; - int precision=64; - int degree=10; + int MaxIter; // Vector? RealD tolerance; // Vector? - RealD MaxIter; // Vector? - OneFlavourRationalParams (RealD lo,RealD hi,int precision=64,int degree = 10); + int degree=10; + int precision=64; + + OneFlavourRationalParams (RealD _lo,RealD _hi,int _maxit,RealD tol=1.0e-8,int _degree = 10,int _precision=64) : + lo(_lo), hi(_hi), MaxIter(_maxit), tolerance(tol), degree(_degree), precision(_precision) + {}; }; }} diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h index 3e62a509..fa6c3eaf 100644 --- a/lib/qcd/action/Actions.h +++ b/lib/qcd/action/Actions.h @@ -159,8 +159,8 @@ typedef DomainWallFermion GparityDomainWallFermionD; //Todo: RHMC #include -//#include -//#include +#include +#include //#include #endif diff --git a/lib/qcd/action/pseudofermion/OneFlavourRationalRatio.h b/lib/qcd/action/pseudofermion/OneFlavourRationalRatio.h new file mode 100644 index 00000000..1144af17 --- /dev/null +++ b/lib/qcd/action/pseudofermion/OneFlavourRationalRatio.h @@ -0,0 +1,226 @@ +#ifndef QCD_PSEUDOFERMION_ONE_FLAVOUR_RATIONAL_RATIO_H +#define QCD_PSEUDOFERMION_ONE_FLAVOUR_RATIONAL_RATIO_H + +namespace Grid{ + namespace QCD{ + + /////////////////////////////////////// + // One flavour rational + /////////////////////////////////////// + + // S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi + // + // Here P/Q \sim R_{1/4} ~ (V^dagV)^{1/4} + // Here N/D \sim R_{-1/2} ~ (M^dagM)^{-1/2} + + template + class OneFlavourRatioRationalPseudoFermionAction : public Action { + public: + + INHERIT_IMPL_TYPES(Impl); + + typedef OneFlavourRationalParams Params; + Params param; + + MultiShiftFunction PowerHalf ; + MultiShiftFunction PowerNegHalf; + MultiShiftFunction PowerQuarter; + MultiShiftFunction PowerNegQuarter; + + private: + + FermionOperator & NumOp;// the basic operator + FermionOperator & DenOp;// the basic operator + FermionField Phi; // the pseudo fermion field for this trajectory + + public: + + OneFlavourRatioRationalPseudoFermionAction(FermionOperator &_NumOp, + FermionOperator &_DenOp, + Params & p + ) : NumOp(_NumOp), DenOp(_DenOp), Phi(_NumOp.FermionGrid()), param(p) + { + AlgRemez remez(param.lo,param.hi,param.precision); + + // MdagM^(+- 1/2) + std::cout< sig^2 = 0.5. + // + // So eta should be of width sig = 1/sqrt(2). + + RealD scale = std::sqrt(0.5); + + FermionField tmp(NumOp.FermionGrid()); + FermionField eta(NumOp.FermionGrid()); + + gaussian(pRNG,eta); + + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + // MdagM^1/4 eta + MdagMLinearOperator ,FermionField> MdagM(DenOp); + ConjugateGradientMultiShift msCG_M(param.MaxIter,PowerQuarter); + msCG_M(MdagM,eta,tmp); + + // VdagV^-1/4 MdagM^1/4 eta + MdagMLinearOperator ,FermionField> VdagV(NumOp); + ConjugateGradientMultiShift msCG_V(param.MaxIter,PowerNegQuarter); + msCG_V(VdagV,tmp,Phi); + + Phi=Phi*scale; + + }; + + ////////////////////////////////////////////////////// + // S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi + ////////////////////////////////////////////////////// + virtual RealD S(const GaugeField &U) { + + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + FermionField X(NumOp.FermionGrid()); + FermionField Y(NumOp.FermionGrid()); + + // VdagV^1/4 Phi + MdagMLinearOperator ,FermionField> VdagV(NumOp); + ConjugateGradientMultiShift msCG_V(param.MaxIter,PowerQuarter); + msCG_V(VdagV,Phi,X); + + // MdagM^-1/4 VdagV^1/4 Phi + MdagMLinearOperator ,FermionField> MdagM(DenOp); + ConjugateGradientMultiShift msCG_M(param.MaxIter,PowerNegQuarter); + msCG_M(MdagM,X,Y); + + // Phidag VdagV^1/4 MdagM^-1/4 MdagM^-1/4 VdagV^1/4 Phi + RealD action = norm2(Y); + + return action; + }; + + // S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi + // + // Here, M is some 5D operator and V is the Pauli-Villars field + // N and D makeup the rat. poly of the M term and P and & makeup the rat.poly of the denom term + // + // Need + // dS_f/dU = chi^dag d[P/Q] N/D P/Q chi + // + chi^dag P/Q d[N/D] P/Q chi + // + chi^dag P/Q N/D d[P/Q] chi + // + // P/Q is expressed as partial fraction expansion: + // + // a0 + \sum_k ak/(V^dagV + bk) + // + // d[P/Q] is then + // + // \sum_k -ak [V^dagV+bk]^{-1} [ dV^dag V + V^dag dV ] [V^dag V + bk]^{-1} + // + // and similar for N/D. + // + // Need + // MpvPhi_k = [Vdag V + bk]^{-1} chi + // MpvPhi = {a0 + \sum_k ak [Vdag V + bk]^{-1} }chi + // + // MfMpvPhi_k = [MdagM+bk]^{-1} MpvPhi + // MfMpvPhi = {a0 + \sum_k ak [Mdag M + bk]^{-1} } MpvPhi + // + // MpvMfMpvPhi_k = [Vdag V + bk]^{-1} MfMpvchi + // + + virtual void deriv(const GaugeField &U,GaugeField & dSdU) { + + const int n_f = PowerNegHalf.poles.size(); + const int n_pv = PowerQuarter.poles.size(); + + std::vector MpvPhi_k (n_pv,NumOp.FermionGrid()); + std::vector MpvMfMpvPhi_k(n_pv,NumOp.FermionGrid()); + std::vector MfMpvPhi_k (n_f,NumOp.FermionGrid()); + + FermionField MpvPhi(NumOp.FermionGrid()); + FermionField MfMpvPhi(NumOp.FermionGrid()); + FermionField MpvMfMpvPhi(NumOp.FermionGrid()); + FermionField Y(NumOp.FermionGrid()); + + GaugeField tmp(NumOp.GaugeGrid()); + + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + MdagMLinearOperator ,FermionField> MdagM(DenOp); + MdagMLinearOperator ,FermionField> VdagV(NumOp); + + ConjugateGradientMultiShift msCG_V(param.MaxIter,PowerQuarter); + ConjugateGradientMultiShift msCG_M(param.MaxIter,PowerNegHalf); + + msCG_V(VdagV,Phi,MpvPhi_k,MpvPhi); + msCG_M(MdagM,MpvPhi,MfMpvPhi_k,MfMpvPhi); + msCG_V(VdagV,MfMpvPhi,MpvMfMpvPhi_k,MpvMfMpvPhi); + + RealD ak; + + dSdU = zero; + + // With these building blocks + // + // dS/dU = + // \sum_k -ak MfMpvPhi_k^dag [ dM^dag M + M^dag dM ] MfMpvPhi_k (1) + // + \sum_k -ak MpvMfMpvPhi_k^\dag [ dV^dag V + V^dag dV ] MpvPhi_k (2) + // -ak MpvPhi_k^dag [ dV^dag V + V^dag dV ] MpvMfMpvPhi_k (3) + + //(1) + for(int k=0;k PCop(FermOp); - FermOp.ImportGauge(U); PCop.MpcDag(etaOdd,PhiOdd); diff --git a/tests/Make.inc b/tests/Make.inc index 125778fe..40ae694f 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_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 +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_rhmc_EOWilson1p1 Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio 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 @@ -146,6 +146,18 @@ Test_remez_SOURCES=Test_remez.cc Test_remez_LDADD=-lGrid +Test_rhmc_EOWilson1p1_SOURCES=Test_rhmc_EOWilson1p1.cc +Test_rhmc_EOWilson1p1_LDADD=-lGrid + + +Test_rhmc_Wilson1p1_SOURCES=Test_rhmc_Wilson1p1.cc +Test_rhmc_Wilson1p1_LDADD=-lGrid + + +Test_rhmc_WilsonRatio_SOURCES=Test_rhmc_WilsonRatio.cc +Test_rhmc_WilsonRatio_LDADD=-lGrid + + Test_rng_SOURCES=Test_rng.cc Test_rng_LDADD=-lGrid diff --git a/tests/Test_rhmc_WilsonRatio.cc b/tests/Test_rhmc_WilsonRatio.cc new file mode 100644 index 00000000..9b27fa70 --- /dev/null +++ b/tests/Test_rhmc_WilsonRatio.cc @@ -0,0 +1,58 @@ +#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); + + // erange,maxiter,resid,npoly + OneFlavourRationalParams Params(1.0e-2,64.0,1000,1.0e-6,6); + OneFlavourRatioRationalPseudoFermionAction WilsonNf1a(NumOp,DenOp,Params); + OneFlavourRatioRationalPseudoFermionAction WilsonNf1b(NumOp,DenOp,Params); + + //Collect actions + ActionLevel Level1; + Level1.push_back(&WilsonNf1a); + Level1.push_back(&WilsonNf1b); + 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); + +}