Merge branch 'develop' of https://github.com/paboyle/Grid into feature/Lanczos

patch to compile with AVX512 for SkyLake Xeon processor using GCC7.2.…

.gitignore

@@ -88,12 +88,18 @@ Thumbs.db
 ###################
 build*/*
 
+# bootstrap #
+#############
+*.tar.bz2*
+
 # IDE related files #
 #####################
 *.xcodeproj/*
 build.sh
 .vscode
 *.code-workspace
+.ctags
+tags
 
 # Eigen source #
 ################

bootstrap.sh

@@ -1,9 +1,16 @@
 #!/usr/bin/env bash
 
-EIGEN_URL='http://bitbucket.org/eigen/eigen/get/3.3.3.tar.bz2'
+EIGEN_SRC='3.3.3.tar.bz2'
+EIGEN_URL="http://bitbucket.org/eigen/eigen/get/${EIGEN_SRC}"
 
-echo "-- deploying Eigen source..."
-wget ${EIGEN_URL} --no-check-certificate && ./scripts/update_eigen.sh `basename ${EIGEN_URL}` && rm `basename ${EIGEN_URL}`
+if [ -f ${EIGEN_SRC} ]; then
+  echo "-- skip deploying Eigen source..."
+else
+  echo "-- deploying Eigen source..."
+  wget ${EIGEN_URL} --no-check-certificate
+  ./scripts/update_eigen.sh `basename ${EIGEN_URL}`
+  #rm `basename ${EIGEN_URL}`
+fi
 
 echo '-- generating Make.inc files...'
 ./scripts/filelist

lib/algorithms/Algorithms.h

@@ -49,6 +49,8 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include <Grid/algorithms/iterative/BlockConjugateGradient.h>
 #include <Grid/algorithms/iterative/ConjugateGradientReliableUpdate.h>
 #include <Grid/algorithms/iterative/ImplicitlyRestartedLanczos.h>
+#include <Grid/algorithms/iterative/ImplicitlyRestartedLanczosCJ.h>
+#include <Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h>
 #include <Grid/algorithms/CoarsenedMatrix.h>
 #include <Grid/algorithms/FFT.h>
 

lib/algorithms/LinearOperator.h

@@ -373,6 +373,75 @@ namespace Grid {
     };
     template<class Matrix,class Field> using SchurStagOperator = SchurStaggeredOperator<Matrix,Field>;
 
+#if 0
+  // This is specific to (Z)mobius fermions
+  template<class Matrix, class Field>
+    class KappaSimilarityTransform {
+  public:
+//    INHERIT_IMPL_TYPES(Matrix);
+    typedef typename Matrix::Coeff_t                     Coeff_t;
+    std::vector<Coeff_t> kappa, kappaDag, kappaInv, kappaInvDag;
+
+    KappaSimilarityTransform (Matrix &zmob) {
+      for (int i=0;i<(int)zmob.bs.size();i++) {
+	Coeff_t k = 1.0 / ( 2.0 * (zmob.bs[i] *(4 - zmob.M5) + 1.0) );
+	kappa.push_back( k );
+	kappaDag.push_back( conj(k) );
+	kappaInv.push_back( 1.0 / k );
+	kappaInvDag.push_back( 1.0 / conj(k) );
+      }
+    }
+
+  template<typename vobj>
+    void sscale(const Lattice<vobj>& in, Lattice<vobj>& out, Coeff_t* s) {
+    GridBase *grid=out._grid;
+    out.checkerboard = in.checkerboard;
+    assert(grid->_simd_layout[0] == 1); // should be fine for ZMobius for now
+    int Ls = grid->_rdimensions[0];
+    parallel_for(int ss=0;ss<grid->oSites();ss++){
+      vobj tmp = s[ss % Ls]*in._odata[ss];
+      vstream(out._odata[ss],tmp);
+    }
+  }
+
+  RealD sscale_norm(const Field& in, Field& out, Coeff_t* s) {
+    sscale(in,out,s);
+    return norm2(out);
+  }
+
+  virtual RealD M       (const Field& in, Field& out) { return sscale_norm(in,out,&kappa[0]);   }
+  virtual RealD MDag    (const Field& in, Field& out) { return sscale_norm(in,out,&kappaDag[0]);}
+  virtual RealD MInv    (const Field& in, Field& out) { return sscale_norm(in,out,&kappaInv[0]);}
+  virtual RealD MInvDag (const Field& in, Field& out) { return sscale_norm(in,out,&kappaInvDag[0]);}
+
+  };
+
+  template<class Matrix,class Field>
+    class SchurDiagTwoKappaOperator :  public SchurOperatorBase<Field> {
+  public:
+    KappaSimilarityTransform<Matrix, Field> _S;
+    SchurDiagTwoOperator<Matrix, Field> _Mat;
+
+    SchurDiagTwoKappaOperator (Matrix &Mat): _S(Mat), _Mat(Mat) {};
+
+    virtual  RealD Mpc      (const Field &in, Field &out) {
+      Field tmp(in._grid);
+
+      _S.MInv(in,out);
+      _Mat.Mpc(out,tmp);
+      return _S.M(tmp,out);
+
+    }
+    virtual  RealD MpcDag   (const Field &in, Field &out){
+      Field tmp(in._grid);
+
+      _S.MDag(in,out);
+      _Mat.MpcDag(out,tmp);
+      return _S.MInvDag(tmp,out);
+    }
+  };
+#endif
+
 
     /////////////////////////////////////////////////////////////
     // Base classes for functions of operators

lib/algorithms/approx/Chebyshev.h

@@ -54,10 +54,16 @@ struct ChebyParams : Serializable {
 
   public:
     void csv(std::ostream &out){
-      RealD diff = hi-lo;
+
+#if 0
       RealD delta = (hi-lo)*1.0e-9;
       for (RealD x=lo; x<hi; x+=delta) {
 	delta*=1.1;
+#else
+	RealD diff = hi-lo;
+      //for (RealD x=lo-0.2*diff; x<hi+0.2*diff; x+=(hi-lo)/1000) {
+      for (RealD x=lo-0.2*diff; x<hi+0.2*diff; x+=diff/1000.0) { // ypj [note] divide by float
+#endif
 	RealD f = approx(x);
 	out<< x<<" "<<f<<std::endl;
       }
@@ -89,7 +95,7 @@ struct ChebyParams : Serializable {
       
       if(order < 2) exit(-1);
       Coeffs.resize(order);
-      Coeffs.assign(0.,order);
+      Coeffs.assign(order,0.);  
       Coeffs[order-1] = 1.;
     };
 

lib/algorithms/densematrix/DenseMatrix.h

@@ -0,0 +1,137 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/iterative/DenseMatrix.h
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: paboyle <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 */
+#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(DenseVector<T > & mat, int N) { 
+  mat.resize(N);
+}
+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 "Householder.h"
+#include "Francis.h"
+
+#endif
+

lib/algorithms/densematrix/Francis.h

@@ -0,0 +1,525 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/iterative/Francis.h
+
+    Copyright (C) 2015
+
+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 */
+#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

lib/algorithms/densematrix/Householder.h

@@ -0,0 +1,242 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/iterative/Householder.h
+
+    Copyright (C) 2015
+
+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 */
+#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

lib/algorithms/iterative/ConjugateGradientMixedPrec.h

