diff --git a/lib/algorithms/iterative/DenseMatrix.h b/lib/algorithms/iterative/DenseMatrix.h
new file mode 100644
index 00000000..2423677d
--- /dev/null
+++ b/lib/algorithms/iterative/DenseMatrix.h
@@ -0,0 +1,106 @@
+#ifndef GRID_DENSE_MATRIX_H
+#define GRID_DENSE_MATRIX_H
+
+namespace Grid {
+    /////////////////////////////////////////////////////////////
+    // Matrix untils
+    /////////////////////////////////////////////////////////////
+
+template<class T> using DenseVector = std::vector<T>;
+template<class T> using DenseMatrix = DenseVector<DenseVector<T> >;
+
+template<class T> void Size(DenseVector<T> & vec, int &N) 
+{ 
+  N= vec.size();
+}
+template<class T> void Size(DenseMatrix<T> & mat, int &N,int &M) 
+{ 
+  N= mat.size();
+  M= mat[0].size();
+}
+
+template<class T> void SizeSquare(DenseMatrix<T> & mat, int &N) 
+{ 
+  int M; Size(mat,N,M);
+  assert(N==M);
+}
+
+template<class T> void Resize(DenseMatrix<T > & mat, int N, int M) { 
+  mat.resize(N);
+  for(int i=0;i<N;i++){
+    mat[i].resize(M);
+  }
+}
+template<class T> void Fill(DenseMatrix<T> & mat, T&val) { 
+  int N,M;
+  Size(mat,N,M);
+  for(int i=0;i<N;i++){
+  for(int j=0;j<M;j++){
+    mat[i][j] = val;
+  }}
+}
+
+/** Transpose of a matrix **/
+template<class T> DenseMatrix<T> Transpose(DenseMatrix<T> & mat){
+  int N,M;
+  Size(mat,N,M);
+  DenseMatrix<T> C; Resize(C,M,N);
+  for(int i=0;i<M;i++){
+  for(int j=0;j<N;j++){
+    C[i][j] = mat[j][i];
+  }} 
+  return C;
+}
+/** Set DenseMatrix to unit matrix **/
+template<class T> void Unity(DenseMatrix<T> &A){
+  int N;  SizeSquare(A,N);
+  for(int i=0;i<N;i++){
+    for(int j=0;j<N;j++){
+      if ( i==j ) A[i][j] = 1;
+      else        A[i][j] = 0;
+    } 
+  } 
+}
+
+/** Add C * I to matrix **/
+template<class T>
+void PlusUnit(DenseMatrix<T> & A,T c){
+  int dim;  SizeSquare(A,dim);
+  for(int i=0;i<dim;i++){A[i][i] = A[i][i] + c;} 
+}
+
+/** return the Hermitian conjugate of matrix **/
+template<class T>
+DenseMatrix<T> HermitianConj(DenseMatrix<T> &mat){
+
+  int dim; SizeSquare(mat,dim);
+
+  DenseMatrix<T> C; Resize(C,dim,dim);
+
+  for(int i=0;i<dim;i++){
+    for(int j=0;j<dim;j++){
+      C[i][j] = conj(mat[j][i]);
+    } 
+  } 
+  return C;
+}
+/**Get a square submatrix**/
+template <class T>
+DenseMatrix<T> GetSubMtx(DenseMatrix<T> &A,int row_st, int row_end, int col_st, int col_end)
+{
+  DenseMatrix<T> H; Resize(H,row_end - row_st,col_end-col_st);
+
+  for(int i = row_st; i<row_end; i++){
+  for(int j = col_st; j<col_end; j++){
+    H[i-row_st][j-col_st]=A[i][j];
+  }}
+  return H;
+}
+
+}
+
+#include <algorithms/iterative/Householder.h>
+#include <algorithms/iterative/Francis.h>
+
+#endif
+
diff --git a/lib/algorithms/iterative/EigenSort.h b/lib/algorithms/iterative/EigenSort.h
new file mode 100644
index 00000000..c036110b
--- /dev/null
+++ b/lib/algorithms/iterative/EigenSort.h
@@ -0,0 +1,52 @@
+#ifndef GRID_EIGENSORT_H
+#define GRID_EIGENSORT_H
+
+
+namespace Grid {
+    /////////////////////////////////////////////////////////////
+    // Eigen sorter to begin with
+    /////////////////////////////////////////////////////////////
+
+template<class Field>
+class SortEigen {
+ private:
+  
+  static bool less_lmd(RealD left,RealD right){
+    return fabs(left) < fabs(right);
+  }  
+  static bool less_pair(std::pair<RealD,Field>& left,
+		 std::pair<RealD,Field>& right){
+    return fabs(left.first) < fabs(right.first);
+  }  
+  
+ public:
+
+  void push(DenseVector<RealD>& lmd,
+	    DenseVector<Field>& evec,int N) {
+
+    DenseVector<std::pair<RealD, Field> > emod;
+    typename DenseVector<std::pair<RealD, Field> >::iterator it;
+    
+    for(int i=0;i<lmd.size();++i){
+      emod.push_back(std::pair<RealD,Field>(lmd[i],evec[i]));
+    }
+
+    partial_sort(emod.begin(),emod.begin()+N,emod.end(),less_pair);
+
+    it=emod.begin();
+    for(int i=0;i<N;++i){
+      lmd[i]=it->first;
+      evec[i]=it->second;
+      ++it;
+    }
+  }
+  void push(DenseVector<RealD>& lmd,int N) {
+    std::partial_sort(lmd.begin(),lmd.begin()+N,lmd.end(),less_lmd);
+  }
+  bool saturated(RealD lmd, RealD thrs) {
+    return fabs(lmd) > fabs(thrs);
+  }
+};
+
+}
+#endif
diff --git a/lib/algorithms/iterative/Francis.h b/lib/algorithms/iterative/Francis.h
new file mode 100644
index 00000000..405193d2
--- /dev/null
+++ b/lib/algorithms/iterative/Francis.h
@@ -0,0 +1,498 @@
+#ifndef FRANCIS_H
+#define FRANCIS_H
+
+#include <cstdlib>
+#include <string>
+#include <cmath>
+#include <iostream>
+#include <sstream>
+#include <stdexcept>
+#include <fstream>
+#include <complex>
+#include <algorithm>
+
+//#include <timer.h>
+//#include <lapacke.h>
+//#include <Eigen/Dense>
+
+namespace Grid {
+
+template <class T> int SymmEigensystem(DenseMatrix<T > &Ain, DenseVector<T> &evals, DenseMatrix<T> &evecs, RealD small);
+template <class T> int     Eigensystem(DenseMatrix<T > &Ain, DenseVector<T> &evals, DenseMatrix<T> &evecs, RealD small);
+
+/**
+  Find the eigenvalues of an upper hessenberg matrix using the Francis QR algorithm.
+H =
+      x  x  x  x  x  x  x  x  x
+      x  x  x  x  x  x  x  x  x
+      0  x  x  x  x  x  x  x  x
+      0  0  x  x  x  x  x  x  x
+      0  0  0  x  x  x  x  x  x
+      0  0  0  0  x  x  x  x  x
+      0  0  0  0  0  x  x  x  x
+      0  0  0  0  0  0  x  x  x
+      0  0  0  0  0  0  0  x  x
+Factorization is P T P^H where T is upper triangular (mod cc blocks) and P is orthagonal/unitary.
+**/
+template <class T>
+int QReigensystem(DenseMatrix<T> &Hin, DenseVector<T> &evals, DenseMatrix<T> &evecs, RealD small)
+{
+  DenseMatrix<T> H = Hin; 
+
+  int N ; SizeSquare(H,N);
+  int M = N;
+
+  Fill(evals,0);
+  Fill(evecs,0);
+
+  T s,t,x=0,y=0,z=0;
+  T u,d;
+  T apd,amd,bc;
+  DenseVector<T> p(N,0);
+  T nrm = Norm(H);    ///DenseMatrix Norm
+  int n, m;
+  int e = 0;
+  int it = 0;
+  int tot_it = 0;
+  int l = 0;
+  int r = 0;
+  DenseMatrix<T> P; Resize(P,N,N); Unity(P);
+  DenseVector<int> trows(N,0);
+
+  /// Check if the matrix is really hessenberg, if not abort
+  RealD sth = 0;
+  for(int j=0;j<N;j++){
+    for(int i=j+2;i<N;i++){
+      sth = abs(H[i][j]);
+      if(sth > small){
+	std::cout << "Non hessenberg H = " << sth << " > " << small << std::endl;
+	exit(1);
+      }
+    }
+  }
+
+  do{
+    std::cout << "Francis QR Step N = " << N << std::endl;
+    /** Check for convergence
+      x  x  x  x  x
+      0  x  x  x  x
+      0  0  x  x  x
+      0  0  x  x  x
+      0  0  0  0  x
+      for this matrix l = 4
+     **/
+    do{
+      l = Chop_subdiag(H,nrm,e,small);
+      r = 0;    ///May have converged on more than one eval
+      ///Single eval
+      if(l == N-1){
+        evals[e] = H[l][l];
+        N--; e++; r++; it = 0;
+      }
+      ///RealD eval
+      if(l == N-2){
+        trows[l+1] = 1;    ///Needed for UTSolve
+        apd = H[l][l] + H[l+1][l+1];
+        amd = H[l][l] - H[l+1][l+1];
+        bc =  (T)4.0*H[l+1][l]*H[l][l+1];
+        evals[e]   = (T)0.5*( apd + sqrt(amd*amd + bc) );
+        evals[e+1] = (T)0.5*( apd - sqrt(amd*amd + bc) );
+        N-=2; e+=2; r++; it = 0;
+      }
+    } while(r>0);
+
+    if(N ==0) break;
+
+    DenseVector<T > ck; Resize(ck,3);
+    DenseVector<T> v;   Resize(v,3);
+
+    for(int m = N-3; m >= l; m--){
+      ///Starting vector essentially random shift.
+      if(it%10 == 0 && N >= 3 && it > 0){
+        s = (T)1.618033989*( abs( H[N-1][N-2] ) + abs( H[N-2][N-3] ) );
+        t = (T)0.618033989*( abs( H[N-1][N-2] ) + abs( H[N-2][N-3] ) );
+        x = H[m][m]*H[m][m] + H[m][m+1]*H[m+1][m] - s*H[m][m] + t;
+        y = H[m+1][m]*(H[m][m] + H[m+1][m+1] - s);
+        z = H[m+1][m]*H[m+2][m+1];
+      }
+      ///Starting vector implicit Q theorem
+      else{
+        s = (H[N-2][N-2] + H[N-1][N-1]);
+        t = (H[N-2][N-2]*H[N-1][N-1] - H[N-2][N-1]*H[N-1][N-2]);
+        x = H[m][m]*H[m][m] + H[m][m+1]*H[m+1][m] - s*H[m][m] + t;
+        y = H[m+1][m]*(H[m][m] + H[m+1][m+1] - s);
+        z = H[m+1][m]*H[m+2][m+1];
+      }
+      ck[0] = x; ck[1] = y; ck[2] = z;
+
+      if(m == l) break;
+
+      /** Some stupid thing from numerical recipies, seems to work**/
+      // PAB.. for heaven's sake quote page, purpose, evidence it works.
+      //       what sort of comment is that!?!?!?
+      u=abs(H[m][m-1])*(abs(y)+abs(z));
+      d=abs(x)*(abs(H[m-1][m-1])+abs(H[m][m])+abs(H[m+1][m+1]));
+      if ((T)abs(u+d) == (T)abs(d) ){
+	l = m; break;
+      }
+
+      //if (u < small){l = m; break;}
+    }
+    if(it > 100000){
+     std::cout << "QReigensystem: bugger it got stuck after 100000 iterations" << std::endl;
+     std::cout << "got " << e << " evals " << l << " " << N << std::endl;
+      exit(1);
+    }
+    normalize(ck);    ///Normalization cancels in PHP anyway
+    T beta;
+    Householder_vector<T >(ck, 0, 2, v, beta);
+    Householder_mult<T >(H,v,beta,0,l,l+2,0);
+    Householder_mult<T >(H,v,beta,0,l,l+2,1);
+    ///Accumulate eigenvector
+    Householder_mult<T >(P,v,beta,0,l,l+2,1);
+    int sw = 0;      ///Are we on the last row?
+    for(int k=l;k<N-2;k++){
+      x = H[k+1][k];
+      y = H[k+2][k];
+      z = (T)0.0;
+      if(k+3 <= N-1){
+	z = H[k+3][k];
+      } else{
+	sw = 1; 
+	v[2] = (T)0.0;
+      }
+      ck[0] = x; ck[1] = y; ck[2] = z;
+      normalize(ck);
+      Householder_vector<T >(ck, 0, 2-sw, v, beta);
+      Householder_mult<T >(H,v, beta,0,k+1,k+3-sw,0);
+      Householder_mult<T >(H,v, beta,0,k+1,k+3-sw,1);
+      ///Accumulate eigenvector
+      Householder_mult<T >(P,v, beta,0,k+1,k+3-sw,1);
+    }
+    it++;
+    tot_it++;
+  }while(N > 1);
+  N = evals.size();
+  ///Annoying - UT solves in reverse order;
+  DenseVector<T> tmp; Resize(tmp,N);
+  for(int i=0;i<N;i++){
+    tmp[i] = evals[N-i-1];
+  } 
+  evals = tmp;
+  UTeigenvectors(H, trows, evals, evecs);
+  for(int i=0;i<evals.size();i++){evecs[i] = P*evecs[i]; normalize(evecs[i]);}
+  return tot_it;
+}
+
+template <class T>
+int my_Wilkinson(DenseMatrix<T> &Hin, DenseVector<T> &evals, DenseMatrix<T> &evecs, RealD small)
+{
+  /**
+  Find the eigenvalues of an upper Hessenberg matrix using the Wilkinson QR algorithm.
+  H =
+  x  x  0  0  0  0
+  x  x  x  0  0  0
+  0  x  x  x  0  0
+  0  0  x  x  x  0
+  0  0  0  x  x  x
+  0  0  0  0  x  x
+  Factorization is P T P^H where T is upper triangular (mod cc blocks) and P is orthagonal/unitary.  **/
+  return my_Wilkinson(Hin, evals, evecs, small, small);
+}
+
+template <class T>
+int my_Wilkinson(DenseMatrix<T> &Hin, DenseVector<T> &evals, DenseMatrix<T> &evecs, RealD small, RealD tol)
+{
+  int N; SizeSquare(Hin,N);
+  int M = N;
+
+  ///I don't want to modify the input but matricies must be passed by reference
+  //Scale a matrix by its "norm"
+  //RealD Hnorm = abs( Hin.LargestDiag() ); H =  H*(1.0/Hnorm);
+  DenseMatrix<T> H;  H = Hin;
+  
+  RealD Hnorm = abs(Norm(Hin));
+  H = H * (1.0 / Hnorm);
+
+  // TODO use openmp and memset
+  Fill(evals,0);
+  Fill(evecs,0);
+
+  T s, t, x = 0, y = 0, z = 0;
+  T u, d;
+  T apd, amd, bc;
+  DenseVector<T> p; Resize(p,N); Fill(p,0);
+
+  T nrm = Norm(H);    ///DenseMatrix Norm
+  int n, m;
+  int e = 0;
+  int it = 0;
+  int tot_it = 0;
+  int l = 0;
+  int r = 0;
+  DenseMatrix<T> P; Resize(P,N,N);
+  Unity(P);
+  DenseVector<int> trows(N, 0);
+  /// Check if the matrix is really symm tridiag
+  RealD sth = 0;
+  for(int j = 0; j < N; ++j)
+  {
+    for(int i = j + 2; i < N; ++i)
+    {
+      if(abs(H[i][j]) > tol || abs(H[j][i]) > tol)
+      {
+	std::cout << "Non Tridiagonal H(" << i << ","<< j << ") = |" << Real( real( H[j][i] ) ) << "| > " << tol << std::endl;
+	std::cout << "Warning tridiagonalize and call again" << std::endl;
+        // exit(1); // see what is going on
+        //return;
+      }
+    }
+  }
+
+  do{
+    do{
+      //Jasper
+      //Check if the subdiagonal term is small enough (<small)
+      //if true then it is converged.
+      //check start from H.dim - e - 1
+      //How to deal with more than 2 are converged?
+      //What if Chop_symm_subdiag return something int the middle?
+      //--------------
+      l = Chop_symm_subdiag(H,nrm, e, small);
+      r = 0;    ///May have converged on more than one eval
+      //Jasper
+      //In this case
+      // x  x  0  0  0  0
+      // x  x  x  0  0  0
+      // 0  x  x  x  0  0
+      // 0  0  x  x  x  0
+      // 0  0  0  x  x  0
+      // 0  0  0  0  0  x  <- l
+      //--------------
+      ///Single eval
+      if(l == N - 1)
+      {
+        evals[e] = H[l][l];
+        N--;
+        e++;
+        r++;
+        it = 0;
+      }
+      //Jasper
+      // x  x  0  0  0  0
+      // x  x  x  0  0  0
+      // 0  x  x  x  0  0
+      // 0  0  x  x  0  0
+      // 0  0  0  0  x  x  <- l
+      // 0  0  0  0  x  x
+      //--------------
+      ///RealD eval
+      if(l == N - 2)
+      {
+        trows[l + 1] = 1;    ///Needed for UTSolve
+        apd = H[l][l] + H[l + 1][ l + 1];
+        amd = H[l][l] - H[l + 1][l + 1];
+        bc =  (T) 4.0 * H[l + 1][l] * H[l][l + 1];
+        evals[e] = (T) 0.5 * (apd + sqrt(amd * amd + bc));
+        evals[e + 1] = (T) 0.5 * (apd - sqrt(amd * amd + bc));
+        N -= 2;
+        e += 2;
+        r++;
+        it = 0;
+      }
+    }while(r > 0);
+    //Jasper
+    //Already converged
+    //--------------
+    if(N == 0) break;
+
+    DenseVector<T> ck,v; Resize(ck,2); Resize(v,2);
+
+    for(int m = N - 3; m >= l; m--)
+    {
+      ///Starting vector essentially random shift.
+      if(it%10 == 0 && N >= 3 && it > 0)
+      {
+        t = abs(H[N - 1][N - 2]) + abs(H[N - 2][N - 3]);
+        x = H[m][m] - t;
+        z = H[m + 1][m];
+      } else {
+      ///Starting vector implicit Q theorem
+        d = (H[N - 2][N - 2] - H[N - 1][N - 1]) * (T) 0.5;
+        t =  H[N - 1][N - 1] - H[N - 1][N - 2] * H[N - 1][N - 2] 
+	  / (d + sign(d) * sqrt(d * d + H[N - 1][N - 2] * H[N - 1][N - 2]));
+        x = H[m][m] - t;
+        z = H[m + 1][m];
+      }
+      //Jasper
+      //why it is here????
+      //-----------------------
+      if(m == l)
+        break;
+
+      u = abs(H[m][m - 1]) * (abs(y) + abs(z));
+      d = abs(x) * (abs(H[m - 1][m - 1]) + abs(H[m][m]) + abs(H[m + 1][m + 1]));
+      if ((T)abs(u + d) == (T)abs(d))
+      {
+        l = m;
+        break;
+      }
+    }
+    //Jasper
+    if(it > 1000000)
+    {
+      std::cout << "Wilkinson: bugger it got stuck after 100000 iterations" << std::endl;
+      std::cout << "got " << e << " evals " << l << " " << N << std::endl;
+      exit(1);
+    }
+    //
+    T s, c;
+    Givens_calc<T>(x, z, c, s);
+    Givens_mult<T>(H, l, l + 1, c, -s, 0);
+    Givens_mult<T>(H, l, l + 1, c,  s, 1);
+    Givens_mult<T>(P, l, l + 1, c,  s, 1);
+    //
+    for(int k = l; k < N - 2; ++k)
+    {
+      x = H.A[k + 1][k];
+      z = H.A[k + 2][k];
+      Givens_calc<T>(x, z, c, s);
+      Givens_mult<T>(H, k + 1, k + 2, c, -s, 0);
+      Givens_mult<T>(H, k + 1, k + 2, c,  s, 1);
+      Givens_mult<T>(P, k + 1, k + 2, c,  s, 1);
+    }
+    it++;
+    tot_it++;
+  }while(N > 1);
+
+  N = evals.size();
+  ///Annoying - UT solves in reverse order;
+  DenseVector<T> tmp(N);
+  for(int i = 0; i < N; ++i)
+    tmp[i] = evals[N-i-1];
+  evals = tmp;
+  //
+  UTeigenvectors(H, trows, evals, evecs);
+  //UTSymmEigenvectors(H, trows, evals, evecs);
+  for(int i = 0; i < evals.size(); ++i)
+  {
+    evecs[i] = P * evecs[i];
+    normalize(evecs[i]);
+    evals[i] = evals[i] * Hnorm;
+  }
+  // // FIXME this is to test
+  // Hin.write("evecs3", evecs);
+  // Hin.write("evals3", evals);
+  // // check rsd
+  // for(int i = 0; i < M; i++) {
+  //   vector<T> Aevec = Hin * evecs[i];
+  //   RealD norm2(0.);
+  //   for(int j = 0; j < M; j++) {
+  //     norm2 += (Aevec[j] - evals[i] * evecs[i][j]) * (Aevec[j] - evals[i] * evecs[i][j]);
+  //   }
+  // }
+  return tot_it;
+}
+
+template <class T>
+void Hess(DenseMatrix<T > &A, DenseMatrix<T> &Q, int start){
+
+  /**
+  turn a matrix A =
+  x  x  x  x  x
+  x  x  x  x  x
+  x  x  x  x  x
+  x  x  x  x  x
+  x  x  x  x  x
+  into
+  x  x  x  x  x
+  x  x  x  x  x
+  0  x  x  x  x
+  0  0  x  x  x
+  0  0  0  x  x
+  with householder rotations
+  Slow.
+  */
+  int N ; SizeSquare(A,N);
+  DenseVector<T > p; Resize(p,N); Fill(p,0);
+
+  for(int k=start;k<N-2;k++){
+    //cerr << "hess" << k << std::endl;
+    DenseVector<T > ck,v; Resize(ck,N-k-1); Resize(v,N-k-1);
+    for(int i=k+1;i<N;i++){ck[i-k-1] = A(i,k);}  ///kth column
+    normalize(ck);    ///Normalization cancels in PHP anyway
+    T beta;
+    Householder_vector<T >(ck, 0, ck.size()-1, v, beta);  ///Householder vector
+    Householder_mult<T>(A,v,beta,start,k+1,N-1,0);  ///A -> PA
+    Householder_mult<T >(A,v,beta,start,k+1,N-1,1);  ///PA -> PAP^H
+    ///Accumulate eigenvector
+    Householder_mult<T >(Q,v,beta,start,k+1,N-1,1);  ///Q -> QP^H
+  }
+  /*for(int l=0;l<N-2;l++){
+    for(int k=l+2;k<N;k++){
+    A(0,k,l);
+    }
+    }*/
+}
+
+template <class T>
+void Tri(DenseMatrix<T > &A, DenseMatrix<T> &Q, int start){
+///Tridiagonalize a matrix
+  int N; SizeSquare(A,N);
+  Hess(A,Q,start);
+  /*for(int l=0;l<N-2;l++){
+    for(int k=l+2;k<N;k++){
+    A(0,l,k);
+    }
+    }*/
+}
+
+template <class T>
+void ForceTridiagonal(DenseMatrix<T> &A){
+///Tridiagonalize a matrix
+  int N ; SizeSquare(A,N);
+  for(int l=0;l<N-2;l++){
+    for(int k=l+2;k<N;k++){
+      A[l][k]=0;
+      A[k][l]=0;
+    }
+  }
+}
+
+template <class T>
+int my_SymmEigensystem(DenseMatrix<T > &Ain, DenseVector<T> &evals, DenseVector<DenseVector<T> > &evecs, RealD small){
+  ///Solve a symmetric eigensystem, not necessarily in tridiagonal form
+  int N; SizeSquare(Ain,N);
+  DenseMatrix<T > A; A = Ain;
+  DenseMatrix<T > Q; Resize(Q,N,N); Unity(Q);
+  Tri(A,Q,0);
+  int it = my_Wilkinson<T>(A, evals, evecs, small);
+  for(int k=0;k<N;k++){evecs[k] = Q*evecs[k];}
+  return it;
+}
+
+
+template <class T>
+int Wilkinson(DenseMatrix<T> &Ain, DenseVector<T> &evals, DenseVector<DenseVector<T> > &evecs, RealD small){
+  return my_Wilkinson(Ain, evals, evecs, small);
+}
+
+template <class T>
+int SymmEigensystem(DenseMatrix<T> &Ain, DenseVector<T> &evals, DenseVector<DenseVector<T> > &evecs, RealD small){
+  return my_SymmEigensystem(Ain, evals, evecs, small);
+}
+
+template <class T>
+int Eigensystem(DenseMatrix<T > &Ain, DenseVector<T> &evals, DenseVector<DenseVector<T> > &evecs, RealD small){
+///Solve a general eigensystem, not necessarily in tridiagonal form
+  int N = Ain.dim;
+  DenseMatrix<T > A(N); A = Ain;
+  DenseMatrix<T > Q(N);Q.Unity();
+  Hess(A,Q,0);
+  int it = QReigensystem<T>(A, evals, evecs, small);
+  for(int k=0;k<N;k++){evecs[k] = Q*evecs[k];}
+  return it;
+}
+
+}
+#endif
diff --git a/lib/algorithms/iterative/Householder.h b/lib/algorithms/iterative/Householder.h
new file mode 100644
index 00000000..01cf3453
--- /dev/null
+++ b/lib/algorithms/iterative/Householder.h
@@ -0,0 +1,215 @@
+#ifndef HOUSEHOLDER_H
+#define HOUSEHOLDER_H
+
+#define TIMER(A) std::cout << GridLogMessage << __FUNC__ << " file "<< __FILE__ <<" line " << __LINE__ << std::endl;
+#define ENTER()  std::cout << GridLogMessage << "ENTRY "<<__FUNC__ << " file "<< __FILE__ <<" line " << __LINE__ << std::endl;
+#define LEAVE()  std::cout << GridLogMessage << "EXIT  "<<__FUNC__ << " file "<< __FILE__ <<" line " << __LINE__ << std::endl;
+
+#include <cstdlib>
+#include <string>
+#include <cmath>
+#include <iostream>
+#include <sstream>
+#include <stdexcept>
+#include <fstream>
+#include <complex>
+#include <algorithm>
+
+namespace Grid {
+/** Comparison function for finding the max element in a vector **/
+template <class T> bool cf(T i, T j) { 
+  return abs(i) < abs(j); 
+}
+
+/** 
+	Calculate a real Givens angle 
+ **/
+template <class T> inline void Givens_calc(T y, T z, T &c, T &s){
+
+  RealD mz = (RealD)abs(z);
+  
+  if(mz==0.0){
+    c = 1; s = 0;
+  }
+  if(mz >= (RealD)abs(y)){
+    T t = -y/z;
+    s = (T)1.0 / sqrt ((T)1.0 + t * t);
+    c = s * t;
+  } else {
+    T t = -z/y;
+    c = (T)1.0 / sqrt ((T)1.0 + t * t);
+    s = c * t;
+  }
+}
+
+template <class T> inline void Givens_mult(DenseMatrix<T> &A,  int i, int k, T c, T s, int dir)
+{
+  int q ; SizeSquare(A,q);
+
+  if(dir == 0){
+    for(int j=0;j<q;j++){
+      T nu = A[i][j];
+      T w  = A[k][j];
+      A[i][j] = (c*nu + s*w);
+      A[k][j] = (-s*nu + c*w);
+    }
+  }
+
+  if(dir == 1){
+    for(int j=0;j<q;j++){
+      T nu = A[j][i];
+      T w  = A[j][k];
+      A[j][i] = (c*nu - s*w);
+      A[j][k] = (s*nu + c*w);
+    }
+  }
+}
+
+/**
+	from input = x;
+	Compute the complex Householder vector, v, such that
+	P = (I - b v transpose(v) )
+	b = 2/v.v
+
+	P | x |    | x | k = 0
+	| x |    | 0 | 
+	| x | =  | 0 |
+	| x |    | 0 | j = 3
+	| x |	   | x |
+
+	These are the "Unreduced" Householder vectors.
+
+ **/
+template <class T> inline void Householder_vector(DenseVector<T> input, int k, int j, DenseVector<T> &v, T &beta)
+{
+  int N ; Size(input,N);
+  T m = *max_element(input.begin() + k, input.begin() + j + 1, cf<T> );
+
+  if(abs(m) > 0.0){
+    T alpha = 0;
+
+    for(int i=k; i<j+1; i++){
+      v[i] = input[i]/m;
+      alpha = alpha + v[i]*conj(v[i]);
+    }
+    alpha = sqrt(alpha);
+    beta = (T)1.0/(alpha*(alpha + abs(v[k]) ));
+
+    if(abs(v[k]) > 0.0)  v[k] = v[k] + (v[k]/abs(v[k]))*alpha;
+    else                 v[k] = -alpha;
+  } else{
+    for(int i=k; i<j+1; i++){
+      v[i] = 0.0;
+    } 
+  }
+}
+
+/**
+	from input = x;
+	Compute the complex Householder vector, v, such that
+	P = (I - b v transpose(v) )
+	b = 2/v.v
+
+	Px = alpha*e_dir
+
+	These are the "Unreduced" Householder vectors.
+
+ **/
+
+template <class T> inline void Householder_vector(DenseVector<T> input, int k, int j, int dir, DenseVector<T> &v, T &beta)
+{
+  int N = input.size();
+  T m = *max_element(input.begin() + k, input.begin() + j + 1, cf);
+  
+  if(abs(m) > 0.0){
+    T alpha = 0;
+
+    for(int i=k; i<j+1; i++){
+      v[i] = input[i]/m;
+      alpha = alpha + v[i]*conj(v[i]);
+    }
+    
+    alpha = sqrt(alpha);
+    beta = 1.0/(alpha*(alpha + abs(v[dir]) ));
+	
+    if(abs(v[dir]) > 0.0) v[dir] = v[dir] + (v[dir]/abs(v[dir]))*alpha;
+    else                  v[dir] = -alpha;
+  }else{
+    for(int i=k; i<j+1; i++){
+      v[i] = 0.0;
+    } 
+  }
+}
+
+/**
+	Compute the product PA if trans = 0
+	AP if trans = 1
+	P = (I - b v transpose(v) )
+	b = 2/v.v
+	start at element l of matrix A
+	v is of length j - k + 1 of v are nonzero
+ **/
+
+template <class T> inline void Householder_mult(DenseMatrix<T> &A , DenseVector<T> v, T beta, int l, int k, int j, int trans)
+{
+  int N ; SizeSquare(A,N);
+
+  if(abs(beta) > 0.0){
+    for(int p=l; p<N; p++){
+      T s = 0;
+      if(trans==0){
+	for(int i=k;i<j+1;i++) s += conj(v[i-k])*A[i][p];
+	s *= beta;
+	for(int i=k;i<j+1;i++){ A[i][p] = A[i][p]-s*conj(v[i-k]);}
+      } else {
+	for(int i=k;i<j+1;i++){ s += conj(v[i-k])*A[p][i];}
+	s *= beta;
+	for(int i=k;i<j+1;i++){ A[p][i]=A[p][i]-s*conj(v[i-k]);}
+      }
+    }
+  }
+}
+
+/**
+	Compute the product PA if trans = 0
+	AP if trans = 1
+	P = (I - b v transpose(v) )
+	b = 2/v.v
+	start at element l of matrix A
+	v is of length j - k + 1 of v are nonzero
+	A is tridiagonal
+ **/
+template <class T> inline void Householder_mult_tri(DenseMatrix<T> &A , DenseVector<T> v, T beta, int l, int M, int k, int j, int trans)
+{
+  if(abs(beta) > 0.0){
+
+    int N ; SizeSquare(A,N);
+
+    DenseMatrix<T> tmp; Resize(tmp,N,N); Fill(tmp,0); 
+
+    T s;
+    for(int p=l; p<M; p++){
+      s = 0;
+      if(trans==0){
+	for(int i=k;i<j+1;i++) s = s + conj(v[i-k])*A[i][p];
+      }else{
+	for(int i=k;i<j+1;i++) s = s + v[i-k]*A[p][i];
+      }
+      s = beta*s;
+      if(trans==0){
+	for(int i=k;i<j+1;i++) tmp[i][p] = tmp(i,p) - s*v[i-k];
+      }else{
+	for(int i=k;i<j+1;i++) tmp[p][i] = tmp[p][i] - s*conj(v[i-k]);
+      }
+    }
+    for(int p=l; p<M; p++){
+      if(trans==0){
+	for(int i=k;i<j+1;i++) A[i][p] = A[i][p] + tmp[i][p];
+      }else{
+	for(int i=k;i<j+1;i++) A[p][i] = A[p][i] + tmp[p][i];
+      }
+    }
+  }
+}
+}
+#endif
diff --git a/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h b/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h
index 9d6153da..c71cbdcd 100644
--- a/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h
+++ b/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h
@@ -1,105 +1,144 @@
 #ifndef GRID_IRL_H
 #define GRID_IRL_H
 
