diff --git a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h
index 71aa2891..2de35244 100644
--- a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h
+++ b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h
@@ -33,10 +33,12 @@ Author: Guido Cossu
 
 #include <string.h> //memset
 
-#define clog std::cout << GridLogMessage 
+#define Glog std::cout << GridLogMessage 
 
 namespace Grid {
 
+enum class LanczosType { irbl, rbl };
+
 /////////////////////////////////////////////////////////////
 // Implicitly restarted block lanczos
 /////////////////////////////////////////////////////////////
@@ -53,6 +55,7 @@ private:
   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;
   ////////////////////////////////////
@@ -69,6 +72,7 @@ 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
@@ -76,7 +80,8 @@ public:
                                  int _MaxIter,  // Max iterations
                                  IRLdiagonalisation _diagonalisation = IRLdiagonaliseWithEigen)
    : _Linop(Linop),    _poly(poly),
-      Nstop(_Nstop), Nu(_Nu), Nk(_Nk), Nm(_Nm), 
+      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),
@@ -117,30 +122,45 @@ public:
     }
     normalize(w);
   }
-
+  
   void calc(std::vector<RealD>& eval,  
             std::vector<Field>& evec, 
-            const std::vector<Field>& src, int& Nconv)
+            const std::vector<Field>& src, int& Nconv, LanczosType Impl)
   {
-    std::string fname = std::string(cname+"::calc()"); 
+    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() );
     
-    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;
+    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 ) { 
-      clog << "Diagonalisation is Eigen "<< std::endl;
+      Glog << "Diagonalisation is Eigen "<< std::endl;
     } else {
       abort();
     }
-    clog << std::string(74,'*') << std::endl;
+    Glog << std::string(74,'*') << std::endl;
     
     assert(Nm == evec.size() && Nm == eval.size());
 	
@@ -167,10 +187,10 @@ public:
   
     // set initial vector
     for (int i=0; i<Nu; ++i) {
-      clog << "norm2(src[" << i << "])= "<< norm2(src[i]) << std::endl;
+      Glog << "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;
+      Glog << "norm2(evec[" << i << "])= "<< norm2(evec[i]) << std::endl;
     }
     
     // initial Nblock_k steps
@@ -180,7 +200,7 @@ public:
     int iter;
     for(iter = 0; iter<MaxIter; ++iter){
       
-      clog <<"#Restart iteration = "<< iter << std::endl;
+      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);
       
@@ -194,7 +214,7 @@ public:
       Qt = Eigen::MatrixXcd::Identity(Nm,Nm);
       diagonalize(eval2,lmd2,lme2,Nu,Nm,Nm,Qt,grid);
       _sort.push(eval2,Nm);
-      clog << "#Ritz value before shift: "<< std::endl;
+      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<<"] ";
@@ -203,7 +223,7 @@ public:
       
       //----------------------------------------------------------------------
       if ( Nm>Nk ) {
-        clog <<" #Apply shifted QR transformations "<<std::endl;
+        Glog <<" #Apply shifted QR transformations "<<std::endl;
         //int k2 = Nk+Nu;
         int k2 = Nk;
       
@@ -240,7 +260,7 @@ public:
         Qt = Eigen::MatrixXcd::Identity(Nm,Nm);
         diagonalize(eval2,lmd2,lme2,Nu,Nk,Nm,Qt,grid);
         _sort.push(eval2,Nk);
-        clog << "#Ritz value after shift: "<< std::endl;
+        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<<"] ";
@@ -250,7 +270,7 @@ public:
       //----------------------------------------------------------------------
 
       // Convergence test
-      clog <<" #Convergence test: "<<std::endl;
+      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){
@@ -284,7 +304,7 @@ public:
 	
       }  // i-loop end
       
-      clog <<" #modes converged: "<<Nconv<<std::endl;
+      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]<<"] ";
@@ -296,11 +316,11 @@ public:
 
     } // end of iter loop
     
-    clog << std::string(74,'*') << std::endl;
+    Glog << std::string(74,'*') << std::endl;
     if ( Nconv<Nstop ) {
-      clog << fname + " NOT converged ; Summary :\n";
+      Glog << fname + " NOT converged ; Summary :\n";
     } else {
-      clog << fname + " CONVERGED ; Summary :\n";
+      Glog << fname + " CONVERGED ; Summary :\n";
       // Sort convered eigenpairs.
       eval.resize(Nconv);
       evec.resize(Nconv,grid);
@@ -310,11 +330,185 @@ public:
       }
       _sort.push(eval,evec,Nconv);
     }
-    clog << std::string(74,'*') << std::endl;
-    clog << " -- Iterations  = "<< iter   << "\n";
-    //clog << " -- beta(k)     = "<< beta_k << "\n";
-    clog << " -- Nconv       = "<< Nconv  << "\n";
-    clog << std::string(74,'*') << std::endl;
+    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;
+      for(int k = 0; k<Nr; ++k) B[k]=0.0;
+      for(int j = 0; j<Nr; j+=Nconv_test_interval){
+        Glog <<" #rotation for 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);
+	}
+      }
+      
+      Nconv = 0;
+      for(int i=0; i<Nr; i+=Nconv_test_interval){
+	
+        _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;
+      } 
+
+      (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;
   
   }
 
@@ -412,9 +606,9 @@ private:
     //lme[0][L] = beta;
     
     for (int u=0; u<Nu; ++u) {
-      clog << "norm2(w[" << u << "])= "<< norm2(w[u]) << std::endl;
+      Glog << "norm2(w[" << u << "])= "<< norm2(w[u]) << std::endl;
       for (int k=L+u; k<R; ++k) {
-        clog <<" In block "<< b << ","; 
+        Glog <<" In block "<< b << ","; 
         std::cout <<" beta[" << u << "," << k-L << "] = ";
         std::cout << lme[u][k] << std::endl;
       }
@@ -508,7 +702,7 @@ private:
          int Nu, int Nb, int Nk, int Nm,
          Eigen::MatrixXcd& M)
   {
-    //clog << "unpackHermitBlockTriDiagMatToEigen() begin" << '\n'; 
+    //Glog << "unpackHermitBlockTriDiagMatToEigen() begin" << '\n'; 
     assert( Nk%Nu == 0 && Nm%Nu == 0 );
     assert( Nk <= Nm );
     M = Eigen::MatrixXcd::Zero(Nk,Nk);
@@ -526,7 +720,7 @@ private:
         M(u+(k/Nu)*Nu,k-Nu) = lme[u][k-Nu];
       }
     }
-    //clog << "unpackHermitBlockTriDiagMatToEigen() end" << endl; 
+    //Glog << "unpackHermitBlockTriDiagMatToEigen() end" << endl; 
   }
  
 
@@ -536,7 +730,7 @@ private:
          int Nu, int Nb, int Nk, int Nm,
          Eigen::MatrixXcd& M)
   {
-    //clog << "packHermitBlockTriDiagMatfromEigen() begin" << '\n'; 
+    //Glog << "packHermitBlockTriDiagMatfromEigen() begin" << '\n'; 
     assert( Nk%Nu == 0 && Nm%Nu == 0 );
     assert( Nk <= Nm );
     
@@ -552,7 +746,7 @@ private:
         lme[u][k-Nu] = M(u+(k/Nu)*Nu,k-Nu);
       }
     }
-    //clog << "packHermitBlockTriDiagMatfromEigen() end" << endl; 
+    //Glog << "packHermitBlockTriDiagMatfromEigen() end" << endl; 
   }
 
 
@@ -561,7 +755,7 @@ private:
 		            RealD Dsh,
 		            Eigen::MatrixXcd& Qprod)
   {
-    //clog << "shiftedQRDecompEigen() begin" << '\n'; 
+    //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);
@@ -652,7 +846,7 @@ private:
     //  }
     //}
     
-    //clog << "shiftedQRDecompEigen() end" << endl; 
+    //Glog << "shiftedQRDecompEigen() end" << endl; 
   }
 
   void exampleQRDecompEigen(void)
@@ -672,74 +866,74 @@ private:
     A(2,1) = 24.0;
     A(2,2) = -41.0;
     
-    clog << "matrix A before ColPivHouseholder" << std::endl;
+    Glog << "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';
+        Glog << "A[" << i << "," << j << "] = " << A(i,j) << '\n';
       }
     }
-    clog << std::endl;
+    Glog << std::endl;
 
     Eigen::ColPivHouseholderQR<Eigen::MatrixXd> QRD(A);
     
-    clog << "matrix A after ColPivHouseholder" << std::endl;
+    Glog << "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';
+        Glog << "A[" << i << "," << j << "] = " << A(i,j) << '\n';
       }
     }
-    clog << std::endl;
+    Glog << std::endl;
     
-    clog << "HouseholderQ with sequence lenth = nonzeroPiviots" << 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++ ) {
-        clog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
+        Glog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
       }
     }
-    clog << std::endl;
+    Glog << std::endl;
     
-    clog << "HouseholderQ with sequence lenth = 1" << 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++ ) {
-        clog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
+        Glog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
       }
     }
-    clog << std::endl;
+    Glog << std::endl;
     
-    clog << "HouseholderQ with sequence lenth = 2" << 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++ ) {
-        clog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
+        Glog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
       }
     }
-    clog << std::endl;
+    Glog << std::endl;
     
-    clog << "matrixR" << std::endl;
+    Glog << "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';
+        Glog << "R[" << i << "," << j << "] = " << R(i,j) << '\n';
       }
     }
-    clog << std::endl;
+    Glog << std::endl;
 
-    clog << "rank = " << QRD.rank() << std::endl;
-    clog << "threshold = " << QRD.threshold() << std::endl;
+    Glog << "rank = " << QRD.rank() << std::endl;
+    Glog << "threshold = " << QRD.threshold() << std::endl;
     
-    clog << "matrixP" << std::endl;
+    Glog << "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';
+        Glog << "P[" << i << "," << j << "] = " << P(i,j) << '\n';
       }
     }
-    clog << std::endl;
+    Glog << std::endl;
 
 
-    clog << "QR decomposition without column pivoting" << std::endl;
+    Glog << "QR decomposition without column pivoting" << std::endl;
     
     A(0,0) = 12.0;
     A(0,1) = -51.0;
@@ -751,35 +945,35 @@ private:
     A(2,1) = 24.0;
     A(2,2) = -41.0;
     
-    clog << "matrix A before Householder" << std::endl;
+    Glog << "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';
+        Glog << "A[" << i << "," << j << "] = " << A(i,j) << '\n';
       }
     }
-    clog << std::endl;
+    Glog << std::endl;
     
     Eigen::HouseholderQR<Eigen::MatrixXd> QRDplain(A);
     
-    clog << "HouseholderQ" << std::endl;
+    Glog << "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';
+        Glog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
       }
     }
-    clog << std::endl;
+    Glog << std::endl;
     
-    clog << "matrix A after Householder" << std::endl;
+    Glog << "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';
+        Glog << "A[" << i << "," << j << "] = " << A(i,j) << '\n';
       }
     }
-    clog << std::endl;
+    Glog << std::endl;
   }
 
  };
 }
-#undef clog
+#undef Glog
 #endif
diff --git a/tests/lanczos/Test_dwf_block_lanczos.cc b/tests/lanczos/Test_dwf_block_lanczos.cc
index 141cf7bf..be306846 100644
--- a/tests/lanczos/Test_dwf_block_lanczos.cc
+++ b/tests/lanczos/Test_dwf_block_lanczos.cc
@@ -48,10 +48,13 @@ class CmdJobParams
     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;
     
@@ -61,10 +64,11 @@ class CmdJobParams
 
     CmdJobParams()
       : gaugefile("Hot"),
-        Ls(8), mass(0.01), M5(1.8), mob_b(1.5), 
-        Nu(4), Nk(200), Np(200), Nstop(100), MaxIter(10), resid(1.0e-8), 
+        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);
 };
@@ -155,6 +159,28 @@ void CmdJobParams::Parse(char **argv,int argc)
     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") ){
@@ -193,7 +219,9 @@ void CmdJobParams::Parse(char **argv,int argc)
     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; 
@@ -230,6 +258,7 @@ int main (int argc, char ** argv)
   } else {
     FieldMetaData header;
     NerscIO::readConfiguration(Umu,header,JP.gaugefile);
+    // ypj [fixme] additional checks for the loaded configuration?
   }
   
   for(int mu=0;mu<Nd;mu++){
@@ -238,6 +267,8 @@ int main (int argc, char ** argv)
   
   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
   
@@ -254,10 +285,6 @@ int main (int argc, char ** argv)
   ZMobiusFermionR  Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,JP.omega,1.,0.,params);
   SchurDiagTwoOperator<ZMobiusFermionR,FermionField> HermOp(Ddwf);
 
-  int Nu = JP.Nu;
-  int Nk = JP.Nk;
-  int Nm = Nk+JP.Np;
-
   //std::vector<double> Coeffs { 0.,-1.}; 
   // ypj [note] this may not be supported by some compilers
   std::vector<double> Coeffs({ 0.,-1.}); 
@@ -267,23 +294,23 @@ int main (int argc, char ** argv)
 //  Cheb.csv(std::cout);
   ImplicitlyRestartedBlockLanczos<FermionField> IRBL(HermOp,
                                                      Cheb,
-                                                     JP.Nstop,
-                                                     Nu,Nk,Nm,
+                                                     JP.Nstop, JP.Ntest,
+                                                     JP.Nu, JP.Nk, JP.Nm,
                                                      JP.resid,
                                                      JP.MaxIter);
   
-  std::vector<RealD> eval(Nm);
+  std::vector<RealD> eval(JP.Nm);
   
-  std::vector<FermionField> src(Nu,FrbGrid);
-  for ( int i=0; i<Nu; ++i ) gaussian(RNG5rb,src[i]);
+  std::vector<FermionField> src(JP.Nu,FrbGrid);
+  for ( int i=0; i<JP.Nu; ++i ) gaussian(RNG5rb,src[i]);
   
-  std::vector<FermionField> evec(Nm,FrbGrid);
+  std::vector<FermionField> evec(JP.Nm,FrbGrid);
   for(int i=0;i<1;++i){
-    std::cout << GridLogMessage << i <<" / "<< Nm <<" grid pointer "<< evec[i]._grid << std::endl;
+    std::cout << GridLogMessage << i <<" / "<< JP.Nm <<" grid pointer "<< evec[i]._grid << std::endl;
   };
 
   int Nconv;
-  IRBL.calc(eval,evec,src,Nconv);
+  IRBL.calc(eval,evec,src,Nconv,JP.Impl);
 
 
   Grid_finalize();