ALl codes compile against the new Lanczos call signature

lib/algorithms/LinearOperator.h

@@ -346,6 +346,7 @@ namespace Grid {
       virtual void operator() (const Field &in, Field &out) = 0;
     };
 
+
     /////////////////////////////////////////////////////////////
     // Base classes for Multishift solvers for operators
     /////////////////////////////////////////////////////////////
@@ -368,6 +369,64 @@ namespace Grid {
      };
     */
 
+  ////////////////////////////////////////////////////////////////////////////////////////////
+  // Hermitian operator Linear function and operator function
+  ////////////////////////////////////////////////////////////////////////////////////////////
+    template<class Field>
+      class HermOpOperatorFunction : public OperatorFunction<Field> {
+      void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) {
+	Linop.HermOp(in,out);
+      };
+    };
+
+    template<typename Field>
+      class PlainHermOp : public LinearFunction<Field> {
+    public:
+      LinearOperatorBase<Field> &_Linop;
+      
+      PlainHermOp(LinearOperatorBase<Field>& linop) : _Linop(linop) 
+      {}
+      
+      void operator()(const Field& in, Field& out) {
+	_Linop.HermOp(in,out);
+      }
+    };
+
+    template<typename Field>
+    class FunctionHermOp : public LinearFunction<Field> {
+    public:
+      OperatorFunction<Field>   & _poly;
+      LinearOperatorBase<Field> &_Linop;
+      
+      FunctionHermOp(OperatorFunction<Field> & poly,LinearOperatorBase<Field>& linop) 
+	: _poly(poly), _Linop(linop) {};
+      
+      void operator()(const Field& in, Field& out) {
+	_poly(_Linop,in,out);
+      }
+    };
+
+  template<class Field>
+  class Polynomial : public OperatorFunction<Field> {
+  private:
+    std::vector<RealD> Coeffs;
+  public:
+    Polynomial(std::vector<RealD> &_Coeffs) : Coeffs(_Coeffs) { };
+
+    // Implement the required interface
+    void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) {
+
+      Field AtoN(in._grid);
+      Field Mtmp(in._grid);
+      AtoN = in;
+      out = AtoN*Coeffs[0];
+      for(int n=1;n<Coeffs.size();n++){
+	Mtmp = AtoN;
+	Linop.HermOp(Mtmp,AtoN);
+	out=out+AtoN*Coeffs[n];
+      }
+    };
+  };
 
 }
 

lib/algorithms/approx/Chebyshev.h

@@ -34,41 +34,6 @@ Author: Christoph Lehner <clehner@bnl.gov>
 
 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<RealD> Coeffs;
-  public:
-    Polynomial(std::vector<RealD> &_Coeffs) : Coeffs(_Coeffs) { };
-
-    // Implement the required interface
-    void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) {
-
-      Field AtoN(in._grid);
-      Field Mtmp(in._grid);
-      AtoN = in;
-      out = AtoN*Coeffs[0];
-//            std::cout <<"Poly in " <<norm2(in)<<" size "<< Coeffs.size()<<std::endl;
-//            std::cout <<"Coeffs[0]= "<<Coeffs[0]<< " 0 " <<norm2(out)<<std::endl;
-      for(int n=1;n<Coeffs.size();n++){
-	Mtmp = AtoN;
-	Linop.HermOp(Mtmp,AtoN);
-	out=out+AtoN*Coeffs[n];
-//            std::cout <<"Coeffs "<<n<<"= "<< Coeffs[n]<< " 0 " <<std::endl;
-//		std::cout << n<<" " <<norm2(out)<<std::endl;
-      }
-    };
-  };
 
   ////////////////////////////////////////////////////////////////////////////////////////////
   // Generic Chebyshev approximations

lib/algorithms/iterative/ImplicitlyRestartedLanczos.h

@@ -186,9 +186,9 @@ public:
 			    int _Nk, // sought vecs
 			    int _Nm, // spare vecs
 			    RealD _eresid, // resid in lmdue deficit 
