Addtional routines for Lanczos (SYM2, Chebyshef)..

lib/algorithms/LinearOperator.h

@@ -270,6 +270,35 @@ namespace Grid {
       }
     };
 
+    template<class Matrix,class Field>
+      class SchurDiagTwoOperator :  public SchurOperatorBase<Field> {
+    protected:
+      Matrix &_Mat;
+    public:
+      SchurDiagTwoOperator (Matrix &Mat): _Mat(Mat){};
+
+      virtual  RealD Mpc      (const Field &in, Field &out) {
+	Field tmp(in._grid);
+
+	_Mat.MooeeInv(in,out);
+	_Mat.Meooe(out,tmp);
+	_Mat.MooeeInv(tmp,out);
+	_Mat.Meooe(out,tmp);
+
+	return axpy_norm(out,-1.0,tmp,in);
+      }
+      virtual  RealD MpcDag   (const Field &in, Field &out){
+	Field tmp(in._grid);
+
+	_Mat.MeooeDag(in,out);
+	_Mat.MooeeInvDag(out,tmp);
+	_Mat.MeooeDag(tmp,out);
+	_Mat.MooeeInvDag(out,tmp);
+
+	return axpy_norm(out,-1.0,tmp,in);
+      }
+    };
+
 
     /////////////////////////////////////////////////////////////
     // Base classes for functions of operators

lib/algorithms/approx/Chebyshev.h

@@ -83,7 +83,8 @@ namespace Grid {
 
   public:
     void csv(std::ostream &out){
-      for (RealD x=lo; x<hi; x+=(hi-lo)/1000) {
+	RealD diff = hi-lo;
+      for (RealD x=lo-0.2*diff; x<hi+0.2*diff; x+=(hi-lo)/1000) {
 	RealD f = approx(x);
 	out<< x<<" "<<f<<std::endl;
       }
@@ -100,10 +101,24 @@ namespace Grid {
 
     Chebyshev(){};
     Chebyshev(RealD _lo,RealD _hi,int _order, RealD (* func)(RealD) ) {Init(_lo,_hi,_order,func);};
-    
+    Chebyshev(RealD _lo,RealD _hi,int _order) {Init(_lo,_hi,_order);};
+
     ////////////////////////////////////////////////////////////////////////////////////////////////////
     // c.f. numerical recipes "chebft"/"chebev". This is sec 5.8 "Chebyshev approximation".
     ////////////////////////////////////////////////////////////////////////////////////////////////////
+// CJ: the one we need for Lanczos
+    void Init(RealD _lo,RealD _hi,int _order)
+    {
+      lo=_lo;
+      hi=_hi;
+      order=_order;
+      
+      if(order < 2) exit(-1);
+      Coeffs.resize(order);
+      Coeffs.assign(0.,order);
+      Coeffs[order-1] = 1.;
+    };
+
     void Init(RealD _lo,RealD _hi,int _order, RealD (* func)(RealD))
     {
       lo=_lo;

lib/algorithms/iterative/EigenSort.h

@@ -38,6 +38,8 @@ template<class Field>
 class SortEigen {
  private:
   
+//hacking for testing for now
+#if 0
   static bool less_lmd(RealD left,RealD right){
     return fabs(left) < fabs(right);
   }  
@@ -45,6 +47,15 @@ class SortEigen {
 		 std::pair<RealD,Field>& right){
     return fabs(left.first) < fabs(right.first);
   }  
+#else
+  static bool less_lmd(RealD left,RealD right){
+    return left > right;
+  }  
+  static bool less_pair(std::pair<RealD,Field>& left,
+		 std::pair<RealD,Field>& right){
+    return left.first > (right.first);
+  }  
+#endif
   
  public:
 

lib/algorithms/iterative/ImplicitlyRestartedLanczos.h

@@ -241,20 +241,15 @@ public:
   double evals_tmp[NN];
   double evec_tmp[NN][NN];
   memset(evec_tmp[0],0,sizeof(double)*NN*NN);
-  double AA[NN][NN];
+//  double AA[NN][NN];
   double DD[NN];
   double EE[NN];
-  memset(AA, 0, sizeof(double)*NN*NN);
-  //    memset(ZZ, 0, sizeof(double)*NN*NN);
   for (int i = 0; i< NN; i++)
     for (int j = i - 1; j <= i + 1; j++)
       if ( j < NN && j >= 0 ) {
-//        AA[i][j] = AH(i,j);
         if (i==j) DD[i] = lmd[i];
         if (i==j) evals_tmp[i] = lmd[i];
         if (j==(i-1)) EE[j] = lme[j];
-        //        if (i<20 && j<20)
-//        std:: cout << "AA["<<i<<"]["<<j<<"]="<<AA[i][j]<<endl;
       }
   int evals_found;
   int lwork = ( (18*NN) > (1+4*NN+NN*NN)? (18*NN):(1+4*NN+NN*NN)) ;
@@ -301,7 +296,6 @@ public:
         }
       }
       {
-//        TIMER("QMP_sum_double_array");
 //        QMP_sum_double_array(evals_tmp,NN);
 //        QMP_sum_double_array((double *)evec_tmp,NN*NN);
          _fgrid.GlobalSumVector(evals_tmp,NN);

tests/Test_dwf_lanczos.cc

@@ -31,6 +31,8 @@ using namespace std;
 using namespace Grid;
 using namespace Grid::QCD;
 
+RealD AllZero(RealD x){ return 0.;}
+
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
@@ -64,14 +66,18 @@ int main (int argc, char ** argv)
 
   const int Nstop = 30;
   const int Nk = 40;
-  const int Np = 100;
+  const int Np = 10;
   const int Nm = Nk+Np;
   const int MaxIt= 10000;
   RealD resid = 1.0e-8;
 
   std::vector<double> Coeffs { 0.,1.};
   Polynomial<LatticeFermion> PolyX(Coeffs);
-  ImplicitlyRestartedLanczos<LatticeFermion> IRL(*FGrid,HermOp,PolyX,Nstop,Nk,Nm,resid,MaxIt);
+  Chebyshev<LatticeFermion> Cheb(1,80.,11);
+//  ChebyshevLanczos<LatticeFermion> Cheb(9.,1.,0.,20);
+  Cheb.csv(std::cout);
+//  exit(-24);
+  ImplicitlyRestartedLanczos<LatticeFermion> IRL(*FGrid,HermOp,Cheb,Nstop,Nk,Nm,resid,MaxIt);
 
   
   std::vector<RealD>          eval(Nm);