diff --git a/Grid/algorithms/iterative/Arnoldi.h b/Grid/algorithms/iterative/Arnoldi.h
index 6232681d..a1cb2598 100644
--- a/Grid/algorithms/iterative/Arnoldi.h
+++ b/Grid/algorithms/iterative/Arnoldi.h
@@ -32,6 +32,53 @@ See the full license in the file "LICENSE" in the top level distribution directo
 
 NAMESPACE_BEGIN(Grid); 
 
+//Moved to KrylovSchur
+#if 0
+/**
+<<<<<<< HEAD
+ * Options for which Ritz values to keep in implicit restart.
+ */
+enum RitzFilter {
+  EvalNormSmall,           // Keep evals with smallest norm
+  EvalNormLarge,           // Keep evals with largest norm
+  EvalReSmall,             // Keep evals with smallest real part
+  EvalReLarge              // Keep evals with largest real part
+};
+
+// Select comparison function from RitzFilter
+struct ComplexComparator
+{
+  RitzFilter f;
+  ComplexComparator (RitzFilter _f) : f(_f) {}
+  bool operator()(std::complex<double> z1, std::complex<double> z2) { 
+    switch (f) {
+      RealD tmp1, tmp2;
+      tmp1=std::abs(std::imag(z1));
+      tmp2=std::abs(std::imag(z2));
+      case EvalNormSmall:
+        return std::abs(z1) < std::abs(z2);
+      case EvalNormLarge:
+        return std::abs(z1) > std::abs(z2);
+// Terrible hack
+//        return std::abs(std::real(z1)) < std::abs(std::real(z2));
+//	if ( std::abs(std::real(z1))  >4.) tmp1 +=1.;
+//	if ( std::abs(std::real(z2))  >4.) tmp2 +=1.;
+      case EvalReSmall:
+	  return tmp1 < tmp2;
+//        return std::abs(std::imag(z1)) < std::abs(std::imag(z2));
+      case EvalReLarge:
+	  return tmp1 > tmp2;
+//        return std::abs(std::real(z1)) > std::abs(std::real(z2));
+      default:
+        assert(0);
+    }
+  }
+};
+
+=======
+>>>>>>> 68af1bba67dd62881ead5ab1e54962a5486a0791
+#endif
+
 /**
  * Implementation of the Arnoldi algorithm.
  */