-			    RealD _betastp, // if beta(k) < betastp: converged
 			    int _MaxIter, // Max iterations
-			    int _MinRestart, int _orth_period = 1,
+			    RealD _betastp=0.0, // if beta(k) < betastp: converged
+			    int _MinRestart=1, int _orth_period = 1,
 			    IRLdiagonalisation _diagonalisation= IRLdiagonaliseWithEigen) :
   _HermOp(HermOp),      _HermOpTest(HermOpTest),
     Nstop(_Nstop)  ,      Nk(_Nk),      Nm(_Nm),
@@ -232,7 +232,7 @@ repeat
   →AVK =VKHK +fKe†K † Extend to an M = K + P step factorization AVM = VMHM + fMeM
 until convergence
 */
-  void calc(std::vector<RealD>& eval, std::vector<Field>& evec,  const Field& src, int& Nconv, bool reverse, int SkipTest)
+  void calc(std::vector<RealD>& eval, std::vector<Field>& evec,  const Field& src, int& Nconv, bool reverse=true, int SkipTest=0)
   {
     GridBase *grid = src._grid;
     assert(grid == evec[0]._grid);

tests/lanczos/Test_dwf_compressed_lanczos.cc

@@ -100,19 +100,6 @@ void write_history(char* fn, std::vector<RealD>& hist) {
   fclose(f);
 }
 
-template<typename Field>
-class FunctionHermOp : public LinearFunction<Field> {
-public:
-  OperatorFunction<Field>   & _poly;
-  LinearOperatorBase<Field> &_Linop;
-
-  FunctionHermOp(OperatorFunction<Field> & poly,LinearOperatorBase<Field>& linop) : _poly(poly), _Linop(linop) {
-  }
-
-  void operator()(const Field& in, Field& out) {
-    _poly(_Linop,in,out);
-  }
-};
 
 template<typename Field>
 class CheckpointedLinearFunction : public LinearFunction<Field> {
@@ -268,19 +255,6 @@ public:
   }
 };
 
-template<typename Field>
-class PlainHermOp : public LinearFunction<Field> {
-public:
-  LinearOperatorBase<Field> &_Linop;
-
-  PlainHermOp(LinearOperatorBase<Field>& linop) : _Linop(linop) {
-  }
-
-  void operator()(const Field& in, Field& out) {
-    _Linop.HermOp(in,out);
-  }
-};
-
 template<typename vtype, int N > using CoarseSiteFieldGeneral = iScalar< iVector<vtype, N> >;
 template<int N> using CoarseSiteFieldD = CoarseSiteFieldGeneral< vComplexD, N >;
 template<int N> using CoarseSiteFieldF = CoarseSiteFieldGeneral< vComplexF, N >;
@@ -326,7 +300,7 @@ void CoarseGridLanczos(BlockProjector<Field>& pr,RealD alpha2,RealD beta,int Npo
     Op2 = &Op2plain;
   }
   ProjectedHermOp<CoarseLatticeFermion<Nstop1>,LatticeFermion> Op2nopoly(pr,HermOp);
-  ImplicitlyRestartedLanczos<CoarseLatticeFermion<Nstop1> > IRL2(*Op2,*Op2,Nstop2,Nk2,Nm2,resid2,betastp2,MaxIt,MinRes2);
+  ImplicitlyRestartedLanczos<CoarseLatticeFermion<Nstop1> > IRL2(*Op2,*Op2,Nstop2,Nk2,Nm2,resid2,MaxIt,betastp2,MinRes2);
 
 
   src_coarse = 1.0;
@@ -648,7 +622,7 @@ int main (int argc, char ** argv) {
   }
 
   // First round of Lanczos to get low mode basis
