diff --git a/lib/algorithms/approx/Chebyshev.h b/lib/algorithms/approx/Chebyshev.h index e4fa0b61..37331816 100644 --- a/lib/algorithms/approx/Chebyshev.h +++ b/lib/algorithms/approx/Chebyshev.h @@ -58,7 +58,9 @@ namespace Grid { } }; - + + // c.f. numerical recipes "chebft"/"chebev". This is sec 5.8 "Chebyshev approximation". + // Chebyshev(double _lo,double _hi,int _order, double (* func)(double) ){ lo=_lo; hi=_hi; @@ -77,7 +79,34 @@ namespace Grid { Coeffs[j] = s * 2.0/order; } }; + void JacksonSmooth(void){ + double M=order; + double alpha = M_PI/(M+2); + double lmax = std::cos(alpha); + double sumUsq =0; + std::vector U(M); + std::vector a(M); + std::vector g(M); + for(int n=0;n<=M;n++){ + U[n] = std::sin((n+1)*std::acos(lmax))/std::sin(std::acos(lmax)); + sumUsq += U[n]*U[n]; + } + sumUsq = std::sqrt(sumUsq); + for(int i=1;i<=M;i++){ + a[i] = U[i]/sumUsq; + } + g[0] = 1.0; + for(int m=1;m<=M;m++){ + g[m] = 0; + for(int i=0;i<=M-m;i++){ + g[m]+= a[i]*a[m+i]; + } + } + for(int m=1;m<=M;m++){ + Coeffs[m]*=g[m]; + } + } double approx(double x) // Convenience for plotting the approximation { double Tn; diff --git a/lib/qcd/action/fermion/CayleyFermion5D.cc b/lib/qcd/action/fermion/CayleyFermion5D.cc index 7db67dd7..c18edb7b 100644 --- a/lib/qcd/action/fermion/CayleyFermion5D.cc +++ b/lib/qcd/action/fermion/CayleyFermion5D.cc @@ -17,11 +17,8 @@ namespace QCD { { } - // override multiply - RealD CayleyFermion5D::M (const LatticeFermion &psi, LatticeFermion &chi) + void CayleyFermion5D::Meooe5D (const LatticeFermion &psi, LatticeFermion &Din) { - LatticeFermion Din(psi._grid); - // Assemble Din for(int s=0;s ChebyAccu(0.5,70.0,30,InverseApproximation); Chebyshev Cheby (2.0,70.0,10,InverseApproximation); Chebyshev ChebyAccu(2.0,70.0,10,InverseApproximation); + Cheby.JacksonSmooth(); + ChebyAccu.JacksonSmooth(); _Aggregates.ProjectToSubspace (Csrc,in); _Aggregates.PromoteFromSubspace(Csrc,out); @@ -385,7 +387,7 @@ int main (int argc, char ** argv) LatticeGaugeField Umu(UGrid); NerscField header; - std::string file("./ckpoint_lat.400"); + std::string file("./ckpoint_lat.4000"); readNerscConfiguration(Umu,header,file); // SU3::ColdConfiguration(RNG4,Umu); diff --git a/tests/Test_hmc_WilsonFermionGauge.cc b/tests/Test_hmc_WilsonFermionGauge.cc index f54503af..867a9f96 100644 --- a/tests/Test_hmc_WilsonFermionGauge.cc +++ b/tests/Test_hmc_WilsonFermionGauge.cc @@ -24,13 +24,14 @@ int main (int argc, char ** argv) // simplify template declaration? Strip the lorentz from the second template WilsonGaugeAction Waction(5.6); - Real mass=0.01; + Real mass=-0.77; WilsonFermion FermOp(U,Fine,RBFine,mass); ConjugateGradient CG(1.0e-8,10000); + ConjugateGradient CGmd(1.0e-6,10000); TwoFlavourPseudoFermionAction - Pseudofermion(FermOp,CG,CG,Fine); + Pseudofermion(FermOp,CGmd,CG,Fine); //Collect actions