diff --git a/Grid/algorithms/iterative/KrylovSchur.h b/Grid/algorithms/iterative/KrylovSchur.h
deleted file mode 100644
index 642f6be3..00000000
--- a/Grid/algorithms/iterative/KrylovSchur.h
+++ /dev/null
@@ -1,852 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: ./lib/algorithms/iterative/KrylovSchur.h
-
-Copyright (C) 2015
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-Author: Patrick Oare <poare@bnl.gov>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-/*  END LEGAL */
-#ifndef GRID_KRYLOVSCHUR_H
-#define GRID_KRYLOVSCHUR_H
-
-NAMESPACE_BEGIN(Grid); 
-
-/**
- * Options for which Ritz values to keep in implicit restart. TODO move this and utilities into a new file
- */
-enum RitzFilter {
-  EvalNormSmall,           // Keep evals with smallest norm
-  EvalNormLarge,           // Keep evals with largest norm
-  EvalReSmall,             // Keep evals with smallest real part
-  EvalReLarge,             // Keep evals with largest real part
-  EvalImSmall,             // Keep evals with smallest imaginary part
-  EvalImLarge,             // Keep evals with largest imaginary part
-  EvalImNormSmall,         // Keep evals with smallest |imaginary| part
-  EvalImNormLarge,         // Keep evals with largest |imaginary| part
-};
-
-/** Selects the RitzFilter corresponding to the input string. */
-inline RitzFilter selectRitzFilter(std::string s) {
-  if (s == "EvalNormSmall")   { return EvalNormSmall; }   else
-  if (s == "EvalNormLarge")   { return EvalNormLarge; }   else
-  if (s == "EvalReSmall")     { return EvalReSmall; }     else
-  if (s == "EvalReLarge")     { return EvalReLarge; }     else
-  if (s == "EvalImSmall")     { return EvalImSmall; }     else
-  if (s == "EvalImLarge")     { return EvalImLarge; }     else 
-  if (s == "EvalImNormSmall") { return EvalImNormSmall; } else
-  if (s == "EvalImNormLarge") { return EvalImNormLarge; } else 
-  { assert(0); }
-}
-
-/** Returns a string saying which RitzFilter it is. */
-inline std::string rfToString(RitzFilter RF) {
-  switch (RF) {
-    case EvalNormSmall:
-      return "EvalNormSmall";
-    case EvalNormLarge:
-      return "EvalNormLarge";
-    case EvalReSmall:
-      return "EvalReSmall";
-    case EvalReLarge:
-      return "EvalReLarge";
-    case EvalImSmall:
-      return "EvalImSmall";
-    case EvalImLarge:
-      return "EvalImLarge";
-    case EvalImNormSmall:
-      return "EvalImNormSmall";
-    case EvalImNormLarge:
-      return "EvalImNormLarge";
-    default:
-      assert(0);
-  }
-}
-
-// Select comparison function from RitzFilter
-struct ComplexComparator
-{
-  RitzFilter RF;
-  ComplexComparator (RitzFilter _rf) : RF(_rf) {}
-  bool operator()(std::complex<double> z1, std::complex<double> z2) { 
-    switch (RF) {
-      case EvalNormSmall:
-        return std::abs(z1) < std::abs(z2);
-      case EvalNormLarge:
-        return std::abs(z1) > std::abs(z2);
-      case EvalReSmall:
-        return std::real(z1) < std::real(z2);     // DELETE THE ABS HERE!!!
-      case EvalReLarge:
-        return std::real(z1) > std::real(z2);
-      case EvalImSmall:
-        return std::imag(z1) < std::imag(z2);
-      case EvalImLarge:
-        return std::imag(z1) > std::imag(z2);
-      case EvalImNormSmall:
-        return std::abs(std::imag(z1)) < std::abs(std::imag(z2));
-      case EvalImNormLarge:
-        return std::abs(std::imag(z1)) > std::abs(std::imag(z2));
-      default:
-        assert(0);
-    }
-  }
-};
-
-/**
- * Computes a complex Schur decomposition of a complex matrix A using Eigen's matrix library. The Schur decomposition,
- *    A = Q^\dag S Q
- * factorizes A into a unitary matrix Q and an upper triangular matrix S. The eigenvalues of A lie on the diagonal of the upper triangular matrix S. 
- * The Schur decomposition is not unique: in particular, any ordering of the eigenvalues of A can be used as the diagonal of the matrix S. 
- * This class supports eigenvalue reordering by swapping two adjacent eigenvalues with a unitary transformation. 
- */
-class ComplexSchurDecomposition {
-
-  private:
-
-    typedef Eigen::MatrixXcd CMat;
-
-    CMat A;           // Matrix to decompose, A = Q^\dag S Q
-    CMat Q;           // Unitary matrix Q
-    CMat S;           // Upper triangular matrix S
-
-    // Placeholders for Givens rotation
-    CMat Givens;      // Givens rotation
-    ComplexD s;       // off-diagonal element
-    ComplexD lam1;    // First eval for swap
-    ComplexD lam2;    // Second eval for swap
-    ComplexD phi;     // phase of s
-    RealD r;          // norm of s and lam2 - lam1
-    
-    int Nm;                               // size of matrix problem
-
-    ComplexComparator cCompare;           // function to sort the Schur matrix. 
-  
-  public:
-
-    /**
-     * If the input matrix _A is in Hessenberg form (upper triangular + first subdiagonal non-zero), then the Schur 
-     * decomposition is easier to compute. 
-     */
-    ComplexSchurDecomposition(CMat _A, bool isHess, RitzFilter ritzFilter = EvalReSmall) : A(_A), Nm (_A.rows()), cCompare (ritzFilter)
-    {
-      Eigen::ComplexSchur<CMat> schur (Nm);
-      if (isHess) {
-        schur.computeFromHessenberg(_A, CMat::Identity(Nm, Nm), true);
-      } else {
-        schur.compute(_A, true);
-      }
-      S = schur.matrixT();
-      Q = schur.matrixU().adjoint();      // Eigen computes A = Q S Q^\dag, we want A = Q^\dag S Q
-    }
-
-    // Getters
-    int  getNm()      { return Nm; }         // size of matrix problem
-    CMat getMatrixA() { return A;  }         // matrix for decomposition
-    CMat getMatrixQ() { return Q;  }         // unitary matrix Q
-    CMat getMatrixS() { return S;  }         // Schur matrix (upper triangular) S
-    CMat getRitz()    { return S.diagonal(); }
-
-    /**
-     * Checks the Schur decomposition A = Q^\dag S Q holds for the computed matrices. Returns if the relative 
-     * Frobenius norm || A - Q^\dag S Q || / || A || is less than rtol. 
-     */
-    bool checkDecomposition(RealD rtol = 1e-8) {
-      RealD Anorm = A.norm();
-      if (Anorm < rtol) {
-        std::cout << GridLogMessage << "Zero matrix" << std::endl;
-        return true;
-      }
-
-      std::cout << GridLogDebug << "S = " << std::endl << S << std::endl;
-      std::cout << GridLogDebug << "Q = " << std::endl << Q << std::endl;
-
-      CMat A2 = Q.adjoint() * S * Q;
-      std::cout << GridLogDebug << "Q^dag S Q = " << std::endl << A2 << std::endl;
-      RealD dA = (A - A2).norm() / Anorm;
-      return (dA < rtol);
-    }
-
-    /**
-     * Swaps the components on the diagonal of the Schur matrix at index i with index i + 1. 
-     * Updates the orthogonal matrix Q accordingly. 
-     */
-    void swapEvals(int i) {
-
-      assert(0 <= i && i <= Nm - 1);        // can only swap blocks with upper left index between 0 and Nm - 1
-
-      // get parameters for rotation
-      s    = S(i, i+1);
-      lam1 = S(i, i);
-      lam2 = S(i+1, i+1);
-      phi  = s / std::abs(s);
-      r    = std::sqrt(std::pow(std::abs(s), 2) + std::pow(std::abs(lam2 - lam1), 2));
-
-      // compute Givens rotation corresponding to these parameters
-      Givens = CMat::Identity(Nm, Nm);
-      Givens(i, i)     = std::abs(s) / r;
-      Givens(i+1, i+1) = Givens(i, i);
-      Givens(i, i+1)   = (phi / r) * std::conj(lam2 - lam1);
-      Givens(i+1, i)   = -std::conj(Givens(i, i+1));
-      
-      // rotate Schur matrix and unitary change of basis matrix Q
-      S = Givens * S * Givens.adjoint();
-      Q = Givens * Q;
-
-      return;
-    }
-
-    /**
-     * Reorders a Schur matrix &Schur to have the Ritz values that we would like to keep for 
-     * restart as the first Nk elements on the diagonal. 
-     * 
-     * This algorithm is implemented as Nk iterations of a a reverse bubble sort with comparator compare.
-     * TODO: pass in compare function as an argument, default to compare with <. 
-     */
-    // void schurReorder(int Nk, std::function compare) {
-    void schurReorder(int Nk) {
-      for (int i = 0; i < Nk; i++) {
-        for (int k = 0; k <= Nm - 2; k++) {
-          int idx = Nm - 2 - k;
-          // TODO use RitzFilter enum here
-          // if (std::abs(S(idx, idx)) < std::abs(S(idx+1, idx+1))) {            // sort by largest modulus of eigenvalue
-          // if (std::real(S(idx, idx)) > std::real(S(idx+1, idx+1))) {        // sort by smallest real eigenvalue
-          if ( cCompare(S(idx+1, idx+1), S(idx, idx)) ) {            // sort by largest modulus of eigenvalue
-            swapEvals(idx);
-          }
-        }
-      }
-      return;
-    }
-
-    void schurReorderBlock() {
-      // TODO method stub
-      return;
-    }
-
-};
-
-// template<class Field>
-// inline void writeFile(const Field &field, const std::string &fname) {
-//   emptyUserRecord record;
-//   ScidacWriter WR(field.Grid()->IsBoss());
-//   WR.open(fname);
-//   WR.writeScidacFieldRecord(field, record, 0); // 0 = Lexico
-//   WR.close();
-// }
-
-/**
- * Implementation of the Krylov-Schur algorithm.
- */
-template<class Field> 
-class KrylovSchur {
-
-  private:
-  
-    std::string cname = std::string("KrylovSchur");
-    int MaxIter;   // Max iterations
-    RealD Tolerance;
-    RealD ssq;
-    RealD rtol;
-    int Nm;           // Number of basis vectors to track (equals MaxIter if no restart)
-    int Nk;           // Number of basis vectors to keep every restart (equals -1 if no restart)
-    int Nstop;       // Stop after converging Nstop eigenvectors.
-
-    LinearOperatorBase<Field> &Linop;
-    GridBase *Grid;
-
-    RealD approxLambdaMax;
-    RealD beta_k;
-    Field u;                                // Residual vector perpendicular to Krylov space (u_{k+1} in notes)
-    Eigen::VectorXcd   b;                   // b vector in Schur decomposition (e_{k+1} in Arnoldi).
-    std::vector<Field> basis;               // orthonormal Krylov basis
-    Eigen::MatrixXcd   Rayleigh;            // Rayleigh quotient of size Nbasis (after construction)
-    Eigen::MatrixXcd   Qt;                  // Transpose of basis rotation which projects out high modes.
-
-    Eigen::VectorXcd   evals;               // evals of Rayleigh quotient
-    std::vector<RealD> ritzEstimates;       // corresponding ritz estimates for evals
-    Eigen::MatrixXcd   littleEvecs;         // Nm x Nm evecs matrix
-    std::vector<Field> evecs;               // Vector of evec fields
-
-    RitzFilter ritzFilter;                  // how to sort evals
-
-  public:       
-    RealD *shift;                           // for Harmonic (shift and invert)
-
-    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),
-        evals (0), ritzEstimates (), evecs (), ssq (0.0), rtol (0.0), beta_k (0.0), approxLambdaMax (0.0),shift(NULL)
-    {
-      u = Zero();
-    };
-
-
-    /* Getters */
-    int                 getNk()                { return Nk;            }
-    Eigen::MatrixXcd    getRayleighQuotient()  { return Rayleigh;      }
-    Field               getU()                 { return u;             }
-    std::vector<Field>  getBasis()             { return basis;         }
-    Eigen::VectorXcd    getEvals()             { return evals;         }
-    std::vector<RealD>  getRitzEstimates()     { return ritzEstimates; }
-    std::vector<Field>  getEvecs()             { return evecs;         }
-
-    /**
-     * Runs the Krylov-Schur loop.
-     *   - Runs an Arnoldi step to generate the Rayleigh quotient and Krylov basis.
-     *   - Schur decompose the Rayleigh quotient. 
-     *   - Permutes the Rayleigh quotient according to the eigenvalues. 
-     *   - Truncate the Krylov-Schur expansion. 
-     */
-    void operator()(const Field& v0, int _maxIter, int _Nm, int _Nk, int _Nstop, bool doubleOrthog = true) {
-      RealD shift_=1.;
-      shift = &shift_;
-      if (shift) 
-          std::cout << GridLogMessage << "Shift " << *shift << std::endl;
-      MaxIter = _maxIter;
-      Nm = _Nm; Nk = _Nk;
-      Nstop = _Nstop;
-      
-      ssq = norm2(v0);
-      RealD approxLambdaMax = approxMaxEval(v0);
-      rtol = Tolerance * approxLambdaMax;
-      std::cout << GridLogMessage << "Approximate max eigenvalue: " << approxLambdaMax << std::endl;
-      // rtol = Tolerance;
-
-      b = Eigen::VectorXcd::Zero(Nm);       // start as e_{k+1}
-      b(Nm-1) = 1.0;
-
-      // basis = new std::vector<Field> (Nm, Grid);
-      // evecs.reserve();
-
-      int start = 0;
-      Field startVec = v0;
-      littleEvecs = Eigen::MatrixXcd::Zero(Nm, Nm);
-      for (int i = 0; i < MaxIter; i++) {
-        std::cout << GridLogMessage << "Restart Iteration " << i << std::endl;
-
-        // Perform Arnoldi steps to compute Krylov basis and Rayleigh quotient (Hess)
-        arnoldiIteration(startVec, Nm, start, doubleOrthog);
-        startVec = u;        // original code
-        start = Nk;
-
-        // checkKSDecomposition();
-
-        // Perform a Schur decomposition on Rayleigh
-        // ComplexSchurDecomposition schur (Rayleigh, false);
-	
-        Eigen::MatrixXcd temp = Rayleigh;
-	for (int m=0;m<Nm;m++) temp(m,m) -= *shift;
-        Eigen::MatrixXcd RayleighS = temp.inverse();
-        Eigen::MatrixXcd temp2 = RayleighS*temp;
-	std::cout << GridLogMessage  << "Shift inverse check: shift= "<<*shift<<" "<< temp2 <<std::endl;
-	
-
-        ComplexSchurDecomposition schur (Rayleigh, false, ritzFilter);
-        ComplexSchurDecomposition schurS (RayleighS, false, ritzFilter);
-        std::cout << GridLogDebug << "Schur decomp holds? " << schur.checkDecomposition() << std::endl;
-
-        // Rearrange Schur matrix so wanted evals are on top left (like MATLAB's ordschur)
-        std::cout << GridLogMessage << "Reordering Schur eigenvalues" << std::endl;
-        schur.schurReorder(Nk);
-        std::cout << GridLogMessage << "Shifted Schur eigenvalues" << std::endl;
-        schurS.schurReorder(Nk);
-        Eigen::MatrixXcd Q = schur.getMatrixQ();
-        Qt = Q.adjoint();                           // TODO should Q be real?
-        Eigen::MatrixXcd S = schur.getMatrixS();
-        // std::cout << GridLogDebug << "Schur decomp holds after reorder? " << schur.checkDecomposition() << std::endl;
-
-        std::cout << GridLogMessage << "*** ROTATING TO SCHUR BASIS *** " << std::endl;
-
-        // Rotate Krylov basis, b vector, redefine Rayleigh quotient and evecs, and truncate.
-        Rayleigh = schur.getMatrixS();
-        b = Q * b;            // b^\dag = b^\dag * Q^\dag <==> b = Q*b
-        // basisRotate(basis, Q, 0, Nm, 0, Nm, Nm);
-        // basisRotate(evecs, Q, 0, Nm, 0, Nm, Nm);
-
-        std::vector<Field> basis2; 
-        // basis2.reserve(Nm);
-        // for (int i = start; i < Nm; i++) {
-        //   basis2.emplace_back(Grid);
-        // }
-        constructUR(basis2, basis, Qt, Nm);
-        basis = basis2;
-
-        // std::vector<Field> evecs2;
-        // constructUR(evecs2, evecs, Qt, Nm);
-        // constructRU(evecs2, evecs, Q, Nm);
-        // evecs = evecs2;
-        // littleEvecs = littleEvecs * Q.adjoint();      // TODO try this and see if it works
-        // littleEvecs = Q * littleEvecs;      // TODO try this and see if it works
-        // std::cout << GridLogDebug << "Ritz vectors rotated correctly? " << checkEvecRotation() << std::endl;
-
-        // checkKSDecomposition();
-
-        std::cout << GridLogMessage << "*** TRUNCATING FOR RESTART *** " << std::endl;
-
-        std::cout << GridLogDebug << "Rayleigh before truncation: " << std::endl << Rayleigh << std::endl;
-
-        Eigen::MatrixXcd RayTmp = Rayleigh(Eigen::seqN(0, Nk), Eigen::seqN(0, Nk));
-        Rayleigh = RayTmp;
-
-        std::vector<Field> basisTmp = std::vector<Field> (basis.begin(), basis.begin() + Nk);
-        basis = basisTmp;
-
-        Eigen::VectorXcd btmp = b.head(Nk);
-        b = btmp;
-
-        std::cout << GridLogDebug << "Rayleigh after truncation: " << std::endl << Rayleigh << std::endl;
-
-        checkKSDecomposition();
-
-        // Compute eigensystem of Rayleigh. Note the eigenvectors correspond to the sorted eigenvalues.
-        computeEigensystem(Rayleigh);
-        std::cout << GridLogMessage << "Eigenvalues (first Nk sorted): " << std::endl << evals << std::endl;
-
-        // check convergence and return if needed. 
-        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 " 
-                        << 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;
-
-          // writeEigensystem(path);
-
-          return;
-        }
-      }
-    }
-
-    /**
-     * Constructs the Arnoldi basis for the Krylov space K_n(D, src). (TODO make private)
-     * 
-     * Parameters
-     * ----------
-     * v0 : Field&
-     *  Source to generate Krylov basis. 
-     * Nm : int
-     *  Final size of the basis desired. If the basis becomes complete before a basis of size Nm is constructed 
-     *  (determined by relative tolerance Tolerance), stops iteration there. 
-     * doubleOrthog : bool (default = false)
-     *  Whether to double orthogonalize the basis (for numerical cancellations) or not. 
-     * start        : int (default = 0)
-     *  If non-zero, assumes part of the Arnoldi basis has already been constructed. 
-     */
-    void arnoldiIteration(const Field& v0, int Nm, int start = 0, bool doubleOrthog = true)
-    {
-
-      ComplexD coeff;
-      Field w (Grid);           // A acting on last Krylov vector. 
-
-      // basis.reserve(Nm);
-      // for (int i = start; i < Nm; i++) {
-      //   basis.emplace_back(Grid);
-      // }
-      // basis.assign(Nm, Field(Grid));
-      // basis.resize(Nm);
-      // for (int i = start; i < Nm; i++) {
-      //   basis[i] = Field(Grid);
-      // }
-
-      if (start == 0) {       // initialize everything that we need.
-        RealD v0Norm = 1 / std::sqrt(ssq);
-        basis.push_back(v0Norm * v0);                // normalized source
-        // basis[0] = v0Norm * v0;                // normalized source
-
-        Rayleigh = Eigen::MatrixXcd::Zero(Nm, Nm); // CJ: B in SLEPc
-        u = Zero();                      
-      } else {
-        // assert( start == basis.size() );      // should be starting at the end of basis (start = Nk)
-        std::cout << GridLogMessage << "Resetting Rayleigh and b" << std::endl;
-        Eigen::MatrixXcd RayleighCp = Rayleigh;
-        Rayleigh = Eigen::MatrixXcd::Zero(Nm, Nm);
-        Rayleigh(Eigen::seqN(0, Nk), Eigen::seqN(0, Nk)) = RayleighCp;
-
-        // append b^\dag to Rayleigh, add u to basis
-        Rayleigh(Nk, Eigen::seqN(0, Nk)) = b.adjoint();
-        basis.push_back(u);
-        // basis[start] = u;       // TODO make sure this is correct
-        b = Eigen::VectorXcd::Zero(Nm);
-      }
-
-      // Construct next Arnoldi vector by normalizing w_i = Dv_i - \sum_j v_j h_{ji}
-      for (int i = start; i < Nm; i++) {
-        Linop.Op(basis.back(), w);
-        // Linop.Op(basis[i], w);
-        for (int j = 0; j < basis.size(); j++) {
-          coeff = innerProduct(basis[j], w);       // coeff = h_{ij}. Note that since {vi} is ONB it's OK to subtract it off after. 
-          Rayleigh(j, i) = coeff;
-          w -= coeff * basis[j];
-        }
-
-        if (doubleOrthog) {
-          std::cout << GridLogMessage << "Double orthogonalizing." << std::endl;
-          for (int j = 0; j < basis.size(); j++) {
-            coeff = innerProduct(basis[j], w);      // see if there is any residual component in basis[j] direction
-            Rayleigh(j, i) += coeff;                // if coeff is non-zero, adjust Rayleigh
-            w -= coeff * basis[j];
-          }
-        }
-
-        // add w_i to the pile
-        if (i < Nm - 1) {
-          coeff = std::sqrt(norm2(w));
-          Rayleigh(i+1, i) = coeff;
-          basis.push_back(
-            (1.0/coeff) * w
-          );
-          // basis[i+1] = (1.0/coeff) * w;
-        }
-
-        // after iterations, update u and beta_k = ||u|| before norm
-        u = w;                                // make sure u is normalized
-        beta_k = std::sqrt(norm2(u));         // beta_k = ||f_k|| determines convergence.
-        u = (1/beta_k) * u;
-      }
-
-      b(Nm - 1) = beta_k;
-
-      // std::cout << GridLogMessage << "|f|^2 after Arnoldi step = " << norm2(f) << std::endl;
-      std::cout << GridLogMessage << "beta_k = |u| (before norm) after Arnoldi step = " << beta_k << std::endl;
-      std::cout << GridLogDebug << "Computed Rayleigh quotient = " << std::endl << Rayleigh << std::endl;
-
-      return;
-    }
-
-    /**
-     * Approximates the eigensystem of the linear operator by computing the eigensystem of 
-     * the Rayleigh quotient. Assumes that the Rayleigh quotient has already been constructed (by 
-     * calling the operator() function).
-     * 
-     * Parameters
-     * ----------
-     * Eigen::MatrixXcd& S
-     *  Schur matrix (upper triangular) similar to original Rayleigh quotient.
-     */
-    void computeEigensystem(Eigen::MatrixXcd& S)
-    {
-      std::cout << GridLogMessage << "Computing eigenvalues." << std::endl;
-
-      // evals = S.diagonal();
-      int n = evals.size();       // should be regular Nm
-
-      evecs.clear();
-      // evecs.assign(n, Field(Grid));
-
-      // TODO: is there a faster way to get the eigenvectors of a triangular matrix?
-      // Rayleigh.triangularView<Eigen::Upper> tri;
-
-      Eigen::ComplexEigenSolver<Eigen::MatrixXcd> es;
-      // es.compute(Rayleigh);
-      es.compute(S);
-      evals = es.eigenvalues();
-      littleEvecs = es.eigenvectors();
-
-      // std::cout << GridLogDebug << "Little evecs: " << littleEvecs << std::endl;
-      // std::cout << "Rayleigh diag: " << S.diagonal() << std::endl;
-      // std::cout << "Rayleigh evals: " << evals << std::endl;
-
-      // Convert evecs to lattice fields
-      for (int k = 0; k < n; k++) {
-        Eigen::VectorXcd vec = littleEvecs.col(k);
-        Field tmp (basis[0].Grid());
-        tmp = Zero();
-        for (int j = 0; j < basis.size(); j++) {
-          tmp = tmp + vec[j] * basis[j];
-        }
-        evecs.push_back(tmp);
-        // evecs[k] = tmp;
-      }
-    }
-
-    /**
-     * Approximates the maximum eigenvalue of Linop.Op to normalize the residual and test for convergence. 
-     * 
-     * TODO implement in parent class eventually
-     * 
-     * Parameters
-     * ----------
-     * Field& v0
-     *  Source field to start with. Must have non-zero norm.
-     * int MAX_ITER (default = 50)
-     *  Maximum number of iterations for power approximation. 
-     * 
-     * Returns
-     * -------
-     * RealD lamApprox
-     *  Approximation of largest eigenvalue. 
-     */
-    RealD approxMaxEval(const Field& v0, int MAX_ITER = 50) {
-      assert (norm2(v0) > 1e-8);                        // must have relatively large source norm to start
-      RealD lamApprox = 0.0;
-      RealD denom = 1.0; RealD num = 1.0;
-      Field v0cp (Grid); Field tmp (Grid);
-      v0cp = v0;
-      denom = std::sqrt(norm2(v0cp));
-      for (int i = 0; i < MAX_ITER; i++) {
-        Linop.Op(v0cp, tmp);                               // CAREFUL: do not do Op(tmp, tmp)
-        v0cp = tmp;
-        num = std::sqrt(norm2(v0cp));                      // num = |A^{n+1} v0|
-        lamApprox = num / denom;                           // lam = |A^{n+1} v0| / |A^n v0|
-        std::cout << GridLogDebug << "Approx for max eval: " << lamApprox << std::endl;
-        denom = num;                                       // denom = |A^{n} v0|
-      }
-      return lamApprox;
-    }
-
-    /**
-     * Computes the number of Krylov-Schur eigenvectors that have converged. An eigenvector s is considered converged 
-     * for a tolerance epsilon if 
-     *    r(s) := |\beta e_m^T s| < epsilon
-     * where beta is the norm of f_{m+1}.
-     * 
-     * TODO implement in parent class eventually
-     * 
-     * Parameters
-     * ----------
-     * 
-     * Returns
-     * -------
-     * int : Number of converged eigenvectors.
-     */
-    int converged() {
-      int Nconv = 0;
-      int _Nm = evecs.size();
-      std::cout << GridLogDebug << "b: " << b << std::endl;
-      Field tmp (Grid); Field fullEvec (Grid);
-      ritzEstimates.clear();
-      // ritzEstimates.resize(_Nm);
-      for (int k = 0; k < _Nm; k++) {
-        Eigen::VectorXcd evec_k = littleEvecs.col(k);
-        RealD ritzEstimate = std::abs(b.dot(evec_k));           // b^\dagger s
-        ritzEstimates.push_back(ritzEstimate);
-        // ritzEstimates[k] = ritzEstimate;
-        std::cout << GridLogMessage << "Ritz estimate for evec " << k << " = " << ritzEstimate << std::endl;
-        if (ritzEstimate < rtol) {
-          Nconv++;
-        }
-
-      }
-      // Check that Ritz estimate is explicitly || D (Uy) - lambda (Uy) ||
-      // checkRitzEstimate();
-      return Nconv;
-    }
-
-    /**
-     * Checks the Krylov-Schur decomposition DU = UR + f b^\dag with the last-computed 
-     * U, R, f, and b. 
-     */
-    bool checkKSDecomposition(RealD tol = 1e-8) {
-
-      std::cout << GridLogMessage << "*** CHECKING KRYLOV-SCHUR DECOMPOSITION *** " << std::endl;
-
-      int k = basis.size();         // number of basis vectors, also the size of Rayleigh.
-
-      // rotate basis by Rayleigh to construct UR
-      // std::vector<Field> rotated;
-
-      // std::cout << GridLogDebug << "Rayleigh in KSDecomposition: " << std::endl << Rayleigh << std::endl;
-
-      std::vector<Field> rotated = basis;
-      constructUR(rotated, basis, Rayleigh, k);             // manually rotate
-      // Eigen::MatrixXcd Rt = Rayleigh.adjoint();
-      // basisRotate(rotated, Rt, 0, k, 0, k, k);           // UR
-
-      // TODO: make a new function that I'm positive does what this is doing
-      // just take the basis U = (u1 u2 ... uNm) and form the linear combination UR from R
-      
-      // For each i, form D u(i) and subtract off (US - u b^\dag)(i)
-      RealD delta = 0.0; RealD deltaSum = 0;
-      Field tmp (Grid); tmp = Zero();
-      for (int i = 0; i < k; i++) {
-        Linop.Op(basis[i], tmp);          // tmp = D u(i)
-        delta = norm2(tmp - rotated[i] - u * std::conj(b(i)));
-        delta = delta / norm2(tmp);               // relative tolerance
-        deltaSum += delta;
-
-        // std::cout << GridLogDebug << "Iteration " << i << std::endl;
-        // std::cout << GridLogDebug << "Du = " << norm2(tmp) << std::endl;
-        // std::cout << GridLogDebug << "rotated = " << norm2(rotated[i]) << std::endl;
-        // std::cout << GridLogDebug << "b[i] = " << b(i) << std::endl;
-        std::cout << GridLogMessage << "Deviation in decomp, column " << i << ": " << delta << std::endl;
-      }
-      std::cout << GridLogMessage << "Squared sum of relative deviations in decomposition: " << deltaSum << std::endl;
-
-      // std::cout << "[DEBUG] testing basis rotate" << std::endl;
-      // std::vector<Field> rotated2;
-      // constructUR(rotated2, basis, Rayleigh, k);
-      // for (int i = 0; i < k; i++) {
-      //   std::cout << "rotated[i] - UR[i] = " << norm2(rotated[i] - rotated2[i]) << std::endl;
-      // }
-
-      return deltaSum < tol;
-    }
-
-    /**
-     * Checks the Ritz vector s was rotated correctly by explicitly recomputing the 
-     * eigenvectors of the rotated Rayleigh quotient. 
-     */
-    bool checkRitzRotation(RealD tol = 1e-8) {
-      std::cout << GridLogMessage << "*** CHECKING RITZ VECTOR ROTATION *** " << std::endl;
-
-      Eigen::ComplexEigenSolver<Eigen::MatrixXcd> es;
-      es.compute(Rayleigh);
-      Eigen::MatrixXcd littleEvecs2 = es.eigenvectors();
-      RealD dLittle = (littleEvecs2 - littleEvecs).norm() / littleEvecs.norm();
-      std::cout << GridLogMessage << "|littleEvecs2 - littleEvecs| / |littleEvecs| = " << dLittle << std::endl;
-
-      std::cout << GridLogMessage << "Forming full eigenvectors" << std::endl;
-      RealD delta = 0.0; RealD deltaSum = 0;
-      for (int k = 0; k < evals.size(); k++) {
-        Eigen::VectorXcd vec = littleEvecs.col(k);
-        Field tmpEvec (Grid);
-        tmpEvec = Zero();
-        for (int j = 0; j < basis.size(); j++) {
-          tmpEvec = tmpEvec + vec[j] * basis[j];
-        }
-        delta = norm2(tmpEvec - evecs[k]) / norm2(evecs[k]);
-        std::cout << GridLogMessage << "Deviation in evec " << k << ": " << delta << std::endl;
-        deltaSum += delta;
-      }
-      return deltaSum < tol;
-    }
-
-    /**
-     * Checks the Ritz estimate R(s) is indeed the deviation of a Ritz eigenvector from being a true eigenvector. 
-     */
-    void checkRitzEstimate(RealD tol = 1e-8) {
-      std::cout << GridLogMessage << "*** CHECKING RITZ ESTIMATE *** " << std::endl;
-
-      // The issue was that the Eigen::eigensolver occasionally returned the complex conjugate pairs in the wrong 
-      // order compared to the diagonal, which is how I was reading them out. When this happened, the Ritz estimate would 
-      // be wrong. So, just need to be more careful and actually read out the eigenvalues. 
-
-      Field tmp (Grid);
-      // std::cout << "n evecs: " << evecs.size() << std::endl;
-      for (int k = 0; k < evecs.size(); k++) {
-        tmp = Zero();
-        Linop.Op(evecs[k], tmp);        // D evec
-        RealD ritz = std::sqrt(norm2(tmp - evals[k] * evecs[k]));
-        std::cout << "Ritz estimate " << k << " = " << ritz << std::endl;
-
-        // Checking little Ritz estimate
-        // Eigen::VectorXcd littleEvec = littleEvecs.col(k);
-        // Eigen::VectorXcd dev = Rayleigh * littleEvec - evals[k] * littleEvec;
-        // std::cout << GridLogMessage << "Little Ritz estimate = " << dev.norm() << std::endl;
-      }
-      return;
-    }
-
-    /**
-     * Given a vector of fields U (equivalently, a LxN matrix, where L is the number of degrees of 
-     * freedom on the lattice field) and an NxN matrix R, forms the product UR. 
-     * 
-     * Note that I believe this is equivalent to basisRotate(U, R.adjoint(), 0, N, 0, N, N), but I'm 
-     * not 100% sure (this will be slower and unoptimized though).
-     */
-    void constructUR(std::vector<Field>& UR, std::vector<Field> &U, Eigen::MatrixXcd& R, int N) {
-      Field tmp (Grid);
-
-      UR.clear();
-      // UR.resize(N);
-
-      std::cout << GridLogDebug << "R to rotate by (should be Rayleigh): " << R << std::endl;
-
-      for (int i = 0; i < N; i++) {
-        tmp = Zero();
-        for (int j = 0; j < N; j++) {
-          std::cout << GridLogDebug << "Adding R("<<j<<", "<<i<<") = " << R(j, i) << " to rotated" << std::endl;
-          std::cout << GridLogDebug << "Norm of U[j] is " << norm2(U[j]) << " to rotated" << std::endl;
-          tmp = tmp + U[j] * R(j, i);
-        }
-        std::cout << GridLogDebug << "rotated norm at i = " << i << " is: " << norm2(tmp) << std::endl;
-        UR.push_back(tmp);
-        // UR[i] = tmp;
-      }
-      return;
-    }
-
-    /**
-     * Same as constructUR but for the product order RU. 
-     */
-    void constructRU(std::vector<Field>& RU, std::vector<Field> &U, Eigen::MatrixXcd& R, int N) {
-      Field tmp (Grid);
-      RU.clear();
-      // RU.resize(N);
-      for (int i = 0; i < N; i++) {
-        tmp = Zero();
-        for (int j = 0; j < N; j++) {
-          tmp = tmp + R(i, j) * U[j];
-        }
-        RU.push_back(tmp);
-        // RU[i] = tmp;
-      }
-      return;
-    }
-
-    // void writeEvec(Field& in, std::string const fname){  
-    //   #ifdef HAVE_LIME
-    //     // Ref: https://github.com/paboyle/Grid/blob/feature/scidac-wp1/tests/debug/Test_general_coarse_hdcg_phys48.cc#L111
-    //     std::cout << GridLogMessage << "Writing evec to: " << fname << std::endl;
-    //     Grid::emptyUserRecord record;
-    //     Grid::ScidacWriter WR(in.Grid()->IsBoss());
-    //     WR.open(fname);
-    //     WR.writeScidacFieldRecord(in,record,0); // Lexico
-    //     WR.close();
-    //   #endif
-    //     // What is the appropriate way to throw error?
-    // }
-
-    // /**
-    //  * Writes the eigensystem of a Krylov Schur object to a directory. 
-    //  * 
-    //  * Parameters
-    //  * ----------
-    //  * std::string path
-    //  *    Directory to write to. 
-    //  */
-    // void writeEigensystem(std::string outDir) {
-    //   std::cout << GridLogMessage << "Writing output to directory: " << outDir << std::endl;
-    //   // TODO write a scidac density file so that we can easily integrate with visualization toolkit
-    //   std::string evalPath = outDir + "/evals.txt";
-    //   std::ofstream fEval;
-    //   fEval.open(evalPath);
-    //   for (int i = 0; i < Nk; i++) {
-    //     // write Eigenvalues
-    //     fEval << i << " " << evals(i);
-    //     if (i < Nk - 1) { fEval << "\n"; }
-    //   }
-    //   fEval.close();
-      
-    //   for (int i = 0; i < Nk; i++) {
-    //     std::string fName = outDir + "/evec" + std::to_string(i);
-    //     // writeFile(evecs[i], fName);     // using method from Grid/HMC/ComputeWilsonFlow.cc
-    //     // writeEvec(evecs[i], fName);
-    //   }
-      
-    // }
-
-};
-
-NAMESPACE_END(Grid);
-#endif
diff --git a/examples/Example_krylov_schur.cc b/examples/Example_krylov_schur.cc
deleted file mode 100644
index 54e3c9ca..00000000
--- a/examples/Example_krylov_schur.cc
+++ /dev/null
@@ -1,441 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./tests/Test_padded_cell.cc
-
-    Copyright (C) 2023
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-    /*  END LEGAL */
-
-// copied here from Test_general_coarse_pvdagm.cc
-
-#include <cstdlib>
-
-#include <Grid/Grid.h>
-#include <Grid/lattice/PaddedCell.h>
-#include <Grid/stencil/GeneralLocalStencil.h>
-
-#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
-#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
-#include <Grid/algorithms/iterative/BiCGSTAB.h>
-
-using namespace std;
-using namespace Grid;
-
-// Hermitize a DWF operator by squaring it
-template<class Matrix,class Field>
-class SquaredLinearOperator : public LinearOperatorBase<Field> {
-
-  public:
-  Matrix &_Mat;
-
-  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){};
-
-  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);
-  }
-  void AdjOp     (const Field &in, Field &out){
-    std::cout << "AdjOp: Mdag PV "<<std::endl;
-    Field tmp(in.Grid());
-    _PV.M(in,tmp);
-    _Mat.Mdag(tmp,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 PV PVdag M"<<std::endl;
-    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;
-  }
-};
-
-template<class Field>
-void testSchurFromHess(Arnoldi<Field>& Arn, Field& src, int Nlarge, int Nm, int Nk) {
-
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout << GridLogMessage << "Testing Schur reordering, Nm = " << Nm << ", Nk = " << Nk << std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-
-  std::cout << GridLogMessage << "Running Arnoldi for 1 iteration to get a Hessenberg." << std::endl;
-  Arn(src, 1, Nlarge, Nm, Nlarge);
-  Eigen::MatrixXcd Hess = Arn.getHessenbergMat();
-  std::cout << GridLogMessage << "Hessenberg for use: " << std::endl << Hess << std::endl;
-
-  ComplexSchurDecomposition schur (Hess, true);
-  bool isDecomposed = schur.checkDecomposition();
-  std::cout << "Schur decomp holds? " << isDecomposed << std::endl;
-
-  std::cout << GridLogMessage << "S = " << std::endl << schur.getMatrixS() << std::endl;
-  std::cout << GridLogMessage << "Swapping S(3, 3) with S(4, 4)" << std::endl;
-  schur.swapEvals(3);
-  std::cout << GridLogMessage << "S after swap = " << std::endl << schur.getMatrixS() << std::endl;
-  std::cout << "Schur decomp still holds? " << schur.checkDecomposition() << std::endl;
-
-  // Now move last diagonal element all the way to the front.
-  std::cout << GridLogMessage << "Moving last eval to front. S at start = " << std::endl << schur.getMatrixS() << std::endl;
-  for (int i = 0; i < Nk - 1; i++) {
-    int swapIdx = Nk - 2 - i;
-    schur.swapEvals(swapIdx);
-    std::cout << GridLogMessage << "S after swap of index " << swapIdx << " = " << std::endl << schur.getMatrixS() << std::endl;
-    std::cout << "Schur decomp still holds? " << schur.checkDecomposition() << std::endl;
-  }
-
-  std::cout << GridLogMessage << "Testing Schur reorder" << std::endl;
-  schur.schurReorder(Nk);
-  std::cout << GridLogMessage << "S after reorder = " << std::endl << schur.getMatrixS() << std::endl;
-  std::cout << "Schur decomp still holds? " << schur.checkDecomposition() << std::endl;
-
-}
-
-int main (int argc, char ** argv)
-{
-  Grid_init(&argc,&argv);
-
-  const int Ls=16;
-
-//   GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
-  std::vector<int> lat_size {16, 16, 16, 16};
-  std::cout << "Lattice size: " << lat_size << std::endl;
-  GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(lat_size, 
-								          GridDefaultSimd(Nd,vComplex::Nsimd()),
-								          GridDefaultMpi());
-  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
-
-  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
-  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
-
-  // Construct a coarsened grid
-  // poare TODO: replace this with the following line?
-  Coordinate clatt = lat_size;
-//   Coordinate clatt = GridDefaultLatt();              // [PO] initial line before I edited it
-  for(int d=0;d<clatt.size();d++){
-    clatt[d] = clatt[d]/2;
-    //    clatt[d] = clatt[d]/4;
-  }
-  GridCartesian *Coarse4d =  SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
-  GridCartesian *Coarse5d =  SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
-
-  std::vector<int> seeds4({1,2,3,4});
-  std::vector<int> seeds5({5,6,7,8});
-  std::vector<int> cseeds({5,6,7,8});
-  GridParallelRNG          RNG5(FGrid);   RNG5.SeedFixedIntegers(seeds5);
-  GridParallelRNG          RNG4(UGrid);   RNG4.SeedFixedIntegers(seeds4);
-  GridParallelRNG          CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
-
-  LatticeFermion    src(FGrid); random(RNG5,src);
-  LatticeFermion result(FGrid); result=Zero();
-  LatticeFermion    ref(FGrid); ref=Zero();
-  LatticeFermion    tmp(FGrid);
-  LatticeFermion    err(FGrid);
-  LatticeGaugeField Umu(UGrid);
-
-  FieldMetaData header;
-  std::string file("config");
-//  std::string file("Users/patrickoare/libraries/PETSc-Grid/ckpoint_lat.4000");
-  NerscIO::readConfiguration(Umu,header,file);
-
-  RealD mass=0.01;
-  RealD M5=1.8;
-
-  DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-  DomainWallFermionD Dpv(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0,M5);
-
-  // const int nbasis = 20;            // size of approximate basis for low-mode space
-  const int nbasis = 3;            // size of approximate basis for low-mode space
-  const int cb = 0 ;
-  LatticeFermion prom(FGrid);
-
-  typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
-  typedef LittleDiracOperator::CoarseVector CoarseVector;
-
-  NextToNearestStencilGeometry5D geom(Coarse5d);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-
-  typedef PVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> PVdagM_t;
-  typedef ShiftedPVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> ShiftedPVdagM_t;
-  typedef ShiftedComplexPVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> ShiftedComplexPVdagM_t;
-  PVdagM_t PVdagM(Ddwf, Dpv);
-  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv);
-  SquaredLinearOperator<DomainWallFermionD, LatticeFermionD> Dsq (Ddwf);
-  NonHermitianLinearOperator<DomainWallFermionD, LatticeFermionD> DLinOp (Ddwf);
-
-  // PowerMethod<LatticeFermion> PM; PM(PVdagM, src);
-  int Nm = 10;
-  int Nk = 4; 
-  // int Nk = Nm+1;     // if just running once
-  // int maxIter = 5;
-  // int maxIter = 1;
-//  int maxIter = 5;
-  int maxIter = 100;
-  int Nstop = 4;
-
-  Coordinate origin ({0,0,0,0});
-  auto tmpSrc = peekSite(src, origin);
-  std::cout << "[DEBUG] Source at origin = " <<  tmpSrc << std::endl;
-  LatticeFermion src2 = src;
-
-  // Run KrylovSchur and Arnoldi on a Hermitian matrix
-  std::cout << GridLogMessage << "Runnning Krylov Schur" << std::endl;
-  // KrylovSchur KrySchur (Dsq, FGrid, 1e-8, EvalNormLarge);
-  KrylovSchur KrySchur (Dsq, FGrid, 1e-8,EvalImNormSmall);
-  KrySchur(src, maxIter, Nm, Nk, Nstop);
-
-  /*
-  std::cout << GridLogMessage << "Running Arnoldi" << std::endl;
-  // Arnoldi Arn (Dsq, FGrid, 1e-8);
-  Arnoldi Arn (DLinOp, FGrid, 1e-8);
-  testSchurFromHess<LatticeFermion>(Arn, src, 10, 6, 4);
-
-  Arnoldi Arn2 (DLinOp, FGrid, 1e-8);
-  testSchurFromHess<LatticeFermion>(Arn2, src, 16, 12, 8);
-  */
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage << "Done "<< std::endl;
-
-  Grid_finalize();
-  return 0;
-}
diff --git a/examples/Example_spec_kryschur.cc b/examples/Example_spec_kryschur.cc
index f87d32be..7e70a180 100644
--- a/examples/Example_spec_kryschur.cc
+++ b/examples/Example_spec_kryschur.cc
@@ -65,6 +65,52 @@
 using namespace std;
 using namespace Grid;
 