@@ -60,6 +60,7 @@ namespace Grid {
     }
   
     void operator() (const FieldD &src_d_in, FieldD &sol_d){
+
       TotalInnerIterations = 0;
 	
       GridStopWatch TotalTimer;

lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h

@@ -0,0 +1,979 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Chulwoo Jung
+Author: Yong-Chull Jang <ypj@quark.phy.bnl.gov> 
+Author: Guido Cossu
+
+    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_IRBL_H
+#define GRID_IRBL_H
+
+#include <string.h> //memset
+
+#define Glog std::cout << GridLogMessage 
+
+namespace Grid {
+
+enum class LanczosType { irbl, rbl };
+
+/////////////////////////////////////////////////////////////
+// Implicitly restarted block lanczos
+/////////////////////////////////////////////////////////////
+template<class Field> 
+class ImplicitlyRestartedBlockLanczos {
+
+private:       
+  
+  std::string cname = std::string("ImplicitlyRestartedBlockLanczos");
+  int MaxIter;   // Max iterations
+  int Nstop;     // Number of evecs checked for convergence
+  int Nu;        // Numbeer of vecs in the unit block
+  int Nk;        // Number of converged sought
+  int Nm;        // total number of vectors
+  int Nblock_k;    // Nk/Nu
+  int Nblock_m;    // Nm/Nu
+  int Nconv_test_interval; // Number of skipped vectors when checking a convergence
+  RealD eresid;
+  IRLdiagonalisation diagonalisation;
+  ////////////////////////////////////
+  // Embedded objects
+  ////////////////////////////////////
+           SortEigen<Field> _sort;
+  LinearOperatorBase<Field> &_Linop;
+    OperatorFunction<Field> &_poly;
+
+  /////////////////////////
+  // Constructor
+  /////////////////////////
+public:       
+ ImplicitlyRestartedBlockLanczos(LinearOperatorBase<Field> &Linop, // op
+                                 OperatorFunction<Field> & poly,   // polynomial
+                                 int _Nstop, // really sought vecs
+                                 int _Nconv_test_interval, // conv check interval
+                                 int _Nu,    // vecs in the unit block
+                                 int _Nk,    // sought vecs
+                                 int _Nm,    // total vecs
+                                 RealD _eresid, // resid in lmd deficit 
+                                 int _MaxIter,  // Max iterations
+                                 IRLdiagonalisation _diagonalisation = IRLdiagonaliseWithEigen)
+   : _Linop(Linop),    _poly(poly),
+      Nstop(_Nstop), Nconv_test_interval(_Nconv_test_interval), 
+      Nu(_Nu), Nk(_Nk), Nm(_Nm), 
+      Nblock_m(_Nm/_Nu), Nblock_k(_Nk/_Nu),
+      //eresid(_eresid),  MaxIter(10),
+      eresid(_eresid),  MaxIter(_MaxIter),
+      diagonalisation(_diagonalisation)
+  { assert( (Nk%Nu==0) && (Nm%Nu==0) ); };
+
+  ////////////////////////////////
+  // Helpers
+  ////////////////////////////////
+  static RealD normalize(Field& v) 
+  {
+    RealD nn = norm2(v);
+    nn = sqrt(nn);
+    v = v * (1.0/nn);
+    return nn;
+  }
+  
+  void orthogonalize(Field& w, std::vector<Field>& evec, int k)
+  {
+    typedef typename Field::scalar_type MyComplex;
+    MyComplex ip;
+    
+    for(int j=0; j<k; ++j){
+      ip = innerProduct(evec[j],w); 
+      w = w - ip * evec[j];
+    }
+    normalize(w);
+  }
+  
+  void orthogonalize_blockhead(Field& w, std::vector<Field>& evec, int k, int Nu)
+  {
+    typedef typename Field::scalar_type MyComplex;
+    MyComplex ip;
+    
+    for(int j=0; j<k; ++j){
+      ip = innerProduct(evec[j*Nu],w); 
+      w = w - ip * evec[j*Nu];
+    }
+    normalize(w);
+  }
+  
+  void calc(std::vector<RealD>& eval,  
+            std::vector<Field>& evec, 
+            const std::vector<Field>& src, int& Nconv, LanczosType Impl)
+  {
+    switch (Impl) {
+      case LanczosType::irbl: 
+        calc_irbl(eval,evec,src,Nconv);
+        break;
+      
+      case LanczosType::rbl: 
+        calc_rbl(eval,evec,src,Nconv);
+        break;
+    }
+  }
+
+  void calc_irbl(std::vector<RealD>& eval,  
+                 std::vector<Field>& evec, 
+                 const std::vector<Field>& src, int& Nconv)
+  {
+    std::string fname = std::string(cname+"::calc_irbl()"); 
+    GridBase *grid = evec[0]._grid;
+    assert(grid == src[0]._grid);
+    assert( Nu = src.size() );
+    
+    Glog << std::string(74,'*') << std::endl;
+    Glog << fname + " starting iteration 0 /  "<< MaxIter<< std::endl;
+    Glog << std::string(74,'*') << std::endl;
+    Glog <<" -- seek   Nk    = "<< Nk    <<" vectors"<< std::endl;
+    Glog <<" -- accept Nstop = "<< Nstop <<" vectors"<< std::endl;
+    Glog <<" -- total  Nm    = "<< Nm    <<" vectors"<< std::endl;
+    Glog <<" -- size of eval = "<< eval.size() << std::endl;
+    Glog <<" -- size of evec = "<< evec.size() << std::endl;
+    if ( diagonalisation == IRLdiagonaliseWithEigen ) { 
+      Glog << "Diagonalisation is Eigen "<< std::endl;
+    } else {
+      abort();
+    }
+    Glog << std::string(74,'*') << std::endl;
+    
+    assert(Nm == evec.size() && Nm == eval.size());
+	
+    std::vector<std::vector<ComplexD>> lmd(Nu,std::vector<ComplexD>(Nm,0.0));  
+    std::vector<std::vector<ComplexD>> lme(Nu,std::vector<ComplexD>(Nm,0.0));  
+    std::vector<std::vector<ComplexD>> lmd2(Nu,std::vector<ComplexD>(Nm,0.0));  
+    std::vector<std::vector<ComplexD>> lme2(Nu,std::vector<ComplexD>(Nm,0.0));  
+    std::vector<RealD> eval2(Nm);
+    std::vector<RealD> resid(Nk);
+
+    Eigen::MatrixXcd    Qt = Eigen::MatrixXcd::Zero(Nm,Nm);
+    Eigen::MatrixXcd    Q = Eigen::MatrixXcd::Zero(Nm,Nm);
+
+    std::vector<int>   Iconv(Nm);
+    std::vector<Field>  B(Nm,grid); // waste of space replicating
+    
+    std::vector<Field> f(Nu,grid);
+    std::vector<Field> f_copy(Nu,grid);
+    Field v(grid);
+    
+    Nconv = 0;
+    
+    RealD beta_k;
+  
+    // set initial vector
+    for (int i=0; i<Nu; ++i) {
+      Glog << "norm2(src[" << i << "])= "<< norm2(src[i]) << std::endl;
+      evec[i] = src[i];
+      orthogonalize(evec[i],evec,i);
+      Glog << "norm2(evec[" << i << "])= "<< norm2(evec[i]) << std::endl;
+    }
+    
+    // initial Nblock_k steps
+    for(int b=0; b<Nblock_k; ++b) blockwiseStep(lmd,lme,evec,f,f_copy,b);
+
+    // restarting loop begins
+    int iter;
+    for(iter = 0; iter<MaxIter; ++iter){
+      
+      Glog <<"#Restart iteration = "<< iter << std::endl;
+      // additional (Nblock_m - Nblock_k) steps
+      for(int b=Nblock_k; b<Nblock_m; ++b) blockwiseStep(lmd,lme,evec,f,f_copy,b);
+      
+      // getting eigenvalues
+      for(int u=0; u<Nu; ++u){
+        for(int k=0; k<Nm; ++k){
+          lmd2[u][k] = lmd[u][k];
+          lme2[u][k] = lme[u][k];
+        }
+      }
+      Qt = Eigen::MatrixXcd::Identity(Nm,Nm);
+      diagonalize(eval2,lmd2,lme2,Nu,Nm,Nm,Qt,grid);
+      _sort.push(eval2,Nm);
+      Glog << "#Ritz value before shift: "<< std::endl;
+      for(int i=0; i<Nm; ++i){
+        std::cout.precision(13);
+        std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
+        std::cout << "Rval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[i] << std::endl;
+      }
+      
+      //----------------------------------------------------------------------
+      if ( Nm>Nk ) {
+        Glog <<" #Apply shifted QR transformations "<<std::endl;
+        //int k2 = Nk+Nu;
+        int k2 = Nk;
+      
+        Eigen::MatrixXcd BTDM = Eigen::MatrixXcd::Identity(Nm,Nm);
+        Q = Eigen::MatrixXcd::Identity(Nm,Nm);
+        
+        unpackHermitBlockTriDiagMatToEigen(lmd,lme,Nu,Nblock_m,Nm,Nm,BTDM);
+
+        for(int ip=Nk; ip<Nm; ++ip){ 
+          shiftedQRDecompEigen(BTDM,Nu,Nm,eval2[ip],Q);
+        }
+        
+        packHermitBlockTriDiagMatfromEigen(lmd,lme,Nu,Nblock_m,Nm,Nm,BTDM);
+
+        for(int i=0; i<k2; ++i) B[i] = 0.0;
+        for(int j=0; j<k2; ++j){
+          for(int k=0; k<Nm; ++k){
+            B[j].checkerboard = evec[k].checkerboard;
+            B[j] += evec[k]*Q(k,j);
+          }
+        }
+        for(int i=0; i<k2; ++i) evec[i] = B[i];
+
+        // reconstruct initial vector for additional pole space
+        blockwiseStep(lmd,lme,evec,f,f_copy,Nblock_k-1);
+
+        // getting eigenvalues
+        for(int u=0; u<Nu; ++u){
+          for(int k=0; k<Nm; ++k){
+            lmd2[u][k] = lmd[u][k];
+            lme2[u][k] = lme[u][k];
+          }
+        }
+        Qt = Eigen::MatrixXcd::Identity(Nm,Nm);
+        diagonalize(eval2,lmd2,lme2,Nu,Nk,Nm,Qt,grid);
+        _sort.push(eval2,Nk);
+        Glog << "#Ritz value after shift: "<< std::endl;
+        for(int i=0; i<Nk; ++i){
+          std::cout.precision(13);
+          std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
+          std::cout << "Rval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[i] << std::endl;
+        }
+      }
+      //----------------------------------------------------------------------
+
+      // Convergence test
+      Glog <<" #Convergence test: "<<std::endl;
+      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){
+	  B[j].checkerboard = evec[k].checkerboard;
+	  B[j] += evec[k]*Qt(k,j);
+	}
+      }
+      
+      Nconv = 0;
+      for(int i=0; i<Nk; ++i){
+	
+        _Linop.HermOp(B[i],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);
+        resid[i] = vv;
+	
+	std::cout.precision(13);
+        std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
+	std::cout << "eval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[i];
+	std::cout << "   resid^2 = "<< std::setw(20)<< std::setiosflags(std::ios_base::right)<< vv<< std::endl;
+	
+	// change the criteria as evals are supposed to be sorted, all evals smaller(larger) than Nstop should have converged
+	//if( (vv<eresid*eresid) && (i == Nconv) ){
+	if (vv<eresid*eresid) {
+	  Iconv[Nconv] = i;
+	  ++Nconv;
+	}
+	
+      }  // i-loop end
+      
+      Glog <<" #modes converged: "<<Nconv<<std::endl;
+      for(int i=0; i<Nconv; ++i){
+	std::cout.precision(13);
+        std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<Iconv[i]<<"] ";
+	std::cout << "eval_conv = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[Iconv[i]];
+	std::cout << "   resid^2 = "<< std::setw(20)<< std::setiosflags(std::ios_base::right)<< resid[Iconv[i]]<< std::endl;
+      } 
+
+      if ( Nconv>=Nstop ) break;
+
+    } // end of iter loop
+    
+    Glog << std::string(74,'*') << std::endl;
+    if ( Nconv<Nstop ) {
+      Glog << fname + " NOT converged ; Summary :\n";
+    } else {
+      Glog << fname + " CONVERGED ; Summary :\n";
+      // Sort convered eigenpairs.
+      eval.resize(Nconv);
+      evec.resize(Nconv,grid);
+      for(int i=0; i<Nconv; ++i){
+        eval[i] = eval2[Iconv[i]];
+        evec[i] = B[Iconv[i]];
+      }
+      _sort.push(eval,evec,Nconv);
+    }
+    Glog << std::string(74,'*') << std::endl;
+    Glog << " -- Iterations  = "<< iter   << "\n";
+    //Glog << " -- beta(k)     = "<< beta_k << "\n";
+    Glog << " -- Nconv       = "<< Nconv  << "\n";
+    Glog << std::string(74,'*') << std::endl;
+  
+  }
+  
+  
+  void calc_rbl(std::vector<RealD>& eval,  
+                 std::vector<Field>& evec, 
+                 const std::vector<Field>& src, int& Nconv)
+  {
+    std::string fname = std::string(cname+"::calc_rbl()"); 
+    GridBase *grid = evec[0]._grid;
+    assert(grid == src[0]._grid);
+    assert( Nu = src.size() );
+
+    int Np = (Nm-Nk);
+    if (Np > 0 && MaxIter > 1) Np /= MaxIter;
+    int Nblock_p = Np/Nu;
+    
+    Glog << std::string(74,'*') << std::endl;
+    Glog << fname + " starting iteration 0 /  "<< MaxIter<< std::endl;
+    Glog << std::string(74,'*') << std::endl;
+    Glog <<" -- seek (min) Nk    = "<< Nk    <<" vectors"<< std::endl;
+    Glog <<" -- seek (inc) Np    = "<< Np <<" vectors"<< std::endl;
+    Glog <<" -- seek (max) Nm    = "<< Nm    <<" vectors"<< std::endl;
+    Glog <<" -- accept Nstop     = "<< Nstop <<" vectors"<< std::endl;
+    Glog <<" -- size of eval     = "<< eval.size() << std::endl;
+    Glog <<" -- size of evec     = "<< evec.size() << std::endl;
+    if ( diagonalisation == IRLdiagonaliseWithEigen ) { 
+      Glog << "Diagonalisation is Eigen "<< std::endl;
+    } else {
+      abort();
+    }
+    Glog << std::string(74,'*') << std::endl;
+    
+    assert(Nm == evec.size() && Nm == eval.size());
+	
+    std::vector<std::vector<ComplexD>> lmd(Nu,std::vector<ComplexD>(Nm,0.0));  
+    std::vector<std::vector<ComplexD>> lme(Nu,std::vector<ComplexD>(Nm,0.0));  
+    std::vector<std::vector<ComplexD>> lmd2(Nu,std::vector<ComplexD>(Nm,0.0));  
+    std::vector<std::vector<ComplexD>> lme2(Nu,std::vector<ComplexD>(Nm,0.0));  
+    std::vector<RealD> eval2(Nk);
+    std::vector<RealD> resid(Nm);
+
+    Eigen::MatrixXcd    Qt = Eigen::MatrixXcd::Zero(Nm,Nm);
+    Eigen::MatrixXcd    Q = Eigen::MatrixXcd::Zero(Nm,Nm);
+
+    std::vector<int>   Iconv(Nm);
+    std::vector<Field>  B(Nm,grid); // waste of space replicating
+    
+    std::vector<Field> f(Nu,grid);
+    std::vector<Field> f_copy(Nu,grid);
+    Field v(grid);
+    
+    Nconv = 0;
+    
+    RealD beta_k;
+  
+    // set initial vector
+    for (int i=0; i<Nu; ++i) {
+      Glog << "norm2(src[" << i << "])= "<< norm2(src[i]) << std::endl;
+      evec[i] = src[i];
+      orthogonalize(evec[i],evec,i);
+      Glog << "norm2(evec[" << i << "])= "<< norm2(evec[i]) << std::endl;
+    }
+    
+    // initial Nblock_k steps
+    for(int b=0; b<Nblock_k; ++b) blockwiseStep(lmd,lme,evec,f,f_copy,b);
+
+    // restarting loop begins
+    int iter;
+    int Nblock_l, Nblock_r;
+    int Nl, Nr;
+    int Nconv_guess = 0;
+
+    for(iter = 0; iter<MaxIter; ++iter){
+         
+      Glog <<"#Restart iteration = "<< iter << std::endl;
+      
+      Nblock_l = Nblock_k + iter*Nblock_p;
+      Nblock_r = Nblock_l + Nblock_p;
+      Nl = Nblock_l*Nu;
+      Nr = Nblock_r*Nu;
+      eval2.resize(Nr);
+
+      // additional Nblock_p steps
+      for(int b=Nblock_l; b<Nblock_r; ++b) blockwiseStep(lmd,lme,evec,f,f_copy,b);
+      
+      // getting eigenvalues
+      for(int u=0; u<Nu; ++u){
+        for(int k=0; k<Nr; ++k){
+          lmd2[u][k] = lmd[u][k];
+          lme2[u][k] = lme[u][k];
+        }
+      }
+      Qt = Eigen::MatrixXcd::Identity(Nr,Nr);
+      diagonalize(eval2,lmd2,lme2,Nu,Nr,Nr,Qt,grid);
+      _sort.push(eval2,Nr);
+      Glog << "#Ritz value: "<< std::endl;
+      for(int i=0; i<Nr; ++i){
+        std::cout.precision(13);
+        std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
+        std::cout << "Rval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[i] << std::endl;
+      }
+      
+      // Convergence test
+      Glog <<" #Convergence test: "<<std::endl;
+      Nconv = 0;
+      for(int k = 0; k<Nr; ++k) B[k]=0.0;
+      for(int j = 0; j<Nr; j+=Nconv_test_interval){
+        if ( j/Nconv_test_interval == Nconv ) {
+          Glog <<" #rotation for next check point evec" 
+               << std::setw(4)<< std::setiosflags(std::ios_base::right) 
+               << "["<< j <<"]" <<std::endl;
+          for(int k = 0; k<Nr; ++k){
+            B[j].checkerboard = evec[k].checkerboard;
+            B[j] += evec[k]*Qt(k,j);
+          }
+          
+          _Linop.HermOp(B[j],v);
+          RealD vnum = real(innerProduct(B[j],v)); // HermOp.
+          RealD vden = norm2(B[j]);
+          eval2[j] = vnum/vden;
+          v -= eval2[j]*B[j];
+          RealD vv = norm2(v);
+          resid[j] = vv;
+          
+          std::cout.precision(13);
+          std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<j<<"] ";
+          std::cout << "eval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[j];
+          std::cout << "   resid^2 = "<< std::setw(20)<< std::setiosflags(std::ios_base::right)<< vv<< std::endl;
+          
+          // change the criteria as evals are supposed to be sorted, all evals smaller(larger) than Nstop should have converged
+          //if( (vv<eresid*eresid) && (i == Nconv) ){
+          if (vv<eresid*eresid) {
+            Iconv[Nconv] = j;
+            ++Nconv;
+          }
+        } else {
+          break;
+        }
+      }  // j-loop end
+      
+      Glog <<" #modes converged: "<<Nconv<<std::endl;
+      for(int i=0; i<Nconv; ++i){
+	std::cout.precision(13);
+        std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<Iconv[i]<<"] ";
+	std::cout << "eval_conv = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[Iconv[i]];
+	std::cout << "   resid^2 = "<< std::setw(20)<< std::setiosflags(std::ios_base::right)<< resid[Iconv[i]]<< std::endl;
+      } 
+
+      (Nconv > 0 ) ? Nconv_guess = 1 + (Nconv-1)*Nconv_test_interval : Nconv_guess = 0;
+      if ( Nconv_guess >= Nstop ) break;
+
+    } // end of iter loop
+    
+    Glog << std::string(74,'*') << std::endl;
+    if ( Nconv_guess < Nstop ) {
+      Glog << fname + " NOT converged ; Summary :\n";
+    } else {
+      Glog << fname + " CONVERGED ; Summary :\n";
+      // Sort convered eigenpairs.
+      eval.resize(Nconv);
+      evec.resize(Nconv,grid);
+      for(int i=0; i<Nconv; ++i){
+        eval[i] = eval2[Iconv[i]];
+        evec[i] = B[Iconv[i]];
+      }
+      _sort.push(eval,evec,Nconv);
+    }
+    Glog << std::string(74,'*') << std::endl;
+    Glog << " -- Iterations    = "<< iter   << "\n";
+    //Glog << " -- beta(k)       = "<< beta_k << "\n";
+    Glog << " -- Nconv         = "<< Nconv  << "\n";
+    Glog << " -- Nconv (guess) = "<< Nconv_guess  << "\n";
+    Glog << std::string(74,'*') << std::endl;
+  
+  }
+
+private:
+  void blockwiseStep(std::vector<std::vector<ComplexD>>& lmd,
+	             std::vector<std::vector<ComplexD>>& lme, 
+	             std::vector<Field>& evec,
+	             std::vector<Field>& w, 
+	             std::vector<Field>& w_copy, 
+                     int b)
+  {
+    const RealD tiny = 1.0e-20;
+    
+    int Nu = w.size();
+    int Nm = evec.size();
+    assert( b < Nm/Nu );
+    
+    // converts block index to full indicies for an interval [L,R)
+    int L = Nu*b;
+    int R = Nu*(b+1);
+
+    Real beta;
+    
+    // 3. wk:=Avk−βkv_{k−1}
+    for (int k=L, u=0; k<R; ++k, ++u) {
+      _poly(_Linop,evec[k],w[u]);      
+    }
+    
+    if (b>0) {
+      for (int u=0; u<Nu; ++u) {
+        //for (int k=L-Nu; k<L; ++k) {
+        for (int k=L-Nu+u; k<L; ++k) {
+          w[u] = w[u] - evec[k] * conjugate(lme[u][k]);
+        }
+      }
+    }
+    
+    // 4. αk:=(vk,wk)
+    //for (int u=0; u<Nu; ++u) {
+    //  for (int k=L; k<R; ++k) {
+    //    lmd[u][k] = innerProduct(evec[k],w[u]);  // lmd = transpose of alpha
+    //  }
+    //  lmd[u][L+u] = real(lmd[u][L+u]);  // force diagonal to be real
+    //}
+    for (int u=0; u<Nu; ++u) {
+      for (int k=L+u; k<R; ++k) {
+        lmd[u][k] = innerProduct(evec[k],w[u]);  // lmd = transpose of alpha
+        lmd[k-L][u+L] = conjugate(lmd[u][k]);     // force hermicity
+      }
+      lmd[u][L+u] = real(lmd[u][L+u]);  // force diagonal to be real
+    }
+    
+    // 5. wk:=wk−αkvk
+    for (int u=0; u<Nu; ++u) {
+      for (int k=L; k<R; ++k) {
+        w[u] = w[u] - evec[k]*lmd[u][k];
+      }
+      w_copy[u] = w[u];
+    }
+    
+    // In block version, the steps 6 and 7 in Lanczos construction is
+    // replaced by the QR decomposition of new basis block.
+    // It results block version beta and orthonormal block basis. 
+    // Here, QR decomposition is done by using Gram-Schmidt.
+    for (int u=0; u<Nu; ++u) {
+      for (int k=L; k<R; ++k) {
+        lme[u][k] = 0.0;
+      }
+    }
+
+#if 0
+    beta = normalize(w[0]);
+    for (int u=1; u<Nu; ++u) {
+      //orthogonalize(w[u],w_copy,u);
+      orthogonalize(w[u],w,u);
+    }
+#else
+    // re-orthogonalization for numerical stability
+    for (int u=0; u<Nu; ++u) {
+      orthogonalize(w[u],evec,R);
+    }
+    // QR part
+    for (int u=1; u<Nu; ++u) {
+      orthogonalize(w[u],w,u);
+    }
+#endif
+    
+    for (int u=0; u<Nu; ++u) {
+      //for (int v=0; v<Nu; ++v) {
+      for (int v=u; v<Nu; ++v) {
+        lme[u][L+v] = innerProduct(w[u],w_copy[v]);
+      }
+      lme[u][L+u] = real(lme[u][L+u]);  // force diagonal to be real
+    }
+    //lme[0][L] = beta;
+    
+    for (int u=0; u<Nu; ++u) {
+      Glog << "norm2(w[" << u << "])= "<< norm2(w[u]) << std::endl;
+      for (int k=L+u; k<R; ++k) {
+        Glog <<" In block "<< b << ","; 
+        std::cout <<" beta[" << u << "," << k-L << "] = ";
+        std::cout << lme[u][k] << std::endl;
+      }
+    }
+#if 0    
+    // re-orthogonalization for numerical stability
+    if (b>0) {
+      for (int u=0; u<Nu; ++u) {
+        orthogonalize(w[u],evec,R);
+      }
+      for (int u=1; u<Nu; ++u) {
+        orthogonalize(w[u],w,u);
+      }
+    }
+    //if (b>0) {
+    //  orthogonalize_blockhead(w[0],evec,b,Nu);
+    //  for (int u=1; u<Nu; ++u) {
+    //    orthogonalize(w[u],w,u);
+    //  }
+    //}
+#endif
+
+    if (b < Nm/Nu-1) {
+      for (int u=0; u<Nu; ++u) {
+        evec[R+u] = w[u];
+      }
+    }
+
+  }
+  
+    
+  void diagonalize_Eigen(std::vector<RealD>& eval, 
+                         std::vector<std::vector<ComplexD>>& lmd,
+                         std::vector<std::vector<ComplexD>>& lme, 
+			 int Nu, int Nk, int Nm,
+			 Eigen::MatrixXcd & Qt, // Nm x Nm
+			 GridBase *grid)
+  {
+    assert( Nk%Nu == 0 && Nm%Nu == 0 );
+    assert( Nk <= Nm );
+    Eigen::MatrixXcd BlockTriDiag = Eigen::MatrixXcd::Zero(Nk,Nk);
+    
+    for ( int u=0; u<Nu; ++u ) {
+      for (int k=0; k<Nk; ++k ) {
+        BlockTriDiag(k,u+(k/Nu)*Nu) = lmd[u][k];
+      }
+    }
+    
+    for ( int u=0; u<Nu; ++u ) {
+      for (int k=Nu; k<Nk; ++k ) {
+        BlockTriDiag(k-Nu,u+(k/Nu)*Nu) = conjugate(lme[u][k-Nu]);
+        BlockTriDiag(u+(k/Nu)*Nu,k-Nu) = lme[u][k-Nu];
+      }
+    }
+    //std::cout << BlockTriDiag << std::endl;
+    
+    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXcd> eigensolver(BlockTriDiag);
+
+    for (int i = 0; i < Nk; i++) {
+      eval[Nk-1-i] = eigensolver.eigenvalues()(i);
+    }
+    for (int i = 0; i < Nk; i++) {
+      for (int j = 0; j < Nk; j++) {
+	Qt(j,Nk-1-i) = eigensolver.eigenvectors()(j,i);
+	//Qt(Nk-1-i,j) = eigensolver.eigenvectors()(i,j);
+	//Qt(i,j) = eigensolver.eigenvectors()(i,j);
+      }
+    }
+  }
+
+
+  void diagonalize(std::vector<RealD>& eval, 
+                   std::vector<std::vector<ComplexD>>& lmd, 
+                   std::vector<std::vector<ComplexD>>& lme, 
+		   int Nu, int Nk, int Nm,   
+		   Eigen::MatrixXcd & Qt,
+		   GridBase *grid)
+  {
+    Qt = Eigen::MatrixXcd::Identity(Nm,Nm);
+    if ( diagonalisation == IRLdiagonaliseWithEigen ) { 
+      diagonalize_Eigen(eval,lmd,lme,Nu,Nk,Nm,Qt,grid);
+    } else { 
+      assert(0);
+    }
+  }
+  
+
+  void unpackHermitBlockTriDiagMatToEigen(
+         std::vector<std::vector<ComplexD>>& lmd,  
+         std::vector<std::vector<ComplexD>>& lme,
+         int Nu, int Nb, int Nk, int Nm,
+         Eigen::MatrixXcd& M)
+  {
+    //Glog << "unpackHermitBlockTriDiagMatToEigen() begin" << '\n'; 
+    assert( Nk%Nu == 0 && Nm%Nu == 0 );
+    assert( Nk <= Nm );
+    M = Eigen::MatrixXcd::Zero(Nk,Nk);
+    
+    // rearrange 
+    for ( int u=0; u<Nu; ++u ) {
+      for (int k=0; k<Nk; ++k ) {
+        M(k,u+(k/Nu)*Nu) = lmd[u][k];
+      }
+    }
+
+    for ( int u=0; u<Nu; ++u ) {
+      for (int k=Nu; k<Nk; ++k ) {
+        M(k-Nu,u+(k/Nu)*Nu) = conjugate(lme[u][k-Nu]);
+        M(u+(k/Nu)*Nu,k-Nu) = lme[u][k-Nu];
+      }
+    }
+    //Glog << "unpackHermitBlockTriDiagMatToEigen() end" << endl; 
+  }
+ 
+
+  void packHermitBlockTriDiagMatfromEigen(
+         std::vector<std::vector<ComplexD>>& lmd,
+         std::vector<std::vector<ComplexD>>& lme,
+         int Nu, int Nb, int Nk, int Nm,
+         Eigen::MatrixXcd& M)
+  {
+    //Glog << "packHermitBlockTriDiagMatfromEigen() begin" << '\n'; 
+    assert( Nk%Nu == 0 && Nm%Nu == 0 );
+    assert( Nk <= Nm );
+    
+    // rearrange 
+    for ( int u=0; u<Nu; ++u ) {
+      for (int k=0; k<Nk; ++k ) {
+        lmd[u][k] = M(k,u+(k/Nu)*Nu);
+      }
+    }
+
+    for ( int u=0; u<Nu; ++u ) {
+      for (int k=Nu; k<Nk; ++k ) {
+        lme[u][k-Nu] = M(u+(k/Nu)*Nu,k-Nu);
+      }
+    }
+    //Glog << "packHermitBlockTriDiagMatfromEigen() end" << endl; 
+  }
+
+
+  // assume the input matrix M is a band matrix
+  void shiftedQRDecompEigen(Eigen::MatrixXcd& M, int Nu, int Nm,
+		            RealD Dsh,
+		            Eigen::MatrixXcd& Qprod)
+  {
+    //Glog << "shiftedQRDecompEigen() begin" << '\n'; 
+    Eigen::MatrixXcd Q = Eigen::MatrixXcd::Zero(Nm,Nm);
+    Eigen::MatrixXcd R = Eigen::MatrixXcd::Zero(Nm,Nm);
+    Eigen::MatrixXcd Mtmp = Eigen::MatrixXcd::Zero(Nm,Nm);
+    
+    Mtmp = M;
+    for (int i=0; i<Nm; ++i ) {
+      Mtmp(i,i) = M(i,i) - Dsh;
+    }
+    
+    Eigen::HouseholderQR<Eigen::MatrixXcd> QRD(Mtmp);
+    Q = QRD.householderQ();
+    R = QRD.matrixQR(); // upper triangular part is the R matrix.
+                        // lower triangular part used to represent series
+                        // of Q sequence.
+
+    // equivalent operation of Qprod *= Q
+    //M = Eigen::MatrixXcd::Zero(Nm,Nm);
+    
+    //for (int i=0; i<Nm; ++i) {
+    //  for (int j=0; j<Nm-2*(Nu+1); ++j) {
+    //    for (int k=0; k<2*(Nu+1)+j; ++k) {
+    //      M(i,j) += Qprod(i,k)*Q(k,j);
+    //    }
+    //  }
+    //}
+    //for (int i=0; i<Nm; ++i) {
+    //  for (int j=Nm-2*(Nu+1); j<Nm; ++j) {
+    //    for (int k=0; k<Nm; ++k) {
+    //      M(i,j) += Qprod(i,k)*Q(k,j);
+    //    }
+    //  }
+    //}
+    
+    Mtmp = Eigen::MatrixXcd::Zero(Nm,Nm);
+
+    for (int i=0; i<Nm; ++i) {
+      for (int j=0; j<Nm-(Nu+1); ++j) {
+        for (int k=0; k<Nu+1+j; ++k) {
+          Mtmp(i,j) += Qprod(i,k)*Q(k,j);
+        }
+      }
+    }
+    for (int i=0; i<Nm; ++i) {
+      for (int j=Nm-(Nu+1); j<Nm; ++j) {
+        for (int k=0; k<Nm; ++k) {
+          Mtmp(i,j) += Qprod(i,k)*Q(k,j);
+        }
+      }
+    }
+    
+    //static int ntimes = 2;
+    //for (int j=0; j<Nm-(ntimes*Nu); ++j) {
+    //  for (int i=ntimes*Nu+j; i<Nm; ++i) {
+    //    Mtmp(i,j) = 0.0;
+    //  }
+    //}
+    //ntimes++;
+
+    Qprod = Mtmp;
+     
+    // equivalent operation of M = Q.adjoint()*(M*Q)
+    Mtmp = Eigen::MatrixXcd::Zero(Nm,Nm);
+    
+    for (int a=0, i=0, kmax=0; a<Nu+1; ++a) {
+      for (int j=0; j<Nm-a; ++j) {
+        i = j+a;
+        kmax = (Nu+1)+j;
+        if (kmax > Nm) kmax = Nm;
+        for (int k=i; k<kmax; ++k) { 
+          Mtmp(i,j) += R(i,k)*Q(k,j);
+        }
+        Mtmp(j,i) = conj(Mtmp(i,j));
+      }
+    }
+
+    for (int i=0; i<Nm; ++i) {
+      Mtmp(i,i) = real(Mtmp(i,i)) + Dsh;
+    }
+    
+    M = Mtmp;
+
+    //M = Q.adjoint()*(M*Q);
+    //for (int i=0; i<Nm; ++i) {
+    //  for (int j=0; j<Nm; ++j) {
+    //    if (i==j) M(i,i) = real(M(i,i));
+    //    if (j>i)  M(i,j) = conj(M(j,i));
+    //    if (i-j > Nu || j-i > Nu) M(i,j) = 0.;
+    //  }
+    //}
+    
+    //Glog << "shiftedQRDecompEigen() end" << endl; 
+  }
+
+  void exampleQRDecompEigen(void)
+  {
+    Eigen::MatrixXd A = Eigen::MatrixXd::Zero(3,3);
+    Eigen::MatrixXd Q = Eigen::MatrixXd::Zero(3,3);
+    Eigen::MatrixXd R = Eigen::MatrixXd::Zero(3,3);
+    Eigen::MatrixXd P = Eigen::MatrixXd::Zero(3,3);
+
+    A(0,0) = 12.0;
+    A(0,1) = -51.0;
+    A(0,2) = 4.0;
+    A(1,0) = 6.0;
+    A(1,1) = 167.0;
+    A(1,2) = -68.0;
+    A(2,0) = -4.0;
+    A(2,1) = 24.0;
+    A(2,2) = -41.0;
+    
+    Glog << "matrix A before ColPivHouseholder" << std::endl;
+    for ( int i=0; i<3; i++ ) {
+      for ( int j=0; j<3; j++ ) {
+        Glog << "A[" << i << "," << j << "] = " << A(i,j) << '\n';
+      }
+    }
+    Glog << std::endl;
+
+    Eigen::ColPivHouseholderQR<Eigen::MatrixXd> QRD(A);
+    
+    Glog << "matrix A after ColPivHouseholder" << std::endl;
+    for ( int i=0; i<3; i++ ) {
+      for ( int j=0; j<3; j++ ) {
+        Glog << "A[" << i << "," << j << "] = " << A(i,j) << '\n';
+      }
+    }
+    Glog << std::endl;
+    
+    Glog << "HouseholderQ with sequence lenth = nonzeroPiviots" << std::endl;
+    Q = QRD.householderQ().setLength(QRD.nonzeroPivots());
+    for ( int i=0; i<3; i++ ) {
+      for ( int j=0; j<3; j++ ) {
+        Glog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
+      }
+    }
+    Glog << std::endl;
+    
+    Glog << "HouseholderQ with sequence lenth = 1" << std::endl;
+    Q = QRD.householderQ().setLength(1);
+    for ( int i=0; i<3; i++ ) {
+      for ( int j=0; j<3; j++ ) {
+        Glog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
+      }
+    }
+    Glog << std::endl;
+    
+    Glog << "HouseholderQ with sequence lenth = 2" << std::endl;
+    Q = QRD.householderQ().setLength(2);
+    for ( int i=0; i<3; i++ ) {
+      for ( int j=0; j<3; j++ ) {
+        Glog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
+      }
+    }
+    Glog << std::endl;
+    
+    Glog << "matrixR" << std::endl;
+    R = QRD.matrixR();
+    for ( int i=0; i<3; i++ ) {
+      for ( int j=0; j<3; j++ ) {
+        Glog << "R[" << i << "," << j << "] = " << R(i,j) << '\n';
+      }
+    }
+    Glog << std::endl;
+
+    Glog << "rank = " << QRD.rank() << std::endl;
+    Glog << "threshold = " << QRD.threshold() << std::endl;
+    
+    Glog << "matrixP" << std::endl;
+    P = QRD.colsPermutation();
+    for ( int i=0; i<3; i++ ) {
+      for ( int j=0; j<3; j++ ) {
+        Glog << "P[" << i << "," << j << "] = " << P(i,j) << '\n';
+      }
+    }
+    Glog << std::endl;
+
+
+    Glog << "QR decomposition without column pivoting" << std::endl;
+    
+    A(0,0) = 12.0;
+    A(0,1) = -51.0;
+    A(0,2) = 4.0;
+    A(1,0) = 6.0;
+    A(1,1) = 167.0;
+    A(1,2) = -68.0;
+    A(2,0) = -4.0;
+    A(2,1) = 24.0;
+    A(2,2) = -41.0;
+    
+    Glog << "matrix A before Householder" << std::endl;
+    for ( int i=0; i<3; i++ ) {
+      for ( int j=0; j<3; j++ ) {
+        Glog << "A[" << i << "," << j << "] = " << A(i,j) << '\n';
+      }
+    }
+    Glog << std::endl;
+    
+    Eigen::HouseholderQR<Eigen::MatrixXd> QRDplain(A);
+    
+    Glog << "HouseholderQ" << std::endl;
+    Q = QRDplain.householderQ();
+    for ( int i=0; i<3; i++ ) {
+      for ( int j=0; j<3; j++ ) {
+        Glog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
+      }
+    }
+    Glog << std::endl;
+    
+    Glog << "matrix A after Householder" << std::endl;
+    for ( int i=0; i<3; i++ ) {
+      for ( int j=0; j<3; j++ ) {
+        Glog << "A[" << i << "," << j << "] = " << A(i,j) << '\n';
+      }
+    }
+    Glog << std::endl;
+  }
+
+ };
+}
+#undef Glog
+#endif

lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h.bak

@@ -0,0 +1,835 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Chulwoo Jung
+Author: Yong-Chull Jang <ypj@quark.phy.bnl.gov> 
+Author: Guido Cossu
+
+    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_IRBL_H
+#define GRID_IRBL_H
+
+#include <string.h> //memset
+
+#define clog std::cout << GridLogMessage 
+
+namespace Grid {
+
+/////////////////////////////////////////////////////////////
+// Implicitly restarted block lanczos
+/////////////////////////////////////////////////////////////
+template<class Field> 
+class ImplicitlyRestartedBlockLanczos {
+
+private:       
+  
+  std::string cname = std::string("ImplicitlyRestartedBlockLanczos");
+  int MaxIter;   // Max iterations
+  int Nstop;     // Number of evecs checked for convergence
+  int Nu;        // Numbeer of vecs in the unit block
+  int Nk;        // Number of converged sought
+  int Nm;        // total number of vectors
+  int Nblock_k;    // Nk/Nu
+  int Nblock_m;    // Nm/Nu
+  RealD eresid;
+  IRLdiagonalisation diagonalisation;
+  ////////////////////////////////////
+  // Embedded objects
+  ////////////////////////////////////
+           SortEigen<Field> _sort;
+  LinearOperatorBase<Field> &_Linop;
+    OperatorFunction<Field> &_poly;
+
+  /////////////////////////
+  // Constructor
+  /////////////////////////
+public:       
+ ImplicitlyRestartedBlockLanczos(LinearOperatorBase<Field> &Linop, // op
+                                 OperatorFunction<Field> & poly,   // polynomial
+                                 int _Nstop, // really sought vecs
+                                 int _Nu,    // vecs in the unit block
+                                 int _Nk,    // sought vecs
+                                 int _Nm,    // total vecs
+                                 RealD _eresid, // resid in lmd deficit 
+                                 int _MaxIter,  // Max iterations
+                                 IRLdiagonalisation _diagonalisation = IRLdiagonaliseWithEigen)
+   : _Linop(Linop),    _poly(poly),
+      Nstop(_Nstop), Nu(_Nu), Nk(_Nk), Nm(_Nm), 
+      Nblock_m(_Nm/_Nu), Nblock_k(_Nk/_Nu),
+      //eresid(_eresid),  MaxIter(10),
+      eresid(_eresid),  MaxIter(_MaxIter),
+      diagonalisation(_diagonalisation)
+  { assert( (Nk%Nu==0) && (Nm%Nu==0) ); };
+
+  ////////////////////////////////
+  // Helpers
+  ////////////////////////////////
+  static RealD normalize(Field& v) 
+  {
+    RealD nn = norm2(v);
+    nn = sqrt(nn);
+    v = v * (1.0/nn);
+    return nn;
+  }
+  
+  void orthogonalize(Field& w, std::vector<Field>& evec, int k)
+  {
+    typedef typename Field::scalar_type MyComplex;
+    MyComplex ip;
+    
+    for(int j=0; j<k; ++j){
+      ip = innerProduct(evec[j],w); 
+      w = w - ip * evec[j];
+    }
+    normalize(w);
+  }
+
+/* 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(std::vector<RealD>& eval,  
+            std::vector<Field>& evec, 
+            const std::vector<Field>& src, int& Nconv)
+  {
+    std::string fname = std::string(cname+"::calc()"); 
+    GridBase *grid = evec[0]._grid;
+    assert(grid == src[0]._grid);
+    assert( Nu = src.size() );
+    
+    clog << std::string(74,'*') << std::endl;
+    clog << fname + " starting iteration 0 /  "<< MaxIter<< std::endl;
+    clog << std::string(74,'*') << std::endl;
+    clog <<" -- seek   Nk    = "<< Nk    <<" vectors"<< std::endl;
+    clog <<" -- accept Nstop = "<< Nstop <<" vectors"<< std::endl;
+    clog <<" -- total  Nm    = "<< Nm    <<" vectors"<< std::endl;
+    clog <<" -- size of eval = "<< eval.size() << std::endl;
+    clog <<" -- size of evec = "<< evec.size() << std::endl;
+    if ( diagonalisation == IRLdiagonaliseWithEigen ) { 
+      clog << "Diagonalisation is Eigen "<< std::endl;
+    } else {
+      abort();
+    }
+    clog << std::string(74,'*') << std::endl;
+    
+    assert(Nm == evec.size() && Nm == eval.size());
+	
+    std::vector<std::vector<ComplexD>> lmd(Nu,std::vector<ComplexD>(Nm,0.0));  
+    std::vector<std::vector<ComplexD>> lme(Nu,std::vector<ComplexD>(Nm,0.0));  
+    std::vector<std::vector<ComplexD>> lmd2(Nu,std::vector<ComplexD>(Nm,0.0));  
+    std::vector<std::vector<ComplexD>> lme2(Nu,std::vector<ComplexD>(Nm,0.0));  
+    std::vector<RealD> eval2(Nm);
+
+    Eigen::MatrixXcd    Qt = Eigen::MatrixXcd::Zero(Nm,Nm);
+    Eigen::MatrixXcd    Q = Eigen::MatrixXcd::Zero(Nm,Nm);
+
+    std::vector<int>   Iconv(Nm);
+    std::vector<Field>  B(Nm,grid); // waste of space replicating
+    
+    std::vector<Field> f(Nu,grid);
+    std::vector<Field> f_copy(Nu,grid);
+    Field v(grid);
+    
+    Nconv = 0;
+    
+    RealD beta_k;
+  
+    // set initial vector
+    for (int i=0; i<Nu; ++i) {
+      clog << "norm2(src[" << i << "])= "<< norm2(src[i]) << std::endl;
+      evec[i] = src[i];
+      orthogonalize(evec[i],evec,i);
+      clog << "norm2(evec[" << i << "])= "<< norm2(evec[i]) << std::endl;
+    }
+    
+    // initial Nblock_k steps
+    for(int b=0; b<Nblock_k; ++b) blockwiseStep(lmd,lme,evec,f,f_copy,b);
+
+    // restarting loop begins
+    int iter;
+    for(iter = 0; iter<MaxIter; ++iter){
+      
+      clog <<" **********************"<< std::endl;
+      clog <<" Restart iteration = "<< iter << std::endl;
+      clog <<" **********************"<< std::endl;
+      
+      // additional (Nblock_m - Nblock_k) steps
+      for(int b=Nblock_k; b<Nblock_m; ++b) blockwiseStep(lmd,lme,evec,f,f_copy,b);
+      
+      for(int k=0; k<Nm; ++k) {
+        clog << "ckpt A1: lme[" << k << "] = " << lme[0][k] << '\n';
+      }
+      for(int k=0; k<Nm; ++k) {
+        clog << "ckpt A2: lmd[" << k << "] = " << lmd[0][k] << '\n';
+      }
+      
+      // residual vector
+#if 1 // ypj[fixme] temporary to check a case when block has one vector
+      for ( int i=0; i<Nu; ++i) f_copy[i] = f[i];
+      for ( int i=0; i<Nu; ++i) {
+        f[i] = f_copy[0]*lme[0][Nm-Nu+i]; 
+        for ( int j=1; j<Nu; ++j) { 
+          f[i] += f_copy[j]*lme[j][Nm-Nu+i]; 
+        }
+        //clog << "ckpt C (i= " << i << ")" << '\n';
+        //clog << "norm2(f) = " << norm2(f[i]) << std::endl;
+      }
+#endif
+      
+      // getting eigenvalues
+      for(int u=0; u<Nu; ++u){
+        for(int k=0; k<Nm; ++k){
+          lmd2[u][k] = lmd[u][k];
+          lme2[u][k] = lme[u][k];
+        }
+      }
+      Qt = Eigen::MatrixXcd::Identity(Nm,Nm);
+      diagonalize(eval2,lmd2,lme2,Nu,Nm,Nm,Qt,grid);
+      
+      //for(int k=0; k<Nm; ++k){
+      //  clog << "ckpt D " << '\n';
+      //  clog << "eval2 [" << k << "] = " << eval2[k] << std::endl;
+      //}
+
+      // sorting
+      _sort.push(eval2,Nm);
+      
+      //for(int k=0; k<Nm; ++k){
+      //  clog << "ckpt E " << '\n';
+      //  clog << "eval2 [" << k << "] = " << eval2[k] << std::endl;
+      //}
+
+      // Implicitly shifted QR transformations
+      Eigen::MatrixXcd BTDM = Eigen::MatrixXcd::Identity(Nm,Nm);
+      Q = Eigen::MatrixXcd::Identity(Nm,Nm);
+      
+      unpackHermitBlockTriDiagMatToEigen(lmd,lme,Nu,Nblock_m,Nm,Nm,BTDM);
+
+      for(int ip=Nk; ip<Nm; ++ip){ 
+        clog << "ckpt B1: shift[" << ip << "] = " << eval2[ip] << endl;
+        shiftedQRDecompEigen(BTDM,Nm,eval2[ip],Q);
+      }
+      BTDM = Q.adjoint()*(BTDM*Q);
+      for (int i=0; i<Nm; ++i ) {
+        for (int j=i+1; j<Nm; ++j ) {
+          BTDM(i,j) = BTDM(j,i);
+        }
+        //BTDM(i,i) = real(BTDM(i,i));
+      }
+      
+      packHermitBlockTriDiagMatfromEigen(lmd,lme,Nu,Nblock_m,Nm,Nm,BTDM);
+
+      //for (int i=0; i<Nm; ++i) {
+      //  for (int j=0; j<Nm; ++j) {
+      //    clog << "ckpt G1: M[" << i << "," << j << "] = " << BTDM(i,j) << '\n';
+      //  }
+      //}
+      //for (int i=0; i<Nm; ++i) {
+      //  for (int j=0; j<Nm; ++j) {
+      //    clog << "ckpt G2: Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
+      //  }
+      //}
+      for (int i=0; i<Nm; ++i) {
+        clog << "ckpt C1: lme[" << i << "] = " << lme[0][i] << '\n';
+      }
+      for (int i=0; i<Nm; ++i) {
+        clog << "ckpt C2: lmd[" << i << "] = " << lmd[0][i] << '\n';
+      }
+
+      for(int i=0; i<Nk+Nu; ++i) B[i] = 0.0;
+      for(int j=0; j<Nk+Nu; ++j){
+	for(int k=0; k<Nm; ++k){
+	  B[j].checkerboard = evec[k].checkerboard;
+	  B[j] += evec[k]*Q(k,j);
+	}
+      }
+      for(int i=0; i<Nk+Nu; ++i) { 
+        evec[i] = B[i];
+        //clog << "ckpt F: norm2_evec[= " << i << "]" << norm2(evec[i]) << std::endl;
+      }
+
+#if 1 // ypj[fixme] temporary to check a case when block has one vector
+      // Compressed vector f and beta(k2)
+      f[0] *= Q(Nm-1,Nk-1);
+      f[0] += lme[0][Nk-1] * evec[Nk]; // was commented out
+      std::cout<< GridLogMessage<<"ckpt D1: Q[Nm-1,Nk-1] = "<<Q(Nm-1,Nk-1)<<std::endl;
+      beta_k = norm2(f[0]);
+      beta_k = sqrt(beta_k);
+      std::cout<< GridLogMessage<<"ckpt D2: beta(k) = "<<beta_k<<std::endl;
+      
+      RealD betar = 1.0/beta_k;
+      evec[Nk] = betar * f[0];
+      lme[0][Nk-1] = beta_k;
+#endif
+
+      // Convergence test
+      for(int u=0; u<Nu; ++u){
+        for(int k=0; k<Nm; ++k){
+          lmd2[u][k] = lmd[u][k];
+          lme2[u][k] = lme[u][k];
+        }
+      }
+      Qt = Eigen::MatrixXcd::Identity(Nm,Nm);
+      diagonalize(eval2,lmd2,lme2,Nu,Nk,Nm,Qt,grid);
+      
+      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){
+	  B[j].checkerboard = evec[k].checkerboard;
+	  B[j] += evec[k]*Qt(k,j);
+	}
+      }
+      
+      //for (int i=0; i<Nk; ++i) {
+      //  for (int j=0; j<Nk; ++j) {
+      //    clog << "ckpt H1: R[" << i << "," << j << "] = " << Qt(i,j) << '\n';
+      //  }
+      //}
+      //for (int i=0; i<Nk; ++i) {
+      //  clog << "ckpt H2: eval2[" << i << "] = " << eval2[i] << '\n';
+      //}
+      
+      //for(int j=0; j<Nk; ++j) {
+      //  clog << "ckpt I: norm2_B[ " << j << "]" << norm2(B[j]) << std::endl;
+      //}
+
+      Nconv = 0;
+      for(int i=0; i<Nk; ++i){
+	
+	_Linop.HermOp(B[i],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.precision(13);
+	clog << "[" << 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;
+	
+	// change the criteria as evals are supposed to be sorted, all evals smaller(larger) than Nstop should have converged
+	if( (vv<eresid*eresid) && (i == Nconv) ){
+	//if( (vv<eresid*eresid) ){
+	  Iconv[Nconv] = i;
+	  ++Nconv;
+	}
+	
+      }  // i-loop end
+      
+      clog <<" #modes converged: "<<Nconv<<std::endl;
+
+      if( Nconv>=Nstop ){
+	goto converged;
+      }
+
+    } // end of iter loop
+    
+    clog <<"**************************************************************************"<< std::endl;
+    std::cout<< GridLogError    << fname + " NOT converged.";
+    clog <<"**************************************************************************"<< std::endl;
+    abort();
+	
+  converged:
+    // Sorting
+    eval.resize(Nconv);
+    evec.resize(Nconv,grid);
+    for(int i=0; i<Nconv; ++i){
+      eval[i] = eval2[Iconv[i]];
+      evec[i] = B[Iconv[i]];
+    }
+    _sort.push(eval,evec,Nconv);
+    
+    clog <<"**************************************************************************"<< std::endl;
+    clog << fname + " CONVERGED ; Summary :\n";
+    clog <<"**************************************************************************"<< std::endl;
+    clog << " -- Iterations  = "<< iter   << "\n";
+    clog << " -- beta(k)     = "<< beta_k << "\n";
+    clog << " -- Nconv       = "<< Nconv  << "\n";
+    clog <<"**************************************************************************"<< std::endl;
+  }
+
+private:
+/* 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 blockwiseStep(std::vector<std::vector<ComplexD>>& lmd,
+	             std::vector<std::vector<ComplexD>>& lme, 
+	             std::vector<Field>& evec,
+	             std::vector<Field>& w, 
+	             std::vector<Field>& w_copy, 
+                     int b)
+  {
+    const RealD tiny = 1.0e-20;
+    
+    int Nu = w.size();
+    int Nm = evec.size();
+    assert( b < Nm/Nu );
+    
+    // converts block index to full indicies for an interval [L,R)
+    int L = Nu*b;
+    int R = Nu*(b+1);
+
+    Real beta;
+    
+    // 3. wk:=Avk−βkv_{k−1}
+    for (int k=L, u=0; k<R; ++k, ++u) {
+      _poly(_Linop,evec[k],w[u]);      
+    }
+    
+    if (b>0) {
+      for (int u=0; u<Nu; ++u) {
+        for (int k=L-Nu; k<L; ++k) {
+          w[u] = w[u] - evec[k] * conjugate(lme[u][k]);
+          //clog << "ckpt A (k= " << k+1 << ")" << '\n';
+          //clog << "lme = " << lme[u][k] << '\n';
+          //clog << "lme = " << conjugate(lme[u][k]) << '\n';
+        }
+        //clog << "norm(w) = " << norm2(w[u]) << std::endl;
+      }
+    }
+    
+    // 4. αk:=(vk,wk)
+    for (int u=0; u<Nu; ++u) {
+      for (int k=L; k<R; ++k) {
+        lmd[u][k] = innerProduct(evec[k],w[u]);  // lmd = transpose of alpha
+      }
+      lmd[u][L+u] = real(lmd[u][L+u]);  // force diagonal to be real
+      //clog << "ckpt B (k= " << L+u << ")" << '\n';
+      //clog << "lmd = " << lmd[u][L+u] << std::endl;
+    }
+    
+    // 5. wk:=wk−αkvk
+    for (int u=0; u<Nu; ++u) {
+      for (int k=L; k<R; ++k) {
+        w[u] = w[u] - evec[k]*lmd[u][k];
+      }
+      w_copy[u] = w[u];
+    }
+    
+    // In block version, the steps 6 and 7 in Lanczos construction is
+    // replaced by the QR decomposition of new basis block.
+    // It results block version beta and orthonormal block basis. 
+    // Here, QR decomposition is done by using Gram-Schmidt
+    for (int u=0; u<Nu; ++u) {
+      for (int k=L; k<R; ++k) {
+        lme[u][k] = 0.0;
+      }
+    }
+
+    beta = normalize(w[0]);
+    for (int u=1; u<Nu; ++u) {
+      //orthogonalize(w[u],w_copy,u);
+      orthogonalize(w[u],w,u);
+    }
+    
+    for (int u=0; u<Nu; ++u) {
+      for (int v=0; v<Nu; ++v) {
+        lme[u][L+v] = innerProduct(w[u],w_copy[v]);
+      }
+    }
+    lme[0][L] = beta;
+    
+#if 0
+    for (int u=0; u<Nu; ++u) {
+      for (int k=L+u; k<R; ++k) {
+        if (lme[u][k] < tiny) {
+          clog <<" In block "<< b << ","; 
+          std::cout <<" beta[" << u << "," << k-L << "] = ";
+          std::cout << lme[u][k] << std::endl;
+        }
+      }
+    }
+#else
+    for (int u=0; u<Nu; ++u) {
+      clog << "norm2(w[" << u << "])= "<< norm2(w[u]) << std::endl;
+      for (int k=L+u; k<R; ++k) {
+        clog <<" In block "<< b << ","; 
+        std::cout <<" beta[" << u << "," << k-L << "] = ";
+        std::cout << lme[u][k] << std::endl;
+      }
+    }
+#endif
+    
+    // re-orthogonalization for numerical stability
+    if (b>0) {
+      for (int u=0; u<Nu; ++u) {
+        orthogonalize(w[u],evec,R);
+      }
+    }
+
+    if (b < Nm/Nu-1) {
+      for (int u=0; u<Nu; ++u) {
+        evec[R+u] = w[u];
+      }
+    }
+
+  }
+  
+
+  void diagonalize_Eigen(std::vector<RealD>& eval, 
+                         std::vector<std::vector<ComplexD>>& lmd,
+                         std::vector<std::vector<ComplexD>>& lme, 
+			 int Nu, int Nk, int Nm,
+			 Eigen::MatrixXcd & Qt, // Nm x Nm
+			 GridBase *grid)
+  {
+    assert( Nk%Nu == 0 && Nm%Nu == 0 );
+    assert( Nk <= Nm );
+    Eigen::MatrixXcd BlockTriDiag = Eigen::MatrixXcd::Zero(Nk,Nk);
+    
+    for ( int u=0; u<Nu; ++u ) {
+      for (int k=0; k<Nk; ++k ) {
+        BlockTriDiag(k,u+(k/Nu)*Nu) = lmd[u][k];
+      }
+    }
+    
+    for ( int u=0; u<Nu; ++u ) {
+      for (int k=Nu; k<Nk; ++k ) {
+        BlockTriDiag(k-Nu,u+(k/Nu)*Nu) = conjugate(lme[u][k-Nu]);
+        BlockTriDiag(u+(k/Nu)*Nu,k-Nu) = lme[u][k-Nu];
+      }
+    }
+    //std::cout << BlockTriDiag << std::endl;
+    
+    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXcd> eigensolver(BlockTriDiag);
+
+    for (int i = 0; i < Nk; i++) {
+      eval[Nk-1-i] = eigensolver.eigenvalues()(i);
+    }
+    for (int i = 0; i < Nk; i++) {
+      for (int j = 0; j < Nk; j++) {
+	Qt(j,Nk-1-i) = eigensolver.eigenvectors()(j,i);
+	//Qt(Nk-1-i,j) = eigensolver.eigenvectors()(i,j);
+	//Qt(i,j) = eigensolver.eigenvectors()(i,j);
+      }
+    }
+  }
+
+
+  void diagonalize(std::vector<RealD>& eval, 
+                   std::vector<std::vector<ComplexD>>& lmd, 
+                   std::vector<std::vector<ComplexD>>& lme, 
+		   int Nu, int Nk, int Nm,   
+		   Eigen::MatrixXcd & Qt,
+		   GridBase *grid)
+  {
+    Qt = Eigen::MatrixXcd::Identity(Nm,Nm);
+    if ( diagonalisation == IRLdiagonaliseWithEigen ) { 
+      diagonalize_Eigen(eval,lmd,lme,Nu,Nk,Nm,Qt,grid);
+    } else { 
+      assert(0);
+    }
+  }
+  
+
+  void unpackHermitBlockTriDiagMatToEigen(
+         std::vector<std::vector<ComplexD>>& lmd,  
+         std::vector<std::vector<ComplexD>>& lme,
+         int Nu, int Nb, int Nk, int Nm,
+         Eigen::MatrixXcd& M)
+  {
+    //clog << "unpackHermitBlockTriDiagMatToEigen() begin" << '\n'; 
+    assert( Nk%Nu == 0 && Nm%Nu == 0 );
+    assert( Nk <= Nm );
+    M = Eigen::MatrixXcd::Zero(Nk,Nk);
+    
+    // rearrange 
+    for ( int u=0; u<Nu; ++u ) {
+      for (int k=0; k<Nk; ++k ) {
+        M(k,u+(k/Nu)*Nu) = lmd[u][k];
+      }
+    }
+
+    for ( int u=0; u<Nu; ++u ) {
+      for (int k=Nu; k<Nk; ++k ) {
+        M(k-Nu,u+(k/Nu)*Nu) = conjugate(lme[u][k-Nu]);
+        M(u+(k/Nu)*Nu,k-Nu) = lme[u][k-Nu];
+      }
+    }
+    //clog << "unpackHermitBlockTriDiagMatToEigen() end" << endl; 
+  }
+ 
+
+  void packHermitBlockTriDiagMatfromEigen(
+         std::vector<std::vector<ComplexD>>& lmd,
+         std::vector<std::vector<ComplexD>>& lme,
+         int Nu, int Nb, int Nk, int Nm,
+         Eigen::MatrixXcd& M)
+  {
+    //clog << "packHermitBlockTriDiagMatfromEigen() begin" << '\n'; 
+    assert( Nk%Nu == 0 && Nm%Nu == 0 );
+    assert( Nk <= Nm );
+    
+    // rearrange 
+    for ( int u=0; u<Nu; ++u ) {
+      for (int k=0; k<Nk; ++k ) {
+        lmd[u][k] = M(k,u+(k/Nu)*Nu);
+      }
+    }
+
+    for ( int u=0; u<Nu; ++u ) {
+      for (int k=Nu; k<Nk; ++k ) {
+        lme[u][k-Nu] = M(u+(k/Nu)*Nu,k-Nu);
+      }
+    }
+    //clog << "packHermitBlockTriDiagMatfromEigen() end" << endl; 
+  }
+
+
+//  void shiftedQRDecompEigen(Eigen::MatrixXcd& M, int Nm,
+//		            RealD Dsh,
+//		            Eigen::MatrixXcd& Qprod, int Nk)
+//  {
+//    //clog << "shiftedQRDecompEigen() begin" << '\n'; 
+//    Eigen::MatrixXcd Mtmp = Eigen::MatrixXcd::Zero(Nm,Nm);
+//    Eigen::MatrixXcd Q = Eigen::MatrixXcd::Zero(Nm,Nm);
+//     
+//    Mtmp = M;
+//    for (int i=0; i<Nm; ++i ) {
+//      Mtmp(i,i) = M(i,i) - Dsh;
+//    }
+//    
+//    Eigen::HouseholderQR<Eigen::MatrixXcd> QRD(Mtmp);
+//    Q = QRD.householderQ();
+//
+//    M = Q.adjoint()*(M*Q);
+//#if 0
+//    Qprod *= Q;
+//#else
+//    Mtmp = Qprod*Q;
+//    
+//    Eigen::HouseholderQR<Eigen::MatrixXcd> QRD2(Mtmp);
+//    Qprod = QRD2.householderQ();
+//
+//    Mtmp -= Qprod;
+//    clog << "Frobenius norm ||Qprod(after) - Qprod|| = " << Mtmp.norm() << std::endl;
+//#endif
+//    //clog << "shiftedQRDecompEigen() end" << endl; 
+//  }
+  void shiftedQRDecompEigen(Eigen::MatrixXcd& M, int Nm,
+		            RealD Dsh,
+		            Eigen::MatrixXcd& Qprod)
+  {
+    //clog << "shiftedQRDecompEigen() begin" << '\n'; 
+    Eigen::MatrixXcd Mtmp = Eigen::MatrixXcd::Zero(Nm,Nm);
+    //Eigen::MatrixXcd Qtmp = Eigen::MatrixXcd::Zero(Nm,Nm);
+    
+    Mtmp = Qprod.adjoint()*(M*Qprod);
+    for (int i=0; i<Nm; ++i ) {
+      for (int j=i+1; j<Nm; ++j ) {
+        Mtmp(i,j) = Mtmp(j,i);
+      }
+    }
+    
+    for (int i=0; i<Nm; ++i ) {
+      Mtmp(i,i) -= Dsh;
+      //Mtmp(i,i) = real(Mtmp(i,i)-Dsh);
+    }
+    
+    Eigen::HouseholderQR<Eigen::MatrixXcd> QRD(Mtmp);
+    //Qtmp = Qprod*QRD.householderQ();
+    
+    //Eigen::HouseholderQR<Eigen::MatrixXcd> QRD2(Qtmp);
+    //Qprod = QRD2.householderQ();
+
+    Qprod *= QRD.householderQ();
+    //ComplexD p;
+    //RealD r;
+
+    //r = 0.;
+    //for (int k=0; k<Nm; ++k) r += real(conj(Qprod(k,0))*Qprod(k,0));
+    //r = sqrt(r);
+    //for (int k=0; k<Nm; ++k) Qprod(k,0) /= r;
+    //
+    //for (int i=1; i<Nm; ++i) {
+    //  for (int j=0; j<i; ++j) {
+    //    p = 0.;
+    //    for (int k=0; k<Nm; ++k) {
+    //      p += conj(Qprod(k,j))*Qprod(k,i);
+    //    }
+    //    for (int k=0; k<Nm; ++k) {
+    //      Qprod(k,i) -= p*Qprod(k,j);
+    //    }
+    //  }
+    //  r = 0.;
+    //  for (int k=0; k<Nm; ++k) r += real(conj(Qprod(k,i))*Qprod(k,i));
+    //  r = sqrt(r);
+    //  for (int k=0; k<Nm; ++k) Qprod(k,i) /= r;
+    //}
+
+    //clog << "shiftedQRDecompEigen() end" << endl; 
+  }
+  
+
+  void exampleQRDecompEigen(void)
+  {
+    Eigen::MatrixXd A = Eigen::MatrixXd::Zero(3,3);
+    Eigen::MatrixXd Q = Eigen::MatrixXd::Zero(3,3);
+    Eigen::MatrixXd R = Eigen::MatrixXd::Zero(3,3);
+    Eigen::MatrixXd P = Eigen::MatrixXd::Zero(3,3);
+
+    A(0,0) = 12.0;
+    A(0,1) = -51.0;
+    A(0,2) = 4.0;
+    A(1,0) = 6.0;
+    A(1,1) = 167.0;
+    A(1,2) = -68.0;
+    A(2,0) = -4.0;
+    A(2,1) = 24.0;
+    A(2,2) = -41.0;
+    
+    clog << "matrix A before ColPivHouseholder" << std::endl;
+    for ( int i=0; i<3; i++ ) {
+      for ( int j=0; j<3; j++ ) {
+        clog << "A[" << i << "," << j << "] = " << A(i,j) << '\n';
+      }
+    }
+    clog << std::endl;
+
+    Eigen::ColPivHouseholderQR<Eigen::MatrixXd> QRD(A);
+    
+    clog << "matrix A after ColPivHouseholder" << std::endl;
+    for ( int i=0; i<3; i++ ) {
+      for ( int j=0; j<3; j++ ) {
+        clog << "A[" << i << "," << j << "] = " << A(i,j) << '\n';
+      }
+    }
+    clog << std::endl;
+    
+    clog << "HouseholderQ with sequence lenth = nonzeroPiviots" << std::endl;
+    Q = QRD.householderQ().setLength(QRD.nonzeroPivots());
+    for ( int i=0; i<3; i++ ) {
+      for ( int j=0; j<3; j++ ) {
+        clog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
+      }
+    }
+    clog << std::endl;
+    
+    clog << "HouseholderQ with sequence lenth = 1" << std::endl;
+    Q = QRD.householderQ().setLength(1);
+    for ( int i=0; i<3; i++ ) {
+      for ( int j=0; j<3; j++ ) {
+        clog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
+      }
+    }
+    clog << std::endl;
+    
+    clog << "HouseholderQ with sequence lenth = 2" << std::endl;
+    Q = QRD.householderQ().setLength(2);
+    for ( int i=0; i<3; i++ ) {
+      for ( int j=0; j<3; j++ ) {
+        clog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
+      }
+    }
+    clog << std::endl;
+    
+    clog << "matrixR" << std::endl;
+    R = QRD.matrixR();
+    for ( int i=0; i<3; i++ ) {
+      for ( int j=0; j<3; j++ ) {
+        clog << "R[" << i << "," << j << "] = " << R(i,j) << '\n';
+      }
+    }
+    clog << std::endl;
+
+    clog << "rank = " << QRD.rank() << std::endl;
+    clog << "threshold = " << QRD.threshold() << std::endl;
+    
+    clog << "matrixP" << std::endl;
+    P = QRD.colsPermutation();
+    for ( int i=0; i<3; i++ ) {
+      for ( int j=0; j<3; j++ ) {
+        clog << "P[" << i << "," << j << "] = " << P(i,j) << '\n';
+      }
+    }
+    clog << std::endl;
+
+
+    clog << "QR decomposition without column pivoting" << std::endl;
+    
+    A(0,0) = 12.0;
+    A(0,1) = -51.0;
+    A(0,2) = 4.0;
+    A(1,0) = 6.0;
+    A(1,1) = 167.0;
+    A(1,2) = -68.0;
+    A(2,0) = -4.0;
+    A(2,1) = 24.0;
+    A(2,2) = -41.0;
+    
+    clog << "matrix A before Householder" << std::endl;
+    for ( int i=0; i<3; i++ ) {
+      for ( int j=0; j<3; j++ ) {
+        clog << "A[" << i << "," << j << "] = " << A(i,j) << '\n';
+      }
+    }
+    clog << std::endl;
+    
+    Eigen::HouseholderQR<Eigen::MatrixXd> QRDplain(A);
+    
+    clog << "HouseholderQ" << std::endl;
+    Q = QRDplain.householderQ();
+    for ( int i=0; i<3; i++ ) {
+      for ( int j=0; j<3; j++ ) {
+        clog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
+      }
+    }
+    clog << std::endl;
+    
+    clog << "matrix A after Householder" << std::endl;
+    for ( int i=0; i<3; i++ ) {
+      for ( int j=0; j<3; j++ ) {
+        clog << "A[" << i << "," << j << "] = " << A(i,j) << '\n';
+      }
+    }
+    clog << std::endl;
+  }
+
+ };
+}
+#undef clog
+#endif

lib/algorithms/iterative/ImplicitlyRestartedLanczos.h.bak

@@ -0,0 +1,625 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/iterative/ImplicitlyRestartedLanczos.h
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Chulwoo Jung
+Author: Guido Cossu
+
+    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_IRL_H
+#define GRID_IRL_H
+
+#include <string.h> //memset
+
+namespace Grid {
+
+  enum IRLdiagonalisation { 
+    IRLdiagonaliseWithDSTEGR,
+    IRLdiagonaliseWithQR,
+    IRLdiagonaliseWithEigen
+  };
+
+////////////////////////////////////////////////////////////////////////////////
+// Helper class for sorting the evalues AND evectors by Field
+// Use pointer swizzle on vectors
+////////////////////////////////////////////////////////////////////////////////
+template<class Field>
+class SortEigen {
+ private:
+  static bool less_lmd(RealD left,RealD right){
+    return left > right;
+  }  
+  static bool less_pair(std::pair<RealD,Field const*>& left,
+                        std::pair<RealD,Field const*>& right){
+    return left.first > (right.first);
+  }  
+  
+ public:
+  void push(std::vector<RealD>& lmd,std::vector<Field>& evec,int N) {
+    
+    ////////////////////////////////////////////////////////////////////////
+    // PAB: FIXME: VERY VERY VERY wasteful: takes a copy of the entire vector set.
+    //    : The vector reorder should be done by pointer swizzle somehow
+    ////////////////////////////////////////////////////////////////////////
+    std::vector<Field> cpy(lmd.size(),evec[0]._grid);
+    for(int i=0;i<lmd.size();i++) cpy[i] = evec[i];
+    
+    std::vector<std::pair<RealD, Field const*> > emod(lmd.size());    
+
+    for(int i=0;i<lmd.size();++i)  emod[i] = std::pair<RealD,Field const*>(lmd[i],&cpy[i]);
+
+    partial_sort(emod.begin(),emod.begin()+N,emod.end(),less_pair);
+
+    typename std::vector<std::pair<RealD, Field const*> >::iterator it = emod.begin();
+    for(int i=0;i<N;++i){
+      lmd[i]=it->first;
+      evec[i]=*(it->second);
+      ++it;
+    }
+  }
+  void push(std::vector<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);
+  }
+};
+
+/////////////////////////////////////////////////////////////
+// Implicitly restarted lanczos
+/////////////////////////////////////////////////////////////
+template<class Field> 
+class ImplicitlyRestartedLanczos {
+
+private:       
+
+  int MaxIter;   // Max iterations
+  int Nstop;     // Number of evecs checked for convergence
+  int Nk;        // Number of converged sought
+  int Nm;        // Nm -- total number of vectors
+  RealD eresid;
+  IRLdiagonalisation diagonalisation;
+  ////////////////////////////////////
+  // Embedded objects
+  ////////////////////////////////////
+           SortEigen<Field> _sort;
+  LinearOperatorBase<Field> &_Linop;
+    OperatorFunction<Field> &_poly;
+
+  /////////////////////////
+  // Constructor
+  /////////////////////////
+public:       
+ ImplicitlyRestartedLanczos(LinearOperatorBase<Field> &Linop, // op
+			    OperatorFunction<Field> & poly,   // polynomial
+			    int _Nstop, // really sought vecs
+			    int _Nk,    // sought vecs
+			    int _Nm,    // total vecs
+			    RealD _eresid, // resid in lmd deficit 
+			    int _MaxIter,  // Max iterations
+			    IRLdiagonalisation _diagonalisation= IRLdiagonaliseWithEigen ) :
+    _Linop(Linop),    _poly(poly),
+      Nstop(_Nstop), Nk(_Nk), Nm(_Nm),
+      eresid(_eresid),  MaxIter(_MaxIter),
+      diagonalisation(_diagonalisation)
+      { };
+
+  ////////////////////////////////
+  // Helpers
+  ////////////////////////////////
+  static RealD normalise(Field& v) 
+  {
+    RealD nn = norm2(v);
+    nn = sqrt(nn);
+    v = v * (1.0/nn);
+    return nn;
+  }
+  
+  void orthogonalize(Field& w, std::vector<Field>& evec, int k)
+  {
+    typedef typename Field::scalar_type MyComplex;
+    MyComplex ip;
+    
+    for(int j=0; j<k; ++j){
+      ip = innerProduct(evec[j],w); 
+      w = w - ip * evec[j];
+    }
+    normalise(w);
+  }
+
+/* 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(std::vector<RealD>& eval,  std::vector<Field>& evec, const Field& src, int& Nconv)
+  {
+    
+    GridBase *grid = evec[0]._grid;
+    assert(grid == src._grid);
+    
+    std::cout << GridLogMessage <<"**************************************************************************"<< std::endl;
+    std::cout << GridLogMessage <<" ImplicitlyRestartedLanczos::calc() starting iteration 0 /  "<< MaxIter<< std::endl;
+    std::cout << GridLogMessage <<"**************************************************************************"<< std::endl;
+    std::cout << GridLogMessage <<" -- seek   Nk    = " << Nk    <<" vectors"<< std::endl;
+    std::cout << GridLogMessage <<" -- accept Nstop = " << Nstop <<" vectors"<< std::endl;
+    std::cout << GridLogMessage <<" -- total  Nm    = " << Nm    <<" vectors"<< std::endl;
+    std::cout << GridLogMessage <<" -- size of eval = " << eval.size() << std::endl;
+    std::cout << GridLogMessage <<" -- size of evec = " << evec.size() << std::endl;
+    if ( diagonalisation == IRLdiagonaliseWithDSTEGR ) {
+      std::cout << GridLogMessage << "Diagonalisation is DSTEGR "<<std::endl;
+    } else if ( diagonalisation == IRLdiagonaliseWithQR ) { 
+      std::cout << GridLogMessage << "Diagonalisation is QR "<<std::endl;
+    }  else if ( diagonalisation == IRLdiagonaliseWithEigen ) { 
+      std::cout << GridLogMessage << "Diagonalisation is Eigen "<<std::endl;
+    }
+    std::cout << GridLogMessage <<"**************************************************************************"<< std::endl;
+    
+    assert(Nm == evec.size() && Nm == eval.size());
+	
+    std::vector<RealD> lme(Nm);  
+    std::vector<RealD> lme2(Nm);
+    std::vector<RealD> eval2(Nm);
+
+    Eigen::MatrixXd    Qt = Eigen::MatrixXd::Zero(Nm,Nm);
+
+    std::vector<int>   Iconv(Nm);
+    std::vector<Field>  B(Nm,grid); // waste of space replicating
+    
+    Field f(grid);
+    Field v(grid);
+    
+    int k1 = 1;
+    int k2 = Nk;
+    
+    Nconv = 0;
+    
+    RealD beta_k;
+  
+    // Set initial vector
+    evec[0] = src;
+    std::cout << GridLogMessage <<"norm2(src)= " << norm2(src)<<std::endl;
+    
+    normalise(evec[0]);
+    std::cout << GridLogMessage <<"norm2(evec[0])= " << norm2(evec[0]) <<std::endl;
+    
+    // Initial Nk steps
+    for(int k=0; k<Nk; ++k) step(eval,lme,evec,f,Nm,k);
+    
+    // Restarting loop begins
+    int iter;
+    for(iter = 0; iter<MaxIter; ++iter){
+      
+      std::cout<< GridLogMessage <<" **********************"<< std::endl;
+      std::cout<< GridLogMessage <<" Restart iteration = "<< iter << std::endl;
+      std::cout<< GridLogMessage <<" **********************"<< std::endl;
+      
+      for(int k=Nk; k<Nm; ++k) step(eval,lme,evec,f,Nm,k);
+      
+      f *= lme[Nm-1];
+      
+      // getting eigenvalues
+      for(int k=0; k<Nm; ++k){
+	eval2[k] = eval[k+k1-1];
+	lme2[k] = lme[k+k1-1];
+      }
+      Qt = Eigen::MatrixXd::Identity(Nm,Nm);
+      diagonalize(eval2,lme2,Nm,Nm,Qt,grid);
+
+      // sorting
+      _sort.push(eval2,Nm);
+      
+      // Implicitly shifted QR transformations
+      Qt = Eigen::MatrixXd::Identity(Nm,Nm);
+      for(int ip=k2; ip<Nm; ++ip){ 
+	// Eigen replacement for qr_decomp ???
+	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 j=k1-1; j<k2+1; ++j){
+	for(int k=0; k<Nm; ++k){
+	  B[j].checkerboard = evec[k].checkerboard;
+	  B[j] += Qt(j,k) * evec[k];
+	}
+      }
+      for(int j=k1-1; j<k2+1; ++j) evec[j] = B[j];
+      
+      // Compressed vector f and beta(k2)
+      f *= Qt(k2-1,Nm-1);
+      f += lme[k2-1] * evec[k2];
+      beta_k = norm2(f);
+      beta_k = sqrt(beta_k);
+      std::cout<< GridLogMessage<<" 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<Nm; ++k){    
+	eval2[k] = eval[k];
+	lme2[k] = lme[k];
+      }
+      Qt = Eigen::MatrixXd::Identity(Nm,Nm);
+      diagonalize(eval2,lme2,Nk,Nm,Qt,grid);
+      
+      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){
+	  B[j].checkerboard = evec[k].checkerboard;
+	  B[j] += Qt(j,k) * evec[k];
+	}
+      }
+
+      Nconv = 0;
+      for(int i=0; i<Nk; ++i){
+	
+	_Linop.HermOp(B[i],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.precision(13);
+	std::cout << GridLogMessage << "[" << 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;
+	
+	// change the criteria as evals are supposed to be sorted, all evals smaller(larger) than Nstop should have converged
+	if((vv<eresid*eresid) && (i == Nconv) ){
+	  Iconv[Nconv] = i;
+	  ++Nconv;
+	}
+	
+      }  // i-loop end
+      
+      std::cout<< GridLogMessage <<" #modes converged: "<<Nconv<<std::endl;
+
+      if( Nconv>=Nstop ){
+	goto converged;
+      }
+    } // end of iter loop
+    
+    std::cout << GridLogMessage <<"**************************************************************************"<< std::endl;
+    std::cout<< GridLogError    <<" ImplicitlyRestartedLanczos::calc() NOT converged.";
+    std::cout << GridLogMessage <<"**************************************************************************"<< std::endl;
+    abort();
+	
+  converged:
+    // Sorting
+    eval.resize(Nconv);
+    evec.resize(Nconv,grid);
+    for(int i=0; i<Nconv; ++i){
+      eval[i] = eval2[Iconv[i]];
+      evec[i] = B[Iconv[i]];
+    }
+    _sort.push(eval,evec,Nconv);
+    
+    std::cout << GridLogMessage <<"**************************************************************************"<< std::endl;
+    std::cout << GridLogMessage << "ImplicitlyRestartedLanczos CONVERGED ; Summary :\n";
+    std::cout << GridLogMessage <<"**************************************************************************"<< std::endl;
+    std::cout << GridLogMessage << " -- Iterations  = "<< iter   << "\n";
+    std::cout << GridLogMessage << " -- beta(k)     = "<< beta_k << "\n";
+    std::cout << GridLogMessage << " -- Nconv       = "<< Nconv  << "\n";
+    std::cout << GridLogMessage <<"**************************************************************************"<< std::endl;
+  }
+
+private:
+/* 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(std::vector<RealD>& lmd,
+	    std::vector<RealD>& lme, 
+	    std::vector<Field>& evec,
+	    Field& w,int Nm,int k)
+  {
+    const RealD tiny = 1.0e-20;
+    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
+    
+    lmd[k] = alph;
+    lme[k] = beta;
+    
+    if ( k > 0 ) orthogonalize(w,evec,k); // orthonormalise
+    if ( k < Nm-1) evec[k+1] = w;
+    
+    if ( beta < tiny ) std::cout << GridLogMessage << " beta is tiny "<<beta<<std::endl;
+  }
+      
+  void diagonalize_Eigen(std::vector<RealD>& lmd, std::vector<RealD>& lme, 
+			 int Nk, int Nm,  
+			 Eigen::MatrixXd & Qt, // Nm x Nm
+			 GridBase *grid)
+  {
+    Eigen::MatrixXd TriDiag = Eigen::MatrixXd::Zero(Nk,Nk);
+
+    for(int i=0;i<Nk;i++)   TriDiag(i,i)   = lmd[i];
+    for(int i=0;i<Nk-1;i++) TriDiag(i,i+1) = lme[i];
+    for(int i=0;i<Nk-1;i++) TriDiag(i+1,i) = lme[i];
+    
+    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigensolver(TriDiag);
+
+    for (int i = 0; i < Nk; i++) {
+      lmd[Nk-1-i] = eigensolver.eigenvalues()(i);
+    }
+    for (int i = 0; i < Nk; i++) {
+      for (int j = 0; j < Nk; j++) {
+	Qt(Nk-1-i,j) = eigensolver.eigenvectors()(j,i);
+      }
+    }
+  }
+  ///////////////////////////////////////////////////////////////////////////
+  // File could end here if settle on Eigen ???
+  ///////////////////////////////////////////////////////////////////////////
+
+  void qr_decomp(std::vector<RealD>& lmd,   // Nm 
+		 std::vector<RealD>& lme,   // Nm 
+		 int Nk, int Nm,            // Nk, Nm
+		 Eigen::MatrixXd& Qt,       // Nm x Nm matrix
+		 RealD Dsh, int kmin, int kmax)
+  {
+    int k = kmin-1;
+    RealD x;
+    
+    RealD Fden = 1.0/hypot(lmd[k]-Dsh,lme[k]);
+    RealD c = ( lmd[k] -Dsh) *Fden;
+    RealD s = -lme[k] *Fden;
+      
+    RealD tmpa1 = lmd[k];
+    RealD tmpa2 = lmd[k+1];
+    RealD tmpb  = lme[k];
+
+    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];
+    lme[k+1] = c*lme[k+1];
+      
+    for(int i=0; i<Nk; ++i){
+      RealD Qtmp1 = Qt(k,i);
+      RealD Qtmp2 = Qt(k+1,i);
+      Qt(k,i)  = c*Qtmp1 - s*Qtmp2;
+      Qt(k+1,i)= s*Qtmp1 + c*Qtmp2; 
+    }
+
+    // Givens transformations
+    for(int k = kmin; k < kmax-1; ++k){
+      
+      RealD Fden = 1.0/hypot(x,lme[k-1]);
+      RealD c = lme[k-1]*Fden;
+      RealD s = - x*Fden;
+	
+      RealD tmpa1 = lmd[k];
+      RealD tmpa2 = lmd[k+1];
+      RealD tmpb  = lme[k];
+
+      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;
+      lme[k-1] = c*lme[k-1] -s*x;
+
+      if(k != kmax-2){
+	x = -s*lme[k+1];
+	lme[k+1] = c*lme[k+1];
+      }
+
+      for(int i=0; i<Nk; ++i){
+	RealD Qtmp1 = Qt(k,i);
+	RealD Qtmp2 = Qt(k+1,i);
+	Qt(k,i)     = c*Qtmp1 -s*Qtmp2;
+	Qt(k+1,i)   = s*Qtmp1 +c*Qtmp2;
+      }
+    }
+  }
+
+  void diagonalize(std::vector<RealD>& lmd, std::vector<RealD>& lme, 
+		   int Nk, int Nm,   
+		   Eigen::MatrixXd & Qt,
+		   GridBase *grid)
+  {
+    Qt = Eigen::MatrixXd::Identity(Nm,Nm);
+    if ( diagonalisation == IRLdiagonaliseWithDSTEGR ) {
+      diagonalize_lapack(lmd,lme,Nk,Nm,Qt,grid);
+    } else if ( diagonalisation == IRLdiagonaliseWithQR ) { 
+      diagonalize_QR(lmd,lme,Nk,Nm,Qt,grid);
+    }  else if ( diagonalisation == IRLdiagonaliseWithEigen ) { 
+      diagonalize_Eigen(lmd,lme,Nk,Nm,Qt,grid);
+    } else { 
+      assert(0);
+    }
+  }
+
+#ifdef USE_LAPACK
+void LAPACK_dstegr(char *jobz, char *range, int *n, double *d, double *e,
+                   double *vl, double *vu, int *il, int *iu, double *abstol,
+                   int *m, double *w, double *z, int *ldz, int *isuppz,
+                   double *work, int *lwork, int *iwork, int *liwork,
+                   int *info);
+#endif
+
+void diagonalize_lapack(std::vector<RealD>& lmd,
+			std::vector<RealD>& lme, 
+			int Nk, int Nm,  
+			Eigen::MatrixXd& Qt,
+			GridBase *grid)
+{
+#ifdef USE_LAPACK
+  const int size = Nm;
+  int NN = Nk;
+  double evals_tmp[NN];
+  double evec_tmp[NN][NN];
+  memset(evec_tmp[0],0,sizeof(double)*NN*NN);
+  double DD[NN];
+  double EE[NN];
+  for (int i = 0; i< NN; i++) {
+    for (int j = i - 1; j <= i + 1; j++) {
+      if ( j < NN && j >= 0 ) {
+	if (i==j) DD[i] = lmd[i];
+	if (i==j) evals_tmp[i] = lmd[i];
+	if (j==(i-1)) EE[j] = lme[j];
+      }
+    }
+  }
+  int evals_found;
+  int lwork = ( (18*NN) > (1+4*NN+NN*NN)? (18*NN):(1+4*NN+NN*NN)) ;
+  int liwork =  3+NN*10 ;
+  int iwork[liwork];
+  double work[lwork];
+  int isuppz[2*NN];
+  char jobz = 'V'; // calculate evals & evecs
+  char range = 'I'; // calculate all evals
+  //    char range = 'A'; // calculate all evals
+  char uplo = 'U'; // refer to upper half of original matrix
+  char compz = 'I'; // Compute eigenvectors of tridiagonal matrix
+  int ifail[NN];
+  int info;
+  int total = grid->_Nprocessors;
+  int node  = grid->_processor;
+  int interval = (NN/total)+1;
+  double vl = 0.0, vu = 0.0;
+  int il = interval*node+1 , iu = interval*(node+1);
+  if (iu > NN)  iu=NN;
+  double tol = 0.0;
+  if (1) {
+    memset(evals_tmp,0,sizeof(double)*NN);
+    if ( il <= NN){
+      LAPACK_dstegr(&jobz, &range, &NN,
+		    (double*)DD, (double*)EE,
+		    &vl, &vu, &il, &iu, // these four are ignored if second parameteris 'A'
+		    &tol, // tolerance
+		    &evals_found, evals_tmp, (double*)evec_tmp, &NN,
+		    isuppz,
+		    work, &lwork, iwork, &liwork,
+		    &info);
+      for (int i = iu-1; i>= il-1; i--){
+	evals_tmp[i] = evals_tmp[i - (il-1)];
+	if (il>1) evals_tmp[i-(il-1)]=0.;
+	for (int j = 0; j< NN; j++){
+	  evec_tmp[i][j] = evec_tmp[i - (il-1)][j];
+	  if (il>1) evec_tmp[i-(il-1)][j]=0.;
+	}
+      }
+    }
+    {
+      grid->GlobalSumVector(evals_tmp,NN);
+      grid->GlobalSumVector((double*)evec_tmp,NN*NN);
+    }
+  } 
+  // Safer to sort instead of just reversing it, 
+  // but the document of the routine says evals are sorted in increasing order. 
+  // qr gives evals in decreasing order.
+  for(int i=0;i<NN;i++){
+    lmd [NN-1-i]=evals_tmp[i];
+    for(int j=0;j<NN;j++){
+      Qt((NN-1-i),j)=evec_tmp[i][j];
+    }
+  }
+#else 
+  assert(0);
+#endif
+}
+
+  void diagonalize_QR(std::vector<RealD>& lmd, std::vector<RealD>& lme, 
+		      int Nk, int Nm,   
+		      Eigen::MatrixXd & Qt,
+		      GridBase *grid)
+  {
+    int Niter = 100*Nm;
+    int kmin = 1;
+    int kmax = Nk;
+
+    // (this should be more sophisticated)
+    for(int iter=0; iter<Niter; ++iter){
+      
+      // determination of 2x2 leading submatrix
+      RealD dsub = lmd[kmax-1]-lmd[kmax-2];
+      RealD dd = sqrt(dsub*dsub + 4.0*lme[kmax-2]*lme[kmax-2]);
+      RealD Dsh = 0.5*(lmd[kmax-2]+lmd[kmax-1] +dd*(dsub/fabs(dsub)));
+      // (Dsh: shift)
+	
+      // transformation
+      qr_decomp(lmd,lme,Nk,Nm,Qt,Dsh,kmin,kmax); // Nk, Nm
+	
+      // Convergence criterion (redef of kmin and kamx)
+      for(int j=kmax-1; j>= kmin; --j){
+	RealD dds = fabs(lmd[j-1])+fabs(lmd[j]);
+	if(fabs(lme[j-1])+dds > dds){
+	  kmax = j+1;
+	  goto continued;
+	}
+      }
+      Niter = iter;
+      return;
+
+    continued:
+      for(int j=0; j<kmax-1; ++j){
+	RealD dds = fabs(lmd[j])+fabs(lmd[j+1]);
+	if(fabs(lme[j])+dds > dds){
+	  kmin = j+1;
+	  break;
+	}
+      }
+    }
+    std::cout << GridLogError << "[QL method] Error - Too many iteration: "<<Niter<<"\n";
+    abort();
+  }
+
+ };
+}
+#endif

lib/algorithms/iterative/SimpleLanczos.h

@@ -0,0 +1,930 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/iterative/ImplicitlyRestartedLanczos.h
+
+    Copyright (C) 2015
+
+Author: Chulwoo Jung <chulwoo@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_LANC_H
+#define GRID_LANC_H
+
+#include <string.h>		//memset
+
+#ifdef USE_LAPACK
+#ifdef USE_MKL
+#include<mkl_lapack.h>
+#else
+void LAPACK_dstegr (char *jobz, char *range, int *n, double *d, double *e,
+		    double *vl, double *vu, int *il, int *iu, double *abstol,
+		    int *m, double *w, double *z, int *ldz, int *isuppz,
+		    double *work, int *lwork, int *iwork, int *liwork,
+		    int *info);
+//#include <lapacke/lapacke.h>
+#endif
+#endif
+
+#include <Grid/algorithms/densematrix/DenseMatrix.h>
+
+// eliminate temorary vector in calc()
+#define MEM_SAVE
+
+namespace Grid
+{
+
+  struct Bisection
+  {
+
+#if 0
+    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.;
+	    }
+	}
+    }
+#endif
+
+    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=%0.14e 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;
+	    }
+	  printf ("x0=%0.14e xu=%0.14e k=%d\n", x0, xu, k);
+	  x[k] = (x0 + xu) / 2;
+	}
+    }
+  };
+
+/////////////////////////////////////////////////////////////
+// Implicitly restarted lanczos
+/////////////////////////////////////////////////////////////
+
+
+  template < class Field > class SimpleLanczos
+  {
+
+    const RealD small = 1.0e-16;
+  public:
+    int lock;
+    int get;
+    int Niter;
+    int converged;
+
+    int Nstop;			// Number of evecs checked for convergence
+    int Nk;			// Number of converged sought
+    int Np;			// Np -- Number of spare vecs in kryloc space
+    int Nm;			// Nm -- total number of vectors
+
+
+    RealD OrthoTime;
+
+    RealD eresid;
+
+    SortEigen < Field > _sort;
+
+    LinearOperatorBase < Field > &_Linop;
+
+    OperatorFunction < Field > &_poly;
+
+    /////////////////////////
+    // Constructor
+    /////////////////////////
+    void init (void)
+    {
+    };
+    void Abort (int ff, DenseVector < RealD > &evals,
+		DenseVector < DenseVector < RealD > >&evecs);
+
+    SimpleLanczos (LinearOperatorBase < Field > &Linop,	// op
+		   OperatorFunction < Field > &poly,	// polynmial
+		   int _Nstop,	// sought vecs
+		   int _Nk,	// sought vecs
+		   int _Nm,	// spare vecs
+		   RealD _eresid,	// resid in lmdue deficit 
+		   int _Niter):	// Max iterations
+     
+      _Linop (Linop),
+      _poly (poly),
+      Nstop (_Nstop), Nk (_Nk), Nm (_Nm), eresid (_eresid), 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;
+
+      GridBase *grid = evec[0]._grid;
+      Field w (grid);
+      std::cout << GridLogMessage << "RitzMatrix " << std::endl;
+      for (int i = 0; i < k; i++)
+	{
+	  _Linop.HermOp (evec[i], w);
+//      _poly(_Linop,evec[i],w);
+	  std::cout << GridLogMessage << "[" << 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 << GridLogMessage << "oops" << std::endl;
+		      abort ();
+		    }
+		  else
+		    std::cout << GridLogMessage << " 0 ";
+		}
+	      else
+		{
+		  std::cout << GridLogMessage << " " << in << " ";
+		}
+	    }
+	  std::cout << GridLogMessage << std::endl;
+	}
+    }
+
+    void step (DenseVector < RealD > &lmd,
+	       DenseVector < RealD > &lme,
+	       Field & last, Field & current, Field & next, uint64_t k)
+    {
+      if (lmd.size () <= k)
+	lmd.resize (k + Nm);
+      if (lme.size () <= k)
+	lme.resize (k + Nm);
+
+
+//      _poly(_Linop,current,next );   // 3. wk:=Avk−βkv_{k−1}
+      _Linop.HermOp (current, next);	// 3. wk:=Avk−βkv_{k−1}
+      if (k > 0)
+	{
+	  next -= lme[k - 1] * last;
+	}
+//      std::cout<<GridLogMessage << "<last|next>" << innerProduct(last,next) <<std::endl;
+
+      ComplexD zalph = innerProduct (current, next);	// 4. αk:=(wk,vk)
+      RealD alph = real (zalph);
+
+      next = next - alph * current;	// 5. wk:=wk−αkvk
+//      std::cout<<GridLogMessage << "<current|next>" << innerProduct(current,next) <<std::endl;
+
+      RealD beta = normalise (next);	// 6. βk+1 := ∥wk∥2. If βk+1 = 0 then Stop
+      // 7. vk+1 := wk/βk+1
+//       norm=beta;
+
+      int interval = Nm / 100 + 1;
+      if ((k % interval) == 0)
+	std::
+	  cout << GridLogMessage << k << " : alpha = " << zalph << " beta " <<
+	  beta << std::endl;
+      const RealD tiny = 1.0e-20;
+      if (beta < tiny)
+	{
+	  std::cout << GridLogMessage << " beta is tiny " << beta << std::
+	    endl;
+	}
+      lmd[k] = alph;
+      lme[k] = beta;
+
+    }
+
+    void qr_decomp (DenseVector < RealD > &lmd,
+		    DenseVector < RealD > &lme,
+		    int Nk,
+		    int Nm,
+		    DenseVector < RealD > &Qt, RealD Dsh, int kmin, int kmax)
+    {
+      int k = kmin - 1;
+      RealD x;
+
+      RealD Fden = 1.0 / hypot (lmd[k] - Dsh, lme[k]);
+      RealD c = (lmd[k] - Dsh) * Fden;
+      RealD s = -lme[k] * Fden;
+
+      RealD tmpa1 = lmd[k];
+      RealD tmpa2 = lmd[k + 1];
+      RealD tmpb = lme[k];
+
+      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];
+      lme[k + 1] = c * lme[k + 1];
+
+      for (int i = 0; i < Nk; ++i)
+	{
+	  RealD Qtmp1 = Qt[i + Nm * k];
+	  RealD Qtmp2 = Qt[i + Nm * (k + 1)];
+	  Qt[i + Nm * k] = c * Qtmp1 - s * Qtmp2;
+	  Qt[i + Nm * (k + 1)] = s * Qtmp1 + c * Qtmp2;
+	}
+
+      // Givens transformations
+      for (int k = kmin; k < kmax - 1; ++k)
+	{
+
+	  RealD Fden = 1.0 / hypot (x, lme[k - 1]);
+	  RealD c = lme[k - 1] * Fden;
+	  RealD s = -x * Fden;
+
+	  RealD tmpa1 = lmd[k];
+	  RealD tmpa2 = lmd[k + 1];
+	  RealD tmpb = lme[k];
+
+	  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;
+	  lme[k - 1] = c * lme[k - 1] - s * x;
+
+	  if (k != kmax - 2)
+	    {
+	      x = -s * lme[k + 1];
+	      lme[k + 1] = c * lme[k + 1];
+	    }
+
+	  for (int i = 0; i < Nk; ++i)
+	    {
+	      RealD Qtmp1 = Qt[i + Nm * k];
+	      RealD Qtmp2 = Qt[i + Nm * (k + 1)];
+	      Qt[i + Nm * k] = c * Qtmp1 - s * Qtmp2;
+	      Qt[i + Nm * (k + 1)] = s * Qtmp1 + c * Qtmp2;
+	    }
+	}
+    }
+
+#ifdef USE_LAPACK
+#ifdef USE_MKL
+#define LAPACK_INT MKL_INT
+#else
+#define LAPACK_INT long long
+#endif
+    void diagonalize_lapack (DenseVector < RealD > &lmd, DenseVector < RealD > &lme, int N1,	// all
+			     int N2,	// get
+			     GridBase * grid)
+    {
+      const int size = Nm;
+      LAPACK_INT NN = N1;
+      double evals_tmp[NN];
+      double DD[NN];
+      double EE[NN];
+      for (int i = 0; i < NN; i++)
+	for (int j = i - 1; j <= i + 1; j++)
+	  if (j < NN && j >= 0)
+	    {
+	      if (i == j)
+		DD[i] = lmd[i];
+	      if (i == j)
+		evals_tmp[i] = lmd[i];
+	      if (j == (i - 1))
+		EE[j] = lme[j];
+	    }
+      LAPACK_INT evals_found;
+      LAPACK_INT lwork =
+	((18 * NN) >
+	 (1 + 4 * NN + NN * NN) ? (18 * NN) : (1 + 4 * NN + NN * NN));
+      LAPACK_INT liwork = 3 + NN * 10;
+      LAPACK_INT iwork[liwork];
+      double work[lwork];
+      LAPACK_INT isuppz[2 * NN];
+      char jobz = 'N';		// calculate evals only
+      char range = 'I';		// calculate il-th to iu-th evals
+      //    char range = 'A'; // calculate all evals
+      char uplo = 'U';		// refer to upper half of original matrix
+      char compz = 'I';		// Compute eigenvectors of tridiagonal matrix
+      int ifail[NN];
+      LAPACK_INT info;
+//  int total = QMP_get_number_of_nodes();
+//  int node = QMP_get_node_number();
+//  GridBase *grid = evec[0]._grid;
+      int total = grid->_Nprocessors;
+      int node = grid->_processor;
+      int interval = (NN / total) + 1;
+      double vl = 0.0, vu = 0.0;
+      LAPACK_INT il = interval * node + 1, iu = interval * (node + 1);
+      if (iu > NN)
+	iu = NN;
+      double tol = 0.0;
+      if (1)
+	{
+	  memset (evals_tmp, 0, sizeof (double) * NN);
+	  if (il <= NN)
+	    {
+	      printf ("total=%d node=%d il=%d iu=%d\n", total, node, il, iu);
+#ifdef USE_MKL
+	      dstegr (&jobz, &range, &NN,
+#else
+	      LAPACK_dstegr (&jobz, &range, &NN,
+#endif
+			     (double *) DD, (double *) EE, &vl, &vu, &il, &iu,	// these four are ignored if second parameteris 'A'
+			     &tol,	// tolerance
+			     &evals_found, evals_tmp, (double *) NULL, &NN,
+			     isuppz, work, &lwork, iwork, &liwork, &info);
+	      for (int i = iu - 1; i >= il - 1; i--)
+		{
+		  printf ("node=%d evals_found=%d evals_tmp[%d] = %g\n", node,
+			  evals_found, i - (il - 1), evals_tmp[i - (il - 1)]);
+		  evals_tmp[i] = evals_tmp[i - (il - 1)];
+		  if (il > 1)
+		    evals_tmp[i - (il - 1)] = 0.;
+		}
+	    }
+	  {
+	    grid->GlobalSumVector (evals_tmp, NN);
+	  }
+	}
+// cheating a bit. It is better to sort instead of just reversing it, but the document of the routine says evals are sorted in increasing order. qr gives evals in decreasing order.
+    }
+#undef LAPACK_INT
+#endif
+
+
+    void diagonalize (DenseVector < RealD > &lmd,
+		      DenseVector < RealD > &lme,
+		      int N2, int N1, GridBase * grid)
+    {
+
+#ifdef USE_LAPACK
+      const int check_lapack = 0;	// just use lapack if 0, check against lapack if 1
+
+      if (!check_lapack)
+	return diagonalize_lapack (lmd, lme, N2, N1, grid);
+
+//      diagonalize_lapack(lmd2,lme2,Nm2,Nm,Qt,grid);
+#endif
+    }
+
+#if 1
+    static RealD normalise (Field & v)
+    {
+      RealD nn = norm2 (v);
+      nn = sqrt (nn);
+      v = v * (1.0 / nn);
+      return nn;
+    }
+
+    void orthogonalize (Field & w, DenseVector < Field > &evec, int k)
+    {
+      double t0 = -usecond () / 1e6;
+      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);
+      t0 += usecond () / 1e6;
+      OrthoTime += t0;
+    }
+
+    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;
+    }
+
+
+    void calc (DenseVector < RealD > &eval, const Field & src, int &Nconv)
+    {
+
+      GridBase *grid = src._grid;
+//      assert(grid == src._grid);
+
+      std::
+	cout << GridLogMessage << " -- Nk = " << Nk << " Np = " << Np << std::
+	endl;
+      std::cout << GridLogMessage << " -- Nm = " << Nm << std::endl;
+      std::cout << GridLogMessage << " -- size of eval   = " << eval.
+	size () << std::endl;
+
+//      assert(c.size() && Nm == eval.size());
+
+      DenseVector < RealD > lme (Nm);
+      DenseVector < RealD > lmd (Nm);
+
+
+      Field current (grid);
+      Field last (grid);
+      Field next (grid);
+
+      Nconv = 0;
+
+      RealD beta_k;
+
+      // Set initial vector
+      // (uniform vector) Why not src??
+      //      evec[0] = 1.0;
+      current = src;
+      std::cout << GridLogMessage << "norm2(src)= " << norm2 (src) << std::
+	endl;
+      normalise (current);
+      std::
+	cout << GridLogMessage << "norm2(evec[0])= " << norm2 (current) <<
+	std::endl;
+
+      // Initial Nk steps
+      OrthoTime = 0.;
+      double t0 = usecond () / 1e6;
+      RealD norm;		// sqrt norm of last vector
+
+      uint64_t iter = 0;
+
+      bool initted = false;
+      std::vector < RealD > low (Nstop * 10);
+      std::vector < RealD > high (Nstop * 10);
+      RealD cont = 0.;
+      while (1) {
+	  cont = 0.;
+	  std::vector < RealD > lme2 (Nm);
+	  std::vector < RealD > lmd2 (Nm);
+	  for (uint64_t k = 0; k < Nm; ++k, iter++) {
+	      step (lmd, lme, last, current, next, iter);
+	      last = current;
+	      current = next;
+	    }
+	  double t1 = usecond () / 1e6;
+	  std::cout << GridLogMessage << "IRL::Initial steps: " << t1 -
+	    t0 << "seconds" << std::endl;
+	  t0 = t1;
+	  std::
+	    cout << GridLogMessage << "IRL::Initial steps:OrthoTime " <<
+	    OrthoTime << "seconds" << std::endl;
+
+	  // getting eigenvalues
+	  lmd2.resize (iter + 2);
+	  lme2.resize (iter + 2);
+	  for (uint64_t k = 0; k < iter; ++k) {
+	      lmd2[k + 1] = lmd[k];
+	      lme2[k + 2] = lme[k];
+	    }
+	  t1 = usecond () / 1e6;
+	  std::cout << GridLogMessage << "IRL:: copy: " << t1 -
+	    t0 << "seconds" << std::endl;
+	  t0 = t1;
+	  {
+	    int total = grid->_Nprocessors;
+	    int node = grid->_processor;
+	    int interval = (Nstop / total) + 1;
+	    int iu = (iter + 1) - (interval * node + 1);
+	    int il = (iter + 1) - (interval * (node + 1));
+	    std::vector < RealD > eval2 (iter + 3);
+	    RealD eps2;
+	    Bisection::bisec (lmd2, lme2, iter, il, iu, 1e-16, 1e-10, eval2,
+			      eps2);
+//        diagonalize(eval2,lme2,iter,Nk,grid);
+	    RealD diff = 0.;
+	    for (int i = il; i <= iu; i++) {
+		if (initted)
+		  diff =
+		    fabs (eval2[i] - high[iu-i]) / (fabs (eval2[i]) +
+						      fabs (high[iu-i]));
+		if (initted && (diff > eresid))
+		  cont = 1.;
+		if (initted)
+		  printf ("eval[%d]=%0.14e %0.14e, %0.14e\n", i, eval2[i],
+			  high[iu-i], diff);
+		high[iu-i] = eval2[i];
+	      }
+	    il = (interval * node + 1);
+	    iu = (interval * (node + 1));
+	    Bisection::bisec (lmd2, lme2, iter, il, iu, 1e-16, 1e-10, eval2,
+			      eps2);
+	    for (int i = il; i <= iu; i++) {
+		if (initted)
+		  diff =
+		    fabs (eval2[i] - low[i]) / (fabs (eval2[i]) +
+						fabs (low[i]));
+		if (initted && (diff > eresid))
+		  cont = 1.;
+		if (initted)
+		  printf ("eval[%d]=%0.14e %0.14e, %0.14e\n", i, eval2[i],
+			  low[i], diff);
+		low[i] = eval2[i];
+	      }
+	    t1 = usecond () / 1e6;
+	    std::cout << GridLogMessage << "IRL:: diagonalize: " << t1 -
+	      t0 << "seconds" << std::endl;
+	    t0 = t1;
+	  }
+
+	  for (uint64_t k = 0; k < Nk; ++k) {
+//          eval[k] = eval2[k];
+	    }
+	  if (initted)
+	    {
+	      grid->GlobalSumVector (&cont, 1);
+	      if (cont < 1.) return;
+	    }
+	  initted = true;
+	}
+
+    }
+
+
+
+
+
+
+/**
+   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 = fabs (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
+
+
+  };
+
+}
+#endif

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

lib/qcd/action/fermion/CayleyFermion5D.cc

@@ -444,6 +444,8 @@ void CayleyFermion5D<Impl>::SetCoefficientsInternal(RealD zolo_hi,std::vector<Co
     assert(omega[i]!=Coeff_t(0.0));
     bs[i] = 0.5*(bpc/omega[i] + bmc);
     cs[i] = 0.5*(bpc/omega[i] - bmc);
+    std::cout<<GridLogMessage << "CayleyFermion5D "<<i<<" bs="<<bs[i]<<" cs="<<cs[i]<< std::endl;
+
   }
 
   ////////////////////////////////////////////////////////

lib/qcd/action/fermion/CayleyFermion5Dvec.cc

@@ -470,6 +470,7 @@ void CayleyFermion5D<Impl>::MooeeInternalAsm(const FermionField &psi, FermionFie
 	a0 = a0+incr;
 	a1 = a1+incr;
 	a2 = a2+sizeof(typename Simd::scalar_type);
+	a2 = a2+sizeof(typename Simd::scalar_type); //ypj [debug]
       }}
     {
       int lexa = s1+LLs*site;
@@ -701,7 +702,8 @@ void CayleyFermion5D<Impl>::MooeeInternalZAsm(const FermionField &psi, FermionFi
 	}
 	a0 = a0+incr;
 	a1 = a1+incr;
-	a2 = a2+sizeof(typename Simd::scalar_type);
+	a2 = a2+sizeof(typename Simd::scalar_type); 
+	a2 = a2+sizeof(typename Simd::scalar_type); // ypj [debug]
       }}
     {
       int lexa = s1+LLs*site;

lib/qcd/action/fermion/DomainWallEOFAFermionvec.cc

@@ -475,7 +475,8 @@ namespace QCD {
                         }
                         a0 = a0 + incr;
                         a1 = a1 + incr;
-                        a2 = a2 + sizeof(typename Simd::scalar_type);
+                        a2 = a2 + sizeof(typename Simd::scalar_type); 
+                        a2 = a2 + sizeof(typename Simd::scalar_type); // ypj [debug]
                     }
                 }
 

lib/qcd/action/fermion/Fermion.h

@@ -63,8 +63,8 @@ Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
 #include <Grid/qcd/action/fermion/MobiusFermion.h>
 #include <Grid/qcd/action/fermion/MobiusEOFAFermion.h>
 #include <Grid/qcd/action/fermion/ZMobiusFermion.h>
-#include <Grid/qcd/action/fermion/SchurDiagTwoKappa.h>
 #include <Grid/qcd/action/fermion/ScaledShamirFermion.h>
+//#include <Grid/qcd/action/fermion/SchurDiagTwoKappa.h>
 #include <Grid/qcd/action/fermion/MobiusZolotarevFermion.h>
 #include <Grid/qcd/action/fermion/ShamirZolotarevFermion.h>
 #include <Grid/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h>

lib/qcd/action/fermion/MobiusEOFAFermionvec.cc

@@ -853,7 +853,8 @@ namespace QCD {
 
               a0 = a0 + incr;
               a1 = a1 + incr;
-              a2 = a2 + sizeof(typename Simd::scalar_type);
+              a2 = a2 + sizeof(typename Simd::scalar_type); 
+              a2 = a2 + sizeof(typename Simd::scalar_type); // ypj [debug]
             }
           }
 

lib/qcd/action/fermion/SchurDiagTwoKappa.h

@@ -1,3 +1,4 @@
+#if 1
     /*************************************************************************************
 
     Grid physics library, www.github.com/paboyle/Grid 
@@ -97,6 +98,117 @@ namespace Grid {
     }
   };
 
+#if 0
+  ///////////////////////////////////////////////////////////////////////////////////////////////////////
+  // Copied from DiagTwoSolve
+  ///////////////////////////////////////////////////////////////////////////////////////////////////////
+  template<class Field> class SchurRedBlackDiagTwoSolve {
+  private:
+    OperatorFunction<Field> & _HermitianRBSolver;
+    int CBfactorise;
+  public:
+
+    /////////////////////////////////////////////////////
+    // Wrap the usual normal equations Schur trick
+    /////////////////////////////////////////////////////
+  SchurRedBlackDiagTwoSolve(OperatorFunction<Field> &HermitianRBSolver)  :
+     _HermitianRBSolver(HermitianRBSolver) 
+    { 
+      CBfactorise=0;
+    };
+
+    template<class Matrix>
+      void operator() (Matrix & _Matrix,const Field &in, Field &out){
+
+      // FIXME CGdiagonalMee not implemented virtual function
+      // FIXME use CBfactorise to control schur decomp
+      GridBase *grid = _Matrix.RedBlackGrid();
+      GridBase *fgrid= _Matrix.Grid();
+
+      SchurDiagTwoOperator<Matrix,Field> _HermOpEO(_Matrix);
+ 
+      Field src_e(grid);
+      Field src_o(grid);
+      Field sol_e(grid);
+      Field sol_o(grid);
+      Field   tmp(grid);
+      Field  Mtmp(grid);
+      Field resid(fgrid);
+
+      pickCheckerboard(Even,src_e,in);
+      pickCheckerboard(Odd ,src_o,in);
+      pickCheckerboard(Even,sol_e,out);
+      pickCheckerboard(Odd ,sol_o,out);
+    
+      /////////////////////////////////////////////////////
+      // src_o = Mdag * (source_o - Moe MeeInv source_e)
+      /////////////////////////////////////////////////////
+      _Matrix.MooeeInv(src_e,tmp);     assert(  tmp.checkerboard ==Even);
+      _Matrix.Meooe   (tmp,Mtmp);      assert( Mtmp.checkerboard ==Odd);     
+      tmp=src_o-Mtmp;                  assert(  tmp.checkerboard ==Odd);     
+
+      // get the right MpcDag
+      _HermOpEO.MpcDag(tmp,src_o);     assert(src_o.checkerboard ==Odd);       
+
+      //////////////////////////////////////////////////////////////
+      // Call the red-black solver
+      //////////////////////////////////////////////////////////////
+      std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl;
+//      _HermitianRBSolver(_HermOpEO,src_o,sol_o);  assert(sol_o.checkerboard==Odd);
+      _HermitianRBSolver(_HermOpEO,src_o,tmp);  assert(tmp.checkerboard==Odd);
+      _Matrix.MooeeInv(tmp,sol_o);        assert(  sol_o.checkerboard   ==Odd);
+
+      ///////////////////////////////////////////////////
+      // sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
+      ///////////////////////////////////////////////////
+      _Matrix.Meooe(sol_o,tmp);        assert(  tmp.checkerboard   ==Even);
+      src_e = src_e-tmp;               assert(  src_e.checkerboard ==Even);
+      _Matrix.MooeeInv(src_e,sol_e);   assert(  sol_e.checkerboard ==Even);
+     
+      setCheckerboard(out,sol_e); assert(  sol_e.checkerboard ==Even);
+      setCheckerboard(out,sol_o); assert(  sol_o.checkerboard ==Odd );
+
+      // Verify the unprec residual
+      _Matrix.M(out,resid); 
+      resid = resid-in;
+      RealD ns = norm2(in);
+      RealD nr = norm2(resid);
+
+      std::cout<<GridLogMessage << "SchurRedBlackDiagTwoKappa solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
+    }     
+  };
+#endif
+namespace QCD{
+    //
+    // Determinant is det of middle factor
+    // This assumes Mee is indept of U.
+    //
+    //
+    template<class Impl>
+    class SchurDifferentiableDiagTwo:  public SchurDiagTwoOperator<FermionOperator<Impl>,typename Impl::FermionField> 
+      {
+      public:
+      INHERIT_IMPL_TYPES(Impl);
+
+ 	typedef FermionOperator<Impl> Matrix;
+
+	SchurDifferentiableDiagTwo (Matrix &Mat) : SchurDiagTwoOperator<Matrix,FermionField>(Mat) {};
+    };
+#if 0
+    template<class Impl>
+    class SchurDifferentiableDiagTwoKappa :  public SchurDiagTwoKappaOperator<FermionOperator<Impl>,typename Impl::FermionField> 
+      {
+      public:
+      INHERIT_IMPL_TYPES(Impl);
+
+ 	typedef FermionOperator<Impl> Matrix;
+
+	SchurDifferentiableDiagTwoKappa (Matrix &Mat) : SchurDiagTwoKappaOperator<Matrix,FermionField>(Mat) {};
+    };
+#endif
+}
+
 }
 
 #endif
+#endif

lib/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h

@@ -140,6 +140,7 @@ namespace Grid{
 
     };
 
+
   }
 }
 #endif

lib/qcd/utils/SpaceTimeGrid.cc

@@ -50,6 +50,7 @@ GridCartesian *SpaceTimeGrid::makeFourDimDWFGrid(const std::vector<int> & latt,c
 GridCartesian         *SpaceTimeGrid::makeFiveDimGrid(int Ls,const GridCartesian *FourDimGrid)
 {
   int N4=FourDimGrid->_ndimension;
+  assert(N4==4);
 
   std::vector<int> latt5(1,Ls);
   std::vector<int> simd5(1,1);

lib/simd/Grid_avx512.h

@@ -556,7 +556,7 @@ namespace Optimization {
     v3  = _mm256_add_epi32(v1, v2);
     v1  = _mm256_hadd_epi32(v3, v3);
     v2  = _mm256_hadd_epi32(v1, v1);
-    u1  = _mm256_castsi256_si128(v2);        // upper half
+    u1  = _mm256_castsi256_si128(v2);        // upper half ypj[debug] ; was missing
     u2  = _mm256_extracti128_si256(v2, 1);  // lower half
     ret = _mm_add_epi32(u1, u2);
     return _mm_cvtsi128_si32(ret);

lib/util/Init.cc

@@ -158,6 +158,14 @@ void GridCmdOptionInt(std::string &str,int & val)
   return;
 }
 
+// ypj [add]
+void GridCmdOptionFloat(std::string &str,double & val)
+{
+  std::stringstream ss(str);
+  ss>>val;
+  return;
+}
+
 
 void GridParseLayout(char **argv,int argc,
 		     std::vector<int> &latt,

lib/util/Init.h

@@ -54,6 +54,9 @@ namespace Grid {
   std::string GridCmdVectorIntToString(const std::vector<int> & vec);
   void GridCmdOptionCSL(std::string str,std::vector<std::string> & vec);
   void GridCmdOptionIntVector(std::string &str,std::vector<int> & vec);
+  // ypj [add]
+  void GridCmdOptionInt(std::string &str,int & val);
+  void GridCmdOptionFloat(std::string &str,double & val);
 
 
   void GridParseLayout(char **argv,int argc,

lib/util/Lexicographic.h

@@ -18,6 +18,10 @@ namespace Grid{
 
     static inline void IndexFromCoor (const std::vector<int>& coor,int &index,const std::vector<int> &dims){
       int nd=dims.size();
+      if(nd > coor.size())  {
+	std::cout<< "coor.size "<<coor.size()<<" >dims.size "<<dims.size()<<std::endl; 
+	assert(0);
+	}
       int stride=1;
       index=0;
       for(int d=0;d<nd;d++){

tests/lanczos/Test_dwf_block_lanczos.cc

@@ -0,0 +1,317 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/Test_dwf_block_lanczos.cc
+
+    Copyright (C) 2015
+
+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 */
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+//typedef typename GparityDomainWallFermionR::FermionField FermionField;
+typedef typename ZMobiusFermionR::FermionField FermionField;
+
+RealD AllZero(RealD x){ return 0.;}
+
+class CmdJobParams 
+{
+  public:
+    std::string gaugefile;
+
+    int Ls;
+    double mass;
+    double M5;
+    double mob_b;
+    std::vector<ComplexD> omega;
+    std::vector<Complex> boundary_phase;
+    
+    LanczosType Impl;
+    int Nu;
+    int Nk;
+    int Np;
+    int Nm;
+    int Nstop;
+    int Ntest;
+    int MaxIter;
+    double resid;
+    
+    double low;
+    double high;
+    int order;
+
+    CmdJobParams()
+      : gaugefile("Hot"),
+        Ls(8), mass(0.01), M5(1.8), mob_b(1.5),
+        Impl(LanczosType::irbl),
+        Nu(4), Nk(200), Np(200), Nstop(100), Ntest(1), MaxIter(10), resid(1.0e-8), 
+        low(0.2), high(5.5), order(11)
+    {Nm=Nk+Np;};
+    
+    void Parse(char **argv, int argc);
+};
+
+
+void CmdJobParams::Parse(char **argv,int argc)
+{
+  std::string arg;
+  std::vector<int> vi;
+  double re,im;
+  int expect, idx;
+  std::string vstr;
+  std::ifstream pfile;
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--gconf") ){
+    gaugefile = GridCmdOptionPayload(argv,argv+argc,"--gconf");
+  }
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--phase") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--phase");
+    pfile.open(arg);
+    assert(pfile);
+    expect = 0;
+    while( pfile >> vstr ) {
+      if ( vstr.compare("boundary_phase") == 0 ) {
+        pfile >> vstr;
+        GridCmdOptionInt(vstr,idx);
+        assert(expect==idx);
+        pfile >> vstr;
+        GridCmdOptionFloat(vstr,re);
+        pfile >> vstr;
+        GridCmdOptionFloat(vstr,im);
+        boundary_phase.push_back({re,im});
+        expect++;
+      }
+    }
+    pfile.close();
+  } else {
+    for (int i=0; i<4; ++i) boundary_phase.push_back({1.,0.});
+  }
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--omega") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--omega");
+    pfile.open(arg);
+    assert(pfile);
+    Ls = 0;
+    while( pfile >> vstr ) {
+      if ( vstr.compare("omega") == 0 ) {
+        pfile >> vstr;
+        GridCmdOptionInt(vstr,idx);
+        assert(Ls==idx);
+        pfile >> vstr;
+        GridCmdOptionFloat(vstr,re);
+        pfile >> vstr;
+        GridCmdOptionFloat(vstr,im);
+        omega.push_back({re,im});
+        Ls++;
+      }
+    }
+    pfile.close();
+  } else {
+    if( GridCmdOptionExists(argv,argv+argc,"--Ls") ){
+      arg = GridCmdOptionPayload(argv,argv+argc,"--Ls");
+      GridCmdOptionInt(arg,Ls);
+    }
+  }
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--mass") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--mass");
+    GridCmdOptionFloat(arg,mass);
+  }
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--M5") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--M5");
+    GridCmdOptionFloat(arg,M5);
+  }
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--mob_b") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--mob_b");
+    GridCmdOptionFloat(arg,mob_b);
+  }
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--irbl") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--irbl");
+    GridCmdOptionIntVector(arg,vi);
+    Nu = vi[0];
+    Nk = vi[1];
+    Np = vi[2];
+    Nstop = vi[3];
+    MaxIter = vi[4];
+    // ypj[fixme] mode overriding message is needed.
+    Impl = LanczosType::irbl;
+    Nm = Nk+Np;
+  }
+  
+  // block Lanczos with explicit extension of its dimensions
+  if( GridCmdOptionExists(argv,argv+argc,"--rbl") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--rbl");
+    GridCmdOptionIntVector(arg,vi);
+    Nu = vi[0];
+    Nk = vi[1];
+    Np = vi[2]; // vector space is enlarged by adding Np vectors
+    Nstop = vi[3];
+    MaxIter = vi[4];
+    // ypj[fixme] mode overriding message is needed.
+    Impl = LanczosType::rbl;
+    Nm = Nk+Np*MaxIter;
+  }
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--check_int") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--check_int");
+    GridCmdOptionInt(arg,Ntest);
+  }
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--resid") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--resid");
+    GridCmdOptionFloat(arg,resid);
+  }
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--cheby_l") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--cheby_l");
+    GridCmdOptionFloat(arg,low);
+  }
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--cheby_u") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--cheby_u");
+    GridCmdOptionFloat(arg,high);
+  }
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--cheby_n") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--cheby_n");
+    GridCmdOptionInt(arg,order);
+  }
+  
+  if ( CartesianCommunicator::RankWorld() == 0 ) {
+    std::streamsize ss = std::cout.precision();
+    std::cout << GridLogMessage <<" Gauge Configuration "<< gaugefile << '\n';
+    std::cout.precision(15);
+    for ( int i=0; i<4; ++i ) std::cout << GridLogMessage <<" boundary_phase["<< i << "] = " << boundary_phase[i] << '\n';
+    std::cout.precision(ss);
+    std::cout << GridLogMessage <<" Ls "<< Ls << '\n';
+    std::cout << GridLogMessage <<" mass "<< mass << '\n';
+    std::cout << GridLogMessage <<" M5 "<< M5 << '\n';
+    std::cout << GridLogMessage <<" mob_b "<< mob_b << '\n';
+    std::cout.precision(15);
+    for ( int i=0; i<Ls; ++i ) std::cout << GridLogMessage <<" omega["<< i << "] = " << omega[i] << '\n';
+    std::cout.precision(ss);
+    std::cout << GridLogMessage <<" Nu "<< Nu << '\n'; 
+    std::cout << GridLogMessage <<" Nk "<< Nk << '\n'; 
+    std::cout << GridLogMessage <<" Np "<< Np << '\n'; 
+    std::cout << GridLogMessage <<" Nm "<< Nm << '\n'; 
+    std::cout << GridLogMessage <<" Nstop "<< Nstop << '\n'; 
+    std::cout << GridLogMessage <<" Ntest "<< Ntest << '\n'; 
+    std::cout << GridLogMessage <<" MaxIter "<< MaxIter << '\n'; 
+    std::cout << GridLogMessage <<" resid "<< resid << '\n'; 
+    std::cout << GridLogMessage <<" Cheby Poly "<< low << "," << high << "," << order << std::endl; 
+  }
+}
+
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+  
+  CmdJobParams JP;
+  JP.Parse(argv,argc);
+
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(JP.Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(JP.Ls,UGrid);
+  printf("UGrid=%p UrbGrid=%p FGrid=%p FrbGrid=%p\n",UGrid,UrbGrid,FGrid,FrbGrid);
+
+  std::vector<int> seeds4({1,2,3,4});
+  std::vector<int> seeds5({5,6,7,8});
+  GridParallelRNG          RNG5(FGrid);  RNG5.SeedFixedIntegers(seeds5);
+  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
+  GridParallelRNG          RNG5rb(FrbGrid);  RNG5.SeedFixedIntegers(seeds5); 
+  // ypj [note] why seed RNG5 again? bug? In this case, run with a default seed().
+  //GridParallelRNG          RNG5rb(FrbGrid);  //RNG5rb.SeedFixedIntegers(seeds5);
+
+  LatticeGaugeField Umu(UGrid); 
+  std::vector<LatticeColourMatrix> U(4,UGrid);
+  
+  if ( JP.gaugefile.compare("Hot") == 0 ) {
+    SU3::HotConfiguration(RNG4, Umu);
+  } else {
+    FieldMetaData header;
+    NerscIO::readConfiguration(Umu,header,JP.gaugefile);
+    // ypj [fixme] additional checks for the loaded configuration?
+  }
+  
+  for(int mu=0;mu<Nd;mu++){
+    U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
+  }
+  
+  RealD mass = JP.mass;
+  RealD M5 = JP.M5;
+
+// ypj [fixme] flexible support for a various Fermions
+//  RealD mob_b = JP.mob_b;      // Gparity
+//  std::vector<ComplexD> omega; // ZMobius
+  
+//  GparityMobiusFermionD ::ImplParams params;
+//  std::vector<int> twists({1,1,1,0});
+//  params.twists = twists;
+//  GparityMobiusFermionR  Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,mob_b,mob_b-1.,params);
+//  SchurDiagTwoOperator<GparityMobiusFermionR,FermionField> HermOp(Ddwf);
+
+  //WilsonFermionR::ImplParams params;
+  ZMobiusFermionR::ImplParams params;
+  params.overlapCommsCompute = true;
+  params.boundary_phases = JP.boundary_phase;
+  ZMobiusFermionR  Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,JP.omega,1.,0.,params);
+  SchurDiagTwoOperator<ZMobiusFermionR,FermionField> HermOp(Ddwf);
+
+  //std::vector<double> Coeffs { 0.,-1.}; 
+  // ypj [note] this may not be supported by some compilers
+  std::vector<double> Coeffs({ 0.,-1.}); 
+  Polynomial<FermionField> PolyX(Coeffs);
+  //Chebyshev<FermionField> Cheb(0.2,5.5,11);
+  Chebyshev<FermionField> Cheb(JP.low,JP.high,JP.order);
+//  Cheb.csv(std::cout);
+  ImplicitlyRestartedBlockLanczos<FermionField> IRBL(HermOp,
+                                                     Cheb,
+                                                     JP.Nstop, JP.Ntest,
+                                                     JP.Nu, JP.Nk, JP.Nm,
+                                                     JP.resid,
+                                                     JP.MaxIter);
+  
+  std::vector<RealD> eval(JP.Nm);
+  
+  std::vector<FermionField> src(JP.Nu,FrbGrid);
+  for ( int i=0; i<JP.Nu; ++i ) gaussian(RNG5rb,src[i]);
+  
+  std::vector<FermionField> evec(JP.Nm,FrbGrid);
+  for(int i=0;i<1;++i){
+    std::cout << GridLogMessage << i <<" / "<< JP.Nm <<" grid pointer "<< evec[i]._grid << std::endl;
+  };
+
+  int Nconv;
+  IRBL.calc(eval,evec,src,Nconv,JP.Impl);
+
+
+  Grid_finalize();
+}

