Added flags for explicit eval checks at the end of operator()

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;
       }
     }

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;
       }
     }

Grid/algorithms/iterative/KrylovSchur.h

@@ -338,6 +338,7 @@ class KrylovSchur {
 
   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),
@@ -606,6 +607,20 @@ if (!shift){
           // 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;

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);