diff --git a/Grid/algorithms/iterative/KrylovSchur.h b/Grid/algorithms/iterative/KrylovSchur.h
index 0c660585..9906741b 100644
--- a/Grid/algorithms/iterative/KrylovSchur.h
+++ b/Grid/algorithms/iterative/KrylovSchur.h
@@ -32,6 +32,12 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #ifndef GRID_KRYLOVSCHUR_H
 #define GRID_KRYLOVSCHUR_H
 
+// #include <Grid/Grid.h>
+// #include <Grid/parallelIO/IldgIOtypes.h>
+// #include <Grid/parallelIO/IldgIO.h>
+
+// #include <Grid/serialisation/emptyUserRecord.h>
+
 NAMESPACE_BEGIN(Grid); 
 
 /**
@@ -167,6 +173,15 @@ class ComplexSchurDecomposition {
 
 };
 
+// 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.
  */
@@ -212,6 +227,7 @@ class KrylovSchur {
 
 
     /* Getters */
+    int                 getNk()                { return Nk;       }
     Eigen::MatrixXcd    getRayleighQuotient()  { return Rayleigh; }
     Field               getU()                 { return u;        }
     std::vector<Field>  getBasis()             { return basis;    }
@@ -321,6 +337,9 @@ class KrylovSchur {
                         << i << "." << std::endl;
           basisRotate(evecs, Qt, 0, Nk, 0, Nk, Nm);
           std::cout << GridLogMessage << "Eigenvalues: " << evals << std::endl;
+
+          // writeEigensystem(path);
+
           return;
         }
       }
@@ -413,8 +432,6 @@ class KrylovSchur {
      * the Rayleigh quotient. Assumes that the Rayleigh quotient has already been constructed (by 
      * calling the operator() function).
      * 
-     * TODO implement in parent class eventually.
-     * 
      * Parameters
      * ----------
      * Eigen::MatrixXcd& S
@@ -653,7 +670,49 @@ class KrylovSchur {
       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_spec_kryschur.cc b/examples/Example_spec_kryschur.cc
index 6b3160ca..61c6ab30 100644
--- a/examples/Example_spec_kryschur.cc
+++ b/examples/Example_spec_kryschur.cc
@@ -40,9 +40,57 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
 #include <Grid/algorithms/iterative/BiCGSTAB.h>
 
+#include <Grid/parallelIO/IldgIOtypes.h>
+#include <Grid/parallelIO/IldgIO.h>
+
 using namespace std;
 using namespace Grid;
 
+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
+    std::cout << Grid::GridLogMessage << "Writes 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
+}
+
+/**
+ * Writes the eigensystem of a Krylov Schur object to a directory. 
+ * 
+ * Parameters
+ * ----------
+ * std::string path
+ *    Directory to write to. 
+ */
+template <class Field>
+void writeEigensystem(KrylovSchur<Field> KS, std::string outDir) {
+  int Nk = KS.getNk();
+  std::cout << GridLogMessage << "Writing output to directory: " << outDir << std::endl;
+  
+  // Write evals
+  std::string evalPath = outDir + "/evals.txt";
+  std::ofstream fEval;
+  fEval.open(evalPath);
+  Eigen::VectorXcd evals = KS.getEvals();
+  for (int i = 0; i < Nk; i++) {
+    // write Eigenvalues
+    fEval << i << " " << evals(i);
+    if (i < Nk - 1) { fEval << "\n"; }
+  }
+  fEval.close();
+  
+  // Write evecs
+  std::vector<Field> evecs = KS.getEvecs();
+  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
+  }
+}
+
 // Hermitize a DWF operator by squaring it
 template<class Matrix,class Field>
 class SquaredLinearOperator : public LinearOperatorBase<Field> {
@@ -293,12 +341,14 @@ int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  // Usage : $ ./Example_spec_kryschur <Nm> <Nk> <maxiter> <Nstop>
+  // Usage : $ ./Example_spec_kryschur <Nm> <Nk> <maxiter> <Nstop> <inFile> <outDir>
   // assert (argc == 5);
-  std::string NmStr = argv[1];
-  std::string NkStr = argv[2];
+  std::string NmStr      = argv[1];
+  std::string NkStr      = argv[2];
   std::string maxIterStr = argv[3];
-  std::string NstopStr = argv[4];
+  std::string NstopStr   = argv[4];
+  std::string file       = argv[5];
+  std::string outDir     = argv[6];
 
   //const int Ls=16;
   const int Ls = 8;
@@ -342,10 +392,11 @@ int main (int argc, char ** argv)
 
   FieldMetaData header;
   // std::string file ("/sdcc/u/poare/PETSc-Grid/ckpoint_EODWF_lat.125");
-  std::string file("/Users/patrickoare/libraries/PETSc-Grid/ckpoint_EODWF_lat.125");
+  // std::string file("/Users/patrickoare/libraries/PETSc-Grid/ckpoint_EODWF_lat.125");
   NerscIO::readConfiguration(Umu,header,file);
 
-  RealD mass=0.01;
+  // RealD mass=0.01;
+  RealD mass=0.001;
   RealD M5=1.8;
 
   DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
@@ -397,6 +448,10 @@ int main (int argc, char ** argv)
   std::cout << GridLogMessage << "Krylov Schur eigenvalues: " << std::endl << KrySchur.getEvals() << std::endl;
   //std::cout << GridLogMessage << "Lanczos eigenvalues: " << std::endl << levals << std::endl;
 
+  // KrySchur.writeEigensystem(outDir);
+
+  writeEigensystem(KrySchur, outDir);
+
   std::cout<<GridLogMessage<<std::endl;
   std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
   std::cout<<GridLogMessage<<std::endl;