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); + +}