Poly support for lanczos

lib/algorithms/approx/Chebyshev.h

@@ -9,23 +9,34 @@ namespace Grid {
   ////////////////////////////////////////////////////////////////////////////////////////////
   // Simple general polynomial with user supplied coefficients
   ////////////////////////////////////////////////////////////////////////////////////////////
+  template<class Field>
+  class HermOpOperatorFunction : public OperatorFunction<Field> {
+    void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) {
+      Linop.HermOp(in,out);
+    };
+  };
+
   template<class Field>
   class Polynomial : public OperatorFunction<Field> {
   private:
-    std::vector<double> Coeffs;
+    std::vector<RealD> Coeffs;
   public:
-    Polynomial(std::vector<double> &_Coeffs) : Coeffs(_Coeffs) {};
+    Polynomial(std::vector<RealD> &_Coeffs) : Coeffs(_Coeffs) { };
 
     // Implement the required interface
     void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) {
 
-      Field AtoN = in;
+      Field AtoN(in._grid);
+      Field Mtmp(in._grid);
+      AtoN = in;
       out = AtoN*Coeffs[0];
-
+      //      std::cout <<"Poly in " <<norm2(in)<<std::endl;
+      //      std::cout <<"0 " <<norm2(out)<<std::endl;
       for(int n=1;n<Coeffs.size();n++){
-	Field Mtmp=AtoN;
+	Mtmp = AtoN;
-	Linop.Op(Mtmp,AtoN);
+	Linop.HermOp(Mtmp,AtoN);
 	out=out+AtoN*Coeffs[n];
+	//	std::cout << n<<" " <<norm2(out)<<std::endl;
       }
     };
   };