tests/lanczos/Test_dwf_compressed_lanczos.cc

@@ -22,6 +22,7 @@
 */
 #include <Grid/Grid.h>
 #include <Grid/algorithms/iterative/ImplicitlyRestartedLanczos.h>
+#include <Grid/parallelIO/IldgIO.h>
 /////////////////////////////////////////////////////////////////////////////
 // The following are now decoupled from the Lanczos and deal with grids.
 // Safe to replace functionality

tests/lanczos/Test_dwf_compressed_lanczos_reorg.cc

@@ -33,6 +33,8 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include <Grid/Grid.h>
 #include <Grid/algorithms/iterative/ImplicitlyRestartedLanczos.h>
 #include <Grid/algorithms/iterative/LocalCoherenceLanczos.h>
+#include <Grid/parallelIO/IldgIOtypes.h>
+#include <Grid/parallelIO/IldgIO.h>
 
 using namespace std;
 using namespace Grid;
@@ -57,7 +59,7 @@ public:
   void checkpointFine(std::string evecs_file,std::string evals_file)
   {
     assert(this->subspace.size()==nbasis);
-    emptyUserRecord record;
+    Grid::emptyUserRecord record;
     Grid::QCD::ScidacWriter WR(this->_FineGrid->IsBoss());
     WR.open(evecs_file);
     for(int k=0;k<nbasis;k++) {

tests/lanczos/Test_dwf_lanczos.cc

@@ -75,16 +75,16 @@ int main (int argc, char ** argv)
   SchurDiagTwoOperator<GparityMobiusFermionR,FermionField> HermOp(Ddwf);
 //  SchurDiagMooeeOperator<DomainWallFermionR,LatticeFermion> HermOp(Ddwf);
 
-  const int Nstop = 30;
-  const int Nk = 40;
-  const int Np = 40;
+  const int Nstop = 120;
+  const int Nk = 240;
+  const int Np = 240;
   const int Nm = Nk+Np;
-  const int MaxIt= 10000;
+  const int MaxIt= 10;
   RealD resid = 1.0e-8;
 
   std::vector<double> Coeffs { 0.,-1.};
   Polynomial<FermionField> PolyX(Coeffs);
-  Chebyshev<FermionField> Cheby(0.2,5.,11);
+  Chebyshev<FermionField> Cheby(0.2,5.5,11);
 
   FunctionHermOp<FermionField> OpCheby(Cheby,HermOp);
      PlainHermOp<FermionField> Op     (HermOp);

tests/lanczos/Test_wilson_lanczos.cc

@@ -58,7 +58,7 @@ int main(int argc, char** argv) {
   GridParallelRNG RNG4(UGrid);
   RNG4.SeedFixedIntegers(seeds4);
   GridParallelRNG RNG5rb(FrbGrid);
-  RNG5.SeedFixedIntegers(seeds5);
+  RNG5.SeedFixedIntegers(seeds5); // ypj [note] Does it mean RNG5rb? RNG5rb is never used.
 
   LatticeGaugeField Umu(UGrid);
   SU3::HotConfiguration(RNG4, Umu);

tests/solver/Test_dwf_mrhs_cg.cc

@@ -26,17 +26,18 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     *************************************************************************************/
     /*  END LEGAL */
 #include <Grid/Grid.h>
+#include <Grid/parallelIO/IldgIOtypes.h>
 #include <Grid/algorithms/iterative/BlockConjugateGradient.h>
 
-using namespace std;
+//using namespace std;
 using namespace Grid;
 using namespace Grid::QCD;
 
 int main (int argc, char ** argv)
 {
-  typedef typename DomainWallFermionR::FermionField FermionField; 
-  typedef typename DomainWallFermionR::ComplexField ComplexField; 
-  typename DomainWallFermionR::ImplParams params; 
+  typedef typename Grid::QCD::DomainWallFermionR::FermionField FermionField; 
+  typedef typename Grid::QCD::DomainWallFermionR::ComplexField ComplexField; 
+  typename Grid::QCD::DomainWallFermionR::ImplParams params; 
 
   const int Ls=4;
 
@@ -102,7 +103,7 @@ int main (int argc, char ** argv)
   std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
   std::cout << GridLogMessage << " Writing out in parallel view "<<std::endl;
   std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
-  emptyUserRecord record;
+//  Grid::emptyUserRecord record;
   std::string file("./scratch.scidac");
   std::string filef("./scratch.scidac.ferm");
 
@@ -115,20 +116,20 @@ int main (int argc, char ** argv)
   {
     FGrid->Barrier();
     ScidacWriter _ScidacWriter(FGrid->IsBoss());
-    _ScidacWriter.open(file);
+//    _ScidacWriter.open(file);
     std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
     std::cout << GridLogMessage << " Writing out gauge field "<<std::endl;
     std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
-    _ScidacWriter.writeScidacFieldRecord(Umu,record);
-    _ScidacWriter.close();
+//    _ScidacWriter.writeScidacFieldRecord(Umu,record);
+//    _ScidacWriter.close();
     FGrid->Barrier();
     std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
     std::cout << GridLogMessage << " Reading in gauge field "<<std::endl;
     std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
-    ScidacReader  _ScidacReader;
-    _ScidacReader.open(file);
-    _ScidacReader.readScidacFieldRecord(s_Umu,record);
-    _ScidacReader.close();
+//    ScidacReader  _ScidacReader;
+//    _ScidacReader.open(file);
+//    _ScidacReader.readScidacFieldRecord(s_Umu,record);
+//    _ScidacReader.close();
     FGrid->Barrier();
     std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
     std::cout << GridLogMessage << " Read in gauge field "<<std::endl;
@@ -145,9 +146,9 @@ int main (int argc, char ** argv)
 
       std::stringstream filefn;      filefn << filef << "."<< n;
       ScidacWriter _ScidacWriter(FGrid->IsBoss());
-      _ScidacWriter.open(filefn.str());
-      _ScidacWriter.writeScidacFieldRecord(src[n],record);
-      _ScidacWriter.close();
+//      _ScidacWriter.open(filefn.str());
+//      _ScidacWriter.writeScidacFieldRecord(src[n],record);
+//      _ScidacWriter.close();
     }
       
     FGrid->Barrier();
@@ -159,10 +160,10 @@ int main (int argc, char ** argv)
     for(int n=0;n<nrhs;n++){
       if ( n==me ) { 
 	std::stringstream filefn;	filefn << filef << "."<< n;
-	ScidacReader  _ScidacReader;
-	_ScidacReader.open(filefn.str());
-	_ScidacReader.readScidacFieldRecord(s_src,record);
-	_ScidacReader.close();
+//	ScidacReader  _ScidacReader;
+//	_ScidacReader.open(filefn.str());
+//	_ScidacReader.readScidacFieldRecord(s_src,record);
+//	_ScidacReader.close();
       }
     }
     FGrid->Barrier();