+<<<<<<< HEAD
+namespace Grid {
+
+struct LanczosParameters: Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(LanczosParameters,
+		  		RealD, mass , 
+		  		RealD, mstep , 
+				Integer, Nstop,
+                                Integer, Nk,
+                                Integer, Np,
+                                Integer, ReadEvec,
+	  			RealD, resid,
+	  			RealD, ChebyLow,
+	  			RealD, ChebyHigh,
+	  			Integer, ChebyOrder)
+
+  LanczosParameters() {
+    /////////////////////////////////
+  }
+
+  template <class ReaderClass >
+  LanczosParameters(Reader<ReaderClass> & TheReader){
+    initialize(TheReader);
+  }
+
+  template < class ReaderClass > 
+  void initialize(Reader<ReaderClass> &TheReader){
+    read(TheReader, "HMC", *this);
+  }
+
+
+  void print_parameters() const {
+//    std::cout << GridLogMessage << "[HMC parameters] Trajectories            : " << Trajectories << "\n";
+//    std::cout << GridLogMessage << "[HMC parameters] Start trajectory        : " << StartTrajectory << "\n";
+//    std::cout << GridLogMessage << "[HMC parameters] Metropolis test (on/off): " << std::boolalpha << MetropolisTest << "\n";
+//    std::cout << GridLogMessage << "[HMC parameters] Thermalization trajs    : " << NoMetropolisUntil << "\n";
+//    std::cout << GridLogMessage << "[HMC parameters] Starting type           : " << StartingType << "\n";
+//    MD.print_parameters();
+  }
+  
+};
+
+}
+
+#if 0
+=======
 template <class T> void writeFile(T& in, std::string const fname){  
   #ifdef HAVE_LIME
     // Ref: https://github.com/paboyle/Grid/blob/feature/scidac-wp1/tests/debug/Test_general_coarse_hdcg_phys48.cc#L111
@@ -111,6 +157,7 @@ void writeEigensystem(KrylovSchur<Field> KS, std::string outDir) {
   // }
 }
 
+>>>>>>> 68af1bba67dd62881ead5ab1e54962a5486a0791
 // Hermitize a DWF operator by squaring it
 template<class Matrix,class Field>
 class SquaredLinearOperator : public LinearOperatorBase<Field> {
@@ -256,6 +303,111 @@ ShiftedComplexPVdagMLinearOperator(ComplexD _shift,Matrix &Mat,Matrix &PV): shif
   }
 };
 
+<<<<<<< HEAD
+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
+
+=======
+>>>>>>> 68af1bba67dd62881ead5ab1e54962a5486a0791
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);