@@ -36,15 +47,15 @@ namespace Grid {
   template<class Field>
   class Chebyshev : public OperatorFunction<Field> {
   private:
-    std::vector<double> Coeffs;
+    std::vector<RealD> Coeffs;
     int order;
-    double hi;
+    RealD hi;
-    double lo;
+    RealD lo;
 
   public:
     void csv(std::ostream &out){
-      for (double x=lo; x<hi; x+=(hi-lo)/1000) {
+      for (RealD x=lo; x<hi; x+=(hi-lo)/1000) {
-	double f = approx(x);
+	RealD f = approx(x);
 	out<< x<<" "<<f<<std::endl;
       }
       return;
@@ -53,15 +64,19 @@ namespace Grid {
     // 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){
+      for(RealD x=lo;x<hi;x+=(hi-lo)/50.0){
 	out <<x<<"\t"<<approx(x)<<std::endl;
       }
     };
 
+    Chebyshev(){};
+    Chebyshev(RealD _lo,RealD _hi,int _order, RealD (* func)(RealD) ) {Init(_lo,_hi,_order,func);};
     
+    ////////////////////////////////////////////////////////////////////////////////////////////////////
     // c.f. numerical recipes "chebft"/"chebev". This is sec 5.8 "Chebyshev approximation".
-    //
+    ////////////////////////////////////////////////////////////////////////////////////////////////////
-    Chebyshev(double _lo,double _hi,int _order, double (* func)(double) ){
+    void Init(RealD _lo,RealD _hi,int _order, RealD (* func)(RealD))
+    {
       lo=_lo;
       hi=_hi;
       order=_order;
@@ -69,24 +84,26 @@ namespace Grid {
       if(order < 2) exit(-1);
       Coeffs.resize(order);
       for(int j=0;j<order;j++){
-	double s=0;
+	RealD s=0;
 	for(int k=0;k<order;k++){
-	  double y=std::cos(M_PI*(k+0.5)/order);
+	  RealD y=std::cos(M_PI*(k+0.5)/order);
-	  double x=0.5*(y*(hi-lo)+(hi+lo));
+	  RealD x=0.5*(y*(hi-lo)+(hi+lo));
-	  double f=func(x);
+	  RealD f=func(x);
 	  s=s+f*std::cos( j*M_PI*(k+0.5)/order );
 	}
 	Coeffs[j] = s * 2.0/order;
       }
     };
+
+    
     void JacksonSmooth(void){
-      double M=order;
+      RealD M=order;
-      double alpha = M_PI/(M+2);
+      RealD alpha = M_PI/(M+2);
-      double lmax = std::cos(alpha);
+      RealD lmax = std::cos(alpha);
-      double sumUsq =0;
+      RealD sumUsq =0;
-      std::vector<double> U(M);
+      std::vector<RealD> U(M);
-      std::vector<double> a(M);
+      std::vector<RealD> a(M);
-      std::vector<double> g(M);
+      std::vector<RealD> 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];
@@ -107,18 +124,18 @@ namespace Grid {
 	Coeffs[m]*=g[m];
       }
     }
-    double approx(double x) // Convenience for plotting the approximation
+    RealD approx(RealD x) // Convenience for plotting the approximation
     {
-      double Tn;
+      RealD Tn;
-      double Tnm;
+      RealD Tnm;
-      double Tnp;
+      RealD Tnp;
       
-      double y=( x-0.5*(hi+lo))/(0.5*(hi-lo));
+      RealD y=( x-0.5*(hi+lo))/(0.5*(hi-lo));
       
-      double T0=1;
+      RealD T0=1;
-      double T1=y;
+      RealD T1=y;
       
-      double sum;
+      RealD sum;
       sum = 0.5*Coeffs[0]*T0;
       sum+= Coeffs[1]*T1;
       
@@ -151,8 +168,8 @@ namespace Grid {
 
       std::cout<<GridLogMessage << "Chebyshev ["<<lo<<","<<hi<<"]"<< " order "<<order <<std::endl;
       // Tn=T1 = (xscale M + mscale)in
-      double xscale = 2.0/(hi-lo);
+      RealD xscale = 2.0/(hi-lo);
-      double mscale = -(hi+lo)/(hi-lo);
+      RealD mscale = -(hi+lo)/(hi-lo);
       Linop.HermOp(T0,y);
       T1=y*xscale+in*mscale;
 
@@ -179,5 +196,121 @@ namespace Grid {
   };
 
 
+  template<class Field>
+  class ChebyshevLanczos : public Chebyshev<Field> {
+  private:
+    std::vector<RealD> Coeffs;
+    int order;
+    RealD alpha;
+    RealD beta;
+    RealD mu;
+
+  public:
+    ChebyshevLanczos(RealD _alpha,RealD _beta,RealD _mu,int _order) :
+    alpha(_alpha),
+      beta(_beta),
+          mu(_mu)
+    {
+      order=_order;
+      Coeffs.resize(order);
+      for(int i=0;i<_order;i++){
+	Coeffs[i] = 0.0;
+      }
+      Coeffs[order-1]=1.0;
+    };
+
+    void csv(std::ostream &out){
+      for (RealD x=-1.2*alpha; x<1.2*alpha; x+=(2.0*alpha)/10000) {
+	RealD f = approx(x);
+	out<< x<<" "<<f<<std::endl;
+      }
+      return;
+    }
+
+    RealD approx(RealD xx) // Convenience for plotting the approximation
+    {
+      RealD Tn;
+      RealD Tnm;
+      RealD Tnp;
+      Real aa = alpha * alpha;
+      Real bb = beta  *  beta;
+      
+      RealD x = ( 2.0 * (xx-mu)*(xx-mu) - (aa+bb) ) / (aa-bb);
+
+      RealD y= x;
+      
+      RealD T0=1;
+      RealD T1=y;
+      
+      RealD sum;
+      sum = 0.5*Coeffs[0]*T0;
+      sum+= Coeffs[1]*T1;
+      
+      Tn =T1;
+      Tnm=T0;
+      for(int i=2;i<order;i++){
+	Tnp=2*y*Tn-Tnm;
+	Tnm=Tn;
+	Tn =Tnp;
+	sum+= Tn*Coeffs[i];
+      }
+      return sum;
+    };
+
+    // shift_Multiply in Rudy's code
+    void AminusMuSq(LinearOperatorBase<Field> &Linop, const Field &in, Field &out) 
+    {
+      GridBase *grid=in._grid;
+      Field tmp(grid);
+
+      RealD aa= alpha*alpha;
+      RealD bb= beta * beta;
+
+      Linop.HermOp(in,out);
+      out = out - mu*in;
+
+      Linop.HermOp(out,tmp);
+      tmp = tmp - mu * out;
+
+      out = (2.0/ (aa-bb) ) * tmp -  ((aa+bb)/(aa-bb))*in;
+    };
+    // Implement the required interface
+    void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) {
+
+      GridBase *grid=in._grid;
+
+      int vol=grid->gSites();
+
+      Field T0(grid); T0 = in;  
+      Field T1(grid); 
+      Field T2(grid);
+      Field  y(grid);
+      
+      Field *Tnm = &T0;
+      Field *Tn  = &T1;
+      Field *Tnp = &T2;
+
+      // Tn=T1 = (xscale M )*in
+      AminusMuSq(Linop,T0,T1);
+
+      // sum = .5 c[0] T0 + c[1] T1
+      out = (0.5*Coeffs[0])*T0 + Coeffs[1]*T1;
+      for(int n=2;n<order;n++){
+	
+	AminusMuSq(Linop,*Tn,y);
+
+	*Tnp=2.0*y-(*Tnm);
+
+	out=out+Coeffs[n]* (*Tnp);
+
+	// Cycle pointers to avoid copies
+	Field *swizzle = Tnm;
+	Tnm    =Tn;
+	Tn     =Tnp;
+	Tnp    =swizzle;
+	  
+      }
+    }
+  };
 }
 #endif