From cbec16feddb60c4359681495c76c57d3e56cecc8 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Tue, 21 Jul 2015 13:48:57 +0900
Subject: [PATCH] More info
---
lib/algorithms/approx/Chebyshev.h | 46 ++++++++++++++++---------------
1 file changed, 24 insertions(+), 22 deletions(-)
diff --git a/lib/algorithms/approx/Chebyshev.h b/lib/algorithms/approx/Chebyshev.h
index 0dad7179..f629ab08 100644
--- a/lib/algorithms/approx/Chebyshev.h
+++ b/lib/algorithms/approx/Chebyshev.h
@@ -50,6 +50,15 @@ namespace Grid {
return;
}
+ // Convenience for plotting the approximation
+ void PlotApprox(std::ostream &out) {
+ out<<"Polynomial approx ["<<lo<<","<<hi<<"]"<<std::endl;
+ for(double x=lo;x<hi;x+=(hi-lo)/50.0){
+ out <<x<<"\t"<<approx(x)<<std::endl;
+ }
+ };
+
+
Chebyshev(double _lo,double _hi,int _order, double (* func)(double) ){
lo=_lo;
hi=_hi;
@@ -95,46 +104,39 @@ namespace Grid {
return sum;
};
- // Convenience for plotting the approximation
- void PlotApprox(std::ostream &out) {
- out<<"Polynomial approx ["<<lo<<","<<hi<<"]"<<std::endl;
- for(double x=lo;x<hi;x+=(hi-lo)/50.0){
- out <<x<<"\t"<<approx(x)<<std::endl;
- }
- };
-
- // Implement the required interface; could require Lattice base class
+ // Implement the required interface
void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) {
- Field T0 = in;
- Field T1 = T0; // Field T1(T0._grid); more efficient but hardwires Lattice class
- Field T2 = T1;
+ GridBase *grid=in._grid;
+
+ int vol=grid->gSites();
+
+ Field T0(grid); T0 = in;
+ Field T1(grid);
+ Field T2(grid);
+ Field y(grid);
- // use a pointer trick to eliminate copies
Field *Tnm = &T0;
Field *Tn = &T1;
Field *Tnp = &T2;
- Field y = in;
-
+
+ std::cout << "Chebyshev ["<<lo<<","<<hi<<"]"<< " order "<<order <<std::endl;
+ // Tn=T1 = (xscale M + mscale)in
double xscale = 2.0/(hi-lo);
double mscale = -(hi+lo)/(hi-lo);
-
- // Tn=T1 = (xscale M + mscale)in
- Linop.Op(T0,y);
-
+ Linop.HermOp(T0,y);
T1=y*xscale+in*mscale;
// sum = .5 c[0] T0 + c[1] T1
out = (0.5*Coeffs[0])*T0 + Coeffs[1]*T1;
-
for(int n=2;n<order;n++){
- Linop.Op(*Tn,y);
+ Linop.HermOp(*Tn,y);
y=xscale*y+mscale*(*Tn);
*Tnp=2.0*y-(*Tnm);
-
+
out=out+Coeffs[n]* (*Tnp);
// Cycle pointers to avoid copies