From 0145685f9658cf6a5f56248b958d804382fe1e77 Mon Sep 17 00:00:00 2001
From: Chulwoo Jung <chulwoo@bnl.gov>
Date: Fri, 18 Aug 2017 01:44:31 -0400
Subject: [PATCH] Added Staggered Type Preconditioned operator

---
 lib/algorithms/LinearOperator.h               | 52 +++++++++++++++++++
 .../iterative/ImplicitlyRestartedLanczosCJ.h  |  8 +--
 2 files changed, 56 insertions(+), 4 deletions(-)

diff --git a/lib/algorithms/LinearOperator.h b/lib/algorithms/LinearOperator.h
index 012ceb5b..d3d46785 100644
--- a/lib/algorithms/LinearOperator.h
+++ b/lib/algorithms/LinearOperator.h
@@ -300,6 +300,58 @@ namespace Grid {
       }
     };
 
+    template<class Matrix,class Field>
+      class SchurStagOperator :  public LinearOperatorBase<Field> {
+    protected:
+      Matrix &_Mat;
+    public:
+      SchurStagOperator (Matrix &Mat): _Mat(Mat){};
+      virtual  RealD Mpc      (const Field &in, Field &out) {
+	Field tmp(in._grid);
+	Field tmp2(in._grid);
+
+	_Mat.Mooee(in,out);
+	_Mat.Mooee(out,tmp);
+
+	_Mat.Meooe(in,out);
+	_Mat.Meooe(out,tmp2);
+
+	return axpy_norm(out,-1.0,tmp2,tmp);
+      }
+      virtual  RealD MpcDag   (const Field &in, Field &out){
+
+	return Mpc(in,out);
+      }
+#if 0
+      virtual void MpcDagMpc(const Field &in, Field &out,RealD &ni,RealD &no) {
+	Field tmp(in._grid);
+	ni=Mpc(in,tmp);
+	no=MpcDag(tmp,out);
+      }
+#endif
+      void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
+	n1 = Mpc(in,out);n2=0.;
+      }
+      void HermOp(const Field &in, Field &out){
+	RealD n1,n2;
+	HermOpAndNorm(in,out,n1,n2);
+      }
+      void Op     (const Field &in, Field &out){
+	Mpc(in,out);
+      }
+      void AdjOp     (const Field &in, Field &out){ 
+	MpcDag(in,out);
+      }
+      // Support for coarsening to a multigrid
+      void OpDiag (const Field &in, Field &out) {
+	assert(0); // must coarsen the unpreconditioned system
+      }
+      void OpDir  (const Field &in, Field &out,int dir,int disp) {
+	assert(0);
+      }
+
+    };
+
 #if 0
   // This is specific to (Z)mobius fermions
   template<class Matrix, class Field>
diff --git a/lib/algorithms/iterative/ImplicitlyRestartedLanczosCJ.h b/lib/algorithms/iterative/ImplicitlyRestartedLanczosCJ.h
index 800bd326..b1438f64 100644
--- a/lib/algorithms/iterative/ImplicitlyRestartedLanczosCJ.h
+++ b/lib/algorithms/iterative/ImplicitlyRestartedLanczosCJ.h
@@ -661,12 +661,12 @@ PARALLEL_FOR_LOOP
 	    RealD vv0 = norm2(v);
 	    eval2[i] = vnum/vden;
 	    v -= eval2[i]*B[i];
-	    RealD vv = norm2(v);
+	    RealD vv = norm2(v)/vv0;
 	    
 	    std::cout.precision(13);
 	    std::cout<<GridLogMessage << "[" << std::setw(3)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
 	    std::cout<<"eval = "<<std::setw(25)<< std::setiosflags(std::ios_base::left)<< eval2[i];
-	    std::cout<<"|H B[i] - eval[i]B[i]|^2 "<< std::setw(25)<< std::setiosflags(std::ios_base::right)<< vv;
+	    std::cout<<"|H B[i] - eval[i]B[i]|^2/|H B[i]|^2 "<< std::setw(25)<< std::setiosflags(std::ios_base::right)<< vv;
 	    std::cout<<" "<< vnum/(sqrt(vden)*sqrt(vv0)) << std::endl;
 	    
 	// change the criteria as evals are supposed to be sorted, all evals smaller(larger) than Nstop should have converged
@@ -696,12 +696,12 @@ void FinalCheck( int _Nk,
 		    RealD vv0 = norm2(v);
 		    RealD eval2 = vnum/vden;
 		    v -= eval2*evec[j];
-		    RealD vv = norm2(v);
+		    RealD vv = norm2(v)/vv0;
 	    
 		    std::cout.precision(13);
 		    std::cout<<GridLogMessage << "[" << std::setw(3)<< std::setiosflags(std::ios_base::right) <<j<<"] ";
 		    std::cout<<"eval = "<<std::setw(25)<< std::setiosflags(std::ios_base::left)<< eval2;
-		    std::cout<<"|H B[i] - eval[i]B[i]|^2 "<< std::setw(25)<< std::setiosflags(std::ios_base::right)<< vv;
+	    	std::cout<<"|H B[i] - eval[i]B[i]|^2/|H B[i]|^2 "<< std::setw(25)<< std::setiosflags(std::ios_base::right)<< vv;
 		    std::cout<<" "<< vnum/(sqrt(vden)*sqrt(vv0)) << std::endl;
 			eval[j] = eval2;