From 89c4e9b16863fcd684404cc405d26676aefc5331 Mon Sep 17 00:00:00 2001
From: Yong-Chull Jang <ypj@quark.phy.bnl.gov>
Date: Thu, 21 Dec 2017 23:13:39 -0500
Subject: [PATCH] first complete version of IRBL; requires practical test and
 clean up

---
 .../ImplicitlyRestartedBlockLanczos.h         | 31 ++++++++++---------
 tests/lanczos/Test_dwf_block_lanczos.cc       |  6 ++--
 tests/lanczos/Test_dwf_lanczos.cc             |  8 ++---
 3 files changed, 24 insertions(+), 21 deletions(-)

diff --git a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h
index 8d7a41c4..694f14fc 100644
--- a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h
+++ b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h
@@ -198,18 +198,7 @@ until convergence
       //  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){
@@ -293,7 +282,19 @@ until convergence
         //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
+      // 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
@@ -305,6 +306,8 @@ until convergence
       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);
 #endif
 
       // Convergence test
@@ -393,7 +396,7 @@ until convergence
     clog << fname + " CONVERGED ; Summary :\n";
     clog <<"**************************************************************************"<< std::endl;
     clog << " -- Iterations  = "<< iter   << "\n";
-    clog << " -- beta(k)     = "<< beta_k << "\n";
+    //clog << " -- beta(k)     = "<< beta_k << "\n";
     clog << " -- Nconv       = "<< Nconv  << "\n";
     clog <<"**************************************************************************"<< std::endl;
   }
diff --git a/tests/lanczos/Test_dwf_block_lanczos.cc b/tests/lanczos/Test_dwf_block_lanczos.cc
index 7050846a..92e3e289 100644
--- a/tests/lanczos/Test_dwf_block_lanczos.cc
+++ b/tests/lanczos/Test_dwf_block_lanczos.cc
@@ -77,10 +77,10 @@ int main (int argc, char ** argv)
   SchurDiagTwoOperator<GparityMobiusFermionR,FermionField> HermOp(Ddwf);
 //  SchurDiagMooeeOperator<DomainWallFermionR,LatticeFermion> HermOp(Ddwf);
 
-  const int Nstop = 50;
+  const int Nstop = 120;
   const int Nu = 1;
-  const int Nk = 200;
-  const int Np = 200;
+  const int Nk = 240;
+  const int Np = 240;
   const int Nm = Nk+Np;
   const int MaxIt= 10;
   RealD resid = 1.0e-8;
diff --git a/tests/lanczos/Test_dwf_lanczos.cc b/tests/lanczos/Test_dwf_lanczos.cc
index 6c65608d..58219007 100644
--- a/tests/lanczos/Test_dwf_lanczos.cc
+++ b/tests/lanczos/Test_dwf_lanczos.cc
@@ -75,11 +75,11 @@ int main (int argc, char ** argv)
   SchurDiagTwoOperator<GparityMobiusFermionR,FermionField> HermOp(Ddwf);
 //  SchurDiagMooeeOperator<DomainWallFermionR,LatticeFermion> HermOp(Ddwf);
 
-  const int Nstop = 50;
-  const int Nk = 200;
-  const int Np = 200;
+  const int Nstop = 120;
+  const int Nk = 240;
+  const int Np = 240;
   const int Nm = Nk+Np;
-  const int MaxIt= 100;
+  const int MaxIt= 10;
   RealD resid = 1.0e-8;
 
   std::vector<double> Coeffs { 0.,-1.};