diff --git a/.gitignore b/.gitignore
index dc59879f..a2118a2d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -88,12 +88,18 @@ Thumbs.db
 ###################
 build*/*
 
+# bootstrap #
+#############
+*.tar.bz2*
+
 # IDE related files #
 #####################
 *.xcodeproj/*
 build.sh
 .vscode
 *.code-workspace
+.ctags
+tags
 
 # Eigen source #
 ################
diff --git a/bootstrap.sh b/bootstrap.sh
index f0a3a7f1..1510e3a1 100755
--- a/bootstrap.sh
+++ b/bootstrap.sh
@@ -1,10 +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}`
-#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
diff --git a/lib/algorithms/Algorithms.h b/lib/algorithms/Algorithms.h
index 09748241..58f6e1d4 100644
--- a/lib/algorithms/Algorithms.h
+++ b/lib/algorithms/Algorithms.h
@@ -50,6 +50,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #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>
 
diff --git a/lib/algorithms/approx/Chebyshev.h b/lib/algorithms/approx/Chebyshev.h
index 11d5c33f..b3364a7c 100644
--- a/lib/algorithms/approx/Chebyshev.h
+++ b/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(order,0.);
+      Coeffs.assign(order,0.);  
       Coeffs[order-1] = 1.;
     };
 
diff --git a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h
new file mode 100644
index 00000000..0ea4e09d
--- /dev/null
+++ b/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
diff --git a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h.bak b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h.bak
new file mode 100644
index 00000000..234d4c5c
--- /dev/null
+++ b/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
diff --git a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h.bak2 b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h.bak2
new file mode 100644
index 00000000..ad8a9c6d
--- /dev/null
+++ b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h.bak2
@@ -0,0 +1,1041 @@
+    /*************************************************************************************
+
+    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);
+  }
+  
+  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);
+  }
+
+/* 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';
+      //}
+      
+
+      
+      // 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;
+        //for (int i=0; i<Nm; ++i) {
+        //  for (int j=0; j<Nm; ++j) {
+        //    clog << "ckpt B2: M(" << ip << ")[" << i << "," << j << "] = " << BTDM(i,j) << '\n';
+        //  }
+        //}
+        //for (int i=0; i<Nm; ++i) {
+        //  for (int j=0; j<Nm; ++j) {
+        //    clog << "ckpt B3: Q(" << ip << ")[" << i << "," << j << "] = " << Q(i,j) << '\n';
+        //  }
+        //}
+        shiftedQRDecompEigen(BTDM,Nu,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) = conj(BTDM(j,i));
+      //  }
+      //  //BTDM(i,i) = real(BTDM(i,i));
+      //}
+      //for (int i=0; i<Nm; ++i) {
+      //  for (int j=0; j<Nm; ++j) {
+      //    clog << "ckpt B2: M(" << Nm-Nk << ")[" << i << "," << j << "] = " << BTDM(i,j) << '\n';
+      //  }
+      //}
+      //for (int i=0; i<Nm; ++i) {
+      //  for (int j=0; j<Nm; ++j) {
+      //    clog << "ckpt B3: Q(" << Nm-Nk << ")[" << i << "," << j << "] = " << Q(i,j) << '\n';
+      //  }
+      //}
+      
+      packHermitBlockTriDiagMatfromEigen(lmd,lme,Nu,Nblock_m,Nm,Nm,BTDM);
+
+      //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;
+      }
+
+      // residual vector
+#if 0 // 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;
+      }
+      
+      // 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;
+#else
+      blockwiseStep(lmd,lme,evec,f,f_copy,Nblock_k-1);
+      // set initial vector
+      //for (int i=0; i<Nu; ++i) {
+      //  orthogonalize(evec[i],evec,i);
+      //}
+      //// initial Nblock_k steps
+      //for(int b=0; b<Nblock_k; ++b) blockwiseStep(lmd,lme,evec,f,f_copy,b);
+#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 k=0; k<Nk; ++k) {
+      //  orthogonalize(B[k],B,k);
+      //}
+      
+      //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) {
+        for (int k=L-Nu+u; 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
+    //}
+    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
+    }
+    //for (int u=0; u<Nu; ++u) {
+    //  for (int k=L; k<L+u+1; ++k) {
+    //    clog << "ckpt B: alphaT[" << u+L << "," << k << "] = " << lmd[u][k] 
+    //         << ", alphaT - conj(alphaT) = " << lmd[u][k] - conjugate(lmd[k-L][u+L]) << 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;
+      }
+    }
+
+#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) {
+      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;
+      }
+    }
+#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)
+  {
+    //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; 
+  }
+
+#if 0
+  void shiftedQRDecompEigen(Eigen::MatrixXcd& M, int Nu, int Nm,
+		            RealD Dsh,
+		            Eigen::MatrixXcd& Qprod)
+  {
+    //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();
+    for (int j=0; j<Nm-2*Nu; ++j ) {
+      for (int i=2*Nu+j; i<Nm; ++i ) {
+        Q(i,j) = 0.;
+      }
+    }
+
+    M = Q.adjoint()*(M*Q);
+    for (int i=0; i<Nm; ++i ) {
+      for (int j=i+Nu; j<Nm; ++j ) {
+        M(i,j) = conj(M(j,i));
+      }
+    }
+    for (int i=0; i<Nm; ++i ) {
+      M(i,i) = real(M(i,i));
+    }
+
+    Qprod *= Q;
+    //clog << "shiftedQRDecompEigen() end" << endl; 
+  }
+#endif
+#if 1
+  // assume the input matrix M is a band matrix
+  void shiftedQRDecompEigen(Eigen::MatrixXcd& M, int Nu, int Nm,
+		            RealD Dsh,
+		            Eigen::MatrixXcd& Qprod)
+  {
+    //clog << "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) {
+      //for (int j=0; j<Nm-2*Nu; ++j) {
+      //  for (int k=0; k<j+2*Nu; ++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 j=Nm-2*Nu; 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)
+    //M = 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) { 
+    //      M(i,j) += R(i,k)*Q(k,j);
+    //    }
+    //    M(j,i) = conj(M(i,j));
+    //  }
+    //}
+
+    //for (int i=0; i<Nm; ++i) {
+    //  M(i,i) = real(M(i,i));
+    //}
+    
+    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.;
+      }
+    }
+    
+    //clog << "shiftedQRDecompEigen() end" << endl; 
+  }
+#endif
+#if 0
+  void shiftedQRDecompEigen(Eigen::MatrixXcd& M, int Nu, int Nm,
+		            RealD Dsh,
+		            Eigen::MatrixXcd& Qprod)
+  {
+    //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);
+    for (int i=0; i<Nm; ++i ) {
+      for (int j=i+1; j<Nm; ++j ) {
+        M(i,j) = conj(M(j,i));
+      }
+    }
+    for (int i=0; i<Nm; ++i ) {
+      M(i,i) = real(M(i,i));
+    }
+
+#if 1
+    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; 
+  }
+#endif
+#if 0
+  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; 
+  }
+#endif 
+
+  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
diff --git a/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h b/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h
index 787cf15a..9e711716 100644
--- a/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h
+++ b/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h
@@ -35,6 +35,8 @@ Author: Christoph Lehner <clehner@bnl.gov>
 //#include <zlib.h>
 #include <sys/stat.h>
 
+#define clog std::cout << GridLogMessage 
+
 namespace Grid { 
 
   ////////////////////////////////////////////////////////
@@ -401,7 +403,16 @@ until convergence
 
       std::cout<<GridLogIRL <<" running "<<Nm-Nk <<" steps: "<<std::endl;
       for(int k=Nk; k<Nm; ++k) step(eval,lme,evec,f,Nm,k);
+      for(int k=0; k<Nm; ++k) {
+        clog << "ckpt A1: lme[" << k << "] = " << lme[k] << '\n';
+      }
+      for(int k=0; k<Nm; ++k) {
+        clog << "ckpt A2: lmd[" << k << "] = " << eval[k] << '\n';
+      }
+      
       f *= lme[Nm-1];
+      //clog << "ckpt C " << '\n';
+      //clog << "norm2(f) = " << norm2(f) << std::endl;
 
       std::cout<<GridLogIRL <<" "<<Nm-Nk <<" steps done "<<std::endl;
       std::cout<<GridLogIRL <<"Initial steps:OrthoTime "<<OrthoTime<< "seconds"<<std::endl;
@@ -416,8 +427,12 @@ until convergence
       Qt = Eigen::MatrixXd::Identity(Nm,Nm);
       diagonalize(eval2,lme2,Nm,Nm,Qt,grid);
       std::cout<<GridLogIRL <<" diagonalized "<<std::endl;
-
+      
       //////////////////////////////////
+      //  clog << "ckpt D " << '\n';
+      //  clog << "eval2 [" << k << "] = " << eval2[k] << std::endl;
+      //}
+
       // sorting
       //////////////////////////////////
       eval2_copy = eval2;
@@ -434,24 +449,44 @@ until convergence
       }
 
       //////////////////////////////////
+      //  clog << "ckpt E " << '\n';
+      //  clog << "eval2 [" << k << "] = " << eval2[k] << std::endl;
+      //}
+      
       // Implicitly shifted QR transformations
       //////////////////////////////////
       Qt = Eigen::MatrixXd::Identity(Nm,Nm);
       for(int ip=k2; ip<Nm; ++ip){ 
 	QR_decomp(eval,lme,Nm,Nm,Qt,eval2[ip],k1,Nm);
+        clog << "ckpt B1: shift[" << ip << "] = " << eval2[ip] << std::endl;
       }
       std::cout<<GridLogIRL <<"QR decomposed "<<std::endl;
 
       assert(k2<Nm);      assert(k2<Nm);      assert(k1>0);
+      //  for (int j=0; j<Nm; ++j) {
+      //    clog << "ckpt G2: Q[" << i << "," << j << "] = " << Qt(j,i) << '\n';
+      //  }
+      //}
+      for (int i=0; i<Nm; ++i) {
+        clog << "ckpt C1: lme[" << i << "] = " << lme[i] << '\n';
+      }
+      for (int i=0; i<Nm; ++i) {
+        clog << "ckpt C2: lmd[" << i << "] = " << eval[i] << '\n';
+      }
+    
 
       basisRotate(evec,Qt,k1-1,k2+1,0,Nm,Nm); /// big constraint on the basis
       std::cout<<GridLogIRL <<"basisRotated  by Qt"<<std::endl;
+        evec[j] = B[j];
+        //clog << "ckpt F: norm2_evec[ " << j << "]" << norm2(evec[j]) << std::endl;
+      }
       
       ////////////////////////////////////////////////////
       // Compressed vector f and beta(k2)
       ////////////////////////////////////////////////////
       f *= Qt(k2-1,Nm-1);
-      f += lme[k2-1] * evec[k2];
+      f += lme[k2-1] * evec[k2]; // was commented out
+      std::cout<< GridLogMessage<<"ckpt D1: Q[Nm-1,Nk-1] = "<<Qt(Nk-1,Nm-1)<<std::endl;
       beta_k = norm2(f);
       beta_k = sqrt(beta_k);
       std::cout<<GridLogIRL<<" beta(k) = "<<beta_k<<std::endl;
@@ -471,6 +506,19 @@ until convergence
       diagonalize(eval2,lme2,Nk,Nm,Qt,grid);
       std::cout<<GridLogIRL <<" Diagonalized "<<std::endl;
 	  
+      //for (int i=0; i<Nk; ++i) {
+      //  for (int j=0; j<Nk; ++j) {
+      //    clog << "ckpt H1: R[" << i << "," << j << "] = " << Qt(j,i) << '\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;
       if (iter >= MinRestart) {
 
@@ -574,10 +622,17 @@ until convergence
 
     _PolyOp(evec_k,w);    std::cout<<GridLogIRL << "PolyOp" <<std::endl;
 
-    if(k>0) w -= lme[k-1] * evec[k-1];
+    if(k>0) { 
+      w -= lme[k-1] * evec[k-1];
+      //clog << "ckpt A (k= " << k << ")" << '\n';
+      //clog << "lme = " << lme[k-1] << '\n';
+      //clog << "norm(w) = " << norm2(w) << std::endl;
+    }
 
     ComplexD zalph = innerProduct(evec_k,w); // 4. αk:=(wk,vk)
     RealD     alph = real(zalph);
+    //clog << "ckpt B (k= " << k << ")" << '\n';
+    //clog << "lmd = " << alph << std::endl;
 
     w = w - alph * evec_k;// 5. wk:=wk−αkvk
 
@@ -838,4 +893,5 @@ void diagonalize_QR(std::vector<RealD>& lmd, std::vector<RealD>& lme,
 }
 };
 }
+#undef clog
 #endif
diff --git a/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h.bak b/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h.bak
new file mode 100644
index 00000000..a8723f32
--- /dev/null
+++ b/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
diff --git a/lib/qcd/action/fermion/CayleyFermion5Dvec.cc b/lib/qcd/action/fermion/CayleyFermion5Dvec.cc
index a7aa6b37..b7412696 100644
--- a/lib/qcd/action/fermion/CayleyFermion5Dvec.cc
+++ b/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;
@@ -702,6 +703,7 @@ 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); // ypj [debug]
       }}
     {
       int lexa = s1+LLs*site;
diff --git a/lib/qcd/action/fermion/DomainWallEOFAFermionvec.cc b/lib/qcd/action/fermion/DomainWallEOFAFermionvec.cc
index 68fd4fec..a3c0fb36 100644
--- a/lib/qcd/action/fermion/DomainWallEOFAFermionvec.cc
+++ b/lib/qcd/action/fermion/DomainWallEOFAFermionvec.cc
@@ -476,6 +476,7 @@ namespace QCD {
                         a0 = a0 + incr;
                         a1 = a1 + incr;
                         a2 = a2 + sizeof(typename Simd::scalar_type); 
+                        a2 = a2 + sizeof(typename Simd::scalar_type); // ypj [debug]
                     }
                 }
 
diff --git a/lib/qcd/action/fermion/MobiusEOFAFermionvec.cc b/lib/qcd/action/fermion/MobiusEOFAFermionvec.cc
index b7980660..a115969c 100644
--- a/lib/qcd/action/fermion/MobiusEOFAFermionvec.cc
+++ b/lib/qcd/action/fermion/MobiusEOFAFermionvec.cc
@@ -854,6 +854,7 @@ namespace QCD {
               a0 = a0 + incr;
               a1 = a1 + incr;
               a2 = a2 + sizeof(typename Simd::scalar_type); 
+              a2 = a2 + sizeof(typename Simd::scalar_type); // ypj [debug]
             }
           }
 
diff --git a/lib/simd/Grid_avx512.h b/lib/simd/Grid_avx512.h
index 4185b7d5..10ea6dc6 100644
--- a/lib/simd/Grid_avx512.h
+++ b/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);
diff --git a/lib/util/Init.cc b/lib/util/Init.cc
index fb3d7a1e..4ac9367b 100644
--- a/lib/util/Init.cc
+++ b/lib/util/Init.cc
@@ -157,6 +157,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,
diff --git a/lib/util/Init.h b/lib/util/Init.h
index 3da00742..e71c95c6 100644
--- a/lib/util/Init.h
+++ b/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,
diff --git a/tests/lanczos/Test_dwf_block_lanczos.cc b/tests/lanczos/Test_dwf_block_lanczos.cc
new file mode 100644
index 00000000..be306846
--- /dev/null
+++ b/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<ComplexD> 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();
+}
diff --git a/tests/lanczos/Test_dwf_lanczos.cc b/tests/lanczos/Test_dwf_lanczos.cc
index b1e205cf..8f00d68b 100644
--- a/tests/lanczos/Test_dwf_lanczos.cc
+++ b/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> Cheb(0.2,5.5,11);
 
   FunctionHermOp<FermionField> OpCheby(Cheby,HermOp);
      PlainHermOp<FermionField> Op     (HermOp);
diff --git a/tests/lanczos/Test_wilson_lanczos.cc b/tests/lanczos/Test_wilson_lanczos.cc
index eabc86d7..d876490d 100644
--- a/tests/lanczos/Test_wilson_lanczos.cc
+++ b/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);