Added explicit shift before pulling

examples/Example_krylov_schur.cc

@@ -105,255 +105,47 @@ template <class T> void writeFile(T& in, std::string const fname){
 typedef WilsonFermionD WilsonOp;
 typedef typename WilsonFermionD::FermionField FermionField;
 
-
-#if 0
-// Hermitize a DWF operator by squaring it
 template<class Matrix,class Field>
-class SquaredLinearOperator : public LinearOperatorBase<Field> {
+class InvertNonHermitianLinearOperator : public LinearOperatorBase<Field> {
-
-  public:
   Matrix &_Mat;
-
+  RealD _stp;
-  public:
-    SquaredLinearOperator(Matrix &Mat): _Mat(Mat) {};
-
-    void OpDiag (const Field &in, Field &out) {    assert(0);  }
-    void OpDir  (const Field &in, Field &out,int dir,int disp) {    assert(0);  }
-    void OpDirAll  (const Field &in, std::vector<Field> &out){    assert(0);  };
-    void Op     (const Field &in, Field &out){
-      // std::cout << "Op is overloaded as HermOp" << std::endl;
-      HermOp(in, out);
-    }
-    void AdjOp     (const Field &in, Field &out){
-      HermOp(in, out);
-    }
-    void _Op     (const Field &in, Field &out){
-      // std::cout << "Op: M "<<std::endl;
-      _Mat.M(in, out);
-    }
-    void _AdjOp     (const Field &in, Field &out){
-      // std::cout << "AdjOp: Mdag "<<std::endl;
-      _Mat.Mdag(in, out);
-    }
-    void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){    assert(0);  }
-    void HermOp(const Field &in, Field &out){
-      // std::cout << "HermOp: Mdag M Mdag M"<<std::endl;
-      Field tmp(in.Grid());
-      _Op(in,tmp);
-      _AdjOp(tmp,out);
-    }
-};
-
-template<class Matrix,class Field>
-class PVdagMLinearOperator : public LinearOperatorBase<Field> {
-  Matrix &_Mat;
-  Matrix &_PV;
 public:
-  PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV){};
+  InvertNonHermitianLinearOperator(Matrix &Mat,RealD stp=1e-8): _Mat(Mat),_stp(stp){};
-
+  // Support for coarsening to a multigrid
-  void OpDiag (const Field &in, Field &out) {    assert(0);  }
+  void OpDiag (const Field &in, Field &out) {
-  void OpDir  (const Field &in, Field &out,int dir,int disp) {    assert(0);  }
+//    _Mat.Mdiag(in,out);
-  void OpDirAll  (const Field &in, std::vector<Field> &out){    assert(0);  };
+//    out = out + shift*in;
+    assert(0);
+  }
+  void OpDir  (const Field &in, Field &out,int dir,int disp) {
+//    _Mat.Mdir(in,out,dir,disp);
+    assert(0);
+  }
+  void OpDirAll  (const Field &in, std::vector<Field> &out){
+//    _Mat.MdirAll(in,out);
+    assert(0);
+  };
   void Op     (const Field &in, Field &out){
-    std::cout << "Op: PVdag M "<<std::endl;
+//    _Mat.M(in,out);
-    Field tmp(in.Grid());
+//  RealD mass=-shift;
-    _Mat.M(in,tmp);
+//  WilsonCloverFermionD Dw(Umu, Grid, RBGrid, mass, csw_r, csw_t);
-    _PV.Mdag(tmp,out);
+  NonHermitianLinearOperator<Matrix,Field> HermOp(_Mat);
+  BiCGSTAB<Field> CG(_stp,10000);
+  CG(HermOp,in,out);
+//    out = out + shift * in;
   }
   void AdjOp     (const Field &in, Field &out){
-    std::cout << "AdjOp: Mdag PV "<<std::endl;
+    _Mat.Mdag(in,out);
-    Field tmp(in.Grid());
+//    out = out + shift * in;
-    _PV.M(in,tmp);
+  }
-    _Mat.Mdag(tmp,out);
+  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
+    assert(0);
   }
-  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){    assert(0);  }
   void HermOp(const Field &in, Field &out){
-    std::cout << "HermOp: Mdag PV PVdag M"<<std::endl;
+    assert(0);
-    Field tmp(in.Grid());
-    //    _Mat.M(in,tmp);
-    //    _PV.Mdag(tmp,out);
-    //    _PV.M(out,tmp);
-    //    _Mat.Mdag(tmp,out);
-    Op(in,tmp);
-    AdjOp(tmp,out);
-    //    std::cout << "HermOp done "<<norm2(out)<<std::endl;
   }
 };
 