-  ImplicitlyRestartedLanczos<LatticeFermion> IRL1(Op1,Op1test,Nstop1,Nk1,Nm1,resid1,betastp1,MaxIt,MinRes1);
+  ImplicitlyRestartedLanczos<LatticeFermion> IRL1(Op1,Op1test,Nstop1,Nk1,Nm1,resid1,MaxIt,betastp1,MinRes1);
   int Nconv;
 
   char tag[1024];

tests/lanczos/Test_dwf_lanczos.cc

@@ -84,11 +84,12 @@ int main (int argc, char ** argv)
 
   std::vector<double> Coeffs { 0.,-1.};
   Polynomial<FermionField> PolyX(Coeffs);
-  Chebyshev<FermionField> Cheb(0.2,5.,11);
+  Chebyshev<FermionField> Cheby(0.2,5.,11);
-//  ChebyshevLanczos<LatticeFermion> Cheb(9.,1.,0.,20);
+
-//  Cheb.csv(std::cout);
+  FunctionHermOp<FermionField> OpCheby(Cheby,HermOp);
-//  exit(-24);
+     PlainHermOp<FermionField> Op     (HermOp);
-  ImplicitlyRestartedLanczos<FermionField> IRL(HermOp,Cheb,Nstop,Nk,Nm,resid,MaxIt);
+
+  ImplicitlyRestartedLanczos<FermionField> IRL(OpCheby,Op,Nstop,Nk,Nm,resid,MaxIt);
 
   
   std::vector<RealD>          eval(Nm);

tests/lanczos/Test_synthetic_lanczos.cc

@@ -119,12 +119,13 @@ int main (int argc, char ** argv)
   RealD beta  = 0.1;
   RealD mu    = 0.0;
   int order = 11;
-  ChebyshevLanczos<LatticeComplex> Cheby(alpha,beta,mu,order);
+  Chebyshev<LatticeComplex> Cheby(alpha,beta,order);
   std::ofstream file("cheby.dat");
   Cheby.csv(file);
 
-  HermOpOperatorFunction<LatticeComplex> X;
   DumbOperator<LatticeComplex> HermOp(grid);
+  FunctionHermOp<LatticeComplex> OpCheby(Cheby,HermOp);
+     PlainHermOp<LatticeComplex> Op(HermOp);
 
   const int Nk = 40;
   const int Nm = 80;
@@ -133,8 +134,9 @@ int main (int argc, char ** argv)
   int Nconv;
   RealD eresid = 1.0e-6;
 
-  ImplicitlyRestartedLanczos<LatticeComplex> IRL(HermOp,X,Nk,Nk,Nm,eresid,Nit);
+
-  ImplicitlyRestartedLanczos<LatticeComplex> ChebyIRL(HermOp,Cheby,Nk,Nk,Nm,eresid,Nit);
+  ImplicitlyRestartedLanczos<LatticeComplex> IRL(Op,Op,Nk,Nk,Nm,eresid,Nit);
+  ImplicitlyRestartedLanczos<LatticeComplex> ChebyIRL(OpCheby,Op,Nk,Nk,Nm,eresid,Nit);
 
   LatticeComplex src(grid); gaussian(RNG,src);
   {

tests/lanczos/Test_wilson_lanczos.cc

@@ -86,9 +86,12 @@ int main(int argc, char** argv) {
 
   std::vector<double> Coeffs{0, 1.};
   Polynomial<FermionField> PolyX(Coeffs);
-  Chebyshev<FermionField> Cheb(0.0, 10., 12);
+  Chebyshev<FermionField> Cheby(0.0, 10., 12);
-  ImplicitlyRestartedLanczos<FermionField> IRL(HermOp, PolyX, Nstop, Nk, Nm,
+
-                                               resid, MaxIt);
+  FunctionHermOp<FermionField> OpCheby(Cheby,HermOp);
+     PlainHermOp<FermionField> Op     (HermOp);
+
+  ImplicitlyRestartedLanczos<FermionField> IRL(OpCheby, Op, Nstop, Nk, Nm, resid, MaxIt);
 
   std::vector<RealD> eval(Nm);
   FermionField src(FGrid);