portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-04 03:05:55 +01:00
Code Releases Activity

More info

Browse Source
This commit is contained in:
Peter Boyle 2015-07-21 13:48:57 +09:00
parent 59baa15d9f
commit 01f9b1f6a4
1 changed files with 24 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
lib/algorithms/approx/Chebyshev.h
View File

@ -50,6 +50,15 @@ namespace Grid {
return; 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) ){ Chebyshev(double _lo,double _hi,int _order, double (* func)(double) ){
lo=_lo; lo=_lo;
hi=_hi; hi=_hi;
@ -95,46 +104,39 @@ namespace Grid {
return sum; return sum;
}; };
// Convenience for plotting the approximation // Implement the required interface
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
void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) { void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) {
Field T0 = in; GridBase *grid=in._grid;
Field T1 = T0; // Field T1(T0._grid); more efficient but hardwires Lattice class
Field T2 = T1; 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 *Tnm = &T0;
Field *Tn = &T1; Field *Tn = &T1;
Field *Tnp = &T2; 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 xscale = 2.0/(hi-lo);
double mscale = -(hi+lo)/(hi-lo); double mscale = -(hi+lo)/(hi-lo);
Linop.HermOp(T0,y);
// Tn=T1 = (xscale M + mscale)in
Linop.Op(T0,y);
T1=y*xscale+in*mscale; T1=y*xscale+in*mscale;
// sum = .5 c[0] T0 + c[1] T1 // sum = .5 c[0] T0 + c[1] T1
out = (0.5*Coeffs[0])*T0 + Coeffs[1]*T1; out = (0.5*Coeffs[0])*T0 + Coeffs[1]*T1;
for(int n=2;n<order;n++){ for(int n=2;n<order;n++){
Linop.Op(*Tn,y); Linop.HermOp(*Tn,y);
y=xscale*y+mscale*(*Tn); y=xscale*y+mscale*(*Tn);
*Tnp=2.0*y-(*Tnm); *Tnp=2.0*y-(*Tnm);
out=out+Coeffs[n]* (*Tnp); out=out+Coeffs[n]* (*Tnp);
// Cycle pointers to avoid copies // Cycle pointers to avoid copies