-template<class Matrix,class Field>
-class ShiftedPVdagMLinearOperator : public LinearOperatorBase<Field> {
-  Matrix &_Mat;
-  Matrix &_PV;
-  RealD shift;
-public:
-  ShiftedPVdagMLinearOperator(RealD _shift,Matrix &Mat,Matrix &PV): shift(_shift),_Mat(Mat),_PV(PV){};
-
-  void OpDiag (const Field &in, Field &out) {    assert(0);  }
-  void OpDir  (const Field &in, Field &out,int dir,int disp) {    assert(0);  }
-  void OpDirAll  (const Field &in, std::vector<Field> &out){    assert(0);  };
-  void Op     (const Field &in, Field &out){
-    std::cout << "Op: PVdag M "<<std::endl;
-    Field tmp(in.Grid());
-    _Mat.M(in,tmp);
-    _PV.Mdag(tmp,out);
-    out = out + shift * in;
-  }
-  void AdjOp     (const Field &in, Field &out){
-    std::cout << "AdjOp: Mdag PV "<<std::endl;
-    Field tmp(in.Grid());
-    _PV.M(tmp,out);
-    _Mat.Mdag(in,tmp);
-    out = out + shift * in;
-  }
-  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){    assert(0);  }
-  void HermOp(const Field &in, Field &out){
-    std::cout << "HermOp: Mdag PV PVdag M"<<std::endl;
-    Field tmp(in.Grid());
-    Op(in,tmp);
-    AdjOp(tmp,out);
-  }
-};
-
-template<class Matrix, class Field>
-class ShiftedComplexPVdagMLinearOperator : public LinearOperatorBase<Field> {
-  Matrix &_Mat;
-  Matrix &_PV;
-  ComplexD shift;
-public:
-ShiftedComplexPVdagMLinearOperator(ComplexD _shift,Matrix &Mat,Matrix &PV): shift(_shift),_Mat(Mat),_PV(PV){};
-
-  void OpDiag (const Field &in, Field &out) {    assert(0);  }
-  void OpDir  (const Field &in, Field &out,int dir,int disp) {    assert(0);  }
-  void OpDirAll  (const Field &in, std::vector<Field> &out){    assert(0);  };
-  void Op     (const Field &in, Field &out){
-    std::cout << "Op: PVdag M "<<std::endl;
-    Field tmp(in.Grid());
-    _Mat.M(in,tmp);
-    _PV.Mdag(tmp,out);
-    out = out + shift * in;
-  }
-  void AdjOp     (const Field &in, Field &out){
-    std::cout << "AdjOp: Mdag PV "<<std::endl;
-    Field tmp(in.Grid());
-    _PV.M(tmp,out);
-    _Mat.Mdag(in,tmp);
-    out = out + shift * in;
-  }
-  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){    assert(0);  }
-  void HermOp(const Field &in, Field &out){
-    std::cout << "HermOp: Mdag PV PVdag M"<<std::endl;
-    Field tmp(in.Grid());
-    Op(in,tmp);
-    AdjOp(tmp,out);
-  }
-  
-  void resetShift(ComplexD newShift) {
-    shift = newShift;
-  }
-};
-
-template<class Fobj,class CComplex,int nbasis>
-class MGPreconditioner : public LinearFunction< Lattice<Fobj> > {
-public:
-  using LinearFunction<Lattice<Fobj> >::operator();
-
-  typedef Aggregation<Fobj,CComplex,nbasis> Aggregates;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::FineField    FineField;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseVector CoarseVector;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseMatrix CoarseMatrix;
-  typedef LinearOperatorBase<FineField>                            FineOperator;
-  typedef LinearFunction    <FineField>                            FineSmoother;
-  typedef LinearOperatorBase<CoarseVector>                         CoarseOperator;
-  typedef LinearFunction    <CoarseVector>                         CoarseSolver;
-  Aggregates     & _Aggregates;
-  FineOperator   & _FineOperator;
-  FineSmoother   & _PreSmoother;
-  FineSmoother   & _PostSmoother;
-  CoarseOperator & _CoarseOperator;
-  CoarseSolver   & _CoarseSolve;
-
-  int    level;  void Level(int lv) {level = lv; };
-
-  MGPreconditioner(Aggregates &Agg,
-		   FineOperator &Fine,
-		   FineSmoother &PreSmoother,
-		   FineSmoother &PostSmoother,
-		   CoarseOperator &CoarseOperator_,
-		   CoarseSolver &CoarseSolve_)
-    : _Aggregates(Agg),
-      _FineOperator(Fine),
-      _PreSmoother(PreSmoother),
-      _PostSmoother(PostSmoother),
-      _CoarseOperator(CoarseOperator_),
-      _CoarseSolve(CoarseSolve_),
-      level(1)  {  }
-
-  virtual void operator()(const FineField &in, FineField & out) 
-  {
-    GridBase *CoarseGrid = _Aggregates.CoarseGrid;
-    //    auto CoarseGrid = _CoarseOperator.Grid();
-    CoarseVector Csrc(CoarseGrid);
-    CoarseVector Csol(CoarseGrid);
-    FineField vec1(in.Grid());
-    FineField vec2(in.Grid());
-
-    std::cout<<GridLogMessage << "Calling PreSmoother " <<std::endl;
-
-    //    std::cout<<GridLogMessage << "Calling PreSmoother input residual "<<norm2(in) <<std::endl;
-    double t;
-    // Fine Smoother
-    //    out = in;
-    out = Zero();
-    t=-usecond();
-    _PreSmoother(in,out);
-    t+=usecond();
-
-    std::cout<<GridLogMessage << "PreSmoother took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Update the residual
-    _FineOperator.Op(out,vec1);  sub(vec1, in ,vec1);   
-    //    std::cout<<GridLogMessage <<"Residual-1 now " <<norm2(vec1)<<std::endl;
-
-    // Fine to Coarse 
-    t=-usecond();
-    _Aggregates.ProjectToSubspace  (Csrc,vec1);
-    t+=usecond();
-    std::cout<<GridLogMessage << "Project to coarse took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Coarse correction
-    t=-usecond();
-    Csol = Zero();
-    _CoarseSolve(Csrc,Csol);
-    //Csol=Zero();
-    t+=usecond();
-    std::cout<<GridLogMessage << "Coarse solve took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Coarse to Fine
-    t=-usecond();  
-    //    _CoarseOperator.PromoteFromSubspace(_Aggregates,Csol,vec1);
-    _Aggregates.PromoteFromSubspace(Csol,vec1); 
-    add(out,out,vec1);
-    t+=usecond();
-    std::cout<<GridLogMessage << "Promote to this level took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Residual
-    _FineOperator.Op(out,vec1);  sub(vec1 ,in , vec1);  
-    //    std::cout<<GridLogMessage <<"Residual-2 now " <<norm2(vec1)<<std::endl;
-
-    // Fine Smoother
-    t=-usecond();
-    //    vec2=vec1;
-    vec2=Zero();
-    _PostSmoother(vec1,vec2);
-    t+=usecond();
-    std::cout<<GridLogMessage << "PostSmoother took "<< t/1000.0<< "ms" <<std::endl;
-
-    add( out,out,vec2);
-    std::cout<<GridLogMessage << "Done " <<std::endl;
-  }
-};
-#endif
-
 template<class Field>
 void testSchurFromHess(Arnoldi<Field>& Arn, Field& src, int Nlarge, int Nm, int Nk) {
 
@@ -500,6 +292,7 @@ int main (int argc, char ** argv)
 
 
   NonHermitianLinearOperator<WilsonOp,FermionField> Dwilson(WilsonOperator); /// <-----
+  InvertNonHermitianLinearOperator<WilsonOp,FermionField> Iwilson(WilsonOperator); /// <-----
   MdagMLinearOperator<WilsonOp,FermionField> HermOp(WilsonOperator); /// <-----
   Gamma5HermitianLinearOperator <WilsonOp,LatticeFermion> HermOp2(WilsonOperator); /// <----
 
@@ -538,11 +331,14 @@ int main (int argc, char ** argv)
 //  KrySchur(src, maxIter, Nm, Nk, Nstop);
 //  KrylovSchur KrySchur (HermOp2, UGrid, resid,EvalNormSmall);
 //  Hacked, really EvalImagSmall
-  KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall);
+#if 0
-//  BlockKrylovSchur KrySchur (HermOp2, UGrid, Nu, resid,EvalNormSmall);
+    RealD shift=1.5;
-  RealD shift=2.5;
+    KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall);
-  KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift);
+    KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift);
-//  KrySchur(src[0], maxIter, Nm, Nk, Nstop);
+#else
+    KrylovSchur KrySchur (Iwilson, UGrid, resid,EvalImNormSmall);
+    KrySchur(src[0], maxIter, Nm, Nk, Nstop);
+#endif
   std::cout << GridLogMessage << "evec.size= " << KrySchur.evecs.size()<< std::endl;
 
   src[0]=KrySchur.evecs[0];