+#include <algorithms/iterative/DenseMatrix.h>
+#include <algorithms/iterative/EigenSort.h>
+
 namespace Grid {
 
-    /////////////////////////////////////////////////////////////
-    // Base classes for iterative processes based on operators
-    // single input vec, single output vec.
-    /////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////
+// Implicitly restarted lanczos
+/////////////////////////////////////////////////////////////
 
-  template<class Field> 
+
+template<class Field> 
     class ImplicitlyRestartedLanczos {
-public:       
 
+    const RealD small = 1.0e-16;
+public:       
+    int lock;
+    int converged;
     int Niter;
-    int Nk;
-    int Np;
-    RealD enorm;
-    RealD vthr;
-    
+
+    int Nk;      // Number of converged sought
+    int Np;      // Np -- Number of spare vecs in kryloc space
+    int Nm;      // Nm -- total number of vectors
+
+    RealD eresid;
+
+    SortEigen<Field> _sort;
+
     LinearOperatorBase<Field> &_Linop;
+
     OperatorFunction<Field>   &_poly;
 
-    ImplicitlyRestartedLanczos(
-			       LinearOperatorBase<Field> &Linop,
-			       OperatorFunction<Field> & poly,
-			       int _Nk,
-			       int _Np,
-			       RealD _enorm,
-			       RealD _vthrs,
-			       int _Niter) :
-    _Linop(Linop),
+    /////////////////////////
+    // Constructor
+    /////////////////////////
+    void init(void){};
+    void Abort(int ff, DenseVector<RealD> &evals,  DenseVector<DenseVector<RealD> > &evecs);
+
+    ImplicitlyRestartedLanczos(LinearOperatorBase<Field> &Linop, // op
+			       OperatorFunction<Field> & poly,   // polynmial
+			       int _Nk, // sought vecs
+			       int _Nm, // spare vecs
+			       RealD _eresid, // resid in lmdue deficit 
+			       int _Niter) : // Max iterations
+      _Linop(Linop),
       _poly(poly),
       Nk(_Nk),
-      Np(_Np),
-      enorm(_enorm),
-      vthr(_vthrs)
+      Nm(_Nm),
+      eresid(_eresid),
+      Niter(_Niter)
     { 
-      vthr=_vthrs;
-      Niter=_Niter;
+      Np = Nm-Nk; assert(Np>0);
     };
 
+    /////////////////////////
+    // Sanity checked this routine (step) against Saad.
+    /////////////////////////
+    void RitzMatrix(DenseVector<Field>& evec,int k){
+      if(1) return;
 
-    void step(Vector<RealD>& lmda,
-	      Vector<RealD>& lmdb, 
-	      Vector<Field>& evec,
-	      Field& f,int Nm,int k)
-    {
-
-      assert( k< Nm );
-      
-      w = opr_->mult(evec[k]);
-
-      if(k==0){  // Initial step
-    
-	RealD wnorm= w*w;
-	std::cout<<"wnorm="<<wnorm<<std::endl;
-
-	RealD alph = evec[k] * w;
-	w -= alph * evec[k];
-	lmd[k] = alph;
-
-	RealD beta = w * w;
-	beta = sqrt(beta);
-	RealD betar = 1.0/beta;
-
-	evec[k+1] = betar * w;
-	lme[k] = beta;
-
-      } else {   // Iteration step
-
-	w -= lme[k-1] * evec[k-1];
-	
-	RealD alph = evec[k] * w;
-	w -= alph * evec[k];
-
-	RealD beta = w * w;
-	beta = sqrt(beta);
-	RealD betar = 1.0/beta;
-	w *= betar;
-	
-	lmd[k] = alph;
-	lme[k] = beta;
-	orthogonalize(w,evec,k);
-	
-	if(k < Nm-1) evec[k+1] = w;
+      GridBase *grid = evec[0]._grid;
+      Field w(grid);
+      std::cout << "RitzMatrix "<<std::endl;
+      for(int i=0;i<k;i++){
+	_poly(_Linop,evec[i],w);
+	std::cout << "["<<i<<"] ";
+	for(int j=0;j<k;j++){
+	  ComplexD in = innerProduct(evec[j],w);
+	  if ( fabs((double)i-j)>1 ) { 
+	    if (abs(in) >1.0e-9 )  { 
+	      std::cout<<"oops"<<std::endl;
+	      abort();
+	    } else 
+	      std::cout << " 0 ";
+	  } else { 
+	    std::cout << " "<<in<<" ";
+	  }
+	}
+	std::cout << std::endl;
       }
     }
 
-    void qr_decomp(Vector<RealD>& lmda,
-		   Vector<RealD>& lmdb,
+/* Saad PP. 195
+1. Choose an initial vector v1 of 2-norm unity. Set β1 ≡ 0, v0 ≡ 0
+2. For k = 1,2,...,m Do:
+3. wk:=Avk−βkv_{k−1}      
+4. αk:=(wk,vk)       // 
+5. wk:=wk−αkvk       // wk orthog vk 
+6. βk+1 := ∥wk∥2. If βk+1 = 0 then Stop
+7. vk+1 := wk/βk+1
+8. EndDo
+ */
+    void step(DenseVector<RealD>& lmd,
+	      DenseVector<RealD>& lme, 
+	      DenseVector<Field>& evec,
+	      Field& w,int Nm,int k)
+    {
+      assert( k< Nm );
+      
+      _poly(_Linop,evec[k],w);      // 3. wk:=Avk−βkv_{k−1}
+      if(k>0){
+	w -= lme[k-1] * evec[k-1];
+      }    
+
+      ComplexD zalph = innerProduct(evec[k],w); // 4. αk:=(wk,vk)
+      RealD     alph = real(zalph);
+
+      w = w - alph * evec[k];// 5. wk:=wk−αkvk
+
+      RealD beta = normalise(w); // 6. βk+1 := ∥wk∥2. If βk+1 = 0 then Stop
+                                 // 7. vk+1 := wk/βk+1
+
+      const RealD tiny = 1.0e-20;
+      if ( beta < tiny ) { 
+	std::cout << " beta is tiny "<<beta<<std::endl;
+      }
+      lmd[k] = alph;
+      lme[k]  = beta;
+
+      if (k>0) { 
+	orthogonalize(w,evec,k); // orthonormalise
+      }
+      
+      if(k < Nm-1) evec[k+1] = w;
+    }
+
+    void qr_decomp(DenseVector<RealD>& lmd,
+		   DenseVector<RealD>& lme,
 		   int Nk,
 		   int Nm,
-		   Vector<RealD>& Qt,
-		   RealD Dsft, 
+		   DenseVector<RealD>& Qt,
+		   RealD Dsh, 
 		   int kmin,
 		   int kmax)
     {
       int k = kmin-1;
       RealD x;
 
-      RealD Fden = 1.0/sqrt((lmd[k]-Dsh)*(lmd[k]-Dsh) +lme[k]*lme[k]);
+      RealD Fden = 1.0/hypot(lmd[k]-Dsh,lme[k]);
       RealD c = ( lmd[k] -Dsh) *Fden;
       RealD s = -lme[k] *Fden;
       
@@ -110,7 +149,7 @@ public:
       lmd[k]   = c*c*tmpa1 +s*s*tmpa2 -2.0*c*s*tmpb;
       lmd[k+1] = s*s*tmpa1 +c*c*tmpa2 +2.0*c*s*tmpb;
       lme[k]   = c*s*(tmpa1-tmpa2) +(c*c-s*s)*tmpb;
-      x        = -s*lme[k+1];
+      x        =-s*lme[k+1];
       lme[k+1] = c*lme[k+1];
       
       for(int i=0; i<Nk; ++i){
@@ -122,7 +161,8 @@ public:
 
       // Givens transformations
       for(int k = kmin; k < kmax-1; ++k){
-	RealD Fden = 1.0/sqrt( x*x +lme[k-1]*lme[k-1]);
+
+	RealD Fden = 1.0/hypot(x,lme[k-1]);
 	RealD c = lme[k-1]*Fden;
 	RealD s = - x*Fden;
 	
@@ -149,11 +189,11 @@ public:
       }
     }
 
-    void diagonalize(Vector<RealD>& lmda,
-		     Vector<RealD>& lmdb, 
+    void diagonalize(DenseVector<RealD>& lmd,
+		     DenseVector<RealD>& lme, 
 		     int Nm2,
 		     int Nm,
-		     Vector<RealD>& Qt)
+		     DenseVector<RealD>& Qt)
     {
       int Niter = 100*Nm;
       int kmin = 1;
@@ -195,102 +235,134 @@ public:
       abort();
     }
 
-    void orthogonalize(Field& w,
-		       const Vector<Field>& evec,
-		       int k)
+    static RealD normalise(Field& v) 
     {
-  // Schmidt orthogonalization                                   
-                                                
-      size_t size = w.size();
-      assert(size%2 ==0);
-
-      std::slice re(0,size/2,2);
-      std::slice im(1,size/2,2);
-
-      for(int j=0; j<k; ++j){
-	RealD prdr = evec[j]*w;
-	RealD prdi = evec[j].im_prod(w);
-	
-	valarray<RealD> evr(evec[j][re]);
-	valarray<RealD> evi(evec[j][im]);
-
-	w.add(re, -prdr*evr +prdi*evi);
-	w.add(im, -prdr*evi -prdi*evr);
-      }
+      RealD nn = norm2(v);
+      nn = sqrt(nn);
+      v = v * (1.0/nn);
+      return nn;
     }
 
-    void calc(Vector<RealD>& lmd,
-	      Vector<Field>& evec,
-	      const Field& b,
-	      int& Nsbt,
+    void orthogonalize(Field& w,
+		       DenseVector<Field>& evec,
+		       int k)
+    {
+      typedef typename Field::scalar_type MyComplex;
+      MyComplex ip;
+
+      if ( 0 ) {
+	for(int j=0; j<k; ++j){
+	  normalise(evec[j]);
+	  for(int i=0;i<j;i++){
+	    ip = innerProduct(evec[i],evec[j]); // are the evecs normalised? ; this assumes so.
+	    evec[j] = evec[j] - ip *evec[i];
+	  }
+	}
+      }
+
+      for(int j=0; j<k; ++j){
+	ip = innerProduct(evec[j],w); // are the evecs normalised? ; this assumes so.
+	w = w - ip * evec[j];
+      }
+      normalise(w);
+    }
+
+    void setUnit_Qt(int Nm, DenseVector<RealD> &Qt) {
+      for(int i=0; i<Qt.size(); ++i) Qt[i] = 0.0;
+      for(int k=0; k<Nm; ++k) Qt[k + k*Nm] = 1.0;
+    }
+
+/* Rudy Arthur's thesis pp.137
+------------------------
+Require: M > K P = M − K †
+Compute the factorization AVM = VM HM + fM eM 
+repeat
+  Q=I
+  for i = 1,...,P do
+    QiRi =HM −θiI Q = QQi
+    H M = Q †i H M Q i
+  end for
+  βK =HM(K+1,K) σK =Q(M,K)
+  r=vK+1βK +rσK
+  VK =VM(1:M)Q(1:M,1:K)
+  HK =HM(1:K,1:K)
+  →AVK =VKHK +fKe†K † Extend to an M = K + P step factorization AVM = VMHM + fMeM
+until convergence
+ */
+    void calc(DenseVector<RealD>& eval,
+	      DenseVector<Field>& evec,
+	      const Field& src,
 	      int& Nconv)
       {
-	const size_t fsize = evec[0].size();
-  
-	Nconv = -1;
-	Nsbt = 0;
-	int Nm = Nk_+Np_;
 
-	std::cout << " -- Nk = " << Nk_ << " Np = "<< Np_ << endl;
-	std::cout << " -- Nm = " << Nm << endl;
-	std::cout << " -- size of lmd   = " << lmd.size() << endl;
-	std::cout << " -- size of evec  = " << evec.size() << endl;
+	GridBase *grid = evec[0]._grid;
+
+	std::cout << " -- Nk = " << Nk << " Np = "<< Np << std::endl;
+	std::cout << " -- Nm = " << Nm << std::endl;
+	std::cout << " -- size of eval   = " << eval.size() << std::endl;
+	std::cout << " -- size of evec  = " << evec.size() << std::endl;
 	
-	assert(Nm < evec.size() && Nm < lmd.size());
+	assert(Nm == evec.size() && Nm == eval.size());
 	
-	vector<RealD> lme(Nm);
-	vector<RealD> lmd2(Nm);
-	vector<RealD> lme2(Nm);
-	vector<RealD> Qt(Nm*Nm);
-	vector<int>    Iconv(Nm);
+	DenseVector<RealD> lme(Nm);  
+	DenseVector<RealD> lme2(Nm);
+	DenseVector<RealD> eval2(Nm);
+	DenseVector<RealD> Qt(Nm*Nm);
+	DenseVector<int>   Iconv(Nm);
+
+	DenseVector<Field>  B(Nm,grid); // waste of space replicating
 	
-	vector<Field>  B(Nm);
-	for(int k=0; k<Nm; ++k) B[k].resize(fsize);
-	
-	Field f(fsize);
-	Field v(fsize);
+	Field f(grid);
+	Field v(grid);
   
 	int k1 = 1;
-	int k2 = Nk_;
-	int kconv = 0;
-	
-	int Kdis  = 0;
-	int Kthrs = 0;
+	int k2 = Nk;
+
+	Nconv = 0;
+
 	RealD beta_k;
   
 	// Set initial vector
-	evec[0] = 1.0;
-	RealD vnorm = evec[0]*evec[0];
-	evec[0] = 1.0/sqrt(vnorm);
-	// (uniform vector)
+	// (uniform vector) Why not src??
+	//	evec[0] = 1.0;
+	evec[0] = src;
+	normalise(evec[0]);
 	
 	// Initial Nk steps
-	for(int k=0; k<k2; ++k) step(lmd,lme,evec,f,Nm,k);
-	
-	// Restarting loop begins
-	for(int iter = 0; iter<Niter_; ++iter){
-	  std::cout<<"\n iteration = "<< iter << endl;
+	for(int k=0; k<Nk; ++k) step(eval,lme,evec,f,Nm,k);
+	RitzMatrix(evec,Nk);
 
-	  int Nm2 = Nm - kconv;
-	  for(int k=k2; k<Nm; ++k) step(lmd,lme,evec,f,Nm,k);
+	// Restarting loop begins
+	for(int iter = 0; iter<Niter; ++iter){
+
+	  std::cout<<"\n Restart iteration = "<< iter << std::endl;
+
+	  // 
+	  // Rudy does a sort first which looks very different. Getting fed up with sorting out the algo defs.
+	  // We loop over 
+	  //
+	  for(int k=Nk; k<Nm; ++k) step(eval,lme,evec,f,Nm,k);
 	  f *= lme[Nm-1];
+
+	  RitzMatrix(evec,k2);
 	  
 	  // getting eigenvalues
-	  for(int k=0; k<Nm2; ++k){
-	    lmd2[k] = lmd[k+k1-1];
+	  for(int k=0; k<Nm; ++k){
+	    eval2[k] = eval[k+k1-1];
 	    lme2[k] = lme[k+k1-1];
 	  }
 	  setUnit_Qt(Nm,Qt);
-	  diagonalize(lmd2,lme2,Nm2,Nm,Qt);
+	  diagonalize(eval2,lme2,Nm,Nm,Qt);
+
 	  // sorting
-	  sort_->push(lmd2,Nm);
+	  _sort.push(eval2,Nm);
 	  
 	  // Implicitly shifted QR transformations
 	  setUnit_Qt(Nm,Qt);
 	  for(int ip=k2; ip<Nm; ++ip) 
-	    qr_decomp(lmd,lme,Nm,Nm,Qt,lmd2[ip],k1,Nm);
+	    qr_decomp(eval,lme,Nm,Nm,Qt,eval2[ip],k1,Nm);
     
-	  for(int i=0; i<(Nk_+1); ++i) B[i] = 0.0;
+	  for(int i=0; i<(Nk+1); ++i) B[i] = 0.0;
 	  
 	  for(int j=k1-1; j<k2+1; ++j){
 	    for(int k=0; k<Nm; ++k){
@@ -302,62 +374,58 @@ public:
 	  // Compressed vector f and beta(k2)
 	  f *= Qt[Nm-1+Nm*(k2-1)];
 	  f += lme[k2-1] * evec[k2];
-	  beta_k = f * f;
+	  beta_k = norm2(f);
 	  beta_k = sqrt(beta_k);
-	  std::cout<<" beta(k) = "<<beta_k<<endl;
+	  std::cout<<" beta(k) = "<<beta_k<<std::endl;
 
 	  RealD betar = 1.0/beta_k;
 	  evec[k2] = betar * f;
 	  lme[k2-1] = beta_k;
 
 	  // Convergence test
-	  for(int k=0; k<Nm2; ++k){    
-	    lmd2[k] = lmd[k];
+	  for(int k=0; k<Nm; ++k){    
+	    eval2[k] = eval[k];
 	    lme2[k] = lme[k];
 	  }
 	  setUnit_Qt(Nm,Qt);
-	  diagonalize(lmd2,lme2,Nk_,Nm,Qt);
+	  diagonalize(eval2,lme2,Nk,Nm,Qt);
 	  
-	  for(int k = 0; k<Nk_; ++k) B[k]=0.0;
+	  for(int k = 0; k<Nk; ++k) B[k]=0.0;
 	  
-	  for(int j = 0; j<Nk_; ++j){
-	    for(int k = 0; k<Nk_; ++k){
+	  for(int j = 0; j<Nk; ++j){
+	    for(int k = 0; k<Nk; ++k){
 	      B[j] += Qt[k+j*Nm] * evec[k];
 	    }
 	  }
-	  Kdis = 0;
-	  Kthrs = 0;
 
-	  std::cout << setiosflags(ios_base::scientific);
-	  for(int i=0; i<Nk_; ++i){
-	    v = opr_->mult(B[i]);
-	    //std::cout<<"vv="<<v*v<<std::endl;
+	  Nconv = 0;
+	  //	  std::cout << std::setiosflags(std::ios_base::scientific);
+	  for(int i=0; i<Nk; ++i){
+
+	    _poly(_Linop,B[i],v);
 	    
-	    RealD vnum = B[i]*v;
-	    RealD vden = B[i]*B[i];
-	    lmd2[i] = vnum/vden;
-	    v -= lmd2[i]*B[i];
-	    RealD vv = v*v;
+	    RealD vnum = real(innerProduct(B[i],v)); // HermOp.
+	    RealD vden = norm2(B[i]);
+	    eval2[i] = vnum/vden;
+	    v -= eval2[i]*B[i];
+	    RealD vv = norm2(v);
 	    
-	    std::cout << " [" << setw(3)<< setiosflags(ios_base::right) <<i<<"] ";
-	    std::cout << setw(25)<< setiosflags(ios_base::left)<< lmd2[i];
-	    std::cout <<"  "<< setw(25)<< setiosflags(ios_base::right)<< vv<< endl;
+	    std::cout << "[" << 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::endl;
 	    
-	    if(vv<enorm_){
-	      Iconv[Kdis] = i;
-	      ++Kdis;
-	      if(sort_->saturated(lmd2[i],vthr)) ++Kthrs;
-	      std::cout<<"Kthrs="<<Kthrs<<endl;
+	    if(vv<eresid*eresid){
+	      Iconv[Nconv] = i;
+	      ++Nconv;
 	    }
+
 	  }  // i-loop end
-	  std::cout << resetiosflags(ios_base::scientific);
+	  //	  std::cout << std::resetiosflags(std::ios_base::scientific);
 
 
-	  std::cout<<" #modes converged: "<<Kdis<<endl;
+	  std::cout<<" #modes converged: "<<Nconv<<std::endl;
 
-	  if(Kthrs > 0){
-	    // (there is a converged eigenvalue larger than Vthrs.)
-	    Nconv = iter;
+	  if( Nconv>=Nk ){
 	    goto converged;
 	  }
 	} // end of iter loop
@@ -368,21 +436,586 @@ public:
       converged:
 	// Sorting
 	
-	lmd.clear();
+	eval.clear();
 	evec.clear();
-	for(int i=0; i<Kdis; ++i){
-	  lmd.push_back(lmd2[Iconv[i]]);
+	for(int i=0; i<Nconv; ++i){
+	  eval.push_back(eval2[Iconv[i]]);
 	  evec.push_back(B[Iconv[i]]);
 	}
-	sort_->push(lmd,evec,Kdis);
-	Nsbt = Kdis - Kthrs;
+	_sort.push(eval,evec,Nconv);
 	
 	std::cout << "\n Converged\n Summary :\n";
 	std::cout << " -- Iterations  = "<< Nconv  << "\n";
 	std::cout << " -- beta(k)     = "<< beta_k << "\n";
-	std::cout << " -- Kdis        = "<< Kdis   << "\n";
-	std::cout << " -- Nsbt        = "<< Nsbt   << "\n";
+	std::cout << " -- Nconv       = "<< Nconv  << "\n";
       }
-  };
+
+
+    // Adapted from Rudy's lanczos factor routine
+    int Lanczos_Factor(int start, int end,  int cont,
+		       DenseVector<Field> & bq, 
+		       Field &bf,
+		       DenseMatrix<RealD> &H){
+      
+      GridBase *grid = bq[0]._grid;
+
+      RealD beta;  
+      RealD sqbt;  
+      RealD alpha;
+
+      for(int i=start;i<Nm;i++){
+	for(int j=start;j<Nm;j++){
+	  H[i][j]=0.0;
+	}
+      }
+
+      std::cout<<"Lanczos_Factor start/end " <<start <<"/"<<end<<std::endl;
+
+      // Starting from scratch, bq[0] contains a random vector and |bq[0]| = 1
+      int first;
+      if(start == 0){
+
+	std::cout << "start == 0\n"; //TESTING
+
+	_poly(_Linop,bq[0],bf);
+
+	alpha = real(innerProduct(bq[0],bf));//alpha =  bq[0]^dag A bq[0]
+
+	std::cout << "alpha = " << alpha << std::endl;
+	
+	bf = bf - alpha * bq[0];  //bf =  A bq[0] - alpha bq[0]
+
+	H[0][0]=alpha;
+
+	std::cout << "Set H(0,0) to " << H[0][0] << std::endl;
+
+	first = 1;
+
+      } else {
+
+	first = start;
+
+      }
+
+      // I think start==0 and cont==zero are the same. Test this
+      // If so I can drop "cont" parameter?
+      if( cont ) assert(start!=0);
+
+      if( start==0 ) assert(cont!=0);
+
+      if( cont){
+
+	beta = 0;sqbt = 0;
+
+	std::cout << "cont is true so setting beta to zero\n";
+
+      }	else {
+
+	beta = norm2(bf);
+	sqbt = sqrt(beta);
+
+	std::cout << "beta = " << beta << std::endl;
+      }
+
+      for(int j=first;j<end;j++){
+
+	std::cout << "Factor j " << j <<std::endl;
+
+	if(cont){ // switches to factoring; understand start!=0 and initial bf value is right.
+	  bq[j] = bf; cont = false;
+	}else{
+	  bq[j] = (1.0/sqbt)*bf ;
+
+	  H[j][j-1]=H[j-1][j] = sqbt;
+	}
+
+	_poly(_Linop,bq[j],bf);
+
+	bf = bf - (1.0/sqbt)*bq[j-1]; 	       //bf = A bq[j] - beta bq[j-1] // PAB this comment was incorrect in beta term??
+
+	alpha = real(innerProduct(bq[j],bf));  //alpha = bq[j]^dag A bq[j]
+
+	bf = bf - alpha*bq[j];                 //bf = A bq[j] - beta bq[j-1] - alpha bq[j]
+	RealD fnorm = norm2(bf);
+
+	RealD bck = sqrt( real( conjugate(alpha)*alpha ) + beta );
+
+	beta = fnorm;
+	sqbt = sqrt(beta);
+	std::cout << "alpha = " << alpha << " fnorm = " << fnorm << '\n';
+
+	///Iterative refinement of orthogonality V = [ bq[0]  bq[1]  ...  bq[M] ]
+#if 0
+	int re = 0;
+	while( re == ref || (sqbt < rho * bck && re < 5) ){
+
+	  //bex = V^dag bf
+	  DenseVector<ComplexD> bex(j+1);
+	  for(int k=0;k<j+1;k++){
+	    bex[k] = innerProduct(bq[k],bf);
+	  }
+	  
+	  zero_fermion(tmp2);
+	  //tmp2 = V s
+	  for(int l=0;l<j+1;l++){
+	    RealD nrm = norm2(bq[l]);
+	    axpy(tmp1,0.0,bq[l],bq[l]); scale(tmp1,bex[l]); 	//tmp1 = V[j] bex[j]
+	    axpy(tmp2,1.0,tmp2,tmp1);					//tmp2 += V[j] bex[j]
+	  }
+
+	  //bf = bf - V V^dag bf.   Subtracting off any component in span { V[j] } 
+	  RealD btc = axpy_norm(bf,-1.0,tmp2,bf);
+	  alpha = alpha + bex[j];	      sqbt = sqrt(real(btc));	      
+	  RealD nmbex = 0;for(int k=0;k<j+1;k++){nmbex = nmbex + real( conjugate(bex[k])*bex[k]  );}
+	  bck = sqrt( nmbex );
+	  re++;
+	}
+	std::cout << "Iteratively refined orthogonality, changes alpha\n";
+	if(re > 1) std::cout << "orthagonality refined " << re << " times" <<std::endl;
+#endif
+	H[j][j]=alpha;
+      }
+
+      return end;
+    }
+
+    void EigenSort(DenseVector<double> evals,
+		   DenseVector<Field>  evecs){
+      int N= evals.size();
+      _sort.push(evals,evecs, evals.size(),N);
+    }
+
+    void ImplicitRestart(int TM, DenseVector<RealD> &evals,  DenseVector<DenseVector<RealD> > &evecs, DenseVector<Field> &bq, Field &bf, int cont)
+    {
+      std::cout << "ImplicitRestart begin. Eigensort starting\n";
+
+      DenseMatrix<RealD> H; Resize(H,Nm,Nm);
+
+      EigenSort(evals, evecs);
+
+      ///Assign shifts
+      int K=Nk;
+      int M=Nm;
+      int P=Np;
+      int converged=0;
+      if(K - converged < 4) P = (M - K-1); //one
+      //      DenseVector<RealD> shifts(P + shift_extra.size());
+      DenseVector<RealD> shifts(P);
+      for(int k = 0; k < P; ++k)
+	shifts[k] = evals[k]; 
+
+      /// Shift to form a new H and q
+      DenseMatrix<RealD> Q; Resize(Q,TM,TM);
+      Unity(Q);
+      Shift(Q, shifts); // H is implicitly passed in in Rudy's Shift routine
+
+      int ff = K;
+
+      /// Shifted H defines a new K step Arnoldi factorization
+      RealD  beta = H[ff][ff-1]; 
+      RealD  sig  = Q[TM - 1][ff - 1];
+      std::cout << "beta = " << beta << " sig = " << real(sig) <<std::endl;
+
+      std::cout << "TM = " << TM << " ";
+      std::cout << norm2(bq[0]) << " -- before" <<std::endl;
+
+      /// q -> q Q
+      times_real(bq, Q, TM);
+
+      std::cout << norm2(bq[0]) << " -- after " << ff <<std::endl;
+      bf =  beta* bq[ff] + sig* bf;
+
+      /// Do the rest of the factorization
+      ff = Lanczos_Factor(ff, M,cont,bq,bf,H);
+      
+      if(ff < M)
+	Abort(ff, evals, evecs);
+    }
+
+///Run the Eigensolver
+    void Run(int cont, DenseVector<Field> &bq, Field &bf, DenseVector<DenseVector<RealD> > & evecs,DenseVector<RealD> &evals)
+    {
+      init();
+
+      int M=Nm;
+
+      DenseMatrix<RealD> H; Resize(H,Nm,Nm);
+      Resize(evals,Nm,Nm);
+      Resize(evecs,Nm);
+
+      int ff = Lanczos_Factor(0, M, cont, bq,bf,H); // 0--M to begin with
+
+      if(ff < M) {
+	std::cout << "Krylov: aborting ff "<<ff <<" "<<M<<std::endl;
+	abort(); // Why would this happen?
+      }
+
+      int itcount = 0;
+      bool stop = false;
+
+      for(int it = 0; it < Niter && (converged < Nk); ++it) {
+
+	std::cout << "Krylov: Iteration --> " << it << std::endl;
+	int lock_num = lock ? converged : 0;
+	DenseVector<RealD> tevals(M - lock_num );
+	DenseMatrix<RealD> tevecs; Resize(tevecs,M - lock_num,M - lock_num);
+	  
+	//check residual of polynominal 
+	TestConv(H,M, tevals, tevecs);
+
+	if(converged >= Nk)
+	    break;
+
+	ImplicitRestart(ff, tevals,tevecs,H);
+      }
+      Wilkinson<RealD>(H, evals, evecs, small); 
+      //      Check();
+
+      std::cout << "Done  "<<std::endl;
+
+    }
+
+   ///H - shift I = QR; H = Q* H Q
+    void Shift(DenseMatrix<RealD> & H,DenseMatrix<RealD> &Q, DenseVector<RealD> shifts) {
+      
+      int P; Size(shifts,P);
+      int M; SizeSquare(Q,M);
+
+      Unity(Q);
+
+      int lock_num = lock ? converged : 0;
+
+      RealD t_Househoulder_vector(0.0);
+      RealD t_Househoulder_mult(0.0);
+
+      for(int i=0;i<P;i++){
+
+	RealD x, y, z;
+	DenseVector<RealD> ck(3), v(3);
+	  
+	x = H[lock_num+0][lock_num+0]-shifts[i];
+	y = H[lock_num+1][lock_num+0];
+	ck[0] = x; ck[1] = y; ck[2] = 0; 
+
+	normalise(ck);	///Normalization cancels in PHP anyway
+	RealD beta;
+
+	Householder_vector<RealD>(ck, 0, 2, v, beta);
+
+	Householder_mult<RealD>(H,v,beta,0,lock_num+0,lock_num+2,0);
+	Householder_mult<RealD>(H,v,beta,0,lock_num+0,lock_num+2,1);
+	///Accumulate eigenvector
+	Householder_mult<RealD>(Q,v,beta,0,lock_num+0,lock_num+2,1);
+	  
+	int sw = 0;
+	for(int k=lock_num+0;k<M-2;k++){
+
+	  x = H[k+1][k]; 
+	  y = H[k+2][k]; 
+	  z = (RealD)0.0;
+	  if(k+3 <= M-1){
+	    z = H[k+3][k];
+	  }else{
+	    sw = 1; v[2] = 0.0;
+	  }
+
+	  ck[0] = x; ck[1] = y; ck[2] = z;
+
+	  normalise(ck);
+
+	  Householder_vector<RealD>(ck, 0, 2-sw, v, beta);
+	  Householder_mult<RealD>(H,v, beta,0,k+1,k+3-sw,0);
+	  Householder_mult<RealD>(H,v, beta,0,k+1,k+3-sw,1);
+	  ///Accumulate eigenvector
+	  Householder_mult<RealD>(Q,v, beta,0,k+1,k+3-sw,1);
+	}
+      }
+    }
+
+    void TestConv(DenseMatrix<RealD> & H,int SS, 
+		  DenseVector<Field> &bq, Field &bf,
+		  DenseVector<RealD> &tevals, DenseVector<DenseVector<RealD> > &tevecs, 
+		  int lock, int converged)
+    {
+      std::cout << "Converged " << converged << " so far." << std::endl;
+      int lock_num = lock ? converged : 0;
+      int M = Nm;
+
+      ///Active Factorization
+      DenseMatrix<RealD> AH; Resize(AH,SS - lock_num,SS - lock_num );
+
+      AH = GetSubMtx(H,lock_num, SS, lock_num, SS);
+
+      int NN=tevals.size();
+      int AHsize=SS-lock_num;
+
+      RealD small=1.0e-16;
+      Wilkinson<RealD>(AH, tevals, tevecs, small);
+
+      EigenSort(tevals, tevecs);
+
+      RealD resid_nrm=  norm2(bf);
+
+      if(!lock) converged = 0;
+
+      for(int i = SS - lock_num - 1; i >= SS - Nk && i >= 0; --i){
+
+	RealD diff = 0;
+	diff = abs(tevecs[i][Nm - 1 - lock_num]) * resid_nrm;
+
+	std::cout << "residual estimate " << SS-1-i << " " << diff << " of (" << tevals[i] << ")" << std::endl;
+
+	if(diff < converged) {
+
+	  if(lock) {
+	    
+	    DenseMatrix<RealD> Q; Resize(Q,M,M);
+	    bool herm = true; 
+
+	    Lock(H, Q, tevals[i], converged, small, SS, herm);
+
+	    times_real(bq, Q, bq.size());
+	    bf = Q[M - 1][M - 1]* bf;
+	    lock_num++;
+	  }
+	  converged++;
+	  std::cout << " converged on eval " << converged << " of " << Nk << std::endl;
+	} else {
+	  break;
+	}
+      }
+      std::cout << "Got " << converged << " so far " <<std::endl;	
+    }
+#if 0
+    ///Check
+    void Check(void) {
+
+      DenseVector<RealD> goodval(get);
+      EigenSort(evals,evecs);
+
+      int NM = Nm;
+      int Nget = this->get;
+      S **V;
+      V = new S* [NM];
+
+      RealD *QZ;
+      QZ = new RealD [NM*NM];
+      for(int i = 0; i < NM; i++){
+	for(int j = 0; j < NM; j++){
+
+	  QZ[i*NM+j] = this->evecs[i][j];
+      
+          int f_size_cb = 24*dop.cbLs*dop.node_cbvol;
+
+	  for(int cb = this->prec; cb < 2; cb++){
+	    for(int i = 0; i < NM; i++){
+	      V[i] = (S*)(this->bq[i][cb]);
+
+	      const int m0 = 4 * 4; // this is new code
+	      assert(m0 % 16 == 0); // see the reason in VtimesQ.C
+
+	      const int row_per_thread = f_size_cb / (bfmarg::threads);
+	      {
+
+		{
+		  DenseVector<RealD> vrow_tmp0(m0*NM);
+		  DenseVector<RealD> vrow_tmp1(m0*NM);
+		  RealD *row_tmp0 = vrow_tmp0.data();
+		  RealD *row_tmp1 = vrow_tmp1.data();
+		  VtimesQ(QZ, NM, V, row_tmp0, row_tmp1, id * row_per_thread, m0, (id + 1) * row_per_thread);
+		}
+	      }
+	    }
+	  }
+	}
+      }
+    }
+#endif
+
+
+/**
+   There is some matrix Q such that for any vector y
+   Q.e_1 = y and Q is unitary.
+**/
+  template<class T>
+  static T orthQ(DenseMatrix<T> &Q, DenseVector<T> y){
+    int N = y.size();	//Matrix Size
+    Fill(Q,0.0);
+    T tau;
+    for(int i=0;i<N;i++){
+      Q[i][0]=y[i];
+    }
+    T sig = conj(y[0])*y[0];
+    T tau0 = abs(sqrt(sig));
+    
+    for(int j=1;j<N;j++){
+      sig += conj(y[j])*y[j]; 
+      tau = abs(sqrt(sig) ); 	
+
+      if(abs(tau0) > 0.0){
+	
+	T gam = conj( (y[j]/tau)/tau0 );
+	for(int k=0;k<=j-1;k++){  
+	  Q[k][j]=-gam*y[k];
+	}
+	Q[j][j]=tau0/tau;
+      } else {
+	Q[j-1][j]=1.0;
+      }
+      tau0 = tau;
+    }
+    return tau;
+  }
+
+/**
+	There is some matrix Q such that for any vector y
+	Q.e_k = y and Q is unitary.
+**/
+  template< class T>
+  static T orthU(DenseMatrix<T> &Q, DenseVector<T> y){
+    T tau = orthQ(Q,y);
+    SL(Q);
+    return tau;
+  }
+
+
+/**
+	Wind up with a matrix with the first con rows untouched
+
+say con = 2
+	Q is such that Qdag H Q has {x, x, val, 0, 0, 0, 0, ...} as 1st colum
+	and the matrix is upper hessenberg
+	and with f and Q appropriately modidied with Q is the arnoldi factorization
+
+**/
+
+template<class T>
+static void Lock(DenseMatrix<T> &H, 	///Hess mtx	
+		 DenseMatrix<T> &Q, 	///Lock Transform
+		 T val, 		///value to be locked
+		 int con, 	///number already locked
+		 RealD small,
+		 int dfg,
+		 bool herm)
+{	
+
+
+  //ForceTridiagonal(H);
+
+  int M = H.dim;
+  DenseVector<T> vec; Resize(vec,M-con);
+
+  DenseMatrix<T> AH; Resize(AH,M-con,M-con);
+  AH = GetSubMtx(H,con, M, con, M);
+
+  DenseMatrix<T> QQ; Resize(QQ,M-con,M-con);
+
+  Unity(Q);   Unity(QQ);
+  
+  DenseVector<T> evals; Resize(evals,M-con);
+  DenseMatrix<T> evecs; Resize(evecs,M-con,M-con);
+
+  Wilkinson<T>(AH, evals, evecs, small);
+
+  int k=0;
+  RealD cold = abs( val - evals[k]); 
+  for(int i=1;i<M-con;i++){
+    RealD cnew = abs( val - evals[i]);
+    if( cnew < cold ){k = i; cold = cnew;}
+  }
+  vec = evecs[k];
+
+  ComplexD tau;
+  orthQ(QQ,vec);
+  //orthQM(QQ,AH,vec);
+
+  AH = Hermitian(QQ)*AH;
+  AH = AH*QQ;
+	
+
+  for(int i=con;i<M;i++){
+    for(int j=con;j<M;j++){
+      Q[i][j]=QQ[i-con][j-con];
+      H[i][j]=AH[i-con][j-con];
+    }
+  }
+
+  for(int j = M-1; j>con+2; j--){
+
+    DenseMatrix<T> U; Resize(U,j-1-con,j-1-con);
+    DenseVector<T> z; Resize(z,j-1-con); 
+    T nm = norm(z); 
+    for(int k = con+0;k<j-1;k++){
+      z[k-con] = conj( H(j,k+1) );
+    }
+    normalise(z);
+
+    RealD tmp = 0;
+    for(int i=0;i<z.size()-1;i++){tmp = tmp + abs(z[i]);}
+
+    if(tmp < small/( (RealD)z.size()-1.0) ){ continue;}	
+
+    tau = orthU(U,z);
+
+    DenseMatrix<T> Hb; Resize(Hb,j-1-con,M);	
+	
+    for(int a = 0;a<M;a++){
+      for(int b = 0;b<j-1-con;b++){
+	T sum = 0;
+	for(int c = 0;c<j-1-con;c++){
+	  sum += H[a][con+1+c]*U[c][b];
+	}//sum += H(a,con+1+c)*U(c,b);}
+	Hb[b][a] = sum;
+      }
+    }
+	
+    for(int k=con+1;k<j;k++){
+      for(int l=0;l<M;l++){
+	H[l][k] = Hb[k-1-con][l];
+      }
+    }//H(Hb[k-1-con][l] , l,k);}}
+
+    DenseMatrix<T> Qb; Resize(Qb,M,M);	
+	
+    for(int a = 0;a<M;a++){
+      for(int b = 0;b<j-1-con;b++){
+	T sum = 0;
+	for(int c = 0;c<j-1-con;c++){
+	  sum += Q[a][con+1+c]*U[c][b];
+	}//sum += Q(a,con+1+c)*U(c,b);}
+	Qb[b][a] = sum;
+      }
+    }
+	
+    for(int k=con+1;k<j;k++){
+      for(int l=0;l<M;l++){
+	Q[l][k] = Qb[k-1-con][l];
+      }
+    }//Q(Qb[k-1-con][l] , l,k);}}
+
+    DenseMatrix<T> Hc; Resize(Hc,M,M);	
+	
+    for(int a = 0;a<j-1-con;a++){
+      for(int b = 0;b<M;b++){
+	T sum = 0;
+	for(int c = 0;c<j-1-con;c++){
+	  sum += conj( U[c][a] )*H[con+1+c][b];
+	}//sum += conj( U(c,a) )*H(con+1+c,b);}
+	Hc[b][a] = sum;
+      }
+    }
+
+    for(int k=0;k<M;k++){
+      for(int l=con+1;l<j;l++){
+	H[l][k] = Hc[k][l-1-con];
+      }
+    }//H(Hc[k][l-1-con] , l,k);}}
+
+  }
+}
+
+
+ };
+
 }
 #endif
diff --git a/lib/algorithms/iterative/Matrix.h b/lib/algorithms/iterative/Matrix.h
new file mode 100644
index 00000000..a455e82e
--- /dev/null
+++ b/lib/algorithms/iterative/Matrix.h
@@ -0,0 +1,426 @@
+#ifndef MATRIX_H
+#define MATRIX_H
+
+#include <cstdlib>
+#include <string>
+#include <cmath>
+#include <vector>
+#include <iostream>
+#include <iomanip>
+#include <complex>
+#include <typeinfo>
+#include <Grid.h>
+
+
+/** Sign function **/
+template <class T> T sign(T p){return ( p/abs(p) );}
+
+/////////////////////////////////////////////////////////////////////////////////////////////////////////
+///////////////////// Hijack STL containers for our wicked means /////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////////////////////////////////////
+template<class T> using Vector = Vector<T>;
+template<class T> using Matrix = Vector<Vector<T> >;
+
+template<class T> void Resize(Vector<T > & vec, int N) { vec.resize(N); }
+
+template<class T> void Resize(Matrix<T > & mat, int N, int M) { 
+  mat.resize(N);
+  for(int i=0;i<N;i++){
+    mat[i].resize(M);
+  }
+}
+template<class T> void Size(Vector<T> & vec, int &N) 
+{ 
+  N= vec.size();
+}
+template<class T> void Size(Matrix<T> & mat, int &N,int &M) 
+{ 
+  N= mat.size();
+  M= mat[0].size();
+}
+template<class T> void SizeSquare(Matrix<T> & mat, int &N) 
+{ 
+  int M; Size(mat,N,M);
+  assert(N==M);
+}
+template<class T> void SizeSame(Matrix<T> & mat1,Matrix<T> &mat2, int &N1,int &M1) 
+{ 
+  int N2,M2;
+  Size(mat1,N1,M1);
+  Size(mat2,N2,M2);
+  assert(N1==N2);
+  assert(M1==M2);
+}
+
+//*****************************************
+//*	(Complex) Vector operations	*
+//*****************************************
+
+/**Conj of a Vector **/
+template <class T> Vector<T> conj(Vector<T> p){
+	Vector<T> q(p.size());
+	for(int i=0;i<p.size();i++){q[i] = conj(p[i]);}
+	return q;
+}
+
+/** Norm of a Vector**/
+template <class T> T norm(Vector<T> p){
+	T sum = 0;
+	for(int i=0;i<p.size();i++){sum = sum + p[i]*conj(p[i]);}
+	return abs(sqrt(sum));
+}
+
+/** Norm squared of a Vector **/
+template <class T> T norm2(Vector<T> p){
+	T sum = 0;
+	for(int i=0;i<p.size();i++){sum = sum + p[i]*conj(p[i]);}
+	return abs((sum));
+}
+
+/** Sum elements of a Vector **/
+template <class T> T trace(Vector<T> p){
+	T sum = 0;
+	for(int i=0;i<p.size();i++){sum = sum + p[i];}
+	return sum;
+}
+
+/** Fill a Vector with constant c **/
+template <class T> void Fill(Vector<T> &p, T c){
+	for(int i=0;i<p.size();i++){p[i] = c;}
+}
+/** Normalize a Vector **/
+template <class T> void normalize(Vector<T> &p){
+	T m = norm(p);
+	if( abs(m) > 0.0) for(int i=0;i<p.size();i++){p[i] /= m;}
+}
+/** Vector by scalar **/
+template <class T, class U> Vector<T> times(Vector<T> p, U s){
+	for(int i=0;i<p.size();i++){p[i] *= s;}
+	return p;
+}
+template <class T, class U> Vector<T> times(U s, Vector<T> p){
+	for(int i=0;i<p.size();i++){p[i] *= s;}
+	return p;
+}
+/** inner product of a and b = conj(a) . b **/
+template <class T> T inner(Vector<T> a, Vector<T> b){
+	T m = 0.;
+	for(int i=0;i<a.size();i++){m = m + conj(a[i])*b[i];}
+	return m;
+}
+/** sum of a and b = a + b **/
+template <class T> Vector<T> add(Vector<T> a, Vector<T> b){
+	Vector<T> m(a.size());
+	for(int i=0;i<a.size();i++){m[i] = a[i] + b[i];}
+	return m;
+}
+/** sum of a and b = a - b **/
+template <class T> Vector<T> sub(Vector<T> a, Vector<T> b){
+	Vector<T> m(a.size());
+	for(int i=0;i<a.size();i++){m[i] = a[i] - b[i];}
+	return m;
+}
+
+/** 
+ *********************************
+ *	Matrices	         *
+ *********************************
+ **/
+
+template<class T> void Fill(Matrix<T> & mat, T&val) { 
+  int N,M;
+  Size(mat,N,M);
+  for(int i=0;i<N;i++){
+  for(int j=0;j<M;j++){
+    mat[i][j] = val;
+  }}
+}
+
+/** Transpose of a matrix **/
+Matrix<T> Transpose(Matrix<T> & mat){
+  int N,M;
+  Size(mat,N,M);
+  Matrix C; Resize(C,M,N);
+  for(int i=0;i<M;i++){
+  for(int j=0;j<N;j++){
+    C[i][j] = mat[j][i];
+  }} 
+  return C;
+}
+/** Set Matrix to unit matrix **/
+template<class T> void Unity(Matrix<T> &mat){
+  int N;  SizeSquare(mat,N);
+  for(int i=0;i<N;i++){
+    for(int j=0;j<N;j++){
+      if ( i==j ) A[i][j] = 1;
+      else        A[i][j] = 0;
+    } 
+  } 
+}
+/** Add C * I to matrix **/
+template<class T>
+void PlusUnit(Matrix<T> & A,T c){
+  int dim;  SizeSquare(A,dim);
+  for(int i=0;i<dim;i++){A[i][i] = A[i][i] + c;} 
+}
+
+/** return the Hermitian conjugate of matrix **/
+Matrix<T> HermitianConj(Matrix<T> &mat){
+
+  int dim; SizeSquare(mat,dim);
+
+  Matrix<T> C; Resize(C,dim,dim);
+
+  for(int i=0;i<dim;i++){
+    for(int j=0;j<dim;j++){
+      C[i][j] = conj(mat[j][i]);
+    } 
+  } 
+  return C;
+}
+
+/** return diagonal entries as a Vector **/
+Vector<T> diag(Matrix<T> &A)
+{
+  int dim; SizeSquare(A,dim);
+  Vector<T> d; Resize(d,dim);
+
+  for(int i=0;i<dim;i++){
+    d[i] = A[i][i];
+  }
+  return d;
+}
+
+/** Left multiply by a Vector **/
+Vector<T> operator *(Vector<T> &B,Matrix<T> &A)
+{
+  int K,M,N; 
+  Size(B,K);
+  Size(A,M,N);
+  assert(K==M);
+  
+  Vector<T> C; Resize(C,N);
+
+  for(int j=0;j<N;j++){
+    T sum = 0.0;
+    for(int i=0;i<M;i++){
+      sum += B[i] * A[i][j];
+    }
+    C[j] =  sum;
+  }
+  return C; 
+}
+
+/** return 1/diagonal entries as a Vector **/
+Vector<T> inv_diag(Matrix<T> & A){
+  int dim; SizeSquare(A,dim);
+  Vector<T> d; Resize(d,dim);
+  for(int i=0;i<dim;i++){
+    d[i] = 1.0/A[i][i];
+  }
+  return d;
+}
+/** Matrix Addition **/
+inline Matrix<T> operator + (Matrix<T> &A,Matrix<T> &B)
+{
+  int N,M  ; SizeSame(A,B,N,M);
+  Matrix C; Resize(C,N,M);
+  for(int i=0;i<N;i++){
+    for(int j=0;j<M;j++){
+      C[i][j] = A[i][j] +  B[i][j];
+    } 
+  } 
+  return C;
+} 
+/** Matrix Subtraction **/
+inline Matrix<T> operator- (Matrix<T> & A,Matrix<T> &B){
+  int N,M  ; SizeSame(A,B,N,M);
+  Matrix C; Resize(C,N,M);
+  for(int i=0;i<N;i++){
+  for(int j=0;j<M;j++){
+    C[i][j] = A[i][j] -  B[i][j];
+  }}
+  return C;
+} 
+
+/** Matrix scalar multiplication **/
+inline Matrix<T> operator* (Matrix<T> & A,T c){
+  int N,M; Size(A,N,M);
+  Matrix C; Resize(C,N,M);
+  for(int i=0;i<N;i++){
+  for(int j=0;j<M;j++){
+    C[i][j] = A[i][j]*c;
+  }} 
+  return C;
+} 
+/** Matrix Matrix multiplication **/
+inline Matrix<T> operator* (Matrix<T> &A,Matrix<T> &B){
+  int K,L,N,M;
+  Size(A,K,L);
+  Size(B,N,M); assert(L==N);
+  Matrix C; Resize(C,K,M);
+
+  for(int i=0;i<K;i++){
+    for(int j=0;j<M;j++){
+      T sum = 0.0;
+      for(int k=0;k<N;k++) sum += A[i][k]*B[k][j];
+      C[i][j] =sum;
+    }
+  }
+  return C; 
+} 
+/** Matrix Vector multiplication **/
+inline Vector<T> operator* (Matrix<T> &A,Vector<T> &B){
+  int M,N,K;
+  Size(A,N,M);
+  Size(B,K); assert(K==M);
+  Vector<T> C; Resize(C,N);
+  for(int i=0;i<N;i++){
+    T sum = 0.0;
+    for(int j=0;j<M;j++) sum += A[i][j]*B[j];
+    C[i] =  sum;
+  }
+  return C; 
+} 
+
+/** Some version of Matrix norm **/
+/*
+inline T Norm(){ // this is not a usual L2 norm
+    T norm = 0;
+    for(int i=0;i<dim;i++){
+      for(int j=0;j<dim;j++){
+	norm += abs(A[i][j]);
+    }}
+    return norm;
+  }
+*/
+
+/** Some version of Matrix norm **/
+template<class T> T LargestDiag(Matrix<T> &A)
+{
+  int dim ; SizeSquare(A,dim); 
+
+  T ld = abs(A[0][0]);
+  for(int i=1;i<dim;i++){
+    T cf = abs(A[i][i]);
+    if(abs(cf) > abs(ld) ){ld = cf;}
+  }
+  return ld;
+}
+
+/** Look for entries on the leading subdiagonal that are smaller than 'small' **/
+template <class T,class U> int Chop_subdiag(Matrix<T> &A,T norm, int offset, U small)
+{
+  int dim; SizeSquare(A,dim);
+  for(int l = dim - 1 - offset; l >= 1; l--) {             		
+    if((U)abs(A[l][l - 1]) < (U)small) {
+      A[l][l-1]=(U)0.0;
+      return l;
+    }
+  }
+  return 0;
+}
+
+/** Look for entries on the leading subdiagonal that are smaller than 'small' **/
+template <class T,class U> int Chop_symm_subdiag(Matrix<T> & A,T norm, int offset, U small) 
+{
+  int dim; SizeSquare(A,dim);
+  for(int l = dim - 1 - offset; l >= 1; l--) {
+    if((U)abs(A[l][l - 1]) < (U)small) {
+      A[l][l - 1] = (U)0.0;
+      A[l - 1][l] = (U)0.0;
+      return l;
+    }
+  }
+  return 0;
+}
+/**Assign a submatrix to a larger one**/
+template<class T>
+void AssignSubMtx(Matrix<T> & A,int row_st, int row_end, int col_st, int col_end, Matrix<T> &S)
+{
+  for(int i = row_st; i<row_end; i++){
+    for(int j = col_st; j<col_end; j++){
+      A[i][j] = S[i - row_st][j - col_st];
+    }
+  }
+}
+
+/**Get a square submatrix**/
+template <class T>
+Matrix<T> GetSubMtx(Matrix<T> &A,int row_st, int row_end, int col_st, int col_end)
+{
+  Matrix<T> H; Resize(row_end - row_st,col_end-col_st);
+
+  for(int i = row_st; i<row_end; i++){
+  for(int j = col_st; j<col_end; j++){
+    H[i-row_st][j-col_st]=A[i][j];
+  }}
+  return H;
+}
+  
+ /**Assign a submatrix to a larger one NB remember Vector Vectors are transposes of the matricies they represent**/
+template<class T>
+void AssignSubMtx(Matrix<T> & A,int row_st, int row_end, int col_st, int col_end, Matrix<T> &S)
+{
+  for(int i = row_st; i<row_end; i++){
+  for(int j = col_st; j<col_end; j++){
+    A[i][j] = S[i - row_st][j - col_st];
+  }}
+}
+  
+/** compute b_i A_ij b_j **/ // surprised no Conj
+template<class T> T proj(Matrix<T> A, Vector<T> B){
+  int dim; SizeSquare(A,dim);
+  int dimB; Size(B,dimB);
+  assert(dimB==dim);
+  T C = 0;
+  for(int i=0;i<dim;i++){
+    T sum = 0.0;
+    for(int j=0;j<dim;j++){
+      sum += A[i][j]*B[j];
+    }
+    C +=  B[i]*sum; // No conj?
+  }
+  return C; 
+}
+
+
+/*
+ *************************************************************
+ *
+ * Matrix Vector products
+ *
+ *************************************************************
+ */
+// Instead make a linop and call my CG;
+
+/// q -> q Q
+template <class T,class Fermion> void times(Vector<Fermion> &q, Matrix<T> &Q)
+{
+  int M; SizeSquare(Q,M);
+  int N; Size(q,N); 
+  assert(M==N);
+
+  times(q,Q,N);
+}
+
+/// q -> q Q
+template <class T> void times(multi1d<LatticeFermion> &q, Matrix<T> &Q, int N)
+{
+  GridBase *grid = q[0]._grid;
+  int M; SizeSquare(Q,M);
+  int K; Size(q,K); 
+  assert(N<M);
+  assert(N<K);
+  Vector<Fermion> S(N,grid );
+  for(int j=0;j<N;j++){
+    S[j] = zero;
+    for(int k=0;k<N;k++){
+      S[j] = S[j] +  q[k]* Q[k][j]; 
+    }
+  }
+  for(int j=0;j<q.size();j++){
+    q[j] = S[j];
+  }
+}
+#endif
diff --git a/lib/algorithms/iterative/bisec.c b/lib/algorithms/iterative/bisec.c
new file mode 100644
index 00000000..a5117aee
--- /dev/null
+++ b/lib/algorithms/iterative/bisec.c
@@ -0,0 +1,122 @@
+#include <math.h>
+#include <stdlib.h>
+#include <vector>
+
+struct Bisection {
+
+static void get_eig2(int row_num,std::vector<RealD> &ALPHA,std::vector<RealD> &BETA, std::vector<RealD> & eig)
+{
+  int i,j;
+  std::vector<RealD> evec1(row_num+3);
+  std::vector<RealD> evec2(row_num+3);
+  RealD eps2;
+  ALPHA[1]=0.;
+  BETHA[1]=0.;
+  for(i=0;i<row_num-1;i++) {
+    ALPHA[i+1] = A[i*(row_num+1)].real();
+    BETHA[i+2] = A[i*(row_num+1)+1].real();
+  }
+  ALPHA[row_num] = A[(row_num-1)*(row_num+1)].real();
+  bisec(ALPHA,BETHA,row_num,1,row_num,1e-10,1e-10,evec1,eps2);
+  bisec(ALPHA,BETHA,row_num,1,row_num,1e-16,1e-16,evec2,eps2);
+
+  // Do we really need to sort here?
+  int begin=1;
+  int end = row_num;
+  int swapped=1;
+  while(swapped) {
+    swapped=0;
+    for(i=begin;i<end;i++){
+      if(mag(evec2[i])>mag(evec2[i+1]))	{
+	swap(evec2+i,evec2+i+1);
+	swapped=1;
+      }
+    }
+    end--;
+    for(i=end-1;i>=begin;i--){
+      if(mag(evec2[i])>mag(evec2[i+1]))	{
+	swap(evec2+i,evec2+i+1);
+	swapped=1;
+      }
+    }
+    begin++;
+  }
+
+  for(i=0;i<row_num;i++){
+    for(j=0;j<row_num;j++) {
+      if(i==j) H[i*row_num+j]=evec2[i+1];
+      else H[i*row_num+j]=0.;
+    }
+  }
+}
+
+static void bisec(std::vector<RealD> &c,   
+		  std::vector<RealD> &b,
+		  int n,
+		  int m1,
+		  int m2,
+		  RealD eps1,
+		  RealD relfeh,
+		  std::vector<RealD> &x,
+		  RealD &eps2)
+{
+  std::vector<RealD> wu(n+2);
+
+  RealD h,q,x1,xu,x0,xmin,xmax; 
+  int i,a,k;
+
+  b[1]=0.0;
+  xmin=c[n]-fabs(b[n]);
+  xmax=c[n]+fabs(b[n]);
+  for(i=1;i<n;i++){
+    h=fabs(b[i])+fabs(b[i+1]);
+    if(c[i]+h>xmax) xmax= c[i]+h;
+    if(c[i]-h<xmin) xmin= c[i]-h;
+  }
+  xmax *=2.;
+
+  eps2=relfeh*((xmin+xmax)>0.0 ? xmax : -xmin);
+  if(eps1<=0.0) eps1=eps2;
+  eps2=0.5*eps1+7.0*(eps2);
+  x0=xmax;
+  for(i=m1;i<=m2;i++){
+    x[i]=xmax;
+    wu[i]=xmin;
+  }
+
+  for(k=m2;k>=m1;k--){
+    xu=xmin;
+    i=k;
+    do{
+      if(xu<wu[i]){
+	xu=wu[i];
+	i=m1-1;
+      }
+      i--;
+    }while(i>=m1);
+    if(x0>x[k]) x0=x[k];
+    while((x0-xu)>2*relfeh*(fabs(xu)+fabs(x0))+eps1){
+      x1=(xu+x0)/2;
+
+      a=0;
+      q=1.0;
+      for(i=1;i<=n;i++){
+	q=c[i]-x1-((q!=0.0)? b[i]*b[i]/q:fabs(b[i])/relfeh);
+	if(q<0) a++;
+      }
+      //			printf("x1=%e a=%d\n",x1,a);
+      if(a<k){
+	if(a<m1){
+	  xu=x1;
+	  wu[m1]=x1;
+	}else {
+	  xu=x1;
+	  wu[a+1]=x1;
+	  if(x[a]>x1) x[a]=x1;
+	}
+      }else x0=x1;
+    }
+    x[k]=(x0+xu)/2;
+  }
+}
+}
diff --git a/lib/algorithms/iterative/get_eig.c b/lib/algorithms/iterative/get_eig.c
new file mode 100644
index 00000000..d3f5a12f
--- /dev/null
+++ b/lib/algorithms/iterative/get_eig.c
@@ -0,0 +1 @@
+
diff --git a/tests/Test_synthetic_lanczos.cc b/tests/Test_synthetic_lanczos.cc
new file mode 100644
index 00000000..c15ca3a7
--- /dev/null
+++ b/tests/Test_synthetic_lanczos.cc
@@ -0,0 +1,123 @@
+#include <fenv.h>
+#include <Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+static int
+FEenableexcept (unsigned int excepts)
+{
+  static fenv_t fenv;
+  unsigned int new_excepts = excepts & FE_ALL_EXCEPT,
+    old_excepts;  // previous masks
+
+  if ( fegetenv (&fenv) ) return -1;
+  old_excepts = fenv.__control & FE_ALL_EXCEPT;
+
+  // unmask
+  fenv.__control &= ~new_excepts;
+  fenv.__mxcsr   &= ~(new_excepts << 7);
+
+  return ( fesetenv (&fenv) ? -1 : old_excepts );
+}
+
+
+template<class Field> class DumbOperator  : public LinearOperatorBase<Field> {
+public:
+  LatticeComplex scale;
+
+  DumbOperator(GridBase *grid)    : scale(grid)
+  {
+    GridParallelRNG  pRNG(grid);  
+    std::vector<int> seeds({5,6,7,8});
+    pRNG.SeedFixedIntegers(seeds);
+
+    random(pRNG,scale);
+
+    scale = exp(-real(scale)*6.0);
+    std::cout << " True matrix \n"<< scale <<std::endl;
+  }
+
+  // Support for coarsening to a multigrid
+  void OpDiag (const Field &in, Field &out) {};
+  void OpDir  (const Field &in, Field &out,int dir,int disp){};
+
+  void Op     (const Field &in, Field &out){
+    out = scale * in;
+  }
+  void AdjOp  (const Field &in, Field &out){
+    out = scale * in;
+  }
+  void HermOp(const Field &in, Field &out){
+    double n1, n2;
+    HermOpAndNorm(in,out,n1,n2);
+  }
+  void HermOpAndNorm(const Field &in, Field &out,double &n1,double &n2){
+    ComplexD dot;
+
+    out = scale * in;
+
+    dot= innerProduct(in,out);
+    n1=real(dot);
+
+    dot = innerProduct(out,out);
+    n2=real(dot);
+  }
+};
+
+
+int main (int argc, char ** argv)
+{
+
+  FEenableexcept(FE_ALL_EXCEPT & ~FE_INEXACT); 
+
+  Grid_init(&argc,&argv);
+
+  GridCartesian *grid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), 
+						       GridDefaultSimd(Nd,vComplex::Nsimd()),
+						       GridDefaultMpi());
+
+  GridParallelRNG  RNG(grid);  
+  std::vector<int> seeds({1,2,3,4});
+  RNG.SeedFixedIntegers(seeds);
+
+
+  RealD alpha = 1.0;
+  RealD beta  = 0.03;
+  RealD mu    = 0.0;
+  int order = 11;
+  ChebyshevLanczos<LatticeComplex> Cheby(alpha,beta,mu,order);
+
+  std::ofstream file("pooh.dat");
+  Cheby.csv(file);
+
+  HermOpOperatorFunction<LatticeComplex> X;
+  DumbOperator<LatticeComplex> HermOp(grid);
+
+  const int Nk = 40;
+  const int Nm = 80;
+  const int Nit= 10000;
+
+  int Nconv;
+  RealD eresid = 1.0e-8;
+
+  ImplicitlyRestartedLanczos<LatticeComplex> IRL(HermOp,X,Nk,Nm,eresid,Nit);
+
+  ImplicitlyRestartedLanczos<LatticeComplex> ChebyIRL(HermOp,Cheby,Nk,Nm,eresid,Nit);
+
+  LatticeComplex src(grid); gaussian(RNG,src);
+  {
+    std::vector<RealD>          eval(Nm);
+    std::vector<LatticeComplex> evec(Nm,grid);
+    IRL.calc(eval,evec,src, Nconv);
+  }
+  
+  {
+    std::vector<RealD>          eval(Nm);
+    std::vector<LatticeComplex> evec(Nm,grid);
+    ChebyIRL.calc(eval,evec,src, Nconv);
+  }
+
+  Grid_finalize();
+}