diff --git a/Grid/algorithms/iterative/BlockKrylovSchur.h b/Grid/algorithms/iterative/BlockKrylovSchur.h
index c48c8b360..6c599f612 100644
--- a/Grid/algorithms/iterative/BlockKrylovSchur.h
+++ b/Grid/algorithms/iterative/BlockKrylovSchur.h
@@ -123,6 +123,7 @@ protected:
 
 public:
   std::vector<Field> evecs;
+  bool doEvalCheck = false;
 
   //--------------------------------------------------------------------
   // Constructor
@@ -252,6 +253,21 @@ public:
         std::cout << GridLogMessage << "BlockKrylovSchur: done after " << iter
                   << " restarts, " << Nconv << " converged." << std::endl;
         std::cout << GridLogMessage << "Eigenvalues: " << evals.transpose() << std::endl;
+
+        if (doEvalCheck) {
+          Field w(Grid_);
+          for (int k = 0; k < (int)evecs.size(); k++) {
+            Linop.Op(evecs[k], w);
+            ComplexD eval_est = toStdCmplx(innerProduct(evecs[k], w));
+            w -= eval_est * evecs[k];
+            RealD res = std::sqrt(norm2(w));
+            std::cout << GridLogMessage << "BlockKrylovSchur: evec[" << k << "]"
+                      << "  eval_reported = " << evals[k]
+                      << "  eval_est = " << eval_est
+                      << "  || A v - eval_est * v || = " << res << std::endl;
+          }
+        }
+
         return;
       }
     }
diff --git a/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h b/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h
index 9a179c1ff..1b405f975 100644
--- a/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h
+++ b/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h
@@ -125,6 +125,7 @@ class HarmonicBlockKrylovSchur {
 
 public:
   std::vector<Field>  evecs;
+  bool doEvalCheck = false;
 
   //--------------------------------------------------------------------
   // Constructor
@@ -253,6 +254,21 @@ public:
                   << "HarmonicBlockKrylovSchur: done after " << iter
                   << " restarts, " << Nconv << " converged." << std::endl;
         std::cout << GridLogMessage << "Eigenvalues: " << evals.transpose() << std::endl;
+
+        if (doEvalCheck) {
+          Field w(Grid_);
+          for (int k = 0; k < (int)evecs.size(); k++) {
+            Linop.Op(evecs[k], w);
+            ComplexD eval_est = toStdCmplx(innerProduct(evecs[k], w));
+            w -= eval_est * evecs[k];
+            RealD res = std::sqrt(norm2(w));
+            std::cout << GridLogMessage << "HarmonicBlockKrylovSchur: evec[" << k << "]"
+                      << "  eval_reported = " << evals[k]
+                      << "  eval_est = " << eval_est
+                      << "  || A v - eval_est * v || = " << res << std::endl;
+          }
+        }
+
         return;
       }
     }
diff --git a/Grid/algorithms/iterative/KrylovSchur.h b/Grid/algorithms/iterative/KrylovSchur.h
index 0ebf29d24..bd97e8529 100644
--- a/Grid/algorithms/iterative/KrylovSchur.h
+++ b/Grid/algorithms/iterative/KrylovSchur.h
@@ -336,8 +336,9 @@ class KrylovSchur {
 
     RitzFilter ritzFilter;                  // how to sort evals
 
-  public:       
+  public:
     RealD *shift;
+    bool doEvalCheck = false;
 
     KrylovSchur(LinearOperatorBase<Field> &_Linop, GridBase *_Grid, RealD _Tolerance, RitzFilter filter = EvalReSmall)
       : Linop(_Linop), Grid(_Grid), Tolerance(_Tolerance), ritzFilter(filter), u(_Grid), MaxIter(-1), Nm(-1), Nk(-1), Nstop (-1),
@@ -601,11 +602,25 @@ if (!shift){
         int Nconv = converged();
         std::cout << GridLogMessage << "Number of evecs converged: " << Nconv << std::endl;
         if (Nconv >= Nstop || i == MaxIter - 1) {
-          std::cout << GridLogMessage << "Converged with " << Nconv << " / " << Nstop << " eigenvectors on iteration " 
+          std::cout << GridLogMessage << "Converged with " << Nconv << " / " << Nstop << " eigenvectors on iteration "
                         << i << "." << std::endl;
           // basisRotate(evecs, Qt, 0, Nk, 0, Nk, Nm);      // Think this might have been the issue
           std::cout << GridLogMessage << "Eigenvalues: " << evals << std::endl;
 
+          if (doEvalCheck) {
+            Field w(Grid);
+            for (int k = 0; k < (int)evecs.size(); k++) {
+              Linop.Op(evecs[k], w);
+              ComplexD eval_est = toStdCmplx(innerProduct(evecs[k], w));
+              w -= eval_est * evecs[k];
+              RealD res = std::sqrt(norm2(w));
+              std::cout << GridLogMessage << "KrylovSchur: evec[" << k << "]"
+                        << "  eval_reported = " << evals[k]
+                        << "  eval_est = " << eval_est
+                        << "  || A v - eval_est * v || = " << res << std::endl;
+            }
+          }
+
           // writeEigensystem(path);
 
           return;
diff --git a/examples/Example_krylov_schur.cc b/examples/Example_krylov_schur.cc
index 991b204ce..08a1e3b04 100644
--- a/examples/Example_krylov_schur.cc
+++ b/examples/Example_krylov_schur.cc
@@ -335,7 +335,7 @@ int main (int argc, char ** argv)
   // Run KrylovSchur and Arnoldi on a Hermitian matrix
   std::cout << GridLogMessage << "Running Krylov Schur" << std::endl;
     RealD shift=LanParams.shift;
-#if 0
+#if 1
     KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall);
 //    KrySchur(src[0], maxIter, Nm, Nk, Nstop);
     KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift);