From fbe1209f7e94a8a4fb0615820c2f176838ad9c20 Mon Sep 17 00:00:00 2001
From: Yong-Chull Jang <integration.field@gmail.com>
Date: Wed, 31 Jan 2018 12:10:24 -0500
Subject: [PATCH] count converged eigenvalues not assuming candidates are
 sorted

---
 .../iterative/ImplicitlyRestartedBlockLanczos.h      | 12 +++++++++++-
 1 file changed, 11 insertions(+), 1 deletion(-)

diff --git a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h
index 933f2882..70c9afeb 100644
--- a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h
+++ b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h
@@ -166,6 +166,7 @@ until convergence
     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);
@@ -267,6 +268,7 @@ until convergence
 	eval2[i] = vnum/vden;
 	v -= eval2[i]*B[i];
 	RealD vv = norm2(v);
+        resid[i] = vv;
 	
 	std::cout.precision(13);
 	clog << "[" << std::setw(3)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
@@ -274,7 +276,8 @@ until convergence
 	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) && (i == Nconv) ){
+	if (vv<eresid*eresid) {
 	  Iconv[Nconv] = i;
 	  ++Nconv;
 	}
@@ -287,6 +290,13 @@ until convergence
 	goto converged;
       }
       
+      for(int i=0; i<Nconv; ++i){
+	std::cout.precision(13);
+	clog << "[" << std::setw(3)<< std::setiosflags(std::ios_base::right) <<Iconv[i]<<"] ";
+	std::cout << "eval_conv = "<<std::setw(25)<< std::setiosflags(std::ios_base::left)<< eval2[Iconv[i]];
+	std::cout << " |H B[i] - eval_conv[i]B[i]|^2 "<< std::setw(25)<< std::setiosflags(std::ios_base::right)<< resid[Iconv[i]]<< std::endl;
+      } 
+      
       if ( iter < MaxIter-1 ) {
         if ( Nu == 1 ) {
           // reconstruct initial vector for additional pole space